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;
66: if (idx == 0) {
67: return 0;
68: } else if (idx == 1) { /* dam break */
69: DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
70: } else if (idx == 2) { /* stagnation in a corner */
71: DMDASetUniformCoordinates(da, -1.0, 1.0, 0.0, 1.0, -1.0, 1.0);
72: } else if (idx == 3) { /* nautilis */
73: DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
74: } else if (idx == 4) DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
76: DMGetCoordinateDM(da, &cda);
77: DMGetCoordinates(da, &Gcoords);
79: VecGetArray(Gcoords, &XX);
80: DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz);
81: DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
82: 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: VecRestoreArray(Gcoords, &XX);
193: return 0;
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;
206: DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
207: DMGetCoordinateDM(da, &cda);
208: DMGetCoordinates(da, &Gcoords);
210: DMDAVecGetArrayRead(cda, Gcoords, &XX);
211: 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: DMDAVecRestoreArrayRead(cda, Gcoords, &XX);
248: return 0;
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;
261: DMGetCoordinateDM(da, &cda);
262: DMGetCoordinates(da, &Gcoords);
264: DMDAVecGetArrayRead(cda, Gcoords, &XX);
265: DMDAVecGetArray(da, field, &FF);
267: 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: DMDAVecRestoreArray(da, field, &FF);
274: DMDAVecRestoreArrayRead(cda, Gcoords, &XX);
275: return 0;
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;
288: DMGetCoordinateDM(da, &cda);
289: DMGetCoordinates(da, &Gcoords);
291: DMDAVecGetArrayRead(cda, Gcoords, &XX);
292: DMDAVecGetArray(da, field, &FF);
294: 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: DMDAVecRestoreArray(da, field, &FF);
306: DMDAVecRestoreArrayRead(cda, Gcoords, &XX);
307: return 0;
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;
321: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac);
322: DMSetFromOptions(dac);
323: DMSetUp(dac);
325: DMRefine(dac, MPI_COMM_NULL, &daf);
326: DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
327: Mx--;
329: DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
330: DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
332: {
333: DM cdaf, cdac;
334: Vec coordsc, coordsf;
336: DMGetCoordinateDM(dac, &cdac);
337: DMGetCoordinateDM(daf, &cdaf);
339: DMGetCoordinates(dac, &coordsc);
340: DMGetCoordinates(daf, &coordsf);
342: DMCreateInterpolation(cdac, cdaf, &II, &scale);
343: MatInterpolate(II, coordsc, coordsf);
344: MatDestroy(&II);
345: VecDestroy(&scale);
346: }
348: DMCreateInterpolation(dac, daf, &INTERP, NULL);
350: DMCreateGlobalVector(dac, &ac);
351: VecSet(ac, 66.99);
353: DMCreateGlobalVector(daf, &af);
355: MatMult(INTERP, ac, af);
357: {
358: Vec afexact;
359: PetscReal nrm;
360: PetscInt N;
362: DMCreateGlobalVector(daf, &afexact);
363: VecSet(afexact, 66.99);
364: VecAXPY(afexact, -1.0, af); /* af <= af - afinterp */
365: VecNorm(afexact, NORM_2, &nrm);
366: VecGetSize(afexact, &N);
367: PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N)));
368: VecDestroy(&afexact);
369: }
371: PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL);
372: if (output) {
373: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv);
374: VecView(ac, vv);
375: PetscViewerDestroy(&vv);
377: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv);
378: VecView(af, vv);
379: PetscViewerDestroy(&vv);
380: }
382: MatDestroy(&INTERP);
383: DMDestroy(&dac);
384: DMDestroy(&daf);
385: VecDestroy(&ac);
386: VecDestroy(&af);
387: return 0;
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;
401: 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: DMSetFromOptions(dac);
403: DMSetUp(dac);
405: DMRefine(dac, MPI_COMM_NULL, &daf);
406: DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
407: Mx--;
408: My--;
410: DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE);
411: DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE);
413: /* apply conformal mappings */
414: map_id = 0;
415: PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL);
416: if (map_id >= 1) DAApplyConformalMapping(dac, map_id);
418: {
419: DM cdaf, cdac;
420: Vec coordsc, coordsf;
422: DMGetCoordinateDM(dac, &cdac);
423: DMGetCoordinateDM(daf, &cdaf);
425: DMGetCoordinates(dac, &coordsc);
426: DMGetCoordinates(daf, &coordsf);
428: DMCreateInterpolation(cdac, cdaf, &II, &scale);
429: MatInterpolate(II, coordsc, coordsf);
430: MatDestroy(&II);
431: VecDestroy(&scale);
432: }
434: DMCreateInterpolation(dac, daf, &INTERP, NULL);
436: DMCreateGlobalVector(dac, &ac);
437: DADefineXLinearField2D(dac, ac);
439: DMCreateGlobalVector(daf, &af);
440: MatMult(INTERP, ac, af);
442: {
443: Vec afexact;
444: PetscReal nrm;
445: PetscInt N;
447: DMCreateGlobalVector(daf, &afexact);
448: VecZeroEntries(afexact);
449: DADefineXLinearField2D(daf, afexact);
450: VecAXPY(afexact, -1.0, af); /* af <= af - afinterp */
451: VecNorm(afexact, NORM_2, &nrm);
452: VecGetSize(afexact, &N);
453: 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: VecDestroy(&afexact);
455: }
457: PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL);
458: if (output) {
459: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv);
460: VecView(ac, vv);
461: PetscViewerDestroy(&vv);
463: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv);
464: VecView(af, vv);
465: PetscViewerDestroy(&vv);
466: }
468: MatDestroy(&INTERP);
469: DMDestroy(&dac);
470: DMDestroy(&daf);
471: VecDestroy(&ac);
472: VecDestroy(&af);
473: return 0;
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;
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: DMSetFromOptions(dac);
490: DMSetUp(dac);
492: DMRefine(dac, MPI_COMM_NULL, &daf);
493: DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
494: Mx--;
495: My--;
496: Mz--;
498: DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
499: DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0);
501: /* apply trilinear mappings */
502: /*DAApplyTrilinearMapping(dac);*/
503: /* apply conformal mappings */
504: map_id = 0;
505: PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL);
506: if (map_id >= 1) DAApplyConformalMapping(dac, map_id);
508: {
509: DM cdaf, cdac;
510: Vec coordsc, coordsf;
512: DMGetCoordinateDM(dac, &cdac);
513: DMGetCoordinateDM(daf, &cdaf);
515: DMGetCoordinates(dac, &coordsc);
516: DMGetCoordinates(daf, &coordsf);
518: DMCreateInterpolation(cdac, cdaf, &II, &scale);
519: MatInterpolate(II, coordsc, coordsf);
520: MatDestroy(&II);
521: VecDestroy(&scale);
522: }
524: DMCreateInterpolation(dac, daf, &INTERP, NULL);
526: DMCreateGlobalVector(dac, &ac);
527: VecZeroEntries(ac);
528: DADefineXLinearField3D(dac, ac);
530: DMCreateGlobalVector(daf, &af);
531: VecZeroEntries(af);
533: MatMult(INTERP, ac, af);
535: {
536: Vec afexact;
537: PetscReal nrm;
538: PetscInt N;
540: DMCreateGlobalVector(daf, &afexact);
541: VecZeroEntries(afexact);
542: DADefineXLinearField3D(daf, afexact);
543: VecAXPY(afexact, -1.0, af); /* af <= af - afinterp */
544: VecNorm(afexact, NORM_2, &nrm);
545: VecGetSize(afexact, &N);
546: 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: VecDestroy(&afexact);
548: }
550: PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL);
551: if (output) {
552: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv);
553: VecView(ac, vv);
554: PetscViewerDestroy(&vv);
556: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv);
557: VecView(af, vv);
558: PetscViewerDestroy(&vv);
559: }
561: MatDestroy(&INTERP);
562: DMDestroy(&dac);
563: DMDestroy(&daf);
564: VecDestroy(&ac);
565: VecDestroy(&af);
566: return 0;
567: }
569: int main(int argc, char **argv)
570: {
571: PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;
574: PetscInitialize(&argc, &argv, 0, help);
575: PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0);
576: PetscOptionsGetInt(NULL, NULL, "-my", &my, 0);
577: PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0);
578: nl = 1;
579: PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0);
580: dim = 2;
581: PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0);
583: for (l = 0; l < nl; l++) {
584: if (dim == 1) da_test_RefineCoords1D(mx);
585: else if (dim == 2) da_test_RefineCoords2D(mx, my);
586: else if (dim == 3) da_test_RefineCoords3D(mx, my, mz);
587: mx = mx * 2;
588: my = my * 2;
589: mz = mz * 2;
590: }
591: 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*/