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: }