Actual source code: plexgeometry.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/petscfeimpl.h>
3: #include <petscblaslapack.h>
4: #include <petsctime.h>
6: const char *const DMPlexCoordMaps[] = {"none", "rotate", "shear", "flare", "annulus", "shell", "sinusoid", "torus", "unknown", "DMPlexCoordMap", "DM_COORD_MAP_", NULL};
8: /*@
9: DMPlexFindVertices - Try to find DAG points based on their coordinates.
11: Not Collective (provided `DMGetCoordinatesLocalSetUp()` has been already called)
13: Input Parameters:
14: + dm - The `DMPLEX` object
15: . coordinates - The `Vec` of coordinates of the sought points
16: - eps - The tolerance or `PETSC_DEFAULT`
18: Output Parameter:
19: . points - The `IS` of found DAG points or -1
21: Level: intermediate
23: Notes:
24: The length of `Vec` coordinates must be npoints * dim where dim is the spatial dimension returned by `DMGetCoordinateDim()` and npoints is the number of sought points.
26: The output `IS` is living on `PETSC_COMM_SELF` and its length is npoints.
27: Each rank does the search independently.
28: If this rank's local `DMPLEX` portion contains the DAG point corresponding to the i-th tuple of coordinates, the i-th entry of the output `IS` is set to that DAG point, otherwise to -1.
30: The output `IS` must be destroyed by user.
32: The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.
34: Complexity of this function is currently O(mn) with m number of vertices to find and n number of vertices in the local mesh. This could probably be improved if needed.
36: .seealso: `DMPLEX`, `DMPlexCreate()`, `DMGetCoordinatesLocal()`
37: @*/
38: PetscErrorCode DMPlexFindVertices(DM dm, Vec coordinates, PetscReal eps, IS *points)
39: {
40: PetscInt c, cdim, i, j, o, p, vStart, vEnd;
41: PetscInt npoints;
42: const PetscScalar *coord;
43: Vec allCoordsVec;
44: const PetscScalar *allCoords;
45: PetscInt *dagPoints;
47: PetscFunctionBegin;
48: if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
49: PetscCall(DMGetCoordinateDim(dm, &cdim));
50: {
51: PetscInt n;
53: PetscCall(VecGetLocalSize(coordinates, &n));
54: PetscCheck(n % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Given coordinates Vec has local length %" PetscInt_FMT " not divisible by coordinate dimension %" PetscInt_FMT " of given DM", n, cdim);
55: npoints = n / cdim;
56: }
57: PetscCall(DMGetCoordinatesLocal(dm, &allCoordsVec));
58: PetscCall(VecGetArrayRead(allCoordsVec, &allCoords));
59: PetscCall(VecGetArrayRead(coordinates, &coord));
60: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
61: if (PetscDefined(USE_DEBUG)) {
62: /* check coordinate section is consistent with DM dimension */
63: PetscSection cs;
64: PetscInt ndof;
66: PetscCall(DMGetCoordinateSection(dm, &cs));
67: for (p = vStart; p < vEnd; p++) {
68: PetscCall(PetscSectionGetDof(cs, p, &ndof));
69: PetscCheck(ndof == cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %" PetscInt_FMT ": ndof = %" PetscInt_FMT " != %" PetscInt_FMT " = cdim", p, ndof, cdim);
70: }
71: }
72: PetscCall(PetscMalloc1(npoints, &dagPoints));
73: if (eps == 0.0) {
74: for (i = 0, j = 0; i < npoints; i++, j += cdim) {
75: dagPoints[i] = -1;
76: for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
77: for (c = 0; c < cdim; c++) {
78: if (coord[j + c] != allCoords[o + c]) break;
79: }
80: if (c == cdim) {
81: dagPoints[i] = p;
82: break;
83: }
84: }
85: }
86: } else {
87: for (i = 0, j = 0; i < npoints; i++, j += cdim) {
88: PetscReal norm;
90: dagPoints[i] = -1;
91: for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
92: norm = 0.0;
93: for (c = 0; c < cdim; c++) norm += PetscRealPart(PetscSqr(coord[j + c] - allCoords[o + c]));
94: norm = PetscSqrtReal(norm);
95: if (norm <= eps) {
96: dagPoints[i] = p;
97: break;
98: }
99: }
100: }
101: }
102: PetscCall(VecRestoreArrayRead(allCoordsVec, &allCoords));
103: PetscCall(VecRestoreArrayRead(coordinates, &coord));
104: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, dagPoints, PETSC_OWN_POINTER, points));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: #if 0
109: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
110: {
111: const PetscReal p0_x = segmentA[0 * 2 + 0];
112: const PetscReal p0_y = segmentA[0 * 2 + 1];
113: const PetscReal p1_x = segmentA[1 * 2 + 0];
114: const PetscReal p1_y = segmentA[1 * 2 + 1];
115: const PetscReal p2_x = segmentB[0 * 2 + 0];
116: const PetscReal p2_y = segmentB[0 * 2 + 1];
117: const PetscReal p3_x = segmentB[1 * 2 + 0];
118: const PetscReal p3_y = segmentB[1 * 2 + 1];
119: const PetscReal s1_x = p1_x - p0_x;
120: const PetscReal s1_y = p1_y - p0_y;
121: const PetscReal s2_x = p3_x - p2_x;
122: const PetscReal s2_y = p3_y - p2_y;
123: const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);
125: PetscFunctionBegin;
126: *hasIntersection = PETSC_FALSE;
127: /* Non-parallel lines */
128: if (denom != 0.0) {
129: const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
130: const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;
132: if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
133: *hasIntersection = PETSC_TRUE;
134: if (intersection) {
135: intersection[0] = p0_x + (t * s1_x);
136: intersection[1] = p0_y + (t * s1_y);
137: }
138: }
139: }
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
144: static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection)
145: {
146: const PetscReal p0_x = segmentA[0 * 3 + 0];
147: const PetscReal p0_y = segmentA[0 * 3 + 1];
148: const PetscReal p0_z = segmentA[0 * 3 + 2];
149: const PetscReal p1_x = segmentA[1 * 3 + 0];
150: const PetscReal p1_y = segmentA[1 * 3 + 1];
151: const PetscReal p1_z = segmentA[1 * 3 + 2];
152: const PetscReal q0_x = segmentB[0 * 3 + 0];
153: const PetscReal q0_y = segmentB[0 * 3 + 1];
154: const PetscReal q0_z = segmentB[0 * 3 + 2];
155: const PetscReal q1_x = segmentB[1 * 3 + 0];
156: const PetscReal q1_y = segmentB[1 * 3 + 1];
157: const PetscReal q1_z = segmentB[1 * 3 + 2];
158: const PetscReal r0_x = segmentC[0 * 3 + 0];
159: const PetscReal r0_y = segmentC[0 * 3 + 1];
160: const PetscReal r0_z = segmentC[0 * 3 + 2];
161: const PetscReal r1_x = segmentC[1 * 3 + 0];
162: const PetscReal r1_y = segmentC[1 * 3 + 1];
163: const PetscReal r1_z = segmentC[1 * 3 + 2];
164: const PetscReal s0_x = p1_x - p0_x;
165: const PetscReal s0_y = p1_y - p0_y;
166: const PetscReal s0_z = p1_z - p0_z;
167: const PetscReal s1_x = q1_x - q0_x;
168: const PetscReal s1_y = q1_y - q0_y;
169: const PetscReal s1_z = q1_z - q0_z;
170: const PetscReal s2_x = r1_x - r0_x;
171: const PetscReal s2_y = r1_y - r0_y;
172: const PetscReal s2_z = r1_z - r0_z;
173: const PetscReal s3_x = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */
174: const PetscReal s3_y = s1_z * s2_x - s1_x * s2_z;
175: const PetscReal s3_z = s1_x * s2_y - s1_y * s2_x;
176: const PetscReal s4_x = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */
177: const PetscReal s4_y = s0_z * s2_x - s0_x * s2_z;
178: const PetscReal s4_z = s0_x * s2_y - s0_y * s2_x;
179: const PetscReal s5_x = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */
180: const PetscReal s5_y = s1_z * s0_x - s1_x * s0_z;
181: const PetscReal s5_z = s1_x * s0_y - s1_y * s0_x;
182: const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */
184: PetscFunctionBegin;
185: *hasIntersection = PETSC_FALSE;
186: /* Line not parallel to plane */
187: if (denom != 0.0) {
188: const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
189: const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
190: const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;
192: if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) {
193: *hasIntersection = PETSC_TRUE;
194: if (intersection) {
195: intersection[0] = p0_x + (t * s0_x);
196: intersection[1] = p0_y + (t * s0_y);
197: intersection[2] = p0_z + (t * s0_z);
198: }
199: }
200: }
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
203: #endif
205: static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Coords_Internal(DM dm, PetscInt dim, PetscInt cdim, const PetscScalar coords[], const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
206: {
207: PetscReal d[4]; // distance of vertices to the plane
208: PetscReal dp; // distance from origin to the plane
209: PetscInt n = 0;
211: PetscFunctionBegin;
212: if (pos) *pos = PETSC_FALSE;
213: if (Nint) *Nint = 0;
214: if (PetscDefined(USE_DEBUG)) {
215: PetscReal mag = DMPlex_NormD_Internal(cdim, normal);
216: PetscCheck(PetscAbsReal(mag - (PetscReal)1.0) < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Normal vector is not normalized: %g", (double)mag);
217: }
219: dp = DMPlex_DotRealD_Internal(cdim, normal, p);
220: for (PetscInt v = 0; v < dim + 1; ++v) {
221: // d[v] is positive, zero, or negative if vertex i is above, on, or below the plane
222: #if PetscDefined(USE_COMPLEX)
223: PetscReal c[4];
224: for (PetscInt i = 0; i < cdim; ++i) c[i] = PetscRealPart(coords[v * cdim + i]);
225: d[v] = DMPlex_DotRealD_Internal(cdim, normal, c);
226: #else
227: d[v] = DMPlex_DotRealD_Internal(cdim, normal, &coords[v * cdim]);
228: #endif
229: d[v] -= dp;
230: }
232: // If all d are positive or negative, no intersection
233: {
234: PetscInt v;
235: for (v = 0; v < dim + 1; ++v)
236: if (d[v] >= 0.) break;
237: if (v == dim + 1) PetscFunctionReturn(PETSC_SUCCESS);
238: for (v = 0; v < dim + 1; ++v)
239: if (d[v] <= 0.) break;
240: if (v == dim + 1) {
241: if (pos) *pos = PETSC_TRUE;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
244: }
246: for (PetscInt v = 0; v < dim + 1; ++v) {
247: // Points with zero distance are automatically added to the list.
248: if (PetscAbsReal(d[v]) < PETSC_MACHINE_EPSILON) {
249: for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = PetscRealPart(coords[v * cdim + i]);
250: ++n;
251: } else {
252: // For each point with nonzero distance, seek another point with opposite sign
253: // and higher index, and compute the intersection of the line between those
254: // points and the plane.
255: for (PetscInt w = v + 1; w < dim + 1; ++w) {
256: if (d[v] * d[w] < 0.) {
257: PetscReal inv_dist = 1. / (d[v] - d[w]);
258: for (PetscInt i = 0; i < cdim; ++i) intPoints[n * cdim + i] = (d[v] * PetscRealPart(coords[w * cdim + i]) - d[w] * PetscRealPart(coords[v * cdim + i])) * inv_dist;
259: ++n;
260: }
261: }
262: }
263: }
264: // TODO order output points if there are 4
265: *Nint = n;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode DMPlexGetPlaneSimplexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
270: {
271: const PetscScalar *array;
272: PetscScalar *coords = NULL;
273: PetscInt numCoords;
274: PetscBool isDG;
275: PetscInt cdim;
277: PetscFunctionBegin;
278: PetscCall(DMGetCoordinateDim(dm, &cdim));
279: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
280: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
281: PetscCheck(numCoords == dim * (dim + 1), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tetrahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * (dim + 1), numCoords);
282: PetscCall(PetscArrayzero(intPoints, dim * (dim + 1)));
284: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, coords, p, normal, pos, Nint, intPoints));
286: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: static PetscErrorCode DMPlexGetPlaneQuadIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
291: {
292: const PetscScalar *array;
293: PetscScalar *coords = NULL;
294: PetscInt numCoords;
295: PetscBool isDG;
296: PetscInt cdim;
297: PetscScalar tcoords[6] = {0., 0., 0., 0., 0., 0.};
298: const PetscInt vertsA[3] = {0, 1, 3};
299: const PetscInt vertsB[3] = {1, 2, 3};
300: PetscInt NintA, NintB;
302: PetscFunctionBegin;
303: PetscCall(DMGetCoordinateDim(dm, &cdim));
304: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
305: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
306: PetscCheck(numCoords == dim * 4, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 4, numCoords);
307: PetscCall(PetscArrayzero(intPoints, dim * 4));
309: for (PetscInt v = 0; v < 3; ++v)
310: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d];
311: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, intPoints));
312: for (PetscInt v = 0; v < 3; ++v)
313: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d];
314: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[NintA * cdim]));
315: *Nint = NintA + NintB;
317: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: static PetscErrorCode DMPlexGetPlaneHexIntersection_Internal(DM dm, PetscInt dim, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
322: {
323: const PetscScalar *array;
324: PetscScalar *coords = NULL;
325: PetscInt numCoords;
326: PetscBool isDG;
327: PetscInt cdim;
328: PetscScalar tcoords[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
329: // We split using the (2, 4) main diagonal, so all tets contain those vertices
330: const PetscInt vertsA[4] = {0, 1, 2, 4};
331: const PetscInt vertsB[4] = {0, 2, 3, 4};
332: const PetscInt vertsC[4] = {1, 7, 2, 4};
333: const PetscInt vertsD[4] = {2, 7, 6, 4};
334: const PetscInt vertsE[4] = {3, 5, 4, 2};
335: const PetscInt vertsF[4] = {4, 5, 6, 2};
336: PetscInt NintA, NintB, NintC, NintD, NintE, NintF, Nsum = 0;
338: PetscFunctionBegin;
339: PetscCall(DMGetCoordinateDim(dm, &cdim));
340: PetscCheck(cdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM has coordinates in %" PetscInt_FMT "D instead of %" PetscInt_FMT "D", cdim, dim);
341: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
342: PetscCheck(numCoords == dim * 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Hexahedron should have %" PetscInt_FMT " coordinates, not %" PetscInt_FMT, dim * 8, numCoords);
343: PetscCall(PetscArrayzero(intPoints, dim * 18));
345: for (PetscInt v = 0; v < 4; ++v)
346: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsA[v] * cdim + d];
347: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintA, &intPoints[Nsum * cdim]));
348: Nsum += NintA;
349: for (PetscInt v = 0; v < 4; ++v)
350: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsB[v] * cdim + d];
351: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintB, &intPoints[Nsum * cdim]));
352: Nsum += NintB;
353: for (PetscInt v = 0; v < 4; ++v)
354: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsC[v] * cdim + d];
355: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintC, &intPoints[Nsum * cdim]));
356: Nsum += NintC;
357: for (PetscInt v = 0; v < 4; ++v)
358: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsD[v] * cdim + d];
359: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintD, &intPoints[Nsum * cdim]));
360: Nsum += NintD;
361: for (PetscInt v = 0; v < 4; ++v)
362: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsE[v] * cdim + d];
363: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintE, &intPoints[Nsum * cdim]));
364: Nsum += NintE;
365: for (PetscInt v = 0; v < 4; ++v)
366: for (PetscInt d = 0; d < cdim; ++d) tcoords[v * cdim + d] = coords[vertsF[v] * cdim + d];
367: PetscCall(DMPlexGetPlaneSimplexIntersection_Coords_Internal(dm, dim, cdim, tcoords, p, normal, pos, &NintF, &intPoints[Nsum * cdim]));
368: Nsum += NintF;
369: *Nint = Nsum;
371: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /*
376: DMPlexGetPlaneCellIntersection_Internal - Finds the intersection of a plane with a cell
378: Not collective
380: Input Parameters:
381: + dm - the DM
382: . c - the mesh point
383: . p - a point on the plane.
384: - normal - a normal vector to the plane, must be normalized
386: Output Parameters:
387: . pos - `PETSC_TRUE` is the cell is on the positive side of the plane, `PETSC_FALSE` on the negative side
388: + Nint - the number of intersection points, in [0, 4]
389: - intPoints - the coordinates of the intersection points, should be length at least 12
391: Note: The `pos` argument is only meaningful if the number of intersections is 0. The algorithmic idea comes from https://github.com/chrisk314/tet-plane-intersection.
393: Level: developer
395: .seealso:
396: @*/
397: static PetscErrorCode DMPlexGetPlaneCellIntersection_Internal(DM dm, PetscInt c, const PetscReal p[], const PetscReal normal[], PetscBool *pos, PetscInt *Nint, PetscReal intPoints[])
398: {
399: DMPolytopeType ct;
401: PetscFunctionBegin;
402: PetscCall(DMPlexGetCellType(dm, c, &ct));
403: switch (ct) {
404: case DM_POLYTOPE_SEGMENT:
405: case DM_POLYTOPE_TRIANGLE:
406: case DM_POLYTOPE_TETRAHEDRON:
407: PetscCall(DMPlexGetPlaneSimplexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
408: break;
409: case DM_POLYTOPE_QUADRILATERAL:
410: PetscCall(DMPlexGetPlaneQuadIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
411: break;
412: case DM_POLYTOPE_HEXAHEDRON:
413: PetscCall(DMPlexGetPlaneHexIntersection_Internal(dm, DMPolytopeTypeGetDim(ct), c, p, normal, pos, Nint, intPoints));
414: break;
415: default:
416: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No plane intersection for cell %" PetscInt_FMT " with type %s", c, DMPolytopeTypes[ct]);
417: }
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
422: {
423: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
424: const PetscReal x = PetscRealPart(point[0]);
425: PetscReal v0, J, invJ, detJ;
426: PetscReal xi;
428: PetscFunctionBegin;
429: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
430: xi = invJ * (x - v0);
432: if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c;
433: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
438: {
439: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
440: PetscReal xi[2] = {0., 0.};
441: PetscReal x[3], v0[3], J[9], invJ[9], detJ;
442: PetscInt embedDim;
444: PetscFunctionBegin;
445: PetscCall(DMGetCoordinateDim(dm, &embedDim));
446: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
447: for (PetscInt j = 0; j < embedDim; ++j) x[j] = PetscRealPart(point[j]);
448: for (PetscInt i = 0; i < 2; ++i) {
449: for (PetscInt j = 0; j < embedDim; ++j) xi[i] += invJ[i * embedDim + j] * (x[j] - v0[j]);
450: }
451: if ((xi[0] >= -eps) && (xi[1] >= -eps) && (xi[0] + xi[1] <= 2.0 + eps)) *cell = c;
452: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
457: {
458: const PetscInt embedDim = 2;
459: PetscReal x = PetscRealPart(point[0]);
460: PetscReal y = PetscRealPart(point[1]);
461: PetscReal v0[2], J[4], invJ[4], detJ;
462: PetscReal xi, eta, r;
464: PetscFunctionBegin;
465: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
466: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
467: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);
469: xi = PetscMax(xi, 0.0);
470: eta = PetscMax(eta, 0.0);
471: if (xi + eta > 2.0) {
472: r = (xi + eta) / 2.0;
473: xi /= r;
474: eta /= r;
475: }
476: cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0];
477: cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1];
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: // This is the ray-casting, or even-odd algorithm: https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule
482: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Linear_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
483: {
484: const PetscScalar *array;
485: PetscScalar *coords = NULL;
486: const PetscInt faces[8] = {0, 1, 1, 2, 2, 3, 3, 0};
487: PetscReal x = PetscRealPart(point[0]);
488: PetscReal y = PetscRealPart(point[1]);
489: PetscInt crossings = 0, numCoords, embedDim;
490: PetscBool isDG;
492: PetscFunctionBegin;
493: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
494: embedDim = numCoords / 4;
495: PetscCheck(!(numCoords % 4), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
496: // Treat linear quads as Monge surfaces, so we just locate on the projection to x-y (could instead project to 2D)
497: for (PetscInt f = 0; f < 4; ++f) {
498: PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * embedDim + 0]);
499: PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * embedDim + 1]);
500: PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * embedDim + 0]);
501: PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * embedDim + 1]);
503: if ((x == x_j) && (y == y_j)) {
504: // point is a corner
505: crossings = 1;
506: break;
507: }
508: if ((y_j > y) != (y_i > y)) {
509: PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j);
510: if (slope == 0) {
511: // point is a corner
512: crossings = 1;
513: break;
514: }
515: if ((slope < 0) != (y_i < y_j)) ++crossings;
516: }
517: }
518: if (crossings % 2) *cell = c;
519: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
520: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
525: {
526: DM cdm;
527: PetscInt degree, dimR, dimC;
528: PetscFE fe;
529: PetscClassId id;
530: PetscSpace sp;
531: PetscReal pointR[3], ref[3], error;
532: Vec coords;
533: PetscBool found = PETSC_FALSE;
535: PetscFunctionBegin;
536: PetscCall(DMGetDimension(dm, &dimR));
537: PetscCall(DMGetCoordinateDM(dm, &cdm));
538: PetscCall(DMGetDimension(cdm, &dimC));
539: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
540: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
541: if (id != PETSCFE_CLASSID) degree = 1;
542: else {
543: PetscCall(PetscFEGetBasisSpace(fe, &sp));
544: PetscCall(PetscSpaceGetDegree(sp, °ree, NULL));
545: }
546: if (degree == 1) {
547: /* Use simple location method for linear elements*/
548: PetscCall(DMPlexLocatePoint_Quad_2D_Linear_Internal(dm, point, c, cell));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
551: /* Otherwise, we have to solve for the real to reference coordinates */
552: PetscCall(DMGetCoordinatesLocal(dm, &coords));
553: error = PETSC_SQRT_MACHINE_EPSILON;
554: for (PetscInt d = 0; d < dimC; d++) pointR[d] = PetscRealPart(point[d]);
555: PetscCall(DMPlexCoordinatesToReference_FE(cdm, fe, c, 1, pointR, ref, coords, dimC, dimR, 10, &error));
556: if (error < PETSC_SQRT_MACHINE_EPSILON) found = PETSC_TRUE;
557: if ((ref[0] > 1.0 + PETSC_SMALL) || (ref[0] < -1.0 - PETSC_SMALL) || (ref[1] > 1.0 + PETSC_SMALL) || (ref[1] < -1.0 - PETSC_SMALL)) found = PETSC_FALSE;
558: if (PetscDefined(USE_DEBUG) && found) {
559: PetscReal real[3], inverseError = 0, normPoint = DMPlex_NormD_Internal(dimC, pointR);
561: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
562: PetscCall(DMPlexReferenceToCoordinates_FE(cdm, fe, c, 1, ref, real, coords, dimC, dimR));
563: inverseError = DMPlex_DistRealD_Internal(dimC, real, pointR);
564: if (inverseError > PETSC_SQRT_MACHINE_EPSILON * normPoint) found = PETSC_FALSE;
565: if (!found) PetscCall(PetscInfo(dm, "Point (%g, %g, %g) != Mapped Ref Coords (%g, %g, %g) with error %g\n", (double)pointR[0], (double)pointR[1], (double)pointR[2], (double)real[0], (double)real[1], (double)real[2], (double)inverseError));
566: }
567: if (found) *cell = c;
568: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
573: {
574: const PetscInt embedDim = 3;
575: const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
576: PetscReal v0[3], J[9], invJ[9], detJ;
577: PetscReal x = PetscRealPart(point[0]);
578: PetscReal y = PetscRealPart(point[1]);
579: PetscReal z = PetscRealPart(point[2]);
580: PetscReal xi, eta, zeta;
582: PetscFunctionBegin;
583: PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
584: xi = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
585: eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
586: zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);
588: if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
589: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: static PetscErrorCode DMPlexLocatePoint_Hex_3D_Linear_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
594: {
595: const PetscScalar *array;
596: PetscScalar *coords = NULL;
597: const PetscInt faces[24] = {0, 3, 2, 1, 5, 4, 7, 6, 3, 0, 4, 5, 1, 2, 6, 7, 3, 5, 6, 2, 0, 1, 7, 4};
598: PetscBool found = PETSC_TRUE;
599: PetscInt numCoords;
600: PetscBool isDG;
602: PetscFunctionBegin;
603: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
604: PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
605: for (PetscInt f = 0; f < 6; ++f) {
606: /* Check the point is under plane */
607: /* Get face normal */
608: PetscReal v_i[3];
609: PetscReal v_j[3];
610: PetscReal normal[3];
611: PetscReal pp[3];
612: PetscReal dot;
614: v_i[0] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
615: v_i[1] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
616: v_i[2] = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
617: v_j[0] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
618: v_j[1] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
619: v_j[2] = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
620: normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
621: normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
622: normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
623: pp[0] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
624: pp[1] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
625: pp[2] = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
626: dot = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];
628: /* Check that projected point is in face (2D location problem) */
629: if (dot < 0.0) {
630: found = PETSC_FALSE;
631: break;
632: }
633: }
634: if (found) *cell = c;
635: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
636: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: static PetscErrorCode DMPlexLocatePoint_Hex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
641: {
642: DM cdm;
643: PetscInt degree, dimR, dimC;
644: PetscFE fe;
645: PetscClassId id;
646: PetscSpace sp;
647: PetscReal pointR[3], ref[3], error;
648: Vec coords;
649: PetscBool found = PETSC_FALSE;
651: PetscFunctionBegin;
652: PetscCall(DMGetDimension(dm, &dimR));
653: PetscCall(DMGetCoordinateDM(dm, &cdm));
654: PetscCall(DMGetDimension(cdm, &dimC));
655: PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
656: PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
657: if (id != PETSCFE_CLASSID) degree = 1;
658: else {
659: PetscCall(PetscFEGetBasisSpace(fe, &sp));
660: PetscCall(PetscSpaceGetDegree(sp, °ree, NULL));
661: }
662: if (degree == 1) {
663: /* Use simple location method for linear elements*/
664: PetscCall(DMPlexLocatePoint_Hex_3D_Linear_Internal(dm, point, c, cell));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
667: /* Otherwise, we have to solve for the real to reference coordinates */
668: PetscCall(DMGetCoordinatesLocal(dm, &coords));
669: error = PETSC_SQRT_MACHINE_EPSILON;
670: for (PetscInt d = 0; d < dimC; d++) pointR[d] = PetscRealPart(point[d]);
671: PetscCall(DMPlexCoordinatesToReference_FE(cdm, fe, c, 1, pointR, ref, coords, dimC, dimR, 10, &error));
672: if (error < PETSC_SQRT_MACHINE_EPSILON) found = PETSC_TRUE;
673: if ((ref[0] > 1.0 + PETSC_SMALL) || (ref[0] < -1.0 - PETSC_SMALL) || (ref[1] > 1.0 + PETSC_SMALL) || (ref[1] < -1.0 - PETSC_SMALL) || (ref[2] > 1.0 + PETSC_SMALL) || (ref[2] < -1.0 - PETSC_SMALL)) found = PETSC_FALSE;
674: if (PetscDefined(USE_DEBUG) && found) {
675: PetscReal real[3], inverseError = 0, normPoint = DMPlex_NormD_Internal(dimC, pointR);
677: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
678: PetscCall(DMPlexReferenceToCoordinates_FE(cdm, fe, c, 1, ref, real, coords, dimC, dimR));
679: inverseError = DMPlex_DistRealD_Internal(dimC, real, pointR);
680: if (inverseError > PETSC_SQRT_MACHINE_EPSILON * normPoint) found = PETSC_FALSE;
681: if (!found) PetscCall(PetscInfo(dm, "Point (%g, %g, %g) != Mapped Ref Coords (%g, %g, %g) with error %g\n", (double)pointR[0], (double)pointR[1], (double)pointR[2], (double)real[0], (double)real[1], (double)real[2], (double)inverseError));
682: }
683: if (found) *cell = c;
684: else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
689: {
690: PetscInt d;
692: PetscFunctionBegin;
693: box->dim = dim;
694: for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.;
695: PetscFunctionReturn(PETSC_SUCCESS);
696: }
698: /*@
699: PetscGridHashCreate - Create a `PetscGridHash` for spatially locating points in a mesh.
701: Collective
703: Input Parameters:
704: + comm - the MPI communicator
705: . dim - the spatial dimension
706: - point - an initial point used to seed the bounding box, or `NULL` for a zero-initialized box
708: Output Parameter:
709: . box - the newly created `PetscGridHash`
711: Level: developer
713: .seealso: `DMPLEX`, `PetscGridHash`, `PetscGridHashEnlarge()`, `PetscGridHashDestroy()`
714: @*/
715: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
716: {
717: PetscFunctionBegin;
718: PetscCall(PetscCalloc1(1, box));
719: PetscCall(PetscGridHashInitialize_Internal(*box, dim, point));
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: /*@
724: PetscGridHashEnlarge - Enlarge the bounding box of a `PetscGridHash` to include a new point.
726: Not Collective
728: Input Parameters:
729: + box - the `PetscGridHash`
730: - point - the point whose coordinates extend the box's lower and upper bounds
732: Level: developer
734: .seealso: `DMPLEX`, `PetscGridHash`, `PetscGridHashCreate()`, `PetscGridHashDestroy()`
735: @*/
736: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
737: {
738: PetscInt d;
740: PetscFunctionBegin;
741: for (d = 0; d < box->dim; ++d) {
742: box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
743: box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
744: }
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: static PetscErrorCode DMPlexCreateGridHash(DM dm, PetscGridHash *box)
749: {
750: Vec coordinates;
751: const PetscScalar *a;
752: PetscInt cdim, cStart, cEnd;
754: PetscFunctionBegin;
755: PetscCall(DMGetCoordinateDim(dm, &cdim));
756: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
757: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
759: PetscCall(VecGetArrayRead(coordinates, &a));
760: PetscCall(PetscGridHashCreate(PetscObjectComm((PetscObject)dm), cdim, a, box));
761: PetscCall(VecRestoreArrayRead(coordinates, &a));
762: for (PetscInt c = cStart; c < cEnd; ++c) {
763: const PetscScalar *array;
764: PetscScalar *coords = NULL;
765: PetscInt numCoords;
766: PetscBool isDG;
768: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
769: for (PetscInt i = 0; i < numCoords / cdim; ++i) PetscCall(PetscGridHashEnlarge(*box, &coords[i * cdim]));
770: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: /*@C
776: PetscGridHashSetGrid - Divide the grid into boxes
778: Not Collective
780: Input Parameters:
781: + box - The grid hash object
782: . n - The number of boxes in each dimension, may use `PETSC_DETERMINE` for the entries
783: - h - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE`, if not needed you can pass in `NULL`
785: Level: developer
787: .seealso: `DMPLEX`, `PetscGridHashCreate()`
788: @*/
789: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
790: {
791: PetscInt d;
793: PetscFunctionBegin;
794: PetscAssertPointer(n, 2);
795: if (h) PetscAssertPointer(h, 3);
796: for (d = 0; d < box->dim; ++d) {
797: box->extent[d] = box->upper[d] - box->lower[d];
798: if (n[d] == PETSC_DETERMINE) {
799: PetscCheck(h, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Missing h");
800: box->h[d] = h[d];
801: box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
802: } else {
803: box->n[d] = n[d];
804: box->h[d] = box->extent[d] / n[d];
805: }
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: /*@C
811: PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point
813: Not Collective
815: Input Parameters:
816: + box - The grid hash object
817: . numPoints - The number of input points
818: - points - The input point coordinates
820: Output Parameters:
821: + dboxes - An array of `numPoints` x `dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
822: - boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
824: Level: developer
826: Note:
827: This only guarantees that a box contains a point, not that a cell does.
829: .seealso: `DMPLEX`, `PetscGridHashCreate()`
830: @*/
831: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
832: {
833: const PetscReal *lower = box->lower;
834: const PetscReal *upper = box->upper;
835: const PetscReal *h = box->h;
836: const PetscInt *n = box->n;
837: const PetscInt dim = box->dim;
838: PetscInt d, p;
840: PetscFunctionBegin;
841: for (p = 0; p < numPoints; ++p) {
842: for (d = 0; d < dim; ++d) {
843: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
845: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
846: if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
847: PetscCheck(dbox >= 0 && dbox < n[d], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %" PetscInt_FMT " (%g, %g, %g) is outside of our bounding box (%g, %g, %g) - (%g, %g, %g)", p, (double)PetscRealPart(points[p * dim + 0]), dim > 1 ? (double)PetscRealPart(points[p * dim + 1]) : 0.0, dim > 2 ? (double)PetscRealPart(points[p * dim + 2]) : 0.0, (double)lower[0], (double)lower[1], (double)lower[2], (double)upper[0], (double)upper[1], (double)upper[2]);
848: dboxes[p * dim + d] = dbox;
849: }
850: if (boxes)
851: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
852: }
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: /*
857: PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point
859: Not Collective
861: Input Parameters:
862: + box - The grid hash object
863: . cellSection - The PetscSection mapping cells to boxes
864: . numPoints - The number of input points
865: - points - The input point coordinates
867: Output Parameters:
868: + dboxes - An array of `numPoints`*`dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
869: . boxes - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
870: - found - Flag indicating if point was located within a box
872: Level: developer
874: Note:
875: This does an additional check that a cell actually contains the point, and found is `PETSC_FALSE` if no cell does. Thus, this function requires that `cellSection` is already constructed.
877: .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()`
878: */
879: static PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
880: {
881: const PetscReal *lower = box->lower;
882: const PetscReal *upper = box->upper;
883: const PetscReal *h = box->h;
884: const PetscInt *n = box->n;
885: const PetscInt dim = box->dim;
886: PetscInt bStart, bEnd, d, p;
888: PetscFunctionBegin;
890: *found = PETSC_FALSE;
891: PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd));
892: for (p = 0; p < numPoints; ++p) {
893: for (d = 0; d < dim; ++d) {
894: PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);
896: if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
897: if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS);
898: dboxes[p * dim + d] = dbox;
899: }
900: if (boxes)
901: for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
902: // It is possible for a box to overlap no grid cells
903: if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS);
904: }
905: *found = PETSC_TRUE;
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }
909: /*@
910: PetscGridHashDestroy - Destroy a `PetscGridHash` and free its resources.
912: Collective
914: Input Parameter:
915: . box - the `PetscGridHash` to destroy; set to `NULL` on return
917: Level: developer
919: .seealso: `DMPLEX`, `PetscGridHash`, `PetscGridHashCreate()`, `PetscGridHashEnlarge()`
920: @*/
921: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
922: {
923: PetscFunctionBegin;
924: if (*box) {
925: PetscCall(PetscSectionDestroy(&(*box)->cellSection));
926: PetscCall(ISDestroy(&(*box)->cells));
927: PetscCall(DMLabelDestroy(&(*box)->cellsSparse));
928: }
929: PetscCall(PetscFree(*box));
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
934: {
935: DMPolytopeType ct;
937: PetscFunctionBegin;
938: PetscCall(DMPlexGetCellType(dm, cellStart, &ct));
939: switch (ct) {
940: case DM_POLYTOPE_SEGMENT:
941: PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell));
942: break;
943: case DM_POLYTOPE_TRIANGLE:
944: PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell));
945: break;
946: case DM_POLYTOPE_QUADRILATERAL:
947: PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell));
948: break;
949: case DM_POLYTOPE_TETRAHEDRON:
950: PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell));
951: break;
952: case DM_POLYTOPE_HEXAHEDRON:
953: PetscCall(DMPlexLocatePoint_Hex_3D_Internal(dm, point, cellStart, cell));
954: break;
955: default:
956: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: /*
962: DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
963: */
964: static PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
965: {
966: DMPolytopeType ct;
968: PetscFunctionBegin;
969: PetscCall(DMPlexGetCellType(dm, cell, &ct));
970: switch (ct) {
971: case DM_POLYTOPE_TRIANGLE:
972: PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint));
973: break;
974: #if 0
975: case DM_POLYTOPE_QUADRILATERAL:
976: PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break;
977: case DM_POLYTOPE_TETRAHEDRON:
978: PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break;
979: case DM_POLYTOPE_HEXAHEDRON:
980: PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break;
981: #endif
982: default:
983: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
984: }
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: /*
989: DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX`
991: Collective
993: Input Parameter:
994: . dm - The `DMPLEX`
996: Output Parameter:
997: . localBox - The grid hash object
999: Level: developer
1001: Notes:
1002: How do we determine all boxes intersecting a given cell?
1004: 1) Get convex body enclosing cell. We will use a box called the box-hull.
1006: 2) Get smallest brick of boxes enclosing the box-hull
1008: 3) Each box is composed of 6 planes, 3 lower and 3 upper. We loop over dimensions, and
1009: for each new plane determine whether the cell is on the negative side, positive side, or intersects it.
1011: a) If the cell is on the negative side of the lower planes, it is not in the box
1013: b) If the cell is on the positive side of the upper planes, it is not in the box
1015: c) If there is no intersection, it is in the box
1017: d) If any intersection point is within the box limits, it is in the box
1019: .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
1020: */
1021: static PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
1022: {
1023: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
1024: PetscGridHash lbox;
1025: PetscSF sf;
1026: const PetscInt *leaves;
1027: PetscInt *dboxes, *boxes;
1028: PetscInt cdim, cStart, cEnd, Nl = -1;
1029: PetscBool flg;
1031: PetscFunctionBegin;
1032: PetscCall(DMGetCoordinateDim(dm, &cdim));
1033: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1034: PetscCall(DMPlexCreateGridHash(dm, &lbox));
1035: {
1036: PetscInt n[3], d = 3;
1038: PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &d, &flg));
1039: if (flg) {
1040: for (PetscInt i = d; i < cdim; ++i) n[i] = n[d - 1];
1041: } else {
1042: for (PetscInt i = 0; i < cdim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / cdim) * 0.8));
1043: }
1044: PetscCall(PetscGridHashSetGrid(lbox, n, NULL));
1045: if (debug)
1046: PetscCall(PetscPrintf(PETSC_COMM_SELF, "GridHash:\n (%g, %g, %g) -- (%g, %g, %g)\n n %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n h %g %g %g\n", (double)lbox->lower[0], (double)lbox->lower[1], cdim > 2 ? (double)lbox->lower[2] : 0.,
1047: (double)lbox->upper[0], (double)lbox->upper[1], cdim > 2 ? (double)lbox->upper[2] : 0, n[0], n[1], cdim > 2 ? n[2] : 0, (double)lbox->h[0], (double)lbox->h[1], cdim > 2 ? (double)lbox->h[2] : 0.));
1048: }
1050: PetscCall(DMGetPointSF(dm, &sf));
1051: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
1052: Nl = PetscMax(Nl, 0);
1053: PetscCall(PetscCalloc2(16 * cdim, &dboxes, 16, &boxes));
1055: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse));
1056: PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd));
1057: for (PetscInt c = cStart; c < cEnd; ++c) {
1058: PetscReal intPoints[6 * 6 * 6 * 3];
1059: const PetscScalar *array;
1060: PetscScalar *coords = NULL;
1061: const PetscReal *h = lbox->h;
1062: PetscReal normal[9] = {1., 0., 0., 0., 1., 0., 0., 0., 1.};
1063: PetscReal *lowerIntPoints[3] = {&intPoints[0 * 6 * 6 * 3], &intPoints[1 * 6 * 6 * 3], &intPoints[2 * 6 * 6 * 3]};
1064: PetscReal *upperIntPoints[3] = {&intPoints[3 * 6 * 6 * 3], &intPoints[4 * 6 * 6 * 3], &intPoints[5 * 6 * 6 * 3]};
1065: PetscReal lp[3], up[3], *tmp;
1066: PetscInt numCoords, idx, dlim[6], lowerInt[3], upperInt[3];
1067: PetscBool isDG, lower[3], upper[3];
1069: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
1070: if (idx >= 0) continue;
1071: // Get grid of boxes containing the cell
1072: PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
1073: PetscCall(PetscGridHashGetEnclosingBox(lbox, numCoords / cdim, coords, dboxes, boxes));
1074: PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
1075: for (PetscInt d = 0; d < cdim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
1076: for (PetscInt d = cdim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
1077: for (PetscInt e = 1; e < numCoords / cdim; ++e) {
1078: for (PetscInt d = 0; d < cdim; ++d) {
1079: dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * cdim + d]);
1080: dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * cdim + d]);
1081: }
1082: }
1083: if (debug > 4) {
1084: for (PetscInt d = 0; d < cdim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " direction %" PetscInt_FMT " box limits %" PetscInt_FMT "--%" PetscInt_FMT "\n", c, d, dlim[d * 2 + 0], dlim[d * 2 + 1]));
1085: }
1086: // Initialize with lower planes for first box
1087: for (PetscInt d = 0; d < cdim; ++d) {
1088: lp[d] = lbox->lower[d] + dlim[d * 2 + 0] * h[d];
1089: up[d] = lp[d] + h[d];
1090: }
1091: for (PetscInt d = 0; d < cdim; ++d) {
1092: PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, lp, &normal[d * 3], &lower[d], &lowerInt[d], lowerIntPoints[d]));
1093: if (debug > 4) {
1094: if (!lowerInt[d])
1095: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) does not intersect %s\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lower[d] ? "positive" : "negative"));
1096: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " lower direction %" PetscInt_FMT " (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, d, (double)lp[0], (double)lp[1], cdim > 2 ? (double)lp[2] : 0., lowerInt[d]));
1097: }
1098: }
1099: // Loop over grid
1100: for (PetscInt k = dlim[2 * 2 + 0]; k <= dlim[2 * 2 + 1]; ++k, lp[2] = up[2], up[2] += h[2]) {
1101: if (cdim > 2) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 2], &upper[2], &upperInt[2], upperIntPoints[2]));
1102: if (cdim > 2 && debug > 4) {
1103: if (!upperInt[2]) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[2] ? "positive" : "negative"));
1104: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 2 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[2]));
1105: }
1106: for (PetscInt j = dlim[1 * 2 + 0]; j <= dlim[1 * 2 + 1]; ++j, lp[1] = up[1], up[1] += h[1]) {
1107: if (cdim > 1) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 1], &upper[1], &upperInt[1], upperIntPoints[1]));
1108: if (cdim > 1 && debug > 4) {
1109: if (!upperInt[1])
1110: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[1] ? "positive" : "negative"));
1111: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 1 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[1]));
1112: }
1113: for (PetscInt i = dlim[0 * 2 + 0]; i <= dlim[0 * 2 + 1]; ++i, lp[0] = up[0], up[0] += h[0]) {
1114: const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i;
1115: PetscBool excNeg = PETSC_TRUE;
1116: PetscBool excPos = PETSC_TRUE;
1117: PetscInt NlInt = 0;
1118: PetscInt NuInt = 0;
1120: PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 0], &upper[0], &upperInt[0], upperIntPoints[0]));
1121: if (debug > 4) {
1122: if (!upperInt[0])
1123: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) does not intersect %s\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upper[0] ? "positive" : "negative"));
1124: else PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " upper direction 0 (%g, %g, %g) intersects %" PetscInt_FMT " times\n", c, (double)up[0], (double)up[1], cdim > 2 ? (double)up[2] : 0., upperInt[0]));
1125: }
1126: for (PetscInt d = 0; d < cdim; ++d) {
1127: NlInt += lowerInt[d];
1128: NuInt += upperInt[d];
1129: }
1130: // If there is no intersection...
1131: if (!NlInt && !NuInt) {
1132: // If the cell is on the negative side of the lower planes, it is not in the box
1133: for (PetscInt d = 0; d < cdim; ++d)
1134: if (lower[d]) {
1135: excNeg = PETSC_FALSE;
1136: break;
1137: }
1138: // If the cell is on the positive side of the upper planes, it is not in the box
1139: for (PetscInt d = 0; d < cdim; ++d)
1140: if (!upper[d]) {
1141: excPos = PETSC_FALSE;
1142: break;
1143: }
1144: if (excNeg || excPos) {
1145: if (debug && excNeg) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the negative side of the lower plane\n", c));
1146: if (debug && excPos) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is on the positive side of the upper plane\n", c));
1147: continue;
1148: }
1149: // Otherwise it is in the box
1150: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " is contained in box %" PetscInt_FMT "\n", c, box));
1151: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1152: continue;
1153: }
1154: /*
1155: If any intersection point is within the box limits, it is in the box
1156: We need to have tolerances here since intersection point calculations can introduce errors
1157: Initialize a count to track which planes have intersection outside the box.
1158: if two adjacent planes have intersection points upper and lower all outside the box, look
1159: first at if another plane has intersection points outside the box, if so, it is inside the cell
1160: look next if no intersection points exist on the other planes, and check if the planes are on the
1161: outside of the intersection points but on opposite ends. If so, the box cuts through the cell.
1162: */
1163: PetscInt outsideCount[6] = {0, 0, 0, 0, 0, 0};
1164: for (PetscInt plane = 0; plane < cdim; ++plane) {
1165: for (PetscInt ip = 0; ip < lowerInt[plane]; ++ip) {
1166: PetscInt d;
1168: for (d = 0; d < cdim; ++d) {
1169: if ((lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (lowerIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1170: if (lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) outsideCount[d]++; // The lower point is to the left of this box, and we count it
1171: break;
1172: }
1173: }
1174: if (d == cdim) {
1175: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected lower plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1176: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1177: goto end;
1178: }
1179: }
1180: for (PetscInt ip = 0; ip < upperInt[plane]; ++ip) {
1181: PetscInt d;
1183: for (d = 0; d < cdim; ++d) {
1184: if ((upperIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1185: if (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL)) outsideCount[cdim + d]++; // The upper point is to the right of this box, and we count it
1186: break;
1187: }
1188: }
1189: if (d == cdim) {
1190: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Cell %" PetscInt_FMT " intersected upper plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1191: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1192: goto end;
1193: }
1194: }
1195: }
1196: /*
1197: Check the planes with intersections
1198: in 2D, check if the square falls in the middle of a cell
1199: ie all four planes have intersection points outside of the box
1200: You do not want to be doing this, because it means your grid hashing is finer than your grid,
1201: but we should still support it I guess
1202: */
1203: if (cdim == 2) {
1204: PetscInt nIntersects = 0;
1205: for (PetscInt d = 0; d < cdim; ++d) nIntersects += (outsideCount[d] + outsideCount[d + cdim]);
1206: // if the count adds up to 8, that means each plane has 2 external intersections and thus it is in the cell
1207: if (nIntersects == 8) {
1208: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1209: goto end;
1210: }
1211: }
1212: /*
1213: In 3 dimensions, if two adjacent planes have at least 3 intersections outside the cell in the appropriate direction,
1214: we then check the 3rd planar dimension. If a plane falls between intersection points, the cell belongs to that box.
1215: If the planes are on opposite sides of the intersection points, the cell belongs to that box and it passes through the cell.
1216: */
1217: if (cdim == 3) {
1218: PetscInt faces[3] = {0, 0, 0}, checkInternalFace = 0;
1219: // Find two adjacent planes with at least 3 intersection points in the upper and lower
1220: // if the third plane has 3 intersection points or more, a pyramid base is formed on that plane and it is in the cell
1221: for (PetscInt d = 0; d < cdim; ++d)
1222: if (outsideCount[d] >= 3 && outsideCount[cdim + d] >= 3) {
1223: faces[d]++;
1224: checkInternalFace++;
1225: }
1226: if (checkInternalFace == 3) {
1227: // All planes have 3 intersection points, add it.
1228: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1229: goto end;
1230: }
1231: // Gross, figure out which adjacent faces have at least 3 points
1232: PetscInt nonIntersectingFace = -1;
1233: if (faces[0] == faces[1]) nonIntersectingFace = 2;
1234: if (faces[0] == faces[2]) nonIntersectingFace = 1;
1235: if (faces[1] == faces[2]) nonIntersectingFace = 0;
1236: if (nonIntersectingFace >= 0) {
1237: for (PetscInt plane = 0; plane < cdim; ++plane) {
1238: if (!lowerInt[nonIntersectingFace] && !upperInt[nonIntersectingFace]) continue;
1239: // If we have 2 adjacent sides with pyramids of intersection outside of them, and there is a point between the end caps at all, it must be between the two non intersecting ends, and the box is inside the cell.
1240: for (PetscInt ip = 0; ip < lowerInt[nonIntersectingFace]; ++ip) {
1241: if (lowerIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || lowerIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1242: }
1243: for (PetscInt ip = 0; ip < upperInt[nonIntersectingFace]; ++ip) {
1244: if (upperIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || upperIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1245: }
1246: goto end;
1247: }
1248: // The points are within the bonds of the non intersecting planes, add it.
1249: setpoint:
1250: PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1251: goto end;
1252: }
1253: }
1254: end:
1255: lower[0] = upper[0];
1256: lowerInt[0] = upperInt[0];
1257: tmp = lowerIntPoints[0];
1258: lowerIntPoints[0] = upperIntPoints[0];
1259: upperIntPoints[0] = tmp;
1260: }
1261: lp[0] = lbox->lower[0] + dlim[0 * 2 + 0] * h[0];
1262: up[0] = lp[0] + h[0];
1263: lower[1] = upper[1];
1264: lowerInt[1] = upperInt[1];
1265: tmp = lowerIntPoints[1];
1266: lowerIntPoints[1] = upperIntPoints[1];
1267: upperIntPoints[1] = tmp;
1268: }
1269: lp[1] = lbox->lower[1] + dlim[1 * 2 + 0] * h[1];
1270: up[1] = lp[1] + h[1];
1271: lower[2] = upper[2];
1272: lowerInt[2] = upperInt[2];
1273: tmp = lowerIntPoints[2];
1274: lowerIntPoints[2] = upperIntPoints[2];
1275: upperIntPoints[2] = tmp;
1276: }
1277: }
1278: PetscCall(PetscFree2(dboxes, boxes));
1280: if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF));
1281: PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells));
1282: PetscCall(DMLabelDestroy(&lbox->cellsSparse));
1283: *localBox = lbox;
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
1288: {
1289: PetscInt debug = ((DM_Plex *)dm->data)->printLocate;
1290: DM_Plex *mesh = (DM_Plex *)dm->data;
1291: PetscBool hash = mesh->useHashLocation, reuse = PETSC_FALSE;
1292: PetscInt bs, numPoints, numFound, *found = NULL;
1293: PetscInt cdim, Nl = 0, cStart, cEnd, numCells;
1294: PetscSF sf;
1295: const PetscInt *leaves;
1296: const PetscInt *boxCells;
1297: PetscSFNode *cells;
1298: PetscScalar *a;
1299: PetscMPIInt result;
1300: PetscLogDouble t0, t1;
1301: PetscReal gmin[3], gmax[3];
1302: PetscInt terminating_query_type[] = {0, 0, 0};
1303: PetscMPIInt rank;
1305: PetscFunctionBegin;
1306: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1307: PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0));
1308: PetscCall(PetscTime(&t0));
1309: PetscCheck(ltype != DM_POINTLOCATION_NEAREST || hash, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it.");
1310: PetscCall(DMGetCoordinateDim(dm, &cdim));
1311: PetscCall(VecGetBlockSize(v, &bs));
1312: PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result));
1313: PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
1314: // We ignore extra coordinates
1315: PetscCheck(bs >= cdim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %" PetscInt_FMT " must be the mesh coordinate dimension %" PetscInt_FMT, bs, cdim);
1316: PetscCall(DMGetCoordinatesLocalSetUp(dm));
1317: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1318: PetscCall(DMGetPointSF(dm, &sf));
1319: if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
1320: Nl = PetscMax(Nl, 0);
1321: PetscCall(VecGetLocalSize(v, &numPoints));
1322: PetscCall(VecGetArray(v, &a));
1323: numPoints /= bs;
1324: {
1325: const PetscSFNode *sf_cells;
1327: PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells));
1328: if (sf_cells) {
1329: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n"));
1330: cells = (PetscSFNode *)sf_cells;
1331: reuse = PETSC_TRUE;
1332: } else {
1333: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n"));
1334: PetscCall(PetscMalloc1(numPoints, &cells));
1335: /* initialize cells if created */
1336: for (PetscInt p = 0; p < numPoints; p++) {
1337: cells[p].rank = 0;
1338: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1339: }
1340: }
1341: }
1342: PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1343: if (hash) {
1344: if (!mesh->lbox) {
1345: PetscCall(PetscInfo(dm, "Initializing grid hashing\n"));
1346: PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox));
1347: }
1348: /* Designate the local box for each point */
1349: /* Send points to correct process */
1350: /* Search cells that lie in each subbox */
1351: /* Should we bin points before doing search? */
1352: PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells));
1353: }
1354: numFound = 0;
1355: for (PetscInt p = 0; p < numPoints; ++p) {
1356: const PetscScalar *point = &a[p * bs];
1357: PetscInt dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
1358: PetscBool point_outside_domain = PETSC_FALSE;
1360: /* check bounding box of domain */
1361: for (PetscInt d = 0; d < cdim; d++) {
1362: if (PetscRealPart(point[d]) < gmin[d]) {
1363: point_outside_domain = PETSC_TRUE;
1364: break;
1365: }
1366: if (PetscRealPart(point[d]) > gmax[d]) {
1367: point_outside_domain = PETSC_TRUE;
1368: break;
1369: }
1370: }
1371: if (point_outside_domain) {
1372: cells[p].rank = 0;
1373: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1374: terminating_query_type[0]++;
1375: continue;
1376: }
1378: /* check initial values in cells[].index - abort early if found */
1379: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1380: PetscInt c = cells[p].index;
1382: cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1383: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, c, &cell));
1384: if (cell >= 0) {
1385: cells[p].rank = 0;
1386: cells[p].index = cell;
1387: numFound++;
1388: }
1389: }
1390: if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1391: terminating_query_type[1]++;
1392: continue;
1393: }
1395: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]Checking point %" PetscInt_FMT " (%.2g, %.2g, %.2g)\n", rank, p, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), cdim > 2 ? (double)PetscRealPart(point[2]) : 0.));
1396: if (hash) {
1397: PetscBool found_box;
1399: /* allow for case that point is outside box - abort early */
1400: PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box));
1401: if (found_box) {
1402: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", rank, bin, dbin[0], dbin[1], cdim > 2 ? dbin[2] : 0));
1403: /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
1404: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1405: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1406: for (PetscInt c = cellOffset; c < cellOffset + numCells; ++c) {
1407: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] Checking for point in cell %" PetscInt_FMT "\n", rank, boxCells[c]));
1408: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, boxCells[c], &cell));
1409: if (cell >= 0) {
1410: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] FOUND in cell %" PetscInt_FMT "\n", rank, cell));
1411: cells[p].rank = 0;
1412: cells[p].index = cell;
1413: numFound++;
1414: terminating_query_type[2]++;
1415: break;
1416: }
1417: }
1418: }
1419: } else {
1420: PetscBool found = PETSC_FALSE;
1421: for (PetscInt c = cStart; c < cEnd; ++c) {
1422: PetscInt idx;
1424: PetscCall(PetscFindInt(c, Nl, leaves, &idx));
1425: if (idx >= 0) continue;
1426: PetscCall(DMPlexLocatePoint_Internal(dm, cdim, point, c, &cell));
1427: if (cell >= 0) {
1428: cells[p].rank = 0;
1429: cells[p].index = cell;
1430: numFound++;
1431: terminating_query_type[2]++;
1432: found = PETSC_TRUE;
1433: break;
1434: }
1435: }
1436: if (!found) terminating_query_type[0]++;
1437: }
1438: }
1439: if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells));
1440: if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
1441: for (PetscInt p = 0; p < numPoints; p++) {
1442: const PetscScalar *point = &a[p * bs];
1443: PetscReal cpoint[3] = {0, 0, 0}, diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
1444: PetscInt dbin[3] = {-1, -1, -1}, bin, cellOffset, bestc = -1;
1446: if (cells[p].index < 0) {
1447: PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin));
1448: PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1449: PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1450: for (PetscInt c = cellOffset; c < cellOffset + numCells; ++c) {
1451: PetscCall(DMPlexClosestPoint_Internal(dm, cdim, point, boxCells[c], cpoint));
1452: for (PetscInt d = 0; d < cdim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
1453: dist = DMPlex_NormD_Internal(cdim, diff);
1454: if (dist < distMax) {
1455: for (PetscInt d = 0; d < cdim; ++d) best[d] = cpoint[d];
1456: bestc = boxCells[c];
1457: distMax = dist;
1458: }
1459: }
1460: if (distMax < PETSC_MAX_REAL) {
1461: ++numFound;
1462: cells[p].rank = 0;
1463: cells[p].index = bestc;
1464: for (PetscInt d = 0; d < cdim; ++d) a[p * bs + d] = best[d];
1465: }
1466: }
1467: }
1468: }
1469: /* This code is only be relevant when interfaced to parallel point location */
1470: /* Check for highest numbered proc that claims a point (do we care?) */
1471: if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
1472: PetscCall(PetscMalloc1(numFound, &found));
1473: numFound = 0;
1474: for (PetscInt p = 0; p < numPoints; p++) {
1475: if (cells[p].rank >= 0 && cells[p].index >= 0) {
1476: if (numFound < p) cells[numFound] = cells[p];
1477: found[numFound++] = p;
1478: }
1479: }
1480: }
1481: PetscCall(VecRestoreArray(v, &a));
1482: if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER));
1483: PetscCall(PetscTime(&t1));
1484: if (hash) {
1485: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [hash]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2]));
1486: } else {
1487: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] terminating_query_type : %" PetscInt_FMT " [outside domain] : %" PetscInt_FMT " [inside initial cell] : %" PetscInt_FMT " [brute-force]\n", terminating_query_type[0], terminating_query_type[1], terminating_query_type[2]));
1488: }
1489: PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] npoints %" PetscInt_FMT " : time(rank0) %1.2e (sec): points/sec %1.4e\n", numPoints, t1 - t0, numPoints / (t1 - t0)));
1490: PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0));
1491: PetscFunctionReturn(PETSC_SUCCESS);
1492: }
1494: /*@
1495: DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates
1497: Not Collective
1499: Input/Output Parameter:
1500: . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x, an array of size 4, last two entries are unchanged
1502: Output Parameter:
1503: . R - The rotation which accomplishes the projection, array of size 4
1505: Level: developer
1507: .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1508: @*/
1509: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
1510: {
1511: const PetscReal x = PetscRealPart(coords[2] - coords[0]);
1512: const PetscReal y = PetscRealPart(coords[3] - coords[1]);
1513: const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;
1515: PetscFunctionBegin;
1516: R[0] = c;
1517: R[1] = -s;
1518: R[2] = s;
1519: R[3] = c;
1520: coords[0] = 0.0;
1521: coords[1] = r;
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: /*@
1526: DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates
1528: Not Collective
1530: Input/Output Parameter:
1531: . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z, an array of size 6, the other entries are unchanged
1533: Output Parameter:
1534: . R - The rotation which accomplishes the projection, an array of size 9
1536: Level: developer
1538: Note:
1539: This uses the basis completion described by Frisvad {cite}`frisvad2012building`
1541: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1542: @*/
1543: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1544: {
1545: PetscReal x = PetscRealPart(coords[3] - coords[0]);
1546: PetscReal y = PetscRealPart(coords[4] - coords[1]);
1547: PetscReal z = PetscRealPart(coords[5] - coords[2]);
1548: PetscReal r = PetscSqrtReal(x * x + y * y + z * z);
1549: PetscReal rinv = 1. / r;
1551: PetscFunctionBegin;
1552: x *= rinv;
1553: y *= rinv;
1554: z *= rinv;
1555: if (x > 0.) {
1556: PetscReal inv1pX = 1. / (1. + x);
1558: R[0] = x;
1559: R[1] = -y;
1560: R[2] = -z;
1561: R[3] = y;
1562: R[4] = 1. - y * y * inv1pX;
1563: R[5] = -y * z * inv1pX;
1564: R[6] = z;
1565: R[7] = -y * z * inv1pX;
1566: R[8] = 1. - z * z * inv1pX;
1567: } else {
1568: PetscReal inv1mX = 1. / (1. - x);
1570: R[0] = x;
1571: R[1] = z;
1572: R[2] = y;
1573: R[3] = y;
1574: R[4] = -y * z * inv1mX;
1575: R[5] = 1. - y * y * inv1mX;
1576: R[6] = z;
1577: R[7] = 1. - z * z * inv1mX;
1578: R[8] = -y * z * inv1mX;
1579: }
1580: coords[0] = 0.0;
1581: coords[1] = r;
1582: coords[2] = 0.0;
1583: PetscFunctionReturn(PETSC_SUCCESS);
1584: }
1586: /*@
1587: DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
1588: plane. The normal is defined by positive orientation of the first 3 points.
1590: Not Collective
1592: Input Parameter:
1593: . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
1595: Input/Output Parameter:
1596: . coords - The interlaced coordinates of each coplanar 3D point; on output the first
1597: 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
1599: Output Parameter:
1600: . R - 3x3 row-major rotation matrix whose columns are the tangent basis [t1, t2, n]. Multiplying by R^T transforms from original frame to tangent frame.
1602: Level: developer
1604: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1605: @*/
1606: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1607: {
1608: PetscReal x1[3], x2[3], n[3], c[3], norm;
1609: const PetscInt dim = 3;
1610: PetscInt d;
1612: PetscFunctionBegin;
1613: /* 0) Calculate normal vector */
1614: for (d = 0; d < dim; ++d) {
1615: x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1616: x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1617: }
1618: // n = x1 \otimes x2
1619: n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1620: n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1621: n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1622: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1623: for (d = 0; d < dim; d++) n[d] /= norm;
1624: norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1625: for (d = 0; d < dim; d++) x1[d] /= norm;
1626: // x2 = n \otimes x1
1627: x2[0] = n[1] * x1[2] - n[2] * x1[1];
1628: x2[1] = n[2] * x1[0] - n[0] * x1[2];
1629: x2[2] = n[0] * x1[1] - n[1] * x1[0];
1630: for (d = 0; d < dim; d++) {
1631: R[d * dim + 0] = x1[d];
1632: R[d * dim + 1] = x2[d];
1633: R[d * dim + 2] = n[d];
1634: c[d] = PetscRealPart(coords[0 * dim + d]);
1635: }
1636: for (PetscInt p = 0; p < coordSize / dim; p++) {
1637: PetscReal y[3];
1638: for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1639: for (d = 0; d < 2; d++) coords[p * 2 + d] = R[0 * dim + d] * y[0] + R[1 * dim + d] * y[1] + R[2 * dim + d] * y[2];
1640: }
1641: PetscFunctionReturn(PETSC_SUCCESS);
1642: }
1644: PETSC_UNUSED static inline void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1645: {
1646: /* Signed volume is 1/2 the determinant
1648: | 1 1 1 |
1649: | x0 x1 x2 |
1650: | y0 y1 y2 |
1652: but if x0,y0 is the origin, we have
1654: | x1 x2 |
1655: | y1 y2 |
1656: */
1657: const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1658: const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1659: PetscReal M[4], detM;
1660: M[0] = x1;
1661: M[1] = x2;
1662: M[2] = y1;
1663: M[3] = y2;
1664: DMPlex_Det2D_Internal(&detM, M);
1665: *vol = 0.5 * detM;
1666: (void)PetscLogFlops(5.0);
1667: }
1669: PETSC_UNUSED static inline void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1670: {
1671: /* Signed volume is 1/6th of the determinant
1673: | 1 1 1 1 |
1674: | x0 x1 x2 x3 |
1675: | y0 y1 y2 y3 |
1676: | z0 z1 z2 z3 |
1678: but if x0,y0,z0 is the origin, we have
1680: | x1 x2 x3 |
1681: | y1 y2 y3 |
1682: | z1 z2 z3 |
1683: */
1684: const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1685: const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1686: const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1687: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1688: PetscReal M[9], detM;
1689: M[0] = x1;
1690: M[1] = x2;
1691: M[2] = x3;
1692: M[3] = y1;
1693: M[4] = y2;
1694: M[5] = y3;
1695: M[6] = z1;
1696: M[7] = z2;
1697: M[8] = z3;
1698: DMPlex_Det3D_Internal(&detM, M);
1699: *vol = -onesixth * detM;
1700: (void)PetscLogFlops(10.0);
1701: }
1703: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1704: {
1705: const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1706: DMPlex_Det3D_Internal(vol, coords);
1707: *vol *= -onesixth;
1708: }
1710: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1711: {
1712: PetscSection coordSection;
1713: Vec coordinates;
1714: const PetscScalar *coords;
1715: PetscInt dim, d, off;
1717: PetscFunctionBegin;
1718: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1719: PetscCall(DMGetCoordinateSection(dm, &coordSection));
1720: PetscCall(PetscSectionGetDof(coordSection, e, &dim));
1721: if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
1722: PetscCall(PetscSectionGetOffset(coordSection, e, &off));
1723: PetscCall(VecGetArrayRead(coordinates, &coords));
1724: if (v0) {
1725: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1726: }
1727: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1728: *detJ = 1.;
1729: if (J) {
1730: for (d = 0; d < dim * dim; d++) J[d] = 0.;
1731: for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1732: if (invJ) {
1733: for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1734: for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1735: }
1736: }
1737: PetscFunctionReturn(PETSC_SUCCESS);
1738: }
1740: /*@C
1741: DMPlexGetCellCoordinates - Get coordinates for a cell, taking into account periodicity
1743: Not Collective
1745: Input Parameters:
1746: + dm - The `DMPLEX`
1747: - cell - The cell number
1749: Output Parameters:
1750: + isDG - Using cellwise coordinates
1751: . Nc - The number of coordinates
1752: . array - The coordinate array
1753: - coords - The cell coordinates
1755: Level: developer
1757: .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1758: @*/
1759: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1760: {
1761: DM cdm;
1762: Vec coordinates;
1763: PetscSection cs;
1764: const PetscScalar *ccoords;
1765: PetscInt pStart, pEnd;
1767: PetscFunctionBeginHot;
1768: *isDG = PETSC_FALSE;
1769: *Nc = 0;
1770: *array = NULL;
1771: *coords = NULL;
1772: /* Check for cellwise coordinates */
1773: PetscCall(DMGetCellCoordinateSection(dm, &cs));
1774: if (!cs) goto cg;
1775: /* Check that the cell exists in the cellwise section */
1776: PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd));
1777: if (cell < pStart || cell >= pEnd) goto cg;
1778: /* Check for cellwise coordinates for this cell */
1779: PetscCall(PetscSectionGetDof(cs, cell, Nc));
1780: if (!*Nc) goto cg;
1781: /* Check for cellwise coordinates */
1782: PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates));
1783: if (!coordinates) goto cg;
1784: /* Get cellwise coordinates */
1785: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1786: PetscCall(VecGetArrayRead(coordinates, array));
1787: PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords));
1788: PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1789: PetscCall(PetscArraycpy(*coords, ccoords, *Nc));
1790: PetscCall(VecRestoreArrayRead(coordinates, array));
1791: *isDG = PETSC_TRUE;
1792: PetscFunctionReturn(PETSC_SUCCESS);
1793: cg:
1794: /* Use continuous coordinates */
1795: PetscCall(DMGetCoordinateDM(dm, &cdm));
1796: PetscCall(DMGetCoordinateSection(dm, &cs));
1797: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1798: PetscCall(DMPlexVecGetOrientedClosure(cdm, cs, PETSC_FALSE, coordinates, cell, 0, Nc, coords));
1799: PetscFunctionReturn(PETSC_SUCCESS);
1800: }
1802: /*@C
1803: DMPlexRestoreCellCoordinates - Get coordinates for a cell, taking into account periodicity
1805: Not Collective
1807: Input Parameters:
1808: + dm - The `DMPLEX`
1809: - cell - The cell number
1811: Output Parameters:
1812: + isDG - Using cellwise coordinates
1813: . Nc - The number of coordinates
1814: . array - The coordinate array
1815: - coords - The cell coordinates
1817: Level: developer
1819: .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1820: @*/
1821: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1822: {
1823: DM cdm;
1824: PetscSection cs;
1825: Vec coordinates;
1827: PetscFunctionBeginHot;
1828: if (*isDG) {
1829: PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1830: PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1831: } else {
1832: PetscCall(DMGetCoordinateDM(dm, &cdm));
1833: PetscCall(DMGetCoordinateSection(dm, &cs));
1834: PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1835: PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, coords));
1836: }
1837: PetscFunctionReturn(PETSC_SUCCESS);
1838: }
1840: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1841: {
1842: const PetscScalar *array;
1843: PetscScalar *coords = NULL;
1844: PetscInt numCoords, d;
1845: PetscBool isDG;
1847: PetscFunctionBegin;
1848: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1849: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1850: *detJ = 0.0;
1851: if (numCoords == 6) {
1852: const PetscInt dim = 3;
1853: PetscReal R[9], J0;
1855: if (v0) {
1856: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1857: }
1858: PetscCall(DMPlexComputeProjection3Dto1D(coords, R));
1859: if (J) {
1860: J0 = 0.5 * PetscRealPart(coords[1]);
1861: J[0] = R[0] * J0;
1862: J[1] = R[1];
1863: J[2] = R[2];
1864: J[3] = R[3] * J0;
1865: J[4] = R[4];
1866: J[5] = R[5];
1867: J[6] = R[6] * J0;
1868: J[7] = R[7];
1869: J[8] = R[8];
1870: DMPlex_Det3D_Internal(detJ, J);
1871: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1872: }
1873: } else if (numCoords == 4) {
1874: const PetscInt dim = 2;
1875: PetscReal R[4], J0;
1877: if (v0) {
1878: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1879: }
1880: PetscCall(DMPlexComputeProjection2Dto1D(coords, R));
1881: if (J) {
1882: J0 = 0.5 * PetscRealPart(coords[1]);
1883: J[0] = R[0] * J0;
1884: J[1] = R[1];
1885: J[2] = R[2] * J0;
1886: J[3] = R[3];
1887: DMPlex_Det2D_Internal(detJ, J);
1888: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1889: }
1890: } else if (numCoords == 2) {
1891: const PetscInt dim = 1;
1893: if (v0) {
1894: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1895: }
1896: if (J) {
1897: J[0] = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1898: *detJ = J[0];
1899: PetscCall(PetscLogFlops(2.0));
1900: if (invJ) {
1901: invJ[0] = 1.0 / J[0];
1902: PetscCall(PetscLogFlops(1.0));
1903: }
1904: }
1905: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for segment %" PetscInt_FMT " is %" PetscInt_FMT " != 2 or 4 or 6", e, numCoords);
1906: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1907: PetscFunctionReturn(PETSC_SUCCESS);
1908: }
1910: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1911: {
1912: const PetscScalar *array;
1913: PetscScalar *coords = NULL;
1914: PetscInt numCoords, d;
1915: PetscBool isDG;
1917: PetscFunctionBegin;
1918: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1919: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1920: *detJ = 0.0;
1921: if (numCoords == 9) {
1922: const PetscInt dim = 3;
1923: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1925: if (v0) {
1926: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1927: }
1928: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1929: if (J) {
1930: const PetscInt pdim = 2;
1932: for (d = 0; d < pdim; d++) {
1933: for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1934: }
1935: PetscCall(PetscLogFlops(8.0));
1936: DMPlex_Det3D_Internal(detJ, J0);
1937: for (d = 0; d < dim; d++) {
1938: for (PetscInt f = 0; f < dim; f++) {
1939: J[d * dim + f] = 0.0;
1940: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1941: }
1942: }
1943: PetscCall(PetscLogFlops(18.0));
1944: }
1945: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1946: } else if (numCoords == 6) {
1947: const PetscInt dim = 2;
1949: if (v0) {
1950: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1951: }
1952: if (J) {
1953: for (d = 0; d < dim; d++) {
1954: for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1955: }
1956: PetscCall(PetscLogFlops(8.0));
1957: DMPlex_Det2D_Internal(detJ, J);
1958: }
1959: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1960: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1961: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1962: PetscFunctionReturn(PETSC_SUCCESS);
1963: }
1965: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1966: {
1967: const PetscScalar *array;
1968: PetscScalar *coords = NULL;
1969: PetscInt numCoords, d;
1970: PetscBool isDG;
1972: PetscFunctionBegin;
1973: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1974: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1975: if (!Nq) {
1976: PetscInt vorder[4] = {0, 1, 2, 3};
1978: if (isTensor) {
1979: vorder[2] = 3;
1980: vorder[3] = 2;
1981: }
1982: *detJ = 0.0;
1983: if (numCoords == 12) {
1984: const PetscInt dim = 3;
1985: PetscReal R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
1987: if (v) {
1988: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1989: }
1990: PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1991: if (J) {
1992: const PetscInt pdim = 2;
1994: for (d = 0; d < pdim; d++) {
1995: J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1996: J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1997: }
1998: PetscCall(PetscLogFlops(8.0));
1999: DMPlex_Det3D_Internal(detJ, J0);
2000: for (d = 0; d < dim; d++) {
2001: for (PetscInt f = 0; f < dim; f++) {
2002: J[d * dim + f] = 0.0;
2003: for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
2004: }
2005: }
2006: PetscCall(PetscLogFlops(18.0));
2007: }
2008: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2009: } else if (numCoords == 8) {
2010: const PetscInt dim = 2;
2012: if (v) {
2013: for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2014: }
2015: if (J) {
2016: for (d = 0; d < dim; d++) {
2017: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
2018: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
2019: }
2020: PetscCall(PetscLogFlops(8.0));
2021: DMPlex_Det2D_Internal(detJ, J);
2022: }
2023: if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
2024: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
2025: } else {
2026: const PetscInt Nv = 4;
2027: const PetscInt dimR = 2;
2028: PetscInt zToPlex[4] = {0, 1, 3, 2};
2029: PetscReal zOrder[12];
2030: PetscReal zCoeff[12];
2031: PetscInt i, j, k, l, dim;
2033: if (isTensor) {
2034: zToPlex[2] = 2;
2035: zToPlex[3] = 3;
2036: }
2037: if (numCoords == 12) {
2038: dim = 3;
2039: } else if (numCoords == 8) {
2040: dim = 2;
2041: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
2042: for (i = 0; i < Nv; i++) {
2043: PetscInt zi = zToPlex[i];
2045: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
2046: }
2047: for (j = 0; j < dim; j++) {
2048: /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
2049: \phi^0 = (1 - xi - eta + xi eta) --> 1 = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
2050: \phi^1 = (1 + xi - eta - xi eta) --> xi = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
2051: \phi^2 = (1 - xi + eta - xi eta) --> eta = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
2052: \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
2053: */
2054: zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2055: zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2056: zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2057: zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
2058: }
2059: for (i = 0; i < Nq; i++) {
2060: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];
2062: if (v) {
2063: PetscReal extPoint[4];
2065: extPoint[0] = 1.;
2066: extPoint[1] = xi;
2067: extPoint[2] = eta;
2068: extPoint[3] = xi * eta;
2069: for (j = 0; j < dim; j++) {
2070: PetscReal val = 0.;
2072: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
2073: v[i * dim + j] = val;
2074: }
2075: }
2076: if (J) {
2077: PetscReal extJ[8];
2079: extJ[0] = 0.;
2080: extJ[1] = 0.;
2081: extJ[2] = 1.;
2082: extJ[3] = 0.;
2083: extJ[4] = 0.;
2084: extJ[5] = 1.;
2085: extJ[6] = eta;
2086: extJ[7] = xi;
2087: for (j = 0; j < dim; j++) {
2088: for (k = 0; k < dimR; k++) {
2089: PetscReal val = 0.;
2091: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
2092: J[i * dim * dim + dim * j + k] = val;
2093: }
2094: }
2095: if (dim == 3) { /* put the cross product in the third component of the Jacobian */
2096: PetscReal x, y, z;
2097: PetscReal *iJ = &J[i * dim * dim];
2098: PetscReal norm;
2100: x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
2101: y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
2102: z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
2103: norm = PetscSqrtReal(x * x + y * y + z * z);
2104: iJ[2] = x / norm;
2105: iJ[5] = y / norm;
2106: iJ[8] = z / norm;
2107: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2108: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2109: } else {
2110: DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
2111: if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2112: }
2113: }
2114: }
2115: }
2116: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }
2120: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2121: {
2122: const PetscScalar *array;
2123: PetscScalar *coords = NULL;
2124: const PetscInt dim = 3;
2125: PetscInt numCoords, d;
2126: PetscBool isDG;
2128: PetscFunctionBegin;
2129: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2130: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2131: *detJ = 0.0;
2132: if (v0) {
2133: for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
2134: }
2135: if (J) {
2136: for (d = 0; d < dim; d++) {
2137: /* I orient with outward face normals */
2138: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2139: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2140: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2141: }
2142: PetscCall(PetscLogFlops(18.0));
2143: DMPlex_Det3D_Internal(detJ, J);
2144: }
2145: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2146: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2147: PetscFunctionReturn(PETSC_SUCCESS);
2148: }
2150: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2151: {
2152: const PetscScalar *array;
2153: PetscScalar *coords = NULL;
2154: const PetscInt dim = 3;
2155: PetscInt numCoords;
2156: PetscBool isDG;
2158: PetscFunctionBegin;
2159: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2160: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2161: if (!Nq) {
2162: *detJ = 0.0;
2163: if (v) {
2164: for (PetscInt d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2165: }
2166: if (J) {
2167: for (PetscInt d = 0; d < dim; d++) {
2168: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2169: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2170: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2171: }
2172: PetscCall(PetscLogFlops(18.0));
2173: DMPlex_Det3D_Internal(detJ, J);
2174: }
2175: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2176: } else {
2177: const PetscInt Nv = 8;
2178: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2179: const PetscInt dim = 3;
2180: const PetscInt dimR = 3;
2181: PetscReal zOrder[24];
2182: PetscReal zCoeff[24];
2183: PetscInt i, j, k, l;
2185: for (i = 0; i < Nv; i++) {
2186: PetscInt zi = zToPlex[i];
2188: for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
2189: }
2190: for (j = 0; j < dim; j++) {
2191: zCoeff[dim * 0 + j] = 0.125 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2192: zCoeff[dim * 1 + j] = 0.125 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2193: zCoeff[dim * 2 + j] = 0.125 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2194: zCoeff[dim * 3 + j] = 0.125 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2195: zCoeff[dim * 4 + j] = 0.125 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] + zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2196: zCoeff[dim * 5 + j] = 0.125 * (+zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] + zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2197: zCoeff[dim * 6 + j] = 0.125 * (+zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] - zOrder[dim * 3 + j] - zOrder[dim * 4 + j] - zOrder[dim * 5 + j] + zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2198: zCoeff[dim * 7 + j] = 0.125 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] - zOrder[dim * 3 + j] + zOrder[dim * 4 + j] - zOrder[dim * 5 + j] - zOrder[dim * 6 + j] + zOrder[dim * 7 + j]);
2199: }
2200: for (i = 0; i < Nq; i++) {
2201: PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];
2203: if (v) {
2204: PetscReal extPoint[8];
2206: extPoint[0] = 1.;
2207: extPoint[1] = xi;
2208: extPoint[2] = eta;
2209: extPoint[3] = xi * eta;
2210: extPoint[4] = theta;
2211: extPoint[5] = theta * xi;
2212: extPoint[6] = theta * eta;
2213: extPoint[7] = theta * eta * xi;
2214: for (j = 0; j < dim; j++) {
2215: PetscReal val = 0.;
2217: for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
2218: v[i * dim + j] = val;
2219: }
2220: }
2221: if (J) {
2222: PetscReal extJ[24];
2224: extJ[0] = 0.;
2225: extJ[1] = 0.;
2226: extJ[2] = 0.;
2227: extJ[3] = 1.;
2228: extJ[4] = 0.;
2229: extJ[5] = 0.;
2230: extJ[6] = 0.;
2231: extJ[7] = 1.;
2232: extJ[8] = 0.;
2233: extJ[9] = eta;
2234: extJ[10] = xi;
2235: extJ[11] = 0.;
2236: extJ[12] = 0.;
2237: extJ[13] = 0.;
2238: extJ[14] = 1.;
2239: extJ[15] = theta;
2240: extJ[16] = 0.;
2241: extJ[17] = xi;
2242: extJ[18] = 0.;
2243: extJ[19] = theta;
2244: extJ[20] = eta;
2245: extJ[21] = theta * eta;
2246: extJ[22] = theta * xi;
2247: extJ[23] = eta * xi;
2249: for (j = 0; j < dim; j++) {
2250: for (k = 0; k < dimR; k++) {
2251: PetscReal val = 0.;
2253: for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
2254: J[i * dim * dim + dim * j + k] = val;
2255: }
2256: }
2257: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2258: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2259: }
2260: }
2261: }
2262: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2263: PetscFunctionReturn(PETSC_SUCCESS);
2264: }
2266: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2267: {
2268: const PetscScalar *array;
2269: PetscScalar *coords = NULL;
2270: const PetscInt dim = 3;
2271: PetscInt numCoords;
2272: PetscBool isDG;
2274: PetscFunctionBegin;
2275: PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2276: PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2277: if (!Nq) {
2278: /* Assume that the map to the reference is affine */
2279: *detJ = 0.0;
2280: if (v) {
2281: for (PetscInt d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2282: }
2283: if (J) {
2284: for (PetscInt d = 0; d < dim; d++) {
2285: J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2286: J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2287: J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2288: }
2289: PetscCall(PetscLogFlops(18.0));
2290: DMPlex_Det3D_Internal(detJ, J);
2291: }
2292: if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2293: } else {
2294: const PetscInt dim = 3;
2295: const PetscInt dimR = 3;
2296: const PetscInt Nv = 6;
2297: PetscReal verts[18];
2298: PetscReal coeff[18];
2299: PetscInt i, j, k, l;
2301: for (i = 0; i < Nv; ++i)
2302: for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
2303: for (j = 0; j < dim; ++j) {
2304: /* Check for triangle,
2305: phi^0 = -1/2 (xi + eta) chi^0 = delta(-1, -1) x(xi) = \sum_k x_k phi^k(xi) = \sum_k chi^k(x) phi^k(xi)
2306: phi^1 = 1/2 (1 + xi) chi^1 = delta( 1, -1) y(xi) = \sum_k y_k phi^k(xi) = \sum_k chi^k(y) phi^k(xi)
2307: phi^2 = 1/2 (1 + eta) chi^2 = delta(-1, 1)
2309: phi^0 + phi^1 + phi^2 = 1 coef_1 = 1/2 ( chi^1 + chi^2)
2310: -phi^0 + phi^1 - phi^2 = xi coef_xi = 1/2 (-chi^0 + chi^1)
2311: -phi^0 - phi^1 + phi^2 = eta coef_eta = 1/2 (-chi^0 + chi^2)
2313: < chi_0 chi_1 chi_2> A / 1 1 1 \ / phi_0 \ <chi> I <phi>^T so we need the inverse transpose
2314: | -1 1 -1 | | phi_1 | =
2315: \ -1 -1 1 / \ phi_2 /
2317: Check phi^0: 1/2 (phi^0 chi^1 + phi^0 chi^2 + phi^0 chi^0 - phi^0 chi^1 + phi^0 chi^0 - phi^0 chi^2) = phi^0 chi^0
2318: */
2319: /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
2320: \phi^0 = 1/4 ( -xi - eta + xi zeta + eta zeta) --> / 1 1 1 1 1 1 \ 1
2321: \phi^1 = 1/4 (1 + eta - zeta - eta zeta) --> | -1 1 -1 -1 -1 1 | eta
2322: \phi^2 = 1/4 (1 + xi - zeta - xi zeta) --> | -1 -1 1 -1 1 -1 | xi
2323: \phi^3 = 1/4 ( -xi - eta - xi zeta - eta zeta) --> | -1 -1 -1 1 1 1 | zeta
2324: \phi^4 = 1/4 (1 + xi + zeta + xi zeta) --> | 1 1 -1 -1 1 -1 | xi zeta
2325: \phi^5 = 1/4 (1 + eta + zeta + eta zeta) --> \ 1 -1 1 -1 -1 1 / eta zeta
2326: 1/4 / 0 1 1 0 1 1 \
2327: | -1 1 0 -1 0 1 |
2328: | -1 0 1 -1 1 0 |
2329: | 0 -1 -1 0 1 1 |
2330: | 1 0 -1 -1 1 0 |
2331: \ 1 -1 0 -1 0 1 /
2332: */
2333: coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2334: coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2335: coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2336: coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2337: coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2338: coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2339: /* For reference prism:
2340: {0, 0, 0}
2341: {0, 1, 0}
2342: {1, 0, 0}
2343: {0, 0, 1}
2344: {0, 0, 0}
2345: {0, 0, 0}
2346: */
2347: }
2348: for (i = 0; i < Nq; ++i) {
2349: const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];
2351: if (v) {
2352: PetscReal extPoint[6];
2354: extPoint[0] = 1.;
2355: extPoint[1] = eta;
2356: extPoint[2] = xi;
2357: extPoint[3] = zeta;
2358: extPoint[4] = xi * zeta;
2359: extPoint[5] = eta * zeta;
2360: for (PetscInt c = 0; c < dim; ++c) {
2361: PetscReal val = 0.;
2363: for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
2364: v[i * dim + c] = val;
2365: }
2366: }
2367: if (J) {
2368: PetscReal extJ[18];
2370: extJ[0] = 0.;
2371: extJ[1] = 0.;
2372: extJ[2] = 0.;
2373: extJ[3] = 0.;
2374: extJ[4] = 1.;
2375: extJ[5] = 0.;
2376: extJ[6] = 1.;
2377: extJ[7] = 0.;
2378: extJ[8] = 0.;
2379: extJ[9] = 0.;
2380: extJ[10] = 0.;
2381: extJ[11] = 1.;
2382: extJ[12] = zeta;
2383: extJ[13] = 0.;
2384: extJ[14] = xi;
2385: extJ[15] = 0.;
2386: extJ[16] = zeta;
2387: extJ[17] = eta;
2389: for (j = 0; j < dim; j++) {
2390: for (k = 0; k < dimR; k++) {
2391: PetscReal val = 0.;
2393: for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
2394: J[i * dim * dim + dim * j + k] = val;
2395: }
2396: }
2397: DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2398: if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2399: }
2400: }
2401: }
2402: PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2403: PetscFunctionReturn(PETSC_SUCCESS);
2404: }
2406: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2407: {
2408: DMPolytopeType ct;
2409: PetscInt depth, dim, coordDim, coneSize, i;
2410: PetscInt Nq = 0;
2411: const PetscReal *points = NULL;
2412: DMLabel depthLabel;
2413: PetscReal xi0[3] = {-1., -1., -1.}, v0[3], J0[9], detJ0;
2414: PetscBool isAffine = PETSC_TRUE;
2416: PetscFunctionBegin;
2417: PetscCall(DMPlexGetDepth(dm, &depth));
2418: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
2419: PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
2420: PetscCall(DMLabelGetValue(depthLabel, cell, &dim));
2421: if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim));
2422: PetscCall(DMGetCoordinateDim(dm, &coordDim));
2423: PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim);
2424: if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL));
2425: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2426: switch (ct) {
2427: case DM_POLYTOPE_POINT:
2428: PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ));
2429: isAffine = PETSC_FALSE;
2430: break;
2431: case DM_POLYTOPE_SEGMENT:
2432: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2433: if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2434: else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ));
2435: break;
2436: case DM_POLYTOPE_TRIANGLE:
2437: if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2438: else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ));
2439: break;
2440: case DM_POLYTOPE_QUADRILATERAL:
2441: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ));
2442: isAffine = PETSC_FALSE;
2443: break;
2444: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2445: PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ));
2446: isAffine = PETSC_FALSE;
2447: break;
2448: case DM_POLYTOPE_TETRAHEDRON:
2449: if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2450: else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ));
2451: break;
2452: case DM_POLYTOPE_HEXAHEDRON:
2453: PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2454: isAffine = PETSC_FALSE;
2455: break;
2456: case DM_POLYTOPE_TRI_PRISM:
2457: PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2458: isAffine = PETSC_FALSE;
2459: break;
2460: default:
2461: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[PetscMax(0, PetscMin(ct, DM_NUM_POLYTOPES))]);
2462: }
2463: if (isAffine && Nq) {
2464: if (v) {
2465: for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
2466: }
2467: if (detJ) {
2468: for (i = 0; i < Nq; i++) detJ[i] = detJ0;
2469: }
2470: if (J) {
2471: PetscInt k;
2473: for (i = 0, k = 0; i < Nq; i++) {
2474: for (PetscInt j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
2475: }
2476: }
2477: if (invJ) {
2478: PetscInt k;
2479: switch (coordDim) {
2480: case 0:
2481: break;
2482: case 1:
2483: invJ[0] = 1. / J0[0];
2484: break;
2485: case 2:
2486: DMPlex_Invert2D_Internal(invJ, J0, detJ0);
2487: break;
2488: case 3:
2489: DMPlex_Invert3D_Internal(invJ, J0, detJ0);
2490: break;
2491: }
2492: for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
2493: for (PetscInt j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
2494: }
2495: }
2496: }
2497: PetscFunctionReturn(PETSC_SUCCESS);
2498: }
2500: /*@C
2501: DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell
2503: Collective
2505: Input Parameters:
2506: + dm - the `DMPLEX`
2507: - cell - the cell
2509: Output Parameters:
2510: + v0 - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
2511: . J - the Jacobian of the transform from the reference element
2512: . invJ - the inverse of the Jacobian
2513: - detJ - the Jacobian determinant
2515: Level: advanced
2517: .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2518: @*/
2519: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2520: {
2521: PetscFunctionBegin;
2522: PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ));
2523: PetscFunctionReturn(PETSC_SUCCESS);
2524: }
2526: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2527: {
2528: const PetscScalar *array;
2529: PetscScalar *coords = NULL;
2530: PetscInt numCoords;
2531: PetscBool isDG;
2532: PetscQuadrature feQuad;
2533: const PetscReal *quadPoints;
2534: PetscTabulation T;
2535: PetscInt dim, cdim, pdim, qdim, Nq, q;
2537: PetscFunctionBegin;
2538: PetscCall(DMGetDimension(dm, &dim));
2539: PetscCall(DMGetCoordinateDim(dm, &cdim));
2540: PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2541: if (!quad) { /* use the first point of the first functional of the dual space */
2542: PetscDualSpace dsp;
2544: PetscCall(PetscFEGetDualSpace(fe, &dsp));
2545: PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad));
2546: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2547: Nq = 1;
2548: } else {
2549: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2550: }
2551: PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim);
2552: PetscCall(PetscFEGetDimension(fe, &pdim));
2553: PetscCall(PetscFEGetQuadrature(fe, &feQuad));
2554: if (feQuad == quad) {
2555: PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T));
2556: PetscCheck(numCoords == pdim * cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %" PetscInt_FMT " coordinates for point %" PetscInt_FMT " != %" PetscInt_FMT "*%" PetscInt_FMT, numCoords, point, pdim, cdim);
2557: } else {
2558: PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T));
2559: }
2560: {
2561: const PetscReal *basis = T->T[0];
2562: const PetscReal *basisDer = T->T[1];
2563: PetscReal detJt;
2565: PetscAssert(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np);
2566: PetscAssert(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb);
2567: PetscAssert(cdim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->Nc);
2568: PetscAssert(dim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->cdim);
2569: if (v) {
2570: PetscCall(PetscArrayzero(v, Nq * cdim));
2571: for (q = 0; q < Nq; ++q) {
2572: PetscInt i, k;
2574: for (k = 0; k < pdim; ++k) {
2575: const PetscInt vertex = k / cdim;
2576: for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2577: }
2578: PetscCall(PetscLogFlops(2.0 * pdim * cdim));
2579: }
2580: }
2581: if (J) {
2582: PetscCall(PetscArrayzero(J, Nq * cdim * cdim));
2583: for (q = 0; q < Nq; ++q) {
2584: PetscInt i, j, k, c, r;
2586: /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2587: for (k = 0; k < pdim; ++k) {
2588: const PetscInt vertex = k / cdim;
2589: for (j = 0; j < dim; ++j) {
2590: for (i = 0; i < cdim; ++i) J[(q * cdim + i) * cdim + j] += basisDer[((q * pdim + k) * cdim + i) * dim + j] * PetscRealPart(coords[vertex * cdim + i]);
2591: }
2592: }
2593: PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim));
2594: if (cdim > dim) {
2595: for (c = dim; c < cdim; ++c)
2596: for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2597: }
2598: if (!detJ && !invJ) continue;
2599: detJt = 0.;
2600: switch (cdim) {
2601: case 3:
2602: DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2603: if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2604: break;
2605: case 2:
2606: DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2607: if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2608: break;
2609: case 1:
2610: detJt = J[q * cdim * dim];
2611: if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2612: }
2613: if (detJ) detJ[q] = detJt;
2614: }
2615: } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
2616: }
2617: if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T));
2618: PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2619: PetscFunctionReturn(PETSC_SUCCESS);
2620: }
2622: /*@C
2623: DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell
2625: Collective
2627: Input Parameters:
2628: + dm - the `DMPLEX`
2629: . cell - the cell
2630: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated. If `quad` is `NULL`, geometry will be
2631: evaluated at the first vertex of the reference element
2633: Output Parameters:
2634: + v - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element. This is a
2635: one-dimensional array of size $cdim * Nq$ where $cdim$ is the dimension of the `DM` coordinate space and $Nq$ is the number of quadrature points
2636: . J - the Jacobian of the transform from the reference element at each quadrature point. This is a one-dimensional array of size $Nq * cdim * cdim$ containing
2637: each Jacobian in column-major order.
2638: . invJ - the inverse of the Jacobian at each quadrature point. This is a one-dimensional array of size $Nq * cdim * cdim$ containing
2639: each inverse Jacobian in column-major order.
2640: - detJ - the Jacobian determinant at each quadrature point. This is a one-dimensional array of size $Nq$.
2642: Level: advanced
2644: Note:
2645: Implicit cell geometry must be used when the topological mesh dimension is not equal to the coordinate dimension, for instance for embedded manifolds.
2647: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2648: @*/
2649: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal detJ[])
2650: {
2651: DM cdm;
2652: PetscFE fe = NULL;
2653: PetscInt dim, cdim;
2655: PetscFunctionBegin;
2656: PetscAssertPointer(detJ, 7);
2657: PetscCall(DMGetDimension(dm, &dim));
2658: PetscCall(DMGetCoordinateDim(dm, &cdim));
2659: PetscCall(DMGetCoordinateDM(dm, &cdm));
2660: if (cdm) {
2661: PetscSpace sp;
2662: PetscClassId id;
2663: PetscDS prob;
2664: PetscObject disc;
2665: PetscInt Nf, spdim, qdim;
2667: PetscCall(DMGetNumFields(cdm, &Nf));
2668: if (Nf) {
2669: PetscCall(DMGetDS(cdm, &prob));
2670: PetscCall(PetscDSGetDiscretization(prob, 0, &disc));
2671: PetscCall(PetscObjectGetClassId(disc, &id));
2672: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2673: if (fe && quad) {
2674: PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, NULL, NULL, NULL));
2675: PetscCall(PetscFEGetBasisSpace(fe, &sp));
2676: PetscCall(PetscSpaceGetNumVariables(sp, &spdim));
2677: if (qdim != spdim) fe = NULL;
2678: }
2679: }
2680: }
2681: if (!fe || (dim != cdim)) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ));
2682: else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ));
2683: PetscFunctionReturn(PETSC_SUCCESS);
2684: }
2686: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2687: {
2688: PetscSection coordSection;
2689: Vec coordinates;
2690: const PetscScalar *coords = NULL;
2691: PetscInt d, dof, off;
2693: PetscFunctionBegin;
2694: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2695: PetscCall(DMGetCoordinateSection(dm, &coordSection));
2696: PetscCall(VecGetArrayRead(coordinates, &coords));
2698: /* for a point the centroid is just the coord */
2699: if (centroid) {
2700: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2701: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2702: for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2703: }
2704: if (normal) {
2705: const PetscInt *support, *cones;
2706: PetscInt supportSize;
2707: PetscReal norm, sign;
2709: /* compute the norm based upon the support centroids */
2710: PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize));
2711: PetscCall(DMPlexGetSupport(dm, cell, &support));
2712: PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL));
2714: /* Take the normal from the centroid of the support to the vertex*/
2715: PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2716: PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2717: for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]);
2719: /* Determine the sign of the normal based upon its location in the support */
2720: PetscCall(DMPlexGetCone(dm, support[0], &cones));
2721: sign = cones[0] == cell ? 1.0 : -1.0;
2723: norm = DMPlex_NormD_Internal(dim, normal);
2724: for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2725: }
2726: if (vol) *vol = 1.0;
2727: PetscCall(VecRestoreArrayRead(coordinates, &coords));
2728: PetscFunctionReturn(PETSC_SUCCESS);
2729: }
2731: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2732: {
2733: const PetscScalar *array;
2734: PetscScalar *coords = NULL;
2735: PetscInt cdim, coordSize, d;
2736: PetscBool isDG;
2738: PetscFunctionBegin;
2739: PetscCall(DMGetCoordinateDim(dm, &cdim));
2740: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2741: PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2);
2742: if (centroid) {
2743: for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2744: }
2745: if (normal) {
2746: PetscReal norm;
2748: switch (cdim) {
2749: case 3:
2750: normal[2] = 0.; /* fall through */
2751: case 2:
2752: normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2753: normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2754: break;
2755: case 1:
2756: normal[0] = 1.0;
2757: break;
2758: default:
2759: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2760: }
2761: norm = DMPlex_NormD_Internal(cdim, normal);
2762: for (d = 0; d < cdim; ++d) normal[d] /= norm;
2763: }
2764: if (vol) {
2765: *vol = 0.0;
2766: for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2767: *vol = PetscSqrtReal(*vol);
2768: }
2769: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2770: PetscFunctionReturn(PETSC_SUCCESS);
2771: }
2773: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2774: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2775: {
2776: DMPolytopeType ct;
2777: const PetscScalar *array;
2778: PetscScalar *coords = NULL;
2779: PetscInt coordSize;
2780: PetscBool isDG;
2781: PetscInt fv[4] = {0, 1, 2, 3};
2782: PetscInt cdim, numCorners, p, d;
2784: PetscFunctionBegin;
2785: /* Must check for hybrid cells because prisms have a different orientation scheme */
2786: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2787: switch (ct) {
2788: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2789: fv[2] = 3;
2790: fv[3] = 2;
2791: break;
2792: default:
2793: break;
2794: }
2795: PetscCall(DMGetCoordinateDim(dm, &cdim));
2796: PetscCall(DMPlexGetConeSize(dm, cell, &numCorners));
2797: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2798: {
2799: PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;
2801: for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2802: for (p = 0; p < numCorners - 2; ++p) {
2803: PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2804: for (d = 0; d < cdim; d++) {
2805: e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2806: e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2807: }
2808: const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2809: const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2810: const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2811: const PetscReal a = PetscSqrtReal(dx * dx + dy * dy + dz * dz);
2813: n[0] += dx;
2814: n[1] += dy;
2815: n[2] += dz;
2816: for (d = 0; d < cdim; d++) c[d] += a * PetscRealPart(origin[d] + coords[cdim * fv[p + 1] + d] + coords[cdim * fv[p + 2] + d]) / 3.;
2817: }
2818: norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2819: // Allow zero volume cells
2820: if (norm != 0) {
2821: n[0] /= norm;
2822: n[1] /= norm;
2823: n[2] /= norm;
2824: c[0] /= norm;
2825: c[1] /= norm;
2826: c[2] /= norm;
2827: }
2828: if (vol) *vol = 0.5 * norm;
2829: if (centroid)
2830: for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2831: if (normal)
2832: for (d = 0; d < cdim; ++d) normal[d] = n[d];
2833: }
2834: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2835: PetscFunctionReturn(PETSC_SUCCESS);
2836: }
2838: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2839: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2840: {
2841: DMPolytopeType ct;
2842: const PetscScalar *array;
2843: PetscScalar *coords = NULL;
2844: PetscInt coordSize;
2845: PetscBool isDG;
2846: PetscReal vsum = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2847: const PetscInt order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2848: const PetscInt *cone, *faceSizes, *faces;
2849: const DMPolytopeType *faceTypes;
2850: PetscBool isHybrid = PETSC_FALSE;
2851: PetscInt numFaces, f, fOff = 0, p, d;
2853: PetscFunctionBegin;
2854: PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim);
2855: /* Must check for hybrid cells because prisms have a different orientation scheme */
2856: PetscCall(DMPlexGetCellType(dm, cell, &ct));
2857: switch (ct) {
2858: case DM_POLYTOPE_POINT_PRISM_TENSOR:
2859: case DM_POLYTOPE_SEG_PRISM_TENSOR:
2860: case DM_POLYTOPE_TRI_PRISM_TENSOR:
2861: case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2862: isHybrid = PETSC_TRUE;
2863: default:
2864: break;
2865: }
2867: if (centroid)
2868: for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2869: PetscCall(DMPlexGetCone(dm, cell, &cone));
2871: // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates
2872: PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2873: PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2874: for (f = 0; f < numFaces; ++f) {
2875: PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
2877: // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2878: // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2879: // so that all tetrahedra have positive volume.
2880: if (f == 0)
2881: for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2882: switch (faceTypes[f]) {
2883: case DM_POLYTOPE_TRIANGLE:
2884: for (d = 0; d < dim; ++d) {
2885: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2886: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2887: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2888: }
2889: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2890: if (flip) vtmp = -vtmp;
2891: vsum += vtmp;
2892: if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2893: for (d = 0; d < dim; ++d) {
2894: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2895: }
2896: }
2897: break;
2898: case DM_POLYTOPE_QUADRILATERAL:
2899: case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2900: PetscInt fv[4] = {0, 1, 2, 3};
2902: /* Side faces for hybrid cells are stored as tensor products */
2903: if (isHybrid && f > 1) {
2904: fv[2] = 3;
2905: fv[3] = 2;
2906: }
2907: /* DO FOR PYRAMID */
2908: /* First tet */
2909: for (d = 0; d < dim; ++d) {
2910: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2911: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2912: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2913: }
2914: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2915: if (flip) vtmp = -vtmp;
2916: vsum += vtmp;
2917: if (centroid) {
2918: for (d = 0; d < dim; ++d) {
2919: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2920: }
2921: }
2922: /* Second tet */
2923: for (d = 0; d < dim; ++d) {
2924: coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2925: coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2926: coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2927: }
2928: Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2929: if (flip) vtmp = -vtmp;
2930: vsum += vtmp;
2931: if (centroid) {
2932: for (d = 0; d < dim; ++d) {
2933: for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2934: }
2935: }
2936: break;
2937: }
2938: default:
2939: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2940: }
2941: fOff += faceSizes[f];
2942: }
2943: PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2944: PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2945: if (vol) *vol = PetscAbsReal(vsum);
2946: if (normal)
2947: for (d = 0; d < dim; ++d) normal[d] = 0.0;
2948: if (centroid)
2949: for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2950: PetscFunctionReturn(PETSC_SUCCESS);
2951: }
2953: /*@C
2954: DMPlexComputeCellGeometryFVM - Compute the volume for a given cell
2956: Collective
2958: Input Parameters:
2959: + dm - the `DMPLEX`
2960: - cell - the cell
2962: Output Parameters:
2963: + vol - the cell volume
2964: . centroid - the cell centroid
2965: - normal - the cell normal, if appropriate
2967: Level: advanced
2969: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2970: @*/
2971: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2972: {
2973: PetscInt depth, dim;
2975: PetscFunctionBegin;
2976: PetscCall(DMPlexGetDepth(dm, &depth));
2977: PetscCall(DMGetDimension(dm, &dim));
2978: PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
2979: switch (depth) {
2980: case 0:
2981: PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal));
2982: break;
2983: case 1:
2984: PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal));
2985: break;
2986: case 2:
2987: PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal));
2988: break;
2989: case 3:
2990: PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal));
2991: break;
2992: default:
2993: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2994: }
2995: PetscFunctionReturn(PETSC_SUCCESS);
2996: }
2998: /*@
2999: DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method
3001: Input Parameter:
3002: . dm - The `DMPLEX`
3004: Output Parameters:
3005: + cellgeom - A `Vec` of `PetscFVCellGeom` data
3006: - facegeom - A `Vec` of `PetscFVFaceGeom` data
3008: Level: developer
3010: .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom`
3011: @*/
3012: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
3013: {
3014: DM dmFace, dmCell;
3015: DMLabel ghostLabel;
3016: PetscSection sectionFace, sectionCell;
3017: PetscSection coordSection;
3018: Vec coordinates;
3019: PetscScalar *fgeom, *cgeom;
3020: PetscReal minradius;
3021: PetscInt dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;
3023: PetscFunctionBegin;
3024: PetscCall(DMGetDimension(dm, &dim));
3025: PetscCall(DMGetCoordinateSection(dm, &coordSection));
3026: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
3027: /* Make cell centroids and volumes */
3028: PetscCall(DMClone(dm, &dmCell));
3029: PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection));
3030: PetscCall(DMSetCoordinatesLocal(dmCell, coordinates));
3031: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell));
3032: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3033: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3034: PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
3035: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar))));
3036: PetscCall(PetscSectionSetUp(sectionCell));
3037: PetscCall(DMSetLocalSection(dmCell, sectionCell));
3038: PetscCall(PetscSectionDestroy(§ionCell));
3039: PetscCall(DMCreateLocalVector(dmCell, cellgeom));
3040: if (cEndInterior < 0) cEndInterior = cEnd;
3041: PetscCall(VecGetArray(*cellgeom, &cgeom));
3042: for (c = cStart; c < cEndInterior; ++c) {
3043: PetscFVCellGeom *cg;
3045: PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg));
3046: PetscCall(PetscArrayzero(cg, 1));
3047: PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL));
3048: }
3049: /* Compute face normals and minimum cell radius */
3050: PetscCall(DMClone(dm, &dmFace));
3051: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionFace));
3052: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
3053: PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd));
3054: for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar))));
3055: PetscCall(PetscSectionSetUp(sectionFace));
3056: PetscCall(DMSetLocalSection(dmFace, sectionFace));
3057: PetscCall(PetscSectionDestroy(§ionFace));
3058: PetscCall(DMCreateLocalVector(dmFace, facegeom));
3059: PetscCall(VecGetArray(*facegeom, &fgeom));
3060: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3061: minradius = PETSC_MAX_REAL;
3062: for (f = fStart; f < fEnd; ++f) {
3063: PetscFVFaceGeom *fg;
3064: PetscReal area;
3065: const PetscInt *cells;
3066: PetscInt ncells, ghost = -1, d, numChildren;
3068: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3069: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3070: PetscCall(DMPlexGetSupport(dm, f, &cells));
3071: PetscCall(DMPlexGetSupportSize(dm, f, &ncells));
3072: /* It is possible to get a face with no support when using partition overlap */
3073: if (!ncells || ghost >= 0 || numChildren) continue;
3074: PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg));
3075: PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal));
3076: for (d = 0; d < dim; ++d) fg->normal[d] *= area;
3077: /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
3078: {
3079: PetscFVCellGeom *cL, *cR;
3080: PetscReal *lcentroid, *rcentroid;
3081: PetscReal l[3], r[3], v[3];
3083: PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL));
3084: lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
3085: if (ncells > 1) {
3086: PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR));
3087: rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
3088: } else {
3089: rcentroid = fg->centroid;
3090: }
3091: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l));
3092: PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r));
3093: DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
3094: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
3095: for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
3096: }
3097: if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
3098: PetscCheck(dim != 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g) v (%g,%g)", f, (double)fg->normal[0], (double)fg->normal[1], (double)v[0], (double)v[1]);
3099: PetscCheck(dim != 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double)fg->normal[0], (double)fg->normal[1], (double)fg->normal[2], (double)v[0], (double)v[1], (double)v[2]);
3100: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f);
3101: }
3102: if (cells[0] < cEndInterior) {
3103: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
3104: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
3105: }
3106: if (ncells > 1 && cells[1] < cEndInterior) {
3107: DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
3108: minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
3109: }
3110: }
3111: }
3112: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &minradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
3113: PetscCall(DMPlexSetMinRadius(dm, minradius));
3114: /* Compute centroids of ghost cells */
3115: for (c = cEndInterior; c < cEnd; ++c) {
3116: PetscFVFaceGeom *fg;
3117: const PetscInt *cone, *support;
3118: PetscInt coneSize, supportSize, s;
3120: PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize));
3121: PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize);
3122: PetscCall(DMPlexGetCone(dmCell, c, &cone));
3123: PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize));
3124: PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize);
3125: PetscCall(DMPlexGetSupport(dmCell, cone[0], &support));
3126: PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg));
3127: for (s = 0; s < 2; ++s) {
3128: /* Reflect ghost centroid across plane of face */
3129: if (support[s] == c) {
3130: PetscFVCellGeom *ci;
3131: PetscFVCellGeom *cg;
3132: PetscReal c2f[3], a;
3134: PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci));
3135: DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
3136: a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
3137: PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg));
3138: DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid);
3139: cg->volume = ci->volume;
3140: }
3141: }
3142: }
3143: PetscCall(VecRestoreArray(*facegeom, &fgeom));
3144: PetscCall(VecRestoreArray(*cellgeom, &cgeom));
3145: PetscCall(DMDestroy(&dmCell));
3146: PetscCall(DMDestroy(&dmFace));
3147: PetscFunctionReturn(PETSC_SUCCESS);
3148: }
3150: /*@
3151: DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face
3153: Not Collective
3155: Input Parameter:
3156: . dm - the `DMPLEX`
3158: Output Parameter:
3159: . minradius - the minimum cell radius
3161: Level: developer
3163: .seealso: `DMPLEX`, `DMGetCoordinates()`
3164: @*/
3165: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
3166: {
3167: PetscFunctionBegin;
3169: PetscAssertPointer(minradius, 2);
3170: *minradius = ((DM_Plex *)dm->data)->minradius;
3171: PetscFunctionReturn(PETSC_SUCCESS);
3172: }
3174: /*@
3175: DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face
3177: Logically Collective
3179: Input Parameters:
3180: + dm - the `DMPLEX`
3181: - minradius - the minimum cell radius
3183: Level: developer
3185: .seealso: `DMPLEX`, `DMSetCoordinates()`
3186: @*/
3187: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
3188: {
3189: PetscFunctionBegin;
3191: ((DM_Plex *)dm->data)->minradius = minradius;
3192: PetscFunctionReturn(PETSC_SUCCESS);
3193: }
3195: /*@C
3196: DMPlexGetCoordinateMap - Returns the function used to map coordinates of newly generated mesh points
3198: Not Collective
3200: Input Parameter:
3201: . dm - the `DMPLEX`
3203: Output Parameter:
3204: . coordFunc - the mapping function
3206: Level: developer
3208: Note:
3209: This function maps from the generated coordinate for the new point to the actual coordinate. Thus it is only practical for manifolds with a nice analytical definition that you can get to from any starting point, like a sphere,
3211: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMPlexSetCoordinateMap()`, `PetscPointFn`
3212: @*/
3213: PetscErrorCode DMPlexGetCoordinateMap(DM dm, PetscPointFn **coordFunc)
3214: {
3215: PetscFunctionBegin;
3217: PetscAssertPointer(coordFunc, 2);
3218: *coordFunc = ((DM_Plex *)dm->data)->coordFunc;
3219: PetscFunctionReturn(PETSC_SUCCESS);
3220: }
3222: /*@C
3223: DMPlexSetCoordinateMap - Sets the function used to map coordinates of newly generated mesh points
3225: Logically Collective
3227: Input Parameters:
3228: + dm - the `DMPLEX`
3229: - coordFunc - the mapping function
3231: Level: developer
3233: Note:
3234: This function maps from the generated coordinate for the new point to the actual coordinate. Thus it is only practical for manifolds with a nice analytical definition that you can get to from any starting point, like a sphere,
3236: .seealso: `DMPLEX`, `DMSetCoordinates()`, `DMPlexGetCoordinateMap()`, `PetscPointFn`
3237: @*/
3238: PetscErrorCode DMPlexSetCoordinateMap(DM dm, PetscPointFn *coordFunc)
3239: {
3240: PetscFunctionBegin;
3242: ((DM_Plex *)dm->data)->coordFunc = coordFunc;
3243: PetscFunctionReturn(PETSC_SUCCESS);
3244: }
3246: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3247: {
3248: DMLabel ghostLabel;
3249: PetscScalar *dx, *grad, **gref;
3250: PetscInt dim, cStart, cEnd, c, cEndInterior, maxNumFaces;
3252: PetscFunctionBegin;
3253: PetscCall(DMGetDimension(dm, &dim));
3254: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3255: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3256: cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
3257: PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL));
3258: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3259: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3260: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3261: for (c = cStart; c < cEndInterior; c++) {
3262: const PetscInt *faces;
3263: PetscInt numFaces, usedFaces, f, d;
3264: PetscFVCellGeom *cg;
3265: PetscBool boundary;
3266: PetscInt ghost;
3268: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
3269: PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3270: if (ghost >= 0) continue;
3272: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3273: PetscCall(DMPlexGetConeSize(dm, c, &numFaces));
3274: PetscCall(DMPlexGetCone(dm, c, &faces));
3275: PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces);
3276: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3277: PetscFVCellGeom *cg1;
3278: PetscFVFaceGeom *fg;
3279: const PetscInt *fcells;
3280: PetscInt ncell, side;
3282: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3283: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3284: if ((ghost >= 0) || boundary) continue;
3285: PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
3286: side = (c != fcells[0]); /* c is on left=0 or right=1 of face */
3287: ncell = fcells[!side]; /* the neighbor */
3288: PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg));
3289: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3290: for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
3291: gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
3292: }
3293: PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
3294: PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad));
3295: for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3296: PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3297: PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3298: if ((ghost >= 0) || boundary) continue;
3299: for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
3300: ++usedFaces;
3301: }
3302: }
3303: PetscCall(PetscFree3(dx, grad, gref));
3304: PetscFunctionReturn(PETSC_SUCCESS);
3305: }
3307: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3308: {
3309: DMLabel ghostLabel;
3310: PetscScalar *dx, *grad, **gref;
3311: PetscInt dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
3312: PetscSection neighSec;
3313: PetscInt (*neighbors)[2];
3314: PetscInt *counter;
3316: PetscFunctionBegin;
3317: PetscCall(DMGetDimension(dm, &dim));
3318: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3319: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3320: if (cEndInterior < 0) cEndInterior = cEnd;
3321: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec));
3322: PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior));
3323: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
3324: PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3325: for (f = fStart; f < fEnd; f++) {
3326: const PetscInt *fcells;
3327: PetscBool boundary;
3328: PetscInt ghost = -1;
3329: PetscInt numChildren, numCells, c;
3331: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3332: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3333: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3334: if ((ghost >= 0) || boundary || numChildren) continue;
3335: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3336: if (numCells == 2) {
3337: PetscCall(DMPlexGetSupport(dm, f, &fcells));
3338: for (c = 0; c < 2; c++) {
3339: PetscInt cell = fcells[c];
3341: if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1));
3342: }
3343: }
3344: }
3345: PetscCall(PetscSectionSetUp(neighSec));
3346: PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces));
3347: PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3348: nStart = 0;
3349: PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd));
3350: PetscCall(PetscMalloc1(nEnd - nStart, &neighbors));
3351: PetscCall(PetscCalloc1(cEndInterior - cStart, &counter));
3352: for (f = fStart; f < fEnd; f++) {
3353: const PetscInt *fcells;
3354: PetscBool boundary;
3355: PetscInt ghost = -1;
3356: PetscInt numChildren, numCells, c;
3358: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3359: PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3360: PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3361: if ((ghost >= 0) || boundary || numChildren) continue;
3362: PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3363: if (numCells == 2) {
3364: PetscCall(DMPlexGetSupport(dm, f, &fcells));
3365: for (c = 0; c < 2; c++) {
3366: PetscInt cell = fcells[c], off;
3368: if (cell >= cStart && cell < cEndInterior) {
3369: PetscCall(PetscSectionGetOffset(neighSec, cell, &off));
3370: off += counter[cell - cStart]++;
3371: neighbors[off][0] = f;
3372: neighbors[off][1] = fcells[1 - c];
3373: }
3374: }
3375: }
3376: }
3377: PetscCall(PetscFree(counter));
3378: PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3379: for (c = cStart; c < cEndInterior; c++) {
3380: PetscInt numFaces, f, d, off, ghost = -1;
3381: PetscFVCellGeom *cg;
3383: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3384: PetscCall(PetscSectionGetDof(neighSec, c, &numFaces));
3385: PetscCall(PetscSectionGetOffset(neighSec, c, &off));
3387: // do not attempt to compute a gradient reconstruction stencil in a ghost cell. It will never be used
3388: if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3389: if (ghost >= 0) continue;
3391: PetscCheck(numFaces >= dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cell %" PetscInt_FMT " has only %" PetscInt_FMT " faces, not enough for gradient reconstruction", c, numFaces);
3392: for (f = 0; f < numFaces; ++f) {
3393: PetscFVCellGeom *cg1;
3394: PetscFVFaceGeom *fg;
3395: const PetscInt *fcells;
3396: PetscInt ncell, side, nface;
3398: nface = neighbors[off + f][0];
3399: ncell = neighbors[off + f][1];
3400: PetscCall(DMPlexGetSupport(dm, nface, &fcells));
3401: side = (c != fcells[0]);
3402: PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg));
3403: PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3404: for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d];
3405: gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
3406: }
3407: PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad));
3408: for (f = 0; f < numFaces; ++f) {
3409: for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d];
3410: }
3411: }
3412: PetscCall(PetscFree3(dx, grad, gref));
3413: PetscCall(PetscSectionDestroy(&neighSec));
3414: PetscCall(PetscFree(neighbors));
3415: PetscFunctionReturn(PETSC_SUCCESS);
3416: }
3418: /*@
3419: DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data
3421: Collective
3423: Input Parameters:
3424: + dm - The `DMPLEX`
3425: . fvm - The `PetscFV`
3426: - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()`
3428: Input/Output Parameter:
3429: . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output
3430: the geometric factors for gradient calculation are inserted
3432: Output Parameter:
3433: . dmGrad - The `DM` describing the layout of gradient data
3435: Level: developer
3437: .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
3438: @*/
3439: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
3440: {
3441: DM dmFace, dmCell;
3442: PetscScalar *fgeom, *cgeom;
3443: PetscSection sectionGrad, parentSection;
3444: PetscInt dim, pdim, cStart, cEnd, cEndInterior, c;
3446: PetscFunctionBegin;
3447: PetscCall(DMGetDimension(dm, &dim));
3448: PetscCall(PetscFVGetNumComponents(fvm, &pdim));
3449: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3450: PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3451: /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
3452: PetscCall(VecGetDM(faceGeometry, &dmFace));
3453: PetscCall(VecGetDM(cellGeometry, &dmCell));
3454: PetscCall(VecGetArray(faceGeometry, &fgeom));
3455: PetscCall(VecGetArray(cellGeometry, &cgeom));
3456: PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL));
3457: if (!parentSection) {
3458: PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3459: } else {
3460: PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3461: }
3462: PetscCall(VecRestoreArray(faceGeometry, &fgeom));
3463: PetscCall(VecRestoreArray(cellGeometry, &cgeom));
3464: /* Create storage for gradients */
3465: PetscCall(DMClone(dm, dmGrad));
3466: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionGrad));
3467: PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd));
3468: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim));
3469: PetscCall(PetscSectionSetUp(sectionGrad));
3470: PetscCall(DMSetLocalSection(*dmGrad, sectionGrad));
3471: PetscCall(PetscSectionDestroy(§ionGrad));
3472: PetscFunctionReturn(PETSC_SUCCESS);
3473: }
3475: /*@
3476: DMPlexGetDataFVM - Retrieve precomputed cell geometry
3478: Collective
3480: Input Parameters:
3481: + dm - The `DM`
3482: - fv - The `PetscFV`
3484: Output Parameters:
3485: + cellgeom - The cell geometry
3486: . facegeom - The face geometry
3487: - gradDM - The gradient matrices
3489: Level: developer
3491: .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()`
3492: @*/
3493: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
3494: {
3495: PetscObject cellgeomobj, facegeomobj;
3497: PetscFunctionBegin;
3498: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3499: if (!cellgeomobj) {
3500: Vec cellgeomInt, facegeomInt;
3502: PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt));
3503: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt));
3504: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt));
3505: PetscCall(VecDestroy(&cellgeomInt));
3506: PetscCall(VecDestroy(&facegeomInt));
3507: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3508: }
3509: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj));
3510: if (cellgeom) *cellgeom = (Vec)cellgeomobj;
3511: if (facegeom) *facegeom = (Vec)facegeomobj;
3512: if (gradDM) {
3513: PetscObject gradobj;
3514: PetscBool computeGradients;
3516: PetscCall(PetscFVGetComputeGradients(fv, &computeGradients));
3517: if (!computeGradients) {
3518: *gradDM = NULL;
3519: PetscFunctionReturn(PETSC_SUCCESS);
3520: }
3521: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3522: if (!gradobj) {
3523: DM dmGradInt;
3525: PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt));
3526: PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt));
3527: PetscCall(DMDestroy(&dmGradInt));
3528: PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3529: }
3530: *gradDM = (DM)gradobj;
3531: }
3532: PetscFunctionReturn(PETSC_SUCCESS);
3533: }
3535: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
3536: {
3537: PetscInt l, m;
3539: PetscFunctionBeginHot;
3540: if (dimC == dimR && dimR <= 3) {
3541: /* invert Jacobian, multiply */
3542: PetscScalar det, idet;
3544: switch (dimR) {
3545: case 1:
3546: invJ[0] = 1. / J[0];
3547: break;
3548: case 2:
3549: det = J[0] * J[3] - J[1] * J[2];
3550: idet = 1. / det;
3551: invJ[0] = J[3] * idet;
3552: invJ[1] = -J[1] * idet;
3553: invJ[2] = -J[2] * idet;
3554: invJ[3] = J[0] * idet;
3555: break;
3556: case 3: {
3557: invJ[0] = J[4] * J[8] - J[5] * J[7];
3558: invJ[1] = J[2] * J[7] - J[1] * J[8];
3559: invJ[2] = J[1] * J[5] - J[2] * J[4];
3560: det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
3561: idet = 1. / det;
3562: invJ[0] *= idet;
3563: invJ[1] *= idet;
3564: invJ[2] *= idet;
3565: invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
3566: invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
3567: invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
3568: invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
3569: invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
3570: invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
3571: } break;
3572: }
3573: for (l = 0; l < dimR; l++) {
3574: for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
3575: }
3576: } else {
3577: char transpose = PetscDefined(USE_COMPLEX) ? 'C' : 'T';
3578: PetscBLASInt m, n, one = 1, worksize;
3580: PetscCall(PetscBLASIntCast(dimR, &m));
3581: PetscCall(PetscBLASIntCast(dimC, &n));
3582: PetscCall(PetscBLASIntCast(dimC * dimC, &worksize));
3583: for (l = 0; l < dimC; l++) invJ[l] = resNeg[l];
3585: PetscCallLAPACKInfo("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info));
3586: for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3587: }
3588: PetscFunctionReturn(PETSC_SUCCESS);
3589: }
3591: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3592: {
3593: PetscInt coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
3594: PetscScalar *coordsScalar = NULL;
3595: PetscReal *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
3596: PetscScalar *J, *invJ, *work;
3598: PetscFunctionBegin;
3600: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3601: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3602: PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3603: PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3604: cellCoords = &cellData[0];
3605: cellCoeffs = &cellData[coordSize];
3606: extJ = &cellData[2 * coordSize];
3607: resNeg = &cellData[2 * coordSize + dimR];
3608: invJ = &J[dimR * dimC];
3609: work = &J[2 * dimR * dimC];
3610: if (dimR == 2) {
3611: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3613: for (i = 0; i < 4; i++) {
3614: PetscInt plexI = zToPlex[i];
3616: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3617: }
3618: } else if (dimR == 3) {
3619: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3621: for (i = 0; i < 8; i++) {
3622: PetscInt plexI = zToPlex[i];
3624: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3625: }
3626: } else {
3627: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3628: }
3629: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3630: for (i = 0; i < dimR; i++) {
3631: PetscReal *swap;
3633: for (j = 0; j < (numV / 2); j++) {
3634: for (k = 0; k < dimC; k++) {
3635: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3636: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3637: }
3638: }
3640: if (i < dimR - 1) {
3641: swap = cellCoeffs;
3642: cellCoeffs = cellCoords;
3643: cellCoords = swap;
3644: }
3645: }
3646: PetscCall(PetscArrayzero(refCoords, numPoints * dimR));
3647: for (j = 0; j < numPoints; j++) {
3648: for (i = 0; i < maxIts; i++) {
3649: PetscReal *guess = &refCoords[dimR * j];
3651: /* compute -residual and Jacobian */
3652: for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3653: for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3654: for (k = 0; k < numV; k++) {
3655: PetscReal extCoord = 1.;
3656: for (l = 0; l < dimR; l++) {
3657: PetscReal coord = guess[l];
3658: PetscInt dep = (k & (1 << l)) >> l;
3660: extCoord *= dep * coord + !dep;
3661: extJ[l] = dep;
3663: for (m = 0; m < dimR; m++) {
3664: PetscReal coord = guess[m];
3665: PetscInt dep = ((k & (1 << m)) >> m) && (m != l);
3666: PetscReal mult = dep * coord + !dep;
3668: extJ[l] *= mult;
3669: }
3670: }
3671: for (l = 0; l < dimC; l++) {
3672: PetscReal coeff = cellCoeffs[dimC * k + l];
3674: resNeg[l] -= coeff * extCoord;
3675: for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3676: }
3677: }
3678: if (0 && PetscDefined(USE_DEBUG)) {
3679: PetscReal maxAbs = 0.;
3681: for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3682: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3683: }
3685: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess));
3686: }
3687: }
3688: PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3689: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3690: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3691: PetscFunctionReturn(PETSC_SUCCESS);
3692: }
3694: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3695: {
3696: PetscInt coordSize, i, j, k, l, numV = (1 << dimR);
3697: PetscScalar *coordsScalar = NULL;
3698: PetscReal *cellData, *cellCoords, *cellCoeffs;
3700: PetscFunctionBegin;
3702: PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3703: PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3704: PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3705: cellCoords = &cellData[0];
3706: cellCoeffs = &cellData[coordSize];
3707: if (dimR == 2) {
3708: const PetscInt zToPlex[4] = {0, 1, 3, 2};
3710: for (i = 0; i < 4; i++) {
3711: PetscInt plexI = zToPlex[i];
3713: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3714: }
3715: } else if (dimR == 3) {
3716: const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
3718: for (i = 0; i < 8; i++) {
3719: PetscInt plexI = zToPlex[i];
3721: for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3722: }
3723: } else {
3724: for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3725: }
3726: /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3727: for (i = 0; i < dimR; i++) {
3728: PetscReal *swap;
3730: for (j = 0; j < (numV / 2); j++) {
3731: for (k = 0; k < dimC; k++) {
3732: cellCoeffs[dimC * j + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3733: cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3734: }
3735: }
3737: if (i < dimR - 1) {
3738: swap = cellCoeffs;
3739: cellCoeffs = cellCoords;
3740: cellCoords = swap;
3741: }
3742: }
3743: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_BEGIN("-Warray-bounds")
3744: PetscCall(PetscArrayzero(realCoords, numPoints * dimC));
3745: PETSC_PRAGMA_DIAGNOSTIC_IGNORED_END()
3746: for (j = 0; j < numPoints; j++) {
3747: const PetscReal *guess = &refCoords[dimR * j];
3748: PetscReal *mapped = &realCoords[dimC * j];
3750: for (k = 0; k < numV; k++) {
3751: PetscReal extCoord = 1.;
3752: for (l = 0; l < dimR; l++) {
3753: PetscReal coord = guess[l];
3754: PetscInt dep = (k & (1 << l)) >> l;
3756: extCoord *= dep * coord + !dep;
3757: }
3758: for (l = 0; l < dimC; l++) {
3759: PetscReal coeff = cellCoeffs[dimC * k + l];
3761: mapped[l] += coeff * extCoord;
3762: }
3763: }
3764: }
3765: PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3766: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3767: PetscFunctionReturn(PETSC_SUCCESS);
3768: }
3770: PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR, PetscInt maxIter, PetscReal *tol)
3771: {
3772: PetscInt numComp, pdim, i, j, k, l, m, coordSize;
3773: PetscScalar *nodes = NULL;
3774: PetscReal *invV, *modes;
3775: PetscReal *B, *D, *resNeg;
3776: PetscScalar *J, *invJ, *work;
3777: PetscReal tolerance = tol == NULL ? 0.0 : *tol;
3779: PetscFunctionBegin;
3780: PetscCall(PetscFEGetDimension(fe, &pdim));
3781: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3782: PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc);
3783: /* we shouldn't apply inverse closure permutation, if one exists */
3784: PetscCall(DMPlexVecGetOrientedClosure(dm, NULL, PETSC_FALSE, coords, cell, 0, &coordSize, &nodes));
3785: /* convert nodes to values in the stable evaluation basis */
3786: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3787: invV = fe->invV;
3788: for (i = 0; i < pdim; ++i) {
3789: modes[i] = 0.;
3790: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3791: }
3792: PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3793: D = &B[pdim * Nc];
3794: resNeg = &D[pdim * Nc * dimR];
3795: PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3796: invJ = &J[Nc * dimR];
3797: work = &invJ[Nc * dimR];
3798: for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.;
3799: for (j = 0; j < numPoints; j++) {
3800: PetscReal normPoint = DMPlex_NormD_Internal(Nc, &realCoords[j * Nc]);
3801: normPoint = normPoint > PETSC_SMALL ? normPoint : 1.0;
3802: for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
3803: PetscReal *guess = &refCoords[j * dimR], error = 0;
3804: PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL));
3805: for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k];
3806: for (k = 0; k < Nc * dimR; k++) J[k] = 0.;
3807: for (k = 0; k < pdim; k++) {
3808: for (l = 0; l < Nc; l++) {
3809: resNeg[l] -= modes[k] * B[k * Nc + l];
3810: for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
3811: }
3812: }
3813: if (0 && PetscDefined(USE_DEBUG)) {
3814: PetscReal maxAbs = 0.;
3816: for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3817: PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3818: }
3819: error = DMPlex_NormD_Internal(Nc, resNeg);
3820: if (error < tolerance * normPoint) {
3821: if (tol) *tol = error / normPoint;
3822: break;
3823: }
3824: PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess));
3825: }
3826: }
3827: PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3828: PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3829: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3830: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3831: PetscFunctionReturn(PETSC_SUCCESS);
3832: }
3834: /* TODO: TOBY please fix this for Nc > 1 */
3835: PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3836: {
3837: PetscInt numComp, pdim, i, j, k, l, coordSize;
3838: PetscScalar *nodes = NULL;
3839: PetscReal *invV, *modes;
3840: PetscReal *B;
3842: PetscFunctionBegin;
3843: PetscCall(PetscFEGetDimension(fe, &pdim));
3844: PetscCall(PetscFEGetNumComponents(fe, &numComp));
3845: PetscCheck(numComp == Nc, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "coordinate discretization must have as many components (%" PetscInt_FMT ") as embedding dimension (!= %" PetscInt_FMT ")", numComp, Nc);
3846: /* we shouldn't apply inverse closure permutation, if one exists */
3847: PetscCall(DMPlexVecGetOrientedClosure(dm, NULL, PETSC_FALSE, coords, cell, 0, &coordSize, &nodes));
3848: /* convert nodes to values in the stable evaluation basis */
3849: PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3850: invV = fe->invV;
3851: for (i = 0; i < pdim; ++i) {
3852: modes[i] = 0.;
3853: for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3854: }
3855: PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3856: PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL));
3857: for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3858: for (j = 0; j < numPoints; j++) {
3859: PetscReal *mapped = &realCoords[j * Nc];
3861: for (k = 0; k < pdim; k++) {
3862: for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3863: }
3864: }
3865: PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3866: PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3867: PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3868: PetscFunctionReturn(PETSC_SUCCESS);
3869: }
3871: /*@
3872: DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element
3873: using a single element map.
3875: Not Collective
3877: Input Parameters:
3878: + dm - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3879: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3880: as a multilinear map for tensor-product elements
3881: . cell - the cell whose map is used.
3882: . numPoints - the number of points to locate
3883: - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3885: Output Parameter:
3886: . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`)
3888: Level: intermediate
3890: Notes:
3891: This inversion will be accurate inside the reference element, but may be inaccurate for
3892: mappings that do not extend uniquely outside the reference cell (e.g, most non-affine maps)
3894: .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()`
3895: @*/
3896: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3897: {
3898: PetscInt dimC, dimR, depth, i, cellHeight, height;
3899: DMPolytopeType ct;
3900: DM coordDM = NULL;
3901: Vec coords;
3902: PetscFE fe = NULL;
3904: PetscFunctionBegin;
3906: PetscCall(DMGetDimension(dm, &dimR));
3907: PetscCall(DMGetCoordinateDim(dm, &dimC));
3908: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3909: PetscCall(DMPlexGetDepth(dm, &depth));
3910: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3911: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3912: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
3913: if (coordDM) {
3914: PetscInt coordFields;
3916: PetscCall(DMGetNumFields(coordDM, &coordFields));
3917: if (coordFields) {
3918: PetscClassId id;
3919: PetscObject disc;
3921: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3922: PetscCall(PetscObjectGetClassId(disc, &id));
3923: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3924: }
3925: }
3926: PetscCall(DMPlexGetCellType(dm, cell, &ct));
3927: PetscCall(DMPlexGetPointHeight(dm, cell, &height));
3928: PetscCheck(height == cellHeight, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in a cell, height = %" PetscInt_FMT, cell, height);
3929: PetscCheck(!DMPolytopeTypeIsHybrid(ct) && ct != DM_POLYTOPE_FV_GHOST, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " is unsupported cell type %s", cell, DMPolytopeTypes[ct]);
3930: if (!fe) { /* implicit discretization: affine or multilinear */
3931: PetscInt coneSize;
3932: PetscBool isSimplex, isTensor;
3934: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3935: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3936: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3937: if (isSimplex) {
3938: PetscReal detJ, *v0, *J, *invJ;
3940: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3941: J = &v0[dimC];
3942: invJ = &J[dimC * dimC];
3943: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ));
3944: for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3945: const PetscReal x0[3] = {-1., -1., -1.};
3947: CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3948: }
3949: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3950: } else if (isTensor) {
3951: PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3952: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3953: } else {
3954: PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR, 7, NULL));
3955: }
3956: PetscFunctionReturn(PETSC_SUCCESS);
3957: }
3959: /*@
3960: DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the mesh for a single element map.
3962: Not Collective
3964: Input Parameters:
3965: + dm - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3966: implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3967: as a multilinear map for tensor-product elements
3968: . cell - the cell whose map is used.
3969: . numPoints - the number of points to locate
3970: - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`)
3972: Output Parameter:
3973: . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)
3975: Level: intermediate
3977: .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()`
3978: @*/
3979: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3980: {
3981: PetscInt dimC, dimR, depth, i, cellHeight, height;
3982: DMPolytopeType ct;
3983: DM coordDM = NULL;
3984: Vec coords;
3985: PetscFE fe = NULL;
3987: PetscFunctionBegin;
3989: PetscCall(DMGetDimension(dm, &dimR));
3990: PetscCall(DMGetCoordinateDim(dm, &dimC));
3991: if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3992: PetscCall(DMPlexGetDepth(dm, &depth));
3993: PetscCall(DMGetCoordinatesLocal(dm, &coords));
3994: PetscCall(DMGetCoordinateDM(dm, &coordDM));
3995: PetscCall(DMPlexGetVTKCellHeight(dm, &cellHeight));
3996: if (coordDM) {
3997: PetscInt coordFields;
3999: PetscCall(DMGetNumFields(coordDM, &coordFields));
4000: if (coordFields) {
4001: PetscClassId id;
4002: PetscObject disc;
4004: PetscCall(DMGetField(coordDM, 0, NULL, &disc));
4005: PetscCall(PetscObjectGetClassId(disc, &id));
4006: if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
4007: }
4008: }
4009: PetscCall(DMPlexGetCellType(dm, cell, &ct));
4010: PetscCall(DMPlexGetPointHeight(dm, cell, &height));
4011: PetscCheck(height == cellHeight, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in a cell, height = %" PetscInt_FMT, cell, height);
4012: PetscCheck(!DMPolytopeTypeIsHybrid(ct) && ct != DM_POLYTOPE_FV_GHOST, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " is unsupported cell type %s", cell, DMPolytopeTypes[ct]);
4013: if (!fe) { /* implicit discretization: affine or multilinear */
4014: PetscInt coneSize;
4015: PetscBool isSimplex, isTensor;
4017: PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
4018: isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
4019: isTensor = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
4020: if (isSimplex) {
4021: PetscReal detJ, *v0, *J;
4023: PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
4024: J = &v0[dimC];
4025: PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ));
4026: for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
4027: const PetscReal xi0[3] = {-1., -1., -1.};
4029: CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
4030: }
4031: PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
4032: } else if (isTensor) {
4033: PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
4034: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
4035: } else {
4036: PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
4037: }
4038: PetscFunctionReturn(PETSC_SUCCESS);
4039: }
4041: void coordMap_identity(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
4042: {
4043: const PetscInt Nc = uOff[1] - uOff[0];
4044: PetscInt c;
4046: for (c = 0; c < Nc; ++c) f0[c] = u[c];
4047: }
4049: /* Constants are
4050: center location
4051: axis vector
4052: rotation angle
4053: */
4054: void coordMap_rotate(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
4055: {
4056: const PetscInt Nc = PetscMin(uOff[1] - uOff[0], 3);
4057: const PetscScalar *center = constants;
4058: const PetscScalar *k = &constants[Nc]; // The rotation axis k
4059: const PetscReal theta = PetscRealPart(constants[Nc * 2]);
4060: const PetscReal ct = PetscCosReal(theta);
4061: const PetscReal st = PetscSinReal(theta);
4062: PetscReal v[3];
4064: // Translate to coordinate with center at the origin
4065: for (PetscInt d = 0; d < Nc; ++d) v[d] = PetscRealPart(u[d] - center[d]);
4066: switch (dim) {
4067: case 2:
4068: /* Rotate point, axis is ignored in 2D
4069: / ct -st \
4070: \ st ct / */
4071: f0[0] = ct * v[0] - st * v[1];
4072: f0[1] = st * v[0] + ct * v[1];
4073: break;
4074: case 3:
4075: /* / ct + k_x^2 (1 - ct) & k_x k_y (1 - ct) - k_z st & k_x k_z (1 - ct) + k_y st \
4076: | k_y k_x (1 - ct) + k_z st & ct + k_y^2 (1 - ct) & k_y k_z (1 - ct) - k_x st |
4077: \ k_z k_x (1 - ct) - k_y st & k_z k_y (1 - ct) + k_x st & ct + k_z^2 (1 - ct) / */
4078: f0[0] = (ct + k[0] * k[0] * (1. - ct)) * v[0] + (k[0] * k[1] * (1. - ct) - k[2] * st) * v[1] + (k[0] * k[2] * (1. - ct) + k[1] * st) * v[2];
4079: f0[1] = (k[1] * k[0] * (1. - ct) + k[2] * st) * v[0] + (ct + k[1] * k[1] * (1. - ct)) * v[1] + (k[1] * k[2] * (1. - ct) - k[0] * st) * v[2];
4080: f0[2] = (k[2] * k[0] * (1. - ct) - k[1] * st) * v[0] + (k[2] * k[1] * (1. - ct) + k[0] * st) * v[1] + (ct + k[2] * k[2] * (1. - ct)) * v[2];
4081: break;
4082: default:
4083: for (PetscInt d = 0; d < Nc; ++d) f0[d] = v[d];
4084: }
4085: // Translate back to original coordinates
4086: for (PetscInt d = 0; d < Nc; ++d) f0[d] += center[d];
4087: }
4089: /* Shear applies the transformation, assuming we fix z,
4090: / 1 0 m_0 \
4091: | 0 1 m_1 |
4092: \ 0 0 1 /
4093: */
4094: void coordMap_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
4095: {
4096: const PetscInt Nc = uOff[1] - uOff[0];
4097: const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
4098: PetscInt c;
4100: for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax];
4101: }
4103: /* Flare applies the transformation, assuming we fix x_f,
4105: x_i = x_i * alpha_i x_f
4106: */
4107: void coordMap_flare(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
4108: {
4109: const PetscInt Nc = uOff[1] - uOff[0];
4110: const PetscInt cf = (PetscInt)PetscRealPart(constants[0]);
4112: for (PetscInt c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]);
4113: }
4115: /*
4116: We would like to map the unit square to a quarter of the annulus between circles of radius 1 and 2. We start by mapping the straight sections, which
4117: will correspond to the top and bottom of our square. So
4119: (0,0)--(1,0) ==> (1,0)--(2,0) Just a shift of (1,0)
4120: (0,1)--(1,1) ==> (0,1)--(0,2) Switch x and y
4122: So it looks like we want to map each layer in y to a ray, so x is the radius and y is the angle:
4124: (x, y) ==> (x+1, \pi/2 y) in (r', \theta') space
4125: ==> ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
4126: */
4127: void coordMap_annulus(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4128: {
4129: const PetscReal ri = PetscRealPart(constants[0]);
4130: const PetscReal ro = PetscRealPart(constants[1]);
4132: xp[0] = (x[0] * (ro - ri) + ri) * PetscCosReal(0.5 * PETSC_PI * x[1]);
4133: xp[1] = (x[0] * (ro - ri) + ri) * PetscSinReal(0.5 * PETSC_PI * x[1]);
4134: }
4136: /*
4137: We would like to map the unit cube to a hemisphere of the spherical shell between balls of radius 1 and 2. We want to map the bottom surface onto the
4138: lower hemisphere and the upper surface onto the top, letting z be the radius.
4140: (x, y) ==> ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x) in (r', \theta', \phi') space
4141: ==> ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
4142: */
4143: void coordMap_shell(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4144: {
4145: const PetscReal pi4 = PETSC_PI / 4.0;
4146: const PetscReal ri = PetscRealPart(constants[0]);
4147: const PetscReal ro = PetscRealPart(constants[1]);
4148: const PetscReal rp = (x[2] + 1) * 0.5 * (ro - ri) + ri;
4149: const PetscReal phip = PetscAtan2Real(x[1], x[0]);
4150: const PetscReal thetap = 0.5 * PETSC_PI * (1.0 - ((((phip <= pi4) && (phip >= -pi4)) || ((phip >= 3.0 * pi4) || (phip <= -3.0 * pi4))) ? PetscAbsReal(x[0]) : PetscAbsReal(x[1])));
4152: xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
4153: xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
4154: xp[2] = rp * PetscSinReal(thetap);
4155: }
4157: void coordMap_sinusoid(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4158: {
4159: const PetscReal c = PetscRealPart(constants[0]);
4160: const PetscReal m = PetscRealPart(constants[1]);
4161: const PetscReal n = PetscRealPart(constants[2]);
4163: xp[0] = x[0];
4164: xp[1] = x[1];
4165: if (dim > 2) xp[2] = c * PetscCosReal(2. * m * PETSC_PI * x[0]) * PetscCosReal(2. * n * PETSC_PI * x[1]);
4166: }
4168: /* This function maps the cylinder [0, r] x [0, 1] along z to the torus revolved around z with radius R,
4169: x' = (R + y) cos(2 pi z)
4170: y' = (R + y) sin(2 pi z)
4171: z' = x
4172: */
4173: void coordMap_torus(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xp[])
4174: {
4175: const PetscReal R = PetscRealPart(constants[0]);
4177: xp[0] = (R + x[1]) * PetscCosReal(2 * PETSC_PI * x[2]);
4178: xp[1] = (R + x[1]) * PetscSinReal(2 * PETSC_PI * x[2]);
4179: xp[2] = x[0];
4180: }
4182: /*@C
4183: DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates.
4185: Not Collective
4187: Input Parameters:
4188: + dm - The `DM`
4189: . time - The time
4190: - func - The function transforming current coordinates to new coordinates
4192: Calling sequence of `func`:
4193: + dim - The spatial dimension
4194: . Nf - The number of input fields (here 1)
4195: . NfAux - The number of input auxiliary fields
4196: . uOff - The offset of the coordinates in u[] (here 0)
4197: . uOff_x - The offset of the coordinates in u_x[] (here 0)
4198: . u - The coordinate values at this point in space
4199: . u_t - The coordinate time derivative at this point in space (here `NULL`)
4200: . u_x - The coordinate derivatives at this point in space
4201: . aOff - The offset of each auxiliary field in u[]
4202: . aOff_x - The offset of each auxiliary field in u_x[]
4203: . a - The auxiliary field values at this point in space
4204: . a_t - The auxiliary field time derivative at this point in space (or `NULL`)
4205: . a_x - The auxiliary field derivatives at this point in space
4206: . t - The current time
4207: . x - The coordinates of this point (here not used)
4208: . numConstants - The number of constants
4209: . constants - The value of each constant
4210: - f - The new coordinates at this point in space
4212: Level: intermediate
4214: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
4215: @*/
4216: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, void (*func)(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]))
4217: {
4218: DM cdm;
4219: PetscDS cds;
4220: DMField cf;
4221: PetscObject obj;
4222: PetscClassId id;
4223: Vec lCoords, tmpCoords;
4225: PetscFunctionBegin;
4226: if (!func) PetscCall(DMPlexGetCoordinateMap(dm, &func));
4227: PetscCall(DMGetCoordinateDM(dm, &cdm));
4228: PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
4229: PetscCall(DMGetDS(cdm, &cds));
4230: PetscCall(PetscDSGetDiscretization(cds, 0, &obj));
4231: PetscCall(PetscObjectGetClassId(obj, &id));
4232: if (id != PETSCFE_CLASSID) {
4233: PetscSection cSection;
4234: const PetscScalar *constants;
4235: PetscScalar *coords, f[16];
4236: PetscInt dim, cdim, Nc, vStart, vEnd;
4238: PetscCall(DMGetDimension(dm, &dim));
4239: PetscCall(DMGetCoordinateDim(dm, &cdim));
4240: PetscCheck(cdim <= 16, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Affine version of DMPlexRemapGeometry is currently limited to dimensions <= 16, not %" PetscInt_FMT, cdim);
4241: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
4242: PetscCall(DMGetCoordinateSection(dm, &cSection));
4243: PetscCall(PetscDSGetConstants(cds, &Nc, &constants));
4244: PetscCall(VecGetArrayWrite(lCoords, &coords));
4245: for (PetscInt v = vStart; v < vEnd; ++v) {
4246: PetscInt uOff[2] = {0, cdim};
4247: PetscInt off;
4249: PetscCall(PetscSectionGetOffset(cSection, v, &off));
4250: (*func)(dim, 1, 0, uOff, NULL, &coords[off], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, Nc, constants, f);
4251: for (PetscInt c = 0; c < cdim; ++c) coords[off + c] = f[c];
4252: }
4253: PetscCall(VecRestoreArrayWrite(lCoords, &coords));
4254: } else {
4255: PetscCall(DMGetLocalVector(cdm, &tmpCoords));
4256: PetscCall(VecCopy(lCoords, tmpCoords));
4257: /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
4258: PetscCall(DMGetCoordinateField(dm, &cf));
4259: cdm->coordinates[0].field = cf;
4260: PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords));
4261: cdm->coordinates[0].field = NULL;
4262: PetscCall(DMRestoreLocalVector(cdm, &tmpCoords));
4263: PetscCall(DMSetCoordinatesLocal(dm, lCoords));
4264: }
4265: PetscFunctionReturn(PETSC_SUCCESS);
4266: }
4268: /*@
4269: DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.
4271: Not Collective
4273: Input Parameters:
4274: + dm - The `DMPLEX`
4275: . direction - The shear coordinate direction, e.g. `DM_X` is the x-axis
4276: - multipliers - The multiplier m for each direction which is not the shear direction
4278: Level: intermediate
4280: .seealso: `DMPLEX`, `DMPlexRemapGeometry()`, `DMDirection`, `DM_X`, `DM_Y`, `DM_Z`
4281: @*/
4282: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
4283: {
4284: DM cdm;
4285: PetscDS cds;
4286: PetscScalar *moduli;
4287: const PetscInt dir = (PetscInt)direction;
4288: PetscInt dE, d, e;
4290: PetscFunctionBegin;
4291: PetscCall(DMGetCoordinateDM(dm, &cdm));
4292: PetscCall(DMGetCoordinateDim(dm, &dE));
4293: PetscCall(PetscMalloc1(dE + 1, &moduli));
4294: moduli[0] = dir;
4295: for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
4296: PetscCall(DMGetDS(cdm, &cds));
4297: PetscCall(PetscDSSetConstants(cds, dE + 1, moduli));
4298: PetscCall(DMPlexRemapGeometry(dm, 0.0, coordMap_shear));
4299: PetscCall(PetscFree(moduli));
4300: PetscFunctionReturn(PETSC_SUCCESS);
4301: }