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, ¢roidExT, 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*/