Actual source code: ex36.c
petsc-3.4.2 2013-07-02
2: static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";
4: #include <petscdmda.h>
6: typedef struct _n_CCmplx CCmplx;
7: struct _n_CCmplx {
8: double real;
9: double imag;
10: };
12: CCmplx CCmplxPow(CCmplx a,double n)
13: {
14: CCmplx b;
15: double r,theta;
16: r = sqrt(a.real*a.real + a.imag*a.imag);
17: theta = atan2(a.imag,a.real);
18: b.real = pow(r,n) * cos(n*theta);
19: b.imag = pow(r,n) * sin(n*theta);
20: return b;
21: }
22: CCmplx CCmplxExp(CCmplx a)
23: {
24: CCmplx b;
25: b.real = exp(a.real) * cos(a.imag);
26: b.imag = exp(a.real) * sin(a.imag);
27: return b;
28: }
29: CCmplx CCmplxSqrt(CCmplx a)
30: {
31: CCmplx b;
32: double r,theta;
33: r = sqrt(a.real*a.real + a.imag*a.imag);
34: theta = atan2(a.imag,a.real);
35: b.real = sqrt(r) * cos(0.5*theta);
36: b.imag = sqrt(r) * sin(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: }
57: PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx)
58: {
60: PetscInt i,n;
61: PetscInt sx,nx,sy,ny,sz,nz,dim;
62: Vec Gcoords;
63: PetscScalar *XX;
64: PetscScalar xx,yy,zz;
65: DM cda;
68: if (idx==0) {
69: return(0);
70: } else if (idx==1) { /* dam break */
71: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
72: } else if (idx==2) { /* stagnation in a corner */
73: DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0);
74: } else if (idx==3) { /* nautilis */
75: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
76: } else if (idx==4) {
77: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
78: }
80: DMGetCoordinateDM(da,&cda);
81: DMGetCoordinates(da,&Gcoords);
83: VecGetArray(Gcoords,&XX);
84: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
85: DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
86: VecGetLocalSize(Gcoords,&n);
87: n = n / dim;
89: for (i=0; i<n; i++) {
90: if ((dim==3) && (idx!=2)) {
91: PetscScalar Ni[8];
92: PetscScalar xi = XX[dim*i];
93: PetscScalar eta = XX[dim*i+1];
94: PetscScalar zeta = XX[dim*i+2];
95: PetscScalar xn[] = {-1.0,1.0,-1.0,1.0, -1.0,1.0,-1.0,1.0 };
96: PetscScalar yn[] = {-1.0,-1.0,1.0,1.0, -1.0,-1.0,1.0,1.0 };
97: PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0, 0.1,4.0,0.2,1.0 };
98: PetscInt p;
100: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
101: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
102: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
103: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
105: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
106: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
107: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
108: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
110: xx = yy = zz = 0.0;
111: for (p=0; p<8; p++) {
112: xx += Ni[p]*xn[p];
113: yy += Ni[p]*yn[p];
114: zz += Ni[p]*zn[p];
115: }
116: XX[dim*i] = xx;
117: XX[dim*i+1] = yy;
118: XX[dim*i+2] = zz;
119: }
121: if (idx==1) {
122: CCmplx zeta,t1,t2;
124: xx = XX[dim*i] - 0.8;
125: yy = XX[dim*i+1] + 1.5;
127: zeta.real = PetscRealPart(xx);
128: zeta.imag = PetscRealPart(yy);
130: t1 = CCmplxPow(zeta,-1.0);
131: t2 = CCmplxAdd(zeta,t1);
133: XX[dim*i] = CCmplxRe(t2);
134: XX[dim*i+1] = CCmplxIm(t2);
135: } else if (idx==2) {
136: CCmplx zeta,t1;
138: xx = XX[dim*i];
139: yy = XX[dim*i+1];
140: zeta.real = PetscRealPart(xx);
141: zeta.imag = PetscRealPart(yy);
143: t1 = CCmplxSqrt(zeta);
144: XX[dim*i] = CCmplxRe(t1);
145: XX[dim*i+1] = CCmplxIm(t1);
146: } else if (idx==3) {
147: CCmplx zeta,t1,t2;
149: xx = XX[dim*i] - 0.8;
150: yy = XX[dim*i+1] + 1.5;
152: zeta.real = PetscRealPart(xx);
153: zeta.imag = PetscRealPart(yy);
154: t1 = CCmplxPow(zeta,-1.0);
155: t2 = CCmplxAdd(zeta,t1);
156: XX[dim*i] = CCmplxRe(t2);
157: XX[dim*i+1] = CCmplxIm(t2);
159: xx = XX[dim*i];
160: yy = XX[dim*i+1];
161: zeta.real = PetscRealPart(xx);
162: zeta.imag = PetscRealPart(yy);
163: t1 = CCmplxExp(zeta);
164: XX[dim*i] = CCmplxRe(t1);
165: XX[dim*i+1] = CCmplxIm(t1);
167: xx = XX[dim*i] + 0.4;
168: yy = XX[dim*i+1];
169: zeta.real = PetscRealPart(xx);
170: zeta.imag = PetscRealPart(yy);
171: t1 = CCmplxPow(zeta,2.0);
172: XX[dim*i] = CCmplxRe(t1);
173: XX[dim*i+1] = CCmplxIm(t1);
174: } else if (idx==4) {
175: PetscScalar Ni[4];
176: PetscScalar xi = XX[dim*i];
177: PetscScalar eta = XX[dim*i+1];
178: PetscScalar xn[] = {0.0,2.0,0.2,3.5};
179: PetscScalar yn[] = {-1.3,0.0,2.0,4.0};
180: PetscInt p;
182: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
183: Ni[1] = 0.25*(1.0+xi)*(1.0-eta);
184: Ni[2] = 0.25*(1.0-xi)*(1.0+eta);
185: Ni[3] = 0.25*(1.0+xi)*(1.0+eta);
187: xx = yy = 0.0;
188: for (p=0; p<4; p++) {
189: xx += Ni[p]*xn[p];
190: yy += Ni[p]*yn[p];
191: }
192: XX[dim*i] = xx;
193: XX[dim*i+1] = yy;
194: }
195: }
196: VecRestoreArray(Gcoords,&XX);
197: return(0);
198: }
202: PetscErrorCode DAApplyTrilinearMapping(DM da)
203: {
205: PetscInt i,j,k;
206: PetscInt sx,nx,sy,ny,sz,nz;
207: Vec Gcoords;
208: DMDACoor3d ***XX;
209: PetscScalar xx,yy,zz;
210: DM cda;
213: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
214: DMGetCoordinateDM(da,&cda);
215: DMGetCoordinates(da,&Gcoords);
217: DMDAVecGetArray(cda,Gcoords,&XX);
218: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
220: for (i=sx; i<sx+nx; i++) {
221: for (j=sy; j<sy+ny; j++) {
222: for (k=sz; k<sz+nz; k++) {
223: PetscScalar Ni[8];
224: PetscScalar xi = XX[k][j][i].x;
225: PetscScalar eta = XX[k][j][i].y;
226: PetscScalar zeta = XX[k][j][i].z;
227: PetscScalar xn[] = {0.0,2.0,0.2,3.5, 0.0,2.1,0.23,3.125 };
228: PetscScalar yn[] = {-1.3,0.0,2.0,4.0, -1.45,-0.1,2.24,3.79 };
229: PetscScalar zn[] = {0.0,0.3,-0.1,0.123, 0.956,1.32,1.12,0.798 };
230: PetscInt p;
232: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
233: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
234: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
235: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
237: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
238: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
239: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
240: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
242: xx = yy = zz = 0.0;
243: for (p=0; p<8; p++) {
244: xx += Ni[p]*xn[p];
245: yy += Ni[p]*yn[p];
246: zz += Ni[p]*zn[p];
247: }
248: XX[k][j][i].x = xx;
249: XX[k][j][i].y = yy;
250: XX[k][j][i].z = zz;
251: }
252: }
253: }
254: DMDAVecRestoreArray(cda,Gcoords,&XX);
255: return(0);
256: }
260: PetscErrorCode DADefineXLinearField2D(DM da,Vec field)
261: {
263: PetscInt i,j;
264: PetscInt sx,nx,sy,ny;
265: Vec Gcoords;
266: DMDACoor2d **XX;
267: PetscScalar **FF;
268: DM cda;
271: DMGetCoordinateDM(da,&cda);
272: DMGetCoordinates(da,&Gcoords);
274: DMDAVecGetArray(cda,Gcoords,&XX);
275: DMDAVecGetArray(da,field,&FF);
277: DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);
279: for (i=sx; i<sx+nx; i++) {
280: for (j=sy; j<sy+ny; j++) {
281: 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;
282: }
283: }
285: DMDAVecRestoreArray(da,field,&FF);
286: DMDAVecRestoreArray(cda,Gcoords,&XX);
287: return(0);
288: }
292: PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
293: {
295: PetscInt i,j,k;
296: PetscInt sx,nx,sy,ny,sz,nz;
297: Vec Gcoords;
298: DMDACoor3d ***XX;
299: PetscScalar ***FF;
300: DM cda;
303: DMGetCoordinateDM(da,&cda);
304: DMGetCoordinates(da,&Gcoords);
306: DMDAVecGetArray(cda,Gcoords,&XX);
307: DMDAVecGetArray(da,field,&FF);
309: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
311: for (k=sz; k<sz+nz; k++) {
312: for (j=sy; j<sy+ny; j++) {
313: for (i=sx; i<sx+nx; i++) {
314: FF[k][j][i] = 10.0
315: + 4.05 * XX[k][j][i].x
316: + 5.50 * XX[k][j][i].y
317: + 1.33 * XX[k][j][i].z
318: + 2.03 * XX[k][j][i].x * XX[k][j][i].y
319: + 0.03 * XX[k][j][i].x * XX[k][j][i].z
320: + 0.83 * XX[k][j][i].y * XX[k][j][i].z
321: + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
322: }
323: }
324: }
326: DMDAVecRestoreArray(da,field,&FF);
327: DMDAVecRestoreArray(cda,Gcoords,&XX);
328: return(0);
329: }
333: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
334: {
336: DM dac,daf;
337: PetscViewer vv;
338: Vec ac,af;
339: PetscInt Mx;
340: Mat II,INTERP;
341: Vec scale;
342: PetscBool output = PETSC_FALSE;
345: DMDACreate1d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE,
346: mx+1,
347: 1, /* 1 dof */
348: 1, /* stencil = 1 */
349: NULL,
350: &dac);
351: DMSetFromOptions(dac);
353: DMRefine(dac,MPI_COMM_NULL,&daf);
354: DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
355: Mx--;
357: DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
358: DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
360: {
361: DM cdaf,cdac;
362: Vec coordsc,coordsf;
364: DMGetCoordinateDM(dac,&cdac);
365: DMGetCoordinateDM(daf,&cdaf);
367: DMGetCoordinates(dac,&coordsc);
368: DMGetCoordinates(daf,&coordsf);
370: DMCreateInterpolation(cdac,cdaf,&II,&scale);
371: MatInterpolate(II,coordsc,coordsf);
372: MatDestroy(&II);
373: VecDestroy(&scale);
374: }
376: DMCreateInterpolation(dac,daf,&INTERP,NULL);
378: DMCreateGlobalVector(dac,&ac);
379: VecSet(ac,66.99);
381: DMCreateGlobalVector(daf,&af);
383: MatMult(INTERP,ac, af);
385: {
386: Vec afexact;
387: PetscReal nrm;
388: PetscInt N;
390: DMCreateGlobalVector(daf,&afexact);
391: VecSet(afexact,66.99);
392: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
393: VecNorm(afexact,NORM_2,&nrm);
394: VecGetSize(afexact,&N);
395: PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,nrm/sqrt((PetscReal)N));
396: VecDestroy(&afexact);
397: }
399: PetscOptionsGetBool(NULL,"-output",&output,NULL);
400: if (output) {
401: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);
402: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
403: DMView(dac, vv);
404: VecView(ac, vv);
405: PetscViewerDestroy(&vv);
407: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);
408: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
409: DMView(daf, vv);
410: VecView(af, vv);
411: PetscViewerDestroy(&vv);
412: }
414: MatDestroy(&INTERP);
415: DMDestroy(&dac);
416: DMDestroy(&daf);
417: VecDestroy(&ac);
418: VecDestroy(&af);
419: return(0);
420: }
424: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
425: {
427: DM dac,daf;
428: PetscViewer vv;
429: Vec ac,af;
430: PetscInt map_id,Mx,My;
431: Mat II,INTERP;
432: Vec scale;
433: PetscBool output = PETSC_FALSE;
436: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,
437: mx+1, my+1,
438: PETSC_DECIDE, PETSC_DECIDE,
439: 1, /* 1 dof */
440: 1, /* stencil = 1 */
441: NULL, NULL,
442: &dac);
443: DMSetFromOptions(dac);
445: DMRefine(dac,MPI_COMM_NULL,&daf);
446: DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
447: Mx--; My--;
449: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
450: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
452: /* apply conformal mappings */
453: map_id = 0;
454: PetscOptionsGetInt(NULL,"-cmap", &map_id,NULL);
455: if (map_id >= 1) {
456: DAApplyConformalMapping(dac,map_id);
457: }
459: {
460: DM cdaf,cdac;
461: Vec coordsc,coordsf;
463: DMGetCoordinateDM(dac,&cdac);
464: DMGetCoordinateDM(daf,&cdaf);
466: DMGetCoordinates(dac,&coordsc);
467: DMGetCoordinates(daf,&coordsf);
469: DMCreateInterpolation(cdac,cdaf,&II,&scale);
470: MatInterpolate(II,coordsc,coordsf);
471: MatDestroy(&II);
472: VecDestroy(&scale);
473: }
476: DMCreateInterpolation(dac,daf,&INTERP,NULL);
478: DMCreateGlobalVector(dac,&ac);
479: DADefineXLinearField2D(dac,ac);
481: DMCreateGlobalVector(daf,&af);
482: MatMult(INTERP,ac, af);
484: {
485: Vec afexact;
486: PetscReal nrm;
487: PetscInt N;
489: DMCreateGlobalVector(daf,&afexact);
490: VecZeroEntries(afexact);
491: DADefineXLinearField2D(daf,afexact);
492: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
493: VecNorm(afexact,NORM_2,&nrm);
494: VecGetSize(afexact,&N);
495: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,nrm/sqrt((PetscReal)N));
496: VecDestroy(&afexact);
497: }
499: PetscOptionsGetBool(NULL,"-output",&output,NULL);
500: if (output) {
501: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);
502: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
503: DMView(dac, vv);
504: VecView(ac, vv);
505: PetscViewerDestroy(&vv);
507: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);
508: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
509: DMView(daf, vv);
510: VecView(af, vv);
511: PetscViewerDestroy(&vv);
512: }
514: MatDestroy(&INTERP);
515: DMDestroy(&dac);
516: DMDestroy(&daf);
517: VecDestroy(&ac);
518: VecDestroy(&af);
519: return(0);
520: }
524: PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
525: {
527: DM dac,daf;
528: PetscViewer vv;
529: Vec ac,af;
530: PetscInt map_id,Mx,My,Mz;
531: Mat II,INTERP;
532: Vec scale;
533: PetscBool output = PETSC_FALSE;
536: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
537: mx+1, my+1,mz+1,
538: PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
539: 1, /* 1 dof */
540: 1, /* stencil = 1 */
541: NULL,NULL,NULL,
542: &dac);
543: DMSetFromOptions(dac);
545: DMRefine(dac,MPI_COMM_NULL,&daf);
546: DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);
547: Mx--; My--; Mz--;
549: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
550: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);
552: /* apply trilinear mappings */
553: /*DAApplyTrilinearMapping(dac);*/
554: /* apply conformal mappings */
555: map_id = 0;
556: PetscOptionsGetInt(NULL,"-cmap", &map_id,NULL);
557: if (map_id >= 1) {
558: DAApplyConformalMapping(dac,map_id);
559: }
561: {
562: DM cdaf,cdac;
563: Vec coordsc,coordsf;
565: DMGetCoordinateDM(dac,&cdac);
566: DMGetCoordinateDM(daf,&cdaf);
568: DMGetCoordinates(dac,&coordsc);
569: DMGetCoordinates(daf,&coordsf);
571: DMCreateInterpolation(cdac,cdaf,&II,&scale);
572: MatInterpolate(II,coordsc,coordsf);
573: MatDestroy(&II);
574: VecDestroy(&scale);
575: }
577: DMCreateInterpolation(dac,daf,&INTERP,NULL);
579: DMCreateGlobalVector(dac,&ac);
580: VecZeroEntries(ac);
581: DADefineXLinearField3D(dac,ac);
583: DMCreateGlobalVector(daf,&af);
584: VecZeroEntries(af);
586: MatMult(INTERP,ac, af);
588: {
589: Vec afexact;
590: PetscReal nrm;
591: PetscInt N;
593: DMCreateGlobalVector(daf,&afexact);
594: VecZeroEntries(afexact);
595: DADefineXLinearField3D(daf,afexact);
596: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
597: VecNorm(afexact,NORM_2,&nrm);
598: VecGetSize(afexact,&N);
599: 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,nrm/sqrt((PetscReal)N));
600: VecDestroy(&afexact);
601: }
603: PetscOptionsGetBool(NULL,"-output",&output,NULL);
604: if (output) {
605: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);
606: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
607: DMView(dac, vv);
608: VecView(ac, vv);
609: PetscViewerDestroy(&vv);
611: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);
612: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
613: DMView(daf, vv);
614: VecView(af, vv);
615: PetscViewerDestroy(&vv);
616: }
618: MatDestroy(&INTERP);
619: DMDestroy(&dac);
620: DMDestroy(&daf);
621: VecDestroy(&ac);
622: VecDestroy(&af);
623: return(0);
624: }
628: int main(int argc,char **argv)
629: {
631: PetscInt mx,my,mz,l,nl,dim;
633: PetscInitialize(&argc,&argv,0,help);
635: mx = my = mz = 2;
636: PetscOptionsGetInt(NULL,"-mx", &mx, 0);
637: PetscOptionsGetInt(NULL,"-my", &my, 0);
638: PetscOptionsGetInt(NULL,"-mz", &mz, 0);
639: nl = 1;
640: PetscOptionsGetInt(NULL,"-nl", &nl, 0);
641: dim = 2;
642: PetscOptionsGetInt(NULL,"-dim", &dim, 0);
644: for (l=0; l<nl; l++) {
645: if (dim==1) da_test_RefineCoords1D(mx);
646: else if (dim==2) da_test_RefineCoords2D(mx,my);
647: else if (dim==3) da_test_RefineCoords3D(mx,my,mz);
648: mx = mx * 2;
649: my = my * 2;
650: mz = mz * 2;
651: }
653: PetscFinalize();
654: return 0;
655: }