Actual source code: ex36.c

  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,PETSC_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:     DMGetInterpolation(cdac,cdaf,&II,&scale);
374:     MatInterpolate(II,coordsc,coordsf);
375:     MatDestroy(&II);
376:     VecDestroy(&scale);
377:   }
378: 
379:   DMGetInterpolation(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,PETSC_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:     DMGetInterpolation(cdac,cdaf,&II,&scale);
474:     MatInterpolate(II,coordsc,coordsf);
475:     MatDestroy(&II);
476:     VecDestroy(&scale);
477:   }
478: 
479: 
480:   DMGetInterpolation(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,PETSC_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:     DMGetInterpolation(cdac,cdaf,&II,&scale);
577:     MatInterpolate(II,coordsc,coordsf);
578:     MatDestroy(&II);
579:     VecDestroy(&scale);
580:   }
581: 
582:   DMGetInterpolation(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: }