Actual source code: ex36.c
2: static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
7: typedef struct _n_CCmplx CCmplx;
8: struct _n_CCmplx {
9: PetscReal real;
10: PetscReal imag;
11: };
13: CCmplx CCmplxPow(CCmplx a, PetscReal n)
14: {
15: CCmplx b;
16: PetscReal r, theta;
17: r = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
18: theta = PetscAtan2Real(a.imag, a.real);
19: b.real = PetscPowReal(r, n) * PetscCosReal(n * theta);
20: b.imag = PetscPowReal(r, n) * PetscSinReal(n * theta);
21: return b;
22: }
23: CCmplx CCmplxExp(CCmplx a)
24: {
25: CCmplx b;
26: b.real = PetscExpReal(a.real) * PetscCosReal(a.imag);
27: b.imag = PetscExpReal(a.real) * PetscSinReal(a.imag);
28: return b;
29: }
30: CCmplx CCmplxSqrt(CCmplx a)
31: {
32: CCmplx b;
33: PetscReal r, theta;
34: r = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
35: theta = PetscAtan2Real(a.imag, a.real);
36: b.real = PetscSqrtReal(r) * PetscCosReal(0.5 * theta);
37: b.imag = PetscSqrtReal(r) * PetscSinReal(0.5 * theta);
38: return b;
39: }
40: CCmplx CCmplxAdd(CCmplx a, CCmplx c)
41: {
42: CCmplx b;
43: b.real = a.real + c.real;
44: b.imag = a.imag + c.imag;
45: return b;
46: }
47: PetscScalar CCmplxRe(CCmplx a)
48: {
49: return (PetscScalar)a.real;
50: }
51: PetscScalar CCmplxIm(CCmplx a)
52: {
53: return (PetscScalar)a.imag;
54: }
56: PetscErrorCode DAApplyConformalMapping(DM da, PetscInt idx)
57: {
58: PetscInt i, n;
59: PetscInt sx, nx, sy, ny, sz, nz, dim;
60: Vec Gcoords;
61: PetscScalar *XX;
62: PetscScalar xx, yy, zz;
63: DM cda;
65: PetscFunctionBeginUser;
66: if (idx == 0) {
67: PetscFunctionReturn(PETSC_SUCCESS);
68: } else if (idx == 1) { /* dam break */
69: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
70: } else if (idx == 2) { /* stagnation in a corner */
71: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, 0.0, 1.0, -1.0, 1.0));
72: } else if (idx == 3) { /* nautilis */
73: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
74: } else if (idx == 4) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
76: PetscCall(DMGetCoordinateDM(da, &cda));
77: PetscCall(DMGetCoordinates(da, &Gcoords));
79: PetscCall(VecGetArray(Gcoords, &XX));
80: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
81: PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
82: PetscCall(VecGetLocalSize(Gcoords, &n));
83: n = n / dim;
85: for (i = 0; i < n; i++) {
86: if ((dim == 3) && (idx != 2)) {
87: PetscScalar Ni[8];
88: PetscScalar xi = XX[dim * i];
89: PetscScalar eta = XX[dim * i + 1];
90: PetscScalar zeta = XX[dim * i + 2];
91: PetscScalar xn[] = {-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0};
92: PetscScalar yn[] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0};
93: PetscScalar zn[] = {-0.1, -4.0, -0.2, -1.0, 0.1, 4.0, 0.2, 1.0};
94: PetscInt p;
96: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
97: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
98: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
99: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
101: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
102: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
103: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
104: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
106: xx = yy = zz = 0.0;
107: for (p = 0; p < 8; p++) {
108: xx += Ni[p] * xn[p];
109: yy += Ni[p] * yn[p];
110: zz += Ni[p] * zn[p];
111: }
112: XX[dim * i] = xx;
113: XX[dim * i + 1] = yy;
114: XX[dim * i + 2] = zz;
115: }
117: if (idx == 1) {
118: CCmplx zeta, t1, t2;
120: xx = XX[dim * i] - 0.8;
121: yy = XX[dim * i + 1] + 1.5;
123: zeta.real = PetscRealPart(xx);
124: zeta.imag = PetscRealPart(yy);
126: t1 = CCmplxPow(zeta, -1.0);
127: t2 = CCmplxAdd(zeta, t1);
129: XX[dim * i] = CCmplxRe(t2);
130: XX[dim * i + 1] = CCmplxIm(t2);
131: } else if (idx == 2) {
132: CCmplx zeta, t1;
134: xx = XX[dim * i];
135: yy = XX[dim * i + 1];
136: zeta.real = PetscRealPart(xx);
137: zeta.imag = PetscRealPart(yy);
139: t1 = CCmplxSqrt(zeta);
140: XX[dim * i] = CCmplxRe(t1);
141: XX[dim * i + 1] = CCmplxIm(t1);
142: } else if (idx == 3) {
143: CCmplx zeta, t1, t2;
145: xx = XX[dim * i] - 0.8;
146: yy = XX[dim * i + 1] + 1.5;
148: zeta.real = PetscRealPart(xx);
149: zeta.imag = PetscRealPart(yy);
150: t1 = CCmplxPow(zeta, -1.0);
151: t2 = CCmplxAdd(zeta, t1);
152: XX[dim * i] = CCmplxRe(t2);
153: XX[dim * i + 1] = CCmplxIm(t2);
155: xx = XX[dim * i];
156: yy = XX[dim * i + 1];
157: zeta.real = PetscRealPart(xx);
158: zeta.imag = PetscRealPart(yy);
159: t1 = CCmplxExp(zeta);
160: XX[dim * i] = CCmplxRe(t1);
161: XX[dim * i + 1] = CCmplxIm(t1);
163: xx = XX[dim * i] + 0.4;
164: yy = XX[dim * i + 1];
165: zeta.real = PetscRealPart(xx);
166: zeta.imag = PetscRealPart(yy);
167: t1 = CCmplxPow(zeta, 2.0);
168: XX[dim * i] = CCmplxRe(t1);
169: XX[dim * i + 1] = CCmplxIm(t1);
170: } else if (idx == 4) {
171: PetscScalar Ni[4];
172: PetscScalar xi = XX[dim * i];
173: PetscScalar eta = XX[dim * i + 1];
174: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5};
175: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0};
176: PetscInt p;
178: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
179: Ni[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
180: Ni[2] = 0.25 * (1.0 - xi) * (1.0 + eta);
181: Ni[3] = 0.25 * (1.0 + xi) * (1.0 + eta);
183: xx = yy = 0.0;
184: for (p = 0; p < 4; p++) {
185: xx += Ni[p] * xn[p];
186: yy += Ni[p] * yn[p];
187: }
188: XX[dim * i] = xx;
189: XX[dim * i + 1] = yy;
190: }
191: }
192: PetscCall(VecRestoreArray(Gcoords, &XX));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PetscErrorCode DAApplyTrilinearMapping(DM da)
197: {
198: PetscInt i, j, k;
199: PetscInt sx, nx, sy, ny, sz, nz;
200: Vec Gcoords;
201: DMDACoor3d ***XX;
202: PetscScalar xx, yy, zz;
203: DM cda;
205: PetscFunctionBeginUser;
206: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
207: PetscCall(DMGetCoordinateDM(da, &cda));
208: PetscCall(DMGetCoordinates(da, &Gcoords));
210: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
211: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
213: for (i = sx; i < sx + nx; i++) {
214: for (j = sy; j < sy + ny; j++) {
215: for (k = sz; k < sz + nz; k++) {
216: PetscScalar Ni[8];
217: PetscScalar xi = XX[k][j][i].x;
218: PetscScalar eta = XX[k][j][i].y;
219: PetscScalar zeta = XX[k][j][i].z;
220: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5, 0.0, 2.1, 0.23, 3.125};
221: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0, -1.45, -0.1, 2.24, 3.79};
222: PetscScalar zn[] = {0.0, 0.3, -0.1, 0.123, 0.956, 1.32, 1.12, 0.798};
223: PetscInt p;
225: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
226: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
227: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
228: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
230: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
231: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
232: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
233: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
235: xx = yy = zz = 0.0;
236: for (p = 0; p < 8; p++) {
237: xx += Ni[p] * xn[p];
238: yy += Ni[p] * yn[p];
239: zz += Ni[p] * zn[p];
240: }
241: XX[k][j][i].x = xx;
242: XX[k][j][i].y = yy;
243: XX[k][j][i].z = zz;
244: }
245: }
246: }
247: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PetscErrorCode DADefineXLinearField2D(DM da, Vec field)
252: {
253: PetscInt i, j;
254: PetscInt sx, nx, sy, ny;
255: Vec Gcoords;
256: DMDACoor2d **XX;
257: PetscScalar **FF;
258: DM cda;
260: PetscFunctionBeginUser;
261: PetscCall(DMGetCoordinateDM(da, &cda));
262: PetscCall(DMGetCoordinates(da, &Gcoords));
264: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
265: PetscCall(DMDAVecGetArray(da, field, &FF));
267: PetscCall(DMDAGetCorners(da, &sx, &sy, 0, &nx, &ny, 0));
269: for (i = sx; i < sx + nx; i++) {
270: 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;
271: }
273: PetscCall(DMDAVecRestoreArray(da, field, &FF));
274: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: PetscErrorCode DADefineXLinearField3D(DM da, Vec field)
279: {
280: PetscInt i, j, k;
281: PetscInt sx, nx, sy, ny, sz, nz;
282: Vec Gcoords;
283: DMDACoor3d ***XX;
284: PetscScalar ***FF;
285: DM cda;
287: PetscFunctionBeginUser;
288: PetscCall(DMGetCoordinateDM(da, &cda));
289: PetscCall(DMGetCoordinates(da, &Gcoords));
291: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
292: PetscCall(DMDAVecGetArray(da, field, &FF));
294: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
296: for (k = sz; k < sz + nz; k++) {
297: for (j = sy; j < sy + ny; j++) {
298: for (i = sx; i < sx + nx; i++) {
299: 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 +
300: 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
301: }
302: }
303: }
305: PetscCall(DMDAVecRestoreArray(da, field, &FF));
306: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
311: {
312: DM dac, daf;
313: PetscViewer vv;
314: Vec ac, af;
315: PetscInt Mx;
316: Mat II, INTERP;
317: Vec scale;
318: PetscBool output = PETSC_FALSE;
320: PetscFunctionBeginUser;
321: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
322: PetscCall(DMSetFromOptions(dac));
323: PetscCall(DMSetUp(dac));
325: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
326: PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
327: Mx--;
329: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
330: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
332: {
333: DM cdaf, cdac;
334: Vec coordsc, coordsf;
336: PetscCall(DMGetCoordinateDM(dac, &cdac));
337: PetscCall(DMGetCoordinateDM(daf, &cdaf));
339: PetscCall(DMGetCoordinates(dac, &coordsc));
340: PetscCall(DMGetCoordinates(daf, &coordsf));
342: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
343: PetscCall(MatInterpolate(II, coordsc, coordsf));
344: PetscCall(MatDestroy(&II));
345: PetscCall(VecDestroy(&scale));
346: }
348: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
350: PetscCall(DMCreateGlobalVector(dac, &ac));
351: PetscCall(VecSet(ac, 66.99));
353: PetscCall(DMCreateGlobalVector(daf, &af));
355: PetscCall(MatMult(INTERP, ac, af));
357: {
358: Vec afexact;
359: PetscReal nrm;
360: PetscInt N;
362: PetscCall(DMCreateGlobalVector(daf, &afexact));
363: PetscCall(VecSet(afexact, 66.99));
364: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
365: PetscCall(VecNorm(afexact, NORM_2, &nrm));
366: PetscCall(VecGetSize(afexact, &N));
367: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
368: PetscCall(VecDestroy(&afexact));
369: }
371: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
372: if (output) {
373: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
374: PetscCall(VecView(ac, vv));
375: PetscCall(PetscViewerDestroy(&vv));
377: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
378: PetscCall(VecView(af, vv));
379: PetscCall(PetscViewerDestroy(&vv));
380: }
382: PetscCall(MatDestroy(&INTERP));
383: PetscCall(DMDestroy(&dac));
384: PetscCall(DMDestroy(&daf));
385: PetscCall(VecDestroy(&ac));
386: PetscCall(VecDestroy(&af));
387: PetscFunctionReturn(PETSC_SUCCESS);
388: }
390: PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
391: {
392: DM dac, daf;
393: PetscViewer vv;
394: Vec ac, af;
395: PetscInt map_id, Mx, My;
396: Mat II, INTERP;
397: Vec scale;
398: PetscBool output = PETSC_FALSE;
400: PetscFunctionBeginUser;
401: 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));
402: PetscCall(DMSetFromOptions(dac));
403: PetscCall(DMSetUp(dac));
405: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
406: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
407: Mx--;
408: My--;
410: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
411: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
413: /* apply conformal mappings */
414: map_id = 0;
415: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
416: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
418: {
419: DM cdaf, cdac;
420: Vec coordsc, coordsf;
422: PetscCall(DMGetCoordinateDM(dac, &cdac));
423: PetscCall(DMGetCoordinateDM(daf, &cdaf));
425: PetscCall(DMGetCoordinates(dac, &coordsc));
426: PetscCall(DMGetCoordinates(daf, &coordsf));
428: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
429: PetscCall(MatInterpolate(II, coordsc, coordsf));
430: PetscCall(MatDestroy(&II));
431: PetscCall(VecDestroy(&scale));
432: }
434: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
436: PetscCall(DMCreateGlobalVector(dac, &ac));
437: PetscCall(DADefineXLinearField2D(dac, ac));
439: PetscCall(DMCreateGlobalVector(daf, &af));
440: PetscCall(MatMult(INTERP, ac, af));
442: {
443: Vec afexact;
444: PetscReal nrm;
445: PetscInt N;
447: PetscCall(DMCreateGlobalVector(daf, &afexact));
448: PetscCall(VecZeroEntries(afexact));
449: PetscCall(DADefineXLinearField2D(daf, afexact));
450: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
451: PetscCall(VecNorm(afexact, NORM_2, &nrm));
452: PetscCall(VecGetSize(afexact, &N));
453: 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))));
454: PetscCall(VecDestroy(&afexact));
455: }
457: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
458: if (output) {
459: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
460: PetscCall(VecView(ac, vv));
461: PetscCall(PetscViewerDestroy(&vv));
463: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
464: PetscCall(VecView(af, vv));
465: PetscCall(PetscViewerDestroy(&vv));
466: }
468: PetscCall(MatDestroy(&INTERP));
469: PetscCall(DMDestroy(&dac));
470: PetscCall(DMDestroy(&daf));
471: PetscCall(VecDestroy(&ac));
472: PetscCall(VecDestroy(&af));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
477: {
478: DM dac, daf;
479: PetscViewer vv;
480: Vec ac, af;
481: PetscInt map_id, Mx, My, Mz;
482: Mat II, INTERP;
483: Vec scale;
484: PetscBool output = PETSC_FALSE;
486: PetscFunctionBeginUser;
487: 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 */
488: 1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
489: PetscCall(DMSetFromOptions(dac));
490: PetscCall(DMSetUp(dac));
492: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
493: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
494: Mx--;
495: My--;
496: Mz--;
498: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
499: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
501: /* apply trilinear mappings */
502: /*PetscCall(DAApplyTrilinearMapping(dac));*/
503: /* apply conformal mappings */
504: map_id = 0;
505: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
506: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
508: {
509: DM cdaf, cdac;
510: Vec coordsc, coordsf;
512: PetscCall(DMGetCoordinateDM(dac, &cdac));
513: PetscCall(DMGetCoordinateDM(daf, &cdaf));
515: PetscCall(DMGetCoordinates(dac, &coordsc));
516: PetscCall(DMGetCoordinates(daf, &coordsf));
518: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
519: PetscCall(MatInterpolate(II, coordsc, coordsf));
520: PetscCall(MatDestroy(&II));
521: PetscCall(VecDestroy(&scale));
522: }
524: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
526: PetscCall(DMCreateGlobalVector(dac, &ac));
527: PetscCall(VecZeroEntries(ac));
528: PetscCall(DADefineXLinearField3D(dac, ac));
530: PetscCall(DMCreateGlobalVector(daf, &af));
531: PetscCall(VecZeroEntries(af));
533: PetscCall(MatMult(INTERP, ac, af));
535: {
536: Vec afexact;
537: PetscReal nrm;
538: PetscInt N;
540: PetscCall(DMCreateGlobalVector(daf, &afexact));
541: PetscCall(VecZeroEntries(afexact));
542: PetscCall(DADefineXLinearField3D(daf, afexact));
543: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
544: PetscCall(VecNorm(afexact, NORM_2, &nrm));
545: PetscCall(VecGetSize(afexact, &N));
546: 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))));
547: PetscCall(VecDestroy(&afexact));
548: }
550: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
551: if (output) {
552: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
553: PetscCall(VecView(ac, vv));
554: PetscCall(PetscViewerDestroy(&vv));
556: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
557: PetscCall(VecView(af, vv));
558: PetscCall(PetscViewerDestroy(&vv));
559: }
561: PetscCall(MatDestroy(&INTERP));
562: PetscCall(DMDestroy(&dac));
563: PetscCall(DMDestroy(&daf));
564: PetscCall(VecDestroy(&ac));
565: PetscCall(VecDestroy(&af));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: int main(int argc, char **argv)
570: {
571: PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;
573: PetscFunctionBeginUser;
574: PetscCall(PetscInitialize(&argc, &argv, 0, help));
575: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
576: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
577: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
578: nl = 1;
579: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
580: dim = 2;
581: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));
583: for (l = 0; l < nl; l++) {
584: if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
585: else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
586: else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
587: mx = mx * 2;
588: my = my * 2;
589: mz = mz * 2;
590: }
591: PetscCall(PetscFinalize());
592: return 0;
593: }
595: /*TEST
597: test:
598: suffix: 1d
599: args: -mx 10 -nl 6 -dim 1
601: test:
602: suffix: 2d
603: output_file: output/ex36_2d.out
604: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}
606: test:
607: suffix: 2dp1
608: nsize: 8
609: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
610: timeoutfactor: 2
612: test:
613: suffix: 2dp2
614: nsize: 8
615: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
616: timeoutfactor: 2
618: test:
619: suffix: 3d
620: args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
622: test:
623: suffix: 3dp1
624: nsize: 32
625: 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
627: TEST*/