Actual source code: ex36.c


  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  = PetscAtan2Real(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  = PetscAtan2Real(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: }

 56: PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx)
 57: {
 58:   PetscInt       i,n;
 59:   PetscInt       sx,nx,sy,ny,sz,nz,dim;
 60:   Vec            Gcoords;
 61:   PetscScalar    *XX;
 62:   PetscScalar    xx,yy,zz;
 63:   DM             cda;

 66:   if (idx==0) {
 67:     return 0;
 68:   } else if (idx==1) { /* dam break */
 69:     DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
 70:   } else if (idx==2) { /* stagnation in a corner */
 71:     DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0);
 72:   } else if (idx==3) { /* nautilis */
 73:     DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
 74:   } else if (idx==4) {
 75:     DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
 76:   }

 78:   DMGetCoordinateDM(da,&cda);
 79:   DMGetCoordinates(da,&Gcoords);

 81:   VecGetArray(Gcoords,&XX);
 82:   DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);
 83:   DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);
 84:   VecGetLocalSize(Gcoords,&n);
 85:   n    = n / dim;

 87:   for (i=0; i<n; i++) {
 88:     if ((dim==3) && (idx!=2)) {
 89:       PetscScalar Ni[8];
 90:       PetscScalar xi   = XX[dim*i];
 91:       PetscScalar eta  = XX[dim*i+1];
 92:       PetscScalar zeta = XX[dim*i+2];
 93:       PetscScalar xn[] = {-1.0,1.0,-1.0,1.0,   -1.0,1.0,-1.0,1.0  };
 94:       PetscScalar yn[] = {-1.0,-1.0,1.0,1.0,   -1.0,-1.0,1.0,1.0  };
 95:       PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0,  0.1,4.0,0.2,1.0  };
 96:       PetscInt    p;

 98:       Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
 99:       Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
100:       Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
101:       Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);

103:       Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
104:       Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
105:       Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
106:       Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);

108:       xx = yy = zz = 0.0;
109:       for (p=0; p<8; p++) {
110:         xx += Ni[p]*xn[p];
111:         yy += Ni[p]*yn[p];
112:         zz += Ni[p]*zn[p];
113:       }
114:       XX[dim*i]   = xx;
115:       XX[dim*i+1] = yy;
116:       XX[dim*i+2] = zz;
117:     }

119:     if (idx==1) {
120:       CCmplx zeta,t1,t2;

122:       xx = XX[dim*i]   - 0.8;
123:       yy = XX[dim*i+1] + 1.5;

125:       zeta.real = PetscRealPart(xx);
126:       zeta.imag = PetscRealPart(yy);

128:       t1 = CCmplxPow(zeta,-1.0);
129:       t2 = CCmplxAdd(zeta,t1);

131:       XX[dim*i]   = CCmplxRe(t2);
132:       XX[dim*i+1] = CCmplxIm(t2);
133:     } else if (idx==2) {
134:       CCmplx zeta,t1;

136:       xx = XX[dim*i];
137:       yy = XX[dim*i+1];
138:       zeta.real = PetscRealPart(xx);
139:       zeta.imag = PetscRealPart(yy);

141:       t1 = CCmplxSqrt(zeta);
142:       XX[dim*i]   = CCmplxRe(t1);
143:       XX[dim*i+1] = CCmplxIm(t1);
144:     } else if (idx==3) {
145:       CCmplx zeta,t1,t2;

147:       xx = XX[dim*i]   - 0.8;
148:       yy = XX[dim*i+1] + 1.5;

150:       zeta.real   = PetscRealPart(xx);
151:       zeta.imag   = PetscRealPart(yy);
152:       t1          = CCmplxPow(zeta,-1.0);
153:       t2          = CCmplxAdd(zeta,t1);
154:       XX[dim*i]   = CCmplxRe(t2);
155:       XX[dim*i+1] = CCmplxIm(t2);

157:       xx          = XX[dim*i];
158:       yy          = XX[dim*i+1];
159:       zeta.real   = PetscRealPart(xx);
160:       zeta.imag   = PetscRealPart(yy);
161:       t1          = CCmplxExp(zeta);
162:       XX[dim*i]   = CCmplxRe(t1);
163:       XX[dim*i+1] = CCmplxIm(t1);

165:       xx          = XX[dim*i] + 0.4;
166:       yy          = XX[dim*i+1];
167:       zeta.real   = PetscRealPart(xx);
168:       zeta.imag   = PetscRealPart(yy);
169:       t1          = CCmplxPow(zeta,2.0);
170:       XX[dim*i]   = CCmplxRe(t1);
171:       XX[dim*i+1] = CCmplxIm(t1);
172:     } else if (idx==4) {
173:       PetscScalar Ni[4];
174:       PetscScalar xi   = XX[dim*i];
175:       PetscScalar eta  = XX[dim*i+1];
176:       PetscScalar xn[] = {0.0,2.0,0.2,3.5};
177:       PetscScalar yn[] = {-1.3,0.0,2.0,4.0};
178:       PetscInt    p;

180:       Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
181:       Ni[1] = 0.25*(1.0+xi)*(1.0-eta);
182:       Ni[2] = 0.25*(1.0-xi)*(1.0+eta);
183:       Ni[3] = 0.25*(1.0+xi)*(1.0+eta);

185:       xx = yy = 0.0;
186:       for (p=0; p<4; p++) {
187:         xx += Ni[p]*xn[p];
188:         yy += Ni[p]*yn[p];
189:       }
190:       XX[dim*i]   = xx;
191:       XX[dim*i+1] = yy;
192:     }
193:   }
194:   VecRestoreArray(Gcoords,&XX);
195:   return 0;
196: }

198: PetscErrorCode DAApplyTrilinearMapping(DM da)
199: {
200:   PetscInt       i,j,k;
201:   PetscInt       sx,nx,sy,ny,sz,nz;
202:   Vec            Gcoords;
203:   DMDACoor3d     ***XX;
204:   PetscScalar    xx,yy,zz;
205:   DM             cda;

208:   DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);
209:   DMGetCoordinateDM(da,&cda);
210:   DMGetCoordinates(da,&Gcoords);

212:   DMDAVecGetArrayRead(cda,Gcoords,&XX);
213:   DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);

215:   for (i=sx; i<sx+nx; i++) {
216:     for (j=sy; j<sy+ny; j++) {
217:       for (k=sz; k<sz+nz; k++) {
218:         PetscScalar Ni[8];
219:         PetscScalar xi   = XX[k][j][i].x;
220:         PetscScalar eta  = XX[k][j][i].y;
221:         PetscScalar zeta = XX[k][j][i].z;
222:         PetscScalar xn[] = {0.0,2.0,0.2,3.5,   0.0,2.1,0.23,3.125  };
223:         PetscScalar yn[] = {-1.3,0.0,2.0,4.0,  -1.45,-0.1,2.24,3.79  };
224:         PetscScalar zn[] = {0.0,0.3,-0.1,0.123,  0.956,1.32,1.12,0.798  };
225:         PetscInt    p;

227:         Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
228:         Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
229:         Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
230:         Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);

232:         Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
233:         Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
234:         Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
235:         Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);

237:         xx = yy = zz = 0.0;
238:         for (p=0; p<8; p++) {
239:           xx += Ni[p]*xn[p];
240:           yy += Ni[p]*yn[p];
241:           zz += Ni[p]*zn[p];
242:         }
243:         XX[k][j][i].x = xx;
244:         XX[k][j][i].y = yy;
245:         XX[k][j][i].z = zz;
246:       }
247:     }
248:   }
249:   DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
250:   return 0;
251: }

253: PetscErrorCode DADefineXLinearField2D(DM da,Vec field)
254: {
255:   PetscInt       i,j;
256:   PetscInt       sx,nx,sy,ny;
257:   Vec            Gcoords;
258:   DMDACoor2d     **XX;
259:   PetscScalar    **FF;
260:   DM             cda;

263:   DMGetCoordinateDM(da,&cda);
264:   DMGetCoordinates(da,&Gcoords);

266:   DMDAVecGetArrayRead(cda,Gcoords,&XX);
267:   DMDAVecGetArray(da,field,&FF);

269:   DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);

271:   for (i=sx; i<sx+nx; i++) {
272:     for (j=sy; j<sy+ny; j++) {
273:       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;
274:     }
275:   }

277:   DMDAVecRestoreArray(da,field,&FF);
278:   DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
279:   return 0;
280: }

282: PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
283: {
284:   PetscInt       i,j,k;
285:   PetscInt       sx,nx,sy,ny,sz,nz;
286:   Vec            Gcoords;
287:   DMDACoor3d     ***XX;
288:   PetscScalar    ***FF;
289:   DM             cda;

292:   DMGetCoordinateDM(da,&cda);
293:   DMGetCoordinates(da,&Gcoords);

295:   DMDAVecGetArrayRead(cda,Gcoords,&XX);
296:   DMDAVecGetArray(da,field,&FF);

298:   DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);

300:   for (k=sz; k<sz+nz; k++) {
301:     for (j=sy; j<sy+ny; j++) {
302:       for (i=sx; i<sx+nx; i++) {
303:         FF[k][j][i] = 10.0
304:                 + 4.05 * XX[k][j][i].x
305:                 + 5.50 * XX[k][j][i].y
306:                 + 1.33 * XX[k][j][i].z
307:                 + 2.03 * XX[k][j][i].x * XX[k][j][i].y
308:                 + 0.03 * XX[k][j][i].x * XX[k][j][i].z
309:                 + 0.83 * XX[k][j][i].y * XX[k][j][i].z
310:                 + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
311:       }
312:     }
313:   }

315:   DMDAVecRestoreArray(da,field,&FF);
316:   DMDAVecRestoreArrayRead(cda,Gcoords,&XX);
317:   return 0;
318: }

320: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
321: {
322:   DM             dac,daf;
323:   PetscViewer    vv;
324:   Vec            ac,af;
325:   PetscInt       Mx;
326:   Mat            II,INTERP;
327:   Vec            scale;
328:   PetscBool      output = PETSC_FALSE;

331:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,mx+1,1, /* 1 dof */1, /* stencil = 1 */NULL,&dac);
332:   DMSetFromOptions(dac);
333:   DMSetUp(dac);

335:   DMRefine(dac,MPI_COMM_NULL,&daf);
336:   DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);
337:   Mx--;

339:   DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);
340:   DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);

342:   {
343:     DM  cdaf,cdac;
344:     Vec coordsc,coordsf;

346:     DMGetCoordinateDM(dac,&cdac);
347:     DMGetCoordinateDM(daf,&cdaf);

349:     DMGetCoordinates(dac,&coordsc);
350:     DMGetCoordinates(daf,&coordsf);

352:     DMCreateInterpolation(cdac,cdaf,&II,&scale);
353:     MatInterpolate(II,coordsc,coordsf);
354:     MatDestroy(&II);
355:     VecDestroy(&scale);
356:   }

358:   DMCreateInterpolation(dac,daf,&INTERP,NULL);

360:   DMCreateGlobalVector(dac,&ac);
361:   VecSet(ac,66.99);

363:   DMCreateGlobalVector(daf,&af);

365:   MatMult(INTERP,ac, af);

367:   {
368:     Vec       afexact;
369:     PetscReal nrm;
370:     PetscInt  N;

372:     DMCreateGlobalVector(daf,&afexact);
373:     VecSet(afexact,66.99);
374:     VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
375:     VecNorm(afexact,NORM_2,&nrm);
376:     VecGetSize(afexact,&N);
377:     PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)));
378:     VecDestroy(&afexact);
379:   }

381:   PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
382:   if (output) {
383:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv);
384:     VecView(ac, vv);
385:     PetscViewerDestroy(&vv);

387:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv);
388:     VecView(af, vv);
389:     PetscViewerDestroy(&vv);
390:   }

392:   MatDestroy(&INTERP);
393:   DMDestroy(&dac);
394:   DMDestroy(&daf);
395:   VecDestroy(&ac);
396:   VecDestroy(&af);
397:   return 0;
398: }

400: PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
401: {
402:   DM             dac,daf;
403:   PetscViewer    vv;
404:   Vec            ac,af;
405:   PetscInt       map_id,Mx,My;
406:   Mat            II,INTERP;
407:   Vec            scale;
408:   PetscBool      output = PETSC_FALSE;

411:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE, PETSC_DECIDE,1, /* 1 dof */1, /* stencil = 1 */NULL, NULL,&dac);
412:   DMSetFromOptions(dac);
413:   DMSetUp(dac);

415:   DMRefine(dac,MPI_COMM_NULL,&daf);
416:   DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);
417:   Mx--; My--;

419:   DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);
420:   DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);

422:   /* apply conformal mappings */
423:   map_id = 0;
424:   PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
425:   if (map_id >= 1) {
426:     DAApplyConformalMapping(dac,map_id);
427:   }

429:   {
430:     DM  cdaf,cdac;
431:     Vec coordsc,coordsf;

433:     DMGetCoordinateDM(dac,&cdac);
434:     DMGetCoordinateDM(daf,&cdaf);

436:     DMGetCoordinates(dac,&coordsc);
437:     DMGetCoordinates(daf,&coordsf);

439:     DMCreateInterpolation(cdac,cdaf,&II,&scale);
440:     MatInterpolate(II,coordsc,coordsf);
441:     MatDestroy(&II);
442:     VecDestroy(&scale);
443:   }

445:   DMCreateInterpolation(dac,daf,&INTERP,NULL);

447:   DMCreateGlobalVector(dac,&ac);
448:   DADefineXLinearField2D(dac,ac);

450:   DMCreateGlobalVector(daf,&af);
451:   MatMult(INTERP,ac, af);

453:   {
454:     Vec       afexact;
455:     PetscReal nrm;
456:     PetscInt  N;

458:     DMCreateGlobalVector(daf,&afexact);
459:     VecZeroEntries(afexact);
460:     DADefineXLinearField2D(daf,afexact);
461:     VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
462:     VecNorm(afexact,NORM_2,&nrm);
463:     VecGetSize(afexact,&N);
464:     PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)));
465:     VecDestroy(&afexact);
466:   }

468:   PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
469:   if (output) {
470:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv);
471:     VecView(ac, vv);
472:     PetscViewerDestroy(&vv);

474:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv);
475:     VecView(af, vv);
476:     PetscViewerDestroy(&vv);
477:   }

479:   MatDestroy(&INTERP);
480:   DMDestroy(&dac);
481:   DMDestroy(&daf);
482:   VecDestroy(&ac);
483:   VecDestroy(&af);
484:   return 0;
485: }

487: PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
488: {
490:   DM             dac,daf;
491:   PetscViewer    vv;
492:   Vec            ac,af;
493:   PetscInt       map_id,Mx,My,Mz;
494:   Mat            II,INTERP;
495:   Vec            scale;
496:   PetscBool      output = PETSC_FALSE;

499:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1, my+1,mz+1,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE,1, /* 1 dof */
500:                       1, /* stencil = 1 */NULL,NULL,NULL,&dac);
501:   DMSetFromOptions(dac);
502:   DMSetUp(dac);

504:   DMRefine(dac,MPI_COMM_NULL,&daf);
505:   DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);
506:   Mx--; My--; Mz--;

508:   DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);
509:   DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);

511:   /* apply trilinear mappings */
512:   /*DAApplyTrilinearMapping(dac);*/
513:   /* apply conformal mappings */
514:   map_id = 0;
515:   PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);
516:   if (map_id >= 1) {
517:     DAApplyConformalMapping(dac,map_id);
518:   }

520:   {
521:     DM  cdaf,cdac;
522:     Vec coordsc,coordsf;

524:     DMGetCoordinateDM(dac,&cdac);
525:     DMGetCoordinateDM(daf,&cdaf);

527:     DMGetCoordinates(dac,&coordsc);
528:     DMGetCoordinates(daf,&coordsf);

530:     DMCreateInterpolation(cdac,cdaf,&II,&scale);
531:     MatInterpolate(II,coordsc,coordsf);
532:     MatDestroy(&II);
533:     VecDestroy(&scale);
534:   }

536:   DMCreateInterpolation(dac,daf,&INTERP,NULL);

538:   DMCreateGlobalVector(dac,&ac);
539:   VecZeroEntries(ac);
540:   DADefineXLinearField3D(dac,ac);

542:   DMCreateGlobalVector(daf,&af);
543:   VecZeroEntries(af);

545:   MatMult(INTERP,ac, af);

547:   {
548:     Vec       afexact;
549:     PetscReal nrm;
550:     PetscInt  N;

552:     DMCreateGlobalVector(daf,&afexact);
553:     VecZeroEntries(afexact);
554:     DADefineXLinearField3D(daf,afexact);
555:     VecAXPY(afexact,-1.0,af); /* af <= af - afinterp */
556:     VecNorm(afexact,NORM_2,&nrm);
557:     VecGetSize(afexact,&N);
558:     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)));
559:     VecDestroy(&afexact);
560:   }

562:   PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);
563:   if (output) {
564:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv);
565:     VecView(ac, vv);
566:     PetscViewerDestroy(&vv);

568:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv);
569:     VecView(af, vv);
570:     PetscViewerDestroy(&vv);
571:   }

573:   MatDestroy(&INTERP);
574:   DMDestroy(&dac);
575:   DMDestroy(&daf);
576:   VecDestroy(&ac);
577:   VecDestroy(&af);
578:   return 0;
579: }

581: int main(int argc,char **argv)
582: {
583:   PetscInt       mx = 2,my = 2,mz = 2,l,nl,dim;

585:   PetscInitialize(&argc,&argv,0,help);
586:   PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);
587:   PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);
588:   PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);
589:   nl = 1;
590:   PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0);
591:   dim = 2;
592:   PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0);

594:   for (l=0; l<nl; l++) {
595:     if (dim==1) {
596:       da_test_RefineCoords1D(mx);
597:     } else if (dim==2) {
598:       da_test_RefineCoords2D(mx,my);
599:     } else if (dim==3) {
600:       da_test_RefineCoords3D(mx,my,mz);
601:     }
602:     mx = mx * 2;
603:     my = my * 2;
604:     mz = mz * 2;
605:   }
606:   PetscFinalize();
607:   return 0;
608: }

610: /*TEST

612:    test:
613:       suffix: 1d
614:       args: -mx 10 -nl 6 -dim 1

616:    test:
617:       suffix: 2d
618:       output_file: output/ex36_2d.out
619:       args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}

621:    test:
622:       suffix: 2dp1
623:       nsize: 8
624:       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
625:       timeoutfactor: 2

627:    test:
628:       suffix: 2dp2
629:       nsize: 8
630:       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
631:       timeoutfactor: 2

633:    test:
634:       suffix: 3d
635:       args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3

637:    test:
638:       suffix: 3dp1
639:       nsize: 32
640:       args: -mx 5 -my 5 -mz 5 -nl 3 -dim 3 -cmap 1 -da_refine_x 1 -da_refine_y 3 -da_refine_z 4

642: TEST*/