Actual source code: ex36.c
petsc-3.3-p6 2013-02-11
2: static char help[] = "Checks the functionality of DMDAGetInterpolation on deformed grids.\n\n";
4: #include <petsc.h>
5: #include <petscvec.h>
6: #include <petscmat.h>
7: #include <petscdmda.h>
9: typedef struct _n_CCmplx CCmplx;
10: struct _n_CCmplx {
11: double real;
12: double imag;
13: };
15: CCmplx CCmplxPow(CCmplx a,double n)
16: {
17: CCmplx b;
18: double r,theta;
19: r = sqrt(a.real*a.real + a.imag*a.imag);
20: theta = atan2(a.imag,a.real);
21: b.real = pow(r,n) * cos(n*theta);
22: b.imag = pow(r,n) * sin(n*theta);
23: return b;
24: }
25: CCmplx CCmplxExp(CCmplx a)
26: {
27: CCmplx b;
28: b.real = exp(a.real) * cos(a.imag);
29: b.imag = exp(a.real) * sin(a.imag);
30: return b;
31: }
32: CCmplx CCmplxSqrt(CCmplx a)
33: {
34: CCmplx b;
35: double r,theta;
36: r = sqrt(a.real*a.real + a.imag*a.imag);
37: theta = atan2(a.imag,a.real);
38: b.real = sqrt(r) * cos(0.5*theta);
39: b.imag = sqrt(r) * sin(0.5*theta);
40: return b;
41: }
42: CCmplx CCmplxAdd(CCmplx a,CCmplx c)
43: {
44: CCmplx b;
45: b.real = a.real +c.real;
46: b.imag = a.imag +c.imag;
47: return b;
48: }
49: PetscScalar CCmplxRe(CCmplx a)
50: {
51: return (PetscScalar)a.real;
52: }
53: PetscScalar CCmplxIm(CCmplx a)
54: {
55: return (PetscScalar)a.imag;
56: }
60: PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx)
61: {
63: PetscInt i,n;
64: PetscInt sx,nx,sy,ny,sz,nz,dim;
65: Vec Gcoords;
66: PetscScalar *XX;
67: PetscScalar xx,yy,zz;
68: DM cda;
69:
71: if (idx==0) {
72: return(0);
73: } else if (idx==1) { /* dam break */
74: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
75: } else if (idx==2) { /* stagnation in a corner */
76: DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0 );
77: } else if (idx==3) { /* nautilis */
78: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
79: } else if (idx==4) {
80: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
81: }
82:
83: DMDAGetCoordinateDA(da,&cda);
84: DMDAGetCoordinates(da,&Gcoords);
85:
86: VecGetArray(Gcoords,&XX);
87: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
88: DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
89: VecGetLocalSize(Gcoords,&n);
90: n = n / dim;
91:
92: for (i=0; i<n; i++) {
93: if ( (dim==3) && (idx!=2) ) {
94: PetscScalar Ni[8];
95: PetscScalar xi = XX[dim*i ];
96: PetscScalar eta = XX[dim*i+1];
97: PetscScalar zeta = XX[dim*i+2];
98: PetscScalar xn[] = {-1.0,1.0,-1.0,1.0, -1.0,1.0,-1.0,1.0 };
99: PetscScalar yn[] = {-1.0,-1.0,1.0,1.0, -1.0,-1.0,1.0,1.0 };
100: PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0, 0.1,4.0,0.2,1.0 };
101: PetscInt p;
102:
103: Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
104: Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
105: Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
106: Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
107:
108: Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
109: Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
110: Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
111: Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
112: xx = yy = zz = 0.0;
113: for (p=0; p<8; p++ ) {
114: xx += Ni[p]*xn[p];
115: yy += Ni[p]*yn[p];
116: zz += Ni[p]*zn[p];
117: }
118: XX[dim*i ] = xx;
119: XX[dim*i+1] = yy;
120: XX[dim*i+2] = zz;
121: }
122:
123: if (idx==1) {
124: CCmplx zeta,t1,t2;
125:
126: xx = XX[dim*i ] - 0.8;
127: yy = XX[dim*i+1] + 1.5;
128: zeta.real = PetscRealPart(xx);
129: zeta.imag = PetscRealPart(yy);
130:
131: t1 = CCmplxPow(zeta,-1.0);
132: 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;
137:
138: xx = XX[dim*i ];
139: yy = XX[dim*i+1];
140: zeta.real = PetscRealPart(xx);
141: zeta.imag = PetscRealPart(yy);
142:
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;
148:
149: xx = XX[dim*i ] - 0.8;
150: yy = XX[dim*i+1] + 1.5;
151:
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);
158:
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);
166:
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: }
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;
182:
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);
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:
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;
212:
214: DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0 );
215: DMDAGetCoordinateDA(da,&cda);
216: DMDAGetCoordinates(da,&Gcoords);
217:
218: DMDAVecGetArray(cda,Gcoords,&XX);
219: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
220:
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;
232:
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);
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:
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;
270:
272: DMDAGetCoordinateDA(da,&cda);
273: DMDAGetCoordinates(da,&Gcoords);
274:
275: DMDAVecGetArray(cda,Gcoords,&XX);
276: DMDAVecGetArray(da,field,&FF);
277:
278: DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);
279:
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: }
285:
286: DMDAVecRestoreArray(da,field,&FF);
287: DMDAVecRestoreArray(cda,Gcoords,&XX);
288:
289: return(0);
290: }
294: PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
295: {
297: PetscInt i,j,k;
298: PetscInt sx,nx,sy,ny,sz,nz;
299: Vec Gcoords;
300: DMDACoor3d ***XX;
301: PetscScalar ***FF;
302: DM cda;
303:
305: DMDAGetCoordinateDA(da,&cda);
306: DMDAGetCoordinates(da,&Gcoords);
307:
308: DMDAVecGetArray(cda,Gcoords,&XX);
309: DMDAVecGetArray(da,field,&FF);
310:
311: DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
312:
313: for (k=sz; k<sz+nz; k++) {
314: for (j=sy; j<sy+ny; j++ ) {
315: for (i=sx; i<sx+nx; i++) {
316: FF[k][j][i] = 10.0
317: + 4.05 * XX[k][j][i].x
318: + 5.50 * XX[k][j][i].y
319: + 1.33 * XX[k][j][i].z
320: + 2.03 * XX[k][j][i].x * XX[k][j][i].y
321: + 0.03 * XX[k][j][i].x * XX[k][j][i].z
322: + 0.83 * XX[k][j][i].y * XX[k][j][i].z
323: + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
324: }
325: }
326: }
327:
328: DMDAVecRestoreArray(da,field,&FF);
329: DMDAVecRestoreArray(cda,Gcoords,&XX);
330:
331: return(0);
332: }
336: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
337: {
339: DM dac,daf;
340: PetscViewer vv;
341: Vec ac,af;
342: PetscInt Mx;
343: Mat II,INTERP;
344: Vec scale;
345: PetscBool output = PETSC_FALSE;
346:
348: DMDACreate1d( PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE,
349: mx+1,
350: 1, /* 1 dof */
351: 1, /* stencil = 1 */
352: PETSC_NULL,
353: &dac );
354: DMSetFromOptions(dac);
355:
356: DMRefine(dac,MPI_COMM_NULL,&daf);
357: DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
358: Mx--;
359:
360: DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE );
361: DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE );
362:
363: {
364: DM cdaf,cdac;
365: Vec coordsc,coordsf;
366:
367: DMDAGetCoordinateDA(dac,&cdac);
368: DMDAGetCoordinateDA(daf,&cdaf);
369:
370: DMDAGetCoordinates(dac,&coordsc);
371: DMDAGetCoordinates(daf,&coordsf);
372:
373: DMCreateInterpolation(cdac,cdaf,&II,&scale);
374: MatInterpolate(II,coordsc,coordsf);
375: MatDestroy(&II);
376: VecDestroy(&scale);
377: }
378:
379: DMCreateInterpolation(dac,daf,&INTERP,PETSC_NULL);
380:
381: DMCreateGlobalVector(dac,&ac);
382: VecSet(ac,66.99);
383:
384: DMCreateGlobalVector(daf,&af);
385:
386: MatMult(INTERP,ac, af);
387:
388: {
389: Vec afexact;
390: PetscReal nrm;
391: PetscInt N;
392:
393: DMCreateGlobalVector(daf,&afexact);
394: VecSet(afexact,66.99);
395: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
396: VecNorm(afexact,NORM_2,&nrm);
397: VecGetSize(afexact,&N);
398: PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,nrm/sqrt((PetscReal)N) );
399: VecDestroy(&afexact);
400: }
401:
402: PetscOptionsGetBool(PETSC_NULL,"-output",&output,PETSC_NULL);
403: if (output) {
404: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);
405: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
406: DMView(dac, vv);
407: VecView(ac, vv);
408: PetscViewerDestroy(&vv);
409:
410: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);
411: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
412: DMView(daf, vv);
413: VecView(af, vv);
414: PetscViewerDestroy(&vv);
415: }
416:
417: MatDestroy(&INTERP);
418: DMDestroy(&dac);
419: DMDestroy(&daf);
420: VecDestroy(&ac);
421: VecDestroy(&af);
422:
423: return(0);
424: }
428: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
429: {
431: DM dac,daf;
432: PetscViewer vv;
433: Vec ac,af;
434: PetscInt map_id,Mx,My;
435: Mat II,INTERP;
436: Vec scale;
437: PetscBool output = PETSC_FALSE;
438:
440: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_STENCIL_BOX,
441: mx+1, my+1,
442: PETSC_DECIDE, PETSC_DECIDE,
443: 1, /* 1 dof */
444: 1, /* stencil = 1 */
445: PETSC_NULL, PETSC_NULL,
446: &dac );
447: DMSetFromOptions(dac);
448:
449: DMRefine(dac,MPI_COMM_NULL,&daf);
450: DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
451: Mx--; My--;
452:
453: DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE );
454: DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE );
455:
456: /* apply conformal mappings */
457: map_id = 0;
458: PetscOptionsGetInt( PETSC_NULL,"-cmap", &map_id,PETSC_NULL );
459: if( map_id >= 1 ) {
460: DAApplyConformalMapping(dac,map_id);
461: }
462:
463: {
464: DM cdaf,cdac;
465: Vec coordsc,coordsf;
466:
467: DMDAGetCoordinateDA(dac,&cdac);
468: DMDAGetCoordinateDA(daf,&cdaf);
469:
470: DMDAGetCoordinates(dac,&coordsc);
471: DMDAGetCoordinates(daf,&coordsf);
472:
473: DMCreateInterpolation(cdac,cdaf,&II,&scale);
474: MatInterpolate(II,coordsc,coordsf);
475: MatDestroy(&II);
476: VecDestroy(&scale);
477: }
478:
479:
480: DMCreateInterpolation(dac,daf,&INTERP,PETSC_NULL);
481:
482: DMCreateGlobalVector(dac,&ac);
483: DADefineXLinearField2D(dac,ac);
484:
485: DMCreateGlobalVector(daf,&af);
486: MatMult(INTERP,ac, af);
487:
488: {
489: Vec afexact;
490: PetscReal nrm;
491: PetscInt N;
492:
493: DMCreateGlobalVector(daf,&afexact);
494: VecZeroEntries(afexact);
495: DADefineXLinearField2D(daf,afexact);
496: VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
497: VecNorm(afexact,NORM_2,&nrm);
498: VecGetSize(afexact,&N);
499: PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,nrm/sqrt((PetscReal)N) );
500: VecDestroy(&afexact);
501: }
502:
503: PetscOptionsGetBool(PETSC_NULL,"-output",&output,PETSC_NULL);
504: if (output) {
505: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);
506: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
507: DMView(dac, vv);
508: VecView(ac, vv);
509: PetscViewerDestroy(&vv);
510:
511: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);
512: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
513: DMView(daf, vv);
514: VecView(af, vv);
515: PetscViewerDestroy(&vv);
516: }
517:
518: MatDestroy(&INTERP);
519: DMDestroy(&dac);
520: DMDestroy(&daf);
521: VecDestroy(&ac);
522: VecDestroy(&af);
523:
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;
539:
541: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_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: PETSC_NULL,PETSC_NULL,PETSC_NULL,
547: &dac );
548: DMSetFromOptions(dac);
549:
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--;
553:
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 );
556:
557: /* apply trilinear mappings */
558: /*DAApplyTrilinearMapping(dac);*/
559: /* apply conformal mappings */
560: map_id = 0;
561: PetscOptionsGetInt( PETSC_NULL,"-cmap", &map_id,PETSC_NULL );
562: if( map_id >= 1 ) {
563: DAApplyConformalMapping(dac,map_id);
564: }
565:
566: {
567: DM cdaf,cdac;
568: Vec coordsc,coordsf;
569:
570: DMDAGetCoordinateDA(dac,&cdac);
571: DMDAGetCoordinateDA(daf,&cdaf);
572:
573: DMDAGetCoordinates(dac,&coordsc);
574: DMDAGetCoordinates(daf,&coordsf);
575:
576: DMCreateInterpolation(cdac,cdaf,&II,&scale);
577: MatInterpolate(II,coordsc,coordsf);
578: MatDestroy(&II);
579: VecDestroy(&scale);
580: }
581:
582: DMCreateInterpolation(dac,daf,&INTERP,PETSC_NULL);
583:
584: DMCreateGlobalVector(dac,&ac);
585: VecZeroEntries(ac);
586: DADefineXLinearField3D(dac,ac);
587:
588: DMCreateGlobalVector(daf,&af);
589: VecZeroEntries(af);
590:
591: MatMult(INTERP,ac, af);
592:
593: {
594: Vec afexact;
595: PetscReal nrm;
596: PetscInt N;
597:
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,nrm/sqrt((PetscReal)N) );
605: VecDestroy(&afexact);
606: }
607:
608: PetscOptionsGetBool(PETSC_NULL,"-output",&output,PETSC_NULL);
609: if (output) {
610: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);
611: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
612: DMView(dac, vv);
613: VecView(ac, vv);
614: PetscViewerDestroy(&vv);
615:
616: PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);
617: PetscViewerSetFormat(vv, PETSC_VIEWER_ASCII_VTK);
618: DMView(daf, vv);
619: VecView(af, vv);
620: PetscViewerDestroy(&vv);
621: }
622:
623: MatDestroy(&INTERP);
624: DMDestroy(&dac);
625: DMDestroy(&daf);
626: VecDestroy(&ac);
627: VecDestroy(&af);
628:
629: return(0);
630: }
634: int main(int argc,char **argv)
635: {
637: PetscInt mx,my,mz,l,nl,dim;
638:
639: PetscInitialize(&argc,&argv,0,help);
641: mx = my = mz = 2;
642: PetscOptionsGetInt(PETSC_NULL,"-mx", &mx, 0 );
643: PetscOptionsGetInt(PETSC_NULL,"-my", &my, 0 );
644: PetscOptionsGetInt(PETSC_NULL,"-mz", &mz, 0 );
645: nl = 1;
646: PetscOptionsGetInt(PETSC_NULL,"-nl", &nl, 0 );
647: dim = 2;
648: PetscOptionsGetInt(PETSC_NULL,"-dim", &dim, 0 );
649:
650: for (l=0; l<nl; l++) {
651: if (dim==1) { da_test_RefineCoords1D(mx); }
652: else if (dim==2) { da_test_RefineCoords2D(mx,my); }
653: else if (dim==3) { da_test_RefineCoords3D(mx,my,mz); }
654: mx = mx * 2;
655: my = my * 2;
656: mz = mz * 2;
657: }
658:
659: PetscFinalize();
660: return 0;
661: }