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 (PetscScalar)a.real;
49: }
50: PetscScalar CCmplxIm(CCmplx a)
51: {
52: return (PetscScalar)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) {
66: PetscFunctionReturn(PETSC_SUCCESS);
67: } else if (idx == 1) { /* dam break */
68: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
69: } else if (idx == 2) { /* stagnation in a corner */
70: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, 0.0, 1.0, -1.0, 1.0));
71: } else if (idx == 3) { /* nautilis */
72: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
73: } else if (idx == 4) PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
75: PetscCall(DMGetCoordinateDM(da, &cda));
76: PetscCall(DMGetCoordinates(da, &Gcoords));
78: PetscCall(VecGetArray(Gcoords, &XX));
79: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
80: PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
81: PetscCall(VecGetLocalSize(Gcoords, &n));
82: n = n / dim;
84: for (i = 0; i < n; i++) {
85: if ((dim == 3) && (idx != 2)) {
86: PetscScalar Ni[8];
87: PetscScalar xi = XX[dim * i];
88: PetscScalar eta = XX[dim * i + 1];
89: PetscScalar zeta = XX[dim * i + 2];
90: PetscScalar xn[] = {-1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0};
91: PetscScalar yn[] = {-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0};
92: PetscScalar zn[] = {-0.1, -4.0, -0.2, -1.0, 0.1, 4.0, 0.2, 1.0};
93: PetscInt p;
95: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
96: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
97: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
98: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
100: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
101: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
102: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
103: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
105: xx = yy = zz = 0.0;
106: for (p = 0; p < 8; p++) {
107: xx += Ni[p] * xn[p];
108: yy += Ni[p] * yn[p];
109: zz += Ni[p] * zn[p];
110: }
111: XX[dim * i] = xx;
112: XX[dim * i + 1] = yy;
113: XX[dim * i + 2] = zz;
114: }
116: if (idx == 1) {
117: CCmplx zeta, t1, t2;
119: xx = XX[dim * i] - 0.8;
120: yy = XX[dim * i + 1] + 1.5;
122: zeta.real = PetscRealPart(xx);
123: zeta.imag = PetscRealPart(yy);
125: t1 = CCmplxPow(zeta, -1.0);
126: t2 = CCmplxAdd(zeta, t1);
128: XX[dim * i] = CCmplxRe(t2);
129: XX[dim * i + 1] = CCmplxIm(t2);
130: } else if (idx == 2) {
131: CCmplx zeta, t1;
133: xx = XX[dim * i];
134: yy = XX[dim * i + 1];
135: zeta.real = PetscRealPart(xx);
136: zeta.imag = PetscRealPart(yy);
138: t1 = CCmplxSqrt(zeta);
139: XX[dim * i] = CCmplxRe(t1);
140: XX[dim * i + 1] = CCmplxIm(t1);
141: } else if (idx == 3) {
142: CCmplx zeta, t1, t2;
144: xx = XX[dim * i] - 0.8;
145: yy = XX[dim * i + 1] + 1.5;
147: zeta.real = PetscRealPart(xx);
148: zeta.imag = PetscRealPart(yy);
149: t1 = CCmplxPow(zeta, -1.0);
150: t2 = CCmplxAdd(zeta, t1);
151: XX[dim * i] = CCmplxRe(t2);
152: XX[dim * i + 1] = CCmplxIm(t2);
154: xx = XX[dim * i];
155: yy = XX[dim * i + 1];
156: zeta.real = PetscRealPart(xx);
157: zeta.imag = PetscRealPart(yy);
158: t1 = CCmplxExp(zeta);
159: XX[dim * i] = CCmplxRe(t1);
160: XX[dim * i + 1] = CCmplxIm(t1);
162: xx = XX[dim * i] + 0.4;
163: yy = XX[dim * i + 1];
164: zeta.real = PetscRealPart(xx);
165: zeta.imag = PetscRealPart(yy);
166: t1 = CCmplxPow(zeta, 2.0);
167: XX[dim * i] = CCmplxRe(t1);
168: XX[dim * i + 1] = CCmplxIm(t1);
169: } else if (idx == 4) {
170: PetscScalar Ni[4];
171: PetscScalar xi = XX[dim * i];
172: PetscScalar eta = XX[dim * i + 1];
173: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5};
174: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0};
175: PetscInt p;
177: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
178: Ni[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
179: Ni[2] = 0.25 * (1.0 - xi) * (1.0 + eta);
180: Ni[3] = 0.25 * (1.0 + xi) * (1.0 + eta);
182: xx = yy = 0.0;
183: for (p = 0; p < 4; p++) {
184: xx += Ni[p] * xn[p];
185: yy += Ni[p] * yn[p];
186: }
187: XX[dim * i] = xx;
188: XX[dim * i + 1] = yy;
189: }
190: }
191: PetscCall(VecRestoreArray(Gcoords, &XX));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PetscErrorCode DAApplyTrilinearMapping(DM da)
196: {
197: PetscInt i, j, k;
198: PetscInt sx, nx, sy, ny, sz, nz;
199: Vec Gcoords;
200: DMDACoor3d ***XX;
201: PetscScalar xx, yy, zz;
202: DM cda;
204: PetscFunctionBeginUser;
205: PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
206: PetscCall(DMGetCoordinateDM(da, &cda));
207: PetscCall(DMGetCoordinates(da, &Gcoords));
209: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
210: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
212: for (i = sx; i < sx + nx; i++) {
213: for (j = sy; j < sy + ny; j++) {
214: for (k = sz; k < sz + nz; k++) {
215: PetscScalar Ni[8];
216: PetscScalar xi = XX[k][j][i].x;
217: PetscScalar eta = XX[k][j][i].y;
218: PetscScalar zeta = XX[k][j][i].z;
219: PetscScalar xn[] = {0.0, 2.0, 0.2, 3.5, 0.0, 2.1, 0.23, 3.125};
220: PetscScalar yn[] = {-1.3, 0.0, 2.0, 4.0, -1.45, -0.1, 2.24, 3.79};
221: PetscScalar zn[] = {0.0, 0.3, -0.1, 0.123, 0.956, 1.32, 1.12, 0.798};
222: PetscInt p;
224: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
225: Ni[1] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
226: Ni[2] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
227: Ni[3] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
229: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
230: Ni[5] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
231: Ni[6] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
232: Ni[7] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
234: xx = yy = zz = 0.0;
235: for (p = 0; p < 8; p++) {
236: xx += Ni[p] * xn[p];
237: yy += Ni[p] * yn[p];
238: zz += Ni[p] * zn[p];
239: }
240: XX[k][j][i].x = xx;
241: XX[k][j][i].y = yy;
242: XX[k][j][i].z = zz;
243: }
244: }
245: }
246: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: PetscErrorCode DADefineXLinearField2D(DM da, Vec field)
251: {
252: PetscInt i, j;
253: PetscInt sx, nx, sy, ny;
254: Vec Gcoords;
255: DMDACoor2d **XX;
256: PetscScalar **FF;
257: DM cda;
259: PetscFunctionBeginUser;
260: PetscCall(DMGetCoordinateDM(da, &cda));
261: PetscCall(DMGetCoordinates(da, &Gcoords));
263: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
264: PetscCall(DMDAVecGetArray(da, field, &FF));
266: PetscCall(DMDAGetCorners(da, &sx, &sy, 0, &nx, &ny, 0));
268: for (i = sx; i < sx + nx; i++) {
269: 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;
270: }
272: PetscCall(DMDAVecRestoreArray(da, field, &FF));
273: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode DADefineXLinearField3D(DM da, Vec field)
278: {
279: PetscInt i, j, k;
280: PetscInt sx, nx, sy, ny, sz, nz;
281: Vec Gcoords;
282: DMDACoor3d ***XX;
283: PetscScalar ***FF;
284: DM cda;
286: PetscFunctionBeginUser;
287: PetscCall(DMGetCoordinateDM(da, &cda));
288: PetscCall(DMGetCoordinates(da, &Gcoords));
290: PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
291: PetscCall(DMDAVecGetArray(da, field, &FF));
293: PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));
295: for (k = sz; k < sz + nz; k++) {
296: for (j = sy; j < sy + ny; j++) {
297: for (i = sx; i < sx + nx; i++) {
298: 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 +
299: 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
300: }
301: }
302: }
304: PetscCall(DMDAVecRestoreArray(da, field, &FF));
305: PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
310: {
311: DM dac, daf;
312: PetscViewer vv;
313: Vec ac, af;
314: PetscInt Mx;
315: Mat II, INTERP;
316: Vec scale;
317: PetscBool output = PETSC_FALSE;
319: PetscFunctionBeginUser;
320: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
321: PetscCall(DMSetFromOptions(dac));
322: PetscCall(DMSetUp(dac));
324: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
325: PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
326: Mx--;
328: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
329: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
331: {
332: DM cdaf, cdac;
333: Vec coordsc, coordsf;
335: PetscCall(DMGetCoordinateDM(dac, &cdac));
336: PetscCall(DMGetCoordinateDM(daf, &cdaf));
338: PetscCall(DMGetCoordinates(dac, &coordsc));
339: PetscCall(DMGetCoordinates(daf, &coordsf));
341: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
342: PetscCall(MatInterpolate(II, coordsc, coordsf));
343: PetscCall(MatDestroy(&II));
344: PetscCall(VecDestroy(&scale));
345: }
347: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
349: PetscCall(DMCreateGlobalVector(dac, &ac));
350: PetscCall(VecSet(ac, 66.99));
352: PetscCall(DMCreateGlobalVector(daf, &af));
354: PetscCall(MatMult(INTERP, ac, af));
356: {
357: Vec afexact;
358: PetscReal nrm;
359: PetscInt N;
361: PetscCall(DMCreateGlobalVector(daf, &afexact));
362: PetscCall(VecSet(afexact, 66.99));
363: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
364: PetscCall(VecNorm(afexact, NORM_2, &nrm));
365: PetscCall(VecGetSize(afexact, &N));
366: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
367: PetscCall(VecDestroy(&afexact));
368: }
370: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
371: if (output) {
372: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
373: PetscCall(VecView(ac, vv));
374: PetscCall(PetscViewerDestroy(&vv));
376: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
377: PetscCall(VecView(af, vv));
378: PetscCall(PetscViewerDestroy(&vv));
379: }
381: PetscCall(MatDestroy(&INTERP));
382: PetscCall(DMDestroy(&dac));
383: PetscCall(DMDestroy(&daf));
384: PetscCall(VecDestroy(&ac));
385: PetscCall(VecDestroy(&af));
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
390: {
391: DM dac, daf;
392: PetscViewer vv;
393: Vec ac, af;
394: PetscInt map_id, Mx, My;
395: Mat II, INTERP;
396: Vec scale;
397: PetscBool output = PETSC_FALSE;
399: PetscFunctionBeginUser;
400: 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));
401: PetscCall(DMSetFromOptions(dac));
402: PetscCall(DMSetUp(dac));
404: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
405: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
406: Mx--;
407: My--;
409: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
410: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
412: /* apply conformal mappings */
413: map_id = 0;
414: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
415: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
417: {
418: DM cdaf, cdac;
419: Vec coordsc, coordsf;
421: PetscCall(DMGetCoordinateDM(dac, &cdac));
422: PetscCall(DMGetCoordinateDM(daf, &cdaf));
424: PetscCall(DMGetCoordinates(dac, &coordsc));
425: PetscCall(DMGetCoordinates(daf, &coordsf));
427: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
428: PetscCall(MatInterpolate(II, coordsc, coordsf));
429: PetscCall(MatDestroy(&II));
430: PetscCall(VecDestroy(&scale));
431: }
433: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
435: PetscCall(DMCreateGlobalVector(dac, &ac));
436: PetscCall(DADefineXLinearField2D(dac, ac));
438: PetscCall(DMCreateGlobalVector(daf, &af));
439: PetscCall(MatMult(INTERP, ac, af));
441: {
442: Vec afexact;
443: PetscReal nrm;
444: PetscInt N;
446: PetscCall(DMCreateGlobalVector(daf, &afexact));
447: PetscCall(VecZeroEntries(afexact));
448: PetscCall(DADefineXLinearField2D(daf, afexact));
449: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
450: PetscCall(VecNorm(afexact, NORM_2, &nrm));
451: PetscCall(VecGetSize(afexact, &N));
452: 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))));
453: PetscCall(VecDestroy(&afexact));
454: }
456: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
457: if (output) {
458: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
459: PetscCall(VecView(ac, vv));
460: PetscCall(PetscViewerDestroy(&vv));
462: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
463: PetscCall(VecView(af, vv));
464: PetscCall(PetscViewerDestroy(&vv));
465: }
467: PetscCall(MatDestroy(&INTERP));
468: PetscCall(DMDestroy(&dac));
469: PetscCall(DMDestroy(&daf));
470: PetscCall(VecDestroy(&ac));
471: PetscCall(VecDestroy(&af));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
476: {
477: DM dac, daf;
478: PetscViewer vv;
479: Vec ac, af;
480: PetscInt map_id, Mx, My, Mz;
481: Mat II, INTERP;
482: Vec scale;
483: PetscBool output = PETSC_FALSE;
485: PetscFunctionBeginUser;
486: 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 */
487: 1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
488: PetscCall(DMSetFromOptions(dac));
489: PetscCall(DMSetUp(dac));
491: PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
492: PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
493: Mx--;
494: My--;
495: Mz--;
497: PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
498: PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
500: /* apply trilinear mappings */
501: /*PetscCall(DAApplyTrilinearMapping(dac));*/
502: /* apply conformal mappings */
503: map_id = 0;
504: PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
505: if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));
507: {
508: DM cdaf, cdac;
509: Vec coordsc, coordsf;
511: PetscCall(DMGetCoordinateDM(dac, &cdac));
512: PetscCall(DMGetCoordinateDM(daf, &cdaf));
514: PetscCall(DMGetCoordinates(dac, &coordsc));
515: PetscCall(DMGetCoordinates(daf, &coordsf));
517: PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
518: PetscCall(MatInterpolate(II, coordsc, coordsf));
519: PetscCall(MatDestroy(&II));
520: PetscCall(VecDestroy(&scale));
521: }
523: PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));
525: PetscCall(DMCreateGlobalVector(dac, &ac));
526: PetscCall(VecZeroEntries(ac));
527: PetscCall(DADefineXLinearField3D(dac, ac));
529: PetscCall(DMCreateGlobalVector(daf, &af));
530: PetscCall(VecZeroEntries(af));
532: PetscCall(MatMult(INTERP, ac, af));
534: {
535: Vec afexact;
536: PetscReal nrm;
537: PetscInt N;
539: PetscCall(DMCreateGlobalVector(daf, &afexact));
540: PetscCall(VecZeroEntries(afexact));
541: PetscCall(DADefineXLinearField3D(daf, afexact));
542: PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
543: PetscCall(VecNorm(afexact, NORM_2, &nrm));
544: PetscCall(VecGetSize(afexact, &N));
545: 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))));
546: PetscCall(VecDestroy(&afexact));
547: }
549: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
550: if (output) {
551: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
552: PetscCall(VecView(ac, vv));
553: PetscCall(PetscViewerDestroy(&vv));
555: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
556: PetscCall(VecView(af, vv));
557: PetscCall(PetscViewerDestroy(&vv));
558: }
560: PetscCall(MatDestroy(&INTERP));
561: PetscCall(DMDestroy(&dac));
562: PetscCall(DMDestroy(&daf));
563: PetscCall(VecDestroy(&ac));
564: PetscCall(VecDestroy(&af));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: int main(int argc, char **argv)
569: {
570: PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;
572: PetscFunctionBeginUser;
573: PetscCall(PetscInitialize(&argc, &argv, 0, help));
574: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
575: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
576: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
577: nl = 1;
578: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
579: dim = 2;
580: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));
582: for (l = 0; l < nl; l++) {
583: if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
584: else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
585: else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
586: mx = mx * 2;
587: my = my * 2;
588: mz = mz * 2;
589: }
590: PetscCall(PetscFinalize());
591: return 0;
592: }
594: /*TEST
596: test:
597: suffix: 1d
598: args: -mx 10 -nl 6 -dim 1
600: test:
601: suffix: 2d
602: output_file: output/ex36_2d.out
603: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}
605: test:
606: suffix: 2dp1
607: nsize: 8
608: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
609: timeoutfactor: 2
611: test:
612: suffix: 2dp2
613: nsize: 8
614: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
615: timeoutfactor: 2
617: test:
618: suffix: 3d
619: args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
621: test:
622: suffix: 3dp1
623: nsize: 32
624: 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
626: TEST*/