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: /*@
  7:   DMPlexFindVertices - Try to find DAG points based on their coordinates.

  9:   Not Collective (provided `DMGetCoordinatesLocalSetUp()` has been already called)

 11:   Input Parameters:
 12: + dm - The `DMPLEX` object
 13: . coordinates - The `Vec` of coordinates of the sought points
 14: - eps - The tolerance or `PETSC_DEFAULT`

 16:   Output Parameter:
 17: . points - The `IS` of found DAG points or -1

 19:   Level: intermediate

 21:   Notes:
 22:   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.

 24:   The output `IS` is living on `PETSC_COMM_SELF` and its length is npoints.
 25:   Each rank does the search independently.
 26:   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.

 28:   The output `IS` must be destroyed by user.

 30:   The tolerance is interpreted as the maximum Euclidean (L2) distance of the sought point from the specified coordinates.

 32:   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.

 34: .seealso: `DMPLEX`, `DMPlexCreate()`, `DMGetCoordinatesLocal()`
 35: @*/
 36: PetscErrorCode DMPlexFindVertices(DM dm, Vec coordinates, PetscReal eps, IS *points)
 37: {
 38:   PetscInt           c, cdim, i, j, o, p, vStart, vEnd;
 39:   PetscInt           npoints;
 40:   const PetscScalar *coord;
 41:   Vec                allCoordsVec;
 42:   const PetscScalar *allCoords;
 43:   PetscInt          *dagPoints;

 45:   PetscFunctionBegin;
 46:   if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
 47:   PetscCall(DMGetCoordinateDim(dm, &cdim));
 48:   {
 49:     PetscInt n;

 51:     PetscCall(VecGetLocalSize(coordinates, &n));
 52:     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);
 53:     npoints = n / cdim;
 54:   }
 55:   PetscCall(DMGetCoordinatesLocal(dm, &allCoordsVec));
 56:   PetscCall(VecGetArrayRead(allCoordsVec, &allCoords));
 57:   PetscCall(VecGetArrayRead(coordinates, &coord));
 58:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
 59:   if (PetscDefined(USE_DEBUG)) {
 60:     /* check coordinate section is consistent with DM dimension */
 61:     PetscSection cs;
 62:     PetscInt     ndof;

 64:     PetscCall(DMGetCoordinateSection(dm, &cs));
 65:     for (p = vStart; p < vEnd; p++) {
 66:       PetscCall(PetscSectionGetDof(cs, p, &ndof));
 67:       PetscCheck(ndof == cdim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %" PetscInt_FMT ": ndof = %" PetscInt_FMT " != %" PetscInt_FMT " = cdim", p, ndof, cdim);
 68:     }
 69:   }
 70:   PetscCall(PetscMalloc1(npoints, &dagPoints));
 71:   if (eps == 0.0) {
 72:     for (i = 0, j = 0; i < npoints; i++, j += cdim) {
 73:       dagPoints[i] = -1;
 74:       for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
 75:         for (c = 0; c < cdim; c++) {
 76:           if (coord[j + c] != allCoords[o + c]) break;
 77:         }
 78:         if (c == cdim) {
 79:           dagPoints[i] = p;
 80:           break;
 81:         }
 82:       }
 83:     }
 84:   } else {
 85:     for (i = 0, j = 0; i < npoints; i++, j += cdim) {
 86:       PetscReal norm;

 88:       dagPoints[i] = -1;
 89:       for (p = vStart, o = 0; p < vEnd; p++, o += cdim) {
 90:         norm = 0.0;
 91:         for (c = 0; c < cdim; c++) norm += PetscRealPart(PetscSqr(coord[j + c] - allCoords[o + c]));
 92:         norm = PetscSqrtReal(norm);
 93:         if (norm <= eps) {
 94:           dagPoints[i] = p;
 95:           break;
 96:         }
 97:       }
 98:     }
 99:   }
100:   PetscCall(VecRestoreArrayRead(allCoordsVec, &allCoords));
101:   PetscCall(VecRestoreArrayRead(coordinates, &coord));
102:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, npoints, dagPoints, PETSC_OWN_POINTER, points));
103:   PetscFunctionReturn(PETSC_SUCCESS);
104: }

106: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
107: {
108:   const PetscReal p0_x  = segmentA[0 * 2 + 0];
109:   const PetscReal p0_y  = segmentA[0 * 2 + 1];
110:   const PetscReal p1_x  = segmentA[1 * 2 + 0];
111:   const PetscReal p1_y  = segmentA[1 * 2 + 1];
112:   const PetscReal p2_x  = segmentB[0 * 2 + 0];
113:   const PetscReal p2_y  = segmentB[0 * 2 + 1];
114:   const PetscReal p3_x  = segmentB[1 * 2 + 0];
115:   const PetscReal p3_y  = segmentB[1 * 2 + 1];
116:   const PetscReal s1_x  = p1_x - p0_x;
117:   const PetscReal s1_y  = p1_y - p0_y;
118:   const PetscReal s2_x  = p3_x - p2_x;
119:   const PetscReal s2_y  = p3_y - p2_y;
120:   const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

122:   PetscFunctionBegin;
123:   *hasIntersection = PETSC_FALSE;
124:   /* Non-parallel lines */
125:   if (denom != 0.0) {
126:     const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
127:     const PetscReal t = (s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;

129:     if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
130:       *hasIntersection = PETSC_TRUE;
131:       if (intersection) {
132:         intersection[0] = p0_x + (t * s1_x);
133:         intersection[1] = p0_y + (t * s1_y);
134:       }
135:     }
136:   }
137:   PetscFunctionReturn(PETSC_SUCCESS);
138: }

140: /* The plane is segmentB x segmentC: https://en.wikipedia.org/wiki/Line%E2%80%93plane_intersection */
141: static PetscErrorCode DMPlexGetLinePlaneIntersection_3D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], const PetscReal segmentC[], PetscReal intersection[], PetscBool *hasIntersection)
142: {
143:   const PetscReal p0_x  = segmentA[0 * 3 + 0];
144:   const PetscReal p0_y  = segmentA[0 * 3 + 1];
145:   const PetscReal p0_z  = segmentA[0 * 3 + 2];
146:   const PetscReal p1_x  = segmentA[1 * 3 + 0];
147:   const PetscReal p1_y  = segmentA[1 * 3 + 1];
148:   const PetscReal p1_z  = segmentA[1 * 3 + 2];
149:   const PetscReal q0_x  = segmentB[0 * 3 + 0];
150:   const PetscReal q0_y  = segmentB[0 * 3 + 1];
151:   const PetscReal q0_z  = segmentB[0 * 3 + 2];
152:   const PetscReal q1_x  = segmentB[1 * 3 + 0];
153:   const PetscReal q1_y  = segmentB[1 * 3 + 1];
154:   const PetscReal q1_z  = segmentB[1 * 3 + 2];
155:   const PetscReal r0_x  = segmentC[0 * 3 + 0];
156:   const PetscReal r0_y  = segmentC[0 * 3 + 1];
157:   const PetscReal r0_z  = segmentC[0 * 3 + 2];
158:   const PetscReal r1_x  = segmentC[1 * 3 + 0];
159:   const PetscReal r1_y  = segmentC[1 * 3 + 1];
160:   const PetscReal r1_z  = segmentC[1 * 3 + 2];
161:   const PetscReal s0_x  = p1_x - p0_x;
162:   const PetscReal s0_y  = p1_y - p0_y;
163:   const PetscReal s0_z  = p1_z - p0_z;
164:   const PetscReal s1_x  = q1_x - q0_x;
165:   const PetscReal s1_y  = q1_y - q0_y;
166:   const PetscReal s1_z  = q1_z - q0_z;
167:   const PetscReal s2_x  = r1_x - r0_x;
168:   const PetscReal s2_y  = r1_y - r0_y;
169:   const PetscReal s2_z  = r1_z - r0_z;
170:   const PetscReal s3_x  = s1_y * s2_z - s1_z * s2_y; /* s1 x s2 */
171:   const PetscReal s3_y  = s1_z * s2_x - s1_x * s2_z;
172:   const PetscReal s3_z  = s1_x * s2_y - s1_y * s2_x;
173:   const PetscReal s4_x  = s0_y * s2_z - s0_z * s2_y; /* s0 x s2 */
174:   const PetscReal s4_y  = s0_z * s2_x - s0_x * s2_z;
175:   const PetscReal s4_z  = s0_x * s2_y - s0_y * s2_x;
176:   const PetscReal s5_x  = s1_y * s0_z - s1_z * s0_y; /* s1 x s0 */
177:   const PetscReal s5_y  = s1_z * s0_x - s1_x * s0_z;
178:   const PetscReal s5_z  = s1_x * s0_y - s1_y * s0_x;
179:   const PetscReal denom = -(s0_x * s3_x + s0_y * s3_y + s0_z * s3_z); /* -s0 . (s1 x s2) */

181:   PetscFunctionBegin;
182:   *hasIntersection = PETSC_FALSE;
183:   /* Line not parallel to plane */
184:   if (denom != 0.0) {
185:     const PetscReal t = (s3_x * (p0_x - q0_x) + s3_y * (p0_y - q0_y) + s3_z * (p0_z - q0_z)) / denom;
186:     const PetscReal u = (s4_x * (p0_x - q0_x) + s4_y * (p0_y - q0_y) + s4_z * (p0_z - q0_z)) / denom;
187:     const PetscReal v = (s5_x * (p0_x - q0_x) + s5_y * (p0_y - q0_y) + s5_z * (p0_z - q0_z)) / denom;

189:     if (t >= 0 && t <= 1 && u >= 0 && u <= 1 && v >= 0 && v <= 1) {
190:       *hasIntersection = PETSC_TRUE;
191:       if (intersection) {
192:         intersection[0] = p0_x + (t * s0_x);
193:         intersection[1] = p0_y + (t * s0_y);
194:         intersection[2] = p0_z + (t * s0_z);
195:       }
196:     }
197:   }
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: static PetscErrorCode DMPlexLocatePoint_Simplex_1D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
202: {
203:   const PetscReal eps = PETSC_SQRT_MACHINE_EPSILON;
204:   const PetscReal x   = PetscRealPart(point[0]);
205:   PetscReal       v0, J, invJ, detJ;
206:   PetscReal       xi;

208:   PetscFunctionBegin;
209:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ));
210:   xi = invJ * (x - v0);

212:   if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c;
213:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
218: {
219:   const PetscInt  embedDim = 2;
220:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
221:   PetscReal       x        = PetscRealPart(point[0]);
222:   PetscReal       y        = PetscRealPart(point[1]);
223:   PetscReal       v0[2], J[4], invJ[4], detJ;
224:   PetscReal       xi, eta;

226:   PetscFunctionBegin;
227:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
228:   xi  = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
229:   eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);

231:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0 + eps)) *cell = c;
232:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
237: {
238:   const PetscInt embedDim = 2;
239:   PetscReal      x        = PetscRealPart(point[0]);
240:   PetscReal      y        = PetscRealPart(point[1]);
241:   PetscReal      v0[2], J[4], invJ[4], detJ;
242:   PetscReal      xi, eta, r;

244:   PetscFunctionBegin;
245:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
246:   xi  = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
247:   eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);

249:   xi  = PetscMax(xi, 0.0);
250:   eta = PetscMax(eta, 0.0);
251:   if (xi + eta > 2.0) {
252:     r = (xi + eta) / 2.0;
253:     xi /= r;
254:     eta /= r;
255:   }
256:   cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0];
257:   cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1];
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: // This is the ray-casting, or even-odd algorithm: https://en.wikipedia.org/wiki/Even%E2%80%93odd_rule
262: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
263: {
264:   const PetscScalar *array;
265:   PetscScalar       *coords    = NULL;
266:   const PetscInt     faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
267:   PetscReal          x         = PetscRealPart(point[0]);
268:   PetscReal          y         = PetscRealPart(point[1]);
269:   PetscInt           crossings = 0, numCoords, f;
270:   PetscBool          isDG;

272:   PetscFunctionBegin;
273:   PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
274:   PetscCheck(numCoords == 8, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
275:   for (f = 0; f < 4; ++f) {
276:     PetscReal x_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 0]);
277:     PetscReal y_i = PetscRealPart(coords[faces[2 * f + 0] * 2 + 1]);
278:     PetscReal x_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 0]);
279:     PetscReal y_j = PetscRealPart(coords[faces[2 * f + 1] * 2 + 1]);

281:     if ((x == x_j) && (y == y_j)) {
282:       // point is a corner
283:       crossings = 1;
284:       break;
285:     }
286:     if ((y_j > y) != (y_i > y)) {
287:       PetscReal slope = (x - x_j) * (y_i - y_j) - (x_i - x_j) * (y - y_j);
288:       if (slope == 0) {
289:         // point is a corner
290:         crossings = 1;
291:         break;
292:       }
293:       if ((slope < 0) != (y_i < y_j)) ++crossings;
294:     }
295:   }
296:   if (crossings % 2) *cell = c;
297:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
298:   PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

302: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
303: {
304:   const PetscInt  embedDim = 3;
305:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
306:   PetscReal       v0[3], J[9], invJ[9], detJ;
307:   PetscReal       x = PetscRealPart(point[0]);
308:   PetscReal       y = PetscRealPart(point[1]);
309:   PetscReal       z = PetscRealPart(point[2]);
310:   PetscReal       xi, eta, zeta;

312:   PetscFunctionBegin;
313:   PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
314:   xi   = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
315:   eta  = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
316:   zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);

318:   if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
319:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
320:   PetscFunctionReturn(PETSC_SUCCESS);
321: }

323: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
324: {
325:   const PetscScalar *array;
326:   PetscScalar       *coords    = NULL;
327:   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};
328:   PetscBool          found     = PETSC_TRUE;
329:   PetscInt           numCoords, f;
330:   PetscBool          isDG;

332:   PetscFunctionBegin;
333:   PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
334:   PetscCheck(numCoords == 24, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Quadrilateral should have 8 coordinates, not %" PetscInt_FMT, numCoords);
335:   for (f = 0; f < 6; ++f) {
336:     /* Check the point is under plane */
337:     /*   Get face normal */
338:     PetscReal v_i[3];
339:     PetscReal v_j[3];
340:     PetscReal normal[3];
341:     PetscReal pp[3];
342:     PetscReal dot;

344:     v_i[0]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
345:     v_i[1]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
346:     v_i[2]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
347:     v_j[0]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
348:     v_j[1]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
349:     v_j[2]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
350:     normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
351:     normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
352:     normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
353:     pp[0]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
354:     pp[1]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
355:     pp[2]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
356:     dot       = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];

358:     /* Check that projected point is in face (2D location problem) */
359:     if (dot < 0.0) {
360:       found = PETSC_FALSE;
361:       break;
362:     }
363:   }
364:   if (found) *cell = c;
365:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
366:   PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords));
367:   PetscFunctionReturn(PETSC_SUCCESS);
368: }

370: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
371: {
372:   PetscInt d;

374:   PetscFunctionBegin;
375:   box->dim = dim;
376:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = point ? PetscRealPart(point[d]) : 0.;
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

380: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
381: {
382:   PetscFunctionBegin;
383:   PetscCall(PetscCalloc1(1, box));
384:   PetscCall(PetscGridHashInitialize_Internal(*box, dim, point));
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
389: {
390:   PetscInt d;

392:   PetscFunctionBegin;
393:   for (d = 0; d < box->dim; ++d) {
394:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
395:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
396:   }
397:   PetscFunctionReturn(PETSC_SUCCESS);
398: }

400: /*
401:   PetscGridHashSetGrid - Divide the grid into boxes

403:   Not Collective

405:   Input Parameters:
406: + box - The grid hash object
407: . n   - The number of boxes in each dimension, or `PETSC_DETERMINE`
408: - h   - The box size in each dimension, only used if n[d] == `PETSC_DETERMINE`

410:   Level: developer

412: .seealso: `DMPLEX`, `PetscGridHashCreate()`
413: */
414: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
415: {
416:   PetscInt d;

418:   PetscFunctionBegin;
419:   for (d = 0; d < box->dim; ++d) {
420:     box->extent[d] = box->upper[d] - box->lower[d];
421:     if (n[d] == PETSC_DETERMINE) {
422:       box->h[d] = h[d];
423:       box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
424:     } else {
425:       box->n[d] = n[d];
426:       box->h[d] = box->extent[d] / n[d];
427:     }
428:   }
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }

432: /*
433:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

435:   Not Collective

437:   Input Parameters:
438: + box       - The grid hash object
439: . numPoints - The number of input points
440: - points    - The input point coordinates

442:   Output Parameters:
443: + dboxes    - An array of numPoints*dim integers expressing the enclosing box as (i_0, i_1, ..., i_dim)
444: - boxes     - An array of numPoints integers expressing the enclosing box as single number, or NULL

446:   Level: developer

448:   Note:
449:   This only guarantees that a box contains a point, not that a cell does.

451: .seealso: `DMPLEX`, `PetscGridHashCreate()`
452: */
453: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
454: {
455:   const PetscReal *lower = box->lower;
456:   const PetscReal *upper = box->upper;
457:   const PetscReal *h     = box->h;
458:   const PetscInt  *n     = box->n;
459:   const PetscInt   dim   = box->dim;
460:   PetscInt         d, p;

462:   PetscFunctionBegin;
463:   for (p = 0; p < numPoints; ++p) {
464:     for (d = 0; d < dim; ++d) {
465:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

467:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
468:       if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
469:       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);
470:       dboxes[p * dim + d] = dbox;
471:     }
472:     if (boxes)
473:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
474:   }
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

478: /*
479:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

481:  Not Collective

483:   Input Parameters:
484: + box         - The grid hash object
485: . cellSection - The PetscSection mapping cells to boxes
486: . numPoints   - The number of input points
487: - points      - The input point coordinates

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

494:   Level: developer

496:   Note:
497:   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.

499: .seealso: `DMPLEX`, `PetscGridHashGetEnclosingBox()`
500: */
501: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
502: {
503:   const PetscReal *lower = box->lower;
504:   const PetscReal *upper = box->upper;
505:   const PetscReal *h     = box->h;
506:   const PetscInt  *n     = box->n;
507:   const PetscInt   dim   = box->dim;
508:   PetscInt         bStart, bEnd, d, p;

510:   PetscFunctionBegin;
512:   *found = PETSC_FALSE;
513:   PetscCall(PetscSectionGetChart(box->cellSection, &bStart, &bEnd));
514:   for (p = 0; p < numPoints; ++p) {
515:     for (d = 0; d < dim; ++d) {
516:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

518:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
519:       if (dbox < 0 || dbox >= n[d]) PetscFunctionReturn(PETSC_SUCCESS);
520:       dboxes[p * dim + d] = dbox;
521:     }
522:     if (boxes)
523:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
524:     // It is possible for a box to overlap no grid cells
525:     if (boxes[p] < bStart || boxes[p] >= bEnd) PetscFunctionReturn(PETSC_SUCCESS);
526:   }
527:   *found = PETSC_TRUE;
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
532: {
533:   PetscFunctionBegin;
534:   if (*box) {
535:     PetscCall(PetscSectionDestroy(&(*box)->cellSection));
536:     PetscCall(ISDestroy(&(*box)->cells));
537:     PetscCall(DMLabelDestroy(&(*box)->cellsSparse));
538:   }
539:   PetscCall(PetscFree(*box));
540:   PetscFunctionReturn(PETSC_SUCCESS);
541: }

543: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
544: {
545:   DMPolytopeType ct;

547:   PetscFunctionBegin;
548:   PetscCall(DMPlexGetCellType(dm, cellStart, &ct));
549:   switch (ct) {
550:   case DM_POLYTOPE_SEGMENT:
551:     PetscCall(DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell));
552:     break;
553:   case DM_POLYTOPE_TRIANGLE:
554:     PetscCall(DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell));
555:     break;
556:   case DM_POLYTOPE_QUADRILATERAL:
557:     PetscCall(DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell));
558:     break;
559:   case DM_POLYTOPE_TETRAHEDRON:
560:     PetscCall(DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell));
561:     break;
562:   case DM_POLYTOPE_HEXAHEDRON:
563:     PetscCall(DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell));
564:     break;
565:   default:
566:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
567:   }
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

571: /*
572:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
573: */
574: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
575: {
576:   DMPolytopeType ct;

578:   PetscFunctionBegin;
579:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
580:   switch (ct) {
581:   case DM_POLYTOPE_TRIANGLE:
582:     PetscCall(DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint));
583:     break;
584: #if 0
585:     case DM_POLYTOPE_QUADRILATERAL:
586:     PetscCall(DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint));break;
587:     case DM_POLYTOPE_TETRAHEDRON:
588:     PetscCall(DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint));break;
589:     case DM_POLYTOPE_HEXAHEDRON:
590:     PetscCall(DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint));break;
591: #endif
592:   default:
593:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
594:   }
595:   PetscFunctionReturn(PETSC_SUCCESS);
596: }

598: /*
599:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the `DMPLEX`

601:   Collective

603:   Input Parameter:
604: . dm - The `DMPLEX`

606:   Output Parameter:
607: . localBox - The grid hash object

609:   Level: developer

611: .seealso: `DMPLEX`, `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
612: */
613: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
614: {
615:   PetscInt           debug = ((DM_Plex *)dm->data)->printLocate;
616:   MPI_Comm           comm;
617:   PetscGridHash      lbox;
618:   PetscSF            sf;
619:   Vec                coordinates;
620:   PetscSection       coordSection;
621:   Vec                coordsLocal;
622:   const PetscScalar *coords;
623:   PetscScalar       *edgeCoords;
624:   PetscInt          *dboxes, *boxes;
625:   const PetscInt    *leaves;
626:   PetscInt           n[3] = {2, 2, 2};
627:   PetscInt           dim, N, Nl = 0, maxConeSize, cStart, cEnd, c, eStart, eEnd, i;
628:   PetscBool          flg;

630:   PetscFunctionBegin;
631:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
632:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
633:   PetscCall(DMGetCoordinateDim(dm, &dim));
634:   PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, NULL));
635:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
636:   PetscCall(VecGetLocalSize(coordinates, &N));
637:   PetscCall(VecGetArrayRead(coordinates, &coords));
638:   PetscCall(PetscGridHashCreate(comm, dim, coords, &lbox));
639:   for (i = 0; i < N; i += dim) PetscCall(PetscGridHashEnlarge(lbox, &coords[i]));
640:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
641:   c = dim;
642:   PetscCall(PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &c, &flg));
643:   if (flg) {
644:     for (i = c; i < dim; ++i) n[i] = n[c - 1];
645:   } else {
646:     for (i = 0; i < dim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / dim) * 0.8));
647:   }
648:   PetscCall(PetscGridHashSetGrid(lbox, n, NULL));
649:   if (debug)
650:     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], (double)lbox->lower[2], (double)lbox->upper[0],
651:                           (double)lbox->upper[1], (double)lbox->upper[2], n[0], n[1], n[2], (double)lbox->h[0], (double)lbox->h[1], (double)lbox->h[2]));
652: #if 0
653:   /* Could define a custom reduction to merge these */
654:   PetscCall(MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm));
655:   PetscCall(MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm));
656: #endif
657:   /* Is there a reason to snap the local bounding box to a division of the global box? */
658:   /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
659:   /* Create label */
660:   PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
661:   if (dim < 2) eStart = eEnd = -1;
662:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse));
663:   PetscCall(DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd));
664:   /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
665:   PetscCall(DMGetCoordinatesLocal(dm, &coordsLocal));
666:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
667:   PetscCall(DMGetPointSF(dm, &sf));
668:   if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
669:   Nl = PetscMax(Nl, 0);
670:   PetscCall(PetscCalloc3(16 * dim, &dboxes, 16, &boxes, PetscPowInt(maxConeSize, dim) * dim, &edgeCoords));
671:   for (c = cStart; c < cEnd; ++c) {
672:     const PetscReal *h       = lbox->h;
673:     PetscScalar     *ccoords = NULL;
674:     PetscInt         csize   = 0;
675:     PetscInt        *closure = NULL;
676:     PetscInt         Ncl, cl, Ne = 0, idx;
677:     PetscScalar      point[3];
678:     PetscInt         dlim[6], d, e, i, j, k;

680:     PetscCall(PetscFindInt(c, Nl, leaves, &idx));
681:     if (idx >= 0) continue;
682:     /* Get all edges in cell */
683:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure));
684:     for (cl = 0; cl < Ncl * 2; ++cl) {
685:       if ((closure[cl] >= eStart) && (closure[cl] < eEnd)) {
686:         PetscScalar *ecoords = &edgeCoords[Ne * dim * 2];
687:         PetscInt     ecsize  = dim * 2;

689:         PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, closure[cl], &ecsize, &ecoords));
690:         PetscCheck(ecsize == dim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Got %" PetscInt_FMT " coords for edge, instead of %" PetscInt_FMT, ecsize, dim * 2);
691:         ++Ne;
692:       }
693:     }
694:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure));
695:     /* Find boxes enclosing each vertex */
696:     PetscCall(DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords));
697:     PetscCall(PetscGridHashGetEnclosingBox(lbox, csize / dim, ccoords, dboxes, boxes));
698:     /* Mark cells containing the vertices */
699:     for (e = 0; e < csize / dim; ++e) {
700:       if (debug)
701:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "Cell %" PetscInt_FMT " has vertex (%g, %g, %g) in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", c, (double)PetscRealPart(ccoords[e * dim + 0]), dim > 1 ? (double)PetscRealPart(ccoords[e * dim + 1]) : 0., dim > 2 ? (double)PetscRealPart(ccoords[e * dim + 2]) : 0., boxes[e], dboxes[e * dim + 0], dim > 1 ? dboxes[e * dim + 1] : -1, dim > 2 ? dboxes[e * dim + 2] : -1));
702:       PetscCall(DMLabelSetValue(lbox->cellsSparse, c, boxes[e]));
703:     }
704:     /* Get grid of boxes containing these */
705:     for (d = 0; d < dim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
706:     for (d = dim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
707:     for (e = 1; e < dim + 1; ++e) {
708:       for (d = 0; d < dim; ++d) {
709:         dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * dim + d]);
710:         dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * dim + d]);
711:       }
712:     }
713:     /* Check for intersection of box with cell */
714:     for (k = dlim[2 * 2 + 0], point[2] = lbox->lower[2] + k * h[2]; k <= dlim[2 * 2 + 1]; ++k, point[2] += h[2]) {
715:       for (j = dlim[1 * 2 + 0], point[1] = lbox->lower[1] + j * h[1]; j <= dlim[1 * 2 + 1]; ++j, point[1] += h[1]) {
716:         for (i = dlim[0 * 2 + 0], point[0] = lbox->lower[0] + i * h[0]; i <= dlim[0 * 2 + 1]; ++i, point[0] += h[0]) {
717:           const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i;
718:           PetscScalar    cpoint[3];
719:           PetscInt       cell, edge, ii, jj, kk;

721:           if (debug)
722:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "Box %" PetscInt_FMT ": (%.2g, %.2g, %.2g) -- (%.2g, %.2g, %.2g)\n", box, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), (double)PetscRealPart(point[2]), (double)PetscRealPart(point[0] + h[0]), (double)PetscRealPart(point[1] + h[1]), (double)PetscRealPart(point[2] + h[2])));
723:           /* Check whether cell contains any vertex of this subbox TODO vectorize this */
724:           for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
725:             for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
726:               for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
727:                 PetscCall(DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell));
728:                 if (cell >= 0) {
729:                   if (debug)
730:                     PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " contains vertex (%.2g, %.2g, %.2g) of box %" PetscInt_FMT "\n", c, (double)PetscRealPart(cpoint[0]), (double)PetscRealPart(cpoint[1]), (double)PetscRealPart(cpoint[2]), box));
731:                   PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
732:                   jj = kk = 2;
733:                   break;
734:                 }
735:               }
736:             }
737:           }
738:           /* Check whether cell edge intersects any face of these subboxes TODO vectorize this */
739:           for (edge = 0; edge < Ne; ++edge) {
740:             PetscReal segA[6] = {0., 0., 0., 0., 0., 0.};
741:             PetscReal segB[6] = {0., 0., 0., 0., 0., 0.};
742:             PetscReal segC[6] = {0., 0., 0., 0., 0., 0.};

744:             PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected dim %" PetscInt_FMT " > 3", dim);
745:             for (d = 0; d < dim * 2; ++d) segA[d] = PetscRealPart(edgeCoords[edge * dim * 2 + d]);
746:             /* 1D: (x) -- (x+h)               0 -- 1
747:                2D: (x,   y)   -- (x,   y+h)   (0, 0) -- (0, 1)
748:                    (x+h, y)   -- (x+h, y+h)   (1, 0) -- (1, 1)
749:                    (x,   y)   -- (x+h, y)     (0, 0) -- (1, 0)
750:                    (x,   y+h) -- (x+h, y+h)   (0, 1) -- (1, 1)
751:                3D: (x,   y,   z)   -- (x,   y+h, z),   (x,   y,   z)   -- (x,   y,   z+h) (0, 0, 0) -- (0, 1, 0), (0, 0, 0) -- (0, 0, 1)
752:                    (x+h, y,   z)   -- (x+h, y+h, z),   (x+h, y,   z)   -- (x+h, y,   z+h) (1, 0, 0) -- (1, 1, 0), (1, 0, 0) -- (1, 0, 1)
753:                    (x,   y,   z)   -- (x+h, y,   z),   (x,   y,   z)   -- (x,   y,   z+h) (0, 0, 0) -- (1, 0, 0), (0, 0, 0) -- (0, 0, 1)
754:                    (x,   y+h, z)   -- (x+h, y+h, z),   (x,   y+h, z)   -- (x,   y+h, z+h) (0, 1, 0) -- (1, 1, 0), (0, 1, 0) -- (0, 1, 1)
755:                    (x,   y,   z)   -- (x+h, y,   z),   (x,   y,   z)   -- (x,   y+h, z)   (0, 0, 0) -- (1, 0, 0), (0, 0, 0) -- (0, 1, 0)
756:                    (x,   y,   z+h) -- (x+h, y,   z+h), (x,   y,   z+h) -- (x,   y+h, z+h) (0, 0, 1) -- (1, 0, 1), (0, 0, 1) -- (0, 1, 1)
757:              */
758:             /* Loop over faces with normal in direction d */
759:             for (d = 0; d < dim; ++d) {
760:               PetscBool intersects = PETSC_FALSE;
761:               PetscInt  e          = (d + 1) % dim;
762:               PetscInt  f          = (d + 2) % dim;

764:               /* There are two faces in each dimension */
765:               for (ii = 0; ii < 2; ++ii) {
766:                 segB[d]       = PetscRealPart(point[d] + ii * h[d]);
767:                 segB[dim + d] = PetscRealPart(point[d] + ii * h[d]);
768:                 segC[d]       = PetscRealPart(point[d] + ii * h[d]);
769:                 segC[dim + d] = PetscRealPart(point[d] + ii * h[d]);
770:                 if (dim > 1) {
771:                   segB[e]       = PetscRealPart(point[e] + 0 * h[e]);
772:                   segB[dim + e] = PetscRealPart(point[e] + 1 * h[e]);
773:                   segC[e]       = PetscRealPart(point[e] + 0 * h[e]);
774:                   segC[dim + e] = PetscRealPart(point[e] + 0 * h[e]);
775:                 }
776:                 if (dim > 2) {
777:                   segB[f]       = PetscRealPart(point[f] + 0 * h[f]);
778:                   segB[dim + f] = PetscRealPart(point[f] + 0 * h[f]);
779:                   segC[f]       = PetscRealPart(point[f] + 0 * h[f]);
780:                   segC[dim + f] = PetscRealPart(point[f] + 1 * h[f]);
781:                 }
782:                 if (dim == 2) {
783:                   PetscCall(DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects));
784:                 } else if (dim == 3) {
785:                   PetscCall(DMPlexGetLinePlaneIntersection_3D_Internal(segA, segB, segC, NULL, &intersects));
786:                 }
787:                 if (intersects) {
788:                   if (debug)
789:                     PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Cell %" PetscInt_FMT " edge %" PetscInt_FMT " (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g) intersects box %" PetscInt_FMT ", face (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g) (%.2g, %.2g, %.2g)--(%.2g, %.2g, %.2g)\n", c, edge, (double)segA[0], (double)segA[1], (double)segA[2], (double)segA[3], (double)segA[4], (double)segA[5], box, (double)segB[0], (double)segB[1], (double)segB[2], (double)segB[3], (double)segB[4], (double)segB[5], (double)segC[0], (double)segC[1], (double)segC[2], (double)segC[3], (double)segC[4], (double)segC[5]));
790:                   PetscCall(DMLabelSetValue(lbox->cellsSparse, c, box));
791:                   edge = Ne;
792:                   break;
793:                 }
794:               }
795:             }
796:           }
797:         }
798:       }
799:     }
800:     PetscCall(DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords));
801:   }
802:   PetscCall(PetscFree3(dboxes, boxes, edgeCoords));
803:   if (debug) PetscCall(DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF));
804:   PetscCall(DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells));
805:   PetscCall(DMLabelDestroy(&lbox->cellsSparse));
806:   *localBox = lbox;
807:   PetscFunctionReturn(PETSC_SUCCESS);
808: }

810: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
811: {
812:   PetscInt        debug = ((DM_Plex *)dm->data)->printLocate;
813:   DM_Plex        *mesh  = (DM_Plex *)dm->data;
814:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
815:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
816:   PetscInt        dim, Nl = 0, cStart, cEnd, numCells, c, d;
817:   PetscSF         sf;
818:   const PetscInt *leaves;
819:   const PetscInt *boxCells;
820:   PetscSFNode    *cells;
821:   PetscScalar    *a;
822:   PetscMPIInt     result;
823:   PetscLogDouble  t0, t1;
824:   PetscReal       gmin[3], gmax[3];
825:   PetscInt        terminating_query_type[] = {0, 0, 0};

827:   PetscFunctionBegin;
828:   PetscCall(PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0));
829:   PetscCall(PetscTime(&t0));
830:   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.");
831:   PetscCall(DMGetCoordinateDim(dm, &dim));
832:   PetscCall(VecGetBlockSize(v, &bs));
833:   PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result));
834:   PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PetscObjectComm((PetscObject)cellSF), PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
835:   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);
836:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
837:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
838:   PetscCall(DMGetPointSF(dm, &sf));
839:   if (sf) PetscCall(PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL));
840:   Nl = PetscMax(Nl, 0);
841:   PetscCall(VecGetLocalSize(v, &numPoints));
842:   PetscCall(VecGetArray(v, &a));
843:   numPoints /= bs;
844:   {
845:     const PetscSFNode *sf_cells;

847:     PetscCall(PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells));
848:     if (sf_cells) {
849:       PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n"));
850:       cells = (PetscSFNode *)sf_cells;
851:       reuse = PETSC_TRUE;
852:     } else {
853:       PetscCall(PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n"));
854:       PetscCall(PetscMalloc1(numPoints, &cells));
855:       /* initialize cells if created */
856:       for (p = 0; p < numPoints; p++) {
857:         cells[p].rank  = 0;
858:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
859:       }
860:     }
861:   }
862:   PetscCall(DMGetBoundingBox(dm, gmin, gmax));
863:   if (hash) {
864:     if (!mesh->lbox) {
865:       PetscCall(PetscInfo(dm, "Initializing grid hashing\n"));
866:       PetscCall(DMPlexComputeGridHash_Internal(dm, &mesh->lbox));
867:     }
868:     /* Designate the local box for each point */
869:     /* Send points to correct process */
870:     /* Search cells that lie in each subbox */
871:     /*   Should we bin points before doing search? */
872:     PetscCall(ISGetIndices(mesh->lbox->cells, &boxCells));
873:   }
874:   for (p = 0, numFound = 0; p < numPoints; ++p) {
875:     const PetscScalar *point   = &a[p * bs];
876:     PetscInt           dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
877:     PetscBool          point_outside_domain = PETSC_FALSE;

879:     /* check bounding box of domain */
880:     for (d = 0; d < dim; d++) {
881:       if (PetscRealPart(point[d]) < gmin[d]) {
882:         point_outside_domain = PETSC_TRUE;
883:         break;
884:       }
885:       if (PetscRealPart(point[d]) > gmax[d]) {
886:         point_outside_domain = PETSC_TRUE;
887:         break;
888:       }
889:     }
890:     if (point_outside_domain) {
891:       cells[p].rank  = 0;
892:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
893:       terminating_query_type[0]++;
894:       continue;
895:     }

897:     /* check initial values in cells[].index - abort early if found */
898:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
899:       c              = cells[p].index;
900:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
901:       PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
902:       if (cell >= 0) {
903:         cells[p].rank  = 0;
904:         cells[p].index = cell;
905:         numFound++;
906:       }
907:     }
908:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
909:       terminating_query_type[1]++;
910:       continue;
911:     }

913:     if (hash) {
914:       PetscBool found_box;

916:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Checking point %" PetscInt_FMT " (%.2g, %.2g, %.2g)\n", p, (double)PetscRealPart(point[0]), (double)PetscRealPart(point[1]), (double)PetscRealPart(point[2])));
917:       /* allow for case that point is outside box - abort early */
918:       PetscCall(PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box));
919:       if (found_box) {
920:         if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", bin, dbin[0], dbin[1], dbin[2]));
921:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
922:         PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
923:         PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
924:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
925:           if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Checking for point in cell %" PetscInt_FMT "\n", boxCells[c]));
926:           PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell));
927:           if (cell >= 0) {
928:             if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "      FOUND in cell %" PetscInt_FMT "\n", cell));
929:             cells[p].rank  = 0;
930:             cells[p].index = cell;
931:             numFound++;
932:             terminating_query_type[2]++;
933:             break;
934:           }
935:         }
936:       }
937:     } else {
938:       for (c = cStart; c < cEnd; ++c) {
939:         PetscInt idx;

941:         PetscCall(PetscFindInt(c, Nl, leaves, &idx));
942:         if (idx >= 0) continue;
943:         PetscCall(DMPlexLocatePoint_Internal(dm, dim, point, c, &cell));
944:         if (cell >= 0) {
945:           cells[p].rank  = 0;
946:           cells[p].index = cell;
947:           numFound++;
948:           terminating_query_type[2]++;
949:           break;
950:         }
951:       }
952:     }
953:   }
954:   if (hash) PetscCall(ISRestoreIndices(mesh->lbox->cells, &boxCells));
955:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
956:     for (p = 0; p < numPoints; p++) {
957:       const PetscScalar *point = &a[p * bs];
958:       PetscReal          cpoint[3], diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
959:       PetscInt           dbin[3] = {-1, -1, -1}, bin, cellOffset, d, bestc = -1;

961:       if (cells[p].index < 0) {
962:         PetscCall(PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin));
963:         PetscCall(PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells));
964:         PetscCall(PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset));
965:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
966:           PetscCall(DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint));
967:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
968:           dist = DMPlex_NormD_Internal(dim, diff);
969:           if (dist < distMax) {
970:             for (d = 0; d < dim; ++d) best[d] = cpoint[d];
971:             bestc   = boxCells[c];
972:             distMax = dist;
973:           }
974:         }
975:         if (distMax < PETSC_MAX_REAL) {
976:           ++numFound;
977:           cells[p].rank  = 0;
978:           cells[p].index = bestc;
979:           for (d = 0; d < dim; ++d) a[p * bs + d] = best[d];
980:         }
981:       }
982:     }
983:   }
984:   /* This code is only be relevant when interfaced to parallel point location */
985:   /* Check for highest numbered proc that claims a point (do we care?) */
986:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
987:     PetscCall(PetscMalloc1(numFound, &found));
988:     for (p = 0, numFound = 0; p < numPoints; p++) {
989:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
990:         if (numFound < p) cells[numFound] = cells[p];
991:         found[numFound++] = p;
992:       }
993:     }
994:   }
995:   PetscCall(VecRestoreArray(v, &a));
996:   if (!reuse) PetscCall(PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER));
997:   PetscCall(PetscTime(&t1));
998:   if (hash) {
999:     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]));
1000:   } else {
1001:     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]));
1002:   }
1003:   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))));
1004:   PetscCall(PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0));
1005:   PetscFunctionReturn(PETSC_SUCCESS);
1006: }

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

1011:   Not Collective

1013:   Input/Output Parameter:
1014: . coords - The coordinates of a segment, on output the new y-coordinate, and 0 for x

1016:   Output Parameter:
1017: . R - The rotation which accomplishes the projection

1019:   Level: developer

1021: .seealso: `DMPLEX`, `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1022: @*/
1023: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
1024: {
1025:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
1026:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
1027:   const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;

1029:   PetscFunctionBegin;
1030:   R[0]      = c;
1031:   R[1]      = -s;
1032:   R[2]      = s;
1033:   R[3]      = c;
1034:   coords[0] = 0.0;
1035:   coords[1] = r;
1036:   PetscFunctionReturn(PETSC_SUCCESS);
1037: }

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

1042:   Not Collective

1044:   Input/Output Parameter:
1045: . coords - The coordinates of a segment; on output, the new y-coordinate, and 0 for x and z

1047:   Output Parameter:
1048: . R - The rotation which accomplishes the projection

1050:   Note:
1051:   This uses the basis completion described by Frisvad in http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html, DOI:10.1080/2165347X.2012.689606

1053:   Level: developer

1055: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1056: @*/
1057: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1058: {
1059:   PetscReal x    = PetscRealPart(coords[3] - coords[0]);
1060:   PetscReal y    = PetscRealPart(coords[4] - coords[1]);
1061:   PetscReal z    = PetscRealPart(coords[5] - coords[2]);
1062:   PetscReal r    = PetscSqrtReal(x * x + y * y + z * z);
1063:   PetscReal rinv = 1. / r;
1064:   PetscFunctionBegin;

1066:   x *= rinv;
1067:   y *= rinv;
1068:   z *= rinv;
1069:   if (x > 0.) {
1070:     PetscReal inv1pX = 1. / (1. + x);

1072:     R[0] = x;
1073:     R[1] = -y;
1074:     R[2] = -z;
1075:     R[3] = y;
1076:     R[4] = 1. - y * y * inv1pX;
1077:     R[5] = -y * z * inv1pX;
1078:     R[6] = z;
1079:     R[7] = -y * z * inv1pX;
1080:     R[8] = 1. - z * z * inv1pX;
1081:   } else {
1082:     PetscReal inv1mX = 1. / (1. - x);

1084:     R[0] = x;
1085:     R[1] = z;
1086:     R[2] = y;
1087:     R[3] = y;
1088:     R[4] = -y * z * inv1mX;
1089:     R[5] = 1. - y * y * inv1mX;
1090:     R[6] = z;
1091:     R[7] = 1. - z * z * inv1mX;
1092:     R[8] = -y * z * inv1mX;
1093:   }
1094:   coords[0] = 0.0;
1095:   coords[1] = r;
1096:   PetscFunctionReturn(PETSC_SUCCESS);
1097: }

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

1103:   Not Collective

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

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

1112:   Output Parameter:
1113: . 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.

1115:   Level: developer

1117: .seealso: `DMPLEX`, `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1118: @*/
1119: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1120: {
1121:   PetscReal      x1[3], x2[3], n[3], c[3], norm;
1122:   const PetscInt dim = 3;
1123:   PetscInt       d, p;

1125:   PetscFunctionBegin;
1126:   /* 0) Calculate normal vector */
1127:   for (d = 0; d < dim; ++d) {
1128:     x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1129:     x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1130:   }
1131:   // n = x1 \otimes x2
1132:   n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1133:   n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1134:   n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1135:   norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1136:   for (d = 0; d < dim; d++) n[d] /= norm;
1137:   norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1138:   for (d = 0; d < dim; d++) x1[d] /= norm;
1139:   // x2 = n \otimes x1
1140:   x2[0] = n[1] * x1[2] - n[2] * x1[1];
1141:   x2[1] = n[2] * x1[0] - n[0] * x1[2];
1142:   x2[2] = n[0] * x1[1] - n[1] * x1[0];
1143:   for (d = 0; d < dim; d++) {
1144:     R[d * dim + 0] = x1[d];
1145:     R[d * dim + 1] = x2[d];
1146:     R[d * dim + 2] = n[d];
1147:     c[d]           = PetscRealPart(coords[0 * dim + d]);
1148:   }
1149:   for (p = 0; p < coordSize / dim; p++) {
1150:     PetscReal y[3];
1151:     for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1152:     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];
1153:   }
1154:   PetscFunctionReturn(PETSC_SUCCESS);
1155: }

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

1161:    |  1  1  1 |
1162:    | x0 x1 x2 |
1163:    | y0 y1 y2 |

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

1167:    | x1 x2 |
1168:    | y1 y2 |
1169:   */
1170:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1171:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1172:   PetscReal       M[4], detM;
1173:   M[0] = x1;
1174:   M[1] = x2;
1175:   M[2] = y1;
1176:   M[3] = y2;
1177:   DMPlex_Det2D_Internal(&detM, M);
1178:   *vol = 0.5 * detM;
1179:   (void)PetscLogFlops(5.0);
1180: }

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

1186:    |  1  1  1  1 |
1187:    | x0 x1 x2 x3 |
1188:    | y0 y1 y2 y3 |
1189:    | z0 z1 z2 z3 |

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

1193:    | x1 x2 x3 |
1194:    | y1 y2 y3 |
1195:    | z1 z2 z3 |
1196:   */
1197:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1198:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1199:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1200:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1201:   PetscReal       M[9], detM;
1202:   M[0] = x1;
1203:   M[1] = x2;
1204:   M[2] = x3;
1205:   M[3] = y1;
1206:   M[4] = y2;
1207:   M[5] = y3;
1208:   M[6] = z1;
1209:   M[7] = z2;
1210:   M[8] = z3;
1211:   DMPlex_Det3D_Internal(&detM, M);
1212:   *vol = -onesixth * detM;
1213:   (void)PetscLogFlops(10.0);
1214: }

1216: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1217: {
1218:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1219:   DMPlex_Det3D_Internal(vol, coords);
1220:   *vol *= -onesixth;
1221: }

1223: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1224: {
1225:   PetscSection       coordSection;
1226:   Vec                coordinates;
1227:   const PetscScalar *coords;
1228:   PetscInt           dim, d, off;

1230:   PetscFunctionBegin;
1231:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1232:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
1233:   PetscCall(PetscSectionGetDof(coordSection, e, &dim));
1234:   if (!dim) PetscFunctionReturn(PETSC_SUCCESS);
1235:   PetscCall(PetscSectionGetOffset(coordSection, e, &off));
1236:   PetscCall(VecGetArrayRead(coordinates, &coords));
1237:   if (v0) {
1238:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1239:   }
1240:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
1241:   *detJ = 1.;
1242:   if (J) {
1243:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1244:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1245:     if (invJ) {
1246:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1247:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1248:     }
1249:   }
1250:   PetscFunctionReturn(PETSC_SUCCESS);
1251: }

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

1256:   Not Collective

1258:   Input Parameters:
1259: + dm   - The `DMPLEX`
1260: - cell - The cell number

1262:   Output Parameters:
1263: + isDG   - Using cellwise coordinates
1264: . Nc     - The number of coordinates
1265: . array  - The coordinate array
1266: - coords - The cell coordinates

1268:   Level: developer

1270: .seealso: `DMPLEX`, `DMPlexRestoreCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1271: @*/
1272: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1273: {
1274:   DM                 cdm;
1275:   Vec                coordinates;
1276:   PetscSection       cs;
1277:   const PetscScalar *ccoords;
1278:   PetscInt           pStart, pEnd;

1280:   PetscFunctionBeginHot;
1281:   *isDG   = PETSC_FALSE;
1282:   *Nc     = 0;
1283:   *array  = NULL;
1284:   *coords = NULL;
1285:   /* Check for cellwise coordinates */
1286:   PetscCall(DMGetCellCoordinateSection(dm, &cs));
1287:   if (!cs) goto cg;
1288:   /* Check that the cell exists in the cellwise section */
1289:   PetscCall(PetscSectionGetChart(cs, &pStart, &pEnd));
1290:   if (cell < pStart || cell >= pEnd) goto cg;
1291:   /* Check for cellwise coordinates for this cell */
1292:   PetscCall(PetscSectionGetDof(cs, cell, Nc));
1293:   if (!*Nc) goto cg;
1294:   /* Check for cellwise coordinates */
1295:   PetscCall(DMGetCellCoordinatesLocalNoncollective(dm, &coordinates));
1296:   if (!coordinates) goto cg;
1297:   /* Get cellwise coordinates */
1298:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1299:   PetscCall(VecGetArrayRead(coordinates, array));
1300:   PetscCall(DMPlexPointLocalRead(cdm, cell, *array, &ccoords));
1301:   PetscCall(DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1302:   PetscCall(PetscArraycpy(*coords, ccoords, *Nc));
1303:   PetscCall(VecRestoreArrayRead(coordinates, array));
1304:   *isDG = PETSC_TRUE;
1305:   PetscFunctionReturn(PETSC_SUCCESS);
1306: cg:
1307:   /* Use continuous coordinates */
1308:   PetscCall(DMGetCoordinateDM(dm, &cdm));
1309:   PetscCall(DMGetCoordinateSection(dm, &cs));
1310:   PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1311:   PetscCall(DMPlexVecGetClosure(cdm, cs, coordinates, cell, Nc, coords));
1312:   PetscFunctionReturn(PETSC_SUCCESS);
1313: }

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

1318:   Not Collective

1320:   Input Parameters:
1321: + dm   - The `DMPLEX`
1322: - cell - The cell number

1324:   Output Parameters:
1325: + isDG   - Using cellwise coordinates
1326: . Nc     - The number of coordinates
1327: . array  - The coordinate array
1328: - coords - The cell coordinates

1330:   Level: developer

1332: .seealso: `DMPLEX`, `DMPlexGetCellCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCellCoordinatesLocal()`
1333: @*/
1334: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1335: {
1336:   DM           cdm;
1337:   PetscSection cs;
1338:   Vec          coordinates;

1340:   PetscFunctionBeginHot;
1341:   if (*isDG) {
1342:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
1343:     PetscCall(DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords));
1344:   } else {
1345:     PetscCall(DMGetCoordinateDM(dm, &cdm));
1346:     PetscCall(DMGetCoordinateSection(dm, &cs));
1347:     PetscCall(DMGetCoordinatesLocalNoncollective(dm, &coordinates));
1348:     PetscCall(DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, (PetscScalar **)coords));
1349:   }
1350:   PetscFunctionReturn(PETSC_SUCCESS);
1351: }

1353: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1354: {
1355:   const PetscScalar *array;
1356:   PetscScalar       *coords = NULL;
1357:   PetscInt           numCoords, d;
1358:   PetscBool          isDG;

1360:   PetscFunctionBegin;
1361:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1362:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1363:   *detJ = 0.0;
1364:   if (numCoords == 6) {
1365:     const PetscInt dim = 3;
1366:     PetscReal      R[9], J0;

1368:     if (v0) {
1369:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1370:     }
1371:     PetscCall(DMPlexComputeProjection3Dto1D(coords, R));
1372:     if (J) {
1373:       J0   = 0.5 * PetscRealPart(coords[1]);
1374:       J[0] = R[0] * J0;
1375:       J[1] = R[1];
1376:       J[2] = R[2];
1377:       J[3] = R[3] * J0;
1378:       J[4] = R[4];
1379:       J[5] = R[5];
1380:       J[6] = R[6] * J0;
1381:       J[7] = R[7];
1382:       J[8] = R[8];
1383:       DMPlex_Det3D_Internal(detJ, J);
1384:       if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1385:     }
1386:   } else if (numCoords == 4) {
1387:     const PetscInt dim = 2;
1388:     PetscReal      R[4], J0;

1390:     if (v0) {
1391:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1392:     }
1393:     PetscCall(DMPlexComputeProjection2Dto1D(coords, R));
1394:     if (J) {
1395:       J0   = 0.5 * PetscRealPart(coords[1]);
1396:       J[0] = R[0] * J0;
1397:       J[1] = R[1];
1398:       J[2] = R[2] * J0;
1399:       J[3] = R[3];
1400:       DMPlex_Det2D_Internal(detJ, J);
1401:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1402:     }
1403:   } else if (numCoords == 2) {
1404:     const PetscInt dim = 1;

1406:     if (v0) {
1407:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1408:     }
1409:     if (J) {
1410:       J[0]  = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1411:       *detJ = J[0];
1412:       PetscCall(PetscLogFlops(2.0));
1413:       if (invJ) {
1414:         invJ[0] = 1.0 / J[0];
1415:         PetscCall(PetscLogFlops(1.0));
1416:       }
1417:     }
1418:   } 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);
1419:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1420:   PetscFunctionReturn(PETSC_SUCCESS);
1421: }

1423: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1424: {
1425:   const PetscScalar *array;
1426:   PetscScalar       *coords = NULL;
1427:   PetscInt           numCoords, d;
1428:   PetscBool          isDG;

1430:   PetscFunctionBegin;
1431:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1432:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1433:   *detJ = 0.0;
1434:   if (numCoords == 9) {
1435:     const PetscInt dim = 3;
1436:     PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1438:     if (v0) {
1439:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1440:     }
1441:     PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1442:     if (J) {
1443:       const PetscInt pdim = 2;

1445:       for (d = 0; d < pdim; d++) {
1446:         for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1447:       }
1448:       PetscCall(PetscLogFlops(8.0));
1449:       DMPlex_Det3D_Internal(detJ, J0);
1450:       for (d = 0; d < dim; d++) {
1451:         for (PetscInt f = 0; f < dim; f++) {
1452:           J[d * dim + f] = 0.0;
1453:           for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1454:         }
1455:       }
1456:       PetscCall(PetscLogFlops(18.0));
1457:     }
1458:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1459:   } else if (numCoords == 6) {
1460:     const PetscInt dim = 2;

1462:     if (v0) {
1463:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1464:     }
1465:     if (J) {
1466:       for (d = 0; d < dim; d++) {
1467:         for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1468:       }
1469:       PetscCall(PetscLogFlops(8.0));
1470:       DMPlex_Det2D_Internal(detJ, J);
1471:     }
1472:     if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1473:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1474:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1475:   PetscFunctionReturn(PETSC_SUCCESS);
1476: }

1478: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1479: {
1480:   const PetscScalar *array;
1481:   PetscScalar       *coords = NULL;
1482:   PetscInt           numCoords, d;
1483:   PetscBool          isDG;

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

1491:     if (isTensor) {
1492:       vorder[2] = 3;
1493:       vorder[3] = 2;
1494:     }
1495:     *detJ = 0.0;
1496:     if (numCoords == 12) {
1497:       const PetscInt dim = 3;
1498:       PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1500:       if (v) {
1501:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1502:       }
1503:       PetscCall(DMPlexComputeProjection3Dto2D(numCoords, coords, R));
1504:       if (J) {
1505:         const PetscInt pdim = 2;

1507:         for (d = 0; d < pdim; d++) {
1508:           J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1509:           J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1510:         }
1511:         PetscCall(PetscLogFlops(8.0));
1512:         DMPlex_Det3D_Internal(detJ, J0);
1513:         for (d = 0; d < dim; d++) {
1514:           for (PetscInt f = 0; f < dim; f++) {
1515:             J[d * dim + f] = 0.0;
1516:             for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1517:           }
1518:         }
1519:         PetscCall(PetscLogFlops(18.0));
1520:       }
1521:       if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1522:     } else if (numCoords == 8) {
1523:       const PetscInt dim = 2;

1525:       if (v) {
1526:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1527:       }
1528:       if (J) {
1529:         for (d = 0; d < dim; d++) {
1530:           J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1531:           J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1532:         }
1533:         PetscCall(PetscLogFlops(8.0));
1534:         DMPlex_Det2D_Internal(detJ, J);
1535:       }
1536:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1537:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1538:   } else {
1539:     const PetscInt Nv         = 4;
1540:     const PetscInt dimR       = 2;
1541:     PetscInt       zToPlex[4] = {0, 1, 3, 2};
1542:     PetscReal      zOrder[12];
1543:     PetscReal      zCoeff[12];
1544:     PetscInt       i, j, k, l, dim;

1546:     if (isTensor) {
1547:       zToPlex[2] = 2;
1548:       zToPlex[3] = 3;
1549:     }
1550:     if (numCoords == 12) {
1551:       dim = 3;
1552:     } else if (numCoords == 8) {
1553:       dim = 2;
1554:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1555:     for (i = 0; i < Nv; i++) {
1556:       PetscInt zi = zToPlex[i];

1558:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1559:     }
1560:     for (j = 0; j < dim; j++) {
1561:       /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
1562:            \phi^0 = (1 - xi - eta + xi eta) --> 1      = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
1563:            \phi^1 = (1 + xi - eta - xi eta) --> xi     = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
1564:            \phi^2 = (1 - xi + eta - xi eta) --> eta    = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
1565:            \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
1566:       */
1567:       zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1568:       zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1569:       zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1570:       zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1571:     }
1572:     for (i = 0; i < Nq; i++) {
1573:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1575:       if (v) {
1576:         PetscReal extPoint[4];

1578:         extPoint[0] = 1.;
1579:         extPoint[1] = xi;
1580:         extPoint[2] = eta;
1581:         extPoint[3] = xi * eta;
1582:         for (j = 0; j < dim; j++) {
1583:           PetscReal val = 0.;

1585:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1586:           v[i * dim + j] = val;
1587:         }
1588:       }
1589:       if (J) {
1590:         PetscReal extJ[8];

1592:         extJ[0] = 0.;
1593:         extJ[1] = 0.;
1594:         extJ[2] = 1.;
1595:         extJ[3] = 0.;
1596:         extJ[4] = 0.;
1597:         extJ[5] = 1.;
1598:         extJ[6] = eta;
1599:         extJ[7] = xi;
1600:         for (j = 0; j < dim; j++) {
1601:           for (k = 0; k < dimR; k++) {
1602:             PetscReal val = 0.;

1604:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1605:             J[i * dim * dim + dim * j + k] = val;
1606:           }
1607:         }
1608:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1609:           PetscReal  x, y, z;
1610:           PetscReal *iJ = &J[i * dim * dim];
1611:           PetscReal  norm;

1613:           x     = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1614:           y     = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1615:           z     = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1616:           norm  = PetscSqrtReal(x * x + y * y + z * z);
1617:           iJ[2] = x / norm;
1618:           iJ[5] = y / norm;
1619:           iJ[8] = z / norm;
1620:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1621:           if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1622:         } else {
1623:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1624:           if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1625:         }
1626:       }
1627:     }
1628:   }
1629:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1630:   PetscFunctionReturn(PETSC_SUCCESS);
1631: }

1633: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1634: {
1635:   const PetscScalar *array;
1636:   PetscScalar       *coords = NULL;
1637:   const PetscInt     dim    = 3;
1638:   PetscInt           numCoords, d;
1639:   PetscBool          isDG;

1641:   PetscFunctionBegin;
1642:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1643:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1644:   *detJ = 0.0;
1645:   if (v0) {
1646:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1647:   }
1648:   if (J) {
1649:     for (d = 0; d < dim; d++) {
1650:       /* I orient with outward face normals */
1651:       J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1652:       J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1653:       J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1654:     }
1655:     PetscCall(PetscLogFlops(18.0));
1656:     DMPlex_Det3D_Internal(detJ, J);
1657:   }
1658:   if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1659:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1660:   PetscFunctionReturn(PETSC_SUCCESS);
1661: }

1663: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1664: {
1665:   const PetscScalar *array;
1666:   PetscScalar       *coords = NULL;
1667:   const PetscInt     dim    = 3;
1668:   PetscInt           numCoords, d;
1669:   PetscBool          isDG;

1671:   PetscFunctionBegin;
1672:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1673:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1674:   if (!Nq) {
1675:     *detJ = 0.0;
1676:     if (v) {
1677:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1678:     }
1679:     if (J) {
1680:       for (d = 0; d < dim; d++) {
1681:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1682:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1683:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1684:       }
1685:       PetscCall(PetscLogFlops(18.0));
1686:       DMPlex_Det3D_Internal(detJ, J);
1687:     }
1688:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1689:   } else {
1690:     const PetscInt Nv         = 8;
1691:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1692:     const PetscInt dim        = 3;
1693:     const PetscInt dimR       = 3;
1694:     PetscReal      zOrder[24];
1695:     PetscReal      zCoeff[24];
1696:     PetscInt       i, j, k, l;

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

1701:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1702:     }
1703:     for (j = 0; j < dim; j++) {
1704:       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]);
1705:       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]);
1706:       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]);
1707:       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]);
1708:       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]);
1709:       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]);
1710:       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]);
1711:       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]);
1712:     }
1713:     for (i = 0; i < Nq; i++) {
1714:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

1716:       if (v) {
1717:         PetscReal extPoint[8];

1719:         extPoint[0] = 1.;
1720:         extPoint[1] = xi;
1721:         extPoint[2] = eta;
1722:         extPoint[3] = xi * eta;
1723:         extPoint[4] = theta;
1724:         extPoint[5] = theta * xi;
1725:         extPoint[6] = theta * eta;
1726:         extPoint[7] = theta * eta * xi;
1727:         for (j = 0; j < dim; j++) {
1728:           PetscReal val = 0.;

1730:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1731:           v[i * dim + j] = val;
1732:         }
1733:       }
1734:       if (J) {
1735:         PetscReal extJ[24];

1737:         extJ[0]  = 0.;
1738:         extJ[1]  = 0.;
1739:         extJ[2]  = 0.;
1740:         extJ[3]  = 1.;
1741:         extJ[4]  = 0.;
1742:         extJ[5]  = 0.;
1743:         extJ[6]  = 0.;
1744:         extJ[7]  = 1.;
1745:         extJ[8]  = 0.;
1746:         extJ[9]  = eta;
1747:         extJ[10] = xi;
1748:         extJ[11] = 0.;
1749:         extJ[12] = 0.;
1750:         extJ[13] = 0.;
1751:         extJ[14] = 1.;
1752:         extJ[15] = theta;
1753:         extJ[16] = 0.;
1754:         extJ[17] = xi;
1755:         extJ[18] = 0.;
1756:         extJ[19] = theta;
1757:         extJ[20] = eta;
1758:         extJ[21] = theta * eta;
1759:         extJ[22] = theta * xi;
1760:         extJ[23] = eta * xi;

1762:         for (j = 0; j < dim; j++) {
1763:           for (k = 0; k < dimR; k++) {
1764:             PetscReal val = 0.;

1766:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1767:             J[i * dim * dim + dim * j + k] = val;
1768:           }
1769:         }
1770:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1771:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1772:       }
1773:     }
1774:   }
1775:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1776:   PetscFunctionReturn(PETSC_SUCCESS);
1777: }

1779: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1780: {
1781:   const PetscScalar *array;
1782:   PetscScalar       *coords = NULL;
1783:   const PetscInt     dim    = 3;
1784:   PetscInt           numCoords, d;
1785:   PetscBool          isDG;

1787:   PetscFunctionBegin;
1788:   PetscCall(DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1789:   PetscCheck(!invJ || J, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1790:   if (!Nq) {
1791:     /* Assume that the map to the reference is affine */
1792:     *detJ = 0.0;
1793:     if (v) {
1794:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1795:     }
1796:     if (J) {
1797:       for (d = 0; d < dim; d++) {
1798:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1799:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1800:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1801:       }
1802:       PetscCall(PetscLogFlops(18.0));
1803:       DMPlex_Det3D_Internal(detJ, J);
1804:     }
1805:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1806:   } else {
1807:     const PetscInt dim  = 3;
1808:     const PetscInt dimR = 3;
1809:     const PetscInt Nv   = 6;
1810:     PetscReal      verts[18];
1811:     PetscReal      coeff[18];
1812:     PetscInt       i, j, k, l;

1814:     for (i = 0; i < Nv; ++i)
1815:       for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
1816:     for (j = 0; j < dim; ++j) {
1817:       /* Check for triangle,
1818:            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)
1819:            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)
1820:            phi^2 =  1/2 (1 + eta)   chi^2 = delta(-1,  1)

1822:            phi^0 + phi^1 + phi^2 = 1    coef_1   = 1/2 (         chi^1 + chi^2)
1823:           -phi^0 + phi^1 - phi^2 = xi   coef_xi  = 1/2 (-chi^0 + chi^1)
1824:           -phi^0 - phi^1 + phi^2 = eta  coef_eta = 1/2 (-chi^0         + chi^2)

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

1830:           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
1831:       */
1832:       /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
1833:            \phi^0 = 1/4 (   -xi - eta        + xi zeta + eta zeta) --> /  1  1  1  1  1  1 \ 1
1834:            \phi^1 = 1/4 (1      + eta - zeta           - eta zeta) --> | -1  1 -1 -1 -1  1 | eta
1835:            \phi^2 = 1/4 (1 + xi       - zeta - xi zeta)            --> | -1 -1  1 -1  1 -1 | xi
1836:            \phi^3 = 1/4 (   -xi - eta        - xi zeta - eta zeta) --> | -1 -1 -1  1  1  1 | zeta
1837:            \phi^4 = 1/4 (1 + xi       + zeta + xi zeta)            --> |  1  1 -1 -1  1 -1 | xi zeta
1838:            \phi^5 = 1/4 (1      + eta + zeta           + eta zeta) --> \  1 -1  1 -1 -1  1 / eta zeta
1839:            1/4 /  0  1  1  0  1  1 \
1840:                | -1  1  0 -1  0  1 |
1841:                | -1  0  1 -1  1  0 |
1842:                |  0 -1 -1  0  1  1 |
1843:                |  1  0 -1 -1  1  0 |
1844:                \  1 -1  0 -1  0  1 /
1845:       */
1846:       coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1847:       coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1848:       coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1849:       coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1850:       coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1851:       coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1852:       /* For reference prism:
1853:       {0, 0, 0}
1854:       {0, 1, 0}
1855:       {1, 0, 0}
1856:       {0, 0, 1}
1857:       {0, 0, 0}
1858:       {0, 0, 0}
1859:       */
1860:     }
1861:     for (i = 0; i < Nq; ++i) {
1862:       const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];

1864:       if (v) {
1865:         PetscReal extPoint[6];
1866:         PetscInt  c;

1868:         extPoint[0] = 1.;
1869:         extPoint[1] = eta;
1870:         extPoint[2] = xi;
1871:         extPoint[3] = zeta;
1872:         extPoint[4] = xi * zeta;
1873:         extPoint[5] = eta * zeta;
1874:         for (c = 0; c < dim; ++c) {
1875:           PetscReal val = 0.;

1877:           for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
1878:           v[i * dim + c] = val;
1879:         }
1880:       }
1881:       if (J) {
1882:         PetscReal extJ[18];

1884:         extJ[0]  = 0.;
1885:         extJ[1]  = 0.;
1886:         extJ[2]  = 0.;
1887:         extJ[3]  = 0.;
1888:         extJ[4]  = 1.;
1889:         extJ[5]  = 0.;
1890:         extJ[6]  = 1.;
1891:         extJ[7]  = 0.;
1892:         extJ[8]  = 0.;
1893:         extJ[9]  = 0.;
1894:         extJ[10] = 0.;
1895:         extJ[11] = 1.;
1896:         extJ[12] = zeta;
1897:         extJ[13] = 0.;
1898:         extJ[14] = xi;
1899:         extJ[15] = 0.;
1900:         extJ[16] = zeta;
1901:         extJ[17] = eta;

1903:         for (j = 0; j < dim; j++) {
1904:           for (k = 0; k < dimR; k++) {
1905:             PetscReal val = 0.;

1907:             for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
1908:             J[i * dim * dim + dim * j + k] = val;
1909:           }
1910:         }
1911:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1912:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1913:       }
1914:     }
1915:   }
1916:   PetscCall(DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords));
1917:   PetscFunctionReturn(PETSC_SUCCESS);
1918: }

1920: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1921: {
1922:   DMPolytopeType   ct;
1923:   PetscInt         depth, dim, coordDim, coneSize, i;
1924:   PetscInt         Nq     = 0;
1925:   const PetscReal *points = NULL;
1926:   DMLabel          depthLabel;
1927:   PetscReal        xi0[3]   = {-1., -1., -1.}, v0[3], J0[9], detJ0;
1928:   PetscBool        isAffine = PETSC_TRUE;

1930:   PetscFunctionBegin;
1931:   PetscCall(DMPlexGetDepth(dm, &depth));
1932:   PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
1933:   PetscCall(DMPlexGetDepthLabel(dm, &depthLabel));
1934:   PetscCall(DMLabelGetValue(depthLabel, cell, &dim));
1935:   if (depth == 1 && dim == 1) PetscCall(DMGetDimension(dm, &dim));
1936:   PetscCall(DMGetCoordinateDim(dm, &coordDim));
1937:   PetscCheck(coordDim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %" PetscInt_FMT " > 3", coordDim);
1938:   if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL));
1939:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
1940:   switch (ct) {
1941:   case DM_POLYTOPE_POINT:
1942:     PetscCall(DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ));
1943:     isAffine = PETSC_FALSE;
1944:     break;
1945:   case DM_POLYTOPE_SEGMENT:
1946:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1947:     if (Nq) PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1948:     else PetscCall(DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ));
1949:     break;
1950:   case DM_POLYTOPE_TRIANGLE:
1951:     if (Nq) PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1952:     else PetscCall(DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ));
1953:     break;
1954:   case DM_POLYTOPE_QUADRILATERAL:
1955:     PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ));
1956:     isAffine = PETSC_FALSE;
1957:     break;
1958:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1959:     PetscCall(DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ));
1960:     isAffine = PETSC_FALSE;
1961:     break;
1962:   case DM_POLYTOPE_TETRAHEDRON:
1963:     if (Nq) PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0));
1964:     else PetscCall(DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ));
1965:     break;
1966:   case DM_POLYTOPE_HEXAHEDRON:
1967:     PetscCall(DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
1968:     isAffine = PETSC_FALSE;
1969:     break;
1970:   case DM_POLYTOPE_TRI_PRISM:
1971:     PetscCall(DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ));
1972:     isAffine = PETSC_FALSE;
1973:     break;
1974:   default:
1975:     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))]);
1976:   }
1977:   if (isAffine && Nq) {
1978:     if (v) {
1979:       for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1980:     }
1981:     if (detJ) {
1982:       for (i = 0; i < Nq; i++) detJ[i] = detJ0;
1983:     }
1984:     if (J) {
1985:       PetscInt k;

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

1990:         for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
1991:       }
1992:     }
1993:     if (invJ) {
1994:       PetscInt k;
1995:       switch (coordDim) {
1996:       case 0:
1997:         break;
1998:       case 1:
1999:         invJ[0] = 1. / J0[0];
2000:         break;
2001:       case 2:
2002:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
2003:         break;
2004:       case 3:
2005:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
2006:         break;
2007:       }
2008:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
2009:         PetscInt j;

2011:         for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
2012:       }
2013:     }
2014:   }
2015:   PetscFunctionReturn(PETSC_SUCCESS);
2016: }

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

2021:   Collective

2023:   Input Parameters:
2024: + dm   - the `DMPLEX`
2025: - cell - the cell

2027:   Output Parameters:
2028: + v0   - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
2029: . J    - the Jacobian of the transform from the reference element
2030: . invJ - the inverse of the Jacobian
2031: - detJ - the Jacobian determinant

2033:   Level: advanced

2035: .seealso: `DMPLEX`, `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2036: @*/
2037: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2038: {
2039:   PetscFunctionBegin;
2040:   PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ));
2041:   PetscFunctionReturn(PETSC_SUCCESS);
2042: }

2044: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
2045: {
2046:   const PetscScalar *array;
2047:   PetscScalar       *coords = NULL;
2048:   PetscInt           numCoords;
2049:   PetscBool          isDG;
2050:   PetscQuadrature    feQuad;
2051:   const PetscReal   *quadPoints;
2052:   PetscTabulation    T;
2053:   PetscInt           dim, cdim, pdim, qdim, Nq, q;

2055:   PetscFunctionBegin;
2056:   PetscCall(DMGetDimension(dm, &dim));
2057:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2058:   PetscCall(DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2059:   if (!quad) { /* use the first point of the first functional of the dual space */
2060:     PetscDualSpace dsp;

2062:     PetscCall(PetscFEGetDualSpace(fe, &dsp));
2063:     PetscCall(PetscDualSpaceGetFunctional(dsp, 0, &quad));
2064:     PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2065:     Nq = 1;
2066:   } else {
2067:     PetscCall(PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL));
2068:   }
2069:   PetscCall(PetscFEGetDimension(fe, &pdim));
2070:   PetscCall(PetscFEGetQuadrature(fe, &feQuad));
2071:   if (feQuad == quad) {
2072:     PetscCall(PetscFEGetCellTabulation(fe, J ? 1 : 0, &T));
2073:     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);
2074:   } else {
2075:     PetscCall(PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T));
2076:   }
2077:   PetscCheck(qdim == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %" PetscInt_FMT " != quadrature dimension %" PetscInt_FMT, dim, qdim);
2078:   {
2079:     const PetscReal *basis    = T->T[0];
2080:     const PetscReal *basisDer = T->T[1];
2081:     PetscReal        detJt;

2083: #if defined(PETSC_USE_DEBUG)
2084:     PetscCheck(Nq == T->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Np %" PetscInt_FMT " != %" PetscInt_FMT, Nq, T->Np);
2085:     PetscCheck(pdim == T->Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nb %" PetscInt_FMT " != %" PetscInt_FMT, pdim, T->Nb);
2086:     PetscCheck(dim == T->Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nc %" PetscInt_FMT " != %" PetscInt_FMT, dim, T->Nc);
2087:     PetscCheck(cdim == T->cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "cdim %" PetscInt_FMT " != %" PetscInt_FMT, cdim, T->cdim);
2088: #endif
2089:     if (v) {
2090:       PetscCall(PetscArrayzero(v, Nq * cdim));
2091:       for (q = 0; q < Nq; ++q) {
2092:         PetscInt i, k;

2094:         for (k = 0; k < pdim; ++k) {
2095:           const PetscInt vertex = k / cdim;
2096:           for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2097:         }
2098:         PetscCall(PetscLogFlops(2.0 * pdim * cdim));
2099:       }
2100:     }
2101:     if (J) {
2102:       PetscCall(PetscArrayzero(J, Nq * cdim * cdim));
2103:       for (q = 0; q < Nq; ++q) {
2104:         PetscInt i, j, k, c, r;

2106:         /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2107:         for (k = 0; k < pdim; ++k) {
2108:           const PetscInt vertex = k / cdim;
2109:           for (j = 0; j < dim; ++j) {
2110:             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]);
2111:           }
2112:         }
2113:         PetscCall(PetscLogFlops(2.0 * pdim * dim * cdim));
2114:         if (cdim > dim) {
2115:           for (c = dim; c < cdim; ++c)
2116:             for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2117:         }
2118:         if (!detJ && !invJ) continue;
2119:         detJt = 0.;
2120:         switch (cdim) {
2121:         case 3:
2122:           DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2123:           if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2124:           break;
2125:         case 2:
2126:           DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2127:           if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2128:           break;
2129:         case 1:
2130:           detJt = J[q * cdim * dim];
2131:           if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2132:         }
2133:         if (detJ) detJ[q] = detJt;
2134:       }
2135:     } else PetscCheck(!detJ && !invJ, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
2136:   }
2137:   if (feQuad != quad) PetscCall(PetscTabulationDestroy(&T));
2138:   PetscCall(DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords));
2139:   PetscFunctionReturn(PETSC_SUCCESS);
2140: }

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

2145:   Collective

2147:   Input Parameters:
2148: + dm   - the `DMPLEX`
2149: . cell - the cell
2150: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If `quad` is `NULL`, geometry will be
2151:          evaluated at the first vertex of the reference element

2153:   Output Parameters:
2154: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
2155: . J    - the Jacobian of the transform from the reference element at each quadrature point
2156: . invJ - the inverse of the Jacobian at each quadrature point
2157: - detJ - the Jacobian determinant at each quadrature point

2159:   Level: advanced

2161: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2162: @*/
2163: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2164: {
2165:   DM      cdm;
2166:   PetscFE fe = NULL;

2168:   PetscFunctionBegin;
2170:   PetscCall(DMGetCoordinateDM(dm, &cdm));
2171:   if (cdm) {
2172:     PetscClassId id;
2173:     PetscInt     numFields;
2174:     PetscDS      prob;
2175:     PetscObject  disc;

2177:     PetscCall(DMGetNumFields(cdm, &numFields));
2178:     if (numFields) {
2179:       PetscCall(DMGetDS(cdm, &prob));
2180:       PetscCall(PetscDSGetDiscretization(prob, 0, &disc));
2181:       PetscCall(PetscObjectGetClassId(disc, &id));
2182:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2183:     }
2184:   }
2185:   if (!fe) PetscCall(DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ));
2186:   else PetscCall(DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ));
2187:   PetscFunctionReturn(PETSC_SUCCESS);
2188: }

2190: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2191: {
2192:   PetscSection       coordSection;
2193:   Vec                coordinates;
2194:   const PetscScalar *coords = NULL;
2195:   PetscInt           d, dof, off;

2197:   PetscFunctionBegin;
2198:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2199:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2200:   PetscCall(VecGetArrayRead(coordinates, &coords));

2202:   /* for a point the centroid is just the coord */
2203:   if (centroid) {
2204:     PetscCall(PetscSectionGetDof(coordSection, cell, &dof));
2205:     PetscCall(PetscSectionGetOffset(coordSection, cell, &off));
2206:     for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2207:   }
2208:   if (normal) {
2209:     const PetscInt *support, *cones;
2210:     PetscInt        supportSize;
2211:     PetscReal       norm, sign;

2213:     /* compute the norm based upon the support centroids */
2214:     PetscCall(DMPlexGetSupportSize(dm, cell, &supportSize));
2215:     PetscCall(DMPlexGetSupport(dm, cell, &support));
2216:     PetscCall(DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL));

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

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

2227:     norm = DMPlex_NormD_Internal(dim, normal);
2228:     for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2229:   }
2230:   if (vol) *vol = 1.0;
2231:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
2232:   PetscFunctionReturn(PETSC_SUCCESS);
2233: }

2235: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2236: {
2237:   const PetscScalar *array;
2238:   PetscScalar       *coords = NULL;
2239:   PetscInt           cdim, coordSize, d;
2240:   PetscBool          isDG;

2242:   PetscFunctionBegin;
2243:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2244:   PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2245:   PetscCheck(coordSize == cdim * 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge has %" PetscInt_FMT " coordinates != %" PetscInt_FMT, coordSize, cdim * 2);
2246:   if (centroid) {
2247:     for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2248:   }
2249:   if (normal) {
2250:     PetscReal norm;

2252:     switch (cdim) {
2253:     case 3:
2254:       normal[2] = 0.; /* fall through */
2255:     case 2:
2256:       normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2257:       normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2258:       break;
2259:     case 1:
2260:       normal[0] = 1.0;
2261:       break;
2262:     default:
2263:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2264:     }
2265:     norm = DMPlex_NormD_Internal(cdim, normal);
2266:     for (d = 0; d < cdim; ++d) normal[d] /= norm;
2267:   }
2268:   if (vol) {
2269:     *vol = 0.0;
2270:     for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2271:     *vol = PetscSqrtReal(*vol);
2272:   }
2273:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2274:   PetscFunctionReturn(PETSC_SUCCESS);
2275: }

2277: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2278: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2279: {
2280:   DMPolytopeType     ct;
2281:   const PetscScalar *array;
2282:   PetscScalar       *coords = NULL;
2283:   PetscInt           coordSize;
2284:   PetscBool          isDG;
2285:   PetscInt           fv[4] = {0, 1, 2, 3};
2286:   PetscInt           cdim, numCorners, p, d;

2288:   PetscFunctionBegin;
2289:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2290:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
2291:   switch (ct) {
2292:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2293:     fv[2] = 3;
2294:     fv[3] = 2;
2295:     break;
2296:   default:
2297:     break;
2298:   }
2299:   PetscCall(DMGetCoordinateDim(dm, &cdim));
2300:   PetscCall(DMPlexGetConeSize(dm, cell, &numCorners));
2301:   PetscCall(DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2302:   {
2303:     PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;

2305:     for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2306:     for (p = 0; p < numCorners - 2; ++p) {
2307:       PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2308:       for (d = 0; d < cdim; d++) {
2309:         e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2310:         e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2311:       }
2312:       const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2313:       const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2314:       const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2315:       const PetscReal a  = PetscSqrtReal(dx * dx + dy * dy + dz * dz);

2317:       n[0] += dx;
2318:       n[1] += dy;
2319:       n[2] += dz;
2320:       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.;
2321:     }
2322:     norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2323:     // Allow zero volume cells
2324:     if (norm != 0) {
2325:       n[0] /= norm;
2326:       n[1] /= norm;
2327:       n[2] /= norm;
2328:       c[0] /= norm;
2329:       c[1] /= norm;
2330:       c[2] /= norm;
2331:     }
2332:     if (vol) *vol = 0.5 * norm;
2333:     if (centroid)
2334:       for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2335:     if (normal)
2336:       for (d = 0; d < cdim; ++d) normal[d] = n[d];
2337:   }
2338:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2339:   PetscFunctionReturn(PETSC_SUCCESS);
2340: }

2342: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2343: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2344: {
2345:   DMPolytopeType        ct;
2346:   const PetscScalar    *array;
2347:   PetscScalar          *coords = NULL;
2348:   PetscInt              coordSize;
2349:   PetscBool             isDG;
2350:   PetscReal             vsum      = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2351:   const PetscInt        order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2352:   const PetscInt       *cone, *faceSizes, *faces;
2353:   const DMPolytopeType *faceTypes;
2354:   PetscBool             isHybrid = PETSC_FALSE;
2355:   PetscInt              numFaces, f, fOff = 0, p, d;

2357:   PetscFunctionBegin;
2358:   PetscCheck(dim <= 3, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No support for dim %" PetscInt_FMT " > 3", dim);
2359:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2360:   PetscCall(DMPlexGetCellType(dm, cell, &ct));
2361:   switch (ct) {
2362:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2363:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2364:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
2365:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2366:     isHybrid = PETSC_TRUE;
2367:   default:
2368:     break;
2369:   }

2371:   if (centroid)
2372:     for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2373:   PetscCall(DMPlexGetCone(dm, cell, &cone));

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

2381:     // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2382:     // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2383:     // so that all tetrahedra have positive volume.
2384:     if (f == 0)
2385:       for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2386:     switch (faceTypes[f]) {
2387:     case DM_POLYTOPE_TRIANGLE:
2388:       for (d = 0; d < dim; ++d) {
2389:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2390:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2391:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2392:       }
2393:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2394:       if (flip) vtmp = -vtmp;
2395:       vsum += vtmp;
2396:       if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2397:         for (d = 0; d < dim; ++d) {
2398:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2399:         }
2400:       }
2401:       break;
2402:     case DM_POLYTOPE_QUADRILATERAL:
2403:     case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2404:       PetscInt fv[4] = {0, 1, 2, 3};

2406:       /* Side faces for hybrid cells are are stored as tensor products */
2407:       if (isHybrid && f > 1) {
2408:         fv[2] = 3;
2409:         fv[3] = 2;
2410:       }
2411:       /* DO FOR PYRAMID */
2412:       /* First tet */
2413:       for (d = 0; d < dim; ++d) {
2414:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2415:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2416:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2417:       }
2418:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2419:       if (flip) vtmp = -vtmp;
2420:       vsum += vtmp;
2421:       if (centroid) {
2422:         for (d = 0; d < dim; ++d) {
2423:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2424:         }
2425:       }
2426:       /* Second tet */
2427:       for (d = 0; d < dim; ++d) {
2428:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2429:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2430:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2431:       }
2432:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2433:       if (flip) vtmp = -vtmp;
2434:       vsum += vtmp;
2435:       if (centroid) {
2436:         for (d = 0; d < dim; ++d) {
2437:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2438:         }
2439:       }
2440:       break;
2441:     }
2442:     default:
2443:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2444:     }
2445:     fOff += faceSizes[f];
2446:   }
2447:   PetscCall(DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces));
2448:   PetscCall(DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords));
2449:   if (vol) *vol = PetscAbsReal(vsum);
2450:   if (normal)
2451:     for (d = 0; d < dim; ++d) normal[d] = 0.0;
2452:   if (centroid)
2453:     for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2454:   PetscFunctionReturn(PETSC_SUCCESS);
2455: }

2457: /*@C
2458:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

2460:   Collective

2462:   Input Parameters:
2463: + dm   - the `DMPLEX`
2464: - cell - the cell

2466:   Output Parameters:
2467: + volume   - the cell volume
2468: . centroid - the cell centroid
2469: - normal - the cell normal, if appropriate

2471:   Level: advanced

2473: .seealso: `DMPLEX`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
2474: @*/
2475: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2476: {
2477:   PetscInt depth, dim;

2479:   PetscFunctionBegin;
2480:   PetscCall(DMPlexGetDepth(dm, &depth));
2481:   PetscCall(DMGetDimension(dm, &dim));
2482:   PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2483:   PetscCall(DMPlexGetPointDepth(dm, cell, &depth));
2484:   switch (depth) {
2485:   case 0:
2486:     PetscCall(DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal));
2487:     break;
2488:   case 1:
2489:     PetscCall(DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal));
2490:     break;
2491:   case 2:
2492:     PetscCall(DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal));
2493:     break;
2494:   case 3:
2495:     PetscCall(DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal));
2496:     break;
2497:   default:
2498:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2499:   }
2500:   PetscFunctionReturn(PETSC_SUCCESS);
2501: }

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

2506:   Input Parameter:
2507: . dm - The `DMPLEX`

2509:   Output Parameters:
2510: + cellgeom - A `Vec` of `PetscFVCellGeom` data
2511: - facegeom - A `Vec` of `PetscFVFaceGeom` data

2513:   Level: developer

2515: .seealso: `DMPLEX`, `PetscFVFaceGeom`, `PetscFVCellGeom`
2516: @*/
2517: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2518: {
2519:   DM           dmFace, dmCell;
2520:   DMLabel      ghostLabel;
2521:   PetscSection sectionFace, sectionCell;
2522:   PetscSection coordSection;
2523:   Vec          coordinates;
2524:   PetscScalar *fgeom, *cgeom;
2525:   PetscReal    minradius, gminradius;
2526:   PetscInt     dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

2528:   PetscFunctionBegin;
2529:   PetscCall(DMGetDimension(dm, &dim));
2530:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
2531:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
2532:   /* Make cell centroids and volumes */
2533:   PetscCall(DMClone(dm, &dmCell));
2534:   PetscCall(DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection));
2535:   PetscCall(DMSetCoordinatesLocal(dmCell, coordinates));
2536:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
2537:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2538:   PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2539:   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
2540:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVCellGeom)) / sizeof(PetscScalar))));
2541:   PetscCall(PetscSectionSetUp(sectionCell));
2542:   PetscCall(DMSetLocalSection(dmCell, sectionCell));
2543:   PetscCall(PetscSectionDestroy(&sectionCell));
2544:   PetscCall(DMCreateLocalVector(dmCell, cellgeom));
2545:   if (cEndInterior < 0) cEndInterior = cEnd;
2546:   PetscCall(VecGetArray(*cellgeom, &cgeom));
2547:   for (c = cStart; c < cEndInterior; ++c) {
2548:     PetscFVCellGeom *cg;

2550:     PetscCall(DMPlexPointLocalRef(dmCell, c, cgeom, &cg));
2551:     PetscCall(PetscArrayzero(cg, 1));
2552:     PetscCall(DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL));
2553:   }
2554:   /* Compute face normals and minimum cell radius */
2555:   PetscCall(DMClone(dm, &dmFace));
2556:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionFace));
2557:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2558:   PetscCall(PetscSectionSetChart(sectionFace, fStart, fEnd));
2559:   for (f = fStart; f < fEnd; ++f) PetscCall(PetscSectionSetDof(sectionFace, f, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFVFaceGeom)) / sizeof(PetscScalar))));
2560:   PetscCall(PetscSectionSetUp(sectionFace));
2561:   PetscCall(DMSetLocalSection(dmFace, sectionFace));
2562:   PetscCall(PetscSectionDestroy(&sectionFace));
2563:   PetscCall(DMCreateLocalVector(dmFace, facegeom));
2564:   PetscCall(VecGetArray(*facegeom, &fgeom));
2565:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2566:   minradius = PETSC_MAX_REAL;
2567:   for (f = fStart; f < fEnd; ++f) {
2568:     PetscFVFaceGeom *fg;
2569:     PetscReal        area;
2570:     const PetscInt  *cells;
2571:     PetscInt         ncells, ghost = -1, d, numChildren;

2573:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2574:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2575:     PetscCall(DMPlexGetSupport(dm, f, &cells));
2576:     PetscCall(DMPlexGetSupportSize(dm, f, &ncells));
2577:     /* It is possible to get a face with no support when using partition overlap */
2578:     if (!ncells || ghost >= 0 || numChildren) continue;
2579:     PetscCall(DMPlexPointLocalRef(dmFace, f, fgeom, &fg));
2580:     PetscCall(DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal));
2581:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2582:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2583:     {
2584:       PetscFVCellGeom *cL, *cR;
2585:       PetscReal       *lcentroid, *rcentroid;
2586:       PetscReal        l[3], r[3], v[3];

2588:       PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL));
2589:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2590:       if (ncells > 1) {
2591:         PetscCall(DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR));
2592:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2593:       } else {
2594:         rcentroid = fg->centroid;
2595:       }
2596:       PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l));
2597:       PetscCall(DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r));
2598:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2599:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2600:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2601:       }
2602:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2603:         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]);
2604:         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]);
2605:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Direction for face %" PetscInt_FMT " could not be fixed", f);
2606:       }
2607:       if (cells[0] < cEndInterior) {
2608:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2609:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2610:       }
2611:       if (ncells > 1 && cells[1] < cEndInterior) {
2612:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2613:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2614:       }
2615:     }
2616:   }
2617:   PetscCall(MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
2618:   PetscCall(DMPlexSetMinRadius(dm, gminradius));
2619:   /* Compute centroids of ghost cells */
2620:   for (c = cEndInterior; c < cEnd; ++c) {
2621:     PetscFVFaceGeom *fg;
2622:     const PetscInt  *cone, *support;
2623:     PetscInt         coneSize, supportSize, s;

2625:     PetscCall(DMPlexGetConeSize(dmCell, c, &coneSize));
2626:     PetscCheck(coneSize == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %" PetscInt_FMT " has cone size %" PetscInt_FMT " != 1", c, coneSize);
2627:     PetscCall(DMPlexGetCone(dmCell, c, &cone));
2628:     PetscCall(DMPlexGetSupportSize(dmCell, cone[0], &supportSize));
2629:     PetscCheck(supportSize == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has support size %" PetscInt_FMT " != 2", cone[0], supportSize);
2630:     PetscCall(DMPlexGetSupport(dmCell, cone[0], &support));
2631:     PetscCall(DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg));
2632:     for (s = 0; s < 2; ++s) {
2633:       /* Reflect ghost centroid across plane of face */
2634:       if (support[s] == c) {
2635:         PetscFVCellGeom *ci;
2636:         PetscFVCellGeom *cg;
2637:         PetscReal        c2f[3], a;

2639:         PetscCall(DMPlexPointLocalRead(dmCell, support[(s + 1) % 2], cgeom, &ci));
2640:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2641:         a = DMPlex_DotRealD_Internal(dim, c2f, fg->normal) / DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2642:         PetscCall(DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg));
2643:         DMPlex_WaxpyD_Internal(dim, 2 * a, fg->normal, ci->centroid, cg->centroid);
2644:         cg->volume = ci->volume;
2645:       }
2646:     }
2647:   }
2648:   PetscCall(VecRestoreArray(*facegeom, &fgeom));
2649:   PetscCall(VecRestoreArray(*cellgeom, &cgeom));
2650:   PetscCall(DMDestroy(&dmCell));
2651:   PetscCall(DMDestroy(&dmFace));
2652:   PetscFunctionReturn(PETSC_SUCCESS);
2653: }

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

2658:   Not Collective

2660:   Input Parameter:
2661: . dm - the `DMPLEX`

2663:   Output Parameter:
2664: . minradius - the minimum cell radius

2666:   Level: developer

2668: .seealso: `DMPLEX`, `DMGetCoordinates()`
2669: @*/
2670: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2671: {
2672:   PetscFunctionBegin;
2675:   *minradius = ((DM_Plex *)dm->data)->minradius;
2676:   PetscFunctionReturn(PETSC_SUCCESS);
2677: }

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

2682:   Logically Collective

2684:   Input Parameters:
2685: + dm - the `DMPLEX`
2686: - minradius - the minimum cell radius

2688:   Level: developer

2690: .seealso: `DMPLEX`, `DMSetCoordinates()`
2691: @*/
2692: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2693: {
2694:   PetscFunctionBegin;
2696:   ((DM_Plex *)dm->data)->minradius = minradius;
2697:   PetscFunctionReturn(PETSC_SUCCESS);
2698: }

2700: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2701: {
2702:   DMLabel      ghostLabel;
2703:   PetscScalar *dx, *grad, **gref;
2704:   PetscInt     dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

2706:   PetscFunctionBegin;
2707:   PetscCall(DMGetDimension(dm, &dim));
2708:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2709:   PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2710:   cEndInterior = cEndInterior < 0 ? cEnd : cEndInterior;
2711:   PetscCall(DMPlexGetMaxSizes(dm, &maxNumFaces, NULL));
2712:   PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
2713:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2714:   PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
2715:   for (c = cStart; c < cEndInterior; c++) {
2716:     const PetscInt  *faces;
2717:     PetscInt         numFaces, usedFaces, f, d;
2718:     PetscFVCellGeom *cg;
2719:     PetscBool        boundary;
2720:     PetscInt         ghost;

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

2726:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
2727:     PetscCall(DMPlexGetConeSize(dm, c, &numFaces));
2728:     PetscCall(DMPlexGetCone(dm, c, &faces));
2729:     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);
2730:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2731:       PetscFVCellGeom *cg1;
2732:       PetscFVFaceGeom *fg;
2733:       const PetscInt  *fcells;
2734:       PetscInt         ncell, side;

2736:       PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
2737:       PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
2738:       if ((ghost >= 0) || boundary) continue;
2739:       PetscCall(DMPlexGetSupport(dm, faces[f], &fcells));
2740:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2741:       ncell = fcells[!side];    /* the neighbor */
2742:       PetscCall(DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg));
2743:       PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
2744:       for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
2745:       gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2746:     }
2747:     PetscCheck(usedFaces, PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2748:     PetscCall(PetscFVComputeGradient(fvm, usedFaces, dx, grad));
2749:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2750:       PetscCall(DMLabelGetValue(ghostLabel, faces[f], &ghost));
2751:       PetscCall(DMIsBoundaryPoint(dm, faces[f], &boundary));
2752:       if ((ghost >= 0) || boundary) continue;
2753:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
2754:       ++usedFaces;
2755:     }
2756:   }
2757:   PetscCall(PetscFree3(dx, grad, gref));
2758:   PetscFunctionReturn(PETSC_SUCCESS);
2759: }

2761: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2762: {
2763:   DMLabel      ghostLabel;
2764:   PetscScalar *dx, *grad, **gref;
2765:   PetscInt     dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2766:   PetscSection neighSec;
2767:   PetscInt(*neighbors)[2];
2768:   PetscInt *counter;

2770:   PetscFunctionBegin;
2771:   PetscCall(DMGetDimension(dm, &dim));
2772:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2773:   PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2774:   if (cEndInterior < 0) cEndInterior = cEnd;
2775:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &neighSec));
2776:   PetscCall(PetscSectionSetChart(neighSec, cStart, cEndInterior));
2777:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
2778:   PetscCall(DMGetLabel(dm, "ghost", &ghostLabel));
2779:   for (f = fStart; f < fEnd; f++) {
2780:     const PetscInt *fcells;
2781:     PetscBool       boundary;
2782:     PetscInt        ghost = -1;
2783:     PetscInt        numChildren, numCells, c;

2785:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2786:     PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
2787:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2788:     if ((ghost >= 0) || boundary || numChildren) continue;
2789:     PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
2790:     if (numCells == 2) {
2791:       PetscCall(DMPlexGetSupport(dm, f, &fcells));
2792:       for (c = 0; c < 2; c++) {
2793:         PetscInt cell = fcells[c];

2795:         if (cell >= cStart && cell < cEndInterior) PetscCall(PetscSectionAddDof(neighSec, cell, 1));
2796:       }
2797:     }
2798:   }
2799:   PetscCall(PetscSectionSetUp(neighSec));
2800:   PetscCall(PetscSectionGetMaxDof(neighSec, &maxNumFaces));
2801:   PetscCall(PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces));
2802:   nStart = 0;
2803:   PetscCall(PetscSectionGetStorageSize(neighSec, &nEnd));
2804:   PetscCall(PetscMalloc1((nEnd - nStart), &neighbors));
2805:   PetscCall(PetscCalloc1((cEndInterior - cStart), &counter));
2806:   for (f = fStart; f < fEnd; f++) {
2807:     const PetscInt *fcells;
2808:     PetscBool       boundary;
2809:     PetscInt        ghost = -1;
2810:     PetscInt        numChildren, numCells, c;

2812:     if (ghostLabel) PetscCall(DMLabelGetValue(ghostLabel, f, &ghost));
2813:     PetscCall(DMIsBoundaryPoint(dm, f, &boundary));
2814:     PetscCall(DMPlexGetTreeChildren(dm, f, &numChildren, NULL));
2815:     if ((ghost >= 0) || boundary || numChildren) continue;
2816:     PetscCall(DMPlexGetSupportSize(dm, f, &numCells));
2817:     if (numCells == 2) {
2818:       PetscCall(DMPlexGetSupport(dm, f, &fcells));
2819:       for (c = 0; c < 2; c++) {
2820:         PetscInt cell = fcells[c], off;

2822:         if (cell >= cStart && cell < cEndInterior) {
2823:           PetscCall(PetscSectionGetOffset(neighSec, cell, &off));
2824:           off += counter[cell - cStart]++;
2825:           neighbors[off][0] = f;
2826:           neighbors[off][1] = fcells[1 - c];
2827:         }
2828:       }
2829:     }
2830:   }
2831:   PetscCall(PetscFree(counter));
2832:   PetscCall(PetscMalloc3(maxNumFaces * dim, &dx, maxNumFaces * dim, &grad, maxNumFaces, &gref));
2833:   for (c = cStart; c < cEndInterior; c++) {
2834:     PetscInt         numFaces, f, d, off, ghost = -1;
2835:     PetscFVCellGeom *cg;

2837:     PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
2838:     PetscCall(PetscSectionGetDof(neighSec, c, &numFaces));
2839:     PetscCall(PetscSectionGetOffset(neighSec, c, &off));

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

2845:     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);
2846:     for (f = 0; f < numFaces; ++f) {
2847:       PetscFVCellGeom *cg1;
2848:       PetscFVFaceGeom *fg;
2849:       const PetscInt  *fcells;
2850:       PetscInt         ncell, side, nface;

2852:       nface = neighbors[off + f][0];
2853:       ncell = neighbors[off + f][1];
2854:       PetscCall(DMPlexGetSupport(dm, nface, &fcells));
2855:       side = (c != fcells[0]);
2856:       PetscCall(DMPlexPointLocalRef(dmFace, nface, fgeom, &fg));
2857:       PetscCall(DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1));
2858:       for (d = 0; d < dim; ++d) dx[f * dim + d] = cg1->centroid[d] - cg->centroid[d];
2859:       gref[f] = fg->grad[side]; /* Gradient reconstruction term will go here */
2860:     }
2861:     PetscCall(PetscFVComputeGradient(fvm, numFaces, dx, grad));
2862:     for (f = 0; f < numFaces; ++f) {
2863:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f * dim + d];
2864:     }
2865:   }
2866:   PetscCall(PetscFree3(dx, grad, gref));
2867:   PetscCall(PetscSectionDestroy(&neighSec));
2868:   PetscCall(PetscFree(neighbors));
2869:   PetscFunctionReturn(PETSC_SUCCESS);
2870: }

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

2875:   Collective

2877:   Input Parameters:
2878: + dm  - The `DMPLEX`
2879: . fvm - The `PetscFV`
2880: - cellGeometry - The face geometry from `DMPlexComputeCellGeometryFVM()`

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

2886:   Output Parameter:
2887: . dmGrad - The `DM` describing the layout of gradient data

2889:   Level: developer

2891: .seealso: `DMPLEX`, `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
2892: @*/
2893: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2894: {
2895:   DM           dmFace, dmCell;
2896:   PetscScalar *fgeom, *cgeom;
2897:   PetscSection sectionGrad, parentSection;
2898:   PetscInt     dim, pdim, cStart, cEnd, cEndInterior, c;

2900:   PetscFunctionBegin;
2901:   PetscCall(DMGetDimension(dm, &dim));
2902:   PetscCall(PetscFVGetNumComponents(fvm, &pdim));
2903:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
2904:   PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
2905:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2906:   PetscCall(VecGetDM(faceGeometry, &dmFace));
2907:   PetscCall(VecGetDM(cellGeometry, &dmCell));
2908:   PetscCall(VecGetArray(faceGeometry, &fgeom));
2909:   PetscCall(VecGetArray(cellGeometry, &cgeom));
2910:   PetscCall(DMPlexGetTree(dm, &parentSection, NULL, NULL, NULL, NULL));
2911:   if (!parentSection) {
2912:     PetscCall(BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom));
2913:   } else {
2914:     PetscCall(BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom));
2915:   }
2916:   PetscCall(VecRestoreArray(faceGeometry, &fgeom));
2917:   PetscCall(VecRestoreArray(cellGeometry, &cgeom));
2918:   /* Create storage for gradients */
2919:   PetscCall(DMClone(dm, dmGrad));
2920:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionGrad));
2921:   PetscCall(PetscSectionSetChart(sectionGrad, cStart, cEnd));
2922:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionGrad, c, pdim * dim));
2923:   PetscCall(PetscSectionSetUp(sectionGrad));
2924:   PetscCall(DMSetLocalSection(*dmGrad, sectionGrad));
2925:   PetscCall(PetscSectionDestroy(&sectionGrad));
2926:   PetscFunctionReturn(PETSC_SUCCESS);
2927: }

2929: /*@
2930:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2932:   Collective

2934:   Input Parameters:
2935: + dm  - The `DM`
2936: - fv  - The `PetscFV`

2938:   Output Parameters:
2939: + cellGeometry - The cell geometry
2940: . faceGeometry - The face geometry
2941: - gradDM       - The gradient matrices

2943:   Level: developer

2945: .seealso: `DMPLEX`, `DMPlexComputeGeometryFVM()`
2946: @*/
2947: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2948: {
2949:   PetscObject cellgeomobj, facegeomobj;

2951:   PetscFunctionBegin;
2952:   PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
2953:   if (!cellgeomobj) {
2954:     Vec cellgeomInt, facegeomInt;

2956:     PetscCall(DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt));
2957:     PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt));
2958:     PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt));
2959:     PetscCall(VecDestroy(&cellgeomInt));
2960:     PetscCall(VecDestroy(&facegeomInt));
2961:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj));
2962:   }
2963:   PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj));
2964:   if (cellgeom) *cellgeom = (Vec)cellgeomobj;
2965:   if (facegeom) *facegeom = (Vec)facegeomobj;
2966:   if (gradDM) {
2967:     PetscObject gradobj;
2968:     PetscBool   computeGradients;

2970:     PetscCall(PetscFVGetComputeGradients(fv, &computeGradients));
2971:     if (!computeGradients) {
2972:       *gradDM = NULL;
2973:       PetscFunctionReturn(PETSC_SUCCESS);
2974:     }
2975:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
2976:     if (!gradobj) {
2977:       DM dmGradInt;

2979:       PetscCall(DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt));
2980:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt));
2981:       PetscCall(DMDestroy(&dmGradInt));
2982:       PetscCall(PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj));
2983:     }
2984:     *gradDM = (DM)gradobj;
2985:   }
2986:   PetscFunctionReturn(PETSC_SUCCESS);
2987: }

2989: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2990: {
2991:   PetscInt l, m;

2993:   PetscFunctionBeginHot;
2994:   if (dimC == dimR && dimR <= 3) {
2995:     /* invert Jacobian, multiply */
2996:     PetscScalar det, idet;

2998:     switch (dimR) {
2999:     case 1:
3000:       invJ[0] = 1. / J[0];
3001:       break;
3002:     case 2:
3003:       det     = J[0] * J[3] - J[1] * J[2];
3004:       idet    = 1. / det;
3005:       invJ[0] = J[3] * idet;
3006:       invJ[1] = -J[1] * idet;
3007:       invJ[2] = -J[2] * idet;
3008:       invJ[3] = J[0] * idet;
3009:       break;
3010:     case 3: {
3011:       invJ[0] = J[4] * J[8] - J[5] * J[7];
3012:       invJ[1] = J[2] * J[7] - J[1] * J[8];
3013:       invJ[2] = J[1] * J[5] - J[2] * J[4];
3014:       det     = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
3015:       idet    = 1. / det;
3016:       invJ[0] *= idet;
3017:       invJ[1] *= idet;
3018:       invJ[2] *= idet;
3019:       invJ[3] = idet * (J[5] * J[6] - J[3] * J[8]);
3020:       invJ[4] = idet * (J[0] * J[8] - J[2] * J[6]);
3021:       invJ[5] = idet * (J[2] * J[3] - J[0] * J[5]);
3022:       invJ[6] = idet * (J[3] * J[7] - J[4] * J[6]);
3023:       invJ[7] = idet * (J[1] * J[6] - J[0] * J[7]);
3024:       invJ[8] = idet * (J[0] * J[4] - J[1] * J[3]);
3025:     } break;
3026:     }
3027:     for (l = 0; l < dimR; l++) {
3028:       for (m = 0; m < dimC; m++) guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
3029:     }
3030:   } else {
3031: #if defined(PETSC_USE_COMPLEX)
3032:     char transpose = 'C';
3033: #else
3034:     char transpose = 'T';
3035: #endif
3036:     PetscBLASInt m        = dimR;
3037:     PetscBLASInt n        = dimC;
3038:     PetscBLASInt one      = 1;
3039:     PetscBLASInt worksize = dimR * dimC, info;

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

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

3046:     for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3047:   }
3048:   PetscFunctionReturn(PETSC_SUCCESS);
3049: }

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

3058:   PetscFunctionBegin;
3060:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3061:   PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3062:   PetscCall(DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3063:   PetscCall(DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3064:   cellCoords = &cellData[0];
3065:   cellCoeffs = &cellData[coordSize];
3066:   extJ       = &cellData[2 * coordSize];
3067:   resNeg     = &cellData[2 * coordSize + dimR];
3068:   invJ       = &J[dimR * dimC];
3069:   work       = &J[2 * dimR * dimC];
3070:   if (dimR == 2) {
3071:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

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

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

3084:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3085:     }
3086:   } else {
3087:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3088:   }
3089:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3090:   for (i = 0; i < dimR; i++) {
3091:     PetscReal *swap;

3093:     for (j = 0; j < (numV / 2); j++) {
3094:       for (k = 0; k < dimC; k++) {
3095:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3096:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3097:       }
3098:     }

3100:     if (i < dimR - 1) {
3101:       swap       = cellCoeffs;
3102:       cellCoeffs = cellCoords;
3103:       cellCoords = swap;
3104:     }
3105:   }
3106:   PetscCall(PetscArrayzero(refCoords, numPoints * dimR));
3107:   for (j = 0; j < numPoints; j++) {
3108:     for (i = 0; i < maxIts; i++) {
3109:       PetscReal *guess = &refCoords[dimR * j];

3111:       /* compute -residual and Jacobian */
3112:       for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3113:       for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3114:       for (k = 0; k < numV; k++) {
3115:         PetscReal extCoord = 1.;
3116:         for (l = 0; l < dimR; l++) {
3117:           PetscReal coord = guess[l];
3118:           PetscInt  dep   = (k & (1 << l)) >> l;

3120:           extCoord *= dep * coord + !dep;
3121:           extJ[l] = dep;

3123:           for (m = 0; m < dimR; m++) {
3124:             PetscReal coord = guess[m];
3125:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
3126:             PetscReal mult  = dep * coord + !dep;

3128:             extJ[l] *= mult;
3129:           }
3130:         }
3131:         for (l = 0; l < dimC; l++) {
3132:           PetscReal coeff = cellCoeffs[dimC * k + l];

3134:           resNeg[l] -= coeff * extCoord;
3135:           for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3136:         }
3137:       }
3138:       if (0 && PetscDefined(USE_DEBUG)) {
3139:         PetscReal maxAbs = 0.;

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

3145:       PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess));
3146:     }
3147:   }
3148:   PetscCall(DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J));
3149:   PetscCall(DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData));
3150:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3151:   PetscFunctionReturn(PETSC_SUCCESS);
3152: }

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

3160:   PetscFunctionBegin;
3162:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3163:   PetscCheck(coordSize >= dimC * numV, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expecting at least %" PetscInt_FMT " coordinates, got %" PetscInt_FMT, dimC * (1 << dimR), coordSize);
3164:   PetscCall(DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3165:   cellCoords = &cellData[0];
3166:   cellCoeffs = &cellData[coordSize];
3167:   if (dimR == 2) {
3168:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

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

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

3181:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3182:     }
3183:   } else {
3184:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3185:   }
3186:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3187:   for (i = 0; i < dimR; i++) {
3188:     PetscReal *swap;

3190:     for (j = 0; j < (numV / 2); j++) {
3191:       for (k = 0; k < dimC; k++) {
3192:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3193:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3194:       }
3195:     }

3197:     if (i < dimR - 1) {
3198:       swap       = cellCoeffs;
3199:       cellCoeffs = cellCoords;
3200:       cellCoords = swap;
3201:     }
3202:   }
3203:   PetscCall(PetscArrayzero(realCoords, numPoints * dimC));
3204:   for (j = 0; j < numPoints; j++) {
3205:     const PetscReal *guess  = &refCoords[dimR * j];
3206:     PetscReal       *mapped = &realCoords[dimC * j];

3208:     for (k = 0; k < numV; k++) {
3209:       PetscReal extCoord = 1.;
3210:       for (l = 0; l < dimR; l++) {
3211:         PetscReal coord = guess[l];
3212:         PetscInt  dep   = (k & (1 << l)) >> l;

3214:         extCoord *= dep * coord + !dep;
3215:       }
3216:       for (l = 0; l < dimC; l++) {
3217:         PetscReal coeff = cellCoeffs[dimC * k + l];

3219:         mapped[l] += coeff * extCoord;
3220:       }
3221:     }
3222:   }
3223:   PetscCall(DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData));
3224:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar));
3225:   PetscFunctionReturn(PETSC_SUCCESS);
3226: }

3228: /* TODO: TOBY please fix this for Nc > 1 */
3229: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3230: {
3231:   PetscInt     numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
3232:   PetscScalar *nodes = NULL;
3233:   PetscReal   *invV, *modes;
3234:   PetscReal   *B, *D, *resNeg;
3235:   PetscScalar *J, *invJ, *work;

3237:   PetscFunctionBegin;
3238:   PetscCall(PetscFEGetDimension(fe, &pdim));
3239:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
3240:   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);
3241:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3242:   /* convert nodes to values in the stable evaluation basis */
3243:   PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3244:   invV = fe->invV;
3245:   for (i = 0; i < pdim; ++i) {
3246:     modes[i] = 0.;
3247:     for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3248:   }
3249:   PetscCall(DMGetWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3250:   D      = &B[pdim * Nc];
3251:   resNeg = &D[pdim * Nc * dimR];
3252:   PetscCall(DMGetWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3253:   invJ = &J[Nc * dimR];
3254:   work = &invJ[Nc * dimR];
3255:   for (i = 0; i < numPoints * dimR; i++) refCoords[i] = 0.;
3256:   for (j = 0; j < numPoints; j++) {
3257:     for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
3258:       PetscReal *guess = &refCoords[j * dimR];
3259:       PetscCall(PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL));
3260:       for (k = 0; k < Nc; k++) resNeg[k] = realCoords[j * Nc + k];
3261:       for (k = 0; k < Nc * dimR; k++) J[k] = 0.;
3262:       for (k = 0; k < pdim; k++) {
3263:         for (l = 0; l < Nc; l++) {
3264:           resNeg[l] -= modes[k] * B[k * Nc + l];
3265:           for (m = 0; m < dimR; m++) J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
3266:         }
3267:       }
3268:       if (0 && PetscDefined(USE_DEBUG)) {
3269:         PetscReal maxAbs = 0.;

3271:         for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3272:         PetscCall(PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs));
3273:       }
3274:       PetscCall(DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess));
3275:     }
3276:   }
3277:   PetscCall(DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J));
3278:   PetscCall(DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B));
3279:   PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3280:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3281:   PetscFunctionReturn(PETSC_SUCCESS);
3282: }

3284: /* TODO: TOBY please fix this for Nc > 1 */
3285: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3286: {
3287:   PetscInt     numComp, pdim, i, j, k, l, coordSize;
3288:   PetscScalar *nodes = NULL;
3289:   PetscReal   *invV, *modes;
3290:   PetscReal   *B;

3292:   PetscFunctionBegin;
3293:   PetscCall(PetscFEGetDimension(fe, &pdim));
3294:   PetscCall(PetscFEGetNumComponents(fe, &numComp));
3295:   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);
3296:   PetscCall(DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3297:   /* convert nodes to values in the stable evaluation basis */
3298:   PetscCall(DMGetWorkArray(dm, pdim, MPIU_REAL, &modes));
3299:   invV = fe->invV;
3300:   for (i = 0; i < pdim; ++i) {
3301:     modes[i] = 0.;
3302:     for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3303:   }
3304:   PetscCall(DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3305:   PetscCall(PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL));
3306:   for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3307:   for (j = 0; j < numPoints; j++) {
3308:     PetscReal *mapped = &realCoords[j * Nc];

3310:     for (k = 0; k < pdim; k++) {
3311:       for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3312:     }
3313:   }
3314:   PetscCall(DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B));
3315:   PetscCall(DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes));
3316:   PetscCall(DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes));
3317:   PetscFunctionReturn(PETSC_SUCCESS);
3318: }

3320: /*@
3321:   DMPlexCoordinatesToReference - Pull coordinates back from the mesh to the reference element using a single element
3322:   map.  This inversion will be accurate inside the reference element, but may be inaccurate for mappings that do not
3323:   extend uniquely outside the reference cell (e.g, most non-affine maps)

3325:   Not Collective

3327:   Input Parameters:
3328: + dm         - The mesh, with coordinate maps defined either by a `PetscDS` for the coordinate `DM` (see `DMGetCoordinateDM()`) or
3329:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3330:                as a multilinear map for tensor-product elements
3331: . cell       - the cell whose map is used.
3332: . numPoints  - the number of points to locate
3333: - realCoords - (numPoints x coordinate dimension) array of coordinates (see `DMGetCoordinateDim()`)

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

3338:   Level: intermediate

3340: .seealso: `DMPLEX`, `DMPlexReferenceToCoordinates()`
3341: @*/
3342: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3343: {
3344:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3345:   DM       coordDM = NULL;
3346:   Vec      coords;
3347:   PetscFE  fe = NULL;

3349:   PetscFunctionBegin;
3351:   PetscCall(DMGetDimension(dm, &dimR));
3352:   PetscCall(DMGetCoordinateDim(dm, &dimC));
3353:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3354:   PetscCall(DMPlexGetDepth(dm, &depth));
3355:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
3356:   PetscCall(DMGetCoordinateDM(dm, &coordDM));
3357:   if (coordDM) {
3358:     PetscInt coordFields;

3360:     PetscCall(DMGetNumFields(coordDM, &coordFields));
3361:     if (coordFields) {
3362:       PetscClassId id;
3363:       PetscObject  disc;

3365:       PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3366:       PetscCall(PetscObjectGetClassId(disc, &id));
3367:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3368:     }
3369:   }
3370:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3371:   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);
3372:   if (!fe) { /* implicit discretization: affine or multilinear */
3373:     PetscInt  coneSize;
3374:     PetscBool isSimplex, isTensor;

3376:     PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3377:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3378:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3379:     if (isSimplex) {
3380:       PetscReal detJ, *v0, *J, *invJ;

3382:       PetscCall(DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3383:       J    = &v0[dimC];
3384:       invJ = &J[dimC * dimC];
3385:       PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ));
3386:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3387:         const PetscReal x0[3] = {-1., -1., -1.};

3389:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3390:       }
3391:       PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3392:     } else if (isTensor) {
3393:       PetscCall(DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3394:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3395:   } else {
3396:     PetscCall(DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR));
3397:   }
3398:   PetscFunctionReturn(PETSC_SUCCESS);
3399: }

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

3404:   Not Collective

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

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

3417:    Level: intermediate

3419: .seealso: `DMPLEX`, `DMPlexCoordinatesToReference()`
3420: @*/
3421: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3422: {
3423:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3424:   DM       coordDM = NULL;
3425:   Vec      coords;
3426:   PetscFE  fe = NULL;

3428:   PetscFunctionBegin;
3430:   PetscCall(DMGetDimension(dm, &dimR));
3431:   PetscCall(DMGetCoordinateDim(dm, &dimC));
3432:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) PetscFunctionReturn(PETSC_SUCCESS);
3433:   PetscCall(DMPlexGetDepth(dm, &depth));
3434:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
3435:   PetscCall(DMGetCoordinateDM(dm, &coordDM));
3436:   if (coordDM) {
3437:     PetscInt coordFields;

3439:     PetscCall(DMGetNumFields(coordDM, &coordFields));
3440:     if (coordFields) {
3441:       PetscClassId id;
3442:       PetscObject  disc;

3444:       PetscCall(DMGetField(coordDM, 0, NULL, &disc));
3445:       PetscCall(PetscObjectGetClassId(disc, &id));
3446:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3447:     }
3448:   }
3449:   PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
3450:   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);
3451:   if (!fe) { /* implicit discretization: affine or multilinear */
3452:     PetscInt  coneSize;
3453:     PetscBool isSimplex, isTensor;

3455:     PetscCall(DMPlexGetConeSize(dm, cell, &coneSize));
3456:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3457:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3458:     if (isSimplex) {
3459:       PetscReal detJ, *v0, *J;

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

3467:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3468:       }
3469:       PetscCall(DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0));
3470:     } else if (isTensor) {
3471:       PetscCall(DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3472:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3473:   } else {
3474:     PetscCall(DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR));
3475:   }
3476:   PetscFunctionReturn(PETSC_SUCCESS);
3477: }

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

3482:   Not Collective

3484:   Input Parameters:
3485: + dm      - The `DM`
3486: . time    - The time
3487: - func    - The function transforming current coordinates to new coordaintes

3489:    Calling sequence of `func`:
3490: .vb
3491:    void func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3492:              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3493:              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3494:              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);
3495: .ve
3496: +  dim          - The spatial dimension
3497: .  Nf           - The number of input fields (here 1)
3498: .  NfAux        - The number of input auxiliary fields
3499: .  uOff         - The offset of the coordinates in u[] (here 0)
3500: .  uOff_x       - The offset of the coordinates in u_x[] (here 0)
3501: .  u            - The coordinate values at this point in space
3502: .  u_t          - The coordinate time derivative at this point in space (here `NULL`)
3503: .  u_x          - The coordinate derivatives at this point in space
3504: .  aOff         - The offset of each auxiliary field in u[]
3505: .  aOff_x       - The offset of each auxiliary field in u_x[]
3506: .  a            - The auxiliary field values at this point in space
3507: .  a_t          - The auxiliary field time derivative at this point in space (or `NULL`)
3508: .  a_x          - The auxiliary field derivatives at this point in space
3509: .  t            - The current time
3510: .  x            - The coordinates of this point (here not used)
3511: .  numConstants - The number of constants
3512: .  constants    - The value of each constant
3513: -  f            - The new coordinates at this point in space

3515:   Level: intermediate

3517: .seealso: `DMPLEX`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
3518: @*/
3519: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time, void (*func)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3520: {
3521:   DM      cdm;
3522:   DMField cf;
3523:   Vec     lCoords, tmpCoords;

3525:   PetscFunctionBegin;
3526:   PetscCall(DMGetCoordinateDM(dm, &cdm));
3527:   PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
3528:   PetscCall(DMGetLocalVector(cdm, &tmpCoords));
3529:   PetscCall(VecCopy(lCoords, tmpCoords));
3530:   /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3531:   PetscCall(DMGetCoordinateField(dm, &cf));
3532:   cdm->coordinates[0].field = cf;
3533:   PetscCall(DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords));
3534:   cdm->coordinates[0].field = NULL;
3535:   PetscCall(DMRestoreLocalVector(cdm, &tmpCoords));
3536:   PetscCall(DMSetCoordinatesLocal(dm, lCoords));
3537:   PetscFunctionReturn(PETSC_SUCCESS);
3538: }

3540: /* Shear applies the transformation, assuming we fix z,
3541:   / 1  0  m_0 \
3542:   | 0  1  m_1 |
3543:   \ 0  0   1  /
3544: */
3545: static void f0_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[])
3546: {
3547:   const PetscInt Nc = uOff[1] - uOff[0];
3548:   const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
3549:   PetscInt       c;

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

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

3557:   Not Collective

3559:   Input Parameters:
3560: + dm          - The `DMPLEX`
3561: . direction   - The shear coordinate direction, e.g. 0 is the x-axis
3562: - multipliers - The multiplier m for each direction which is not the shear direction

3564:   Level: intermediate

3566: .seealso: `DMPLEX`, `DMPlexRemapGeometry()`
3567: @*/
3568: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3569: {
3570:   DM             cdm;
3571:   PetscDS        cds;
3572:   PetscObject    obj;
3573:   PetscClassId   id;
3574:   PetscScalar   *moduli;
3575:   const PetscInt dir = (PetscInt)direction;
3576:   PetscInt       dE, d, e;

3578:   PetscFunctionBegin;
3579:   PetscCall(DMGetCoordinateDM(dm, &cdm));
3580:   PetscCall(DMGetCoordinateDim(dm, &dE));
3581:   PetscCall(PetscMalloc1(dE + 1, &moduli));
3582:   moduli[0] = dir;
3583:   for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3584:   PetscCall(DMGetDS(cdm, &cds));
3585:   PetscCall(PetscDSGetDiscretization(cds, 0, &obj));
3586:   PetscCall(PetscObjectGetClassId(obj, &id));
3587:   if (id != PETSCFE_CLASSID) {
3588:     Vec          lCoords;
3589:     PetscSection cSection;
3590:     PetscScalar *coords;
3591:     PetscInt     vStart, vEnd, v;

3593:     PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
3594:     PetscCall(DMGetCoordinateSection(dm, &cSection));
3595:     PetscCall(DMGetCoordinatesLocal(dm, &lCoords));
3596:     PetscCall(VecGetArray(lCoords, &coords));
3597:     for (v = vStart; v < vEnd; ++v) {
3598:       PetscReal ds;
3599:       PetscInt  off, c;

3601:       PetscCall(PetscSectionGetOffset(cSection, v, &off));
3602:       ds = PetscRealPart(coords[off + dir]);
3603:       for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds;
3604:     }
3605:     PetscCall(VecRestoreArray(lCoords, &coords));
3606:   } else {
3607:     PetscCall(PetscDSSetConstants(cds, dE + 1, moduli));
3608:     PetscCall(DMPlexRemapGeometry(dm, 0.0, f0_shear));
3609:   }
3610:   PetscCall(PetscFree(moduli));
3611:   PetscFunctionReturn(PETSC_SUCCESS);
3612: }