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