Actual source code: ex36.c
1: static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
6: typedef struct _n_CCmplx CCmplx;
7: struct _n_CCmplx {
8: PetscReal real;
9: PetscReal imag;
10: };
12: CCmplx CCmplxPow(CCmplx a, PetscReal n)
13: {
14: CCmplx b;
15: PetscReal r, theta;
16: r = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
17: theta = PetscAtan2Real(a.imag, a.real);
18: b.real = PetscPowReal(r, n) * PetscCosReal(n * theta);
19: b.imag = PetscPowReal(r, n) * PetscSinReal(n * theta);
20: return b;
21: }
22: CCmplx CCmplxExp(CCmplx a)
23: {
24: CCmplx b;
25: b.real = PetscExpReal(a.real) * PetscCosReal(a.imag);
26: b.imag = PetscExpReal(a.real) * PetscSinReal(a.imag);
27: return b;
28: }
29: CCmplx CCmplxSqrt(CCmplx a)
30: {
31: CCmplx b;
32: PetscReal r, theta;
33: r = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
34: theta = PetscAtan2Real(a.imag, a.real);
35: b.real = PetscSqrtReal(r) * PetscCosReal(0.5 * theta);
36: b.imag = PetscSqrtReal(r) * PetscSinReal(0.5 * theta);
37: return b;
38: }
39: CCmplx CCmplxAdd(CCmplx a, CCmplx c)
40: {
41: CCmplx b;
42: b.real = a.real + c.real;
43: b.imag = a.imag + c.imag;
44: return b;
45: }
46: PetscScalar CCmplxRe(CCmplx a)
47: {
48: return a.real;
49: }
50: PetscScalar CCmplxIm(CCmplx a)
51: {
52: return a.imag;
53: }
55: PetscErrorCode DAApplyConformalMapping(DM da, PetscInt idx)
56: {
57: PetscInt i, n;
58: PetscInt sx, nx, sy, ny, sz, nz, dim;
59: Vec Gcoords;
60: PetscScalar *XX;
61: PetscScalar xx, yy, zz;
62: DM cda;
64: PetscFunctionBeginUser;
65: if (idx == 0) PetscFunctionReturn(PETSC_SUCCESS);
66: else if (idx == 1) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0)); /* dam break */
67: else if (idx == 2) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, 0.0, 1.0, -1.0, 1.0)); /* stagnation in a corner */
68: else if (idx == 3) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0)); /* nautilis */
69: else if (idx == 4) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
71: PetscCall(DMGetCoordinateDM(da, &cda));
72: PetscCall(DMGetCoordinates(da, &Gcoords));
74: PetscCall(VecGetArray(Gcoords, &XX));
75: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
76: PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
77: PetscCall(VecGetLocalSize(Gcoords, &n));
78: n = n / dim;
80: for (i = 0; i < n; i++) {
81: if (dim == 3 && idx != 2) {
82: PetscScalar Ni[8];
83: PetscScalar xi = XX[dim * i];
84: PetscScalar eta = XX[dim * i + 1];
85: PetscScalar zeta = XX[dim * i + 2];
86: PetscScalar xn[] = {-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0};
87: PetscScalar yn[] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0};
88: PetscScalar zn[] = {-0.1, -4.0, -0.2, -1.0, 0.1, 4.0, 0.2, 1.0};
90: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
91: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
92: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
93: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
95: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
96: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
97: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
98: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
100: xx = yy = zz = 0.0;
101: for (PetscInt p = 0; p < 8; p++) {
102: xx += Ni[p] * xn[p];
103: yy += Ni[p] * yn[p];
104: zz += Ni[p] * zn[p];
105: }
106: XX[dim * i] = xx;
107: XX[dim * i + 1] = yy;
108: XX[dim * i + 2] = zz;
109: }
111: if (idx == 1) {
112: CCmplx zeta, t1, t2;
114: xx = XX[dim * i] - 0.8;
115: yy = XX[dim * i + 1] + 1.5;
117: zeta.real = PetscRealPart(xx);
118: zeta.imag = PetscRealPart(yy);
120: t1 = CCmplxPow(zeta, -1.0);
121: t2 = CCmplxAdd(zeta, t1);
123: XX[dim * i] = CCmplxRe(t2);
124: XX[dim * i + 1] = CCmplxIm(t2);
125: } else if (idx == 2) {
126: CCmplx zeta, t1;
128: xx = XX[dim * i];
129: yy = XX[dim * i + 1];
130: zeta.real = PetscRealPart(xx);
131: zeta.imag = PetscRealPart(yy);
133: t1 = CCmplxSqrt(zeta);
134: XX[dim * i] = CCmplxRe(t1);
135: XX[dim * i + 1] = CCmplxIm(t1);
136: } else if (idx == 3) {
137: CCmplx zeta, t1, t2;
139: xx = XX[dim * i] - 0.8;
140: yy = XX[dim * i + 1] + 1.5;
142: zeta.real = PetscRealPart(xx);
143: zeta.imag = PetscRealPart(yy);
144: t1 = CCmplxPow(zeta, -1.0);
145: t2 = CCmplxAdd(zeta, t1);
146: XX[dim * i] = CCmplxRe(t2);
147: XX[dim * i + 1] = CCmplxIm(t2);
149: xx = XX[dim * i];
150: yy = XX[dim * i + 1];
151: zeta.real = PetscRealPart(xx);
152: zeta.imag = PetscRealPart(yy);
153: t1 = CCmplxExp(zeta);
154: XX[dim * i] = CCmplxRe(t1);
155: XX[dim * i + 1] = CCmplxIm(t1);
157: xx = XX[dim * i] + 0.4;
158: yy = XX[dim * i + 1];
159: zeta.real = PetscRealPart(xx);
160: zeta.imag = PetscRealPart(yy);
161: t1 = CCmplxPow(zeta, 2.0);
162: XX[dim * i] = CCmplxRe(t1);
163: XX[dim * i + 1] = CCmplxIm(t1);
164: } else if (idx == 4) {
165: PetscScalar Ni[4];
166: PetscScalar xi = XX[dim * i];
167: PetscScalar eta = XX[dim * i + 1];
168: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5};
169: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0};
171: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
172: Ni[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
173: Ni[2] = 0.25 * (1.0 - xi) * (1.0 + eta);
174: Ni[3] = 0.25 * (1.0 + xi) * (1.0 + eta);
176: xx = yy = 0.0;
177: for (PetscInt p = 0; p < 4; p++) {
178: xx += Ni[p] * xn[p];
179: yy += Ni[p] * yn[p];
180: }
181: XX[dim * i] = xx;
182: XX[dim * i + 1] = yy;
183: }
184: }
185: PetscCall(VecRestoreArray(Gcoords, &XX));
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: PetscErrorCode DAApplyTrilinearMapping(DM da)
190: {
191: PetscInt i, j, k;
192: PetscInt sx, nx, sy, ny, sz, nz;
193: Vec Gcoords;
194: DMDACoor3d ***XX;
195: PetscScalar xx, yy, zz;
196: DM cda;
198: PetscFunctionBeginUser;
199: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
200: PetscCall(DMGetCoordinateDM(da, &cda));
201: PetscCall(DMGetCoordinates(da, &Gcoords));
203: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
204: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
206: for (i = sx; i < sx + nx; i++) {
207: for (j = sy; j < sy + ny; j++) {
208: for (k = sz; k < sz + nz; k++) {
209: PetscScalar Ni[8];
210: PetscScalar xi = XX[k][j][i].x;
211: PetscScalar eta = XX[k][j][i].y;
212: PetscScalar zeta = XX[k][j][i].z;
213: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5, 0.0, 2.1, 0.23, 3.125};
214: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0, -1.45, -0.1, 2.24, 3.79};
215: PetscScalar zn[] = {0.0, 0.3, -0.1, 0.123, 0.956, 1.32, 1.12, 0.798};
217: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
218: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
219: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
220: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
222: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
223: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
224: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
225: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
227: xx = yy = zz = 0.0;
228: for (PetscInt p = 0; p < 8; p++) {
229: xx += Ni[p] * xn[p];
230: yy += Ni[p] * yn[p];
231: zz += Ni[p] * zn[p];
232: }
233: XX[k][j][i].x = xx;
234: XX[k][j][i].y = yy;
235: XX[k][j][i].z = zz;
236: }
237: }
238: }
239: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: PetscErrorCode DADefineXLinearField2D(DM da, Vec field)
244: {
245: PetscInt i, j;
246: PetscInt sx, nx, sy, ny;
247: Vec Gcoords;
248: DMDACoor2d **XX;
249: PetscScalar **FF;
250: DM cda;
252: PetscFunctionBeginUser;
253: PetscCall(DMGetCoordinateDM(da, &cda));
254: PetscCall(DMGetCoordinates(da, &Gcoords));
256: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
257: PetscCall(DMDAVecGetArray(da, field, &FF));
259: PetscCall(DMDAGetCorners(da, &sx, &sy, 0, &nx, &ny, 0));
261: for (i = sx; i < sx + nx; i++) {
262: for (j = sy; j < sy + ny; j++) FF[j][i] = 10.0 + 3.0 * XX[j][i].x + 5.5 * XX[j][i].y + 8.003 * XX[j][i].x * XX[j][i].y;
263: }
265: PetscCall(DMDAVecRestoreArray(da, field, &FF));
266: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscErrorCode DADefineXLinearField3D(DM da, Vec field)
271: {
272: PetscInt i, j, k;
273: PetscInt sx, nx, sy, ny, sz, nz;
274: Vec Gcoords;
275: DMDACoor3d ***XX;
276: PetscScalar ***FF;
277: DM cda;
279: PetscFunctionBeginUser;
280: PetscCall(DMGetCoordinateDM(da, &cda));
281: PetscCall(DMGetCoordinates(da, &Gcoords));
283: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
284: PetscCall(DMDAVecGetArray(da, field, &FF));
286: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
288: for (k = sz; k < sz + nz; k++) {
289: for (j = sy; j < sy + ny; j++) {
290: for (i = sx; i < sx + nx; i++) {
291: FF[k][j][i] = 10.0 + 4.05 * XX[k][j][i].x + 5.50 * XX[k][j][i].y + 1.33 * XX[k][j][i].z + 2.03 * XX[k][j][i].x * XX[k][j][i].y + 0.03 * XX[k][j][i].x * XX[k][j][i].z + 0.83 * XX[k][j][i].y * XX[k][j][i].z +
292: 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
293: }
294: }
295: }
297: PetscCall(DMDAVecRestoreArray(da, field, &FF));
298: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
303: {
304: DM dac, daf;
305: PetscViewer vv;
306: Vec ac, af;
307: PetscInt Mx;
308: Mat II, INTERP;
309: Vec scale;
310: PetscBool output = PETSC_FALSE;
312: PetscFunctionBeginUser;
313: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
314: PetscCall(DMSetFromOptions(dac));
315: PetscCall(DMSetUp(dac));
317: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
318: PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
319: Mx--;
321: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
322: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
324: {
325: DM cdaf, cdac;
326: Vec coordsc, coordsf;
328: PetscCall(DMGetCoordinateDM(dac, &cdac));
329: PetscCall(DMGetCoordinateDM(daf, &cdaf));
331: PetscCall(DMGetCoordinates(dac, &coordsc));
332: PetscCall(DMGetCoordinates(daf, &coordsf));
334: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
335: PetscCall(MatInterpolate(II, coordsc, coordsf));
336: PetscCall(MatDestroy(&II));
337: PetscCall(VecDestroy(&scale));
338: }
340: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
342: PetscCall(DMCreateGlobalVector(dac, &ac));
343: PetscCall(VecSet(ac, 66.99));
345: PetscCall(DMCreateGlobalVector(daf, &af));
347: PetscCall(MatMult(INTERP, ac, af));
349: {
350: Vec afexact;
351: PetscReal nrm;
352: PetscInt N;
354: PetscCall(DMCreateGlobalVector(daf, &afexact));
355: PetscCall(VecSet(afexact, 66.99));
356: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
357: PetscCall(VecNorm(afexact, NORM_2, &nrm));
358: PetscCall(VecGetSize(afexact, &N));
359: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
360: PetscCall(VecDestroy(&afexact));
361: }
363: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
364: if (output) {
365: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
366: PetscCall(VecView(ac, vv));
367: PetscCall(PetscViewerDestroy(&vv));
369: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
370: PetscCall(VecView(af, vv));
371: PetscCall(PetscViewerDestroy(&vv));
372: }
374: PetscCall(MatDestroy(&INTERP));
375: PetscCall(DMDestroy(&dac));
376: PetscCall(DMDestroy(&daf));
377: PetscCall(VecDestroy(&ac));
378: PetscCall(VecDestroy(&af));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
383: {
384: DM dac, daf;
385: PetscViewer vv;
386: Vec ac, af;
387: PetscInt map_id, Mx, My;
388: Mat II, INTERP;
389: Vec scale;
390: PetscBool output = PETSC_FALSE;
392: PetscFunctionBeginUser;
393: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, NULL, &dac));
394: PetscCall(DMSetFromOptions(dac));
395: PetscCall(DMSetUp(dac));
397: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
398: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
399: Mx--;
400: My--;
402: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
403: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
405: /* apply conformal mappings */
406: map_id = 0;
407: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
408: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
410: {
411: DM cdaf, cdac;
412: Vec coordsc, coordsf;
414: PetscCall(DMGetCoordinateDM(dac, &cdac));
415: PetscCall(DMGetCoordinateDM(daf, &cdaf));
417: PetscCall(DMGetCoordinates(dac, &coordsc));
418: PetscCall(DMGetCoordinates(daf, &coordsf));
420: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
421: PetscCall(MatInterpolate(II, coordsc, coordsf));
422: PetscCall(MatDestroy(&II));
423: PetscCall(VecDestroy(&scale));
424: }
426: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
428: PetscCall(DMCreateGlobalVector(dac, &ac));
429: PetscCall(DADefineXLinearField2D(dac, ac));
431: PetscCall(DMCreateGlobalVector(daf, &af));
432: PetscCall(MatMult(INTERP, ac, af));
434: {
435: Vec afexact;
436: PetscReal nrm;
437: PetscInt N;
439: PetscCall(DMCreateGlobalVector(daf, &afexact));
440: PetscCall(DADefineXLinearField2D(daf, afexact));
441: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
442: PetscCall(VecNorm(afexact, NORM_2, &nrm));
443: PetscCall(VecGetSize(afexact, &N));
444: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, Mx, My, (double)(nrm / PetscSqrtReal((PetscReal)N))));
445: PetscCall(VecDestroy(&afexact));
446: }
448: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
449: if (output) {
450: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
451: PetscCall(VecView(ac, vv));
452: PetscCall(PetscViewerDestroy(&vv));
454: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
455: PetscCall(VecView(af, vv));
456: PetscCall(PetscViewerDestroy(&vv));
457: }
459: PetscCall(MatDestroy(&INTERP));
460: PetscCall(DMDestroy(&dac));
461: PetscCall(DMDestroy(&daf));
462: PetscCall(VecDestroy(&ac));
463: PetscCall(VecDestroy(&af));
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
468: {
469: DM dac, daf;
470: PetscViewer vv;
471: Vec ac, af;
472: PetscInt map_id, Mx, My, Mz;
473: Mat II, INTERP;
474: Vec scale;
475: PetscBool output = PETSC_FALSE;
477: PetscFunctionBeginUser;
478: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, /* 1 dof */
479: 1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
480: PetscCall(DMSetFromOptions(dac));
481: PetscCall(DMSetUp(dac));
483: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
484: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
485: Mx--;
486: My--;
487: Mz--;
489: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
490: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
492: /* apply trilinear mappings */
493: /*PetscCall(DAApplyTrilinearMapping(dac));*/
494: /* apply conformal mappings */
495: map_id = 0;
496: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
497: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
499: {
500: DM cdaf, cdac;
501: Vec coordsc, coordsf;
503: PetscCall(DMGetCoordinateDM(dac, &cdac));
504: PetscCall(DMGetCoordinateDM(daf, &cdaf));
506: PetscCall(DMGetCoordinates(dac, &coordsc));
507: PetscCall(DMGetCoordinates(daf, &coordsf));
509: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
510: PetscCall(MatInterpolate(II, coordsc, coordsf));
511: PetscCall(MatDestroy(&II));
512: PetscCall(VecDestroy(&scale));
513: }
515: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
517: PetscCall(DMCreateGlobalVector(dac, &ac));
518: PetscCall(DADefineXLinearField3D(dac, ac));
520: PetscCall(DMCreateGlobalVector(daf, &af));
522: PetscCall(MatMult(INTERP, ac, af));
524: {
525: Vec afexact;
526: PetscReal nrm;
527: PetscInt N;
529: PetscCall(DMCreateGlobalVector(daf, &afexact));
530: PetscCall(DADefineXLinearField3D(daf, afexact));
531: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
532: PetscCall(VecNorm(afexact, NORM_2, &nrm));
533: PetscCall(VecGetSize(afexact, &N));
534: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, mz, Mx, My, Mz, (double)(nrm / PetscSqrtReal((PetscReal)N))));
535: PetscCall(VecDestroy(&afexact));
536: }
538: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
539: if (output) {
540: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
541: PetscCall(VecView(ac, vv));
542: PetscCall(PetscViewerDestroy(&vv));
544: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
545: PetscCall(VecView(af, vv));
546: PetscCall(PetscViewerDestroy(&vv));
547: }
549: PetscCall(MatDestroy(&INTERP));
550: PetscCall(DMDestroy(&dac));
551: PetscCall(DMDestroy(&daf));
552: PetscCall(VecDestroy(&ac));
553: PetscCall(VecDestroy(&af));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: int main(int argc, char **argv)
558: {
559: PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;
561: PetscFunctionBeginUser;
562: PetscCall(PetscInitialize(&argc, &argv, 0, help));
563: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
564: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
565: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
566: nl = 1;
567: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
568: dim = 2;
569: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));
571: for (l = 0; l < nl; l++) {
572: if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
573: else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
574: else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
575: mx = mx * 2;
576: my = my * 2;
577: mz = mz * 2;
578: }
579: PetscCall(PetscFinalize());
580: return 0;
581: }
583: /*TEST
585: test:
586: suffix: 1d
587: args: -mx 10 -nl 6 -dim 1
589: test:
590: suffix: 2d
591: output_file: output/ex36_2d.out
592: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}
594: test:
595: suffix: 2dp1
596: nsize: 8
597: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
598: timeoutfactor: 2
600: test:
601: suffix: 2dp2
602: nsize: 8
603: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
604: timeoutfactor: 2
606: test:
607: suffix: 3d
608: args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
610: test:
611: suffix: 3dp1
612: nsize: 32
613: args: -mx 5 -my 5 -mz 5 -nl 3 -dim 3 -cmap 1 -da_refine_x 1 -da_refine_y 3 -da_refine_z 4
615: TEST*/