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};
89: PetscInt p;
91: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
92: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
93: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
94: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
96: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
97: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
98: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
99: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
101: xx = yy = zz = 0.0;
102: for (p = 0; p < 8; p++) {
103: xx += Ni[p] * xn[p];
104: yy += Ni[p] * yn[p];
105: zz += Ni[p] * zn[p];
106: }
107: XX[dim * i] = xx;
108: XX[dim * i + 1] = yy;
109: XX[dim * i + 2] = zz;
110: }
112: if (idx == 1) {
113: CCmplx zeta, t1, t2;
115: xx = XX[dim * i] - 0.8;
116: yy = XX[dim * i + 1] + 1.5;
118: zeta.real = PetscRealPart(xx);
119: zeta.imag = PetscRealPart(yy);
121: t1 = CCmplxPow(zeta, -1.0);
122: t2 = CCmplxAdd(zeta, t1);
124: XX[dim * i] = CCmplxRe(t2);
125: XX[dim * i + 1] = CCmplxIm(t2);
126: } else if (idx == 2) {
127: CCmplx zeta, t1;
129: xx = XX[dim * i];
130: yy = XX[dim * i + 1];
131: zeta.real = PetscRealPart(xx);
132: zeta.imag = PetscRealPart(yy);
134: t1 = CCmplxSqrt(zeta);
135: XX[dim * i] = CCmplxRe(t1);
136: XX[dim * i + 1] = CCmplxIm(t1);
137: } else if (idx == 3) {
138: CCmplx zeta, t1, t2;
140: xx = XX[dim * i] - 0.8;
141: yy = XX[dim * i + 1] + 1.5;
143: zeta.real = PetscRealPart(xx);
144: zeta.imag = PetscRealPart(yy);
145: t1 = CCmplxPow(zeta, -1.0);
146: t2 = CCmplxAdd(zeta, t1);
147: XX[dim * i] = CCmplxRe(t2);
148: XX[dim * i + 1] = CCmplxIm(t2);
150: xx = XX[dim * i];
151: yy = XX[dim * i + 1];
152: zeta.real = PetscRealPart(xx);
153: zeta.imag = PetscRealPart(yy);
154: t1 = CCmplxExp(zeta);
155: XX[dim * i] = CCmplxRe(t1);
156: XX[dim * i + 1] = CCmplxIm(t1);
158: xx = XX[dim * i] + 0.4;
159: yy = XX[dim * i + 1];
160: zeta.real = PetscRealPart(xx);
161: zeta.imag = PetscRealPart(yy);
162: t1 = CCmplxPow(zeta, 2.0);
163: XX[dim * i] = CCmplxRe(t1);
164: XX[dim * i + 1] = CCmplxIm(t1);
165: } else if (idx == 4) {
166: PetscScalar Ni[4];
167: PetscScalar xi = XX[dim * i];
168: PetscScalar eta = XX[dim * i + 1];
169: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5};
170: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0};
171: PetscInt p;
173: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
174: Ni[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
175: Ni[2] = 0.25 * (1.0 - xi) * (1.0 + eta);
176: Ni[3] = 0.25 * (1.0 + xi) * (1.0 + eta);
178: xx = yy = 0.0;
179: for (p = 0; p < 4; p++) {
180: xx += Ni[p] * xn[p];
181: yy += Ni[p] * yn[p];
182: }
183: XX[dim * i] = xx;
184: XX[dim * i + 1] = yy;
185: }
186: }
187: PetscCall(VecRestoreArray(Gcoords, &XX));
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode DAApplyTrilinearMapping(DM da)
192: {
193: PetscInt i, j, k;
194: PetscInt sx, nx, sy, ny, sz, nz;
195: Vec Gcoords;
196: DMDACoor3d ***XX;
197: PetscScalar xx, yy, zz;
198: DM cda;
200: PetscFunctionBeginUser;
201: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
202: PetscCall(DMGetCoordinateDM(da, &cda));
203: PetscCall(DMGetCoordinates(da, &Gcoords));
205: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
206: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
208: for (i = sx; i < sx + nx; i++) {
209: for (j = sy; j < sy + ny; j++) {
210: for (k = sz; k < sz + nz; k++) {
211: PetscScalar Ni[8];
212: PetscScalar xi = XX[k][j][i].x;
213: PetscScalar eta = XX[k][j][i].y;
214: PetscScalar zeta = XX[k][j][i].z;
215: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5, 0.0, 2.1, 0.23, 3.125};
216: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0, -1.45, -0.1, 2.24, 3.79};
217: PetscScalar zn[] = {0.0, 0.3, -0.1, 0.123, 0.956, 1.32, 1.12, 0.798};
218: PetscInt p;
220: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
221: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
222: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
223: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
225: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
226: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
227: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
228: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
230: xx = yy = zz = 0.0;
231: for (p = 0; p < 8; p++) {
232: xx += Ni[p] * xn[p];
233: yy += Ni[p] * yn[p];
234: zz += Ni[p] * zn[p];
235: }
236: XX[k][j][i].x = xx;
237: XX[k][j][i].y = yy;
238: XX[k][j][i].z = zz;
239: }
240: }
241: }
242: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: PetscErrorCode DADefineXLinearField2D(DM da, Vec field)
247: {
248: PetscInt i, j;
249: PetscInt sx, nx, sy, ny;
250: Vec Gcoords;
251: DMDACoor2d **XX;
252: PetscScalar **FF;
253: DM cda;
255: PetscFunctionBeginUser;
256: PetscCall(DMGetCoordinateDM(da, &cda));
257: PetscCall(DMGetCoordinates(da, &Gcoords));
259: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
260: PetscCall(DMDAVecGetArray(da, field, &FF));
262: PetscCall(DMDAGetCorners(da, &sx, &sy, 0, &nx, &ny, 0));
264: for (i = sx; i < sx + nx; i++) {
265: 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;
266: }
268: PetscCall(DMDAVecRestoreArray(da, field, &FF));
269: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: PetscErrorCode DADefineXLinearField3D(DM da, Vec field)
274: {
275: PetscInt i, j, k;
276: PetscInt sx, nx, sy, ny, sz, nz;
277: Vec Gcoords;
278: DMDACoor3d ***XX;
279: PetscScalar ***FF;
280: DM cda;
282: PetscFunctionBeginUser;
283: PetscCall(DMGetCoordinateDM(da, &cda));
284: PetscCall(DMGetCoordinates(da, &Gcoords));
286: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
287: PetscCall(DMDAVecGetArray(da, field, &FF));
289: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
291: for (k = sz; k < sz + nz; k++) {
292: for (j = sy; j < sy + ny; j++) {
293: for (i = sx; i < sx + nx; i++) {
294: 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 +
295: 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
296: }
297: }
298: }
300: PetscCall(DMDAVecRestoreArray(da, field, &FF));
301: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
306: {
307: DM dac, daf;
308: PetscViewer vv;
309: Vec ac, af;
310: PetscInt Mx;
311: Mat II, INTERP;
312: Vec scale;
313: PetscBool output = PETSC_FALSE;
315: PetscFunctionBeginUser;
316: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
317: PetscCall(DMSetFromOptions(dac));
318: PetscCall(DMSetUp(dac));
320: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
321: PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
322: Mx--;
324: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
325: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
327: {
328: DM cdaf, cdac;
329: Vec coordsc, coordsf;
331: PetscCall(DMGetCoordinateDM(dac, &cdac));
332: PetscCall(DMGetCoordinateDM(daf, &cdaf));
334: PetscCall(DMGetCoordinates(dac, &coordsc));
335: PetscCall(DMGetCoordinates(daf, &coordsf));
337: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
338: PetscCall(MatInterpolate(II, coordsc, coordsf));
339: PetscCall(MatDestroy(&II));
340: PetscCall(VecDestroy(&scale));
341: }
343: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
345: PetscCall(DMCreateGlobalVector(dac, &ac));
346: PetscCall(VecSet(ac, 66.99));
348: PetscCall(DMCreateGlobalVector(daf, &af));
350: PetscCall(MatMult(INTERP, ac, af));
352: {
353: Vec afexact;
354: PetscReal nrm;
355: PetscInt N;
357: PetscCall(DMCreateGlobalVector(daf, &afexact));
358: PetscCall(VecSet(afexact, 66.99));
359: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
360: PetscCall(VecNorm(afexact, NORM_2, &nrm));
361: PetscCall(VecGetSize(afexact, &N));
362: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
363: PetscCall(VecDestroy(&afexact));
364: }
366: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
367: if (output) {
368: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
369: PetscCall(VecView(ac, vv));
370: PetscCall(PetscViewerDestroy(&vv));
372: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
373: PetscCall(VecView(af, vv));
374: PetscCall(PetscViewerDestroy(&vv));
375: }
377: PetscCall(MatDestroy(&INTERP));
378: PetscCall(DMDestroy(&dac));
379: PetscCall(DMDestroy(&daf));
380: PetscCall(VecDestroy(&ac));
381: PetscCall(VecDestroy(&af));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
386: {
387: DM dac, daf;
388: PetscViewer vv;
389: Vec ac, af;
390: PetscInt map_id, Mx, My;
391: Mat II, INTERP;
392: Vec scale;
393: PetscBool output = PETSC_FALSE;
395: PetscFunctionBeginUser;
396: 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));
397: PetscCall(DMSetFromOptions(dac));
398: PetscCall(DMSetUp(dac));
400: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
401: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
402: Mx--;
403: My--;
405: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
406: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
408: /* apply conformal mappings */
409: map_id = 0;
410: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
411: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
413: {
414: DM cdaf, cdac;
415: Vec coordsc, coordsf;
417: PetscCall(DMGetCoordinateDM(dac, &cdac));
418: PetscCall(DMGetCoordinateDM(daf, &cdaf));
420: PetscCall(DMGetCoordinates(dac, &coordsc));
421: PetscCall(DMGetCoordinates(daf, &coordsf));
423: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
424: PetscCall(MatInterpolate(II, coordsc, coordsf));
425: PetscCall(MatDestroy(&II));
426: PetscCall(VecDestroy(&scale));
427: }
429: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
431: PetscCall(DMCreateGlobalVector(dac, &ac));
432: PetscCall(DADefineXLinearField2D(dac, ac));
434: PetscCall(DMCreateGlobalVector(daf, &af));
435: PetscCall(MatMult(INTERP, ac, af));
437: {
438: Vec afexact;
439: PetscReal nrm;
440: PetscInt N;
442: PetscCall(DMCreateGlobalVector(daf, &afexact));
443: PetscCall(VecZeroEntries(afexact));
444: PetscCall(DADefineXLinearField2D(daf, afexact));
445: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
446: PetscCall(VecNorm(afexact, NORM_2, &nrm));
447: PetscCall(VecGetSize(afexact, &N));
448: 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))));
449: PetscCall(VecDestroy(&afexact));
450: }
452: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
453: if (output) {
454: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
455: PetscCall(VecView(ac, vv));
456: PetscCall(PetscViewerDestroy(&vv));
458: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
459: PetscCall(VecView(af, vv));
460: PetscCall(PetscViewerDestroy(&vv));
461: }
463: PetscCall(MatDestroy(&INTERP));
464: PetscCall(DMDestroy(&dac));
465: PetscCall(DMDestroy(&daf));
466: PetscCall(VecDestroy(&ac));
467: PetscCall(VecDestroy(&af));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
472: {
473: DM dac, daf;
474: PetscViewer vv;
475: Vec ac, af;
476: PetscInt map_id, Mx, My, Mz;
477: Mat II, INTERP;
478: Vec scale;
479: PetscBool output = PETSC_FALSE;
481: PetscFunctionBeginUser;
482: 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 */
483: 1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
484: PetscCall(DMSetFromOptions(dac));
485: PetscCall(DMSetUp(dac));
487: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
488: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
489: Mx--;
490: My--;
491: Mz--;
493: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
494: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
496: /* apply trilinear mappings */
497: /*PetscCall(DAApplyTrilinearMapping(dac));*/
498: /* apply conformal mappings */
499: map_id = 0;
500: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
501: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
503: {
504: DM cdaf, cdac;
505: Vec coordsc, coordsf;
507: PetscCall(DMGetCoordinateDM(dac, &cdac));
508: PetscCall(DMGetCoordinateDM(daf, &cdaf));
510: PetscCall(DMGetCoordinates(dac, &coordsc));
511: PetscCall(DMGetCoordinates(daf, &coordsf));
513: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
514: PetscCall(MatInterpolate(II, coordsc, coordsf));
515: PetscCall(MatDestroy(&II));
516: PetscCall(VecDestroy(&scale));
517: }
519: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
521: PetscCall(DMCreateGlobalVector(dac, &ac));
522: PetscCall(VecZeroEntries(ac));
523: PetscCall(DADefineXLinearField3D(dac, ac));
525: PetscCall(DMCreateGlobalVector(daf, &af));
526: PetscCall(VecZeroEntries(af));
528: PetscCall(MatMult(INTERP, ac, af));
530: {
531: Vec afexact;
532: PetscReal nrm;
533: PetscInt N;
535: PetscCall(DMCreateGlobalVector(daf, &afexact));
536: PetscCall(VecZeroEntries(afexact));
537: PetscCall(DADefineXLinearField3D(daf, afexact));
538: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
539: PetscCall(VecNorm(afexact, NORM_2, &nrm));
540: PetscCall(VecGetSize(afexact, &N));
541: 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))));
542: PetscCall(VecDestroy(&afexact));
543: }
545: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
546: if (output) {
547: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
548: PetscCall(VecView(ac, vv));
549: PetscCall(PetscViewerDestroy(&vv));
551: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
552: PetscCall(VecView(af, vv));
553: PetscCall(PetscViewerDestroy(&vv));
554: }
556: PetscCall(MatDestroy(&INTERP));
557: PetscCall(DMDestroy(&dac));
558: PetscCall(DMDestroy(&daf));
559: PetscCall(VecDestroy(&ac));
560: PetscCall(VecDestroy(&af));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: int main(int argc, char **argv)
565: {
566: PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;
568: PetscFunctionBeginUser;
569: PetscCall(PetscInitialize(&argc, &argv, 0, help));
570: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
571: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
572: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
573: nl = 1;
574: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
575: dim = 2;
576: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));
578: for (l = 0; l < nl; l++) {
579: if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
580: else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
581: else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
582: mx = mx * 2;
583: my = my * 2;
584: mz = mz * 2;
585: }
586: PetscCall(PetscFinalize());
587: return 0;
588: }
590: /*TEST
592: test:
593: suffix: 1d
594: args: -mx 10 -nl 6 -dim 1
596: test:
597: suffix: 2d
598: output_file: output/ex36_2d.out
599: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}
601: test:
602: suffix: 2dp1
603: nsize: 8
604: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
605: timeoutfactor: 2
607: test:
608: suffix: 2dp2
609: nsize: 8
610: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
611: timeoutfactor: 2
613: test:
614: suffix: 3d
615: args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
617: test:
618: suffix: 3dp1
619: nsize: 32
620: 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
622: TEST*/