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", "shear", "flare", "annulus", "shell", "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 defined(PETSC_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 PetscInt  embedDim = 2;
440:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
441:   PetscReal       x        = PetscRealPart(point[0]);
442:   PetscReal       y        = PetscRealPart(point[1]);
443:   PetscReal       v0[2], J[4], invJ[4], detJ;
444:   PetscReal       xi, eta;

446:   PetscFunctionBegin;
447:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
448:   xi  = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
449:   eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);

451:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 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_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, f;
490:   PetscBool          isDG;

492:   PetscFunctionBegin;
493:   PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
494:   PetscCheck(numCoords == 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
495:   for (f = 0; f < 4; ++f) {
496:     PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 0]);
497:     PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 1]);
498:     PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 0]);
499:     PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 1]);

501:     if ((x == x_j) && (y == y_j)) {
502:       // point is a corner
503:       crossings = 1;
504:       break;
505:     }
506:     if ((y_j > y) != (y_i > y)) {
507:       PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j);
508:       if (slope == 0) {
509:         // point is a corner
510:         crossings = 1;
511:         break;
512:       }
513:       if ((slope < 0) != (y_i < y_j)) ++crossings;
514:     }
515:   }
516:   if (crossings % 2) *cell = c;
517:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
518:   PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
523: {
524:   const PetscInt  embedDim = 3;
525:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
526:   PetscReal       v0[3], J[9], invJ[9], detJ;
527:   PetscReal       x = PetscRealPart(point[0]);
528:   PetscReal       y = PetscRealPart(point[1]);
529:   PetscReal       z = PetscRealPart(point[2]);
530:   PetscReal       xi, eta, zeta;

532:   PetscFunctionBegin;
533:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
534:   xi   = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
535:   eta  = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
536:   zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);

538:   if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
539:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
540:   PetscFunctionReturn(PETSC_SUCCESS);
541: }

543: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
544: {
545:   const PetscScalar *array;
546:   PetscScalar       *coords    = NULL;
547:   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};
548:   PetscBool          found     = PETSC_TRUE;
549:   PetscInt           numCoords, f;
550:   PetscBool          isDG;

552:   PetscFunctionBegin;
553:   PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
554:   PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
555:   for (f = 0; f < 6; ++f) {
556:     /* Check the point is under plane */
557:     /*   Get face normal */
558:     PetscReal v_i[3];
559:     PetscReal v_j[3];
560:     PetscReal normal[3];
561:     PetscReal pp[3];
562:     PetscReal dot;

564:     v_i[0]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
565:     v_i[1]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
566:     v_i[2]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
567:     v_j[0]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
568:     v_j[1]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
569:     v_j[2]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
570:     normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
571:     normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
572:     normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
573:     pp[0]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
574:     pp[1]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
575:     pp[2]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
576:     dot       = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];

578:     /* Check that projected point is in face (2D location problem) */
579:     if (dot < 0.0) {
580:       found = PETSC_FALSE;
581:       break;
582:     }
583:   }
584:   if (found) *cell = c;
585:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
586:   PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
587:   PetscFunctionReturn(PETSC_SUCCESS);
588: }

590: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
591: {
592:   PetscInt d;

594:   PetscFunctionBegin;
595:   box->dim = dim;
596:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.;
597:   PetscFunctionReturn(PETSC_SUCCESS);
598: }

600: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
601: {
602:   PetscFunctionBegin;
603:   PetscCall(PetscCalloc1(1, box));
604:   PetscCall(PetscGridHashInitialize_Internal(*box, dim, point));
605:   PetscFunctionReturn(PETSC_SUCCESS);
606: }

608: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
609: {
610:   PetscInt d;

612:   PetscFunctionBegin;
613:   for (d = 0; d < box->dim; ++d) {
614:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
615:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
616:   }
617:   PetscFunctionReturn(PETSC_SUCCESS);
618: }

620: static PetscErrorCode DMPlexCreateGridHash(DM dm, PetscGridHash *box)
621: {
622:   Vec                coordinates;
623:   const PetscScalar *coords;
624:   PetscInt           cdim, N, bs;

626:   PetscFunctionBegin;
627:   PetscCall(DMGetCoordinateDim(dm, &cdim));
628:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
629:   PetscCall(VecGetArrayRead(coordinates, &coords));
630:   PetscCall(VecGetLocalSize(coordinates, &N));
631:   PetscCall(VecGetBlockSize(coordinates, &bs));
632:   PetscCheck(bs == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Coordinate block size %" PetscInt_FMT " != %" PetscInt_FMT " coordinate dimension", bs, cdim);

634:   PetscCall(PetscGridHashCreate(PetscObjectComm((PetscObject)dm), cdim, coords, box));
635:   for (PetscInt i = 0; i < N; i += cdim) PetscCall(PetscGridHashEnlarge(*box, &coords[i]));

637:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

641: /*@C
642:   PetscGridHashSetGrid - Divide the grid into boxes

644:   Not Collective

646:   Input Parameters:
647: + box - The grid hash object
648: . n   - The number of boxes in each dimension, may use `PETSC_DETERMINE` for the entries
649: - h   - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE`, if not needed you can pass in `NULL`

651:   Level: developer

653: .seealso: `DMPLEX`, `PetscGridHashCreate()`
654: @*/
655: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
656: {
657:   PetscInt d;

659:   PetscFunctionBegin;
660:   PetscAssertPointer(n, 2);
661:   if (h) PetscAssertPointer(h, 3);
662:   for (d = 0; d < box->dim; ++d) {
663:     box->extent[d] = box->upper[d] - box->lower[d];
664:     if (n[d] == PETSC_DETERMINE) {
665:       PetscCheck(h, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Missing h");
666:       box->h[d] = h[d];
667:       box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
668:     } else {
669:       box->n[d] = n[d];
670:       box->h[d] = box->extent[d] / n[d];
671:     }
672:   }
673:   PetscFunctionReturn(PETSC_SUCCESS);
674: }

676: /*@C
677:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

679:   Not Collective

681:   Input Parameters:
682: + box       - The grid hash object
683: . numPoints - The number of input points
684: - points    - The input point coordinates

686:   Output Parameters:
687: + dboxes - An array of `numPoints` x `dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
688: - boxes  - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`

690:   Level: developer

692:   Note:
693:   This only guarantees that a box contains a point, not that a cell does.

695: .seealso: `DMPLEX`, `PetscGridHashCreate()`
696: @*/
697: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
698: {
699:   const PetscReal *lower = box->lower;
700:   const PetscReal *upper = box->upper;
701:   const PetscReal *h     = box->h;
702:   const PetscInt  *n     = box->n;
703:   const PetscInt   dim   = box->dim;
704:   PetscInt         d, p;

706:   PetscFunctionBegin;
707:   for (p = 0; p < numPoints; ++p) {
708:     for (d = 0; d < dim; ++d) {
709:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

711:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
712:       if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
713:       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", 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);
714:       dboxes[p * dim + d] = dbox;
715:     }
716:     if (boxes)
717:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
718:   }
719:   PetscFunctionReturn(PETSC_SUCCESS);
720: }

722: /*
723:   PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

725:   Not Collective

727:   Input Parameters:
728: + box         - The grid hash object
729: . cellSection - The PetscSection mapping cells to boxes
730: . numPoints   - The number of input points
731: - points      - The input point coordinates

733:   Output Parameters:
734: + dboxes - An array of `numPoints`*`dim` integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
735: . boxes  - An array of `numPoints` integers expressing the enclosing box as single number, or `NULL`
736: - found  - Flag indicating if point was located within a box

738:   Level: developer

740:   Note:
741:   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.

743: .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()`
744: */
745: static PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
746: {
747:   const PetscReal *lower = box->lower;
748:   const PetscReal *upper = box->upper;
749:   const PetscReal *h     = box->h;
750:   const PetscInt  *n     = box->n;
751:   const PetscInt   dim   = box->dim;
752:   PetscInt         bStart, bEnd, d, p;

754:   PetscFunctionBegin;
756:   *found = PETSC_FALSE;
757:   PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd));
758:   for (p = 0; p < numPoints; ++p) {
759:     for (d = 0; d < dim; ++d) {
760:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

762:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
763:       if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS);
764:       dboxes[p * dim + d] = dbox;
765:     }
766:     if (boxes)
767:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
768:     // It is possible for a box to overlap no grid cells
769:     if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS);
770:   }
771:   *found = PETSC_TRUE;
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

775: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
776: {
777:   PetscFunctionBegin;
778:   if (*box) {
779:     PetscCall(PetscSectionDestroy(&(*box)->cellSection));
780:     PetscCall(ISDestroy(&(*box)->cells));
781:     PetscCall(DMLabelDestroy(&(*box)->cellsSparse));
782:   }
783:   PetscCall(PetscFree(*box));
784:   PetscFunctionReturn(PETSC_SUCCESS);
785: }

787: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
788: {
789:   DMPolytopeType ct;

791:   PetscFunctionBegin;
792:   PetscCall(DMPlexGetCellType(dm, cellStart, &ct));
793:   switch (ct) {
794:   case DM_POLYTOPE_SEGMENT:
795:     PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell));
796:     break;
797:   case DM_POLYTOPE_TRIANGLE:
798:     PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell));
799:     break;
800:   case DM_POLYTOPE_QUADRILATERAL:
801:     PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell));
802:     break;
803:   case DM_POLYTOPE_TETRAHEDRON:
804:     PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell));
805:     break;
806:   case DM_POLYTOPE_HEXAHEDRON:
807:     PetscCall(DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell));
808:     break;
809:   default:
810:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
811:   }
812:   PetscFunctionReturn(PETSC_SUCCESS);
813: }

815: /*
816:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
817: */
818: static PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
819: {
820:   DMPolytopeType ct;

822:   PetscFunctionBegin;
823:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
824:   switch (ct) {
825:   case DM_POLYTOPE_TRIANGLE:
826:     PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint));
827:     break;
828: #if 0
829:     case DM_POLYTOPE_QUADRILATERAL:
830:     PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break;
831:     case DM_POLYTOPE_TETRAHEDRON:
832:     PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break;
833:     case DM_POLYTOPE_HEXAHEDRON:
834:     PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break;
835: #endif
836:   default:
837:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
838:   }
839:   PetscFunctionReturn(PETSC_SUCCESS);
840: }

842: /*
843:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX`

845:   Collective

847:   Input Parameter:
848: . dm - The `DMPLEX`

850:   Output Parameter:
851: . localBox - The grid hash object

853:   Level: developer

855:   Notes:
856:   How do we determine all boxes intersecting a given cell?

858:   1) Get convex body enclosing cell. We will use a box called the box-hull.

860:   2) Get smallest brick of boxes enclosing the box-hull

862:   3) Each box is composed of 6 planes, 3 lower and 3 upper. We loop over dimensions, and
863:      for each new plane determine whether the cell is on the negative side, positive side, or intersects it.

865:      a) If the cell is on the negative side of the lower planes, it is not in the box

867:      b) If the cell is on the positive side of the upper planes, it is not in the box

869:      c) If there is no intersection, it is in the box

871:      d) If any intersection point is within the box limits, it is in the box

873: .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
874: */
875: static PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
876: {
877:   PetscInt        debug = ((DM_Plex *)dm->data)->printLocate;
878:   PetscGridHash   lbox;
879:   PetscSF         sf;
880:   const PetscInt *leaves;
881:   PetscInt       *dboxes, *boxes;
882:   PetscInt        cdim, cStart, cEnd, Nl = -1;
883:   PetscBool       flg;

885:   PetscFunctionBegin;
886:   PetscCall(DMGetCoordinateDim(dm, &cdim));
887:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
888:   PetscCall(DMPlexCreateGridHash(dm, &lbox));
889:   {
890:     PetscInt n[3], d;

892:     PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &d, &flg));
893:     if (flg) {
894:       for (PetscInt i = d; i < cdim; ++i) n[i] = n[d - 1];
895:     } else {
896:       for (PetscInt i = 0; i < cdim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / cdim) * 0.8));
897:     }
898:     PetscCall(PetscGridHashSetGrid(lbox, n, NULL));
899:     if (debug)
900:       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.,
901:                             (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.));
902:   }

904:   PetscCall(DMGetPointSF(dm, &sf));
905:   if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
906:   Nl = PetscMax(Nl, 0);
907:   PetscCall(PetscCalloc2(16 * cdim, &dboxes, 16, &boxes));

909:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse));
910:   PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd));
911:   for (PetscInt c = cStart; c < cEnd; ++c) {
912:     PetscReal          intPoints[6 * 6 * 6 * 3];
913:     const PetscScalar *array;
914:     PetscScalar       *coords            = NULL;
915:     const PetscReal   *h                 = lbox->h;
916:     PetscReal          normal[9]         = {1., 0., 0., 0., 1., 0., 0., 0., 1.};
917:     PetscReal         *lowerIntPoints[3] = {&intPoints[0 * 6 * 6 * 3], &intPoints[1 * 6 * 6 * 3], &intPoints[2 * 6 * 6 * 3]};
918:     PetscReal         *upperIntPoints[3] = {&intPoints[3 * 6 * 6 * 3], &intPoints[4 * 6 * 6 * 3], &intPoints[5 * 6 * 6 * 3]};
919:     PetscReal          lp[3], up[3], *tmp;
920:     PetscInt           numCoords, idx, dlim[6], lowerInt[3], upperInt[3];
921:     PetscBool          isDG, lower[3], upper[3];

923:     PetscCall(PetscFindInt(c, Nl, leaves, &idx));
924:     if (idx >= 0) continue;
925:     // Get grid of boxes containing the cell
926:     PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
927:     PetscCall(PetscGridHashGetEnclosingBox(lbox, numCoords / cdim, coords, dboxes, boxes));
928:     PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
929:     for (PetscInt d = 0; d < cdim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
930:     for (PetscInt d = cdim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
931:     for (PetscInt e = 1; e < numCoords / cdim; ++e) {
932:       for (PetscInt d = 0; d < cdim; ++d) {
933:         dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * cdim + d]);
934:         dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * cdim + d]);
935:       }
936:     }
937:     if (debug > 4) {
938:       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]));
939:     }
940:     // Initialize with lower planes for first box
941:     for (PetscInt d = 0; d < cdim; ++d) {
942:       lp[d] = lbox->lower[d] + dlim[d * 2 + 0] * h[d];
943:       up[d] = lp[d] + h[d];
944:     }
945:     for (PetscInt d = 0; d < cdim; ++d) {
946:       PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, lp, &normal[d * 3], &lower[d], &lowerInt[d], lowerIntPoints[d]));
947:       if (debug > 4) {
948:         if (!lowerInt[d])
949:           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"));
950:         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]));
951:       }
952:     }
953:     // Loop over grid
954:     for (PetscInt k = dlim[2 * 2 + 0]; k <= dlim[2 * 2 + 1]; ++k, lp[2] = up[2], up[2] += h[2]) {
955:       if (cdim > 2) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 2], &upper[2], &upperInt[2], upperIntPoints[2]));
956:       if (cdim > 2 && debug > 4) {
957:         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"));
958:         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]));
959:       }
960:       for (PetscInt j = dlim[1 * 2 + 0]; j <= dlim[1 * 2 + 1]; ++j, lp[1] = up[1], up[1] += h[1]) {
961:         if (cdim > 1) PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 1], &upper[1], &upperInt[1], upperIntPoints[1]));
962:         if (cdim > 1 && debug > 4) {
963:           if (!upperInt[1])
964:             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"));
965:           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]));
966:         }
967:         for (PetscInt i = dlim[0 * 2 + 0]; i <= dlim[0 * 2 + 1]; ++i, lp[0] = up[0], up[0] += h[0]) {
968:           const PetscInt box    = (k * lbox->n[1] + j) * lbox->n[0] + i;
969:           PetscBool      excNeg = PETSC_TRUE;
970:           PetscBool      excPos = PETSC_TRUE;
971:           PetscInt       NlInt  = 0;
972:           PetscInt       NuInt  = 0;

974:           PetscCall(DMPlexGetPlaneCellIntersection_Internal(dm, c, up, &normal[3 * 0], &upper[0], &upperInt[0], upperIntPoints[0]));
975:           if (debug > 4) {
976:             if (!upperInt[0])
977:               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"));
978:             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]));
979:           }
980:           for (PetscInt d = 0; d < cdim; ++d) {
981:             NlInt += lowerInt[d];
982:             NuInt += upperInt[d];
983:           }
984:           // If there is no intersection...
985:           if (!NlInt && !NuInt) {
986:             // If the cell is on the negative side of the lower planes, it is not in the box
987:             for (PetscInt d = 0; d < cdim; ++d)
988:               if (lower[d]) {
989:                 excNeg = PETSC_FALSE;
990:                 break;
991:               }
992:             // If the cell is on the positive side of the upper planes, it is not in the box
993:             for (PetscInt d = 0; d < cdim; ++d)
994:               if (!upper[d]) {
995:                 excPos = PETSC_FALSE;
996:                 break;
997:               }
998:             if (excNeg || excPos) {
999:               if (debug && excNeg) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " is on the negative side of the lower plane\n", c));
1000:               if (debug && excPos) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " is on the positive side of the upper plane\n", c));
1001:               continue;
1002:             }
1003:             // Otherwise it is in the box
1004:             if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " is contained in box %" PetscInt_FMT "\n", c, box));
1005:             PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1006:             continue;
1007:           }
1008:           /*
1009:             If any intersection point is within the box limits, it is in the box
1010:             We need to have tolerances here since intersection point calculations can introduce errors
1011:             Initialize a count to track which planes have intersection outside the box.
1012:             if two adjacent planes have intersection points upper and lower all outside the box, look
1013:             first at if another plane has intersection points outside the box, if so, it is inside the cell
1014:             look next if no intersection points exist on the other planes, and check if the planes are on the
1015:             outside of the intersection points but on opposite ends. If so, the box cuts through the cell.
1016:           */
1017:           PetscInt outsideCount[6] = {0, 0, 0, 0, 0, 0};
1018:           for (PetscInt plane = 0; plane < cdim; ++plane) {
1019:             for (PetscInt ip = 0; ip < lowerInt[plane]; ++ip) {
1020:               PetscInt d;

1022:               for (d = 0; d < cdim; ++d) {
1023:                 if ((lowerIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (lowerIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1024:                   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
1025:                   break;
1026:                 }
1027:               }
1028:               if (d == cdim) {
1029:                 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " intersected lower plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1030:                 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1031:                 goto end;
1032:               }
1033:             }
1034:             for (PetscInt ip = 0; ip < upperInt[plane]; ++ip) {
1035:               PetscInt d;

1037:               for (d = 0; d < cdim; ++d) {
1038:                 if ((upperIntPoints[plane][ip * cdim + d] < (lp[d] - PETSC_SMALL)) || (upperIntPoints[plane][ip * cdim + d] > (up[d] + PETSC_SMALL))) {
1039:                   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
1040:                   break;
1041:                 }
1042:               }
1043:               if (d == cdim) {
1044:                 if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " intersected upper plane %" PetscInt_FMT " of box %" PetscInt_FMT "\n", c, plane, box));
1045:                 PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1046:                 goto end;
1047:               }
1048:             }
1049:           }
1050:           /*
1051:              Check the planes with intersections
1052:              in 2D, check if the square falls in the middle of a cell
1053:              ie all four planes have intersection points outside of the box
1054:              You do not want to be doing this, because it means your grid hashing is finer than your grid,
1055:              but we should still support it I guess
1056:           */
1057:           if (cdim == 2) {
1058:             PetscInt nIntersects = 0;
1059:             for (PetscInt d = 0; d < cdim; ++d) nIntersects += (outsideCount[d] + outsideCount[d + cdim]);
1060:             // if the count adds up to 8, that means each plane has 2 external intersections and thus it is in the cell
1061:             if (nIntersects == 8) {
1062:               PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1063:               goto end;
1064:             }
1065:           }
1066:           /*
1067:              In 3 dimensions, if two adjacent planes have at least 3 intersections outside the cell in the appropriate direction,
1068:              we then check the 3rd planar dimension. If a plane falls between intersection points, the cell belongs to that box.
1069:              If the planes are on opposite sides of the intersection points, the cell belongs to that box and it passes through the cell.
1070:           */
1071:           if (cdim == 3) {
1072:             PetscInt faces[3] = {0, 0, 0}, checkInternalFace = 0;
1073:             // Find two adjacent planes with at least 3 intersection points in the upper and lower
1074:             // if the third plane has 3 intersection points or more, a pyramid base is formed on that plane and it is in the cell
1075:             for (PetscInt d = 0; d < cdim; ++d)
1076:               if (outsideCount[d] >= 3 && outsideCount[cdim + d] >= 3) {
1077:                 faces[d]++;
1078:                 checkInternalFace++;
1079:               }
1080:             if (checkInternalFace == 3) {
1081:               // All planes have 3 intersection points, add it.
1082:               PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1083:               goto end;
1084:             }
1085:             // Gross, figure out which adjacent faces have at least 3 points
1086:             PetscInt nonIntersectingFace = -1;
1087:             if (faces[0] == faces[1]) nonIntersectingFace = 2;
1088:             if (faces[0] == faces[2]) nonIntersectingFace = 1;
1089:             if (faces[1] == faces[2]) nonIntersectingFace = 0;
1090:             if (nonIntersectingFace >= 0) {
1091:               for (PetscInt plane = 0; plane < cdim; ++plane) {
1092:                 if (!lowerInt[nonIntersectingFace] && !upperInt[nonIntersectingFace]) continue;
1093:                 // 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.
1094:                 for (PetscInt ip = 0; ip < lowerInt[nonIntersectingFace]; ++ip) {
1095:                   if (lowerIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || lowerIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1096:                 }
1097:                 for (PetscInt ip = 0; ip < upperInt[nonIntersectingFace]; ++ip) {
1098:                   if (upperIntPoints[plane][ip * cdim + nonIntersectingFace] > lp[nonIntersectingFace] - PETSC_SMALL || upperIntPoints[plane][ip * cdim + nonIntersectingFace] < up[nonIntersectingFace] + PETSC_SMALL) goto setpoint;
1099:                 }
1100:                 goto end;
1101:               }
1102:               // The points are within the bonds of the non intersecting planes, add it.
1103:             setpoint:
1104:               PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
1105:               goto end;
1106:             }
1107:           }
1108:         end:
1109:           lower[0]          = upper[0];
1110:           lowerInt[0]       = upperInt[0];
1111:           tmp               = lowerIntPoints[0];
1112:           lowerIntPoints[0] = upperIntPoints[0];
1113:           upperIntPoints[0] = tmp;
1114:         }
1115:         lp[0]             = lbox->lower[0] + dlim[0 * 2 + 0] * h[0];
1116:         up[0]             = lp[0] + h[0];
1117:         lower[1]          = upper[1];
1118:         lowerInt[1]       = upperInt[1];
1119:         tmp               = lowerIntPoints[1];
1120:         lowerIntPoints[1] = upperIntPoints[1];
1121:         upperIntPoints[1] = tmp;
1122:       }
1123:       lp[1]             = lbox->lower[1] + dlim[1 * 2 + 0] * h[1];
1124:       up[1]             = lp[1] + h[1];
1125:       lower[2]          = upper[2];
1126:       lowerInt[2]       = upperInt[2];
1127:       tmp               = lowerIntPoints[2];
1128:       lowerIntPoints[2] = upperIntPoints[2];
1129:       upperIntPoints[2] = tmp;
1130:     }
1131:   }
1132:   PetscCall(PetscFree2(dboxes, boxes));

1134:   if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF));
1135:   PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells));
1136:   PetscCall(DMLabelDestroy(&lbox->cellsSparse));
1137:   *localBox = lbox;
1138:   PetscFunctionReturn(PETSC_SUCCESS);
1139: }

1141: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
1142: {
1143:   PetscInt        debug = ((DM_Plex *)dm->data)->printLocate;
1144:   DM_Plex        *mesh  = (DM_Plex *)dm->data;
1145:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
1146:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
1147:   PetscInt        dim, Nl = 0, cStart, cEnd, numCells, c, d;
1148:   PetscSF         sf;
1149:   const PetscInt *leaves;
1150:   const PetscInt *boxCells;
1151:   PetscSFNode    *cells;
1152:   PetscScalar    *a;
1153:   PetscMPIInt     result;
1154:   PetscLogDouble  t0, t1;
1155:   PetscReal       gmin[3], gmax[3];
1156:   PetscInt        terminating_query_type[] = {0, 0, 0};
1157:   PetscMPIInt     rank;

1159:   PetscFunctionBegin;
1160:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
1161:   PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0));
1162:   PetscCall(PetscTime(&t0));
1163:   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.");
1164:   PetscCall(DMGetCoordinateDim(dm, &dim));
1165:   PetscCall(VecGetBlockSize(v, &bs));
1166:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result));
1167:   PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
1168:   PetscCheck(bs == dim, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %" PetscInt_FMT " must be the mesh coordinate dimension %" PetscInt_FMT, bs, dim);
1169:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
1170:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1171:   PetscCall(DMGetPointSF(dm, &sf));
1172:   if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
1173:   Nl = PetscMax(Nl, 0);
1174:   PetscCall(VecGetLocalSize(v, &numPoints));
1175:   PetscCall(VecGetArray(v, &a));
1176:   numPoints /= bs;
1177:   {
1178:     const PetscSFNode *sf_cells;

1180:     PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells));
1181:     if (sf_cells) {
1182:       PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n"));
1183:       cells = (PetscSFNode *)sf_cells;
1184:       reuse = PETSC_TRUE;
1185:     } else {
1186:       PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n"));
1187:       PetscCall(PetscMalloc1(numPoints, &cells));
1188:       /* initialize cells if created */
1189:       for (p = 0; p < numPoints; p++) {
1190:         cells[p].rank  = 0;
1191:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1192:       }
1193:     }
1194:   }
1195:   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1196:   if (hash) {
1197:     if (!mesh->lbox) {
1198:       PetscCall(PetscInfo(dm, "Initializing grid hashing\n"));
1199:       PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox));
1200:     }
1201:     /* Designate the local box for each point */
1202:     /* Send points to correct process */
1203:     /* Search cells that lie in each subbox */
1204:     /*   Should we bin points before doing search? */
1205:     PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells));
1206:   }
1207:   for (p = 0, numFound = 0; p < numPoints; ++p) {
1208:     const PetscScalar *point   = &a[p * bs];
1209:     PetscInt           dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
1210:     PetscBool          point_outside_domain = PETSC_FALSE;

1212:     /* check bounding box of domain */
1213:     for (d = 0; d < dim; d++) {
1214:       if (PetscRealPart(point[d]) < gmin[d]) {
1215:         point_outside_domain = PETSC_TRUE;
1216:         break;
1217:       }
1218:       if (PetscRealPart(point[d]) > gmax[d]) {
1219:         point_outside_domain = PETSC_TRUE;
1220:         break;
1221:       }
1222:     }
1223:     if (point_outside_domain) {
1224:       cells[p].rank  = 0;
1225:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1226:       terminating_query_type[0]++;
1227:       continue;
1228:     }

1230:     /* check initial values in cells[].index - abort early if found */
1231:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1232:       c              = cells[p].index;
1233:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
1234:       PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
1235:       if (cell >= 0) {
1236:         cells[p].rank  = 0;
1237:         cells[p].index = cell;
1238:         numFound++;
1239:       }
1240:     }
1241:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
1242:       terminating_query_type[1]++;
1243:       continue;
1244:     }

1246:     if (hash) {
1247:       PetscBool found_box;

1249:       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]), dim > 2 ? (double)PetscRealPart(point[2]) : 0.));
1250:       /* allow for case that point is outside box - abort early */
1251:       PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box));
1252:       if (found_box) {
1253:         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], dim > 2 ? dbin[2] : 0));
1254:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
1255:         PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1256:         PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1257:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
1258:           if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]    Checking for point in cell %" PetscInt_FMT "\n", rank, boxCells[c]));
1259:           PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell));
1260:           if (cell >= 0) {
1261:             if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d]      FOUND in cell %" PetscInt_FMT "\n", rank, cell));
1262:             cells[p].rank  = 0;
1263:             cells[p].index = cell;
1264:             numFound++;
1265:             terminating_query_type[2]++;
1266:             break;
1267:           }
1268:         }
1269:       }
1270:     } else {
1271:       for (c = cStart; c < cEnd; ++c) {
1272:         PetscInt idx;

1274:         PetscCall(PetscFindInt(c, Nl, leaves, &idx));
1275:         if (idx >= 0) continue;
1276:         PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
1277:         if (cell >= 0) {
1278:           cells[p].rank  = 0;
1279:           cells[p].index = cell;
1280:           numFound++;
1281:           terminating_query_type[2]++;
1282:           break;
1283:         }
1284:       }
1285:     }
1286:   }
1287:   if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells));
1288:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
1289:     for (p = 0; p < numPoints; p++) {
1290:       const PetscScalar *point     = &a[p * bs];
1291:       PetscReal          cpoint[3] = {0, 0, 0}, diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
1292:       PetscInt           dbin[3] = {-1, -1, -1}, bin, cellOffset, d, bestc = -1;

1294:       if (cells[p].index < 0) {
1295:         PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin));
1296:         PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
1297:         PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
1298:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
1299:           PetscCall(DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint));
1300:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
1301:           dist = DMPlex_NormD_Internal(dim, diff);
1302:           if (dist < distMax) {
1303:             for (d = 0; d < dim; ++d) best[d] = cpoint[d];
1304:             bestc   = boxCells[c];
1305:             distMax = dist;
1306:           }
1307:         }
1308:         if (distMax < PETSC_MAX_REAL) {
1309:           ++numFound;
1310:           cells[p].rank  = 0;
1311:           cells[p].index = bestc;
1312:           for (d = 0; d < dim; ++d) a[p * bs + d] = best[d];
1313:         }
1314:       }
1315:     }
1316:   }
1317:   /* This code is only be relevant when interfaced to parallel point location */
1318:   /* Check for highest numbered proc that claims a point (do we care?) */
1319:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
1320:     PetscCall(PetscMalloc1(numFound, &found));
1321:     for (p = 0, numFound = 0; p < numPoints; p++) {
1322:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
1323:         if (numFound < p) cells[numFound] = cells[p];
1324:         found[numFound++] = p;
1325:       }
1326:     }
1327:   }
1328:   PetscCall(VecRestoreArray(v, &a));
1329:   if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER));
1330:   PetscCall(PetscTime(&t1));
1331:   if (hash) {
1332:     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]));
1333:   } else {
1334:     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]));
1335:   }
1336:   PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] npoints %" PetscInt_FMT " : time(rank0) %1.2e (sec): points/sec %1.4e\n", numPoints, t1 - t0, (double)((double)numPoints / (t1 - t0))));
1337:   PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0));
1338:   PetscFunctionReturn(PETSC_SUCCESS);
1339: }

1341: /*@C
1342:   DMPlexComputeProjection2Dto1D - Rewrite coordinates to be the 1D projection of the 2D coordinates

1344:   Not Collective

1346:   Input/Output Parameter:
1347: . 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

1349:   Output Parameter:
1350: . R - The rotation which accomplishes the projection, array of size 4

1352:   Level: developer

1354: .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1355: @*/
1356: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
1357: {
1358:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
1359:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
1360:   const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;

1362:   PetscFunctionBegin;
1363:   R[0]      = c;
1364:   R[1]      = -s;
1365:   R[2]      = s;
1366:   R[3]      = c;
1367:   coords[0] = 0.0;
1368:   coords[1] = r;
1369:   PetscFunctionReturn(PETSC_SUCCESS);
1370: }

1372: /*@C
1373:   DMPlexComputeProjection3Dto1D - Rewrite coordinates to be the 1D projection of the 3D coordinates

1375:   Not Collective

1377:   Input/Output Parameter:
1378: . 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

1380:   Output Parameter:
1381: . R - The rotation which accomplishes the projection, an array of size 9

1383:   Level: developer

1385:   Note:
1386:   This uses the basis completion described by Frisvad {cite}`frisvad2012building`

1388: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1389: @*/
1390: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1391: {
1392:   PetscReal x    = PetscRealPart(coords[3] - coords[0]);
1393:   PetscReal y    = PetscRealPart(coords[4] - coords[1]);
1394:   PetscReal z    = PetscRealPart(coords[5] - coords[2]);
1395:   PetscReal r    = PetscSqrtReal(x * x + y * y + z * z);
1396:   PetscReal rinv = 1. / r;

1398:   PetscFunctionBegin;
1399:   x *= rinv;
1400:   y *= rinv;
1401:   z *= rinv;
1402:   if (x > 0.) {
1403:     PetscReal inv1pX = 1. / (1. + x);

1405:     R[0] = x;
1406:     R[1] = -y;
1407:     R[2] = -z;
1408:     R[3] = y;
1409:     R[4] = 1. - y * y * inv1pX;
1410:     R[5] = -y * z * inv1pX;
1411:     R[6] = z;
1412:     R[7] = -y * z * inv1pX;
1413:     R[8] = 1. - z * z * inv1pX;
1414:   } else {
1415:     PetscReal inv1mX = 1. / (1. - x);

1417:     R[0] = x;
1418:     R[1] = z;
1419:     R[2] = y;
1420:     R[3] = y;
1421:     R[4] = -y * z * inv1mX;
1422:     R[5] = 1. - y * y * inv1mX;
1423:     R[6] = z;
1424:     R[7] = 1. - z * z * inv1mX;
1425:     R[8] = -y * z * inv1mX;
1426:   }
1427:   coords[0] = 0.0;
1428:   coords[1] = r;
1429:   PetscFunctionReturn(PETSC_SUCCESS);
1430: }

1432: /*@
1433:   DMPlexComputeProjection3Dto2D - Rewrite coordinates of 3 or more coplanar 3D points to a common 2D basis for the
1434:   plane.  The normal is defined by positive orientation of the first 3 points.

1436:   Not Collective

1438:   Input Parameter:
1439: . coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)

1441:   Input/Output Parameter:
1442: . coords - The interlaced coordinates of each coplanar 3D point; on output the first
1443:            2*coordSize/3 entries contain interlaced 2D points, with the rest undefined

1445:   Output Parameter:
1446: . 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.

1448:   Level: developer

1450: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1451: @*/
1452: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1453: {
1454:   PetscReal      x1[3], x2[3], n[3], c[3], norm;
1455:   const PetscInt dim = 3;
1456:   PetscInt       d, p;

1458:   PetscFunctionBegin;
1459:   /* 0) Calculate normal vector */
1460:   for (d = 0; d < dim; ++d) {
1461:     x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1462:     x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1463:   }
1464:   // n = x1 \otimes x2
1465:   n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1466:   n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1467:   n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1468:   norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1469:   for (d = 0; d < dim; d++) n[d] /= norm;
1470:   norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1471:   for (d = 0; d < dim; d++) x1[d] /= norm;
1472:   // x2 = n \otimes x1
1473:   x2[0] = n[1] * x1[2] - n[2] * x1[1];
1474:   x2[1] = n[2] * x1[0] - n[0] * x1[2];
1475:   x2[2] = n[0] * x1[1] - n[1] * x1[0];
1476:   for (d = 0; d < dim; d++) {
1477:     R[d * dim + 0] = x1[d];
1478:     R[d * dim + 1] = x2[d];
1479:     R[d * dim + 2] = n[d];
1480:     c[d]           = PetscRealPart(coords[0 * dim + d]);
1481:   }
1482:   for (p = 0; p < coordSize / dim; p++) {
1483:     PetscReal y[3];
1484:     for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1485:     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];
1486:   }
1487:   PetscFunctionReturn(PETSC_SUCCESS);
1488: }

1490: PETSC_UNUSED static inline void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1491: {
1492:   /* Signed volume is 1/2 the determinant

1494:    |  1  1  1 |
1495:    | x0 x1 x2 |
1496:    | y0 y1 y2 |

1498:      but if x0,y0 is the origin, we have

1500:    | x1 x2 |
1501:    | y1 y2 |
1502:   */
1503:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1504:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1505:   PetscReal       M[4], detM;
1506:   M[0] = x1;
1507:   M[1] = x2;
1508:   M[2] = y1;
1509:   M[3] = y2;
1510:   DMPlex_Det2D_Internal(&detM, M);
1511:   *vol = 0.5 * detM;
1512:   (void)PetscLogFlops(5.0);
1513: }

1515: PETSC_UNUSED static inline void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1516: {
1517:   /* Signed volume is 1/6th of the determinant

1519:    |  1  1  1  1 |
1520:    | x0 x1 x2 x3 |
1521:    | y0 y1 y2 y3 |
1522:    | z0 z1 z2 z3 |

1524:      but if x0,y0,z0 is the origin, we have

1526:    | x1 x2 x3 |
1527:    | y1 y2 y3 |
1528:    | z1 z2 z3 |
1529:   */
1530:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1531:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1532:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1533:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1534:   PetscReal       M[9], detM;
1535:   M[0] = x1;
1536:   M[1] = x2;
1537:   M[2] = x3;
1538:   M[3] = y1;
1539:   M[4] = y2;
1540:   M[5] = y3;
1541:   M[6] = z1;
1542:   M[7] = z2;
1543:   M[8] = z3;
1544:   DMPlex_Det3D_Internal(&detM, M);
1545:   *vol = -onesixth * detM;
1546:   (void)PetscLogFlops(10.0);
1547: }

1549: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1550: {
1551:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1552:   DMPlex_Det3D_Internal(vol, coords);
1553:   *vol *= -onesixth;
1554: }

1556: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1557: {
1558:   PetscSection       coordSection;
1559:   Vec                coordinates;
1560:   const PetscScalar *coords;
1561:   PetscInt           dim, d, off;

1563:   PetscFunctionBegin;
1564:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1565:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
1566:   PetscCall(PetscSectionGetDof(coordSection, e, &dim));
1567:   if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
1568:   PetscCall(PetscSectionGetOffset(coordSection, e, &off));
1569:   PetscCall(VecGetArrayRead(coordinates, &coords));
1570:   if (v0) {
1571:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1572:   }
1573:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
1574:   *detJ = 1.;
1575:   if (J) {
1576:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1577:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1578:     if (invJ) {
1579:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1580:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1581:     }
1582:   }
1583:   PetscFunctionReturn(PETSC_SUCCESS);
1584: }

1586: /*@C
1587:   DMPlexGetCellCoordinates - Get coordinates for a cell, taking into account periodicity

1589:   Not Collective

1591:   Input Parameters:
1592: + dm   - The `DMPLEX`
1593: - cell - The cell number

1595:   Output Parameters:
1596: + isDG   - Using cellwise coordinates
1597: . Nc     - The number of coordinates
1598: . array  - The coordinate array
1599: - coords - The cell coordinates

1601:   Level: developer

1603: .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1604: @*/
1605: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1606: {
1607:   DM                 cdm;
1608:   Vec                coordinates;
1609:   PetscSection       cs;
1610:   const PetscScalar *ccoords;
1611:   PetscInt           pStart, pEnd;

1613:   PetscFunctionBeginHot;
1614:   *isDG   = PETSC_FALSE;
1615:   *Nc     = 0;
1616:   *array  = NULL;
1617:   *coords = NULL;
1618:   /* Check for cellwise coordinates */
1619:   PetscCall(DMGetCellCoordinateSection(dm, &cs));
1620:   if (!cs) goto cg;
1621:   /* Check that the cell exists in the cellwise section */
1622:   PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd));
1623:   if (cell < pStart || cell >= pEnd) goto cg;
1624:   /* Check for cellwise coordinates for this cell */
1625:   PetscCall(PetscSectionGetDof(cs, cell, Nc));
1626:   if (!*Nc) goto cg;
1627:   /* Check for cellwise coordinates */
1628:   PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates));
1629:   if (!coordinates) goto cg;
1630:   /* Get cellwise coordinates */
1631:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1632:   PetscCall(VecGetArrayRead(coordinates, array));
1633:   PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords));
1634:   PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1635:   PetscCall(PetscArraycpy(*coords, ccoords, *Nc));
1636:   PetscCall(VecRestoreArrayRead(coordinates, array));
1637:   *isDG = PETSC_TRUE;
1638:   PetscFunctionReturn(PETSC_SUCCESS);
1639: cg:
1640:   /* Use continuous coordinates */
1641:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1642:   PetscCall(DMGetCoordinateSection(dm, &cs));
1643:   PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1644:   PetscCall(DMPlexVecGetOrientedClosure_Internal(cdm, cs, PETSC_FALSE, coordinates, cell, 0, Nc, coords));
1645:   PetscFunctionReturn(PETSC_SUCCESS);
1646: }

1648: /*@C
1649:   DMPlexRestoreCellCoordinates - Get coordinates for a cell, taking into account periodicity

1651:   Not Collective

1653:   Input Parameters:
1654: + dm   - The `DMPLEX`
1655: - cell - The cell number

1657:   Output Parameters:
1658: + isDG   - Using cellwise coordinates
1659: . Nc     - The number of coordinates
1660: . array  - The coordinate array
1661: - coords - The cell coordinates

1663:   Level: developer

1665: .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1666: @*/
1667: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1668: {
1669:   DM           cdm;
1670:   PetscSection cs;
1671:   Vec          coordinates;

1673:   PetscFunctionBeginHot;
1674:   if (*isDG) {
1675:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1676:     PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1677:   } else {
1678:     PetscCall(DMGetCoordinateDM(dm, &cdm));
1679:     PetscCall(DMGetCoordinateSection(dm, &cs));
1680:     PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1681:     PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, (PetscScalar **)coords));
1682:   }
1683:   PetscFunctionReturn(PETSC_SUCCESS);
1684: }

1686: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1687: {
1688:   const PetscScalar *array;
1689:   PetscScalar       *coords = NULL;
1690:   PetscInt           numCoords, d;
1691:   PetscBool          isDG;

1693:   PetscFunctionBegin;
1694:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1695:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1696:   *detJ = 0.0;
1697:   if (numCoords == 6) {
1698:     const PetscInt dim = 3;
1699:     PetscReal      R[9], J0;

1701:     if (v0) {
1702:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1703:     }
1704:     PetscCall(DMPlexComputeProjection3Dto1D(coords, R));
1705:     if (J) {
1706:       J0   = 0.5 * PetscRealPart(coords[1]);
1707:       J[0] = R[0] * J0;
1708:       J[1] = R[1];
1709:       J[2] = R[2];
1710:       J[3] = R[3] * J0;
1711:       J[4] = R[4];
1712:       J[5] = R[5];
1713:       J[6] = R[6] * J0;
1714:       J[7] = R[7];
1715:       J[8] = R[8];
1716:       DMPlex_Det3D_Internal(detJ, J);
1717:       if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1718:     }
1719:   } else if (numCoords == 4) {
1720:     const PetscInt dim = 2;
1721:     PetscReal      R[4], J0;

1723:     if (v0) {
1724:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1725:     }
1726:     PetscCall(DMPlexComputeProjection2Dto1D(coords, R));
1727:     if (J) {
1728:       J0   = 0.5 * PetscRealPart(coords[1]);
1729:       J[0] = R[0] * J0;
1730:       J[1] = R[1];
1731:       J[2] = R[2] * J0;
1732:       J[3] = R[3];
1733:       DMPlex_Det2D_Internal(detJ, J);
1734:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1735:     }
1736:   } else if (numCoords == 2) {
1737:     const PetscInt dim = 1;

1739:     if (v0) {
1740:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1741:     }
1742:     if (J) {
1743:       J[0]  = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1744:       *detJ = J[0];
1745:       PetscCall(PetscLogFlops(2.0));
1746:       if (invJ) {
1747:         invJ[0] = 1.0 / J[0];
1748:         PetscCall(PetscLogFlops(1.0));
1749:       }
1750:     }
1751:   } 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);
1752:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1753:   PetscFunctionReturn(PETSC_SUCCESS);
1754: }

1756: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1757: {
1758:   const PetscScalar *array;
1759:   PetscScalar       *coords = NULL;
1760:   PetscInt           numCoords, d;
1761:   PetscBool          isDG;

1763:   PetscFunctionBegin;
1764:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1765:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1766:   *detJ = 0.0;
1767:   if (numCoords == 9) {
1768:     const PetscInt dim = 3;
1769:     PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1771:     if (v0) {
1772:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1773:     }
1774:     PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1775:     if (J) {
1776:       const PetscInt pdim = 2;

1778:       for (d = 0; d < pdim; d++) {
1779:         for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1780:       }
1781:       PetscCall(PetscLogFlops(8.0));
1782:       DMPlex_Det3D_Internal(detJ, J0);
1783:       for (d = 0; d < dim; d++) {
1784:         for (PetscInt f = 0; f < dim; f++) {
1785:           J[d * dim + f] = 0.0;
1786:           for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1787:         }
1788:       }
1789:       PetscCall(PetscLogFlops(18.0));
1790:     }
1791:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1792:   } else if (numCoords == 6) {
1793:     const PetscInt dim = 2;

1795:     if (v0) {
1796:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1797:     }
1798:     if (J) {
1799:       for (d = 0; d < dim; d++) {
1800:         for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1801:       }
1802:       PetscCall(PetscLogFlops(8.0));
1803:       DMPlex_Det2D_Internal(detJ, J);
1804:     }
1805:     if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1806:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1807:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1808:   PetscFunctionReturn(PETSC_SUCCESS);
1809: }

1811: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1812: {
1813:   const PetscScalar *array;
1814:   PetscScalar       *coords = NULL;
1815:   PetscInt           numCoords, d;
1816:   PetscBool          isDG;

1818:   PetscFunctionBegin;
1819:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1820:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1821:   if (!Nq) {
1822:     PetscInt vorder[4] = {0, 1, 2, 3};

1824:     if (isTensor) {
1825:       vorder[2] = 3;
1826:       vorder[3] = 2;
1827:     }
1828:     *detJ = 0.0;
1829:     if (numCoords == 12) {
1830:       const PetscInt dim = 3;
1831:       PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1833:       if (v) {
1834:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1835:       }
1836:       PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1837:       if (J) {
1838:         const PetscInt pdim = 2;

1840:         for (d = 0; d < pdim; d++) {
1841:           J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1842:           J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1843:         }
1844:         PetscCall(PetscLogFlops(8.0));
1845:         DMPlex_Det3D_Internal(detJ, J0);
1846:         for (d = 0; d < dim; d++) {
1847:           for (PetscInt f = 0; f < dim; f++) {
1848:             J[d * dim + f] = 0.0;
1849:             for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1850:           }
1851:         }
1852:         PetscCall(PetscLogFlops(18.0));
1853:       }
1854:       if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1855:     } else if (numCoords == 8) {
1856:       const PetscInt dim = 2;

1858:       if (v) {
1859:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1860:       }
1861:       if (J) {
1862:         for (d = 0; d < dim; d++) {
1863:           J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1864:           J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1865:         }
1866:         PetscCall(PetscLogFlops(8.0));
1867:         DMPlex_Det2D_Internal(detJ, J);
1868:       }
1869:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1870:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1871:   } else {
1872:     const PetscInt Nv         = 4;
1873:     const PetscInt dimR       = 2;
1874:     PetscInt       zToPlex[4] = {0, 1, 3, 2};
1875:     PetscReal      zOrder[12];
1876:     PetscReal      zCoeff[12];
1877:     PetscInt       i, j, k, l, dim;

1879:     if (isTensor) {
1880:       zToPlex[2] = 2;
1881:       zToPlex[3] = 3;
1882:     }
1883:     if (numCoords == 12) {
1884:       dim = 3;
1885:     } else if (numCoords == 8) {
1886:       dim = 2;
1887:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1888:     for (i = 0; i < Nv; i++) {
1889:       PetscInt zi = zToPlex[i];

1891:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1892:     }
1893:     for (j = 0; j < dim; j++) {
1894:       /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
1895:            \phi^0 = (1 - xi - eta + xi eta) --> 1      = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
1896:            \phi^1 = (1 + xi - eta - xi eta) --> xi     = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
1897:            \phi^2 = (1 - xi + eta - xi eta) --> eta    = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
1898:            \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
1899:       */
1900:       zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1901:       zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1902:       zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1903:       zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1904:     }
1905:     for (i = 0; i < Nq; i++) {
1906:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1908:       if (v) {
1909:         PetscReal extPoint[4];

1911:         extPoint[0] = 1.;
1912:         extPoint[1] = xi;
1913:         extPoint[2] = eta;
1914:         extPoint[3] = xi * eta;
1915:         for (j = 0; j < dim; j++) {
1916:           PetscReal val = 0.;

1918:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1919:           v[i * dim + j] = val;
1920:         }
1921:       }
1922:       if (J) {
1923:         PetscReal extJ[8];

1925:         extJ[0] = 0.;
1926:         extJ[1] = 0.;
1927:         extJ[2] = 1.;
1928:         extJ[3] = 0.;
1929:         extJ[4] = 0.;
1930:         extJ[5] = 1.;
1931:         extJ[6] = eta;
1932:         extJ[7] = xi;
1933:         for (j = 0; j < dim; j++) {
1934:           for (k = 0; k < dimR; k++) {
1935:             PetscReal val = 0.;

1937:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1938:             J[i * dim * dim + dim * j + k] = val;
1939:           }
1940:         }
1941:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1942:           PetscReal  x, y, z;
1943:           PetscReal *iJ = &J[i * dim * dim];
1944:           PetscReal  norm;

1946:           x     = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1947:           y     = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1948:           z     = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1949:           norm  = PetscSqrtReal(x * x + y * y + z * z);
1950:           iJ[2] = x / norm;
1951:           iJ[5] = y / norm;
1952:           iJ[8] = z / norm;
1953:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1954:           if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1955:         } else {
1956:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1957:           if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1958:         }
1959:       }
1960:     }
1961:   }
1962:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1963:   PetscFunctionReturn(PETSC_SUCCESS);
1964: }

1966: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1967: {
1968:   const PetscScalar *array;
1969:   PetscScalar       *coords = NULL;
1970:   const PetscInt     dim    = 3;
1971:   PetscInt           numCoords, d;
1972:   PetscBool          isDG;

1974:   PetscFunctionBegin;
1975:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1976:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1977:   *detJ = 0.0;
1978:   if (v0) {
1979:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1980:   }
1981:   if (J) {
1982:     for (d = 0; d < dim; d++) {
1983:       /* I orient with outward face normals */
1984:       J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1985:       J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1986:       J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1987:     }
1988:     PetscCall(PetscLogFlops(18.0));
1989:     DMPlex_Det3D_Internal(detJ, J);
1990:   }
1991:   if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1992:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1993:   PetscFunctionReturn(PETSC_SUCCESS);
1994: }

1996: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1997: {
1998:   const PetscScalar *array;
1999:   PetscScalar       *coords = NULL;
2000:   const PetscInt     dim    = 3;
2001:   PetscInt           numCoords, d;
2002:   PetscBool          isDG;

2004:   PetscFunctionBegin;
2005:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2006:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2007:   if (!Nq) {
2008:     *detJ = 0.0;
2009:     if (v) {
2010:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2011:     }
2012:     if (J) {
2013:       for (d = 0; d < dim; d++) {
2014:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2015:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2016:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2017:       }
2018:       PetscCall(PetscLogFlops(18.0));
2019:       DMPlex_Det3D_Internal(detJ, J);
2020:     }
2021:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2022:   } else {
2023:     const PetscInt Nv         = 8;
2024:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
2025:     const PetscInt dim        = 3;
2026:     const PetscInt dimR       = 3;
2027:     PetscReal      zOrder[24];
2028:     PetscReal      zCoeff[24];
2029:     PetscInt       i, j, k, l;

2031:     for (i = 0; i < Nv; i++) {
2032:       PetscInt zi = zToPlex[i];

2034:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
2035:     }
2036:     for (j = 0; j < dim; j++) {
2037:       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]);
2038:       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]);
2039:       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]);
2040:       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]);
2041:       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]);
2042:       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]);
2043:       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]);
2044:       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]);
2045:     }
2046:     for (i = 0; i < Nq; i++) {
2047:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

2049:       if (v) {
2050:         PetscReal extPoint[8];

2052:         extPoint[0] = 1.;
2053:         extPoint[1] = xi;
2054:         extPoint[2] = eta;
2055:         extPoint[3] = xi * eta;
2056:         extPoint[4] = theta;
2057:         extPoint[5] = theta * xi;
2058:         extPoint[6] = theta * eta;
2059:         extPoint[7] = theta * eta * xi;
2060:         for (j = 0; j < dim; j++) {
2061:           PetscReal val = 0.;

2063:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
2064:           v[i * dim + j] = val;
2065:         }
2066:       }
2067:       if (J) {
2068:         PetscReal extJ[24];

2070:         extJ[0]  = 0.;
2071:         extJ[1]  = 0.;
2072:         extJ[2]  = 0.;
2073:         extJ[3]  = 1.;
2074:         extJ[4]  = 0.;
2075:         extJ[5]  = 0.;
2076:         extJ[6]  = 0.;
2077:         extJ[7]  = 1.;
2078:         extJ[8]  = 0.;
2079:         extJ[9]  = eta;
2080:         extJ[10] = xi;
2081:         extJ[11] = 0.;
2082:         extJ[12] = 0.;
2083:         extJ[13] = 0.;
2084:         extJ[14] = 1.;
2085:         extJ[15] = theta;
2086:         extJ[16] = 0.;
2087:         extJ[17] = xi;
2088:         extJ[18] = 0.;
2089:         extJ[19] = theta;
2090:         extJ[20] = eta;
2091:         extJ[21] = theta * eta;
2092:         extJ[22] = theta * xi;
2093:         extJ[23] = eta * xi;

2095:         for (j = 0; j < dim; j++) {
2096:           for (k = 0; k < dimR; k++) {
2097:             PetscReal val = 0.;

2099:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
2100:             J[i * dim * dim + dim * j + k] = val;
2101:           }
2102:         }
2103:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2104:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2105:       }
2106:     }
2107:   }
2108:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2109:   PetscFunctionReturn(PETSC_SUCCESS);
2110: }

2112: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2113: {
2114:   const PetscScalar *array;
2115:   PetscScalar       *coords = NULL;
2116:   const PetscInt     dim    = 3;
2117:   PetscInt           numCoords, d;
2118:   PetscBool          isDG;

2120:   PetscFunctionBegin;
2121:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2122:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
2123:   if (!Nq) {
2124:     /* Assume that the map to the reference is affine */
2125:     *detJ = 0.0;
2126:     if (v) {
2127:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
2128:     }
2129:     if (J) {
2130:       for (d = 0; d < dim; d++) {
2131:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2132:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2133:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
2134:       }
2135:       PetscCall(PetscLogFlops(18.0));
2136:       DMPlex_Det3D_Internal(detJ, J);
2137:     }
2138:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
2139:   } else {
2140:     const PetscInt dim  = 3;
2141:     const PetscInt dimR = 3;
2142:     const PetscInt Nv   = 6;
2143:     PetscReal      verts[18];
2144:     PetscReal      coeff[18];
2145:     PetscInt       i, j, k, l;

2147:     for (i = 0; i < Nv; ++i)
2148:       for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
2149:     for (j = 0; j < dim; ++j) {
2150:       /* Check for triangle,
2151:            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)
2152:            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)
2153:            phi^2 =  1/2 (1 + eta)   chi^2 = delta(-1,  1)

2155:            phi^0 + phi^1 + phi^2 = 1    coef_1   = 1/2 (         chi^1 + chi^2)
2156:           -phi^0 + phi^1 - phi^2 = xi   coef_xi  = 1/2 (-chi^0 + chi^1)
2157:           -phi^0 - phi^1 + phi^2 = eta  coef_eta = 1/2 (-chi^0         + chi^2)

2159:           < chi_0 chi_1 chi_2> A /  1  1  1 \ / phi_0 \   <chi> I <phi>^T  so we need the inverse transpose
2160:                                  | -1  1 -1 | | phi_1 | =
2161:                                  \ -1 -1  1 / \ phi_2 /

2163:           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
2164:       */
2165:       /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
2166:            \phi^0 = 1/4 (   -xi - eta        + xi zeta + eta zeta) --> /  1  1  1  1  1  1 \ 1
2167:            \phi^1 = 1/4 (1      + eta - zeta           - eta zeta) --> | -1  1 -1 -1 -1  1 | eta
2168:            \phi^2 = 1/4 (1 + xi       - zeta - xi zeta)            --> | -1 -1  1 -1  1 -1 | xi
2169:            \phi^3 = 1/4 (   -xi - eta        - xi zeta - eta zeta) --> | -1 -1 -1  1  1  1 | zeta
2170:            \phi^4 = 1/4 (1 + xi       + zeta + xi zeta)            --> |  1  1 -1 -1  1 -1 | xi zeta
2171:            \phi^5 = 1/4 (1      + eta + zeta           + eta zeta) --> \  1 -1  1 -1 -1  1 / eta zeta
2172:            1/4 /  0  1  1  0  1  1 \
2173:                | -1  1  0 -1  0  1 |
2174:                | -1  0  1 -1  1  0 |
2175:                |  0 -1 -1  0  1  1 |
2176:                |  1  0 -1 -1  1  0 |
2177:                \  1 -1  0 -1  0  1 /
2178:       */
2179:       coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2180:       coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2181:       coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2182:       coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
2183:       coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
2184:       coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
2185:       /* For reference prism:
2186:       {0, 0, 0}
2187:       {0, 1, 0}
2188:       {1, 0, 0}
2189:       {0, 0, 1}
2190:       {0, 0, 0}
2191:       {0, 0, 0}
2192:       */
2193:     }
2194:     for (i = 0; i < Nq; ++i) {
2195:       const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];

2197:       if (v) {
2198:         PetscReal extPoint[6];
2199:         PetscInt  c;

2201:         extPoint[0] = 1.;
2202:         extPoint[1] = eta;
2203:         extPoint[2] = xi;
2204:         extPoint[3] = zeta;
2205:         extPoint[4] = xi * zeta;
2206:         extPoint[5] = eta * zeta;
2207:         for (c = 0; c < dim; ++c) {
2208:           PetscReal val = 0.;

2210:           for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
2211:           v[i * dim + c] = val;
2212:         }
2213:       }
2214:       if (J) {
2215:         PetscReal extJ[18];

2217:         extJ[0]  = 0.;
2218:         extJ[1]  = 0.;
2219:         extJ[2]  = 0.;
2220:         extJ[3]  = 0.;
2221:         extJ[4]  = 1.;
2222:         extJ[5]  = 0.;
2223:         extJ[6]  = 1.;
2224:         extJ[7]  = 0.;
2225:         extJ[8]  = 0.;
2226:         extJ[9]  = 0.;
2227:         extJ[10] = 0.;
2228:         extJ[11] = 1.;
2229:         extJ[12] = zeta;
2230:         extJ[13] = 0.;
2231:         extJ[14] = xi;
2232:         extJ[15] = 0.;
2233:         extJ[16] = zeta;
2234:         extJ[17] = eta;

2236:         for (j = 0; j < dim; j++) {
2237:           for (k = 0; k < dimR; k++) {
2238:             PetscReal val = 0.;

2240:             for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
2241:             J[i * dim * dim + dim * j + k] = val;
2242:           }
2243:         }
2244:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
2245:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
2246:       }
2247:     }
2248:   }
2249:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
2250:   PetscFunctionReturn(PETSC_SUCCESS);
2251: }

2253: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2254: {
2255:   DMPolytopeType   ct;
2256:   PetscInt         depth, dim, coordDim, coneSize, i;
2257:   PetscInt         Nq     = 0;
2258:   const PetscReal *points = NULL;
2259:   DMLabel          depthLabel;
2260:   PetscReal        xi0[3]   = {-1., -1., -1.}, v0[3], J0[9], detJ0;
2261:   PetscBool        isAffine = PETSC_TRUE;

2263:   PetscFunctionBegin;
2264:   PetscCall(DMPlexGetDepth(dm, &depth));
2265:   PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
2266:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
2267:   PetscCall(DMLabelGetValue(depthLabel, cell, &dim));
2268:   if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim));
2269:   PetscCall(DMGetCoordinateDim(dm, &coordDim));
2270:   PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim);
2271:   if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL));
2272:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
2273:   switch (ct) {
2274:   case DM_POLYTOPE_POINT:
2275:     PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ));
2276:     isAffine = PETSC_FALSE;
2277:     break;
2278:   case DM_POLYTOPE_SEGMENT:
2279:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2280:     if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2281:     else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ));
2282:     break;
2283:   case DM_POLYTOPE_TRIANGLE:
2284:     if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2285:     else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ));
2286:     break;
2287:   case DM_POLYTOPE_QUADRILATERAL:
2288:     PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ));
2289:     isAffine = PETSC_FALSE;
2290:     break;
2291:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2292:     PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ));
2293:     isAffine = PETSC_FALSE;
2294:     break;
2295:   case DM_POLYTOPE_TETRAHEDRON:
2296:     if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
2297:     else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ));
2298:     break;
2299:   case DM_POLYTOPE_HEXAHEDRON:
2300:     PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2301:     isAffine = PETSC_FALSE;
2302:     break;
2303:   case DM_POLYTOPE_TRI_PRISM:
2304:     PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
2305:     isAffine = PETSC_FALSE;
2306:     break;
2307:   default:
2308:     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))]);
2309:   }
2310:   if (isAffine && Nq) {
2311:     if (v) {
2312:       for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
2313:     }
2314:     if (detJ) {
2315:       for (i = 0; i < Nq; i++) detJ[i] = detJ0;
2316:     }
2317:     if (J) {
2318:       PetscInt k;

2320:       for (i = 0, k = 0; i < Nq; i++) {
2321:         PetscInt j;

2323:         for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
2324:       }
2325:     }
2326:     if (invJ) {
2327:       PetscInt k;
2328:       switch (coordDim) {
2329:       case 0:
2330:         break;
2331:       case 1:
2332:         invJ[0] = 1. / J0[0];
2333:         break;
2334:       case 2:
2335:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
2336:         break;
2337:       case 3:
2338:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
2339:         break;
2340:       }
2341:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
2342:         PetscInt j;

2344:         for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
2345:       }
2346:     }
2347:   }
2348:   PetscFunctionReturn(PETSC_SUCCESS);
2349: }

2351: /*@C
2352:   DMPlexComputeCellGeometryAffineFEM - Assuming an affine map, compute the Jacobian, inverse Jacobian, and Jacobian determinant for a given cell

2354:   Collective

2356:   Input Parameters:
2357: + dm   - the `DMPLEX`
2358: - cell - the cell

2360:   Output Parameters:
2361: + v0   - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
2362: . J    - the Jacobian of the transform from the reference element
2363: . invJ - the inverse of the Jacobian
2364: - detJ - the Jacobian determinant

2366:   Level: advanced

2368: .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2369: @*/
2370: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2371: {
2372:   PetscFunctionBegin;
2373:   PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ));
2374:   PetscFunctionReturn(PETSC_SUCCESS);
2375: }

2377: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2378: {
2379:   const PetscScalar *array;
2380:   PetscScalar       *coords = NULL;
2381:   PetscInt           numCoords;
2382:   PetscBool          isDG;
2383:   PetscQuadrature    feQuad;
2384:   const PetscReal   *quadPoints;
2385:   PetscTabulation    T;
2386:   PetscInt           dim, cdim, pdim, qdim, Nq, q;

2388:   PetscFunctionBegin;
2389:   PetscCall(DMGetDimension(dm, &dim));
2390:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2391:   PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2392:   if (!quad) { /* use the first point of the first functional of the dual space */
2393:     PetscDualSpace dsp;

2395:     PetscCall(PetscFEGetDualSpace(fe, &dsp));
2396:     PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad));
2397:     PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2398:     Nq = 1;
2399:   } else {
2400:     PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2401:   }
2402:   PetscCall(PetscFEGetDimension(fe, &pdim));
2403:   PetscCall(PetscFEGetQuadrature(fe, &feQuad));
2404:   if (feQuad == quad) {
2405:     PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T));
2406:     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);
2407:   } else {
2408:     PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T));
2409:   }
2410:   PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim);
2411:   {
2412:     const PetscReal *basis    = T->T[0];
2413:     const PetscReal *basisDer = T->T[1];
2414:     PetscReal        detJt;

2416:     PetscAssert(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np);
2417:     PetscAssert(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb);
2418:     PetscAssert(dim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->Nc);
2419:     PetscAssert(cdim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->cdim);
2420:     if (v) {
2421:       PetscCall(PetscArrayzero(v, Nq * cdim));
2422:       for (q = 0; q < Nq; ++q) {
2423:         PetscInt i, k;

2425:         for (k = 0; k < pdim; ++k) {
2426:           const PetscInt vertex = k / cdim;
2427:           for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2428:         }
2429:         PetscCall(PetscLogFlops(2.0 * pdim * cdim));
2430:       }
2431:     }
2432:     if (J) {
2433:       PetscCall(PetscArrayzero(J, Nq * cdim * cdim));
2434:       for (q = 0; q < Nq; ++q) {
2435:         PetscInt i, j, k, c, r;

2437:         /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2438:         for (k = 0; k < pdim; ++k) {
2439:           const PetscInt vertex = k / cdim;
2440:           for (j = 0; j < dim; ++j) {
2441:             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]);
2442:           }
2443:         }
2444:         PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim));
2445:         if (cdim > dim) {
2446:           for (c = dim; c < cdim; ++c)
2447:             for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2448:         }
2449:         if (!detJ && !invJ) continue;
2450:         detJt = 0.;
2451:         switch (cdim) {
2452:         case 3:
2453:           DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2454:           if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2455:           break;
2456:         case 2:
2457:           DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2458:           if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2459:           break;
2460:         case 1:
2461:           detJt = J[q * cdim * dim];
2462:           if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2463:         }
2464:         if (detJ) detJ[q] = detJt;
2465:       }
2466:     } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
2467:   }
2468:   if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T));
2469:   PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2470:   PetscFunctionReturn(PETSC_SUCCESS);
2471: }

2473: /*@C
2474:   DMPlexComputeCellGeometryFEM - Compute the Jacobian, inverse Jacobian, and Jacobian determinant at each quadrature point in the given cell

2476:   Collective

2478:   Input Parameters:
2479: + dm   - the `DMPLEX`
2480: . cell - the cell
2481: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If `quad` is `NULL`, geometry will be
2482:          evaluated at the first vertex of the reference element

2484:   Output Parameters:
2485: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
2486: . J    - the Jacobian of the transform from the reference element at each quadrature point
2487: . invJ - the inverse of the Jacobian at each quadrature point
2488: - detJ - the Jacobian determinant at each quadrature point

2490:   Level: advanced

2492: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2493: @*/
2494: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2495: {
2496:   DM      cdm;
2497:   PetscFE fe = NULL;

2499:   PetscFunctionBegin;
2500:   PetscAssertPointer(detJ, 7);
2501:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2502:   if (cdm) {
2503:     PetscClassId id;
2504:     PetscInt     numFields;
2505:     PetscDS      prob;
2506:     PetscObject  disc;

2508:     PetscCall(DMGetNumFields(cdm, &numFields));
2509:     if (numFields) {
2510:       PetscCall(DMGetDS(cdm, &prob));
2511:       PetscCall(PetscDSGetDiscretization(prob, 0, &disc));
2512:       PetscCall(PetscObjectGetClassId(disc, &id));
2513:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2514:     }
2515:   }
2516:   if (!fe) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ));
2517:   else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ));
2518:   PetscFunctionReturn(PETSC_SUCCESS);
2519: }

2521: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2522: {
2523:   PetscSection       coordSection;
2524:   Vec                coordinates;
2525:   const PetscScalar *coords = NULL;
2526:   PetscInt           d, dof, off;

2528:   PetscFunctionBegin;
2529:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2530:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2531:   PetscCall(VecGetArrayRead(coordinates, &coords));

2533:   /* for a point the centroid is just the coord */
2534:   if (centroid) {
2535:     PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2536:     PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2537:     for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2538:   }
2539:   if (normal) {
2540:     const PetscInt *support, *cones;
2541:     PetscInt        supportSize;
2542:     PetscReal       norm, sign;

2544:     /* compute the norm based upon the support centroids */
2545:     PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize));
2546:     PetscCall(DMPlexGetSupport(dm, cell, &support));
2547:     PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL));

2549:     /* Take the normal from the centroid of the support to the vertex*/
2550:     PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2551:     PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2552:     for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]);

2554:     /* Determine the sign of the normal based upon its location in the support */
2555:     PetscCall(DMPlexGetCone(dm, support[0], &cones));
2556:     sign = cones[0] == cell ? 1.0 : -1.0;

2558:     norm = DMPlex_NormD_Internal(dim, normal);
2559:     for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2560:   }
2561:   if (vol) *vol = 1.0;
2562:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
2563:   PetscFunctionReturn(PETSC_SUCCESS);
2564: }

2566: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2567: {
2568:   const PetscScalar *array;
2569:   PetscScalar       *coords = NULL;
2570:   PetscInt           cdim, coordSize, d;
2571:   PetscBool          isDG;

2573:   PetscFunctionBegin;
2574:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2575:   PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2576:   PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2);
2577:   if (centroid) {
2578:     for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2579:   }
2580:   if (normal) {
2581:     PetscReal norm;

2583:     switch (cdim) {
2584:     case 3:
2585:       normal[2] = 0.; /* fall through */
2586:     case 2:
2587:       normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2588:       normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2589:       break;
2590:     case 1:
2591:       normal[0] = 1.0;
2592:       break;
2593:     default:
2594:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2595:     }
2596:     norm = DMPlex_NormD_Internal(cdim, normal);
2597:     for (d = 0; d < cdim; ++d) normal[d] /= norm;
2598:   }
2599:   if (vol) {
2600:     *vol = 0.0;
2601:     for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2602:     *vol = PetscSqrtReal(*vol);
2603:   }
2604:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2605:   PetscFunctionReturn(PETSC_SUCCESS);
2606: }

2608: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2609: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2610: {
2611:   DMPolytopeType     ct;
2612:   const PetscScalar *array;
2613:   PetscScalar       *coords = NULL;
2614:   PetscInt           coordSize;
2615:   PetscBool          isDG;
2616:   PetscInt           fv[4] = {0, 1, 2, 3};
2617:   PetscInt           cdim, numCorners, p, d;

2619:   PetscFunctionBegin;
2620:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2621:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
2622:   switch (ct) {
2623:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2624:     fv[2] = 3;
2625:     fv[3] = 2;
2626:     break;
2627:   default:
2628:     break;
2629:   }
2630:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2631:   PetscCall(DMPlexGetConeSize(dm, cell, &numCorners));
2632:   PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2633:   {
2634:     PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;

2636:     for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2637:     for (p = 0; p < numCorners - 2; ++p) {
2638:       PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2639:       for (d = 0; d < cdim; d++) {
2640:         e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2641:         e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2642:       }
2643:       const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2644:       const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2645:       const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2646:       const PetscReal a  = PetscSqrtReal(dx * dx + dy * dy + dz * dz);

2648:       n[0] += dx;
2649:       n[1] += dy;
2650:       n[2] += dz;
2651:       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.;
2652:     }
2653:     norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2654:     // Allow zero volume cells
2655:     if (norm != 0) {
2656:       n[0] /= norm;
2657:       n[1] /= norm;
2658:       n[2] /= norm;
2659:       c[0] /= norm;
2660:       c[1] /= norm;
2661:       c[2] /= norm;
2662:     }
2663:     if (vol) *vol = 0.5 * norm;
2664:     if (centroid)
2665:       for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2666:     if (normal)
2667:       for (d = 0; d < cdim; ++d) normal[d] = n[d];
2668:   }
2669:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2670:   PetscFunctionReturn(PETSC_SUCCESS);
2671: }

2673: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2674: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2675: {
2676:   DMPolytopeType        ct;
2677:   const PetscScalar    *array;
2678:   PetscScalar          *coords = NULL;
2679:   PetscInt              coordSize;
2680:   PetscBool             isDG;
2681:   PetscReal             vsum      = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2682:   const PetscInt        order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2683:   const PetscInt       *cone, *faceSizes, *faces;
2684:   const DMPolytopeType *faceTypes;
2685:   PetscBool             isHybrid = PETSC_FALSE;
2686:   PetscInt              numFaces, f, fOff = 0, p, d;

2688:   PetscFunctionBegin;
2689:   PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim);
2690:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2691:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
2692:   switch (ct) {
2693:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2694:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2695:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
2696:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2697:     isHybrid = PETSC_TRUE;
2698:   default:
2699:     break;
2700:   }

2702:   if (centroid)
2703:     for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2704:   PetscCall(DMPlexGetCone(dm, cell, &cone));

2706:   // Using the closure of faces for coordinates does not work in periodic geometries, so we index into the cell coordinates
2707:   PetscCall(DMPlexGetRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2708:   PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2709:   for (f = 0; f < numFaces; ++f) {
2710:     PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */

2712:     // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2713:     // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2714:     // so that all tetrahedra have positive volume.
2715:     if (f == 0)
2716:       for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2717:     switch (faceTypes[f]) {
2718:     case DM_POLYTOPE_TRIANGLE:
2719:       for (d = 0; d < dim; ++d) {
2720:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2721:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2722:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2723:       }
2724:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2725:       if (flip) vtmp = -vtmp;
2726:       vsum += vtmp;
2727:       if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2728:         for (d = 0; d < dim; ++d) {
2729:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2730:         }
2731:       }
2732:       break;
2733:     case DM_POLYTOPE_QUADRILATERAL:
2734:     case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2735:       PetscInt fv[4] = {0, 1, 2, 3};

2737:       /* Side faces for hybrid cells are stored as tensor products */
2738:       if (isHybrid && f > 1) {
2739:         fv[2] = 3;
2740:         fv[3] = 2;
2741:       }
2742:       /* DO FOR PYRAMID */
2743:       /* First tet */
2744:       for (d = 0; d < dim; ++d) {
2745:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2746:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2747:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2748:       }
2749:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2750:       if (flip) vtmp = -vtmp;
2751:       vsum += vtmp;
2752:       if (centroid) {
2753:         for (d = 0; d < dim; ++d) {
2754:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2755:         }
2756:       }
2757:       /* Second tet */
2758:       for (d = 0; d < dim; ++d) {
2759:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2760:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2761:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2762:       }
2763:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2764:       if (flip) vtmp = -vtmp;
2765:       vsum += vtmp;
2766:       if (centroid) {
2767:         for (d = 0; d < dim; ++d) {
2768:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2769:         }
2770:       }
2771:       break;
2772:     }
2773:     default:
2774:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2775:     }
2776:     fOff += faceSizes[f];
2777:   }
2778:   PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2779:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2780:   if (vol) *vol = PetscAbsReal(vsum);
2781:   if (normal)
2782:     for (d = 0; d < dim; ++d) normal[d] = 0.0;
2783:   if (centroid)
2784:     for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2785:   PetscFunctionReturn(PETSC_SUCCESS);
2786: }

2788: /*@C
2789:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

2791:   Collective

2793:   Input Parameters:
2794: + dm   - the `DMPLEX`
2795: - cell - the cell

2797:   Output Parameters:
2798: + vol      - the cell volume
2799: . centroid - the cell centroid
2800: - normal   - the cell normal, if appropriate

2802:   Level: advanced

2804: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2805: @*/
2806: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2807: {
2808:   PetscInt depth, dim;

2810:   PetscFunctionBegin;
2811:   PetscCall(DMPlexGetDepth(dm, &depth));
2812:   PetscCall(DMGetDimension(dm, &dim));
2813:   PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2814:   PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
2815:   switch (depth) {
2816:   case 0:
2817:     PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal));
2818:     break;
2819:   case 1:
2820:     PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal));
2821:     break;
2822:   case 2:
2823:     PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal));
2824:     break;
2825:   case 3:
2826:     PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal));
2827:     break;
2828:   default:
2829:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2830:   }
2831:   PetscFunctionReturn(PETSC_SUCCESS);
2832: }

2834: /*@
2835:   DMPlexComputeGeometryFVM - Computes the cell and face geometry for a finite volume method

2837:   Input Parameter:
2838: . dm - The `DMPLEX`

2840:   Output Parameters:
2841: + cellgeom - A `Vec` of `PetscFVCellGeom` data
2842: - facegeom - A `Vec` of `PetscFVFaceGeom` data

2844:   Level: developer

2846: .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom`
2847: @*/
2848: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2849: {
2850:   DM           dmFace, dmCell;
2851:   DMLabel      ghostLabel;
2852:   PetscSection sectionFace, sectionCell;
2853:   PetscSection coordSection;
2854:   Vec          coordinates;
2855:   PetscScalar *fgeom, *cgeom;
2856:   PetscReal    minradius, gminradius;
2857:   PetscInt     dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

2859:   PetscFunctionBegin;
2860:   PetscCall(DMGetDimension(dm, &dim));
2861:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2862:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2863:   /* Make cell centroids and volumes */
2864:   PetscCall(DMClone(dm, &dmCell));
2865:   PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection));
2866:   PetscCall(DMSetCoordinatesLocal(dmCell, coordinates));
2867:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
2868:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2869:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
2870:   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
2871:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar))));
2872:   PetscCall(PetscSectionSetUp(sectionCell));
2873:   PetscCall(DMSetLocalSection(dmCell, sectionCell));
2874:   PetscCall(PetscSectionDestroy(&sectionCell));
2875:   PetscCall(DMCreateLocalVector(dmCell, cellgeom));
2876:   if (cEndInterior < 0) cEndInterior = cEnd;
2877:   PetscCall(VecGetArray(*cellgeom, &cgeom));
2878:   for (c = cStart; c < cEndInterior; ++c) {
2879:     PetscFVCellGeom *cg;

2881:     PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg));
2882:     PetscCall(PetscArrayzero(cg, 1));
2883:     PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL));
2884:   }
2885:   /* Compute face normals and minimum cell radius */
2886:   PetscCall(DMClone(dm, &dmFace));
2887:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace));
2888:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2889:   PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd));
2890:   for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar))));
2891:   PetscCall(PetscSectionSetUp(sectionFace));
2892:   PetscCall(DMSetLocalSection(dmFace, sectionFace));
2893:   PetscCall(PetscSectionDestroy(&sectionFace));
2894:   PetscCall(DMCreateLocalVector(dmFace, facegeom));
2895:   PetscCall(VecGetArray(*facegeom, &fgeom));
2896:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2897:   minradius = PETSC_MAX_REAL;
2898:   for (f = fStart; f < fEnd; ++f) {
2899:     PetscFVFaceGeom *fg;
2900:     PetscReal        area;
2901:     const PetscInt  *cells;
2902:     PetscInt         ncells, ghost = -1, d, numChildren;

2904:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2905:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2906:     PetscCall(DMPlexGetSupport(dm, f, &cells));
2907:     PetscCall(DMPlexGetSupportSize(dm, f, &ncells));
2908:     /* It is possible to get a face with no support when using partition overlap */
2909:     if (!ncells || ghost >= 0 || numChildren) continue;
2910:     PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg));
2911:     PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal));
2912:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2913:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2914:     {
2915:       PetscFVCellGeom *cL, *cR;
2916:       PetscReal       *lcentroid, *rcentroid;
2917:       PetscReal        l[3], r[3], v[3];

2919:       PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL));
2920:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2921:       if (ncells > 1) {
2922:         PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR));
2923:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2924:       } else {
2925:         rcentroid = fg->centroid;
2926:       }
2927:       PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l));
2928:       PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r));
2929:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2930:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2931:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2932:       }
2933:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2934:         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]);
2935:         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]);
2936:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f);
2937:       }
2938:       if (cells[0] < cEndInterior) {
2939:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2940:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2941:       }
2942:       if (ncells > 1 && cells[1] < cEndInterior) {
2943:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2944:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2945:       }
2946:     }
2947:   }
2948:   PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
2949:   PetscCall(DMPlexSetMinRadius(dm, gminradius));
2950:   /* Compute centroids of ghost cells */
2951:   for (c = cEndInterior; c < cEnd; ++c) {
2952:     PetscFVFaceGeom *fg;
2953:     const PetscInt  *cone, *support;
2954:     PetscInt         coneSize, supportSize, s;

2956:     PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize));
2957:     PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize);
2958:     PetscCall(DMPlexGetCone(dmCell, c, &cone));
2959:     PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize));
2960:     PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize);
2961:     PetscCall(DMPlexGetSupport(dmCell, cone[0], &support));
2962:     PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg));
2963:     for (s = 0; s < 2; ++s) {
2964:       /* Reflect ghost centroid across plane of face */
2965:       if (support[s] == c) {
2966:         PetscFVCellGeom *ci;
2967:         PetscFVCellGeom *cg;
2968:         PetscReal        c2f[3], a;

2970:         PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci));
2971:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2972:         a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2973:         PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg));
2974:         DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid);
2975:         cg->volume = ci->volume;
2976:       }
2977:     }
2978:   }
2979:   PetscCall(VecRestoreArray(*facegeom, &fgeom));
2980:   PetscCall(VecRestoreArray(*cellgeom, &cgeom));
2981:   PetscCall(DMDestroy(&dmCell));
2982:   PetscCall(DMDestroy(&dmFace));
2983:   PetscFunctionReturn(PETSC_SUCCESS);
2984: }

2986: /*@C
2987:   DMPlexGetMinRadius - Returns the minimum distance from any cell centroid to a face

2989:   Not Collective

2991:   Input Parameter:
2992: . dm - the `DMPLEX`

2994:   Output Parameter:
2995: . minradius - the minimum cell radius

2997:   Level: developer

2999: .seealso: `DMPLEX`, `DMGetCoordinates()`
3000: @*/
3001: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
3002: {
3003:   PetscFunctionBegin;
3005:   PetscAssertPointer(minradius, 2);
3006:   *minradius = ((DM_Plex *)dm->data)->minradius;
3007:   PetscFunctionReturn(PETSC_SUCCESS);
3008: }

3010: /*@C
3011:   DMPlexSetMinRadius - Sets the minimum distance from the cell centroid to a face

3013:   Logically Collective

3015:   Input Parameters:
3016: + dm        - the `DMPLEX`
3017: - minradius - the minimum cell radius

3019:   Level: developer

3021: .seealso: `DMPLEX`, `DMSetCoordinates()`
3022: @*/
3023: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
3024: {
3025:   PetscFunctionBegin;
3027:   ((DM_Plex *)dm->data)->minradius = minradius;
3028:   PetscFunctionReturn(PETSC_SUCCESS);
3029: }

3031: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3032: {
3033:   DMLabel      ghostLabel;
3034:   PetscScalar *dx, *grad, **gref;
3035:   PetscInt     dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

3037:   PetscFunctionBegin;
3038:   PetscCall(DMGetDimension(dm, &dim));
3039:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3040:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3041:   cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
3042:   PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL));
3043:   PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3044:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3045:   PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3046:   for (c = cStart; c < cEndInterior; c++) {
3047:     const PetscInt  *faces;
3048:     PetscInt         numFaces, usedFaces, f, d;
3049:     PetscFVCellGeom *cg;
3050:     PetscBool        boundary;
3051:     PetscInt         ghost;

3053:     // do not attempt to compute a gradient reconstruction stencil in a ghost cell.  It will never be used
3054:     PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3055:     if (ghost >= 0) continue;

3057:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3058:     PetscCall(DMPlexGetConeSize(dm, c, &numFaces));
3059:     PetscCall(DMPlexGetCone(dm, c, &faces));
3060:     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);
3061:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3062:       PetscFVCellGeom *cg1;
3063:       PetscFVFaceGeom *fg;
3064:       const PetscInt  *fcells;
3065:       PetscInt         ncell, side;

3067:       PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3068:       PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3069:       if ((ghost >= 0) || boundary) continue;
3070:       PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
3071:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
3072:       ncell = fcells[!side];    /* the neighbor */
3073:       PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg));
3074:       PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3075:       for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
3076:       gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
3077:     }
3078:     PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
3079:     PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad));
3080:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
3081:       PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
3082:       PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
3083:       if ((ghost >= 0) || boundary) continue;
3084:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
3085:       ++usedFaces;
3086:     }
3087:   }
3088:   PetscCall(PetscFree3(dx, grad, gref));
3089:   PetscFunctionReturn(PETSC_SUCCESS);
3090: }

3092: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
3093: {
3094:   DMLabel      ghostLabel;
3095:   PetscScalar *dx, *grad, **gref;
3096:   PetscInt     dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
3097:   PetscSection neighSec;
3098:   PetscInt(*neighbors)[2];
3099:   PetscInt *counter;

3101:   PetscFunctionBegin;
3102:   PetscCall(DMGetDimension(dm, &dim));
3103:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3104:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3105:   if (cEndInterior < 0) cEndInterior = cEnd;
3106:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec));
3107:   PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior));
3108:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
3109:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
3110:   for (f = fStart; f < fEnd; f++) {
3111:     const PetscInt *fcells;
3112:     PetscBool       boundary;
3113:     PetscInt        ghost = -1;
3114:     PetscInt        numChildren, numCells, c;

3116:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3117:     PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3118:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3119:     if ((ghost >= 0) || boundary || numChildren) continue;
3120:     PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3121:     if (numCells == 2) {
3122:       PetscCall(DMPlexGetSupport(dm, f, &fcells));
3123:       for (c = 0; c < 2; c++) {
3124:         PetscInt cell = fcells[c];

3126:         if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1));
3127:       }
3128:     }
3129:   }
3130:   PetscCall(PetscSectionSetUp(neighSec));
3131:   PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces));
3132:   PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
3133:   nStart = 0;
3134:   PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd));
3135:   PetscCall(PetscMalloc1((nEnd - nStart), &neighbors));
3136:   PetscCall(PetscCalloc1((cEndInterior - cStart), &counter));
3137:   for (f = fStart; f < fEnd; f++) {
3138:     const PetscInt *fcells;
3139:     PetscBool       boundary;
3140:     PetscInt        ghost = -1;
3141:     PetscInt        numChildren, numCells, c;

3143:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
3144:     PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
3145:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
3146:     if ((ghost >= 0) || boundary || numChildren) continue;
3147:     PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
3148:     if (numCells == 2) {
3149:       PetscCall(DMPlexGetSupport(dm, f, &fcells));
3150:       for (c = 0; c < 2; c++) {
3151:         PetscInt cell = fcells[c], off;

3153:         if (cell >= cStart && cell < cEndInterior) {
3154:           PetscCall(PetscSectionGetOffset(neighSec, cell, &off));
3155:           off += counter[cell - cStart]++;
3156:           neighbors[off][0] = f;
3157:           neighbors[off][1] = fcells[1 - c];
3158:         }
3159:       }
3160:     }
3161:   }
3162:   PetscCall(PetscFree(counter));
3163:   PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
3164:   for (c = cStart; c < cEndInterior; c++) {
3165:     PetscInt         numFaces, f, d, off, ghost = -1;
3166:     PetscFVCellGeom *cg;

3168:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
3169:     PetscCall(PetscSectionGetDof(neighSec, c, &numFaces));
3170:     PetscCall(PetscSectionGetOffset(neighSec, c, &off));

3172:     // do not attempt to compute a gradient reconstruction stencil in a ghost cell.  It will never be used
3173:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, c, &ghost));
3174:     if (ghost >= 0) continue;

3176:     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);
3177:     for (f = 0; f < numFaces; ++f) {
3178:       PetscFVCellGeom *cg1;
3179:       PetscFVFaceGeom *fg;
3180:       const PetscInt  *fcells;
3181:       PetscInt         ncell, side, nface;

3183:       nface = neighbors[off + f][0];
3184:       ncell = neighbors[off + f][1];
3185:       PetscCall(DMPlexGetSupport(dm, nface, &fcells));
3186:       side = (c != fcells[0]);
3187:       PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg));
3188:       PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
3189:       for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d];
3190:       gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
3191:     }
3192:     PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad));
3193:     for (f = 0; f < numFaces; ++f) {
3194:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d];
3195:     }
3196:   }
3197:   PetscCall(PetscFree3(dx, grad, gref));
3198:   PetscCall(PetscSectionDestroy(&neighSec));
3199:   PetscCall(PetscFree(neighbors));
3200:   PetscFunctionReturn(PETSC_SUCCESS);
3201: }

3203: /*@
3204:   DMPlexComputeGradientFVM - Compute geometric factors for gradient reconstruction, which are stored in the geometry data, and compute layout for gradient data

3206:   Collective

3208:   Input Parameters:
3209: + dm           - The `DMPLEX`
3210: . fvm          - The `PetscFV`
3211: - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()`

3213:   Input/Output Parameter:
3214: . faceGeometry - The face geometry from `DMPlexComputeFaceGeometryFVM()`; on output
3215:                  the geometric factors for gradient calculation are inserted

3217:   Output Parameter:
3218: . dmGrad - The `DM` describing the layout of gradient data

3220:   Level: developer

3222: .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
3223: @*/
3224: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
3225: {
3226:   DM           dmFace, dmCell;
3227:   PetscScalar *fgeom, *cgeom;
3228:   PetscSection sectionGrad, parentSection;
3229:   PetscInt     dim, pdim, cStart, cEnd, cEndInterior, c;

3231:   PetscFunctionBegin;
3232:   PetscCall(DMGetDimension(dm, &dim));
3233:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
3234:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
3235:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));
3236:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
3237:   PetscCall(VecGetDM(faceGeometry, &dmFace));
3238:   PetscCall(VecGetDM(cellGeometry, &dmCell));
3239:   PetscCall(VecGetArray(faceGeometry, &fgeom));
3240:   PetscCall(VecGetArray(cellGeometry, &cgeom));
3241:   PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL));
3242:   if (!parentSection) {
3243:     PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3244:   } else {
3245:     PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom));
3246:   }
3247:   PetscCall(VecRestoreArray(faceGeometry, &fgeom));
3248:   PetscCall(VecRestoreArray(cellGeometry, &cgeom));
3249:   /* Create storage for gradients */
3250:   PetscCall(DMClone(dm, dmGrad));
3251:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionGrad));
3252:   PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd));
3253:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim));
3254:   PetscCall(PetscSectionSetUp(sectionGrad));
3255:   PetscCall(DMSetLocalSection(*dmGrad, sectionGrad));
3256:   PetscCall(PetscSectionDestroy(&sectionGrad));
3257:   PetscFunctionReturn(PETSC_SUCCESS);
3258: }

3260: /*@
3261:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

3263:   Collective

3265:   Input Parameters:
3266: + dm - The `DM`
3267: - fv - The `PetscFV`

3269:   Output Parameters:
3270: + cellgeom - The cell geometry
3271: . facegeom - The face geometry
3272: - gradDM   - The gradient matrices

3274:   Level: developer

3276: .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()`
3277: @*/
3278: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
3279: {
3280:   PetscObject cellgeomobj, facegeomobj;

3282:   PetscFunctionBegin;
3283:   PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3284:   if (!cellgeomobj) {
3285:     Vec cellgeomInt, facegeomInt;

3287:     PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt));
3288:     PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt));
3289:     PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt));
3290:     PetscCall(VecDestroy(&cellgeomInt));
3291:     PetscCall(VecDestroy(&facegeomInt));
3292:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
3293:   }
3294:   PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj));
3295:   if (cellgeom) *cellgeom = (Vec)cellgeomobj;
3296:   if (facegeom) *facegeom = (Vec)facegeomobj;
3297:   if (gradDM) {
3298:     PetscObject gradobj;
3299:     PetscBool   computeGradients;

3301:     PetscCall(PetscFVGetComputeGradients(fv, &computeGradients));
3302:     if (!computeGradients) {
3303:       *gradDM = NULL;
3304:       PetscFunctionReturn(PETSC_SUCCESS);
3305:     }
3306:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3307:     if (!gradobj) {
3308:       DM dmGradInt;

3310:       PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt));
3311:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt));
3312:       PetscCall(DMDestroy(&dmGradInt));
3313:       PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
3314:     }
3315:     *gradDM = (DM)gradobj;
3316:   }
3317:   PetscFunctionReturn(PETSC_SUCCESS);
3318: }

3320: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
3321: {
3322:   PetscInt l, m;

3324:   PetscFunctionBeginHot;
3325:   if (dimC == dimR && dimR <= 3) {
3326:     /* invert Jacobian, multiply */
3327:     PetscScalar det, idet;

3329:     switch (dimR) {
3330:     case 1:
3331:       invJ[0] = 1. / J[0];
3332:       break;
3333:     case 2:
3334:       det     = J[0] * J[3] - J[1] * J[2];
3335:       idet    = 1. / det;
3336:       invJ[0] = J[3] * idet;
3337:       invJ[1] = -J[1] * idet;
3338:       invJ[2] = -J[2] * idet;
3339:       invJ[3] = J[0] * idet;
3340:       break;
3341:     case 3: {
3342:       invJ[0] = J[4] * J[8] - J[5] * J[7];
3343:       invJ[1] = J[2] * J[7] - J[1] * J[8];
3344:       invJ[2] = J[1] * J[5] - J[2] * J[4];
3345:       det     = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
3346:       idet    = 1. / det;
3347:       invJ[0] *= idet;
3348:       invJ[1] *= idet;
3349:       invJ[2] *= idet;
3350:       invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
3351:       invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
3352:       invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
3353:       invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
3354:       invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
3355:       invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
3356:     } break;
3357:     }
3358:     for (l = 0; l < dimR; l++) {
3359:       for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
3360:     }
3361:   } else {
3362: #if defined(PETSC_USE_COMPLEX)
3363:     char transpose = 'C';
3364: #else
3365:     char transpose = 'T';
3366: #endif
3367:     PetscBLASInt m        = dimR;
3368:     PetscBLASInt n        = dimC;
3369:     PetscBLASInt one      = 1;
3370:     PetscBLASInt worksize = dimR * dimC, info;

3372:     for (l = 0; l < dimC; l++) invJ[l] = resNeg[l];

3374:     PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info));
3375:     PetscCheck(info == 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELS");

3377:     for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3378:   }
3379:   PetscFunctionReturn(PETSC_SUCCESS);
3380: }

3382: static PetscErrorCode DMPlexCoordinatesToReference_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3383: {
3384:   PetscInt     coordSize, i, j, k, l, m, maxIts = 7, numV = (1 << dimR);
3385:   PetscScalar *coordsScalar = NULL;
3386:   PetscReal   *cellData, *cellCoords, *cellCoeffs, *extJ, *resNeg;
3387:   PetscScalar *J, *invJ, *work;

3389:   PetscFunctionBegin;
3391:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3392:   PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3393:   PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3394:   PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3395:   cellCoords = &cellData[0];
3396:   cellCoeffs = &cellData[coordSize];
3397:   extJ       = &cellData[2 * coordSize];
3398:   resNeg     = &cellData[2 * coordSize + dimR];
3399:   invJ       = &J[dimR * dimC];
3400:   work       = &J[2 * dimR * dimC];
3401:   if (dimR == 2) {
3402:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

3404:     for (i = 0; i < 4; i++) {
3405:       PetscInt plexI = zToPlex[i];

3407:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3408:     }
3409:   } else if (dimR == 3) {
3410:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

3412:     for (i = 0; i < 8; i++) {
3413:       PetscInt plexI = zToPlex[i];

3415:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3416:     }
3417:   } else {
3418:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3419:   }
3420:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3421:   for (i = 0; i < dimR; i++) {
3422:     PetscReal *swap;

3424:     for (j = 0; j < (numV / 2); j++) {
3425:       for (k = 0; k < dimC; k++) {
3426:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3427:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3428:       }
3429:     }

3431:     if (i < dimR - 1) {
3432:       swap       = cellCoeffs;
3433:       cellCoeffs = cellCoords;
3434:       cellCoords = swap;
3435:     }
3436:   }
3437:   PetscCall(PetscArrayzero(refCoords, numPoints * dimR));
3438:   for (j = 0; j < numPoints; j++) {
3439:     for (i = 0; i < maxIts; i++) {
3440:       PetscReal *guess = &refCoords[dimR * j];

3442:       /* compute -residual and Jacobian */
3443:       for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3444:       for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3445:       for (k = 0; k < numV; k++) {
3446:         PetscReal extCoord = 1.;
3447:         for (l = 0; l < dimR; l++) {
3448:           PetscReal coord = guess[l];
3449:           PetscInt  dep   = (k & (1 << l)) >> l;

3451:           extCoord *= dep * coord + !dep;
3452:           extJ[l] = dep;

3454:           for (m = 0; m < dimR; m++) {
3455:             PetscReal coord = guess[m];
3456:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
3457:             PetscReal mult  = dep * coord + !dep;

3459:             extJ[l] *= mult;
3460:           }
3461:         }
3462:         for (l = 0; l < dimC; l++) {
3463:           PetscReal coeff = cellCoeffs[dimC * k + l];

3465:           resNeg[l] -= coeff * extCoord;
3466:           for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3467:         }
3468:       }
3469:       if (0 && PetscDefined(USE_DEBUG)) {
3470:         PetscReal maxAbs = 0.;

3472:         for (l = 0; l < dimC; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3473:         PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3474:       }

3476:       PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess));
3477:     }
3478:   }
3479:   PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3480:   PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3481:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3482:   PetscFunctionReturn(PETSC_SUCCESS);
3483: }

3485: static PetscErrorCode DMPlexReferenceToCoordinates_Tensor(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt dimC, PetscInt dimR)
3486: {
3487:   PetscInt     coordSize, i, j, k, l, numV = (1 << dimR);
3488:   PetscScalar *coordsScalar = NULL;
3489:   PetscReal   *cellData, *cellCoords, *cellCoeffs;

3491:   PetscFunctionBegin;
3493:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3494:   PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3495:   PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3496:   cellCoords = &cellData[0];
3497:   cellCoeffs = &cellData[coordSize];
3498:   if (dimR == 2) {
3499:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

3501:     for (i = 0; i < 4; i++) {
3502:       PetscInt plexI = zToPlex[i];

3504:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3505:     }
3506:   } else if (dimR == 3) {
3507:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

3509:     for (i = 0; i < 8; i++) {
3510:       PetscInt plexI = zToPlex[i];

3512:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3513:     }
3514:   } else {
3515:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3516:   }
3517:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3518:   for (i = 0; i < dimR; i++) {
3519:     PetscReal *swap;

3521:     for (j = 0; j < (numV / 2); j++) {
3522:       for (k = 0; k < dimC; k++) {
3523:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3524:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3525:       }
3526:     }

3528:     if (i < dimR - 1) {
3529:       swap       = cellCoeffs;
3530:       cellCoeffs = cellCoords;
3531:       cellCoords = swap;
3532:     }
3533:   }
3534:   PetscCall(PetscArrayzero(realCoords, numPoints * dimC));
3535:   for (j = 0; j < numPoints; j++) {
3536:     const PetscReal *guess  = &refCoords[dimR * j];
3537:     PetscReal       *mapped = &realCoords[dimC * j];

3539:     for (k = 0; k < numV; k++) {
3540:       PetscReal extCoord = 1.;
3541:       for (l = 0; l < dimR; l++) {
3542:         PetscReal coord = guess[l];
3543:         PetscInt  dep   = (k & (1 << l)) >> l;

3545:         extCoord *= dep * coord + !dep;
3546:       }
3547:       for (l = 0; l < dimC; l++) {
3548:         PetscReal coeff = cellCoeffs[dimC * k + l];

3550:         mapped[l] += coeff * extCoord;
3551:       }
3552:     }
3553:   }
3554:   PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3555:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3556:   PetscFunctionReturn(PETSC_SUCCESS);
3557: }

3559: /* TODO: TOBY please fix this for Nc > 1 */
3560: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3561: {
3562:   PetscInt     numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
3563:   PetscScalar *nodes = NULL;
3564:   PetscReal   *invV, *modes;
3565:   PetscReal   *B, *D, *resNeg;
3566:   PetscScalar *J, *invJ, *work;

3568:   PetscFunctionBegin;
3569:   PetscCall(PetscFEGetDimension(fe, &pdim));
3570:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
3571:   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);
3572:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3573:   /* convert nodes to values in the stable evaluation basis */
3574:   PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3575:   invV = fe->invV;
3576:   for (i = 0; i < pdim; ++i) {
3577:     modes[i] = 0.;
3578:     for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3579:   }
3580:   PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3581:   D      = &B[pdim * Nc];
3582:   resNeg = &D[pdim * Nc * dimR];
3583:   PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3584:   invJ = &J[Nc * dimR];
3585:   work = &invJ[Nc * dimR];
3586:   for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.;
3587:   for (j = 0; j < numPoints; j++) {
3588:     for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
3589:       PetscReal *guess = &refCoords[j * dimR];
3590:       PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL));
3591:       for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k];
3592:       for (k = 0; k < Nc * dimR; k++) J[k] = 0.;
3593:       for (k = 0; k < pdim; k++) {
3594:         for (l = 0; l < Nc; l++) {
3595:           resNeg[l] -= modes[k] * B[k * Nc + l];
3596:           for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
3597:         }
3598:       }
3599:       if (0 && PetscDefined(USE_DEBUG)) {
3600:         PetscReal maxAbs = 0.;

3602:         for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3603:         PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3604:       }
3605:       PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess));
3606:     }
3607:   }
3608:   PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3609:   PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3610:   PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3611:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3612:   PetscFunctionReturn(PETSC_SUCCESS);
3613: }

3615: /* TODO: TOBY please fix this for Nc > 1 */
3616: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3617: {
3618:   PetscInt     numComp, pdim, i, j, k, l, coordSize;
3619:   PetscScalar *nodes = NULL;
3620:   PetscReal   *invV, *modes;
3621:   PetscReal   *B;

3623:   PetscFunctionBegin;
3624:   PetscCall(PetscFEGetDimension(fe, &pdim));
3625:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
3626:   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);
3627:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3628:   /* convert nodes to values in the stable evaluation basis */
3629:   PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3630:   invV = fe->invV;
3631:   for (i = 0; i < pdim; ++i) {
3632:     modes[i] = 0.;
3633:     for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3634:   }
3635:   PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3636:   PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL));
3637:   for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3638:   for (j = 0; j < numPoints; j++) {
3639:     PetscReal *mapped = &realCoords[j * Nc];

3641:     for (k = 0; k < pdim; k++) {
3642:       for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3643:     }
3644:   }
3645:   PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3646:   PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3647:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3648:   PetscFunctionReturn(PETSC_SUCCESS);
3649: }

3651: /*@
3652:   DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element
3653:   using a single element map.

3655:   Not Collective

3657:   Input Parameters:
3658: + dm         - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3659:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3660:                as a multilinear map for tensor-product elements
3661: . cell       - the cell whose map is used.
3662: . numPoints  - the number of points to locate
3663: - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)

3665:   Output Parameter:
3666: . refCoords - (`numPoints` x `dimension`) array of reference coordinates (see `DMGetDimension()`)

3668:   Level: intermediate

3670:   Notes:
3671:   This inversion will be accurate inside the reference element, but may be inaccurate for
3672:   mappings that do not extend uniquely outside the reference cell (e.g, most non-affine maps)

3674: .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()`
3675: @*/
3676: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3677: {
3678:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3679:   DM       coordDM = NULL;
3680:   Vec      coords;
3681:   PetscFE  fe = NULL;

3683:   PetscFunctionBegin;
3685:   PetscCall(DMGetDimension(dm, &dimR));
3686:   PetscCall(DMGetCoordinateDim(dm, &dimC));
3687:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3688:   PetscCall(DMPlexGetDepth(dm, &depth));
3689:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
3690:   PetscCall(DMGetCoordinateDM(dm, &coordDM));
3691:   if (coordDM) {
3692:     PetscInt coordFields;

3694:     PetscCall(DMGetNumFields(coordDM, &coordFields));
3695:     if (coordFields) {
3696:       PetscClassId id;
3697:       PetscObject  disc;

3699:       PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3700:       PetscCall(PetscObjectGetClassId(disc, &id));
3701:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3702:     }
3703:   }
3704:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3705:   PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd);
3706:   if (!fe) { /* implicit discretization: affine or multilinear */
3707:     PetscInt  coneSize;
3708:     PetscBool isSimplex, isTensor;

3710:     PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3711:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3712:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3713:     if (isSimplex) {
3714:       PetscReal detJ, *v0, *J, *invJ;

3716:       PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3717:       J    = &v0[dimC];
3718:       invJ = &J[dimC * dimC];
3719:       PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ));
3720:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3721:         const PetscReal x0[3] = {-1., -1., -1.};

3723:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3724:       }
3725:       PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3726:     } else if (isTensor) {
3727:       PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3728:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3729:   } else {
3730:     PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3731:   }
3732:   PetscFunctionReturn(PETSC_SUCCESS);
3733: }

3735: /*@
3736:   DMPlexReferenceToCoordinates - Map references coordinates to coordinates in the mesh for a single element map.

3738:   Not Collective

3740:   Input Parameters:
3741: + dm        - The mesh, with coordinate maps defined either by a PetscDS for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3742:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3743:                as a multilinear map for tensor-product elements
3744: . cell      - the cell whose map is used.
3745: . numPoints - the number of points to locate
3746: - refCoords - (numPoints x dimension) array of reference coordinates (see `DMGetDimension()`)

3748:   Output Parameter:
3749: . realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)

3751:   Level: intermediate

3753: .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()`
3754: @*/
3755: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3756: {
3757:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3758:   DM       coordDM = NULL;
3759:   Vec      coords;
3760:   PetscFE  fe = NULL;

3762:   PetscFunctionBegin;
3764:   PetscCall(DMGetDimension(dm, &dimR));
3765:   PetscCall(DMGetCoordinateDim(dm, &dimC));
3766:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3767:   PetscCall(DMPlexGetDepth(dm, &depth));
3768:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
3769:   PetscCall(DMGetCoordinateDM(dm, &coordDM));
3770:   if (coordDM) {
3771:     PetscInt coordFields;

3773:     PetscCall(DMGetNumFields(coordDM, &coordFields));
3774:     if (coordFields) {
3775:       PetscClassId id;
3776:       PetscObject  disc;

3778:       PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3779:       PetscCall(PetscObjectGetClassId(disc, &id));
3780:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3781:     }
3782:   }
3783:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3784:   PetscCheck(cell >= cStart && cell < cEnd, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "point %" PetscInt_FMT " not in cell range [%" PetscInt_FMT ",%" PetscInt_FMT ")", cell, cStart, cEnd);
3785:   if (!fe) { /* implicit discretization: affine or multilinear */
3786:     PetscInt  coneSize;
3787:     PetscBool isSimplex, isTensor;

3789:     PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3790:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3791:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3792:     if (isSimplex) {
3793:       PetscReal detJ, *v0, *J;

3795:       PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3796:       J = &v0[dimC];
3797:       PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ));
3798:       for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3799:         const PetscReal xi0[3] = {-1., -1., -1.};

3801:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3802:       }
3803:       PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3804:     } else if (isTensor) {
3805:       PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3806:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3807:   } else {
3808:     PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3809:   }
3810:   PetscFunctionReturn(PETSC_SUCCESS);
3811: }

3813: 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[])
3814: {
3815:   const PetscInt Nc = uOff[1] - uOff[0];
3816:   PetscInt       c;

3818:   for (c = 0; c < Nc; ++c) f0[c] = u[c];
3819: }

3821: /* Shear applies the transformation, assuming we fix z,
3822:   / 1  0  m_0 \
3823:   | 0  1  m_1 |
3824:   \ 0  0   1  /
3825: */
3826: 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[])
3827: {
3828:   const PetscInt Nc = uOff[1] - uOff[0];
3829:   const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
3830:   PetscInt       c;

3832:   for (c = 0; c < Nc; ++c) coords[c] = u[c] + constants[c + 1] * u[ax];
3833: }

3835: /* Flare applies the transformation, assuming we fix x_f,

3837:    x_i = x_i * alpha_i x_f
3838: */
3839: 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[])
3840: {
3841:   const PetscInt Nc = uOff[1] - uOff[0];
3842:   const PetscInt cf = (PetscInt)PetscRealPart(constants[0]);
3843:   PetscInt       c;

3845:   for (c = 0; c < Nc; ++c) coords[c] = u[c] * (c == cf ? 1.0 : constants[c + 1] * u[cf]);
3846: }

3848: /*
3849:   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
3850:   will correspond to the top and bottom of our square. So

3852:     (0,0)--(1,0)  ==>  (1,0)--(2,0)      Just a shift of (1,0)
3853:     (0,1)--(1,1)  ==>  (0,1)--(0,2)      Switch x and y

3855:   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:

3857:     (x, y)  ==>  (x+1, \pi/2 y)                           in (r', \theta') space
3858:             ==>  ((x+1) cos(\pi/2 y), (x+1) sin(\pi/2 y)) in (x', y') space
3859: */
3860: 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[])
3861: {
3862:   const PetscReal ri = PetscRealPart(constants[0]);
3863:   const PetscReal ro = PetscRealPart(constants[1]);

3865:   xp[0] = (x[0] * (ro - ri) + ri) * PetscCosReal(0.5 * PETSC_PI * x[1]);
3866:   xp[1] = (x[0] * (ro - ri) + ri) * PetscSinReal(0.5 * PETSC_PI * x[1]);
3867: }

3869: /*
3870:   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
3871:   lower hemisphere and the upper surface onto the top, letting z be the radius.

3873:     (x, y)  ==>  ((z+3)/2, \pi/2 (|x| or |y|), arctan y/x)                                                  in (r', \theta', \phi') space
3874:             ==>  ((z+3)/2 \cos(\theta') cos(\phi'), (z+3)/2 \cos(\theta') sin(\phi'), (z+3)/2 sin(\theta')) in (x', y', z') space
3875: */
3876: 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[])
3877: {
3878:   const PetscReal pi4    = PETSC_PI / 4.0;
3879:   const PetscReal ri     = PetscRealPart(constants[0]);
3880:   const PetscReal ro     = PetscRealPart(constants[1]);
3881:   const PetscReal rp     = (x[2] + 1) * 0.5 * (ro - ri) + ri;
3882:   const PetscReal phip   = PetscAtan2Real(x[1], x[0]);
3883:   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])));

3885:   xp[0] = rp * PetscCosReal(thetap) * PetscCosReal(phip);
3886:   xp[1] = rp * PetscCosReal(thetap) * PetscSinReal(phip);
3887:   xp[2] = rp * PetscSinReal(thetap);
3888: }

3890: /*@C
3891:   DMPlexRemapGeometry - This function maps the original `DM` coordinates to new coordinates.

3893:   Not Collective

3895:   Input Parameters:
3896: + dm   - The `DM`
3897: . time - The time
3898: - func - The function transforming current coordinates to new coordinates

3900:   Calling sequence of `func`:
3901: + dim          - The spatial dimension
3902: . Nf           - The number of input fields (here 1)
3903: . NfAux        - The number of input auxiliary fields
3904: . uOff         - The offset of the coordinates in u[] (here 0)
3905: . uOff_x       - The offset of the coordinates in u_x[] (here 0)
3906: . u            - The coordinate values at this point in space
3907: . u_t          - The coordinate time derivative at this point in space (here `NULL`)
3908: . u_x          - The coordinate derivatives at this point in space
3909: . aOff         - The offset of each auxiliary field in u[]
3910: . aOff_x       - The offset of each auxiliary field in u_x[]
3911: . a            - The auxiliary field values at this point in space
3912: . a_t          - The auxiliary field time derivative at this point in space (or `NULL`)
3913: . a_x          - The auxiliary field derivatives at this point in space
3914: . t            - The current time
3915: . x            - The coordinates of this point (here not used)
3916: . numConstants - The number of constants
3917: . constants    - The value of each constant
3918: - f            - The new coordinates at this point in space

3920:   Level: intermediate

3922: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
3923: @*/
3924: 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[]))
3925: {
3926:   DM           cdm;
3927:   PetscDS      cds;
3928:   DMField      cf;
3929:   PetscObject  obj;
3930:   PetscClassId id;
3931:   Vec          lCoords, tmpCoords;

3933:   PetscFunctionBegin;
3934:   PetscCall(DMGetCoordinateDM(dm, &cdm));
3935:   PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
3936:   PetscCall(DMGetDS(cdm, &cds));
3937:   PetscCall(PetscDSGetDiscretization(cds, 0, &obj));
3938:   PetscCall(PetscObjectGetClassId(obj, &id));
3939:   if (id != PETSCFE_CLASSID) {
3940:     PetscSection       cSection;
3941:     const PetscScalar *constants;
3942:     PetscScalar       *coords, f[16];
3943:     PetscInt           dim, cdim, Nc, vStart, vEnd;

3945:     PetscCall(DMGetDimension(dm, &dim));
3946:     PetscCall(DMGetCoordinateDim(dm, &cdim));
3947:     PetscCheck(cdim <= 16, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Affine version of DMPlexRemapGeometry is currently limited to dimensions <= 16, not %" PetscInt_FMT, cdim);
3948:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
3949:     PetscCall(DMGetCoordinateSection(dm, &cSection));
3950:     PetscCall(PetscDSGetConstants(cds, &Nc, &constants));
3951:     PetscCall(VecGetArrayWrite(lCoords, &coords));
3952:     for (PetscInt v = vStart; v < vEnd; ++v) {
3953:       PetscInt uOff[2] = {0, cdim};
3954:       PetscInt off, c;

3956:       PetscCall(PetscSectionGetOffset(cSection, v, &off));
3957:       (*func)(dim, 1, 0, uOff, NULL, &coords[off], NULL, NULL, NULL, NULL, NULL, NULL, NULL, 0.0, NULL, Nc, constants, f);
3958:       for (c = 0; c < cdim; ++c) coords[off + c] = f[c];
3959:     }
3960:     PetscCall(VecRestoreArrayWrite(lCoords, &coords));
3961:   } else {
3962:     PetscCall(DMGetLocalVector(cdm, &tmpCoords));
3963:     PetscCall(VecCopy(lCoords, tmpCoords));
3964:     /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3965:     PetscCall(DMGetCoordinateField(dm, &cf));
3966:     cdm->coordinates[0].field = cf;
3967:     PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords));
3968:     cdm->coordinates[0].field = NULL;
3969:     PetscCall(DMRestoreLocalVector(cdm, &tmpCoords));
3970:     PetscCall(DMSetCoordinatesLocal(dm, lCoords));
3971:   }
3972:   PetscFunctionReturn(PETSC_SUCCESS);
3973: }

3975: /*@C
3976:   DMPlexShearGeometry - This shears the domain, meaning adds a multiple of the shear coordinate to all other coordinates.

3978:   Not Collective

3980:   Input Parameters:
3981: + dm          - The `DMPLEX`
3982: . direction   - The shear coordinate direction, e.g. `DM_X` is the x-axis
3983: - multipliers - The multiplier m for each direction which is not the shear direction

3985:   Level: intermediate

3987: .seealso: `DMPLEX`, `DMPlexRemapGeometry()`, `DMDirection`, `DM_X`, `DM_Y`, `DM_Z`
3988: @*/
3989: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3990: {
3991:   DM             cdm;
3992:   PetscDS        cds;
3993:   PetscScalar   *moduli;
3994:   const PetscInt dir = (PetscInt)direction;
3995:   PetscInt       dE, d, e;

3997:   PetscFunctionBegin;
3998:   PetscCall(DMGetCoordinateDM(dm, &cdm));
3999:   PetscCall(DMGetCoordinateDim(dm, &dE));
4000:   PetscCall(PetscMalloc1(dE + 1, &moduli));
4001:   moduli[0] = dir;
4002:   for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
4003:   PetscCall(DMGetDS(cdm, &cds));
4004:   PetscCall(PetscDSSetConstants(cds, dE + 1, moduli));
4005:   PetscCall(DMPlexRemapGeometry(dm, 0.0, coordMap_shear));
4006:   PetscCall(PetscFree(moduli));
4007:   PetscFunctionReturn(PETSC_SUCCESS);
4008: }