xref: /petsc/src/dm/tests/ex36.c (revision bef158480efac06de457f7a665168877ab3c2fd7) 
1 
2 static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";
3 
4 #include <petscdm.h>
5 #include <petscdmda.h>
6 
7 typedef struct _n_CCmplx CCmplx;
8 struct _n_CCmplx {
9   PetscReal real;
10   PetscReal imag;
11 };
12 
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 }
55 
56 PetscErrorCode DAApplyConformalMapping(DM da,PetscInt idx)
57 {
58   PetscErrorCode ierr;
59   PetscInt       i,n;
60   PetscInt       sx,nx,sy,ny,sz,nz,dim;
61   Vec            Gcoords;
62   PetscScalar    *XX;
63   PetscScalar    xx,yy,zz;
64   DM             cda;
65 
66   PetscFunctionBeginUser;
67   if (idx==0) {
68     PetscFunctionReturn(0);
69   } else if (idx==1) { /* dam break */
70     ierr = DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
71   } else if (idx==2) { /* stagnation in a corner */
72     ierr = DMDASetUniformCoordinates(da, -1.0,1.0, 0.0,1.0, -1.0,1.0);CHKERRQ(ierr);
73   } else if (idx==3) { /* nautilis */
74     ierr = DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
75   } else if (idx==4) {
76     ierr = DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
77   }
78 
79   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
80   ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);
81 
82   ierr = VecGetArray(Gcoords,&XX);CHKERRQ(ierr);
83   ierr = DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);CHKERRQ(ierr);
84   ierr = DMDAGetInfo(da, &dim, 0,0,0, 0,0,0, 0,0,0,0,0,0);CHKERRQ(ierr);
85   ierr = VecGetLocalSize(Gcoords,&n);CHKERRQ(ierr);
86   n    = n / dim;
87 
88   for (i=0; i<n; i++) {
89     if ((dim==3) && (idx!=2)) {
90       PetscScalar Ni[8];
91       PetscScalar xi   = XX[dim*i];
92       PetscScalar eta  = XX[dim*i+1];
93       PetscScalar zeta = XX[dim*i+2];
94       PetscScalar xn[] = {-1.0,1.0,-1.0,1.0,   -1.0,1.0,-1.0,1.0  };
95       PetscScalar yn[] = {-1.0,-1.0,1.0,1.0,   -1.0,-1.0,1.0,1.0  };
96       PetscScalar zn[] = {-0.1,-4.0,-0.2,-1.0,  0.1,4.0,0.2,1.0  };
97       PetscInt    p;
98 
99       Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
100       Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
101       Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
102       Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
103 
104       Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
105       Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
106       Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
107       Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
108 
109       xx = yy = zz = 0.0;
110       for (p=0; p<8; p++) {
111         xx += Ni[p]*xn[p];
112         yy += Ni[p]*yn[p];
113         zz += Ni[p]*zn[p];
114       }
115       XX[dim*i]   = xx;
116       XX[dim*i+1] = yy;
117       XX[dim*i+2] = zz;
118     }
119 
120     if (idx==1) {
121       CCmplx zeta,t1,t2;
122 
123       xx = XX[dim*i]   - 0.8;
124       yy = XX[dim*i+1] + 1.5;
125 
126       zeta.real = PetscRealPart(xx);
127       zeta.imag = PetscRealPart(yy);
128 
129       t1 = CCmplxPow(zeta,-1.0);
130       t2 = CCmplxAdd(zeta,t1);
131 
132       XX[dim*i]   = CCmplxRe(t2);
133       XX[dim*i+1] = CCmplxIm(t2);
134     } else if (idx==2) {
135       CCmplx zeta,t1;
136 
137       xx = XX[dim*i];
138       yy = XX[dim*i+1];
139       zeta.real = PetscRealPart(xx);
140       zeta.imag = PetscRealPart(yy);
141 
142       t1 = CCmplxSqrt(zeta);
143       XX[dim*i]   = CCmplxRe(t1);
144       XX[dim*i+1] = CCmplxIm(t1);
145     } else if (idx==3) {
146       CCmplx zeta,t1,t2;
147 
148       xx = XX[dim*i]   - 0.8;
149       yy = XX[dim*i+1] + 1.5;
150 
151       zeta.real   = PetscRealPart(xx);
152       zeta.imag   = PetscRealPart(yy);
153       t1          = CCmplxPow(zeta,-1.0);
154       t2          = CCmplxAdd(zeta,t1);
155       XX[dim*i]   = CCmplxRe(t2);
156       XX[dim*i+1] = CCmplxIm(t2);
157 
158       xx          = XX[dim*i];
159       yy          = XX[dim*i+1];
160       zeta.real   = PetscRealPart(xx);
161       zeta.imag   = PetscRealPart(yy);
162       t1          = CCmplxExp(zeta);
163       XX[dim*i]   = CCmplxRe(t1);
164       XX[dim*i+1] = CCmplxIm(t1);
165 
166       xx          = XX[dim*i] + 0.4;
167       yy          = XX[dim*i+1];
168       zeta.real   = PetscRealPart(xx);
169       zeta.imag   = PetscRealPart(yy);
170       t1          = CCmplxPow(zeta,2.0);
171       XX[dim*i]   = CCmplxRe(t1);
172       XX[dim*i+1] = CCmplxIm(t1);
173     } else if (idx==4) {
174       PetscScalar Ni[4];
175       PetscScalar xi   = XX[dim*i];
176       PetscScalar eta  = XX[dim*i+1];
177       PetscScalar xn[] = {0.0,2.0,0.2,3.5};
178       PetscScalar yn[] = {-1.3,0.0,2.0,4.0};
179       PetscInt    p;
180 
181       Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
182       Ni[1] = 0.25*(1.0+xi)*(1.0-eta);
183       Ni[2] = 0.25*(1.0-xi)*(1.0+eta);
184       Ni[3] = 0.25*(1.0+xi)*(1.0+eta);
185 
186       xx = yy = 0.0;
187       for (p=0; p<4; p++) {
188         xx += Ni[p]*xn[p];
189         yy += Ni[p]*yn[p];
190       }
191       XX[dim*i]   = xx;
192       XX[dim*i+1] = yy;
193     }
194   }
195   ierr = VecRestoreArray(Gcoords,&XX);CHKERRQ(ierr);
196   PetscFunctionReturn(0);
197 }
198 
199 PetscErrorCode DAApplyTrilinearMapping(DM da)
200 {
201   PetscErrorCode ierr;
202   PetscInt       i,j,k;
203   PetscInt       sx,nx,sy,ny,sz,nz;
204   Vec            Gcoords;
205   DMDACoor3d     ***XX;
206   PetscScalar    xx,yy,zz;
207   DM             cda;
208 
209   PetscFunctionBeginUser;
210   ierr = DMDASetUniformCoordinates(da, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
211   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
212   ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);
213 
214   ierr = DMDAVecGetArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
215   ierr = DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);CHKERRQ(ierr);
216 
217   for (i=sx; i<sx+nx; i++) {
218     for (j=sy; j<sy+ny; j++) {
219       for (k=sz; k<sz+nz; k++) {
220         PetscScalar Ni[8];
221         PetscScalar xi   = XX[k][j][i].x;
222         PetscScalar eta  = XX[k][j][i].y;
223         PetscScalar zeta = XX[k][j][i].z;
224         PetscScalar xn[] = {0.0,2.0,0.2,3.5,   0.0,2.1,0.23,3.125  };
225         PetscScalar yn[] = {-1.3,0.0,2.0,4.0,  -1.45,-0.1,2.24,3.79  };
226         PetscScalar zn[] = {0.0,0.3,-0.1,0.123,  0.956,1.32,1.12,0.798  };
227         PetscInt    p;
228 
229         Ni[0] = 0.125*(1.0-xi)*(1.0-eta)*(1.0-zeta);
230         Ni[1] = 0.125*(1.0+xi)*(1.0-eta)*(1.0-zeta);
231         Ni[2] = 0.125*(1.0-xi)*(1.0+eta)*(1.0-zeta);
232         Ni[3] = 0.125*(1.0+xi)*(1.0+eta)*(1.0-zeta);
233 
234         Ni[4] = 0.125*(1.0-xi)*(1.0-eta)*(1.0+zeta);
235         Ni[5] = 0.125*(1.0+xi)*(1.0-eta)*(1.0+zeta);
236         Ni[6] = 0.125*(1.0-xi)*(1.0+eta)*(1.0+zeta);
237         Ni[7] = 0.125*(1.0+xi)*(1.0+eta)*(1.0+zeta);
238 
239         xx = yy = zz = 0.0;
240         for (p=0; p<8; p++) {
241           xx += Ni[p]*xn[p];
242           yy += Ni[p]*yn[p];
243           zz += Ni[p]*zn[p];
244         }
245         XX[k][j][i].x = xx;
246         XX[k][j][i].y = yy;
247         XX[k][j][i].z = zz;
248       }
249     }
250   }
251   ierr = DMDAVecRestoreArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
252   PetscFunctionReturn(0);
253 }
254 
255 PetscErrorCode DADefineXLinearField2D(DM da,Vec field)
256 {
257   PetscErrorCode ierr;
258   PetscInt       i,j;
259   PetscInt       sx,nx,sy,ny;
260   Vec            Gcoords;
261   DMDACoor2d     **XX;
262   PetscScalar    **FF;
263   DM             cda;
264 
265   PetscFunctionBeginUser;
266   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
267   ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);
268 
269   ierr = DMDAVecGetArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
270   ierr = DMDAVecGetArray(da,field,&FF);CHKERRQ(ierr);
271 
272   ierr = DMDAGetCorners(da,&sx,&sy,0,&nx,&ny,0);CHKERRQ(ierr);
273 
274   for (i=sx; i<sx+nx; i++) {
275     for (j=sy; j<sy+ny; j++) {
276       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;
277     }
278   }
279 
280   ierr = DMDAVecRestoreArray(da,field,&FF);CHKERRQ(ierr);
281   ierr = DMDAVecRestoreArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
282   PetscFunctionReturn(0);
283 }
284 
285 PetscErrorCode DADefineXLinearField3D(DM da,Vec field)
286 {
287   PetscErrorCode ierr;
288   PetscInt       i,j,k;
289   PetscInt       sx,nx,sy,ny,sz,nz;
290   Vec            Gcoords;
291   DMDACoor3d     ***XX;
292   PetscScalar    ***FF;
293   DM             cda;
294 
295   PetscFunctionBeginUser;
296   ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
297   ierr = DMGetCoordinates(da,&Gcoords);CHKERRQ(ierr);
298 
299   ierr = DMDAVecGetArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
300   ierr = DMDAVecGetArray(da,field,&FF);CHKERRQ(ierr);
301 
302   ierr = DMDAGetCorners(da,&sx,&sy,&sz,&nx,&ny,&nz);CHKERRQ(ierr);
303 
304   for (k=sz; k<sz+nz; k++) {
305     for (j=sy; j<sy+ny; j++) {
306       for (i=sx; i<sx+nx; i++) {
307         FF[k][j][i] = 10.0
308                 + 4.05 * XX[k][j][i].x
309                 + 5.50 * XX[k][j][i].y
310                 + 1.33 * XX[k][j][i].z
311                 + 2.03 * XX[k][j][i].x * XX[k][j][i].y
312                 + 0.03 * XX[k][j][i].x * XX[k][j][i].z
313                 + 0.83 * XX[k][j][i].y * XX[k][j][i].z
314                 + 3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
315       }
316     }
317   }
318 
319   ierr = DMDAVecRestoreArray(da,field,&FF);CHKERRQ(ierr);
320   ierr = DMDAVecRestoreArrayRead(cda,Gcoords,&XX);CHKERRQ(ierr);
321   PetscFunctionReturn(0);
322 }
323 
324 PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
325 {
326   PetscErrorCode ierr;
327   DM             dac,daf;
328   PetscViewer    vv;
329   Vec            ac,af;
330   PetscInt       Mx;
331   Mat            II,INTERP;
332   Vec            scale;
333   PetscBool      output = PETSC_FALSE;
334 
335   PetscFunctionBeginUser;
336   ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,mx+1,1, /* 1 dof */1, /* stencil = 1 */NULL,&dac);CHKERRQ(ierr);
337   ierr = DMSetFromOptions(dac);CHKERRQ(ierr);
338   ierr = DMSetUp(dac);CHKERRQ(ierr);
339 
340   ierr = DMRefine(dac,MPI_COMM_NULL,&daf);CHKERRQ(ierr);
341   ierr = DMDAGetInfo(daf,0,&Mx,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
342   Mx--;
343 
344   ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
345   ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
346 
347   {
348     DM  cdaf,cdac;
349     Vec coordsc,coordsf;
350 
351     ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
352     ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
353 
354     ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
355     ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
356 
357     ierr = DMCreateInterpolation(cdac,cdaf,&II,&scale);CHKERRQ(ierr);
358     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
359     ierr = MatDestroy(&II);CHKERRQ(ierr);
360     ierr = VecDestroy(&scale);CHKERRQ(ierr);
361   }
362 
363   ierr = DMCreateInterpolation(dac,daf,&INTERP,NULL);CHKERRQ(ierr);
364 
365   ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
366   ierr = VecSet(ac,66.99);CHKERRQ(ierr);
367 
368   ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
369 
370   ierr = MatMult(INTERP,ac, af);CHKERRQ(ierr);
371 
372   {
373     Vec       afexact;
374     PetscReal nrm;
375     PetscInt  N;
376 
377     ierr = DMCreateGlobalVector(daf,&afexact);CHKERRQ(ierr);
378     ierr = VecSet(afexact,66.99);CHKERRQ(ierr);
379     ierr = VecAXPY(afexact,-1.0,af);CHKERRQ(ierr); /* af <= af - afinterp */
380     ierr = VecNorm(afexact,NORM_2,&nrm);CHKERRQ(ierr);
381     ierr = VecGetSize(afexact,&N);CHKERRQ(ierr);
382     ierr = PetscPrintf(PETSC_COMM_WORLD,"%D=>%D, interp err = %1.4e\n",mx,Mx,(double)(nrm/PetscSqrtReal((PetscReal)N)));CHKERRQ(ierr);
383     ierr = VecDestroy(&afexact);CHKERRQ(ierr);
384   }
385 
386   ierr = PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);CHKERRQ(ierr);
387   if (output) {
388     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtk", &vv);CHKERRQ(ierr);
389     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
390     ierr = DMView(dac, vv);CHKERRQ(ierr);
391     ierr = VecView(ac, vv);CHKERRQ(ierr);
392     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
393     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
394 
395     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtk", &vv);CHKERRQ(ierr);
396     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
397     ierr = DMView(daf, vv);CHKERRQ(ierr);
398     ierr = VecView(af, vv);CHKERRQ(ierr);
399     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
400     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
401   }
402 
403   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
404   ierr = DMDestroy(&dac);CHKERRQ(ierr);
405   ierr = DMDestroy(&daf);CHKERRQ(ierr);
406   ierr = VecDestroy(&ac);CHKERRQ(ierr);
407   ierr = VecDestroy(&af);CHKERRQ(ierr);
408   PetscFunctionReturn(0);
409 }
410 
411 PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
412 {
413   PetscErrorCode ierr;
414   DM             dac,daf;
415   PetscViewer    vv;
416   Vec            ac,af;
417   PetscInt       map_id,Mx,My;
418   Mat            II,INTERP;
419   Vec            scale;
420   PetscBool      output = PETSC_FALSE;
421 
422   PetscFunctionBeginUser;
423   ierr = 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);CHKERRQ(ierr);
424   ierr = DMSetFromOptions(dac);CHKERRQ(ierr);
425   ierr = DMSetUp(dac);CHKERRQ(ierr);
426 
427   ierr = DMRefine(dac,MPI_COMM_NULL,&daf);CHKERRQ(ierr);
428   ierr = DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
429   Mx--; My--;
430 
431   ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
432   ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
433 
434   /* apply conformal mappings */
435   map_id = 0;
436   ierr   = PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);CHKERRQ(ierr);
437   if (map_id >= 1) {
438     ierr = DAApplyConformalMapping(dac,map_id);CHKERRQ(ierr);
439   }
440 
441   {
442     DM  cdaf,cdac;
443     Vec coordsc,coordsf;
444 
445     ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
446     ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
447 
448     ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
449     ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
450 
451     ierr = DMCreateInterpolation(cdac,cdaf,&II,&scale);CHKERRQ(ierr);
452     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
453     ierr = MatDestroy(&II);CHKERRQ(ierr);
454     ierr = VecDestroy(&scale);CHKERRQ(ierr);
455   }
456 
457 
458   ierr = DMCreateInterpolation(dac,daf,&INTERP,NULL);CHKERRQ(ierr);
459 
460   ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
461   ierr = DADefineXLinearField2D(dac,ac);CHKERRQ(ierr);
462 
463   ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
464   ierr = MatMult(INTERP,ac, af);CHKERRQ(ierr);
465 
466   {
467     Vec       afexact;
468     PetscReal nrm;
469     PetscInt  N;
470 
471     ierr = DMCreateGlobalVector(daf,&afexact);CHKERRQ(ierr);
472     ierr = VecZeroEntries(afexact);CHKERRQ(ierr);
473     ierr = DADefineXLinearField2D(daf,afexact);CHKERRQ(ierr);
474     ierr = VecAXPY(afexact,-1.0,af);CHKERRQ(ierr); /* af <= af - afinterp */
475     ierr = VecNorm(afexact,NORM_2,&nrm);CHKERRQ(ierr);
476     ierr = VecGetSize(afexact,&N);CHKERRQ(ierr);
477     ierr = PetscPrintf(PETSC_COMM_WORLD,"[%D x %D]=>[%D x %D], interp err = %1.4e\n",mx,my,Mx,My,(double)(nrm/PetscSqrtReal((PetscReal)N)));CHKERRQ(ierr);
478     ierr = VecDestroy(&afexact);CHKERRQ(ierr);
479   }
480 
481   ierr = PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);CHKERRQ(ierr);
482   if (output) {
483     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtk", &vv);CHKERRQ(ierr);
484     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
485     ierr = DMView(dac, vv);CHKERRQ(ierr);
486     ierr = VecView(ac, vv);CHKERRQ(ierr);
487     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
488     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
489 
490     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtk", &vv);CHKERRQ(ierr);
491     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
492     ierr = DMView(daf, vv);CHKERRQ(ierr);
493     ierr = VecView(af, vv);CHKERRQ(ierr);
494     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
495     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
496   }
497 
498   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
499   ierr = DMDestroy(&dac);CHKERRQ(ierr);
500   ierr = DMDestroy(&daf);CHKERRQ(ierr);
501   ierr = VecDestroy(&ac);CHKERRQ(ierr);
502   ierr = VecDestroy(&af);CHKERRQ(ierr);
503   PetscFunctionReturn(0);
504 }
505 
506 PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
507 {
508   PetscErrorCode ierr;
509   DM             dac,daf;
510   PetscViewer    vv;
511   Vec            ac,af;
512   PetscInt       map_id,Mx,My,Mz;
513   Mat            II,INTERP;
514   Vec            scale;
515   PetscBool      output = PETSC_FALSE;
516 
517   PetscFunctionBeginUser;
518   ierr = 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 */
519                       1, /* stencil = 1 */NULL,NULL,NULL,&dac);CHKERRQ(ierr);
520   ierr = DMSetFromOptions(dac);CHKERRQ(ierr);
521   ierr = DMSetUp(dac);CHKERRQ(ierr);
522 
523   ierr = DMRefine(dac,MPI_COMM_NULL,&daf);CHKERRQ(ierr);
524   ierr = DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
525   Mx--; My--; Mz--;
526 
527   ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
528   ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
529 
530   /* apply trilinear mappings */
531   /*ierr = DAApplyTrilinearMapping(dac);CHKERRQ(ierr);*/
532   /* apply conformal mappings */
533   map_id = 0;
534   ierr   = PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);CHKERRQ(ierr);
535   if (map_id >= 1) {
536     ierr = DAApplyConformalMapping(dac,map_id);CHKERRQ(ierr);
537   }
538 
539   {
540     DM  cdaf,cdac;
541     Vec coordsc,coordsf;
542 
543     ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
544     ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
545 
546     ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
547     ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
548 
549     ierr = DMCreateInterpolation(cdac,cdaf,&II,&scale);CHKERRQ(ierr);
550     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
551     ierr = MatDestroy(&II);CHKERRQ(ierr);
552     ierr = VecDestroy(&scale);CHKERRQ(ierr);
553   }
554 
555   ierr = DMCreateInterpolation(dac,daf,&INTERP,NULL);CHKERRQ(ierr);
556 
557   ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
558   ierr = VecZeroEntries(ac);CHKERRQ(ierr);
559   ierr = DADefineXLinearField3D(dac,ac);CHKERRQ(ierr);
560 
561   ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
562   ierr = VecZeroEntries(af);CHKERRQ(ierr);
563 
564   ierr = MatMult(INTERP,ac, af);CHKERRQ(ierr);
565 
566   {
567     Vec       afexact;
568     PetscReal nrm;
569     PetscInt  N;
570 
571     ierr = DMCreateGlobalVector(daf,&afexact);CHKERRQ(ierr);
572     ierr = VecZeroEntries(afexact);CHKERRQ(ierr);
573     ierr = DADefineXLinearField3D(daf,afexact);CHKERRQ(ierr);
574     ierr = VecAXPY(afexact,-1.0,af);CHKERRQ(ierr); /* af <= af - afinterp */
575     ierr = VecNorm(afexact,NORM_2,&nrm);CHKERRQ(ierr);
576     ierr = VecGetSize(afexact,&N);CHKERRQ(ierr);
577     ierr = 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)));CHKERRQ(ierr);
578     ierr = VecDestroy(&afexact);CHKERRQ(ierr);
579   }
580 
581   PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);CHKERRQ(ierr);
582   if (output) {
583     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtk", &vv);CHKERRQ(ierr);
584     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
585     ierr = DMView(dac, vv);CHKERRQ(ierr);
586     ierr = VecView(ac, vv);CHKERRQ(ierr);
587     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
588     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
589 
590     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtk", &vv);CHKERRQ(ierr);
591     ierr = PetscViewerPushFormat(vv, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
592     ierr = DMView(daf, vv);CHKERRQ(ierr);
593     ierr = VecView(af, vv);CHKERRQ(ierr);
594     ierr = PetscViewerPopFormat(vv);CHKERRQ(ierr);
595     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
596   }
597 
598   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
599   ierr = DMDestroy(&dac);CHKERRQ(ierr);
600   ierr = DMDestroy(&daf);CHKERRQ(ierr);
601   ierr = VecDestroy(&ac);CHKERRQ(ierr);
602   ierr = VecDestroy(&af);CHKERRQ(ierr);
603   PetscFunctionReturn(0);
604 }
605 
606 int main(int argc,char **argv)
607 {
608   PetscErrorCode ierr;
609   PetscInt       mx = 2,my = 2,mz = 2,l,nl,dim;
610 
611   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
612   ierr = PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);CHKERRQ(ierr);
613   ierr = PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);CHKERRQ(ierr);
614   ierr = PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);CHKERRQ(ierr);
615   nl = 1;
616   ierr = PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0);CHKERRQ(ierr);
617   dim = 2;
618   ierr = PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0);CHKERRQ(ierr);
619 
620   for (l=0; l<nl; l++) {
621     if (dim==1) {
622       ierr = da_test_RefineCoords1D(mx);CHKERRQ(ierr);
623     } else if (dim==2) {
624       ierr = da_test_RefineCoords2D(mx,my);CHKERRQ(ierr);
625     } else if (dim==3) {
626       ierr = da_test_RefineCoords3D(mx,my,mz);CHKERRQ(ierr);
627     }
628     mx = mx * 2;
629     my = my * 2;
630     mz = mz * 2;
631   }
632   ierr = PetscFinalize();
633   return ierr;
634 }
635 
636 
637 /*TEST
638 
639    test:
640       suffix: 1d
641       args: -mx 10 -nl 6 -dim 1
642 
643    test:
644       suffix: 2d
645 
646       test:
647          args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 0
648 
649       test:
650          args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 1
651 
652       test:
653          args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 2
654 
655       test:
656          args: -mx 10 -my 10 -nl 6 -dim 2 -cmap 3
657 
658    test:
659       suffix: 2dp1
660       nsize: 8
661       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
662       timeoutfactor: 2
663 
664    test:
665       suffix: 2dp2
666       nsize: 8
667       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
668       timeoutfactor: 2
669 
670    test:
671       suffix: 3d
672       args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
673 
674    test:
675       suffix: 3dp1
676       nsize: 32
677       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
678 
679 TEST*/
680