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:   PetscFunctionBegin;
 24:   PetscCall(DMCreate(comm, dm));
 25:   PetscCall(DMSetType(*dm, DMPLEX));
 26:   PetscCall(DMSetFromOptions(*dm));
 27:   PetscCall(DMSetApplicationContext(*dm, user));
 28:   PetscCall(PetscObjectSetName((PetscObject)*dm, "Input Mesh"));
 29:   PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
 30:   PetscFunctionReturn(PETSC_SUCCESS);
 31: }

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

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

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

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

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

102:   if (options->transform) PetscCall(PetscPrintf(comm, "Using random transforms\n"));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

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

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

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

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

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

151:   PetscFunctionBegin;
152:   PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
153:   for (d = 0; d < spaceDim; ++d) {
154:     if (v0[d] != v0Ex[d]) {
155:       switch (spaceDim) {
156:       case 2:
157:         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]);
158:       case 3:
159:         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]);
160:       default:
161:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid space dimension %" PetscInt_FMT, spaceDim);
162:       }
163:     }
164:   }
165:   for (i = 0; i < spaceDim; ++i) {
166:     for (j = 0; j < spaceDim; ++j) {
167:       PetscCheck(RelativeError(J[i * spaceDim + j], JEx[i * spaceDim + j]) < 10 * PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid J[%" PetscInt_FMT ",%" PetscInt_FMT "]: %g != %g", i, j, (double)J[i * spaceDim + j], (double)JEx[i * spaceDim + j]);
168:       PetscCheck(RelativeError(invJ[i * spaceDim + j], invJEx[i * spaceDim + j]) < 10 * PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid invJ[%" PetscInt_FMT ",%" PetscInt_FMT "]: %g != %g", i, j, (double)invJ[i * spaceDim + j], (double)invJEx[i * spaceDim + j]);
169:     }
170:   }
171:   PetscCheck(RelativeError(detJ, detJEx) < 10 * PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid |J| = %g != %g diff %g", (double)detJ, (double)detJEx, (double)(detJ - detJEx));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

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

181:   PetscFunctionBegin;
182:   PetscCall(DMPlexComputeCellGeometryFVM(dm, cell, volEx ? &vol : NULL, centroidEx ? centroid : NULL, normalEx ? normal : NULL));
183:   for (d = 0; d < spaceDim; ++d) {
184:     if (centroidEx)
185:       PetscCheck(RelativeError(centroid[d], centroidEx[d]) < tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT ", Invalid centroid[%" PetscInt_FMT "]: %g != %g diff %g", cell, d, (double)centroid[d], (double)centroidEx[d], (double)(centroid[d] - centroidEx[d]));
186:     if (normalEx) PetscCheck(RelativeError(normal[d], normalEx[d]) < tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT ", Invalid normal[%" PetscInt_FMT "]: %g != %g", cell, d, (double)normal[d], (double)normalEx[d]);
187:   }
188:   if (volEx) PetscCheck(RelativeError(volEx, vol) < tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT ", Invalid volume = %g != %g diff %g", cell, (double)vol, (double)volEx, (double)(vol - volEx));
189:   PetscFunctionReturn(PETSC_SUCCESS);
190: }

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

200:   PetscFunctionBegin;
201:   PetscCall(DMGetDimension(dm, &dim));
202:   PetscCall(DMGetCoordinateDim(dm, &cdim));
203:   if (dim != cdim) PetscFunctionReturn(PETSC_SUCCESS);
204:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
205:   if (ct == DM_POLYTOPE_TRI_PRISM_TENSOR) PetscFunctionReturn(PETSC_SUCCESS);
206:   PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
207:   PetscCall(DMPlexGetCone(dm, cell, &cone));
208:   PetscCall(DMPlexGetConeOrientation(dm, cell, &ornt));
209:   for (f = 0; f < coneSize; ++f) {
210:     const PetscInt sgn = dim == 1 ? (f == 0 ? -1 : 1) : (ornt[f] < 0 ? -1 : 1);

212:     PetscCall(DMPlexComputeCellGeometryFVM(dm, cone[f], &area, NULL, normal));
213:     for (d = 0; d < cdim; ++d) integral[d] += sgn * area * normal[d];
214:   }
215:   for (d = 0; d < cdim; ++d)
216:     PetscCheck(PetscAbsReal(integral[d]) < tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " Surface integral for component %" PetscInt_FMT ": %g != 0. as it should be for a constant field", cell, d, (double)integral[d]);
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: 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[])
221: {
222:   const PetscInt *cone;
223:   PetscInt        coneSize, c;
224:   PetscInt        dim, depth, cdim;

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

249:     PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
250:     PetscCall(DMGetCoordinateSection(dm, &coordSection));
251:     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords));
252:     PetscCall(PetscMalloc2(coordSize, &origCoords, coordSize, &newCoords));
253:     PetscCall(PetscMalloc5(cdim, &v0ExT, cdim * cdim, &JExT, cdim * cdim, &invJExT, cdim, &centroidExT, cdim, &normalExT));
254:     PetscCall(PetscMalloc2(cdim, &faceCentroidExT, cdim, &faceNormalExT));
255:     for (c = 0; c < coordSize; ++c) origCoords[c] = coords[c];
256:     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords));
257:     numCorners = coordSize / cdim;

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

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

385:       PetscCall(ChangeCoordinates(dm, cdim, newCoords));
386:       if (v0Ex) PetscCall(CheckFEMGeometry(dm, 0, cdim, v0ExT, JExT, invJExT, detJExT));
387:       if (dim == depth && centroidEx) {
388:         PetscCall(CheckFVMGeometry(dm, cell, cdim, centroidExT, normalExT, volExT));
389:         PetscCall(CheckGaussLaw(dm, cell));
390:         if (faceCentroidEx) {
391:           PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
392:           PetscCall(DMPlexGetCone(dm, cell, &cone));
393:           for (c = 0; c < coneSize; ++c) {
394:             PetscInt off = c * cdim;

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

420:             PetscCall(CheckFVMGeometry(dm, cone[c], cdim, faceCentroidExT, faceNormalExT, faceVolExT));
421:           }
422:         }
423:       }
424:     }
425:     PetscCall(PetscRandomDestroy(&r));
426:     PetscCall(PetscRandomDestroy(&ang));
427:     PetscCall(PetscRandomDestroy(&ang2));
428:     PetscCall(PetscFree2(origCoords, newCoords));
429:     PetscCall(PetscFree5(v0ExT, JExT, invJExT, centroidExT, normalExT));
430:     PetscCall(PetscFree2(faceCentroidExT, faceNormalExT));
431:   }
432:   PetscFunctionReturn(PETSC_SUCCESS);
433: }

435: static PetscErrorCode TestTriangle(MPI_Comm comm, PetscBool transform)
436: {
437:   DM dm;

439:   PetscFunctionBegin;
440:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRIANGLE, &dm));
441:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
442:   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
443:   {
444:     PetscReal v0Ex[2]       = {-1.0, -1.0};
445:     PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
446:     PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
447:     PetscReal detJEx        = 1.0;
448:     PetscReal centroidEx[2] = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.)};
449:     PetscReal normalEx[2]   = {0.0, 0.0};
450:     PetscReal volEx         = 2.0;

452:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
453:   }
454:   /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (2.0) */
455:   {
456:     PetscScalar vertexCoords[9] = {-1.0, -1.0, 0.0, 1.0, -1.0, 0.0, -1.0, 1.0, 0.0};
457:     PetscReal   v0Ex[3]         = {-1.0, -1.0, 0.0};
458:     PetscReal   JEx[9]          = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
459:     PetscReal   invJEx[9]       = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
460:     PetscReal   detJEx          = 1.0;
461:     PetscReal   centroidEx[3]   = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), 0.0};
462:     PetscReal   normalEx[3]     = {0.0, 0.0, 1.0};
463:     PetscReal   volEx           = 2.0;

465:     PetscCall(ChangeCoordinates(dm, 3, vertexCoords));
466:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
467:   }
468:   /* Cleanup */
469:   PetscCall(DMDestroy(&dm));
470:   PetscFunctionReturn(PETSC_SUCCESS);
471: }

473: static PetscErrorCode TestQuadrilateral(MPI_Comm comm, PetscBool transform)
474: {
475:   DM dm;

477:   PetscFunctionBegin;
478:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_QUADRILATERAL, &dm));
479:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
480:   /* Check reference geometry: determinant is scaled by reference volume (2.0) */
481:   {
482:     PetscReal v0Ex[2]       = {-1.0, -1.0};
483:     PetscReal JEx[4]        = {1.0, 0.0, 0.0, 1.0};
484:     PetscReal invJEx[4]     = {1.0, 0.0, 0.0, 1.0};
485:     PetscReal detJEx        = 1.0;
486:     PetscReal centroidEx[2] = {0.0, 0.0};
487:     PetscReal normalEx[2]   = {0.0, 0.0};
488:     PetscReal volEx         = 4.0;

490:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
491:   }
492:   /* Move to 3D: Check reference geometry: determinant is scaled by reference volume (4.0) */
493:   {
494:     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};
495:     PetscReal   v0Ex[3]          = {-1.0, -1.0, 0.0};
496:     PetscReal   JEx[9]           = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
497:     PetscReal   invJEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
498:     PetscReal   detJEx           = 1.0;
499:     PetscReal   centroidEx[3]    = {0.0, 0.0, 0.0};
500:     PetscReal   normalEx[3]      = {0.0, 0.0, 1.0};
501:     PetscReal   volEx            = 4.0;

503:     PetscCall(ChangeCoordinates(dm, 3, vertexCoords));
504:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
505:   }
506:   /* Cleanup */
507:   PetscCall(DMDestroy(&dm));
508:   PetscFunctionReturn(PETSC_SUCCESS);
509: }

511: static PetscErrorCode TestTetrahedron(MPI_Comm comm, PetscBool transform)
512: {
513:   DM dm;

515:   PetscFunctionBegin;
516:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TETRAHEDRON, &dm));
517:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
518:   /* Check reference geometry: determinant is scaled by reference volume (4/3) */
519:   {
520:     PetscReal v0Ex[3]       = {-1.0, -1.0, -1.0};
521:     PetscReal JEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
522:     PetscReal invJEx[9]     = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
523:     PetscReal detJEx        = 1.0;
524:     PetscReal centroidEx[3] = {-0.5, -0.5, -0.5};
525:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
526:     PetscReal volEx         = (PetscReal)4.0 / (PetscReal)3.0;

528:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
529:   }
530:   /* Cleanup */
531:   PetscCall(DMDestroy(&dm));
532:   PetscFunctionReturn(PETSC_SUCCESS);
533: }

535: static PetscErrorCode TestHexahedron(MPI_Comm comm, PetscBool transform)
536: {
537:   DM dm;

539:   PetscFunctionBegin;
540:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm));
541:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
542:   /* Check reference geometry: determinant is scaled by reference volume 8.0 */
543:   {
544:     PetscReal v0Ex[3]       = {-1.0, -1.0, -1.0};
545:     PetscReal JEx[9]        = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
546:     PetscReal invJEx[9]     = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
547:     PetscReal detJEx        = 1.0;
548:     PetscReal centroidEx[3] = {0.0, 0.0, 0.0};
549:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
550:     PetscReal volEx         = 8.0;

552:     PetscCall(CheckCell(dm, 0, transform, v0Ex, JEx, invJEx, detJEx, centroidEx, normalEx, volEx, NULL, NULL, NULL));
553:   }
554:   /* Cleanup */
555:   PetscCall(DMDestroy(&dm));
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

559: static PetscErrorCode TestHexahedronCurved(MPI_Comm comm)
560: {
561:   DM          dm;
562:   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};

564:   PetscFunctionBegin;
565:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_HEXAHEDRON, &dm));
566:   PetscCall(ChangeCoordinates(dm, 3, coords));
567:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
568:   {
569:     PetscReal centroidEx[3] = {0.0, 0.0, 0.016803278688524603};
570:     PetscReal normalEx[3]   = {0.0, 0.0, 0.0};
571:     PetscReal volEx         = 8.1333333333333346;

573:     PetscCall(CheckCell(dm, 0, PETSC_FALSE, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, NULL, NULL, NULL));
574:   }
575:   PetscCall(DMDestroy(&dm));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: /* This wedge is a tensor product cell, rather than a normal wedge */
580: static PetscErrorCode TestWedge(MPI_Comm comm, PetscBool transform)
581: {
582:   DM dm;

584:   PetscFunctionBegin;
585:   PetscCall(DMPlexCreateReferenceCell(comm, DM_POLYTOPE_TRI_PRISM_TENSOR, &dm));
586:   PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
587:   /* Check reference geometry: determinant is scaled by reference volume 4.0 */
588:   {
589: #if 0
590:     /* FEM geometry is not functional for wedges */
591:     PetscReal v0Ex[3]   = {-1.0, -1.0, -1.0};
592:     PetscReal JEx[9]    = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
593:     PetscReal invJEx[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
594:     PetscReal detJEx    = 1.0;
595: #endif

597:     {
598:       PetscReal centroidEx[3]      = {-((PetscReal)1.) / ((PetscReal)3.), -((PetscReal)1.) / ((PetscReal)3.), 0.0};
599:       PetscReal normalEx[3]        = {0.0, 0.0, 0.0};
600:       PetscReal volEx              = 4.0;
601:       PetscReal faceVolEx[5]       = {2.0, 2.0, 4.0, PETSC_SQRT2 * 4.0, 4.0};
602:       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};
603:       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};

605:       PetscCall(CheckCell(dm, 0, transform, NULL, NULL, NULL, 0.0, centroidEx, normalEx, volEx, faceCentroidEx, faceNormalEx, faceVolEx));
606:     }
607:   }
608:   /* Cleanup */
609:   PetscCall(DMDestroy(&dm));
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

613: int main(int argc, char **argv)
614: {
615:   AppCtx user;

617:   PetscFunctionBeginUser;
618:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
619:   PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
620:   if (user.runType == RUN_REFERENCE) {
621:     PetscCall(TestTriangle(PETSC_COMM_SELF, user.transform));
622:     PetscCall(TestQuadrilateral(PETSC_COMM_SELF, user.transform));
623:     PetscCall(TestTetrahedron(PETSC_COMM_SELF, user.transform));
624:     PetscCall(TestHexahedron(PETSC_COMM_SELF, user.transform));
625:     PetscCall(TestWedge(PETSC_COMM_SELF, user.transform));
626:   } else if (user.runType == RUN_HEX_CURVED) {
627:     PetscCall(TestHexahedronCurved(PETSC_COMM_SELF));
628:   } else if (user.runType == RUN_FILE) {
629:     PetscInt dim, cStart, cEnd, c;

631:     PetscCall(DMGetDimension(user.dm, &dim));
632:     PetscCall(DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd));
633:     for (c = 0; c < cEnd - cStart; ++c) {
634:       PetscReal *v0       = user.v0 ? &user.v0[c * dim] : NULL;
635:       PetscReal *J        = user.J ? &user.J[c * dim * dim] : NULL;
636:       PetscReal *invJ     = user.invJ ? &user.invJ[c * dim * dim] : NULL;
637:       PetscReal  detJ     = user.detJ ? user.detJ[c] : 0.0;
638:       PetscReal *centroid = user.centroid ? &user.centroid[c * dim] : NULL;
639:       PetscReal *normal   = user.normal ? &user.normal[c * dim] : NULL;
640:       PetscReal  vol      = user.vol ? user.vol[c] : 0.0;

642:       PetscCall(CheckCell(user.dm, c + cStart, PETSC_FALSE, v0, J, invJ, detJ, centroid, normal, vol, NULL, NULL, NULL));
643:     }
644:     PetscCall(PetscFree4(user.v0, user.J, user.invJ, user.detJ));
645:     PetscCall(PetscFree(user.centroid));
646:     PetscCall(PetscFree(user.normal));
647:     PetscCall(PetscFree(user.vol));
648:     PetscCall(DMDestroy(&user.dm));
649:   } else if (user.runType == RUN_DISPLAY) {
650:     DM                 gdm, dmCell;
651:     Vec                cellgeom, facegeom;
652:     const PetscScalar *cgeom;
653:     PetscInt           dim, d, cStart, cEnd, cEndInterior, c;

655:     PetscCall(DMGetCoordinateDim(user.dm, &dim));
656:     PetscCall(DMPlexConstructGhostCells(user.dm, NULL, NULL, &gdm));
657:     if (gdm) {
658:       PetscCall(DMDestroy(&user.dm));
659:       user.dm = gdm;
660:     }
661:     PetscCall(DMPlexComputeGeometryFVM(user.dm, &cellgeom, &facegeom));
662:     PetscCall(DMPlexGetHeightStratum(user.dm, 0, &cStart, &cEnd));
663:     PetscCall(DMPlexGetGhostCellStratum(user.dm, &cEndInterior, NULL));
664:     if (cEndInterior >= 0) cEnd = cEndInterior;
665:     PetscCall(VecGetDM(cellgeom, &dmCell));
666:     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
667:     for (c = 0; c < cEnd - cStart; ++c) {
668:       PetscFVCellGeom *cg;

670:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
671:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %4" PetscInt_FMT ": Centroid (", c));
672:       for (d = 0; d < dim; ++d) {
673:         if (d > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
674:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "%12.2g", (double)cg->centroid[d]));
675:       }
676:       PetscCall(PetscPrintf(PETSC_COMM_SELF, ") Vol %12.2g\n", (double)cg->volume));
677:     }
678:     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
679:     PetscCall(VecDestroy(&cellgeom));
680:     PetscCall(VecDestroy(&facegeom));
681:     PetscCall(DMDestroy(&user.dm));
682:   }
683:   PetscCall(PetscFinalize());
684:   return 0;
685: }

687: /*TEST

689:   test:
690:     suffix: 1
691:     args: -dm_view ascii::ascii_info_detail
692:   test:
693:     suffix: 2
694:     args: -run_type hex_curved
695:   test:
696:     suffix: 3
697:     args: -transform
698:   test:
699:     suffix: 4
700:     requires: exodusii
701:     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
702:   test:
703:     suffix: 5
704:     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
705:   test:
706:     suffix: 6
707:     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
708: TEST*/