Actual source code: ex36.c
petsc-3.5.4 2015-05-23
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: DMDAVecGetArray(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: DMDAVecRestoreArray(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: DMDAVecGetArray(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: DMDAVecRestoreArray(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: DMDAVecGetArray(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: DMDAVecRestoreArray(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,"-output",&output,NULL);
401: if (output) {
402: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);
403: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
404: DMView(dac, vv);
405: VecView(ac, vv);
406: PetscViewerDestroy(&vv);
408: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);
409: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
410: DMView(daf, vv);
411: VecView(af, vv);
412: PetscViewerDestroy(&vv);
413: }
415: MatDestroy(&INTERP);
416: DMDestroy(&dac);
417: DMDestroy(&daf);
418: VecDestroy(&ac);
419: VecDestroy(&af);
420: return(0);
421: }
425: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
426: {
428: DM dac,daf;
429: PetscViewer vv;
430: Vec ac,af;
431: PetscInt map_id,Mx,My;
432: Mat II,INTERP;
433: Vec scale;
434: PetscBool output = PETSC_FALSE;
437: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,
438: mx+1, my+1,
439: PETSC_DECIDE, PETSC_DECIDE,
440: 1, /* 1 dof */
441: 1, /* stencil = 1 */
442: NULL, NULL,
443: &dac);
444: DMSetFromOptions(dac);
446: DMRefine(dac,MPI_COMM_NULL,&daf);
447: DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
448: Mx--; My--;
450: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
451: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
453: /* apply conformal mappings */
454: map_id = 0;
455: PetscOptionsGetInt(NULL,"-cmap", &map_id,NULL);
456: if (map_id >= 1) {
457: DAApplyConformalMapping(dac,map_id);
458: }
460: {
461: DM cdaf,cdac;
462: Vec coordsc,coordsf;
464: DMGetCoordinateDM(dac,&cdac);
465: DMGetCoordinateDM(daf,&cdaf);
467: DMGetCoordinates(dac,&coordsc);
468: DMGetCoordinates(daf,&coordsf);
470: DMCreateInterpolation(cdac,cdaf,&II,&scale);
471: MatInterpolate(II,coordsc,coordsf);
472: MatDestroy(&II);
473: VecDestroy(&scale);
474: }
477: DMCreateInterpolation(dac,daf,&INTERP,NULL);
479: DMCreateGlobalVector(dac,&ac);
480: DADefineXLinearField2D(dac,ac);
482: DMCreateGlobalVector(daf,&af);
483: MatMult(INTERP,ac, af);
485: {
486: Vec afexact;
487: PetscReal nrm;
488: PetscInt N;
490: DMCreateGlobalVector(daf,&afexact);
491: VecZeroEntries(afexact);
492: DADefineXLinearField2D(daf,afexact);
493: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
494: VecNorm(afexact,NORM_2,&nrm);
495: VecGetSize(afexact,&N);
496: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)));
497: VecDestroy(&afexact);
498: }
500: PetscOptionsGetBool(NULL,"-output",&output,NULL);
501: if (output) {
502: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);
503: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
504: DMView(dac, vv);
505: VecView(ac, vv);
506: PetscViewerDestroy(&vv);
508: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);
509: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
510: DMView(daf, vv);
511: VecView(af, vv);
512: PetscViewerDestroy(&vv);
513: }
515: MatDestroy(&INTERP);
516: DMDestroy(&dac);
517: DMDestroy(&daf);
518: VecDestroy(&ac);
519: VecDestroy(&af);
520: return(0);
521: }
525: PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
526: {
528: DM dac,daf;
529: PetscViewer vv;
530: Vec ac,af;
531: PetscInt map_id,Mx,My,Mz;
532: Mat II,INTERP;
533: Vec scale;
534: PetscBool output = PETSC_FALSE;
537: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
538: mx+1, my+1,mz+1,
539: PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,
540: 1, /* 1 dof */
541: 1, /* stencil = 1 */
542: NULL,NULL,NULL,
543: &dac);
544: DMSetFromOptions(dac);
546: DMRefine(dac,MPI_COMM_NULL,&daf);
547: DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);
548: Mx--; My--; Mz--;
550: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
551: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);
553: /* apply trilinear mappings */
554: /*DAApplyTrilinearMapping(dac);*/
555: /* apply conformal mappings */
556: map_id = 0;
557: PetscOptionsGetInt(NULL,"-cmap", &map_id,NULL);
558: if (map_id >= 1) {
559: DAApplyConformalMapping(dac,map_id);
560: }
562: {
563: DM cdaf,cdac;
564: Vec coordsc,coordsf;
566: DMGetCoordinateDM(dac,&cdac);
567: DMGetCoordinateDM(daf,&cdaf);
569: DMGetCoordinates(dac,&coordsc);
570: DMGetCoordinates(daf,&coordsf);
572: DMCreateInterpolation(cdac,cdaf,&II,&scale);
573: MatInterpolate(II,coordsc,coordsf);
574: MatDestroy(&II);
575: VecDestroy(&scale);
576: }
578: DMCreateInterpolation(dac,daf,&INTERP,NULL);
580: DMCreateGlobalVector(dac,&ac);
581: VecZeroEntries(ac);
582: DADefineXLinearField3D(dac,ac);
584: DMCreateGlobalVector(daf,&af);
585: VecZeroEntries(af);
587: MatMult(INTERP,ac, af);
589: {
590: Vec afexact;
591: PetscReal nrm;
592: PetscInt N;
594: DMCreateGlobalVector(daf,&afexact);
595: VecZeroEntries(afexact);
596: DADefineXLinearField3D(daf,afexact);
597: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
598: VecNorm(afexact,NORM_2,&nrm);
599: VecGetSize(afexact,&N);
600: 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)));
601: VecDestroy(&afexact);
602: }
604: PetscOptionsGetBool(NULL,"-output",&output,NULL);
605: if (output) {
606: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);
607: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
608: DMView(dac, vv);
609: VecView(ac, vv);
610: PetscViewerDestroy(&vv);
612: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);
613: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
614: DMView(daf, vv);
615: VecView(af, vv);
616: PetscViewerDestroy(&vv);
617: }
619: MatDestroy(&INTERP);
620: DMDestroy(&dac);
621: DMDestroy(&daf);
622: VecDestroy(&ac);
623: VecDestroy(&af);
624: return(0);
625: }
629: int main(int argc,char **argv)
630: {
632: PetscInt mx,my,mz,l,nl,dim;
634: PetscInitialize(&argc,&argv,0,help);
636: mx = my = mz = 2;
637: PetscOptionsGetInt(NULL,"-mx", &mx, 0);
638: PetscOptionsGetInt(NULL,"-my", &my, 0);
639: PetscOptionsGetInt(NULL,"-mz", &mz, 0);
640: nl = 1;
641: PetscOptionsGetInt(NULL,"-nl", &nl, 0);
642: dim = 2;
643: PetscOptionsGetInt(NULL,"-dim", &dim, 0);
645: for (l=0; l<nl; l++) {
646: if (dim==1) da_test_RefineCoords1D(mx);
647: else if (dim==2) da_test_RefineCoords2D(mx,my);
648: else if (dim==3) da_test_RefineCoords3D(mx,my,mz);
649: mx = mx * 2;
650: my = my * 2;
651: mz = mz * 2;
652: }
654: PetscFinalize();
655: return 0;
656: }