Actual source code: ex36.c

  1: static char help[] = "Checks the functionality of DMGetInterpolation() on deformed grids.\n\n";

  3: #include <petscdm.h>
  4: #include <petscdmda.h>

  6: typedef struct _n_CCmplx CCmplx;
  7: struct _n_CCmplx {
  8:   PetscReal real;
  9:   PetscReal imag;
 10: };

 12: CCmplx CCmplxPow(CCmplx a, PetscReal n)
 13: {
 14:   CCmplx    b;
 15:   PetscReal r, theta;
 16:   r      = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
 17:   theta  = PetscAtan2Real(a.imag, a.real);
 18:   b.real = PetscPowReal(r, n) * PetscCosReal(n * theta);
 19:   b.imag = PetscPowReal(r, n) * PetscSinReal(n * theta);
 20:   return b;
 21: }
 22: CCmplx CCmplxExp(CCmplx a)
 23: {
 24:   CCmplx b;
 25:   b.real = PetscExpReal(a.real) * PetscCosReal(a.imag);
 26:   b.imag = PetscExpReal(a.real) * PetscSinReal(a.imag);
 27:   return b;
 28: }
 29: CCmplx CCmplxSqrt(CCmplx a)
 30: {
 31:   CCmplx    b;
 32:   PetscReal r, theta;
 33:   r      = PetscSqrtReal(a.real * a.real + a.imag * a.imag);
 34:   theta  = PetscAtan2Real(a.imag, a.real);
 35:   b.real = PetscSqrtReal(r) * PetscCosReal(0.5 * theta);
 36:   b.imag = PetscSqrtReal(r) * PetscSinReal(0.5 * theta);
 37:   return b;
 38: }
 39: CCmplx CCmplxAdd(CCmplx a, CCmplx c)
 40: {
 41:   CCmplx b;
 42:   b.real = a.real + c.real;
 43:   b.imag = a.imag + c.imag;
 44:   return b;
 45: }
 46: PetscScalar CCmplxRe(CCmplx a)
 47: {
 48:   return (PetscScalar)a.real;
 49: }
 50: PetscScalar CCmplxIm(CCmplx a)
 51: {
 52:   return (PetscScalar)a.imag;
 53: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

295:   for (k = sz; k < sz + nz; k++) {
296:     for (j = sy; j < sy + ny; j++) {
297:       for (i = sx; i < sx + nx; i++) {
298:         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 +
299:                       3.79 * XX[k][j][i].x * XX[k][j][i].y * XX[k][j][i].z;
300:       }
301:     }
302:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

446:     PetscCall(DMCreateGlobalVector(daf, &afexact));
447:     PetscCall(VecZeroEntries(afexact));
448:     PetscCall(DADefineXLinearField2D(daf, afexact));
449:     PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
450:     PetscCall(VecNorm(afexact, NORM_2, &nrm));
451:     PetscCall(VecGetSize(afexact, &N));
452:     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))));
453:     PetscCall(VecDestroy(&afexact));
454:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

539:     PetscCall(DMCreateGlobalVector(daf, &afexact));
540:     PetscCall(VecZeroEntries(afexact));
541:     PetscCall(DADefineXLinearField3D(daf, afexact));
542:     PetscCall(VecAXPY(afexact, -1.0, af)); /* af <= af - afinterp */
543:     PetscCall(VecNorm(afexact, NORM_2, &nrm));
544:     PetscCall(VecGetSize(afexact, &N));
545:     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))));
546:     PetscCall(VecDestroy(&afexact));
547:   }

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

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

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

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

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

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

594: /*TEST

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

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

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

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

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

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

626: TEST*/