Actual source code: ex36.c
petsc-3.13.6 2020-09-29
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: {
59: PetscInt i,n;
60: PetscInt sx,nx,sy,ny,sz,nz,dim;
61: Vec Gcoords;
62: PetscScalar *XX;
63: PetscScalar xx,yy,zz;
64: DM cda;
67: if (idx==0) {
68: return(0);
69: } else if (idx==1) { /* dam break */
70: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
71: } else if (idx==2) { /* stagnation in a corner */
72: DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0);
73: } else if (idx==3) { /* nautilis */
74: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
75: } else if (idx==4) {
76: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
77: }
79: DMGetCoordinateDM(da,&cda);
80: DMGetCoordinates(da,&Gcoords);
82: VecGetArray(Gcoords,&XX);
83: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
84: DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
85: VecGetLocalSize(Gcoords,&n);
86: n = n / dim;
88: for (i=0; i<n; i++) {
89: if ((dim==3) && (idx!=2)) {
90: PetscScalar Ni[8];
91: PetscScalar xi = XX[dim*i];
92: PetscScalar eta = XX[dim*i+1];
93: PetscScalar zeta = XX[dim*i+2];
94: PetscScalar xn[] = {-1.0,1.0,-1.0,1.0, -1.0,1.0,-1.0,1.0 };
95: PetscScalar yn[] = {-1.0,-1.0,1.0,1.0, -1.0,-1.0,1.0,1.0 };
96: PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0, 0.1,4.0,0.2,1.0 };
97: PetscInt p;
99: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
100: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
101: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
102: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
104: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
105: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
106: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
107: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
109: xx = yy = zz = 0.0;
110: for (p=0; p<8; p++) {
111: xx += Ni[p]*xn[p];
112: yy += Ni[p]*yn[p];
113: zz += Ni[p]*zn[p];
114: }
115: XX[dim*i] = xx;
116: XX[dim*i+1] = yy;
117: XX[dim*i+2] = zz;
118: }
120: if (idx==1) {
121: CCmplx zeta,t1,t2;
123: xx = XX[dim*i] - 0.8;
124: yy = XX[dim*i+1] + 1.5;
126: zeta.real = PetscRealPart(xx);
127: zeta.imag = PetscRealPart(yy);
129: t1 = CCmplxPow(zeta,-1.0);
130: t2 = CCmplxAdd(zeta,t1);
132: XX[dim*i] = CCmplxRe(t2);
133: XX[dim*i+1] = CCmplxIm(t2);
134: } else if (idx==2) {
135: CCmplx zeta,t1;
137: xx = XX[dim*i];
138: yy = XX[dim*i+1];
139: zeta.real = PetscRealPart(xx);
140: zeta.imag = PetscRealPart(yy);
142: t1 = CCmplxSqrt(zeta);
143: XX[dim*i] = CCmplxRe(t1);
144: XX[dim*i+1] = CCmplxIm(t1);
145: } else if (idx==3) {
146: CCmplx zeta,t1,t2;
148: xx = XX[dim*i] - 0.8;
149: yy = XX[dim*i+1] + 1.5;
151: zeta.real = PetscRealPart(xx);
152: zeta.imag = PetscRealPart(yy);
153: t1 = CCmplxPow(zeta,-1.0);
154: t2 = CCmplxAdd(zeta,t1);
155: XX[dim*i] = CCmplxRe(t2);
156: XX[dim*i+1] = CCmplxIm(t2);
158: xx = XX[dim*i];
159: yy = XX[dim*i+1];
160: zeta.real = PetscRealPart(xx);
161: zeta.imag = PetscRealPart(yy);
162: t1 = CCmplxExp(zeta);
163: XX[dim*i] = CCmplxRe(t1);
164: XX[dim*i+1] = CCmplxIm(t1);
166: xx = XX[dim*i] + 0.4;
167: yy = XX[dim*i+1];
168: zeta.real = PetscRealPart(xx);
169: zeta.imag = PetscRealPart(yy);
170: t1 = CCmplxPow(zeta,2.0);
171: XX[dim*i] = CCmplxRe(t1);
172: XX[dim*i+1] = CCmplxIm(t1);
173: } else if (idx==4) {
174: PetscScalar Ni[4];
175: PetscScalar xi = XX[dim*i];
176: PetscScalar eta = XX[dim*i+1];
177: PetscScalar xn[] = {0.0,2.0,0.2,3.5};
178: PetscScalar yn[] = {-1.3,0.0,2.0,4.0};
179: PetscInt p;
181: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
182: Ni[1] = 0.25*(1.0+xi)*(1.0-eta);
183: Ni[2] = 0.25*(1.0-xi)*(1.0+eta);
184: Ni[3] = 0.25*(1.0+xi)*(1.0+eta);
186: xx = yy = 0.0;
187: for (p=0; p<4; p++) {
188: xx += Ni[p]*xn[p];
189: yy += Ni[p]*yn[p];
190: }
191: XX[dim*i] = xx;
192: XX[dim*i+1] = yy;
193: }
194: }
195: VecRestoreArray(Gcoords,&XX);
196: return(0);
197: }
199: PetscErrorCode DAApplyTrilinearMapping(DM da)
200: {
202: PetscInt i,j,k;
203: PetscInt sx,nx,sy,ny,sz,nz;
204: Vec Gcoords;
205: DMDACoor3d ***XX;
206: PetscScalar xx,yy,zz;
207: DM cda;
210: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
211: DMGetCoordinateDM(da,&cda);
212: DMGetCoordinates(da,&Gcoords);
214: DMDAVecGetArrayRead(cda,Gcoords,&XX);
215: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
217: for (i=sx; i<sx+nx; i++) {
218: for (j=sy; j<sy+ny; j++) {
219: for (k=sz; k<sz+nz; k++) {
220: PetscScalar Ni[8];
221: PetscScalar xi = XX[k][j][i].x;
222: PetscScalar eta = XX[k][j][i].y;
223: PetscScalar zeta = XX[k][j][i].z;
224: PetscScalar xn[] = {0.0,2.0,0.2,3.5, 0.0,2.1,0.23,3.125 };
225: PetscScalar yn[] = {-1.3,0.0,2.0,4.0, -1.45,-0.1,2.24,3.79 };
226: PetscScalar zn[] = {0.0,0.3,-0.1,0.123, 0.956,1.32,1.12,0.798 };
227: PetscInt p;
229: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
230: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
231: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
232: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
234: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
235: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
236: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
237: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
239: xx = yy = zz = 0.0;
240: for (p=0; p<8; p++) {
241: xx += Ni[p]*xn[p];
242: yy += Ni[p]*yn[p];
243: zz += Ni[p]*zn[p];
244: }
245: XX[k][j][i].x = xx;
246: XX[k][j][i].y = yy;
247: XX[k][j][i].z = zz;
248: }
249: }
250: }
251: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
252: return(0);
253: }
255: PetscErrorCode DADefineXLinearField2D(DM da,Vec field)
256: {
258: PetscInt i,j;
259: PetscInt sx,nx,sy,ny;
260: Vec Gcoords;
261: DMDACoor2d **XX;
262: PetscScalar **FF;
263: DM cda;
266: DMGetCoordinateDM(da,&cda);
267: DMGetCoordinates(da,&Gcoords);
269: DMDAVecGetArrayRead(cda,Gcoords,&XX);
270: DMDAVecGetArray(da,field,&FF);
272: DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);
274: for (i=sx; i<sx+nx; i++) {
275: for (j=sy; j<sy+ny; j++) {
276: 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;
277: }
278: }
280: DMDAVecRestoreArray(da,field,&FF);
281: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
282: return(0);
283: }
285: PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
286: {
288: PetscInt i,j,k;
289: PetscInt sx,nx,sy,ny,sz,nz;
290: Vec Gcoords;
291: DMDACoor3d ***XX;
292: PetscScalar ***FF;
293: DM cda;
296: DMGetCoordinateDM(da,&cda);
297: DMGetCoordinates(da,&Gcoords);
299: DMDAVecGetArrayRead(cda,Gcoords,&XX);
300: DMDAVecGetArray(da,field,&FF);
302: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
304: for (k=sz; k<sz+nz; k++) {
305: for (j=sy; j<sy+ny; j++) {
306: for (i=sx; i<sx+nx; i++) {
307: FF[k][j][i] = 10.0
308: + 4.05 * XX[k][j][i].x
309: + 5.50 * XX[k][j][i].y
310: + 1.33 * XX[k][j][i].z
311: + 2.03 * XX[k][j][i].x * XX[k][j][i].y
312: + 0.03 * XX[k][j][i].x * XX[k][j][i].z
313: + 0.83 * XX[k][j][i].y * XX[k][j][i].z
314: + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
315: }
316: }
317: }
319: DMDAVecRestoreArray(da,field,&FF);
320: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
321: return(0);
322: }
324: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
325: {
327: DM dac,daf;
328: PetscViewer vv;
329: Vec ac,af;
330: PetscInt Mx;
331: Mat II,INTERP;
332: Vec scale;
333: PetscBool output = PETSC_FALSE;
336: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,mx+1,1, /* 1 dof */1, /* stencil = 1 */NULL,&dac);
337: DMSetFromOptions(dac);
338: DMSetUp(dac);
340: DMRefine(dac,MPI_COMM_NULL,&daf);
341: DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
342: Mx--;
344: DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
345: DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
347: {
348: DM cdaf,cdac;
349: Vec coordsc,coordsf;
351: DMGetCoordinateDM(dac,&cdac);
352: DMGetCoordinateDM(daf,&cdaf);
354: DMGetCoordinates(dac,&coordsc);
355: DMGetCoordinates(daf,&coordsf);
357: DMCreateInterpolation(cdac,cdaf,&II,&scale);
358: MatInterpolate(II,coordsc,coordsf);
359: MatDestroy(&II);
360: VecDestroy(&scale);
361: }
363: DMCreateInterpolation(dac,daf,&INTERP,NULL);
365: DMCreateGlobalVector(dac,&ac);
366: VecSet(ac,66.99);
368: DMCreateGlobalVector(daf,&af);
370: MatMult(INTERP,ac, af);
372: {
373: Vec afexact;
374: PetscReal nrm;
375: PetscInt N;
377: DMCreateGlobalVector(daf,&afexact);
378: VecSet(afexact,66.99);
379: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
380: VecNorm(afexact,NORM_2,&nrm);
381: VecGetSize(afexact,&N);
382: PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)));
383: VecDestroy(&afexact);
384: }
386: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
387: if (output) {
388: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);
389: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
390: DMView(dac, vv);
391: VecView(ac, vv);
392: PetscViewerPopFormat(vv);
393: PetscViewerDestroy(&vv);
395: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);
396: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
397: DMView(daf, vv);
398: VecView(af, vv);
399: PetscViewerPopFormat(vv);
400: PetscViewerDestroy(&vv);
401: }
403: MatDestroy(&INTERP);
404: DMDestroy(&dac);
405: DMDestroy(&daf);
406: VecDestroy(&ac);
407: VecDestroy(&af);
408: return(0);
409: }
411: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
412: {
414: DM dac,daf;
415: PetscViewer vv;
416: Vec ac,af;
417: PetscInt map_id,Mx,My;
418: Mat II,INTERP;
419: Vec scale;
420: PetscBool output = PETSC_FALSE;
423: 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);
424: DMSetFromOptions(dac);
425: DMSetUp(dac);
427: DMRefine(dac,MPI_COMM_NULL,&daf);
428: DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
429: Mx--; My--;
431: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
432: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
434: /* apply conformal mappings */
435: map_id = 0;
436: PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
437: if (map_id >= 1) {
438: DAApplyConformalMapping(dac,map_id);
439: }
441: {
442: DM cdaf,cdac;
443: Vec coordsc,coordsf;
445: DMGetCoordinateDM(dac,&cdac);
446: DMGetCoordinateDM(daf,&cdaf);
448: DMGetCoordinates(dac,&coordsc);
449: DMGetCoordinates(daf,&coordsf);
451: DMCreateInterpolation(cdac,cdaf,&II,&scale);
452: MatInterpolate(II,coordsc,coordsf);
453: MatDestroy(&II);
454: VecDestroy(&scale);
455: }
458: DMCreateInterpolation(dac,daf,&INTERP,NULL);
460: DMCreateGlobalVector(dac,&ac);
461: DADefineXLinearField2D(dac,ac);
463: DMCreateGlobalVector(daf,&af);
464: MatMult(INTERP,ac, af);
466: {
467: Vec afexact;
468: PetscReal nrm;
469: PetscInt N;
471: DMCreateGlobalVector(daf,&afexact);
472: VecZeroEntries(afexact);
473: DADefineXLinearField2D(daf,afexact);
474: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
475: VecNorm(afexact,NORM_2,&nrm);
476: VecGetSize(afexact,&N);
477: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)));
478: VecDestroy(&afexact);
479: }
481: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
482: if (output) {
483: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);
484: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
485: DMView(dac, vv);
486: VecView(ac, vv);
487: PetscViewerPopFormat(vv);
488: PetscViewerDestroy(&vv);
490: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);
491: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
492: DMView(daf, vv);
493: VecView(af, vv);
494: PetscViewerPopFormat(vv);
495: PetscViewerDestroy(&vv);
496: }
498: MatDestroy(&INTERP);
499: DMDestroy(&dac);
500: DMDestroy(&daf);
501: VecDestroy(&ac);
502: VecDestroy(&af);
503: return(0);
504: }
506: PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
507: {
509: DM dac,daf;
510: PetscViewer vv;
511: Vec ac,af;
512: PetscInt map_id,Mx,My,Mz;
513: Mat II,INTERP;
514: Vec scale;
515: PetscBool output = PETSC_FALSE;
518: 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 */
519: 1, /* stencil = 1 */NULL,NULL,NULL,&dac);
520: DMSetFromOptions(dac);
521: DMSetUp(dac);
523: DMRefine(dac,MPI_COMM_NULL,&daf);
524: DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);
525: Mx--; My--; Mz--;
527: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
528: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);
530: /* apply trilinear mappings */
531: /*DAApplyTrilinearMapping(dac);*/
532: /* apply conformal mappings */
533: map_id = 0;
534: PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
535: if (map_id >= 1) {
536: DAApplyConformalMapping(dac,map_id);
537: }
539: {
540: DM cdaf,cdac;
541: Vec coordsc,coordsf;
543: DMGetCoordinateDM(dac,&cdac);
544: DMGetCoordinateDM(daf,&cdaf);
546: DMGetCoordinates(dac,&coordsc);
547: DMGetCoordinates(daf,&coordsf);
549: DMCreateInterpolation(cdac,cdaf,&II,&scale);
550: MatInterpolate(II,coordsc,coordsf);
551: MatDestroy(&II);
552: VecDestroy(&scale);
553: }
555: DMCreateInterpolation(dac,daf,&INTERP,NULL);
557: DMCreateGlobalVector(dac,&ac);
558: VecZeroEntries(ac);
559: DADefineXLinearField3D(dac,ac);
561: DMCreateGlobalVector(daf,&af);
562: VecZeroEntries(af);
564: MatMult(INTERP,ac, af);
566: {
567: Vec afexact;
568: PetscReal nrm;
569: PetscInt N;
571: DMCreateGlobalVector(daf,&afexact);
572: VecZeroEntries(afexact);
573: DADefineXLinearField3D(daf,afexact);
574: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
575: VecNorm(afexact,NORM_2,&nrm);
576: VecGetSize(afexact,&N);
577: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D x %D]=>[%D x %D x %D], interp err = %1.4e\n",mx,my,mz,Mx,My,Mz,(double)(nrm/PetscSqrtReal((PetscReal)N)));
578: VecDestroy(&afexact);
579: }
581: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
582: if (output) {
583: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);
584: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
585: DMView(dac, vv);
586: VecView(ac, vv);
587: PetscViewerPopFormat(vv);
588: PetscViewerDestroy(&vv);
590: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);
591: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
592: DMView(daf, vv);
593: VecView(af, vv);
594: PetscViewerPopFormat(vv);
595: PetscViewerDestroy(&vv);
596: }
598: MatDestroy(&INTERP);
599: DMDestroy(&dac);
600: DMDestroy(&daf);
601: VecDestroy(&ac);
602: VecDestroy(&af);
603: return(0);
604: }
606: int main(int argc,char **argv)
607: {
609: PetscInt mx = 2,my = 2,mz = 2,l,nl,dim;
611: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
612: PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);
613: PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);
614: PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);
615: nl = 1;
616: PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0);
617: dim = 2;
618: PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0);
620: for (l=0; l<nl; l++) {
621: if (dim==1) da_test_RefineCoords1D(mx);
622: else if (dim==2) da_test_RefineCoords2D(mx,my);
623: else if (dim==3) da_test_RefineCoords3D(mx,my,mz);
624: mx = mx * 2;
625: my = my * 2;
626: mz = mz * 2;
627: }
628: PetscFinalize();
629: return ierr;
630: }
633: /*TEST
635: test:
636: suffix: 1d
637: args: -mx 10 -nl 6 -dim 1
639: test:
640: suffix: 2d
642: test:
643: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 0
645: test:
646: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 1
648: test:
649: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 2
651: test:
652: args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 3
654: test:
655: suffix: 2dp1
656: nsize: 8
657: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
658: timeoutfactor: 2
660: test:
661: suffix: 2dp2
662: nsize: 8
663: args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
664: timeoutfactor: 2
666: test:
667: suffix: 3d
668: args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
670: test:
671: suffix: 3dp1
672: nsize: 32
673: 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
675: TEST*/