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;

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

 76:   PetscCall(DMGetCoordinateDM(da, &cda));
 77:   PetscCall(DMGetCoordinates(da, &Gcoords));

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

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

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

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

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

117:     if (idx == 1) {
118:       CCmplx zeta, t1, t2;

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

123:       zeta.real = PetscRealPart(xx);
124:       zeta.imag = PetscRealPart(yy);

126:       t1 = CCmplxPow(zeta, -1.0);
127:       t2 = CCmplxAdd(zeta, t1);

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

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

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

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

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

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

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

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

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

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

205:   PetscFunctionBeginUser;
206:   PetscCall(DMDASetUniformCoordinates(da, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
207:   PetscCall(DMGetCoordinateDM(da, &cda));
208:   PetscCall(DMGetCoordinates(da, &Gcoords));

210:   PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
211:   PetscCall(DMDAGetCorners(da, &sx, &sy, &sz, &nx, &ny, &nz));

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

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

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

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

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

260:   PetscFunctionBeginUser;
261:   PetscCall(DMGetCoordinateDM(da, &cda));
262:   PetscCall(DMGetCoordinates(da, &Gcoords));

264:   PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
265:   PetscCall(DMDAVecGetArray(da, field, &FF));

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

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

273:   PetscCall(DMDAVecRestoreArray(da, field, &FF));
274:   PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

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

287:   PetscFunctionBeginUser;
288:   PetscCall(DMGetCoordinateDM(da, &cda));
289:   PetscCall(DMGetCoordinates(da, &Gcoords));

291:   PetscCall(DMDAVecGetArrayRead(cda, Gcoords, &XX));
292:   PetscCall(DMDAVecGetArray(da, field, &FF));

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

296:   for (k = sz; k < sz + nz; k++) {
297:     for (j = sy; j < sy + ny; j++) {
298:       for (i = sx; i < sx + nx; i++) {
299:         FF[k][j][i] = 10.0 + 4.05 * XX[k][j][i].x + 5.50 * XX[k][j][i].y + 1.33 * XX[k][j][i].z + 2.03 * XX[k][j][i].x * XX[k][j][i].y + 0.03 * XX[k][j][i].x * XX[k][j][i].z + 0.83 * XX[k][j][i].y * XX[k][j][i].z +
300:                       3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
301:       }
302:     }
303:   }

305:   PetscCall(DMDAVecRestoreArray(da, field, &FF));
306:   PetscCall(DMDAVecRestoreArrayRead(cda, Gcoords, &XX));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: PetscErrorCode da_test_RefineCoords1D(PetscInt mx)
311: {
312:   DM          dac, daf;
313:   PetscViewer vv;
314:   Vec         ac, af;
315:   PetscInt    Mx;
316:   Mat         II, INTERP;
317:   Vec         scale;
318:   PetscBool   output = PETSC_FALSE;

320:   PetscFunctionBeginUser;
321:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, mx + 1, 1, /* 1 dof */ 1, /* stencil = 1 */ NULL, &dac));
322:   PetscCall(DMSetFromOptions(dac));
323:   PetscCall(DMSetUp(dac));

325:   PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
326:   PetscCall(DMDAGetInfo(daf, 0, &Mx, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
327:   Mx--;

329:   PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
330:   PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));

332:   {
333:     DM  cdaf, cdac;
334:     Vec coordsc, coordsf;

336:     PetscCall(DMGetCoordinateDM(dac, &cdac));
337:     PetscCall(DMGetCoordinateDM(daf, &cdaf));

339:     PetscCall(DMGetCoordinates(dac, &coordsc));
340:     PetscCall(DMGetCoordinates(daf, &coordsf));

342:     PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
343:     PetscCall(MatInterpolate(II, coordsc, coordsf));
344:     PetscCall(MatDestroy(&II));
345:     PetscCall(VecDestroy(&scale));
346:   }

348:   PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

350:   PetscCall(DMCreateGlobalVector(dac, &ac));
351:   PetscCall(VecSet(ac, 66.99));

353:   PetscCall(DMCreateGlobalVector(daf, &af));

355:   PetscCall(MatMult(INTERP, ac, af));

357:   {
358:     Vec       afexact;
359:     PetscReal nrm;
360:     PetscInt  N;

362:     PetscCall(DMCreateGlobalVector(daf, &afexact));
363:     PetscCall(VecSet(afexact, 66.99));
364:     PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
365:     PetscCall(VecNorm(afexact, NORM_2, &nrm));
366:     PetscCall(VecGetSize(afexact, &N));
367:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%" PetscInt_FMT "=>%" PetscInt_FMT ", interp err = %1.4e\n", mx, Mx, (double)(nrm / PetscSqrtReal((PetscReal)N))));
368:     PetscCall(VecDestroy(&afexact));
369:   }

371:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
372:   if (output) {
373:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_1D.vtr", &vv));
374:     PetscCall(VecView(ac, vv));
375:     PetscCall(PetscViewerDestroy(&vv));

377:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_1D.vtr", &vv));
378:     PetscCall(VecView(af, vv));
379:     PetscCall(PetscViewerDestroy(&vv));
380:   }

382:   PetscCall(MatDestroy(&INTERP));
383:   PetscCall(DMDestroy(&dac));
384:   PetscCall(DMDestroy(&daf));
385:   PetscCall(VecDestroy(&ac));
386:   PetscCall(VecDestroy(&af));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: PetscErrorCode da_test_RefineCoords2D(PetscInt mx, PetscInt my)
391: {
392:   DM          dac, daf;
393:   PetscViewer vv;
394:   Vec         ac, af;
395:   PetscInt    map_id, Mx, My;
396:   Mat         II, INTERP;
397:   Vec         scale;
398:   PetscBool   output = PETSC_FALSE;

400:   PetscFunctionBeginUser;
401:   PetscCall(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));
402:   PetscCall(DMSetFromOptions(dac));
403:   PetscCall(DMSetUp(dac));

405:   PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
406:   PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
407:   Mx--;
408:   My--;

410:   PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));
411:   PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, PETSC_DECIDE, PETSC_DECIDE));

413:   /* apply conformal mappings */
414:   map_id = 0;
415:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
416:   if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));

418:   {
419:     DM  cdaf, cdac;
420:     Vec coordsc, coordsf;

422:     PetscCall(DMGetCoordinateDM(dac, &cdac));
423:     PetscCall(DMGetCoordinateDM(daf, &cdaf));

425:     PetscCall(DMGetCoordinates(dac, &coordsc));
426:     PetscCall(DMGetCoordinates(daf, &coordsf));

428:     PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
429:     PetscCall(MatInterpolate(II, coordsc, coordsf));
430:     PetscCall(MatDestroy(&II));
431:     PetscCall(VecDestroy(&scale));
432:   }

434:   PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

436:   PetscCall(DMCreateGlobalVector(dac, &ac));
437:   PetscCall(DADefineXLinearField2D(dac, ac));

439:   PetscCall(DMCreateGlobalVector(daf, &af));
440:   PetscCall(MatMult(INTERP, ac, af));

442:   {
443:     Vec       afexact;
444:     PetscReal nrm;
445:     PetscInt  N;

447:     PetscCall(DMCreateGlobalVector(daf, &afexact));
448:     PetscCall(VecZeroEntries(afexact));
449:     PetscCall(DADefineXLinearField2D(daf, afexact));
450:     PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
451:     PetscCall(VecNorm(afexact, NORM_2, &nrm));
452:     PetscCall(VecGetSize(afexact, &N));
453:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, Mx, My, (double)(nrm / PetscSqrtReal((PetscReal)N))));
454:     PetscCall(VecDestroy(&afexact));
455:   }

457:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
458:   if (output) {
459:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_2D.vtr", &vv));
460:     PetscCall(VecView(ac, vv));
461:     PetscCall(PetscViewerDestroy(&vv));

463:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_2D.vtr", &vv));
464:     PetscCall(VecView(af, vv));
465:     PetscCall(PetscViewerDestroy(&vv));
466:   }

468:   PetscCall(MatDestroy(&INTERP));
469:   PetscCall(DMDestroy(&dac));
470:   PetscCall(DMDestroy(&daf));
471:   PetscCall(VecDestroy(&ac));
472:   PetscCall(VecDestroy(&af));
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: PetscErrorCode da_test_RefineCoords3D(PetscInt mx, PetscInt my, PetscInt mz)
477: {
478:   DM          dac, daf;
479:   PetscViewer vv;
480:   Vec         ac, af;
481:   PetscInt    map_id, Mx, My, Mz;
482:   Mat         II, INTERP;
483:   Vec         scale;
484:   PetscBool   output = PETSC_FALSE;

486:   PetscFunctionBeginUser;
487:   PetscCall(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 */
488:                          1, /* stencil = 1 */ NULL, NULL, NULL, &dac));
489:   PetscCall(DMSetFromOptions(dac));
490:   PetscCall(DMSetUp(dac));

492:   PetscCall(DMRefine(dac, MPI_COMM_NULL, &daf));
493:   PetscCall(DMDAGetInfo(daf, 0, &Mx, &My, &Mz, 0, 0, 0, 0, 0, 0, 0, 0, 0));
494:   Mx--;
495:   My--;
496:   Mz--;

498:   PetscCall(DMDASetUniformCoordinates(dac, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));
499:   PetscCall(DMDASetUniformCoordinates(daf, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0));

501:   /* apply trilinear mappings */
502:   /*PetscCall(DAApplyTrilinearMapping(dac));*/
503:   /* apply conformal mappings */
504:   map_id = 0;
505:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-cmap", &map_id, NULL));
506:   if (map_id >= 1) PetscCall(DAApplyConformalMapping(dac, map_id));

508:   {
509:     DM  cdaf, cdac;
510:     Vec coordsc, coordsf;

512:     PetscCall(DMGetCoordinateDM(dac, &cdac));
513:     PetscCall(DMGetCoordinateDM(daf, &cdaf));

515:     PetscCall(DMGetCoordinates(dac, &coordsc));
516:     PetscCall(DMGetCoordinates(daf, &coordsf));

518:     PetscCall(DMCreateInterpolation(cdac, cdaf, &II, &scale));
519:     PetscCall(MatInterpolate(II, coordsc, coordsf));
520:     PetscCall(MatDestroy(&II));
521:     PetscCall(VecDestroy(&scale));
522:   }

524:   PetscCall(DMCreateInterpolation(dac, daf, &INTERP, NULL));

526:   PetscCall(DMCreateGlobalVector(dac, &ac));
527:   PetscCall(VecZeroEntries(ac));
528:   PetscCall(DADefineXLinearField3D(dac, ac));

530:   PetscCall(DMCreateGlobalVector(daf, &af));
531:   PetscCall(VecZeroEntries(af));

533:   PetscCall(MatMult(INTERP, ac, af));

535:   {
536:     Vec       afexact;
537:     PetscReal nrm;
538:     PetscInt  N;

540:     PetscCall(DMCreateGlobalVector(daf, &afexact));
541:     PetscCall(VecZeroEntries(afexact));
542:     PetscCall(DADefineXLinearField3D(daf, afexact));
543:     PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
544:     PetscCall(VecNorm(afexact, NORM_2, &nrm));
545:     PetscCall(VecGetSize(afexact, &N));
546:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "]=>[%" PetscInt_FMT " x %" PetscInt_FMT " x %" PetscInt_FMT "], interp err = %1.4e\n", mx, my, mz, Mx, My, Mz, (double)(nrm / PetscSqrtReal((PetscReal)N))));
547:     PetscCall(VecDestroy(&afexact));
548:   }

550:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-output", &output, NULL));
551:   if (output) {
552:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "dac_3D.vtr", &vv));
553:     PetscCall(VecView(ac, vv));
554:     PetscCall(PetscViewerDestroy(&vv));

556:     PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "daf_3D.vtr", &vv));
557:     PetscCall(VecView(af, vv));
558:     PetscCall(PetscViewerDestroy(&vv));
559:   }

561:   PetscCall(MatDestroy(&INTERP));
562:   PetscCall(DMDestroy(&dac));
563:   PetscCall(DMDestroy(&daf));
564:   PetscCall(VecDestroy(&ac));
565:   PetscCall(VecDestroy(&af));
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

569: int main(int argc, char **argv)
570: {
571:   PetscInt mx = 2, my = 2, mz = 2, l, nl, dim;

573:   PetscFunctionBeginUser;
574:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
575:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, 0));
576:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, 0));
577:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, 0));
578:   nl = 1;
579:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nl", &nl, 0));
580:   dim = 2;
581:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, 0));

583:   for (l = 0; l < nl; l++) {
584:     if (dim == 1) PetscCall(da_test_RefineCoords1D(mx));
585:     else if (dim == 2) PetscCall(da_test_RefineCoords2D(mx, my));
586:     else if (dim == 3) PetscCall(da_test_RefineCoords3D(mx, my, mz));
587:     mx = mx * 2;
588:     my = my * 2;
589:     mz = mz * 2;
590:   }
591:   PetscCall(PetscFinalize());
592:   return 0;
593: }

595: /*TEST

597:    test:
598:       suffix: 1d
599:       args: -mx 10 -nl 6 -dim 1

601:    test:
602:       suffix: 2d
603:       output_file: output/ex36_2d.out
604:       args: -mx 10 -my 10 -nl 6 -dim 2 -cmap {{0 1 2 3}}

606:    test:
607:       suffix: 2dp1
608:       nsize: 8
609:       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 4
610:       timeoutfactor: 2

612:    test:
613:       suffix: 2dp2
614:       nsize: 8
615:       args: -mx 10 -my 10 -nl 4 -dim 2 -cmap 3 -da_refine_x 3 -da_refine_y 1
616:       timeoutfactor: 2

618:    test:
619:       suffix: 3d
620:       args: -mx 5 -my 5 -mz 5 -nl 4 -dim 3 -cmap 3

622:    test:
623:       suffix: 3dp1
624:       nsize: 32
625:       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

627: TEST*/