Actual source code: ex36.c
petsc-3.7.7 2017-09-25
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 = atan2(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 = atan2(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: }
58: PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx)
59: {
61: PetscInt i,n;
62: PetscInt sx,nx,sy,ny,sz,nz,dim;
63: Vec Gcoords;
64: PetscScalar *XX;
65: PetscScalar xx,yy,zz;
66: DM cda;
69: if (idx==0) {
70: return(0);
71: } else if (idx==1) { /* dam break */
72: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
73: } else if (idx==2) { /* stagnation in a corner */
74: DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0);
75: } else if (idx==3) { /* nautilis */
76: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
77: } else if (idx==4) {
78: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
79: }
81: DMGetCoordinateDM(da,&cda);
82: DMGetCoordinates(da,&Gcoords);
84: VecGetArray(Gcoords,&XX);
85: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
86: DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
87: VecGetLocalSize(Gcoords,&n);
88: n = n / dim;
90: for (i=0; i<n; i++) {
91: if ((dim==3) && (idx!=2)) {
92: PetscScalar Ni[8];
93: PetscScalar xi = XX[dim*i];
94: PetscScalar eta = XX[dim*i+1];
95: PetscScalar zeta = XX[dim*i+2];
96: PetscScalar xn[] = {-1.0,1.0,-1.0,1.0, -1.0,1.0,-1.0,1.0 };
97: PetscScalar yn[] = {-1.0,-1.0,1.0,1.0, -1.0,-1.0,1.0,1.0 };
98: PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0, 0.1,4.0,0.2,1.0 };
99: PetscInt p;
101: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
102: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
103: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
104: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
106: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
107: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
108: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
109: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
111: xx = yy = zz = 0.0;
112: for (p=0; p<8; p++) {
113: xx += Ni[p]*xn[p];
114: yy += Ni[p]*yn[p];
115: zz += Ni[p]*zn[p];
116: }
117: XX[dim*i] = xx;
118: XX[dim*i+1] = yy;
119: XX[dim*i+2] = zz;
120: }
122: if (idx==1) {
123: CCmplx zeta,t1,t2;
125: xx = XX[dim*i] - 0.8;
126: yy = XX[dim*i+1] + 1.5;
128: zeta.real = PetscRealPart(xx);
129: zeta.imag = PetscRealPart(yy);
131: t1 = CCmplxPow(zeta,-1.0);
132: t2 = CCmplxAdd(zeta,t1);
134: XX[dim*i] = CCmplxRe(t2);
135: XX[dim*i+1] = CCmplxIm(t2);
136: } else if (idx==2) {
137: CCmplx zeta,t1;
139: xx = XX[dim*i];
140: yy = XX[dim*i+1];
141: zeta.real = PetscRealPart(xx);
142: zeta.imag = PetscRealPart(yy);
144: t1 = CCmplxSqrt(zeta);
145: XX[dim*i] = CCmplxRe(t1);
146: XX[dim*i+1] = CCmplxIm(t1);
147: } else if (idx==3) {
148: CCmplx zeta,t1,t2;
150: xx = XX[dim*i] - 0.8;
151: yy = XX[dim*i+1] + 1.5;
153: zeta.real = PetscRealPart(xx);
154: zeta.imag = PetscRealPart(yy);
155: t1 = CCmplxPow(zeta,-1.0);
156: t2 = CCmplxAdd(zeta,t1);
157: XX[dim*i] = CCmplxRe(t2);
158: XX[dim*i+1] = CCmplxIm(t2);
160: xx = XX[dim*i];
161: yy = XX[dim*i+1];
162: zeta.real = PetscRealPart(xx);
163: zeta.imag = PetscRealPart(yy);
164: t1 = CCmplxExp(zeta);
165: XX[dim*i] = CCmplxRe(t1);
166: XX[dim*i+1] = CCmplxIm(t1);
168: xx = XX[dim*i] + 0.4;
169: yy = XX[dim*i+1];
170: zeta.real = PetscRealPart(xx);
171: zeta.imag = PetscRealPart(yy);
172: t1 = CCmplxPow(zeta,2.0);
173: XX[dim*i] = CCmplxRe(t1);
174: XX[dim*i+1] = CCmplxIm(t1);
175: } else if (idx==4) {
176: PetscScalar Ni[4];
177: PetscScalar xi = XX[dim*i];
178: PetscScalar eta = XX[dim*i+1];
179: PetscScalar xn[] = {0.0,2.0,0.2,3.5};
180: PetscScalar yn[] = {-1.3,0.0,2.0,4.0};
181: PetscInt p;
183: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
184: Ni[1] = 0.25*(1.0+xi)*(1.0-eta);
185: Ni[2] = 0.25*(1.0-xi)*(1.0+eta);
186: Ni[3] = 0.25*(1.0+xi)*(1.0+eta);
188: xx = yy = 0.0;
189: for (p=0; p<4; p++) {
190: xx += Ni[p]*xn[p];
191: yy += Ni[p]*yn[p];
192: }
193: XX[dim*i] = xx;
194: XX[dim*i+1] = yy;
195: }
196: }
197: VecRestoreArray(Gcoords,&XX);
198: return(0);
199: }
203: PetscErrorCode DAApplyTrilinearMapping(DM da)
204: {
206: PetscInt i,j,k;
207: PetscInt sx,nx,sy,ny,sz,nz;
208: Vec Gcoords;
209: DMDACoor3d ***XX;
210: PetscScalar xx,yy,zz;
211: DM cda;
214: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
215: DMGetCoordinateDM(da,&cda);
216: DMGetCoordinates(da,&Gcoords);
218: DMDAVecGetArrayRead(cda,Gcoords,&XX);
219: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
221: for (i=sx; i<sx+nx; i++) {
222: for (j=sy; j<sy+ny; j++) {
223: for (k=sz; k<sz+nz; k++) {
224: PetscScalar Ni[8];
225: PetscScalar xi = XX[k][j][i].x;
226: PetscScalar eta = XX[k][j][i].y;
227: PetscScalar zeta = XX[k][j][i].z;
228: PetscScalar xn[] = {0.0,2.0,0.2,3.5, 0.0,2.1,0.23,3.125 };
229: PetscScalar yn[] = {-1.3,0.0,2.0,4.0, -1.45,-0.1,2.24,3.79 };
230: PetscScalar zn[] = {0.0,0.3,-0.1,0.123, 0.956,1.32,1.12,0.798 };
231: PetscInt p;
233: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
234: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
235: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
236: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
238: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
239: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
240: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
241: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
243: xx = yy = zz = 0.0;
244: for (p=0; p<8; p++) {
245: xx += Ni[p]*xn[p];
246: yy += Ni[p]*yn[p];
247: zz += Ni[p]*zn[p];
248: }
249: XX[k][j][i].x = xx;
250: XX[k][j][i].y = yy;
251: XX[k][j][i].z = zz;
252: }
253: }
254: }
255: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
256: return(0);
257: }
261: PetscErrorCode DADefineXLinearField2D(DM da,Vec field)
262: {
264: PetscInt i,j;
265: PetscInt sx,nx,sy,ny;
266: Vec Gcoords;
267: DMDACoor2d **XX;
268: PetscScalar **FF;
269: DM cda;
272: DMGetCoordinateDM(da,&cda);
273: DMGetCoordinates(da,&Gcoords);
275: DMDAVecGetArrayRead(cda,Gcoords,&XX);
276: DMDAVecGetArray(da,field,&FF);
278: DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);
280: for (i=sx; i<sx+nx; i++) {
281: for (j=sy; j<sy+ny; j++) {
282: 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;
283: }
284: }
286: DMDAVecRestoreArray(da,field,&FF);
287: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
288: return(0);
289: }
293: PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
294: {
296: PetscInt i,j,k;
297: PetscInt sx,nx,sy,ny,sz,nz;
298: Vec Gcoords;
299: DMDACoor3d ***XX;
300: PetscScalar ***FF;
301: DM cda;
304: DMGetCoordinateDM(da,&cda);
305: DMGetCoordinates(da,&Gcoords);
307: DMDAVecGetArrayRead(cda,Gcoords,&XX);
308: DMDAVecGetArray(da,field,&FF);
310: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
312: for (k=sz; k<sz+nz; k++) {
313: for (j=sy; j<sy+ny; j++) {
314: for (i=sx; i<sx+nx; i++) {
315: FF[k][j][i] = 10.0
316: + 4.05 * XX[k][j][i].x
317: + 5.50 * XX[k][j][i].y
318: + 1.33 * XX[k][j][i].z
319: + 2.03 * XX[k][j][i].x * XX[k][j][i].y
320: + 0.03 * XX[k][j][i].x * XX[k][j][i].z
321: + 0.83 * XX[k][j][i].y * XX[k][j][i].z
322: + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
323: }
324: }
325: }
327: DMDAVecRestoreArray(da,field,&FF);
328: DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
329: return(0);
330: }
334: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
335: {
337: DM dac,daf;
338: PetscViewer vv;
339: Vec ac,af;
340: PetscInt Mx;
341: Mat II,INTERP;
342: Vec scale;
343: PetscBool output = PETSC_FALSE;
346: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,
347: mx+1,
348: 1, /* 1 dof */
349: 1, /* stencil = 1 */
350: NULL,
351: &dac);
352: DMSetFromOptions(dac);
354: DMRefine(dac,MPI_COMM_NULL,&daf);
355: DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
356: Mx--;
358: DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
359: DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
361: {
362: DM cdaf,cdac;
363: Vec coordsc,coordsf;
365: DMGetCoordinateDM(dac,&cdac);
366: DMGetCoordinateDM(daf,&cdaf);
368: DMGetCoordinates(dac,&coordsc);
369: DMGetCoordinates(daf,&coordsf);
371: DMCreateInterpolation(cdac,cdaf,&II,&scale);
372: MatInterpolate(II,coordsc,coordsf);
373: MatDestroy(&II);
374: VecDestroy(&scale);
375: }
377: DMCreateInterpolation(dac,daf,&INTERP,NULL);
379: DMCreateGlobalVector(dac,&ac);
380: VecSet(ac,66.99);
382: DMCreateGlobalVector(daf,&af);
384: MatMult(INTERP,ac, af);
386: {
387: Vec afexact;
388: PetscReal nrm;
389: PetscInt N;
391: DMCreateGlobalVector(daf,&afexact);
392: VecSet(afexact,66.99);
393: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
394: VecNorm(afexact,NORM_2,&nrm);
395: VecGetSize(afexact,&N);
396: PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)));
397: VecDestroy(&afexact);
398: }
400: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
401: if (output) {
402: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);
403: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
404: DMView(dac, vv);
405: VecView(ac, vv);
406: PetscViewerPopFormat(vv);
407: PetscViewerDestroy(&vv);
409: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);
410: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
411: DMView(daf, vv);
412: VecView(af, vv);
413: PetscViewerPopFormat(vv);
414: PetscViewerDestroy(&vv);
415: }
417: MatDestroy(&INTERP);
418: DMDestroy(&dac);
419: DMDestroy(&daf);
420: VecDestroy(&ac);
421: VecDestroy(&af);
422: return(0);
423: }
427: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
428: {
430: DM dac,daf;
431: PetscViewer vv;
432: Vec ac,af;
433: PetscInt map_id,Mx,My;
434: Mat II,INTERP;
435: Vec scale;
436: PetscBool output = PETSC_FALSE;
439: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
440: mx+1, my+1,
441: PETSC_DECIDE, PETSC_DECIDE,
442: 1, /* 1 dof */
443: 1, /* stencil = 1 */
444: NULL, NULL,
445: &dac);
446: DMSetFromOptions(dac);
448: DMRefine(dac,MPI_COMM_NULL,&daf);
449: DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
450: Mx--; My--;
452: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
453: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
455: /* apply conformal mappings */
456: map_id = 0;
457: PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
458: if (map_id >= 1) {
459: DAApplyConformalMapping(dac,map_id);
460: }
462: {
463: DM cdaf,cdac;
464: Vec coordsc,coordsf;
466: DMGetCoordinateDM(dac,&cdac);
467: DMGetCoordinateDM(daf,&cdaf);
469: DMGetCoordinates(dac,&coordsc);
470: DMGetCoordinates(daf,&coordsf);
472: DMCreateInterpolation(cdac,cdaf,&II,&scale);
473: MatInterpolate(II,coordsc,coordsf);
474: MatDestroy(&II);
475: VecDestroy(&scale);
476: }
479: DMCreateInterpolation(dac,daf,&INTERP,NULL);
481: DMCreateGlobalVector(dac,&ac);
482: DADefineXLinearField2D(dac,ac);
484: DMCreateGlobalVector(daf,&af);
485: MatMult(INTERP,ac, af);
487: {
488: Vec afexact;
489: PetscReal nrm;
490: PetscInt N;
492: DMCreateGlobalVector(daf,&afexact);
493: VecZeroEntries(afexact);
494: DADefineXLinearField2D(daf,afexact);
495: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
496: VecNorm(afexact,NORM_2,&nrm);
497: VecGetSize(afexact,&N);
498: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)));
499: VecDestroy(&afexact);
500: }
502: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
503: if (output) {
504: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);
505: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
506: DMView(dac, vv);
507: VecView(ac, vv);
508: PetscViewerPopFormat(vv);
509: PetscViewerDestroy(&vv);
511: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);
512: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
513: DMView(daf, vv);
514: VecView(af, vv);
515: PetscViewerPopFormat(vv);
516: PetscViewerDestroy(&vv);
517: }
519: MatDestroy(&INTERP);
520: DMDestroy(&dac);
521: DMDestroy(&daf);
522: VecDestroy(&ac);
523: VecDestroy(&af);
524: return(0);
525: }
529: PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
530: {
532: DM dac,daf;
533: PetscViewer vv;
534: Vec ac,af;
535: PetscInt map_id,Mx,My,Mz;
536: Mat II,INTERP;
537: Vec scale;
538: PetscBool output = PETSC_FALSE;
541: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
542: mx+1, my+1,mz+1,
543: PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
544: 1, /* 1 dof */
545: 1, /* stencil = 1 */
546: NULL,NULL,NULL,
547: &dac);
548: DMSetFromOptions(dac);
550: DMRefine(dac,MPI_COMM_NULL,&daf);
551: DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);
552: Mx--; My--; Mz--;
554: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
555: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);
557: /* apply trilinear mappings */
558: /*DAApplyTrilinearMapping(dac);*/
559: /* apply conformal mappings */
560: map_id = 0;
561: PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
562: if (map_id >= 1) {
563: DAApplyConformalMapping(dac,map_id);
564: }
566: {
567: DM cdaf,cdac;
568: Vec coordsc,coordsf;
570: DMGetCoordinateDM(dac,&cdac);
571: DMGetCoordinateDM(daf,&cdaf);
573: DMGetCoordinates(dac,&coordsc);
574: DMGetCoordinates(daf,&coordsf);
576: DMCreateInterpolation(cdac,cdaf,&II,&scale);
577: MatInterpolate(II,coordsc,coordsf);
578: MatDestroy(&II);
579: VecDestroy(&scale);
580: }
582: DMCreateInterpolation(dac,daf,&INTERP,NULL);
584: DMCreateGlobalVector(dac,&ac);
585: VecZeroEntries(ac);
586: DADefineXLinearField3D(dac,ac);
588: DMCreateGlobalVector(daf,&af);
589: VecZeroEntries(af);
591: MatMult(INTERP,ac, af);
593: {
594: Vec afexact;
595: PetscReal nrm;
596: PetscInt N;
598: DMCreateGlobalVector(daf,&afexact);
599: VecZeroEntries(afexact);
600: DADefineXLinearField3D(daf,afexact);
601: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
602: VecNorm(afexact,NORM_2,&nrm);
603: VecGetSize(afexact,&N);
604: 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)));
605: VecDestroy(&afexact);
606: }
608: PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
609: if (output) {
610: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);
611: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
612: DMView(dac, vv);
613: VecView(ac, vv);
614: PetscViewerPopFormat(vv);
615: PetscViewerDestroy(&vv);
617: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);
618: PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);
619: DMView(daf, vv);
620: VecView(af, vv);
621: PetscViewerPopFormat(vv);
622: PetscViewerDestroy(&vv);
623: }
625: MatDestroy(&INTERP);
626: DMDestroy(&dac);
627: DMDestroy(&daf);
628: VecDestroy(&ac);
629: VecDestroy(&af);
630: return(0);
631: }
635: int main(int argc,char **argv)
636: {
638: PetscInt mx,my,mz,l,nl,dim;
640: PetscInitialize(&argc,&argv,0,help);
642: mx = my = mz = 2;
643: PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);
644: PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);
645: PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);
646: nl = 1;
647: PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0);
648: dim = 2;
649: PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0);
651: for (l=0; l<nl; l++) {
652: if (dim==1) da_test_RefineCoords1D(mx);
653: else if (dim==2) da_test_RefineCoords2D(mx,my);
654: else if (dim==3) da_test_RefineCoords3D(mx,my,mz);
655: mx = mx * 2;
656: my = my * 2;
657: mz = mz * 2;
658: }
660: PetscFinalize();
661: return 0;
662: }