xref: /petsc/src/dm/tests/ex36.c (revision 030f984af8d8bb4c203755d35bded3c05b3d83ce) 
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.vtr", &vv);CHKERRQ(ierr);
389     ierr = VecView(ac, vv);CHKERRQ(ierr);
390     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
391 
392     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv);CHKERRQ(ierr);
393     ierr = VecView(af, vv);CHKERRQ(ierr);
394     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
395   }
396 
397   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
398   ierr = DMDestroy(&dac);CHKERRQ(ierr);
399   ierr = DMDestroy(&daf);CHKERRQ(ierr);
400   ierr = VecDestroy(&ac);CHKERRQ(ierr);
401   ierr = VecDestroy(&af);CHKERRQ(ierr);
402   PetscFunctionReturn(0);
403 }
404 
405 PetscErrorCode da_test_RefineCoords2D(PetscInt mx,PetscInt my)
406 {
407   PetscErrorCode ierr;
408   DM             dac,daf;
409   PetscViewer    vv;
410   Vec            ac,af;
411   PetscInt       map_id,Mx,My;
412   Mat            II,INTERP;
413   Vec            scale;
414   PetscBool      output = PETSC_FALSE;
415 
416   PetscFunctionBeginUser;
417   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);
418   ierr = DMSetFromOptions(dac);CHKERRQ(ierr);
419   ierr = DMSetUp(dac);CHKERRQ(ierr);
420 
421   ierr = DMRefine(dac,MPI_COMM_NULL,&daf);CHKERRQ(ierr);
422   ierr = DMDAGetInfo(daf,0,&Mx,&My,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
423   Mx--; My--;
424 
425   ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
426   ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
427 
428   /* apply conformal mappings */
429   map_id = 0;
430   ierr   = PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);CHKERRQ(ierr);
431   if (map_id >= 1) {
432     ierr = DAApplyConformalMapping(dac,map_id);CHKERRQ(ierr);
433   }
434 
435   {
436     DM  cdaf,cdac;
437     Vec coordsc,coordsf;
438 
439     ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
440     ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
441 
442     ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
443     ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
444 
445     ierr = DMCreateInterpolation(cdac,cdaf,&II,&scale);CHKERRQ(ierr);
446     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
447     ierr = MatDestroy(&II);CHKERRQ(ierr);
448     ierr = VecDestroy(&scale);CHKERRQ(ierr);
449   }
450 
451   ierr = DMCreateInterpolation(dac,daf,&INTERP,NULL);CHKERRQ(ierr);
452 
453   ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
454   ierr = DADefineXLinearField2D(dac,ac);CHKERRQ(ierr);
455 
456   ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
457   ierr = MatMult(INTERP,ac, af);CHKERRQ(ierr);
458 
459   {
460     Vec       afexact;
461     PetscReal nrm;
462     PetscInt  N;
463 
464     ierr = DMCreateGlobalVector(daf,&afexact);CHKERRQ(ierr);
465     ierr = VecZeroEntries(afexact);CHKERRQ(ierr);
466     ierr = DADefineXLinearField2D(daf,afexact);CHKERRQ(ierr);
467     ierr = VecAXPY(afexact,-1.0,af);CHKERRQ(ierr); /* af <= af - afinterp */
468     ierr = VecNorm(afexact,NORM_2,&nrm);CHKERRQ(ierr);
469     ierr = VecGetSize(afexact,&N);CHKERRQ(ierr);
470     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);
471     ierr = VecDestroy(&afexact);CHKERRQ(ierr);
472   }
473 
474   ierr = PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);CHKERRQ(ierr);
475   if (output) {
476     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv);CHKERRQ(ierr);
477     ierr = VecView(ac, vv);CHKERRQ(ierr);
478     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
479 
480     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv);CHKERRQ(ierr);
481     ierr = VecView(af, vv);CHKERRQ(ierr);
482     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
483   }
484 
485   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
486   ierr = DMDestroy(&dac);CHKERRQ(ierr);
487   ierr = DMDestroy(&daf);CHKERRQ(ierr);
488   ierr = VecDestroy(&ac);CHKERRQ(ierr);
489   ierr = VecDestroy(&af);CHKERRQ(ierr);
490   PetscFunctionReturn(0);
491 }
492 
493 PetscErrorCode da_test_RefineCoords3D(PetscInt mx,PetscInt my,PetscInt mz)
494 {
495   PetscErrorCode ierr;
496   DM             dac,daf;
497   PetscViewer    vv;
498   Vec            ac,af;
499   PetscInt       map_id,Mx,My,Mz;
500   Mat            II,INTERP;
501   Vec            scale;
502   PetscBool      output = PETSC_FALSE;
503 
504   PetscFunctionBeginUser;
505   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 */
506                       1, /* stencil = 1 */NULL,NULL,NULL,&dac);CHKERRQ(ierr);
507   ierr = DMSetFromOptions(dac);CHKERRQ(ierr);
508   ierr = DMSetUp(dac);CHKERRQ(ierr);
509 
510   ierr = DMRefine(dac,MPI_COMM_NULL,&daf);CHKERRQ(ierr);
511   ierr = DMDAGetInfo(daf,0,&Mx,&My,&Mz,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
512   Mx--; My--; Mz--;
513 
514   ierr = DMDASetUniformCoordinates(dac, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
515   ierr = DMDASetUniformCoordinates(daf, -1.0,1.0, -1.0,1.0, -1.0,1.0);CHKERRQ(ierr);
516 
517   /* apply trilinear mappings */
518   /*ierr = DAApplyTrilinearMapping(dac);CHKERRQ(ierr);*/
519   /* apply conformal mappings */
520   map_id = 0;
521   ierr   = PetscOptionsGetInt(NULL,NULL,"-cmap", &map_id,NULL);CHKERRQ(ierr);
522   if (map_id >= 1) {
523     ierr = DAApplyConformalMapping(dac,map_id);CHKERRQ(ierr);
524   }
525 
526   {
527     DM  cdaf,cdac;
528     Vec coordsc,coordsf;
529 
530     ierr = DMGetCoordinateDM(dac,&cdac);CHKERRQ(ierr);
531     ierr = DMGetCoordinateDM(daf,&cdaf);CHKERRQ(ierr);
532 
533     ierr = DMGetCoordinates(dac,&coordsc);CHKERRQ(ierr);
534     ierr = DMGetCoordinates(daf,&coordsf);CHKERRQ(ierr);
535 
536     ierr = DMCreateInterpolation(cdac,cdaf,&II,&scale);CHKERRQ(ierr);
537     ierr = MatInterpolate(II,coordsc,coordsf);CHKERRQ(ierr);
538     ierr = MatDestroy(&II);CHKERRQ(ierr);
539     ierr = VecDestroy(&scale);CHKERRQ(ierr);
540   }
541 
542   ierr = DMCreateInterpolation(dac,daf,&INTERP,NULL);CHKERRQ(ierr);
543 
544   ierr = DMCreateGlobalVector(dac,&ac);CHKERRQ(ierr);
545   ierr = VecZeroEntries(ac);CHKERRQ(ierr);
546   ierr = DADefineXLinearField3D(dac,ac);CHKERRQ(ierr);
547 
548   ierr = DMCreateGlobalVector(daf,&af);CHKERRQ(ierr);
549   ierr = VecZeroEntries(af);CHKERRQ(ierr);
550 
551   ierr = MatMult(INTERP,ac, af);CHKERRQ(ierr);
552 
553   {
554     Vec       afexact;
555     PetscReal nrm;
556     PetscInt  N;
557 
558     ierr = DMCreateGlobalVector(daf,&afexact);CHKERRQ(ierr);
559     ierr = VecZeroEntries(afexact);CHKERRQ(ierr);
560     ierr = DADefineXLinearField3D(daf,afexact);CHKERRQ(ierr);
561     ierr = VecAXPY(afexact,-1.0,af);CHKERRQ(ierr); /* af <= af - afinterp */
562     ierr = VecNorm(afexact,NORM_2,&nrm);CHKERRQ(ierr);
563     ierr = VecGetSize(afexact,&N);CHKERRQ(ierr);
564     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);
565     ierr = VecDestroy(&afexact);CHKERRQ(ierr);
566   }
567 
568   ierr = PetscOptionsGetBool(NULL,NULL,"-output",&output,NULL);CHKERRQ(ierr);
569   if (output) {
570     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv);CHKERRQ(ierr);
571     ierr = VecView(ac, vv);CHKERRQ(ierr);
572     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
573 
574     ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv);CHKERRQ(ierr);
575     ierr = VecView(af, vv);CHKERRQ(ierr);
576     ierr = PetscViewerDestroy(&vv);CHKERRQ(ierr);
577   }
578 
579   ierr = MatDestroy(&INTERP);CHKERRQ(ierr);
580   ierr = DMDestroy(&dac);CHKERRQ(ierr);
581   ierr = DMDestroy(&daf);CHKERRQ(ierr);
582   ierr = VecDestroy(&ac);CHKERRQ(ierr);
583   ierr = VecDestroy(&af);CHKERRQ(ierr);
584   PetscFunctionReturn(0);
585 }
586 
587 int main(int argc,char **argv)
588 {
589   PetscErrorCode ierr;
590   PetscInt       mx = 2,my = 2,mz = 2,l,nl,dim;
591 
592   ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
593   ierr = PetscOptionsGetInt(NULL,NULL,"-mx", &mx, 0);CHKERRQ(ierr);
594   ierr = PetscOptionsGetInt(NULL,NULL,"-my", &my, 0);CHKERRQ(ierr);
595   ierr = PetscOptionsGetInt(NULL,NULL,"-mz", &mz, 0);CHKERRQ(ierr);
596   nl = 1;
597   ierr = PetscOptionsGetInt(NULL,NULL,"-nl", &nl, 0);CHKERRQ(ierr);
598   dim = 2;
599   ierr = PetscOptionsGetInt(NULL,NULL,"-dim", &dim, 0);CHKERRQ(ierr);
600 
601   for (l=0; l<nl; l++) {
602     if (dim==1) {
603       ierr = da_test_RefineCoords1D(mx);CHKERRQ(ierr);
604     } else if (dim==2) {
605       ierr = da_test_RefineCoords2D(mx,my);CHKERRQ(ierr);
606     } else if (dim==3) {
607       ierr = da_test_RefineCoords3D(mx,my,mz);CHKERRQ(ierr);
608     }
609     mx = mx * 2;
610     my = my * 2;
611     mz = mz * 2;
612   }
613   ierr = PetscFinalize();
614   return ierr;
615 }
616 
617 /*TEST
618 
619    test:
620       suffix: 1d
621       args: -mx 10 -nl 6 -dim 1
622 
623    test:
624       suffix: 2d
625       output_file: output/ex36_2d.out
626       args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}
627 
628    test:
629       suffix: 2dp1
630       nsize: 8
631       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
632       timeoutfactor: 2
633 
634    test:
635       suffix: 2dp2
636       nsize: 8
637       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
638       timeoutfactor: 2
639 
640    test:
641       suffix: 3d
642       args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3
643 
644    test:
645       suffix: 3dp1
646       nsize: 32
647       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
648 
649 TEST*/
650