Actual source code: ex8.c

  1: static char help[] = "Tests for cell geometry\n\n";

  3: #include <petscdmplex.h>

  5: typedef enum {
  6:   RUN_REFERENCE,
  7:   RUN_HEX_CURVED,
  8:   RUN_FILE,
  9:   RUN_DISPLAY
 10: } RunType;

 12: typedef struct {
 13:   DM        dm;
 14:   RunType   runType;   /* Type of mesh to use */
 15:   PetscBool transform; /* Use random coordinate transformations */
 16:   /* Data for input meshes */
 17:   PetscReal *v0, *J, *invJ, *detJ;    /* FEM data */
 18:   PetscReal *centroid, *normal, *vol; /* FVM data */
 19: } AppCtx;

 21: static PetscErrorCode ReadMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 22: {
 23:   DMCreate(comm, dm);
 24:   DMSetType(*dm, DMPLEX);
 25:   DMSetFromOptions(*dm);
 26:   DMSetApplicationContext(*dm, user);
 27:   PetscObjectSetName((PetscObject)*dm, "Input Mesh");
 28:   DMViewFromOptions(*dm, NULL, "-dm_view");
 29:   return 0;
 30: }

 32: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 33: {
 34:   const char *runTypes[4] = {"reference", "hex_curved", "file", "display"};
 35:   PetscInt    run;

 38:   options->runType   = RUN_REFERENCE;
 39:   options->transform = PETSC_FALSE;

 41:   PetscOptionsBegin(comm, "", "Geometry Test Options", "DMPLEX");
 42:   run = options->runType;
 43:   PetscOptionsEList("-run_type", "The run type", "ex8.c", runTypes, 3, runTypes[options->runType], &run, NULL);
 44:   options->runType = (RunType)run;
 45:   PetscOptionsBool("-transform", "Use random transforms", "ex8.c", options->transform, &options->transform, NULL);

 47:   if (options->runType == RUN_FILE) {
 48:     PetscInt  dim, cStart, cEnd, numCells, n;
 49:     PetscBool flg, feFlg;

 51:     ReadMesh(PETSC_COMM_WORLD, options, &options->dm);
 52:     DMGetDimension(options->dm, &dim);
 53:     DMPlexGetHeightStratum(options->dm, 0, &cStart, &cEnd);
 54:     numCells = cEnd - cStart;
 55:     PetscMalloc4(numCells * dim, &options->v0, numCells * dim * dim, &options->J, numCells * dim * dim, &options->invJ, numCells, &options->detJ);
 56:     PetscMalloc1(numCells * dim, &options->centroid);
 57:     PetscMalloc1(numCells * dim, &options->normal);
 58:     PetscMalloc1(numCells, &options->vol);
 59:     n = numCells * dim;
 60:     PetscOptionsRealArray("-v0", "Input v0 for each cell", "ex8.c", options->v0, &n, &feFlg);
 62:     n = numCells * dim * dim;
 63:     PetscOptionsRealArray("-J", "Input Jacobian for each cell", "ex8.c", options->J, &n, &flg);
 65:     n = numCells * dim * dim;
 66:     PetscOptionsRealArray("-invJ", "Input inverse Jacobian for each cell", "ex8.c", options->invJ, &n, &flg);
 68:     n = numCells;
 69:     PetscOptionsRealArray("-detJ", "Input Jacobian determinant for each cell", "ex8.c", options->detJ, &n, &flg);
 71:     n = numCells * dim;
 72:     if (!feFlg) {
 73:       PetscFree4(options->v0, options->J, options->invJ, options->detJ);
 74:       options->v0 = options->J = options->invJ = options->detJ = NULL;
 75:     }
 76:     PetscOptionsRealArray("-centroid", "Input centroid for each cell", "ex8.c", options->centroid, &n, &flg);
 78:     if (!flg) {
 79:       PetscFree(options->centroid);
 80:       options->centroid = NULL;
 81:     }
 82:     n = numCells * dim;
 83:     PetscOptionsRealArray("-normal", "Input normal for each cell", "ex8.c", options->normal, &n, &flg);
 85:     if (!flg) {
 86:       PetscFree(options->normal);
 87:       options->normal = NULL;
 88:     }
 89:     n = numCells;
 90:     PetscOptionsRealArray("-vol", "Input volume for each cell", "ex8.c", options->vol, &n, &flg);
 92:     if (!flg) {
 93:       PetscFree(options->vol);
 94:       options->vol = NULL;
 95:     }
 96:   } else if (options->runType == RUN_DISPLAY) {
 97:     ReadMesh(PETSC_COMM_WORLD, options, &options->dm);
 98:   }
 99:   PetscOptionsEnd();

101:   if (options->transform) PetscPrintf(comm, "Using random transforms\n");
102:   return 0;
103: }

105: static PetscErrorCode ChangeCoordinates(DM dm, PetscInt spaceDim, PetscScalar vertexCoords[])
106: {
107:   PetscSection coordSection;
108:   Vec          coordinates;
109:   PetscScalar *coords;
110:   PetscInt     vStart, vEnd, v, d, coordSize;

112:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
113:   DMGetCoordinateSection(dm, &coordSection);
114:   PetscSectionSetNumFields(coordSection, 1);
115:   PetscSectionSetFieldComponents(coordSection, 0, spaceDim);
116:   PetscSectionSetChart(coordSection, vStart, vEnd);
117:   for (v = vStart; v < vEnd; ++v) {
118:     PetscSectionSetDof(coordSection, v, spaceDim);
119:     PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);
120:   }
121:   PetscSectionSetUp(coordSection);
122:   PetscSectionGetStorageSize(coordSection, &coordSize);
123:   VecCreate(PETSC_COMM_SELF, &coordinates);
124:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
125:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
126:   VecSetFromOptions(coordinates);
127:   VecGetArray(coordinates, &coords);
128:   for (v = vStart; v < vEnd; ++v) {
129:     PetscInt off;

131:     PetscSectionGetOffset(coordSection, v, &off);
132:     for (d = 0; d < spaceDim; ++d) coords[off + d] = vertexCoords[(v - vStart) * spaceDim + d];
133:   }
134:   VecRestoreArray(coordinates, &coords);
135:   DMSetCoordinateDim(dm, spaceDim);
136:   DMSetCoordinatesLocal(dm, coordinates);
137:   VecDestroy(&coordinates);
138:   DMViewFromOptions(dm, NULL, "-dm_view");
139:   return 0;
140: }

142: #define RelativeError(a, b) PetscAbs(a - b) / (1.0 + PetscMax(PetscAbs(a), PetscAbs(b)))

144: static PetscErrorCode CheckFEMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx)
145: {
146:   PetscReal v0[3], J[9], invJ[9], detJ;
147:   PetscInt  d, i, j;

149:   DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ);
150:   for (d = 0; d < spaceDim; ++d) {
151:     if (v0[d] != v0Ex[d]) {
152:       switch (spaceDim) {
153:       case 2:
154:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g) != (%g, %g)", (double)v0[0], (double)v0[1], (double)v0Ex[0], (double)v0Ex[1]);
155:       case 3:
156:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid v0 (%g, %g, %g) != (%g, %g, %g)", (double)v0[0], (double)v0[1], (double)v0[2], (double)v0Ex[0], (double)v0Ex[1], (double)v0Ex[2]);
157:       default:
158:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %" PetscInt_FMT, spaceDim);
159:       }
160:     }
161:   }
162:   for (i = 0; i < spaceDim; ++i) {
163:     for (j = 0; j < spaceDim; ++j) {
166:     }
167:   }
169:   return 0;
170: }

172: static PetscErrorCode CheckFVMGeometry(DM dm, PetscInt cell, PetscInt spaceDim, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx)
173: {
174:   PetscReal tol = PetscMax(10 * PETSC_SMALL, 1e-10);
175:   PetscReal centroid[3], normal[3], vol;
176:   PetscInt  d;

178:   DMPlexComputeCellGeometryFVM(dm, cell, volEx ? &vol : NULL, centroidEx ? centroid : NULL, normalEx ? normal : NULL);
179:   for (d = 0; d < spaceDim; ++d) {
180:     if (centroidEx)
183:   }
185:   return 0;
186: }

188: static PetscErrorCode CheckGaussLaw(DM dm, PetscInt cell)
189: {
190:   DMPolytopeType  ct;
191:   PetscReal       tol = PetscMax(10 * PETSC_SMALL, 1e-10);
192:   PetscReal       normal[3], integral[3] = {0., 0., 0.}, area;
193:   const PetscInt *cone, *ornt;
194:   PetscInt        coneSize, f, dim, cdim, d;

196:   DMGetDimension(dm, &dim);
197:   DMGetCoordinateDim(dm, &cdim);
198:   if (dim != cdim) return 0;
199:   DMPlexGetCellType(dm, cell, &ct);
200:   if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) return 0;
201:   DMPlexGetConeSize(dm, cell, &coneSize);
202:   DMPlexGetCone(dm, cell, &cone);
203:   DMPlexGetConeOrientation(dm, cell, &ornt);
204:   for (f = 0; f < coneSize; ++f) {
205:     const PetscInt sgn = dim == 1 ? (f == 0 ? -1 : 1) : (ornt[f] < 0 ? -1 : 1);

207:     DMPlexComputeCellGeometryFVM(dm, cone[f], &area, NULL, normal);
208:     for (d = 0; d < cdim; ++d) integral[d] += sgn * area * normal[d];
209:   }
210:   for (d = 0; d < cdim; ++d)
212:   return 0;
213: }

215: static PetscErrorCode CheckCell(DM dm, PetscInt cell, PetscBool transform, PetscReal v0Ex[], PetscReal JEx[], PetscReal invJEx[], PetscReal detJEx, PetscReal centroidEx[], PetscReal normalEx[], PetscReal volEx, PetscReal faceCentroidEx[], PetscReal faceNormalEx[], PetscReal faceVolEx[])
216: {
217:   const PetscInt *cone;
218:   PetscInt        coneSize, c;
219:   PetscInt        dim, depth, cdim;

221:   DMPlexGetDepth(dm, &depth);
222:   DMGetDimension(dm, &dim);
223:   DMGetCoordinateDim(dm, &cdim);
224:   if (v0Ex) CheckFEMGeometry(dm, cell, cdim, v0Ex, JEx, invJEx, detJEx);
225:   if (dim == depth && centroidEx) {
226:     CheckFVMGeometry(dm, cell, cdim, centroidEx, normalEx, volEx);
227:     CheckGaussLaw(dm, cell);
228:     if (faceCentroidEx) {
229:       DMPlexGetConeSize(dm, cell, &coneSize);
230:       DMPlexGetCone(dm, cell, &cone);
231:       for (c = 0; c < coneSize; ++c) CheckFVMGeometry(dm, cone[c], dim, &faceCentroidEx[c * dim], &faceNormalEx[c * dim], faceVolEx[c]);
232:     }
233:   }
234:   if (transform) {
235:     Vec          coordinates;
236:     PetscSection coordSection;
237:     PetscScalar *coords = NULL, *origCoords, *newCoords;
238:     PetscReal   *v0ExT, *JExT, *invJExT, detJExT = 0, *centroidExT, *normalExT, volExT = 0;
239:     PetscReal   *faceCentroidExT, *faceNormalExT, faceVolExT;
240:     PetscRandom  r, ang, ang2;
241:     PetscInt     coordSize, numCorners, t;

243:     DMGetCoordinatesLocal(dm, &coordinates);
244:     DMGetCoordinateSection(dm, &coordSection);
245:     DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
246:     PetscMalloc2(coordSize, &origCoords, coordSize, &newCoords);
247:     PetscMalloc5(cdim, &v0ExT, cdim * cdim, &JExT, cdim * cdim, &invJExT, cdim, &centroidExT, cdim, &normalExT);
248:     PetscMalloc2(cdim, &faceCentroidExT, cdim, &faceNormalExT);
249:     for (c = 0; c < coordSize; ++c) origCoords[c] = coords[c];
250:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
251:     numCorners = coordSize / cdim;

253:     PetscRandomCreate(PETSC_COMM_SELF, &r);
254:     PetscRandomSetFromOptions(r);
255:     PetscRandomSetInterval(r, 0.0, 10.0);
256:     PetscRandomCreate(PETSC_COMM_SELF, &ang);
257:     PetscRandomSetFromOptions(ang);
258:     PetscRandomSetInterval(ang, 0.0, 2 * PETSC_PI);
259:     PetscRandomCreate(PETSC_COMM_SELF, &ang2);
260:     PetscRandomSetFromOptions(ang2);
261:     PetscRandomSetInterval(ang2, 0.0, PETSC_PI);
262:     for (t = 0; t < 1; ++t) {
263:       PetscScalar trans[3];
264:       PetscReal   R[9], rot[3], rotM[9];
265:       PetscReal   scale, phi, theta, psi = 0.0, norm;
266:       PetscInt    d, e, f, p;

268:       for (c = 0; c < coordSize; ++c) newCoords[c] = origCoords[c];
269:       PetscRandomGetValueReal(r, &scale);
270:       PetscRandomGetValueReal(ang, &phi);
271:       PetscRandomGetValueReal(ang2, &theta);
272:       for (d = 0; d < cdim; ++d) PetscRandomGetValue(r, &trans[d]);
273:       switch (cdim) {
274:       case 2:
275:         R[0] = PetscCosReal(phi);
276:         R[1] = -PetscSinReal(phi);
277:         R[2] = PetscSinReal(phi);
278:         R[3] = PetscCosReal(phi);
279:         break;
280:       case 3: {
281:         const PetscReal ct = PetscCosReal(theta), st = PetscSinReal(theta);
282:         const PetscReal cp = PetscCosReal(phi), sp = PetscSinReal(phi);
283:         const PetscReal cs = PetscCosReal(psi), ss = PetscSinReal(psi);
284:         R[0] = ct * cs;
285:         R[1] = sp * st * cs - cp * ss;
286:         R[2] = sp * ss + cp * st * cs;
287:         R[3] = ct * ss;
288:         R[4] = cp * cs + sp * st * ss;
289:         R[5] = cp * st * ss - sp * cs;
290:         R[6] = -st;
291:         R[7] = sp * ct;
292:         R[8] = cp * ct;
293:         break;
294:       }
295:       default:
296:         SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid coordinate dimension %" PetscInt_FMT, cdim);
297:       }
298:       if (v0Ex) {
299:         detJExT = detJEx;
300:         for (d = 0; d < cdim; ++d) {
301:           v0ExT[d] = v0Ex[d];
302:           for (e = 0; e < cdim; ++e) {
303:             JExT[d * cdim + e]    = JEx[d * cdim + e];
304:             invJExT[d * cdim + e] = invJEx[d * cdim + e];
305:           }
306:         }
307:         for (d = 0; d < cdim; ++d) {
308:           v0ExT[d] *= scale;
309:           v0ExT[d] += PetscRealPart(trans[d]);
310:           /* Only scale dimensions in the manifold */
311:           for (e = 0; e < dim; ++e) {
312:             JExT[d * cdim + e] *= scale;
313:             invJExT[d * cdim + e] /= scale;
314:           }
315:           if (d < dim) detJExT *= scale;
316:         }
317:         /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */
318:         for (d = 0; d < cdim; ++d) {
319:           for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * v0ExT[e];
320:         }
321:         for (d = 0; d < cdim; ++d) v0ExT[d] = rot[d];
322:         for (d = 0; d < cdim; ++d) {
323:           for (e = 0; e < cdim; ++e) {
324:             for (f = 0, rotM[d * cdim + e] = 0.0; f < cdim; ++f) rotM[d * cdim + e] += R[d * cdim + f] * JExT[f * cdim + e];
325:           }
326:         }
327:         for (d = 0; d < cdim; ++d) {
328:           for (e = 0; e < cdim; ++e) JExT[d * cdim + e] = rotM[d * cdim + e];
329:         }
330:         for (d = 0; d < cdim; ++d) {
331:           for (e = 0; e < cdim; ++e) {
332:             for (f = 0, rotM[d * cdim + e] = 0.0; f < cdim; ++f) rotM[d * cdim + e] += invJExT[d * cdim + f] * R[e * cdim + f];
333:           }
334:         }
335:         for (d = 0; d < cdim; ++d) {
336:           for (e = 0; e < cdim; ++e) invJExT[d * cdim + e] = rotM[d * cdim + e];
337:         }
338:       }
339:       if (centroidEx) {
340:         volExT = volEx;
341:         for (d = 0; d < cdim; ++d) {
342:           centroidExT[d] = centroidEx[d];
343:           normalExT[d]   = normalEx[d];
344:         }
345:         for (d = 0; d < cdim; ++d) {
346:           centroidExT[d] *= scale;
347:           centroidExT[d] += PetscRealPart(trans[d]);
348:           normalExT[d] /= scale;
349:           /* Only scale dimensions in the manifold */
350:           if (d < dim) volExT *= scale;
351:         }
352:         /* Do scaling and translation before rotation, so that we can leave out the normal dimension for lower dimensional manifolds */
353:         for (d = 0; d < cdim; ++d) {
354:           for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * centroidExT[e];
355:         }
356:         for (d = 0; d < cdim; ++d) centroidExT[d] = rot[d];
357:         for (d = 0; d < cdim; ++d) {
358:           for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * normalExT[e];
359:         }
360:         for (d = 0; d < cdim; ++d) normalExT[d] = rot[d];
361:         for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(normalExT[d]);
362:         norm = PetscSqrtReal(norm);
363:         if (norm != 0.)
364:           for (d = 0; d < cdim; ++d) normalExT[d] /= norm;
365:       }
366:       for (d = 0; d < cdim; ++d) {
367:         for (p = 0; p < numCorners; ++p) {
368:           newCoords[p * cdim + d] *= scale;
369:           newCoords[p * cdim + d] += trans[d];
370:         }
371:       }
372:       for (p = 0; p < numCorners; ++p) {
373:         for (d = 0; d < cdim; ++d) {
374:           for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * PetscRealPart(newCoords[p * cdim + e]);
375:         }
376:         for (d = 0; d < cdim; ++d) newCoords[p * cdim + d] = rot[d];
377:       }

379:       ChangeCoordinates(dm, cdim, newCoords);
380:       if (v0Ex) CheckFEMGeometry(dm, 0, cdim, v0ExT, JExT, invJExT, detJExT);
381:       if (dim == depth && centroidEx) {
382:         CheckFVMGeometry(dm, cell, cdim, centroidExT, normalExT, volExT);
383:         CheckGaussLaw(dm, cell);
384:         if (faceCentroidEx) {
385:           DMPlexGetConeSize(dm, cell, &coneSize);
386:           DMPlexGetCone(dm, cell, &cone);
387:           for (c = 0; c < coneSize; ++c) {
388:             PetscInt off = c * cdim;

390:             faceVolExT = faceVolEx[c];
391:             for (d = 0; d < cdim; ++d) {
392:               faceCentroidExT[d] = faceCentroidEx[off + d];
393:               faceNormalExT[d]   = faceNormalEx[off + d];
394:             }
395:             for (d = 0; d < cdim; ++d) {
396:               faceCentroidExT[d] *= scale;
397:               faceCentroidExT[d] += PetscRealPart(trans[d]);
398:               faceNormalExT[d] /= scale;
399:               /* Only scale dimensions in the manifold */
400:               if (d < dim - 1) faceVolExT *= scale;
401:             }
402:             for (d = 0; d < cdim; ++d) {
403:               for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * faceCentroidExT[e];
404:             }
405:             for (d = 0; d < cdim; ++d) faceCentroidExT[d] = rot[d];
406:             for (d = 0; d < cdim; ++d) {
407:               for (e = 0, rot[d] = 0.0; e < cdim; ++e) rot[d] += R[d * cdim + e] * faceNormalExT[e];
408:             }
409:             for (d = 0; d < cdim; ++d) faceNormalExT[d] = rot[d];
410:             for (d = 0, norm = 0.0; d < cdim; ++d) norm += PetscSqr(faceNormalExT[d]);
411:             norm = PetscSqrtReal(norm);
412:             for (d = 0; d < cdim; ++d) faceNormalExT[d] /= norm;

414:             CheckFVMGeometry(dm, cone[c], cdim, faceCentroidExT, faceNormalExT, faceVolExT);
415:           }
416:         }
417:       }
418:     }
419:     PetscRandomDestroy(&r);
420:     PetscRandomDestroy(&ang);
421:     PetscRandomDestroy(&ang2);
422:     PetscFree2(origCoords, newCoords);
423:     PetscFree5(v0ExT, JExT, invJExT, centroidExT, normalExT);
424:     PetscFree2(faceCentroidExT, faceNormalExT);
425:   }
426:   return 0;
427: }

429: static PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool transform)
430: {
431:   DM dm;

433:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRIANGLE, &dm);
434:   DMViewFromOptions(dm, NULL, "-dm_view");
435:   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
436:   {
437:     PetscReal v0Ex[2]       = {-1.0, -1.0};
438:     PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
439:     PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
440:     PetscReal detJEx        = 1.0;
441:     PetscReal centroidEx[2] = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.)};
442:     PetscReal normalEx[2]   = {0.0, 0.0};
443:     PetscReal volEx         = 2.0;

445:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
446:   }
447:   /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */
448:   {
449:     PetscScalar vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0};
450:     PetscReal   v0Ex[3]         = {-1.0, -1.0, 0.0};
451:     PetscReal   JEx[9]          = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
452:     PetscReal   invJEx[9]       = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
453:     PetscReal   detJEx          = 1.0;
454:     PetscReal   centroidEx[3]   = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), 0.0};
455:     PetscReal   normalEx[3]     = {0.0, 0.0, 1.0};
456:     PetscReal   volEx           = 2.0;

458:     ChangeCoordinates(dm, 3, vertexCoords);
459:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
460:   }
461:   /* Cleanup */
462:   DMDestroy(&dm);
463:   return 0;
464: }

466: static PetscErrorCode TestQuadrilateral(MPI_Comm comm, PetscBool transform)
467: {
468:   DM dm;

470:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_QUADRILATERAL, &dm);
471:   DMViewFromOptions(dm, NULL, "-dm_view");
472:   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
473:   {
474:     PetscReal v0Ex[2]       = {-1.0, -1.0};
475:     PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
476:     PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
477:     PetscReal detJEx        = 1.0;
478:     PetscReal centroidEx[2] = {0.0, 0.0};
479:     PetscReal normalEx[2]   = {0.0, 0.0};
480:     PetscReal volEx         = 4.0;

482:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
483:   }
484:   /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (4.0) */
485:   {
486:     PetscScalar vertexCoords[12] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, 1.0, 0.0, -1.0, 1.0, 0.0};
487:     PetscReal   v0Ex[3]          = {-1.0, -1.0, 0.0};
488:     PetscReal   JEx[9]           = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
489:     PetscReal   invJEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
490:     PetscReal   detJEx           = 1.0;
491:     PetscReal   centroidEx[3]    = {0.0, 0.0, 0.0};
492:     PetscReal   normalEx[3]      = {0.0, 0.0, 1.0};
493:     PetscReal   volEx            = 4.0;

495:     ChangeCoordinates(dm, 3, vertexCoords);
496:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
497:   }
498:   /* Cleanup */
499:   DMDestroy(&dm);
500:   return 0;
501: }

503: static PetscErrorCode TestTetrahedron(MPI_Comm comm, PetscBool transform)
504: {
505:   DM dm;

507:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TETRAHEDRON, &dm);
508:   DMViewFromOptions(dm, NULL, "-dm_view");
509:   /* Check reference geometry: determinant is scaled by reference volume (4/3) */
510:   {
511:     PetscReal v0Ex[3]       = {-1.0, -1.0, -1.0};
512:     PetscReal JEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
513:     PetscReal invJEx[9]     = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
514:     PetscReal detJEx        = 1.0;
515:     PetscReal centroidEx[3] = {-0.5, -0.5, -0.5};
516:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
517:     PetscReal volEx         = (PetscReal)4.0 / (PetscReal)3.0;

519:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
520:   }
521:   /* Cleanup */
522:   DMDestroy(&dm);
523:   return 0;
524: }

526: static PetscErrorCode TestHexahedron(MPI_Comm comm, PetscBool transform)
527: {
528:   DM dm;

530:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm);
531:   DMViewFromOptions(dm, NULL, "-dm_view");
532:   /* Check reference geometry: determinant is scaled by reference volume 8.0 */
533:   {
534:     PetscReal v0Ex[3]       = {-1.0, -1.0, -1.0};
535:     PetscReal JEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
536:     PetscReal invJEx[9]     = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
537:     PetscReal detJEx        = 1.0;
538:     PetscReal centroidEx[3] = {0.0, 0.0, 0.0};
539:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
540:     PetscReal volEx         = 8.0;

542:     CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL);
543:   }
544:   /* Cleanup */
545:   DMDestroy(&dm);
546:   return 0;
547: }

549: static PetscErrorCode TestHexahedronCurved(MPI_Comm comm)
550: {
551:   DM          dm;
552:   PetscScalar coords[24] = {-1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.1, 1.0, -1.0, 1.0, 1.0, 1.0, 1.1, -1.0, 1.0, 1.0};

554:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm);
555:   ChangeCoordinates(dm, 3, coords);
556:   DMViewFromOptions(dm, NULL, "-dm_view");
557:   {
558:     PetscReal centroidEx[3] = {0.0, 0.0, 0.016803278688524603};
559:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
560:     PetscReal volEx         = 8.1333333333333346;

562:     CheckCell(dm, 0, PETSC_FALSE, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, NULL, NULL, NULL);
563:   }
564:   DMDestroy(&dm);
565:   return 0;
566: }

568: /* This wedge is a tensor product cell, rather than a normal wedge */
569: static PetscErrorCode TestWedge(MPI_Comm comm, PetscBool transform)
570: {
571:   DM dm;

573:   DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRI_PRISM_TENSOR, &dm);
574:   DMViewFromOptions(dm, NULL, "-dm_view");
575:   /* Check reference geometry: determinant is scaled by reference volume 4.0 */
576:   {
577: #if 0
578:     /* FEM geometry is not functional for wedges */
579:     PetscReal v0Ex[3]   = {-1.0, -1.0, -1.0};
580:     PetscReal JEx[9]    = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
581:     PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
582:     PetscReal detJEx    = 1.0;
583: #endif

585:     {
586:       PetscReal centroidEx[3]      = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), 0.0};
587:       PetscReal normalEx[3]        = {0.0, 0.0, 0.0};
588:       PetscReal volEx              = 4.0;
589:       PetscReal faceVolEx[5]       = {2.0, 2.0, 4.0, PETSC_SQRT2 * 4.0, 4.0};
590:       PetscReal faceNormalEx[15]   = {0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, -1.0, 0.0, PETSC_SQRT2 / 2.0, PETSC_SQRT2 / 2.0, 0.0, -1.0, 0.0, 0.0};
591:       PetscReal faceCentroidEx[15] = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), -1.0, -((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), 1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0};

593:       CheckCell(dm, 0, transform, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, faceCentroidEx, faceNormalEx, faceVolEx);
594:     }
595:   }
596:   /* Cleanup */
597:   DMDestroy(&dm);
598:   return 0;
599: }

601: int main(int argc, char **argv)
602: {
603:   AppCtx user;

606:   PetscInitialize(&argc, &argv, NULL, help);
607:   ProcessOptions(PETSC_COMM_WORLD, &user);
608:   if (user.runType == RUN_REFERENCE) {
609:     TestTriangle(PETSC_COMM_SELF, user.transform);
610:     TestQuadrilateral(PETSC_COMM_SELF, user.transform);
611:     TestTetrahedron(PETSC_COMM_SELF, user.transform);
612:     TestHexahedron(PETSC_COMM_SELF, user.transform);
613:     TestWedge(PETSC_COMM_SELF, user.transform);
614:   } else if (user.runType == RUN_HEX_CURVED) {
615:     TestHexahedronCurved(PETSC_COMM_SELF);
616:   } else if (user.runType == RUN_FILE) {
617:     PetscInt dim, cStart, cEnd, c;

619:     DMGetDimension(user.dm, &dim);
620:     DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd);
621:     for (c = 0; c < cEnd - cStart; ++c) {
622:       PetscReal *v0       = user.v0 ? &user.v0[c * dim] : NULL;
623:       PetscReal *J        = user.J ? &user.J[c * dim * dim] : NULL;
624:       PetscReal *invJ     = user.invJ ? &user.invJ[c * dim * dim] : NULL;
625:       PetscReal  detJ     = user.detJ ? user.detJ[c] : 0.0;
626:       PetscReal *centroid = user.centroid ? &user.centroid[c * dim] : NULL;
627:       PetscReal *normal   = user.normal ? &user.normal[c * dim] : NULL;
628:       PetscReal  vol      = user.vol ? user.vol[c] : 0.0;

630:       CheckCell(user.dm, c + cStart, PETSC_FALSE, v0, J, invJ, detJ, centroid, normal, vol, NULL, NULL, NULL);
631:     }
632:     PetscFree4(user.v0, user.J, user.invJ, user.detJ);
633:     PetscFree(user.centroid);
634:     PetscFree(user.normal);
635:     PetscFree(user.vol);
636:     DMDestroy(&user.dm);
637:   } else if (user.runType == RUN_DISPLAY) {
638:     DM                 gdm, dmCell;
639:     Vec                cellgeom, facegeom;
640:     const PetscScalar *cgeom;
641:     PetscInt           dim, d, cStart, cEnd, cEndInterior, c;

643:     DMGetCoordinateDim(user.dm, &dim);
644:     DMPlexConstructGhostCells(user.dm, NULL, NULL, &gdm);
645:     if (gdm) {
646:       DMDestroy(&user.dm);
647:       user.dm = gdm;
648:     }
649:     DMPlexComputeGeometryFVM(user.dm, &cellgeom, &facegeom);
650:     DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd);
651:     DMPlexGetGhostCellStratum(user.dm, &cEndInterior, NULL);
652:     if (cEndInterior >= 0) cEnd = cEndInterior;
653:     VecGetDM(cellgeom, &dmCell);
654:     VecGetArrayRead(cellgeom, &cgeom);
655:     for (c = 0; c < cEnd - cStart; ++c) {
656:       PetscFVCellGeom *cg;

658:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
659:       PetscPrintf(PETSC_COMM_SELF, "Cell %4" PetscInt_FMT ": Centroid (", c);
660:       for (d = 0; d < dim; ++d) {
661:         if (d > 0) PetscPrintf(PETSC_COMM_SELF, ", ");
662:         PetscPrintf(PETSC_COMM_SELF, "%12.2g", (double)cg->centroid[d]);
663:       }
664:       PetscPrintf(PETSC_COMM_SELF, ") Vol %12.2g\n", (double)cg->volume);
665:     }
666:     VecRestoreArrayRead(cellgeom, &cgeom);
667:     VecDestroy(&cellgeom);
668:     VecDestroy(&facegeom);
669:     DMDestroy(&user.dm);
670:   }
671:   PetscFinalize();
672:   return 0;
673: }

675: /*TEST

677:   test:
678:     suffix: 1
679:     args: -dm_view ascii::ascii_info_detail
680:   test:
681:     suffix: 2
682:     args: -run_type hex_curved
683:   test:
684:     suffix: 3
685:     args: -transform
686:   test:
687:     suffix: 4
688:     requires: exodusii
689:     args: -run_type file -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/simpleblock-100.exo -dm_view ascii::ascii_info_detail -v0 -1.5,-0.5,0.5,-0.5,-0.5,0.5,0.5,-0.5,0.5 -J 0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0,0.0,0.0,0.5,0.0,0.5,0.0,-0.5,0.0,0.0 -invJ 0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0,0.0,0.0,-2.0,0.0,2.0,0.0,2.0,0.0,0.0 -detJ 0.125,0.125,0.125 -centroid -1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0 -normal 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 -vol 1.0,1.0,1.0
690:   test:
691:     suffix: 5
692:     args: -run_type file -dm_plex_dim 3 -dm_plex_simplex 0 -dm_plex_box_faces 3,1,1 -dm_plex_box_lower -1.5,-0.5,-0.5 -dm_plex_box_upper 1.5,0.5,0.5 -dm_view ascii::ascii_info_detail -centroid -1.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0 -normal 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 -vol 1.0,1.0,1.0
693:   test:
694:     suffix: 6
695:     args: -run_type file -dm_plex_dim 1 -dm_plex_simplex 0 -dm_plex_box_faces 3 -dm_plex_box_lower -1.5 -dm_plex_box_upper 1.5 -dm_view ascii::ascii_info_detail -centroid -1.0,0.0,1.0 -vol 1.0,1.0,1.0
696: TEST*/