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 called already)

 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 Parameters:
 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: `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:   if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
 46:   DMGetCoordinateDim(dm, &cdim);
 47:   {
 48:     PetscInt n;

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

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

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

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

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

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

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

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

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

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

205:   DMPlexComputeCellGeometryFEM(dm, c, NULL, &v0, &J, &invJ, &detJ);
206:   xi = invJ * (x - v0);

208:   if ((xi >= -eps) && (xi <= 2. + eps)) *cell = c;
209:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
210:   return 0;
211: }

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

222:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
223:   xi  = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
224:   eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);

226:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0 + eps)) *cell = c;
227:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
228:   return 0;
229: }

231: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
232: {
233:   const PetscInt embedDim = 2;
234:   PetscReal      x        = PetscRealPart(point[0]);
235:   PetscReal      y        = PetscRealPart(point[1]);
236:   PetscReal      v0[2], J[4], invJ[4], detJ;
237:   PetscReal      xi, eta, r;

239:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
240:   xi  = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]);
241:   eta = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]);

243:   xi  = PetscMax(xi, 0.0);
244:   eta = PetscMax(eta, 0.0);
245:   if (xi + eta > 2.0) {
246:     r = (xi + eta) / 2.0;
247:     xi /= r;
248:     eta /= r;
249:   }
250:   cpoint[0] = J[0 * embedDim + 0] * xi + J[0 * embedDim + 1] * eta + v0[0];
251:   cpoint[1] = J[1 * embedDim + 0] * xi + J[1 * embedDim + 1] * eta + v0[1];
252:   return 0;
253: }

255: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
256: {
257:   const PetscScalar *array;
258:   PetscScalar       *coords    = NULL;
259:   const PetscInt     faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
260:   PetscReal          x         = PetscRealPart(point[0]);
261:   PetscReal          y         = PetscRealPart(point[1]);
262:   PetscInt           crossings = 0, numCoords, f;
263:   PetscBool          isDG;

265:   DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords);
267:   for (f = 0; f < 4; ++f) {
268:     PetscReal x_i   = PetscRealPart(coords[faces[2 * f + 0] * 2 + 0]);
269:     PetscReal y_i   = PetscRealPart(coords[faces[2 * f + 0] * 2 + 1]);
270:     PetscReal x_j   = PetscRealPart(coords[faces[2 * f + 1] * 2 + 0]);
271:     PetscReal y_j   = PetscRealPart(coords[faces[2 * f + 1] * 2 + 1]);
272:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
273:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
274:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
275:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
276:     if ((cond1 || cond2) && above) ++crossings;
277:   }
278:   if (crossings % 2) *cell = c;
279:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
280:   DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords);
281:   return 0;
282: }

284: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
285: {
286:   const PetscInt  embedDim = 3;
287:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
288:   PetscReal       v0[3], J[9], invJ[9], detJ;
289:   PetscReal       x = PetscRealPart(point[0]);
290:   PetscReal       y = PetscRealPart(point[1]);
291:   PetscReal       z = PetscRealPart(point[2]);
292:   PetscReal       xi, eta, zeta;

294:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
295:   xi   = invJ[0 * embedDim + 0] * (x - v0[0]) + invJ[0 * embedDim + 1] * (y - v0[1]) + invJ[0 * embedDim + 2] * (z - v0[2]);
296:   eta  = invJ[1 * embedDim + 0] * (x - v0[0]) + invJ[1 * embedDim + 1] * (y - v0[1]) + invJ[1 * embedDim + 2] * (z - v0[2]);
297:   zeta = invJ[2 * embedDim + 0] * (x - v0[0]) + invJ[2 * embedDim + 1] * (y - v0[1]) + invJ[2 * embedDim + 2] * (z - v0[2]);

299:   if ((xi >= -eps) && (eta >= -eps) && (zeta >= -eps) && (xi + eta + zeta <= 2.0 + eps)) *cell = c;
300:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
301:   return 0;
302: }

304: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
305: {
306:   const PetscScalar *array;
307:   PetscScalar       *coords    = NULL;
308:   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};
309:   PetscBool          found     = PETSC_TRUE;
310:   PetscInt           numCoords, f;
311:   PetscBool          isDG;

313:   DMPlexGetCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords);
315:   for (f = 0; f < 6; ++f) {
316:     /* Check the point is under plane */
317:     /*   Get face normal */
318:     PetscReal v_i[3];
319:     PetscReal v_j[3];
320:     PetscReal normal[3];
321:     PetscReal pp[3];
322:     PetscReal dot;

324:     v_i[0]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
325:     v_i[1]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
326:     v_i[2]    = PetscRealPart(coords[faces[f * 4 + 3] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
327:     v_j[0]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 0] - coords[faces[f * 4 + 0] * 3 + 0]);
328:     v_j[1]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 1] - coords[faces[f * 4 + 0] * 3 + 1]);
329:     v_j[2]    = PetscRealPart(coords[faces[f * 4 + 1] * 3 + 2] - coords[faces[f * 4 + 0] * 3 + 2]);
330:     normal[0] = v_i[1] * v_j[2] - v_i[2] * v_j[1];
331:     normal[1] = v_i[2] * v_j[0] - v_i[0] * v_j[2];
332:     normal[2] = v_i[0] * v_j[1] - v_i[1] * v_j[0];
333:     pp[0]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 0] - point[0]);
334:     pp[1]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 1] - point[1]);
335:     pp[2]     = PetscRealPart(coords[faces[f * 4 + 0] * 3 + 2] - point[2]);
336:     dot       = normal[0] * pp[0] + normal[1] * pp[1] + normal[2] * pp[2];

338:     /* Check that projected point is in face (2D location problem) */
339:     if (dot < 0.0) {
340:       found = PETSC_FALSE;
341:       break;
342:     }
343:   }
344:   if (found) *cell = c;
345:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
346:   DMPlexRestoreCellCoordinates(dm, c, &isDG, &numCoords, &array, &coords);
347:   return 0;
348: }

350: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
351: {
352:   PetscInt d;

354:   box->dim = dim;
355:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
356:   return 0;
357: }

359: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
360: {
361:   PetscMalloc1(1, box);
362:   PetscGridHashInitialize_Internal(*box, dim, point);
363:   return 0;
364: }

366: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
367: {
368:   PetscInt d;

370:   for (d = 0; d < box->dim; ++d) {
371:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
372:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
373:   }
374:   return 0;
375: }

377: /*
378:   PetscGridHashSetGrid - Divide the grid into boxes

380:   Not collective

382:   Input Parameters:
383: + box - The grid hash object
384: . n   - The number of boxes in each dimension, or PETSC_DETERMINE
385: - h   - The box size in each dimension, only used if n[d] == PETSC_DETERMINE

387:   Level: developer

389: .seealso: `PetscGridHashCreate()`
390: */
391: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
392: {
393:   PetscInt d;

395:   for (d = 0; d < box->dim; ++d) {
396:     box->extent[d] = box->upper[d] - box->lower[d];
397:     if (n[d] == PETSC_DETERMINE) {
398:       box->h[d] = h[d];
399:       box->n[d] = PetscCeilReal(box->extent[d] / h[d]);
400:     } else {
401:       box->n[d] = n[d];
402:       box->h[d] = box->extent[d] / n[d];
403:     }
404:   }
405:   return 0;
406: }

408: /*
409:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

411:   Not collective

413:   Input Parameters:
414: + box       - The grid hash object
415: . numPoints - The number of input points
416: - points    - The input point coordinates

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

422:   Level: developer

424:   Note:
425:   This only guarantees that a box contains a point, not that a cell does.

427: .seealso: `PetscGridHashCreate()`
428: */
429: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
430: {
431:   const PetscReal *lower = box->lower;
432:   const PetscReal *upper = box->upper;
433:   const PetscReal *h     = box->h;
434:   const PetscInt  *n     = box->n;
435:   const PetscInt   dim   = box->dim;
436:   PetscInt         d, p;

438:   for (p = 0; p < numPoints; ++p) {
439:     for (d = 0; d < dim; ++d) {
440:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

442:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
443:       if (dbox == -1 && PetscAbsReal(PetscRealPart(points[p * dim + d]) - lower[d]) < 1.0e-9) dbox = 0;
445:       dboxes[p * dim + d] = dbox;
446:     }
447:     if (boxes)
448:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
449:   }
450:   return 0;
451: }

453: /*
454:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

456:  Not collective

458:   Input Parameters:
459: + box         - The grid hash object
460: . cellSection - The PetscSection mapping cells to boxes
461: . numPoints   - The number of input points
462: - points      - The input point coordinates

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

469:   Level: developer

471:   Note:
472:   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 the cellSection is already constructed.

474: .seealso: `PetscGridHashGetEnclosingBox()`
475: */
476: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscSection cellSection, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[], PetscBool *found)
477: {
478:   const PetscReal *lower = box->lower;
479:   const PetscReal *upper = box->upper;
480:   const PetscReal *h     = box->h;
481:   const PetscInt  *n     = box->n;
482:   const PetscInt   dim   = box->dim;
483:   PetscInt         bStart, bEnd, d, p;

486:   *found = PETSC_FALSE;
487:   PetscSectionGetChart(box->cellSection, &bStart, &bEnd);
488:   for (p = 0; p < numPoints; ++p) {
489:     for (d = 0; d < dim; ++d) {
490:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p * dim + d]) - lower[d]) / h[d]);

492:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p * dim + d]) - upper[d]) < 1.0e-9) dbox = n[d] - 1;
493:       if (dbox < 0 || dbox >= n[d]) return 0;
494:       dboxes[p * dim + d] = dbox;
495:     }
496:     if (boxes)
497:       for (d = dim - 2, boxes[p] = dboxes[p * dim + dim - 1]; d >= 0; --d) boxes[p] = boxes[p] * n[d] + dboxes[p * dim + d];
498:     // It is possible for a box to overlap no grid cells
499:     if (boxes[p] < bStart || boxes[p] >= bEnd) return 0;
500:   }
501:   *found = PETSC_TRUE;
502:   return 0;
503: }

505: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
506: {
507:   if (*box) {
508:     PetscSectionDestroy(&(*box)->cellSection);
509:     ISDestroy(&(*box)->cells);
510:     DMLabelDestroy(&(*box)->cellsSparse);
511:   }
512:   PetscFree(*box);
513:   return 0;
514: }

516: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
517: {
518:   DMPolytopeType ct;

520:   DMPlexGetCellType(dm, cellStart, &ct);
521:   switch (ct) {
522:   case DM_POLYTOPE_SEGMENT:
523:     DMPlexLocatePoint_Simplex_1D_Internal(dm, point, cellStart, cell);
524:     break;
525:   case DM_POLYTOPE_TRIANGLE:
526:     DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
527:     break;
528:   case DM_POLYTOPE_QUADRILATERAL:
529:     DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell);
530:     break;
531:   case DM_POLYTOPE_TETRAHEDRON:
532:     DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
533:     break;
534:   case DM_POLYTOPE_HEXAHEDRON:
535:     DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
536:     break;
537:   default:
538:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %" PetscInt_FMT " with type %s", cellStart, DMPolytopeTypes[ct]);
539:   }
540:   return 0;
541: }

543: /*
544:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
545: */
546: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
547: {
548:   DMPolytopeType ct;

550:   DMPlexGetCellType(dm, cell, &ct);
551:   switch (ct) {
552:   case DM_POLYTOPE_TRIANGLE:
553:     DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);
554:     break;
555: #if 0
556:     case DM_POLYTOPE_QUADRILATERAL:
557:     DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);break;
558:     case DM_POLYTOPE_TETRAHEDRON:
559:     DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);break;
560:     case DM_POLYTOPE_HEXAHEDRON:
561:     DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);break;
562: #endif
563:   default:
564:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %" PetscInt_FMT " with type %s", cell, DMPolytopeTypes[ct]);
565:   }
566:   return 0;
567: }

569: /*
570:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex

572:   Collective on dm

574:   Input Parameter:
575: . dm - The Plex

577:   Output Parameter:
578: . localBox - The grid hash object

580:   Level: developer

582: .seealso: `PetscGridHashCreate()`, `PetscGridHashGetEnclosingBox()`
583: */
584: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
585: {
586:   PetscInt           debug = ((DM_Plex *)dm->data)->printLocate;
587:   MPI_Comm           comm;
588:   PetscGridHash      lbox;
589:   Vec                coordinates;
590:   PetscSection       coordSection;
591:   Vec                coordsLocal;
592:   const PetscScalar *coords;
593:   PetscScalar       *edgeCoords;
594:   PetscInt          *dboxes, *boxes;
595:   PetscInt           n[3] = {2, 2, 2};
596:   PetscInt           dim, N, maxConeSize, cStart, cEnd, c, eStart, eEnd, i;
597:   PetscBool          flg;

599:   PetscObjectGetComm((PetscObject)dm, &comm);
600:   DMGetCoordinatesLocal(dm, &coordinates);
601:   DMGetCoordinateDim(dm, &dim);
602:   DMPlexGetMaxSizes(dm, &maxConeSize, NULL);
603:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
604:   VecGetLocalSize(coordinates, &N);
605:   VecGetArrayRead(coordinates, &coords);
606:   PetscGridHashCreate(comm, dim, coords, &lbox);
607:   for (i = 0; i < N; i += dim) PetscGridHashEnlarge(lbox, &coords[i]);
608:   VecRestoreArrayRead(coordinates, &coords);
609:   c = dim;
610:   PetscOptionsGetIntArray(NULL, ((PetscObject)dm)->prefix, "-dm_plex_hash_box_faces", n, &c, &flg);
611:   if (flg) {
612:     for (i = c; i < dim; ++i) n[i] = n[c - 1];
613:   } else {
614:     for (i = 0; i < dim; ++i) n[i] = PetscMax(2, PetscFloorReal(PetscPowReal((PetscReal)(cEnd - cStart), 1.0 / dim) * 0.8));
615:   }
616:   PetscGridHashSetGrid(lbox, n, NULL);
617:   if (debug)
618:     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],
619:                           (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]));
620: #if 0
621:   /* Could define a custom reduction to merge these */
622:   MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
623:   MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
624: #endif
625:   /* Is there a reason to snap the local bounding box to a division of the global box? */
626:   /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
627:   /* Create label */
628:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
629:   if (dim < 2) eStart = eEnd = -1;
630:   DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
631:   DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
632:   /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
633:   DMGetCoordinatesLocal(dm, &coordsLocal);
634:   DMGetCoordinateSection(dm, &coordSection);
635:   PetscCalloc3(16 * dim, &dboxes, 16, &boxes, PetscPowInt(maxConeSize, dim) * dim, &edgeCoords);
636:   for (c = cStart; c < cEnd; ++c) {
637:     const PetscReal *h       = lbox->h;
638:     PetscScalar     *ccoords = NULL;
639:     PetscInt         csize   = 0;
640:     PetscInt        *closure = NULL;
641:     PetscInt         Ncl, cl, Ne = 0;
642:     PetscScalar      point[3];
643:     PetscInt         dlim[6], d, e, i, j, k;

645:     /* Get all edges in cell */
646:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
647:     for (cl = 0; cl < Ncl * 2; ++cl) {
648:       if ((closure[cl] >= eStart) && (closure[cl] < eEnd)) {
649:         PetscScalar *ecoords = &edgeCoords[Ne * dim * 2];
650:         PetscInt     ecsize  = dim * 2;

652:         DMPlexVecGetClosure(dm, coordSection, coordsLocal, closure[cl], &ecsize, &ecoords);
654:         ++Ne;
655:       }
656:     }
657:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &Ncl, &closure);
658:     /* Find boxes enclosing each vertex */
659:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
660:     PetscGridHashGetEnclosingBox(lbox, csize / dim, ccoords, dboxes, boxes);
661:     /* Mark cells containing the vertices */
662:     for (e = 0; e < csize / dim; ++e) {
663:       if (debug)
664:         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);
665:       DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);
666:     }
667:     /* Get grid of boxes containing these */
668:     for (d = 0; d < dim; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = dboxes[d];
669:     for (d = dim; d < 3; ++d) dlim[d * 2 + 0] = dlim[d * 2 + 1] = 0;
670:     for (e = 1; e < dim + 1; ++e) {
671:       for (d = 0; d < dim; ++d) {
672:         dlim[d * 2 + 0] = PetscMin(dlim[d * 2 + 0], dboxes[e * dim + d]);
673:         dlim[d * 2 + 1] = PetscMax(dlim[d * 2 + 1], dboxes[e * dim + d]);
674:       }
675:     }
676:     /* Check for intersection of box with cell */
677:     for (k = dlim[2 * 2 + 0], point[2] = lbox->lower[2] + k * h[2]; k <= dlim[2 * 2 + 1]; ++k, point[2] += h[2]) {
678:       for (j = dlim[1 * 2 + 0], point[1] = lbox->lower[1] + j * h[1]; j <= dlim[1 * 2 + 1]; ++j, point[1] += h[1]) {
679:         for (i = dlim[0 * 2 + 0], point[0] = lbox->lower[0] + i * h[0]; i <= dlim[0 * 2 + 1]; ++i, point[0] += h[0]) {
680:           const PetscInt box = (k * lbox->n[1] + j) * lbox->n[0] + i;
681:           PetscScalar    cpoint[3];
682:           PetscInt       cell, edge, ii, jj, kk;

684:           if (debug)
685:             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]));
686:           /* Check whether cell contains any vertex of this subbox TODO vectorize this */
687:           for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
688:             for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
689:               for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {
690:                 DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
691:                 if (cell >= 0) {
692:                   if (debug)
693:                     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);
694:                   DMLabelSetValue(lbox->cellsSparse, c, box);
695:                   jj = kk = 2;
696:                   break;
697:                 }
698:               }
699:             }
700:           }
701:           /* Check whether cell edge intersects any face of these subboxes TODO vectorize this */
702:           for (edge = 0; edge < Ne; ++edge) {
703:             PetscReal segA[6] = {0., 0., 0., 0., 0., 0.};
704:             PetscReal segB[6] = {0., 0., 0., 0., 0., 0.};
705:             PetscReal segC[6] = {0., 0., 0., 0., 0., 0.};

708:             for (d = 0; d < dim * 2; ++d) segA[d] = PetscRealPart(edgeCoords[edge * dim * 2 + d]);
709:             /* 1D: (x) -- (x+h)               0 -- 1
710:                2D: (x,   y)   -- (x,   y+h)   (0, 0) -- (0, 1)
711:                    (x+h, y)   -- (x+h, y+h)   (1, 0) -- (1, 1)
712:                    (x,   y)   -- (x+h, y)     (0, 0) -- (1, 0)
713:                    (x,   y+h) -- (x+h, y+h)   (0, 1) -- (1, 1)
714:                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)
715:                    (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)
716:                    (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)
717:                    (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)
718:                    (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)
719:                    (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)
720:              */
721:             /* Loop over faces with normal in direction d */
722:             for (d = 0; d < dim; ++d) {
723:               PetscBool intersects = PETSC_FALSE;
724:               PetscInt  e          = (d + 1) % dim;
725:               PetscInt  f          = (d + 2) % dim;

727:               /* There are two faces in each dimension */
728:               for (ii = 0; ii < 2; ++ii) {
729:                 segB[d]       = PetscRealPart(point[d] + ii * h[d]);
730:                 segB[dim + d] = PetscRealPart(point[d] + ii * h[d]);
731:                 segC[d]       = PetscRealPart(point[d] + ii * h[d]);
732:                 segC[dim + d] = PetscRealPart(point[d] + ii * h[d]);
733:                 if (dim > 1) {
734:                   segB[e]       = PetscRealPart(point[e] + 0 * h[e]);
735:                   segB[dim + e] = PetscRealPart(point[e] + 1 * h[e]);
736:                   segC[e]       = PetscRealPart(point[e] + 0 * h[e]);
737:                   segC[dim + e] = PetscRealPart(point[e] + 0 * h[e]);
738:                 }
739:                 if (dim > 2) {
740:                   segB[f]       = PetscRealPart(point[f] + 0 * h[f]);
741:                   segB[dim + f] = PetscRealPart(point[f] + 0 * h[f]);
742:                   segC[f]       = PetscRealPart(point[f] + 0 * h[f]);
743:                   segC[dim + f] = PetscRealPart(point[f] + 1 * h[f]);
744:                 }
745:                 if (dim == 2) {
746:                   DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
747:                 } else if (dim == 3) {
748:                   DMPlexGetLinePlaneIntersection_3D_Internal(segA, segB, segC, NULL, &intersects);
749:                 }
750:                 if (intersects) {
751:                   if (debug)
752:                     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]);
753:                   DMLabelSetValue(lbox->cellsSparse, c, box);
754:                   edge = Ne;
755:                   break;
756:                 }
757:               }
758:             }
759:           }
760:         }
761:       }
762:     }
763:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
764:   }
765:   PetscFree3(dboxes, boxes, edgeCoords);
766:   if (debug) DMLabelView(lbox->cellsSparse, PETSC_VIEWER_STDOUT_SELF);
767:   DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
768:   DMLabelDestroy(&lbox->cellsSparse);
769:   *localBox = lbox;
770:   return 0;
771: }

773: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
774: {
775:   PetscInt        debug = ((DM_Plex *)dm->data)->printLocate;
776:   DM_Plex        *mesh  = (DM_Plex *)dm->data;
777:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
778:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
779:   PetscInt        dim, Nl = 0, cStart, cEnd, numCells, c, d;
780:   PetscSF         sf;
781:   const PetscInt *leaves;
782:   const PetscInt *boxCells;
783:   PetscSFNode    *cells;
784:   PetscScalar    *a;
785:   PetscMPIInt     result;
786:   PetscLogDouble  t0, t1;
787:   PetscReal       gmin[3], gmax[3];
788:   PetscInt        terminating_query_type[] = {0, 0, 0};

790:   PetscLogEventBegin(DMPLEX_LocatePoints, 0, 0, 0, 0);
791:   PetscTime(&t0);
793:   DMGetCoordinateDim(dm, &dim);
794:   VecGetBlockSize(v, &bs);
795:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF), PETSC_COMM_SELF, &result);
798:   DMGetCoordinatesLocalSetUp(dm);
799:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
800:   DMGetPointSF(dm, &sf);
801:   if (sf) PetscSFGetGraph(sf, NULL, &Nl, &leaves, NULL);
802:   Nl = PetscMax(Nl, 0);
803:   VecGetLocalSize(v, &numPoints);
804:   VecGetArray(v, &a);
805:   numPoints /= bs;
806:   {
807:     const PetscSFNode *sf_cells;

809:     PetscSFGetGraph(cellSF, NULL, NULL, NULL, &sf_cells);
810:     if (sf_cells) {
811:       PetscInfo(dm, "[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
812:       cells = (PetscSFNode *)sf_cells;
813:       reuse = PETSC_TRUE;
814:     } else {
815:       PetscInfo(dm, "[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
816:       PetscMalloc1(numPoints, &cells);
817:       /* initialize cells if created */
818:       for (p = 0; p < numPoints; p++) {
819:         cells[p].rank  = 0;
820:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
821:       }
822:     }
823:   }
824:   DMGetBoundingBox(dm, gmin, gmax);
825:   if (hash) {
826:     if (!mesh->lbox) {
827:       PetscInfo(dm, "Initializing grid hashing");
828:       DMPlexComputeGridHash_Internal(dm, &mesh->lbox);
829:     }
830:     /* Designate the local box for each point */
831:     /* Send points to correct process */
832:     /* Search cells that lie in each subbox */
833:     /*   Should we bin points before doing search? */
834:     ISGetIndices(mesh->lbox->cells, &boxCells);
835:   }
836:   for (p = 0, numFound = 0; p < numPoints; ++p) {
837:     const PetscScalar *point   = &a[p * bs];
838:     PetscInt           dbin[3] = {-1, -1, -1}, bin, cell = -1, cellOffset;
839:     PetscBool          point_outside_domain = PETSC_FALSE;

841:     /* check bounding box of domain */
842:     for (d = 0; d < dim; d++) {
843:       if (PetscRealPart(point[d]) < gmin[d]) {
844:         point_outside_domain = PETSC_TRUE;
845:         break;
846:       }
847:       if (PetscRealPart(point[d]) > gmax[d]) {
848:         point_outside_domain = PETSC_TRUE;
849:         break;
850:       }
851:     }
852:     if (point_outside_domain) {
853:       cells[p].rank  = 0;
854:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
855:       terminating_query_type[0]++;
856:       continue;
857:     }

859:     /* check initial values in cells[].index - abort early if found */
860:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
861:       c              = cells[p].index;
862:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
863:       DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
864:       if (cell >= 0) {
865:         cells[p].rank  = 0;
866:         cells[p].index = cell;
867:         numFound++;
868:       }
869:     }
870:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
871:       terminating_query_type[1]++;
872:       continue;
873:     }

875:     if (hash) {
876:       PetscBool found_box;

878:       if (debug) 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]));
879:       /* allow for case that point is outside box - abort early */
880:       PetscGridHashGetEnclosingBoxQuery(mesh->lbox, mesh->lbox->cellSection, 1, point, dbin, &bin, &found_box);
881:       if (found_box) {
882:         if (debug) PetscPrintf(PETSC_COMM_SELF, "  Found point in box %" PetscInt_FMT " (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", bin, dbin[0], dbin[1], dbin[2]);
883:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
884:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
885:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
886:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
887:           if (debug) PetscPrintf(PETSC_COMM_SELF, "    Checking for point in cell %" PetscInt_FMT "\n", boxCells[c]);
888:           DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
889:           if (cell >= 0) {
890:             if (debug) PetscPrintf(PETSC_COMM_SELF, "      FOUND in cell %" PetscInt_FMT "\n", cell);
891:             cells[p].rank  = 0;
892:             cells[p].index = cell;
893:             numFound++;
894:             terminating_query_type[2]++;
895:             break;
896:           }
897:         }
898:       }
899:     } else {
900:       for (c = cStart; c < cEnd; ++c) {
901:         PetscInt idx;

903:         PetscFindInt(c, Nl, leaves, &idx);
904:         if (idx >= 0) continue;
905:         DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
906:         if (cell >= 0) {
907:           cells[p].rank  = 0;
908:           cells[p].index = cell;
909:           numFound++;
910:           terminating_query_type[2]++;
911:           break;
912:         }
913:       }
914:     }
915:   }
916:   if (hash) ISRestoreIndices(mesh->lbox->cells, &boxCells);
917:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
918:     for (p = 0; p < numPoints; p++) {
919:       const PetscScalar *point = &a[p * bs];
920:       PetscReal          cpoint[3], diff[3], best[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL}, dist, distMax = PETSC_MAX_REAL;
921:       PetscInt           dbin[3] = {-1, -1, -1}, bin, cellOffset, d, bestc = -1;

923:       if (cells[p].index < 0) {
924:         PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
925:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
926:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
927:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
928:           DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
929:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
930:           dist = DMPlex_NormD_Internal(dim, diff);
931:           if (dist < distMax) {
932:             for (d = 0; d < dim; ++d) best[d] = cpoint[d];
933:             bestc   = boxCells[c];
934:             distMax = dist;
935:           }
936:         }
937:         if (distMax < PETSC_MAX_REAL) {
938:           ++numFound;
939:           cells[p].rank  = 0;
940:           cells[p].index = bestc;
941:           for (d = 0; d < dim; ++d) a[p * bs + d] = best[d];
942:         }
943:       }
944:     }
945:   }
946:   /* This code is only be relevant when interfaced to parallel point location */
947:   /* Check for highest numbered proc that claims a point (do we care?) */
948:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
949:     PetscMalloc1(numFound, &found);
950:     for (p = 0, numFound = 0; p < numPoints; p++) {
951:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
952:         if (numFound < p) cells[numFound] = cells[p];
953:         found[numFound++] = p;
954:       }
955:     }
956:   }
957:   VecRestoreArray(v, &a);
958:   if (!reuse) PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
959:   PetscTime(&t1);
960:   if (hash) {
961:     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]);
962:   } else {
963:     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]);
964:   }
965:   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)));
966:   PetscLogEventEnd(DMPLEX_LocatePoints, 0, 0, 0, 0);
967:   return 0;
968: }

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

973:   Not collective

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

978:   Output Parameter:
979: . R - The rotation which accomplishes the projection

981:   Level: developer

983: .seealso: `DMPlexComputeProjection3Dto1D()`, `DMPlexComputeProjection3Dto2D()`
984: @*/
985: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
986: {
987:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
988:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
989:   const PetscReal r = PetscSqrtReal(x * x + y * y), c = x / r, s = y / r;

991:   R[0]      = c;
992:   R[1]      = -s;
993:   R[2]      = s;
994:   R[3]      = c;
995:   coords[0] = 0.0;
996:   coords[1] = r;
997:   return 0;
998: }

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

1003:   Not collective

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

1008:   Output Parameter:
1009: . R - The rotation which accomplishes the projection

1011:   Note: 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

1013:   Level: developer

1015: .seealso: `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto2D()`
1016: @*/
1017: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
1018: {
1019:   PetscReal x    = PetscRealPart(coords[3] - coords[0]);
1020:   PetscReal y    = PetscRealPart(coords[4] - coords[1]);
1021:   PetscReal z    = PetscRealPart(coords[5] - coords[2]);
1022:   PetscReal r    = PetscSqrtReal(x * x + y * y + z * z);
1023:   PetscReal rinv = 1. / r;

1025:   x *= rinv;
1026:   y *= rinv;
1027:   z *= rinv;
1028:   if (x > 0.) {
1029:     PetscReal inv1pX = 1. / (1. + x);

1031:     R[0] = x;
1032:     R[1] = -y;
1033:     R[2] = -z;
1034:     R[3] = y;
1035:     R[4] = 1. - y * y * inv1pX;
1036:     R[5] = -y * z * inv1pX;
1037:     R[6] = z;
1038:     R[7] = -y * z * inv1pX;
1039:     R[8] = 1. - z * z * inv1pX;
1040:   } else {
1041:     PetscReal inv1mX = 1. / (1. - x);

1043:     R[0] = x;
1044:     R[1] = z;
1045:     R[2] = y;
1046:     R[3] = y;
1047:     R[4] = -y * z * inv1mX;
1048:     R[5] = 1. - y * y * inv1mX;
1049:     R[6] = z;
1050:     R[7] = 1. - z * z * inv1mX;
1051:     R[8] = -y * z * inv1mX;
1052:   }
1053:   coords[0] = 0.0;
1054:   coords[1] = r;
1055:   return 0;
1056: }

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

1062:   Not collective

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

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

1071:   Output Parameter:
1072: . 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.

1074:   Level: developer

1076: .seealso: `DMPlexComputeProjection2Dto1D()`, `DMPlexComputeProjection3Dto1D()`
1077: @*/
1078: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
1079: {
1080:   PetscReal      x1[3], x2[3], n[3], c[3], norm;
1081:   const PetscInt dim = 3;
1082:   PetscInt       d, p;

1084:   /* 0) Calculate normal vector */
1085:   for (d = 0; d < dim; ++d) {
1086:     x1[d] = PetscRealPart(coords[1 * dim + d] - coords[0 * dim + d]);
1087:     x2[d] = PetscRealPart(coords[2 * dim + d] - coords[0 * dim + d]);
1088:   }
1089:   // n = x1 \otimes x2
1090:   n[0] = x1[1] * x2[2] - x1[2] * x2[1];
1091:   n[1] = x1[2] * x2[0] - x1[0] * x2[2];
1092:   n[2] = x1[0] * x2[1] - x1[1] * x2[0];
1093:   norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
1094:   for (d = 0; d < dim; d++) n[d] /= norm;
1095:   norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
1096:   for (d = 0; d < dim; d++) x1[d] /= norm;
1097:   // x2 = n \otimes x1
1098:   x2[0] = n[1] * x1[2] - n[2] * x1[1];
1099:   x2[1] = n[2] * x1[0] - n[0] * x1[2];
1100:   x2[2] = n[0] * x1[1] - n[1] * x1[0];
1101:   for (d = 0; d < dim; d++) {
1102:     R[d * dim + 0] = x1[d];
1103:     R[d * dim + 1] = x2[d];
1104:     R[d * dim + 2] = n[d];
1105:     c[d]           = PetscRealPart(coords[0 * dim + d]);
1106:   }
1107:   for (p = 0; p < coordSize / dim; p++) {
1108:     PetscReal y[3];
1109:     for (d = 0; d < dim; d++) y[d] = PetscRealPart(coords[p * dim + d]) - c[d];
1110:     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];
1111:   }
1112:   return 0;
1113: }

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

1119:    |  1  1  1 |
1120:    | x0 x1 x2 |
1121:    | y0 y1 y2 |

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

1125:    | x1 x2 |
1126:    | y1 y2 |
1127:   */
1128:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1129:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1130:   PetscReal       M[4], detM;
1131:   M[0] = x1;
1132:   M[1] = x2;
1133:   M[2] = y1;
1134:   M[3] = y2;
1135:   DMPlex_Det2D_Internal(&detM, M);
1136:   *vol = 0.5 * detM;
1137:   (void)PetscLogFlops(5.0);
1138: }

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

1144:    |  1  1  1  1 |
1145:    | x0 x1 x2 x3 |
1146:    | y0 y1 y2 y3 |
1147:    | z0 z1 z2 z3 |

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

1151:    | x1 x2 x3 |
1152:    | y1 y2 y3 |
1153:    | z1 z2 z3 |
1154:   */
1155:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4] - coords[1], z1 = coords[5] - coords[2];
1156:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7] - coords[1], z2 = coords[8] - coords[2];
1157:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1158:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1159:   PetscReal       M[9], detM;
1160:   M[0] = x1;
1161:   M[1] = x2;
1162:   M[2] = x3;
1163:   M[3] = y1;
1164:   M[4] = y2;
1165:   M[5] = y3;
1166:   M[6] = z1;
1167:   M[7] = z2;
1168:   M[8] = z3;
1169:   DMPlex_Det3D_Internal(&detM, M);
1170:   *vol = -onesixth * detM;
1171:   (void)PetscLogFlops(10.0);
1172: }

1174: static inline void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1175: {
1176:   const PetscReal onesixth = ((PetscReal)1. / (PetscReal)6.);
1177:   DMPlex_Det3D_Internal(vol, coords);
1178:   *vol *= -onesixth;
1179: }

1181: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1182: {
1183:   PetscSection       coordSection;
1184:   Vec                coordinates;
1185:   const PetscScalar *coords;
1186:   PetscInt           dim, d, off;

1188:   DMGetCoordinatesLocal(dm, &coordinates);
1189:   DMGetCoordinateSection(dm, &coordSection);
1190:   PetscSectionGetDof(coordSection, e, &dim);
1191:   if (!dim) return 0;
1192:   PetscSectionGetOffset(coordSection, e, &off);
1193:   VecGetArrayRead(coordinates, &coords);
1194:   if (v0) {
1195:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);
1196:   }
1197:   VecRestoreArrayRead(coordinates, &coords);
1198:   *detJ = 1.;
1199:   if (J) {
1200:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1201:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1202:     if (invJ) {
1203:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1204:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1205:     }
1206:   }
1207:   return 0;
1208: }

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

1213:   Not collective

1215:   Input Parameters:
1216: + dm   - The DM
1217: - cell - The cell number

1219:   Output Parameters:
1220: + isDG   - Using cellwise coordinates
1221: . Nc     - The number of coordinates
1222: . array  - The coordinate array
1223: - coords - The cell coordinates

1225:   Level: developer

1227: .seealso: DMPlexRestoreCellCoordinates(), DMGetCoordinatesLocal(), DMGetCellCoordinatesLocal()
1228: @*/
1229: PetscErrorCode DMPlexGetCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1230: {
1231:   DM                 cdm;
1232:   Vec                coordinates;
1233:   PetscSection       cs;
1234:   const PetscScalar *ccoords;
1235:   PetscInt           pStart, pEnd;

1238:   *isDG   = PETSC_FALSE;
1239:   *Nc     = 0;
1240:   *array  = NULL;
1241:   *coords = NULL;
1242:   /* Check for cellwise coordinates */
1243:   DMGetCellCoordinateSection(dm, &cs);
1244:   if (!cs) goto cg;
1245:   /* Check that the cell exists in the cellwise section */
1246:   PetscSectionGetChart(cs, &pStart, &pEnd);
1247:   if (cell < pStart || cell >= pEnd) goto cg;
1248:   /* Check for cellwise coordinates for this cell */
1249:   PetscSectionGetDof(cs, cell, Nc);
1250:   if (!*Nc) goto cg;
1251:   /* Check for cellwise coordinates */
1252:   DMGetCellCoordinatesLocalNoncollective(dm, &coordinates);
1253:   if (!coordinates) goto cg;
1254:   /* Get cellwise coordinates */
1255:   DMGetCellCoordinateDM(dm, &cdm);
1256:   VecGetArrayRead(coordinates, array);
1257:   DMPlexPointLocalRead(cdm, cell, *array, &ccoords);
1258:   DMGetWorkArray(cdm, *Nc, MPIU_SCALAR, coords);
1259:   PetscArraycpy(*coords, ccoords, *Nc);
1260:   VecRestoreArrayRead(coordinates, array);
1261:   *isDG = PETSC_TRUE;
1262:   return 0;
1263: cg:
1264:   /* Use continuous coordinates */
1265:   DMGetCoordinateDM(dm, &cdm);
1266:   DMGetCoordinateSection(dm, &cs);
1267:   DMGetCoordinatesLocalNoncollective(dm, &coordinates);
1268:   DMPlexVecGetClosure(cdm, cs, coordinates, cell, Nc, coords);
1269:   return 0;
1270: }

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

1275:   Not collective

1277:   Input Parameters:
1278: + dm   - The DM
1279: - cell - The cell number

1281:   Output Parameters:
1282: + isDG   - Using cellwise coordinates
1283: . Nc     - The number of coordinates
1284: . array  - The coordinate array
1285: - coords - The cell coordinates

1287:   Level: developer

1289: .seealso: DMPlexGetCellCoordinates(), DMGetCoordinatesLocal(), DMGetCellCoordinatesLocal()
1290: @*/
1291: PetscErrorCode DMPlexRestoreCellCoordinates(DM dm, PetscInt cell, PetscBool *isDG, PetscInt *Nc, const PetscScalar *array[], PetscScalar *coords[])
1292: {
1293:   DM           cdm;
1294:   PetscSection cs;
1295:   Vec          coordinates;

1298:   if (*isDG) {
1299:     DMGetCellCoordinateDM(dm, &cdm);
1300:     DMRestoreWorkArray(cdm, *Nc, MPIU_SCALAR, coords);
1301:   } else {
1302:     DMGetCoordinateDM(dm, &cdm);
1303:     DMGetCoordinateSection(dm, &cs);
1304:     DMGetCoordinatesLocalNoncollective(dm, &coordinates);
1305:     DMPlexVecRestoreClosure(cdm, cs, coordinates, cell, Nc, (PetscScalar **)coords);
1306:   }
1307:   return 0;
1308: }

1310: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1311: {
1312:   const PetscScalar *array;
1313:   PetscScalar       *coords = NULL;
1314:   PetscInt           numCoords, d;
1315:   PetscBool          isDG;

1317:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1319:   *detJ = 0.0;
1320:   if (numCoords == 6) {
1321:     const PetscInt dim = 3;
1322:     PetscReal      R[9], J0;

1324:     if (v0) {
1325:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1326:     }
1327:     DMPlexComputeProjection3Dto1D(coords, R);
1328:     if (J) {
1329:       J0   = 0.5 * PetscRealPart(coords[1]);
1330:       J[0] = R[0] * J0;
1331:       J[1] = R[1];
1332:       J[2] = R[2];
1333:       J[3] = R[3] * J0;
1334:       J[4] = R[4];
1335:       J[5] = R[5];
1336:       J[6] = R[6] * J0;
1337:       J[7] = R[7];
1338:       J[8] = R[8];
1339:       DMPlex_Det3D_Internal(detJ, J);
1340:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1341:     }
1342:   } else if (numCoords == 4) {
1343:     const PetscInt dim = 2;
1344:     PetscReal      R[4], J0;

1346:     if (v0) {
1347:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1348:     }
1349:     DMPlexComputeProjection2Dto1D(coords, R);
1350:     if (J) {
1351:       J0   = 0.5 * PetscRealPart(coords[1]);
1352:       J[0] = R[0] * J0;
1353:       J[1] = R[1];
1354:       J[2] = R[2] * J0;
1355:       J[3] = R[3];
1356:       DMPlex_Det2D_Internal(detJ, J);
1357:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1358:     }
1359:   } else if (numCoords == 2) {
1360:     const PetscInt dim = 1;

1362:     if (v0) {
1363:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1364:     }
1365:     if (J) {
1366:       J[0]  = 0.5 * (PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1367:       *detJ = J[0];
1368:       PetscLogFlops(2.0);
1369:       if (invJ) {
1370:         invJ[0] = 1.0 / J[0];
1371:         PetscLogFlops(1.0);
1372:       }
1373:     }
1374:   } 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);
1375:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1376:   return 0;
1377: }

1379: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1380: {
1381:   const PetscScalar *array;
1382:   PetscScalar       *coords = NULL;
1383:   PetscInt           numCoords, d;
1384:   PetscBool          isDG;

1386:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1388:   *detJ = 0.0;
1389:   if (numCoords == 9) {
1390:     const PetscInt dim = 3;
1391:     PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1393:     if (v0) {
1394:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1395:     }
1396:     DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1397:     if (J) {
1398:       const PetscInt pdim = 2;

1400:       for (d = 0; d < pdim; d++) {
1401:         for (PetscInt f = 0; f < pdim; f++) J0[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * pdim + d]) - PetscRealPart(coords[0 * pdim + d]));
1402:       }
1403:       PetscLogFlops(8.0);
1404:       DMPlex_Det3D_Internal(detJ, J0);
1405:       for (d = 0; d < dim; d++) {
1406:         for (PetscInt f = 0; f < dim; f++) {
1407:           J[d * dim + f] = 0.0;
1408:           for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1409:         }
1410:       }
1411:       PetscLogFlops(18.0);
1412:     }
1413:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1414:   } else if (numCoords == 6) {
1415:     const PetscInt dim = 2;

1417:     if (v0) {
1418:       for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1419:     }
1420:     if (J) {
1421:       for (d = 0; d < dim; d++) {
1422:         for (PetscInt f = 0; f < dim; f++) J[d * dim + f] = 0.5 * (PetscRealPart(coords[(f + 1) * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1423:       }
1424:       PetscLogFlops(8.0);
1425:       DMPlex_Det2D_Internal(detJ, J);
1426:     }
1427:     if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1428:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %" PetscInt_FMT " != 6 or 9", numCoords);
1429:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1430:   return 0;
1431: }

1433: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1434: {
1435:   const PetscScalar *array;
1436:   PetscScalar       *coords = NULL;
1437:   PetscInt           numCoords, d;
1438:   PetscBool          isDG;

1440:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1442:   if (!Nq) {
1443:     PetscInt vorder[4] = {0, 1, 2, 3};

1445:     if (isTensor) {
1446:       vorder[2] = 3;
1447:       vorder[3] = 2;
1448:     }
1449:     *detJ = 0.0;
1450:     if (numCoords == 12) {
1451:       const PetscInt dim = 3;
1452:       PetscReal      R[9], J0[9] = {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0};

1454:       if (v) {
1455:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1456:       }
1457:       DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1458:       if (J) {
1459:         const PetscInt pdim = 2;

1461:         for (d = 0; d < pdim; d++) {
1462:           J0[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * pdim + d]) - PetscRealPart(coords[vorder[0] * pdim + d]));
1463:           J0[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[2] * pdim + d]) - PetscRealPart(coords[vorder[1] * pdim + d]));
1464:         }
1465:         PetscLogFlops(8.0);
1466:         DMPlex_Det3D_Internal(detJ, J0);
1467:         for (d = 0; d < dim; d++) {
1468:           for (PetscInt f = 0; f < dim; f++) {
1469:             J[d * dim + f] = 0.0;
1470:             for (PetscInt g = 0; g < dim; g++) J[d * dim + f] += R[d * dim + g] * J0[g * dim + f];
1471:           }
1472:         }
1473:         PetscLogFlops(18.0);
1474:       }
1475:       if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1476:     } else if (numCoords == 8) {
1477:       const PetscInt dim = 2;

1479:       if (v) {
1480:         for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1481:       }
1482:       if (J) {
1483:         for (d = 0; d < dim; d++) {
1484:           J[d * dim + 0] = 0.5 * (PetscRealPart(coords[vorder[1] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1485:           J[d * dim + 1] = 0.5 * (PetscRealPart(coords[vorder[3] * dim + d]) - PetscRealPart(coords[vorder[0] * dim + d]));
1486:         }
1487:         PetscLogFlops(8.0);
1488:         DMPlex_Det2D_Internal(detJ, J);
1489:       }
1490:       if (invJ) DMPlex_Invert2D_Internal(invJ, J, *detJ);
1491:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1492:   } else {
1493:     const PetscInt Nv         = 4;
1494:     const PetscInt dimR       = 2;
1495:     PetscInt       zToPlex[4] = {0, 1, 3, 2};
1496:     PetscReal      zOrder[12];
1497:     PetscReal      zCoeff[12];
1498:     PetscInt       i, j, k, l, dim;

1500:     if (isTensor) {
1501:       zToPlex[2] = 2;
1502:       zToPlex[3] = 3;
1503:     }
1504:     if (numCoords == 12) {
1505:       dim = 3;
1506:     } else if (numCoords == 8) {
1507:       dim = 2;
1508:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %" PetscInt_FMT " != 8 or 12", numCoords);
1509:     for (i = 0; i < Nv; i++) {
1510:       PetscInt zi = zToPlex[i];

1512:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1513:     }
1514:     for (j = 0; j < dim; j++) {
1515:       /* Nodal basis for evaluation at the vertices: (1 \mp xi) (1 \mp eta):
1516:            \phi^0 = (1 - xi - eta + xi eta) --> 1      = 1/4 ( \phi^0 + \phi^1 + \phi^2 + \phi^3)
1517:            \phi^1 = (1 + xi - eta - xi eta) --> xi     = 1/4 (-\phi^0 + \phi^1 - \phi^2 + \phi^3)
1518:            \phi^2 = (1 - xi + eta - xi eta) --> eta    = 1/4 (-\phi^0 - \phi^1 + \phi^2 + \phi^3)
1519:            \phi^3 = (1 + xi + eta + xi eta) --> xi eta = 1/4 ( \phi^0 - \phi^1 - \phi^2 + \phi^3)
1520:       */
1521:       zCoeff[dim * 0 + j] = 0.25 * (zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1522:       zCoeff[dim * 1 + j] = 0.25 * (-zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1523:       zCoeff[dim * 2 + j] = 0.25 * (-zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1524:       zCoeff[dim * 3 + j] = 0.25 * (zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1525:     }
1526:     for (i = 0; i < Nq; i++) {
1527:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1529:       if (v) {
1530:         PetscReal extPoint[4];

1532:         extPoint[0] = 1.;
1533:         extPoint[1] = xi;
1534:         extPoint[2] = eta;
1535:         extPoint[3] = xi * eta;
1536:         for (j = 0; j < dim; j++) {
1537:           PetscReal val = 0.;

1539:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1540:           v[i * dim + j] = val;
1541:         }
1542:       }
1543:       if (J) {
1544:         PetscReal extJ[8];

1546:         extJ[0] = 0.;
1547:         extJ[1] = 0.;
1548:         extJ[2] = 1.;
1549:         extJ[3] = 0.;
1550:         extJ[4] = 0.;
1551:         extJ[5] = 1.;
1552:         extJ[6] = eta;
1553:         extJ[7] = xi;
1554:         for (j = 0; j < dim; j++) {
1555:           for (k = 0; k < dimR; k++) {
1556:             PetscReal val = 0.;

1558:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1559:             J[i * dim * dim + dim * j + k] = val;
1560:           }
1561:         }
1562:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1563:           PetscReal  x, y, z;
1564:           PetscReal *iJ = &J[i * dim * dim];
1565:           PetscReal  norm;

1567:           x     = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1568:           y     = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1569:           z     = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1570:           norm  = PetscSqrtReal(x * x + y * y + z * z);
1571:           iJ[2] = x / norm;
1572:           iJ[5] = y / norm;
1573:           iJ[8] = z / norm;
1574:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1575:           if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1576:         } else {
1577:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1578:           if (invJ) DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1579:         }
1580:       }
1581:     }
1582:   }
1583:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1584:   return 0;
1585: }

1587: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1588: {
1589:   const PetscScalar *array;
1590:   PetscScalar       *coords = NULL;
1591:   const PetscInt     dim    = 3;
1592:   PetscInt           numCoords, d;
1593:   PetscBool          isDG;

1595:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1597:   *detJ = 0.0;
1598:   if (v0) {
1599:     for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);
1600:   }
1601:   if (J) {
1602:     for (d = 0; d < dim; d++) {
1603:       /* I orient with outward face normals */
1604:       J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1605:       J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1606:       J[d * dim + 2] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1607:     }
1608:     PetscLogFlops(18.0);
1609:     DMPlex_Det3D_Internal(detJ, J);
1610:   }
1611:   if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1612:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1613:   return 0;
1614: }

1616: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1617: {
1618:   const PetscScalar *array;
1619:   PetscScalar       *coords = NULL;
1620:   const PetscInt     dim    = 3;
1621:   PetscInt           numCoords, d;
1622:   PetscBool          isDG;

1624:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1626:   if (!Nq) {
1627:     *detJ = 0.0;
1628:     if (v) {
1629:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1630:     }
1631:     if (J) {
1632:       for (d = 0; d < dim; d++) {
1633:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[3 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1634:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1635:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1636:       }
1637:       PetscLogFlops(18.0);
1638:       DMPlex_Det3D_Internal(detJ, J);
1639:     }
1640:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1641:   } else {
1642:     const PetscInt Nv         = 8;
1643:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1644:     const PetscInt dim        = 3;
1645:     const PetscInt dimR       = 3;
1646:     PetscReal      zOrder[24];
1647:     PetscReal      zCoeff[24];
1648:     PetscInt       i, j, k, l;

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

1653:       for (j = 0; j < dim; j++) zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1654:     }
1655:     for (j = 0; j < dim; j++) {
1656:       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]);
1657:       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]);
1658:       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]);
1659:       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]);
1660:       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]);
1661:       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]);
1662:       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]);
1663:       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]);
1664:     }
1665:     for (i = 0; i < Nq; i++) {
1666:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

1668:       if (v) {
1669:         PetscReal extPoint[8];

1671:         extPoint[0] = 1.;
1672:         extPoint[1] = xi;
1673:         extPoint[2] = eta;
1674:         extPoint[3] = xi * eta;
1675:         extPoint[4] = theta;
1676:         extPoint[5] = theta * xi;
1677:         extPoint[6] = theta * eta;
1678:         extPoint[7] = theta * eta * xi;
1679:         for (j = 0; j < dim; j++) {
1680:           PetscReal val = 0.;

1682:           for (k = 0; k < Nv; k++) val += extPoint[k] * zCoeff[dim * k + j];
1683:           v[i * dim + j] = val;
1684:         }
1685:       }
1686:       if (J) {
1687:         PetscReal extJ[24];

1689:         extJ[0]  = 0.;
1690:         extJ[1]  = 0.;
1691:         extJ[2]  = 0.;
1692:         extJ[3]  = 1.;
1693:         extJ[4]  = 0.;
1694:         extJ[5]  = 0.;
1695:         extJ[6]  = 0.;
1696:         extJ[7]  = 1.;
1697:         extJ[8]  = 0.;
1698:         extJ[9]  = eta;
1699:         extJ[10] = xi;
1700:         extJ[11] = 0.;
1701:         extJ[12] = 0.;
1702:         extJ[13] = 0.;
1703:         extJ[14] = 1.;
1704:         extJ[15] = theta;
1705:         extJ[16] = 0.;
1706:         extJ[17] = xi;
1707:         extJ[18] = 0.;
1708:         extJ[19] = theta;
1709:         extJ[20] = eta;
1710:         extJ[21] = theta * eta;
1711:         extJ[22] = theta * xi;
1712:         extJ[23] = eta * xi;

1714:         for (j = 0; j < dim; j++) {
1715:           for (k = 0; k < dimR; k++) {
1716:             PetscReal val = 0.;

1718:             for (l = 0; l < Nv; l++) val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1719:             J[i * dim * dim + dim * j + k] = val;
1720:           }
1721:         }
1722:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1723:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1724:       }
1725:     }
1726:   }
1727:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1728:   return 0;
1729: }

1731: static PetscErrorCode DMPlexComputeTriangularPrismGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1732: {
1733:   const PetscScalar *array;
1734:   PetscScalar       *coords = NULL;
1735:   const PetscInt     dim    = 3;
1736:   PetscInt           numCoords, d;
1737:   PetscBool          isDG;

1739:   DMPlexGetCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1741:   if (!Nq) {
1742:     /* Assume that the map to the reference is affine */
1743:     *detJ = 0.0;
1744:     if (v) {
1745:       for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);
1746:     }
1747:     if (J) {
1748:       for (d = 0; d < dim; d++) {
1749:         J[d * dim + 0] = 0.5 * (PetscRealPart(coords[2 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1750:         J[d * dim + 1] = 0.5 * (PetscRealPart(coords[1 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1751:         J[d * dim + 2] = 0.5 * (PetscRealPart(coords[4 * dim + d]) - PetscRealPart(coords[0 * dim + d]));
1752:       }
1753:       PetscLogFlops(18.0);
1754:       DMPlex_Det3D_Internal(detJ, J);
1755:     }
1756:     if (invJ) DMPlex_Invert3D_Internal(invJ, J, *detJ);
1757:   } else {
1758:     const PetscInt dim  = 3;
1759:     const PetscInt dimR = 3;
1760:     const PetscInt Nv   = 6;
1761:     PetscReal      verts[18];
1762:     PetscReal      coeff[18];
1763:     PetscInt       i, j, k, l;

1765:     for (i = 0; i < Nv; ++i)
1766:       for (j = 0; j < dim; ++j) verts[dim * i + j] = PetscRealPart(coords[dim * i + j]);
1767:     for (j = 0; j < dim; ++j) {
1768:       /* Check for triangle,
1769:            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)
1770:            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)
1771:            phi^2 =  1/2 (1 + eta)   chi^2 = delta(-1,  1)

1773:            phi^0 + phi^1 + phi^2 = 1    coef_1   = 1/2 (         chi^1 + chi^2)
1774:           -phi^0 + phi^1 - phi^2 = xi   coef_xi  = 1/2 (-chi^0 + chi^1)
1775:           -phi^0 - phi^1 + phi^2 = eta  coef_eta = 1/2 (-chi^0         + chi^2)

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

1781:           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
1782:       */
1783:       /* Nodal basis for evaluation at the vertices: {-xi - eta, 1 + xi, 1 + eta} (1 \mp zeta):
1784:            \phi^0 = 1/4 (   -xi - eta        + xi zeta + eta zeta) --> /  1  1  1  1  1  1 \ 1
1785:            \phi^1 = 1/4 (1      + eta - zeta           - eta zeta) --> | -1  1 -1 -1 -1  1 | eta
1786:            \phi^2 = 1/4 (1 + xi       - zeta - xi zeta)            --> | -1 -1  1 -1  1 -1 | xi
1787:            \phi^3 = 1/4 (   -xi - eta        - xi zeta - eta zeta) --> | -1 -1 -1  1  1  1 | zeta
1788:            \phi^4 = 1/4 (1 + xi       + zeta + xi zeta)            --> |  1  1 -1 -1  1 -1 | xi zeta
1789:            \phi^5 = 1/4 (1      + eta + zeta           + eta zeta) --> \  1 -1  1 -1 -1  1 / eta zeta
1790:            1/4 /  0  1  1  0  1  1 \
1791:                | -1  1  0 -1  0  1 |
1792:                | -1  0  1 -1  1  0 |
1793:                |  0 -1 -1  0  1  1 |
1794:                |  1  0 -1 -1  1  0 |
1795:                \  1 -1  0 -1  0  1 /
1796:       */
1797:       coeff[dim * 0 + j] = (1. / 4.) * (verts[dim * 1 + j] + verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1798:       coeff[dim * 1 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1799:       coeff[dim * 2 + j] = (1. / 4.) * (-verts[dim * 0 + j] + verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1800:       coeff[dim * 3 + j] = (1. / 4.) * (-verts[dim * 1 + j] - verts[dim * 2 + j] + verts[dim * 4 + j] + verts[dim * 5 + j]);
1801:       coeff[dim * 4 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 2 + j] - verts[dim * 3 + j] + verts[dim * 4 + j]);
1802:       coeff[dim * 5 + j] = (1. / 4.) * (verts[dim * 0 + j] - verts[dim * 1 + j] - verts[dim * 3 + j] + verts[dim * 5 + j]);
1803:       /* For reference prism:
1804:       {0, 0, 0}
1805:       {0, 1, 0}
1806:       {1, 0, 0}
1807:       {0, 0, 1}
1808:       {0, 0, 0}
1809:       {0, 0, 0}
1810:       */
1811:     }
1812:     for (i = 0; i < Nq; ++i) {
1813:       const PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], zeta = points[dimR * i + 2];

1815:       if (v) {
1816:         PetscReal extPoint[6];
1817:         PetscInt  c;

1819:         extPoint[0] = 1.;
1820:         extPoint[1] = eta;
1821:         extPoint[2] = xi;
1822:         extPoint[3] = zeta;
1823:         extPoint[4] = xi * zeta;
1824:         extPoint[5] = eta * zeta;
1825:         for (c = 0; c < dim; ++c) {
1826:           PetscReal val = 0.;

1828:           for (k = 0; k < Nv; ++k) val += extPoint[k] * coeff[k * dim + c];
1829:           v[i * dim + c] = val;
1830:         }
1831:       }
1832:       if (J) {
1833:         PetscReal extJ[18];

1835:         extJ[0]  = 0.;
1836:         extJ[1]  = 0.;
1837:         extJ[2]  = 0.;
1838:         extJ[3]  = 0.;
1839:         extJ[4]  = 1.;
1840:         extJ[5]  = 0.;
1841:         extJ[6]  = 1.;
1842:         extJ[7]  = 0.;
1843:         extJ[8]  = 0.;
1844:         extJ[9]  = 0.;
1845:         extJ[10] = 0.;
1846:         extJ[11] = 1.;
1847:         extJ[12] = zeta;
1848:         extJ[13] = 0.;
1849:         extJ[14] = xi;
1850:         extJ[15] = 0.;
1851:         extJ[16] = zeta;
1852:         extJ[17] = eta;

1854:         for (j = 0; j < dim; j++) {
1855:           for (k = 0; k < dimR; k++) {
1856:             PetscReal val = 0.;

1858:             for (l = 0; l < Nv; l++) val += coeff[dim * l + j] * extJ[dimR * l + k];
1859:             J[i * dim * dim + dim * j + k] = val;
1860:           }
1861:         }
1862:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1863:         if (invJ) DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);
1864:       }
1865:     }
1866:   }
1867:   DMPlexRestoreCellCoordinates(dm, e, &isDG, &numCoords, &array, &coords);
1868:   return 0;
1869: }

1871: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1872: {
1873:   DMPolytopeType   ct;
1874:   PetscInt         depth, dim, coordDim, coneSize, i;
1875:   PetscInt         Nq     = 0;
1876:   const PetscReal *points = NULL;
1877:   DMLabel          depthLabel;
1878:   PetscReal        xi0[3]   = {-1., -1., -1.}, v0[3], J0[9], detJ0;
1879:   PetscBool        isAffine = PETSC_TRUE;

1881:   DMPlexGetDepth(dm, &depth);
1882:   DMPlexGetConeSize(dm, cell, &coneSize);
1883:   DMPlexGetDepthLabel(dm, &depthLabel);
1884:   DMLabelGetValue(depthLabel, cell, &dim);
1885:   if (depth == 1 && dim == 1) DMGetDimension(dm, &dim);
1886:   DMGetCoordinateDim(dm, &coordDim);
1888:   if (quad) PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);
1889:   DMPlexGetCellType(dm, cell, &ct);
1890:   switch (ct) {
1891:   case DM_POLYTOPE_POINT:
1892:     DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1893:     isAffine = PETSC_FALSE;
1894:     break;
1895:   case DM_POLYTOPE_SEGMENT:
1896:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
1897:     if (Nq) DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1898:     else DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);
1899:     break;
1900:   case DM_POLYTOPE_TRIANGLE:
1901:     if (Nq) DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1902:     else DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);
1903:     break;
1904:   case DM_POLYTOPE_QUADRILATERAL:
1905:     DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);
1906:     isAffine = PETSC_FALSE;
1907:     break;
1908:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
1909:     DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);
1910:     isAffine = PETSC_FALSE;
1911:     break;
1912:   case DM_POLYTOPE_TETRAHEDRON:
1913:     if (Nq) DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1914:     else DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);
1915:     break;
1916:   case DM_POLYTOPE_HEXAHEDRON:
1917:     DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1918:     isAffine = PETSC_FALSE;
1919:     break;
1920:   case DM_POLYTOPE_TRI_PRISM:
1921:     DMPlexComputeTriangularPrismGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1922:     isAffine = PETSC_FALSE;
1923:     break;
1924:   default:
1925:     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))]);
1926:   }
1927:   if (isAffine && Nq) {
1928:     if (v) {
1929:       for (i = 0; i < Nq; i++) CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1930:     }
1931:     if (detJ) {
1932:       for (i = 0; i < Nq; i++) detJ[i] = detJ0;
1933:     }
1934:     if (J) {
1935:       PetscInt k;

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

1940:         for (j = 0; j < coordDim * coordDim; j++, k++) J[k] = J0[j];
1941:       }
1942:     }
1943:     if (invJ) {
1944:       PetscInt k;
1945:       switch (coordDim) {
1946:       case 0:
1947:         break;
1948:       case 1:
1949:         invJ[0] = 1. / J0[0];
1950:         break;
1951:       case 2:
1952:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1953:         break;
1954:       case 3:
1955:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1956:         break;
1957:       }
1958:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1959:         PetscInt j;

1961:         for (j = 0; j < coordDim * coordDim; j++, k++) invJ[k] = invJ[j];
1962:       }
1963:     }
1964:   }
1965:   return 0;
1966: }

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

1971:   Collective on dm

1973:   Input Parameters:
1974: + dm   - the DM
1975: - cell - the cell

1977:   Output Parameters:
1978: + v0   - the translation part of this affine transform, meaning the translation to the origin (not the first vertex of the reference cell)
1979: . J    - the Jacobian of the transform from the reference element
1980: . invJ - the inverse of the Jacobian
1981: - detJ - the Jacobian determinant

1983:   Level: advanced

1985:   Fortran Notes:
1986:   Since it returns arrays, this routine is only available in Fortran 90, and you must
1987:   include petsc.h90 in your code.

1989: .seealso: `DMPlexComputeCellGeometryFEM()`, `DMGetCoordinateSection()`, `DMGetCoordinates()`
1990: @*/
1991: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1992: {
1993:   DMPlexComputeCellGeometryFEM_Implicit(dm, cell, NULL, v0, J, invJ, detJ);
1994:   return 0;
1995: }

1997: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1998: {
1999:   const PetscScalar *array;
2000:   PetscScalar       *coords = NULL;
2001:   PetscInt           numCoords;
2002:   PetscBool          isDG;
2003:   PetscQuadrature    feQuad;
2004:   const PetscReal   *quadPoints;
2005:   PetscTabulation    T;
2006:   PetscInt           dim, cdim, pdim, qdim, Nq, q;

2008:   DMGetDimension(dm, &dim);
2009:   DMGetCoordinateDim(dm, &cdim);
2010:   DMPlexGetCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords);
2011:   if (!quad) { /* use the first point of the first functional of the dual space */
2012:     PetscDualSpace dsp;

2014:     PetscFEGetDualSpace(fe, &dsp);
2015:     PetscDualSpaceGetFunctional(dsp, 0, &quad);
2016:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
2017:     Nq = 1;
2018:   } else {
2019:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
2020:   }
2021:   PetscFEGetDimension(fe, &pdim);
2022:   PetscFEGetQuadrature(fe, &feQuad);
2023:   if (feQuad == quad) {
2024:     PetscFEGetCellTabulation(fe, J ? 1 : 0, &T);
2026:   } else {
2027:     PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
2028:   }
2030:   {
2031:     const PetscReal *basis    = T->T[0];
2032:     const PetscReal *basisDer = T->T[1];
2033:     PetscReal        detJt;

2035: #if defined(PETSC_USE_DEBUG)
2040: #endif
2041:     if (v) {
2042:       PetscArrayzero(v, Nq * cdim);
2043:       for (q = 0; q < Nq; ++q) {
2044:         PetscInt i, k;

2046:         for (k = 0; k < pdim; ++k) {
2047:           const PetscInt vertex = k / cdim;
2048:           for (i = 0; i < cdim; ++i) v[q * cdim + i] += basis[(q * pdim + k) * cdim + i] * PetscRealPart(coords[vertex * cdim + i]);
2049:         }
2050:         PetscLogFlops(2.0 * pdim * cdim);
2051:       }
2052:     }
2053:     if (J) {
2054:       PetscArrayzero(J, Nq * cdim * cdim);
2055:       for (q = 0; q < Nq; ++q) {
2056:         PetscInt i, j, k, c, r;

2058:         /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
2059:         for (k = 0; k < pdim; ++k) {
2060:           const PetscInt vertex = k / cdim;
2061:           for (j = 0; j < dim; ++j) {
2062:             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]);
2063:           }
2064:         }
2065:         PetscLogFlops(2.0 * pdim * dim * cdim);
2066:         if (cdim > dim) {
2067:           for (c = dim; c < cdim; ++c)
2068:             for (r = 0; r < cdim; ++r) J[r * cdim + c] = r == c ? 1.0 : 0.0;
2069:         }
2070:         if (!detJ && !invJ) continue;
2071:         detJt = 0.;
2072:         switch (cdim) {
2073:         case 3:
2074:           DMPlex_Det3D_Internal(&detJt, &J[q * cdim * dim]);
2075:           if (invJ) DMPlex_Invert3D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2076:           break;
2077:         case 2:
2078:           DMPlex_Det2D_Internal(&detJt, &J[q * cdim * dim]);
2079:           if (invJ) DMPlex_Invert2D_Internal(&invJ[q * cdim * dim], &J[q * cdim * dim], detJt);
2080:           break;
2081:         case 1:
2082:           detJt = J[q * cdim * dim];
2083:           if (invJ) invJ[q * cdim * dim] = 1.0 / detJt;
2084:         }
2085:         if (detJ) detJ[q] = detJt;
2086:       }
2088:   }
2089:   if (feQuad != quad) PetscTabulationDestroy(&T);
2090:   DMPlexRestoreCellCoordinates(dm, point, &isDG, &numCoords, &array, &coords);
2091:   return 0;
2092: }

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

2097:   Collective on dm

2099:   Input Parameters:
2100: + dm   - the DM
2101: . cell - the cell
2102: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If quad == NULL, geometry will be
2103:          evaluated at the first vertex of the reference element

2105:   Output Parameters:
2106: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
2107: . J    - the Jacobian of the transform from the reference element at each quadrature point
2108: . invJ - the inverse of the Jacobian at each quadrature point
2109: - detJ - the Jacobian determinant at each quadrature point

2111:   Level: advanced

2113:   Fortran Notes:
2114:   Since it returns arrays, this routine is only available in Fortran 90, and you must
2115:   include petsc.h90 in your code.

2117: .seealso: `DMGetCoordinateSection()`, `DMGetCoordinates()`
2118: @*/
2119: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
2120: {
2121:   DM      cdm;
2122:   PetscFE fe = NULL;

2125:   DMGetCoordinateDM(dm, &cdm);
2126:   if (cdm) {
2127:     PetscClassId id;
2128:     PetscInt     numFields;
2129:     PetscDS      prob;
2130:     PetscObject  disc;

2132:     DMGetNumFields(cdm, &numFields);
2133:     if (numFields) {
2134:       DMGetDS(cdm, &prob);
2135:       PetscDSGetDiscretization(prob, 0, &disc);
2136:       PetscObjectGetClassId(disc, &id);
2137:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
2138:     }
2139:   }
2140:   if (!fe) DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);
2141:   else DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);
2142:   return 0;
2143: }

2145: static PetscErrorCode DMPlexComputeGeometryFVM_0D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2146: {
2147:   PetscSection       coordSection;
2148:   Vec                coordinates;
2149:   const PetscScalar *coords = NULL;
2150:   PetscInt           d, dof, off;

2152:   DMGetCoordinatesLocal(dm, &coordinates);
2153:   DMGetCoordinateSection(dm, &coordSection);
2154:   VecGetArrayRead(coordinates, &coords);

2156:   /* for a point the centroid is just the coord */
2157:   if (centroid) {
2158:     PetscSectionGetDof(coordSection, cell, &dof);
2159:     PetscSectionGetOffset(coordSection, cell, &off);
2160:     for (d = 0; d < dof; d++) centroid[d] = PetscRealPart(coords[off + d]);
2161:   }
2162:   if (normal) {
2163:     const PetscInt *support, *cones;
2164:     PetscInt        supportSize;
2165:     PetscReal       norm, sign;

2167:     /* compute the norm based upon the support centroids */
2168:     DMPlexGetSupportSize(dm, cell, &supportSize);
2169:     DMPlexGetSupport(dm, cell, &support);
2170:     DMPlexComputeCellGeometryFVM(dm, support[0], NULL, normal, NULL);

2172:     /* Take the normal from the centroid of the support to the vertex*/
2173:     PetscSectionGetDof(coordSection, cell, &dof);
2174:     PetscSectionGetOffset(coordSection, cell, &off);
2175:     for (d = 0; d < dof; d++) normal[d] -= PetscRealPart(coords[off + d]);

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

2181:     norm = DMPlex_NormD_Internal(dim, normal);
2182:     for (d = 0; d < dim; ++d) normal[d] /= (norm * sign);
2183:   }
2184:   if (vol) *vol = 1.0;
2185:   VecRestoreArrayRead(coordinates, &coords);
2186:   return 0;
2187: }

2189: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2190: {
2191:   const PetscScalar *array;
2192:   PetscScalar       *coords = NULL;
2193:   PetscInt           cdim, coordSize, d;
2194:   PetscBool          isDG;

2196:   DMGetCoordinateDim(dm, &cdim);
2197:   DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords);
2199:   if (centroid) {
2200:     for (d = 0; d < cdim; ++d) centroid[d] = 0.5 * PetscRealPart(coords[d] + coords[cdim + d]);
2201:   }
2202:   if (normal) {
2203:     PetscReal norm;

2205:     switch (cdim) {
2206:     case 3:
2207:       normal[2] = 0.; /* fall through */
2208:     case 2:
2209:       normal[0] = -PetscRealPart(coords[1] - coords[cdim + 1]);
2210:       normal[1] = PetscRealPart(coords[0] - coords[cdim + 0]);
2211:       break;
2212:     case 1:
2213:       normal[0] = 1.0;
2214:       break;
2215:     default:
2216:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Dimension %" PetscInt_FMT " not supported", cdim);
2217:     }
2218:     norm = DMPlex_NormD_Internal(cdim, normal);
2219:     for (d = 0; d < cdim; ++d) normal[d] /= norm;
2220:   }
2221:   if (vol) {
2222:     *vol = 0.0;
2223:     for (d = 0; d < cdim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - coords[cdim + d]));
2224:     *vol = PetscSqrtReal(*vol);
2225:   }
2226:   DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords);
2227:   return 0;
2228: }

2230: /* Centroid_i = (\sum_n A_n Cn_i) / A */
2231: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2232: {
2233:   DMPolytopeType     ct;
2234:   const PetscScalar *array;
2235:   PetscScalar       *coords = NULL;
2236:   PetscInt           coordSize;
2237:   PetscBool          isDG;
2238:   PetscInt           fv[4] = {0, 1, 2, 3};
2239:   PetscInt           cdim, numCorners, p, d;

2241:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2242:   DMPlexGetCellType(dm, cell, &ct);
2243:   switch (ct) {
2244:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2245:     fv[2] = 3;
2246:     fv[3] = 2;
2247:     break;
2248:   default:
2249:     break;
2250:   }
2251:   DMGetCoordinateDim(dm, &cdim);
2252:   DMPlexGetConeSize(dm, cell, &numCorners);
2253:   DMPlexGetCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords);
2254:   {
2255:     PetscReal c[3] = {0., 0., 0.}, n[3] = {0., 0., 0.}, origin[3] = {0., 0., 0.}, norm;

2257:     for (d = 0; d < cdim; d++) origin[d] = PetscRealPart(coords[d]);
2258:     for (p = 0; p < numCorners - 2; ++p) {
2259:       PetscReal e0[3] = {0., 0., 0.}, e1[3] = {0., 0., 0.};
2260:       for (d = 0; d < cdim; d++) {
2261:         e0[d] = PetscRealPart(coords[cdim * fv[p + 1] + d]) - origin[d];
2262:         e1[d] = PetscRealPart(coords[cdim * fv[p + 2] + d]) - origin[d];
2263:       }
2264:       const PetscReal dx = e0[1] * e1[2] - e0[2] * e1[1];
2265:       const PetscReal dy = e0[2] * e1[0] - e0[0] * e1[2];
2266:       const PetscReal dz = e0[0] * e1[1] - e0[1] * e1[0];
2267:       const PetscReal a  = PetscSqrtReal(dx * dx + dy * dy + dz * dz);

2269:       n[0] += dx;
2270:       n[1] += dy;
2271:       n[2] += dz;
2272:       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.;
2273:     }
2274:     norm = PetscSqrtReal(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2275:     n[0] /= norm;
2276:     n[1] /= norm;
2277:     n[2] /= norm;
2278:     c[0] /= norm;
2279:     c[1] /= norm;
2280:     c[2] /= norm;
2281:     if (vol) *vol = 0.5 * norm;
2282:     if (centroid)
2283:       for (d = 0; d < cdim; ++d) centroid[d] = c[d];
2284:     if (normal)
2285:       for (d = 0; d < cdim; ++d) normal[d] = n[d];
2286:   }
2287:   DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords);
2288:   return 0;
2289: }

2291: /* Centroid_i = (\sum_n V_n Cn_i) / V */
2292: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2293: {
2294:   DMPolytopeType        ct;
2295:   const PetscScalar    *array;
2296:   PetscScalar          *coords = NULL;
2297:   PetscInt              coordSize;
2298:   PetscBool             isDG;
2299:   PetscReal             vsum      = 0.0, vtmp, coordsTmp[3 * 3], origin[3];
2300:   const PetscInt        order[16] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
2301:   const PetscInt       *cone, *faceSizes, *faces;
2302:   const DMPolytopeType *faceTypes;
2303:   PetscBool             isHybrid = PETSC_FALSE;
2304:   PetscInt              numFaces, f, fOff = 0, p, d;

2307:   /* Must check for hybrid cells because prisms have a different orientation scheme */
2308:   DMPlexGetCellType(dm, cell, &ct);
2309:   switch (ct) {
2310:   case DM_POLYTOPE_POINT_PRISM_TENSOR:
2311:   case DM_POLYTOPE_SEG_PRISM_TENSOR:
2312:   case DM_POLYTOPE_TRI_PRISM_TENSOR:
2313:   case DM_POLYTOPE_QUAD_PRISM_TENSOR:
2314:     isHybrid = PETSC_TRUE;
2315:   default:
2316:     break;
2317:   }

2319:   if (centroid)
2320:     for (d = 0; d < dim; ++d) centroid[d] = 0.0;
2321:   DMPlexGetCone(dm, cell, &cone);

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

2329:     // If using zero as the origin vertex for each tetrahedron, an element far from the origin will have positive and
2330:     // negative volumes that nearly cancel, thus incurring rounding error. Here we define origin[] as the first vertex
2331:     // so that all tetrahedra have positive volume.
2332:     if (f == 0)
2333:       for (d = 0; d < dim; d++) origin[d] = PetscRealPart(coords[d]);
2334:     switch (faceTypes[f]) {
2335:     case DM_POLYTOPE_TRIANGLE:
2336:       for (d = 0; d < dim; ++d) {
2337:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + 0] * dim + d]) - origin[d];
2338:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + 1] * dim + d]) - origin[d];
2339:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + 2] * dim + d]) - origin[d];
2340:       }
2341:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2342:       if (flip) vtmp = -vtmp;
2343:       vsum += vtmp;
2344:       if (centroid) { /* Centroid of OABC = (a+b+c)/4 */
2345:         for (d = 0; d < dim; ++d) {
2346:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2347:         }
2348:       }
2349:       break;
2350:     case DM_POLYTOPE_QUADRILATERAL:
2351:     case DM_POLYTOPE_SEG_PRISM_TENSOR: {
2352:       PetscInt fv[4] = {0, 1, 2, 3};

2354:       /* Side faces for hybrid cells are are stored as tensor products */
2355:       if (isHybrid && f > 1) {
2356:         fv[2] = 3;
2357:         fv[3] = 2;
2358:       }
2359:       /* DO FOR PYRAMID */
2360:       /* First tet */
2361:       for (d = 0; d < dim; ++d) {
2362:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[0]] * dim + d]) - origin[d];
2363:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2364:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2365:       }
2366:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2367:       if (flip) vtmp = -vtmp;
2368:       vsum += vtmp;
2369:       if (centroid) {
2370:         for (d = 0; d < dim; ++d) {
2371:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2372:         }
2373:       }
2374:       /* Second tet */
2375:       for (d = 0; d < dim; ++d) {
2376:         coordsTmp[0 * dim + d] = PetscRealPart(coords[faces[fOff + fv[1]] * dim + d]) - origin[d];
2377:         coordsTmp[1 * dim + d] = PetscRealPart(coords[faces[fOff + fv[2]] * dim + d]) - origin[d];
2378:         coordsTmp[2 * dim + d] = PetscRealPart(coords[faces[fOff + fv[3]] * dim + d]) - origin[d];
2379:       }
2380:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
2381:       if (flip) vtmp = -vtmp;
2382:       vsum += vtmp;
2383:       if (centroid) {
2384:         for (d = 0; d < dim; ++d) {
2385:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p * dim + d] * vtmp;
2386:         }
2387:       }
2388:       break;
2389:     }
2390:     default:
2391:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %" PetscInt_FMT " of type %s", cone[f], DMPolytopeTypes[ct]);
2392:     }
2393:     fOff += faceSizes[f];
2394:   }
2395:   DMPlexRestoreRawFaces_Internal(dm, ct, order, &numFaces, &faceTypes, &faceSizes, &faces);
2396:   DMPlexRestoreCellCoordinates(dm, cell, &isDG, &coordSize, &array, &coords);
2397:   if (vol) *vol = PetscAbsReal(vsum);
2398:   if (normal)
2399:     for (d = 0; d < dim; ++d) normal[d] = 0.0;
2400:   if (centroid)
2401:     for (d = 0; d < dim; ++d) centroid[d] = centroid[d] / (vsum * 4) + origin[d];
2402:   return 0;
2403: }

2405: /*@C
2406:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

2408:   Collective on dm

2410:   Input Parameters:
2411: + dm   - the DM
2412: - cell - the cell

2414:   Output Parameters:
2415: + volume   - the cell volume
2416: . centroid - the cell centroid
2417: - normal - the cell normal, if appropriate

2419:   Level: advanced

2421:   Fortran Notes:
2422:   Since it returns arrays, this routine is only available in Fortran 90, and you must
2423:   include petsc.h90 in your code.

2425: .seealso: `DMGetCoordinateSection()`, `DMGetCoordinates()`
2426: @*/
2427: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2428: {
2429:   PetscInt depth, dim;

2431:   DMPlexGetDepth(dm, &depth);
2432:   DMGetDimension(dm, &dim);
2434:   DMPlexGetPointDepth(dm, cell, &depth);
2435:   switch (depth) {
2436:   case 0:
2437:     DMPlexComputeGeometryFVM_0D_Internal(dm, dim, cell, vol, centroid, normal);
2438:     break;
2439:   case 1:
2440:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2441:     break;
2442:   case 2:
2443:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2444:     break;
2445:   case 3:
2446:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2447:     break;
2448:   default:
2449:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %" PetscInt_FMT " (depth %" PetscInt_FMT ") for element geometry computation", dim, depth);
2450:   }
2451:   return 0;
2452: }

2454: /*@
2455:   DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh

2457:   Collective on dm

2459:   Input Parameter:
2460: . dm - The DMPlex

2462:   Output Parameter:
2463: . cellgeom - A vector with the cell geometry data for each cell

2465:   Level: beginner

2467: @*/
2468: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2469: {
2470:   DM           dmCell;
2471:   Vec          coordinates;
2472:   PetscSection coordSection, sectionCell;
2473:   PetscScalar *cgeom;
2474:   PetscInt     cStart, cEnd, c;

2476:   DMClone(dm, &dmCell);
2477:   DMGetCoordinateSection(dm, &coordSection);
2478:   DMGetCoordinatesLocal(dm, &coordinates);
2479:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2480:   DMSetCoordinatesLocal(dmCell, coordinates);
2481:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);
2482:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2483:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2484:   /* TODO This needs to be multiplied by Nq for non-affine */
2485:   for (c = cStart; c < cEnd; ++c) PetscSectionSetDof(sectionCell, c, (PetscInt)PetscCeilReal(((PetscReal)sizeof(PetscFEGeom)) / sizeof(PetscScalar)));
2486:   PetscSectionSetUp(sectionCell);
2487:   DMSetLocalSection(dmCell, sectionCell);
2488:   PetscSectionDestroy(&sectionCell);
2489:   DMCreateLocalVector(dmCell, cellgeom);
2490:   VecGetArray(*cellgeom, &cgeom);
2491:   for (c = cStart; c < cEnd; ++c) {
2492:     PetscFEGeom *cg;

2494:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2495:     PetscArrayzero(cg, 1);
2496:     DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2498:   }
2499:   return 0;
2500: }

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

2505:   Input Parameter:
2506: . dm - The DM

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

2512:   Level: developer

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

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

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

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

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

2623:     DMPlexGetConeSize(dmCell, c, &coneSize);
2625:     DMPlexGetCone(dmCell, c, &cone);
2626:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2628:     DMPlexGetSupport(dmCell, cone[0], &support);
2629:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2630:     for (s = 0; s < 2; ++s) {
2631:       /* Reflect ghost centroid across plane of face */
2632:       if (support[s] == c) {
2633:         PetscFVCellGeom *ci;
2634:         PetscFVCellGeom *cg;
2635:         PetscReal        c2f[3], a;

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

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

2656:   Not collective

2658:   Input Parameter:
2659: . dm - the DM

2661:   Output Parameter:
2662: . minradius - the minimum cell radius

2664:   Level: developer

2666: .seealso: `DMGetCoordinates()`
2667: @*/
2668: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2669: {
2672:   *minradius = ((DM_Plex *)dm->data)->minradius;
2673:   return 0;
2674: }

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

2679:   Logically collective

2681:   Input Parameters:
2682: + dm - the DM
2683: - minradius - the minimum cell radius

2685:   Level: developer

2687: .seealso: `DMSetCoordinates()`
2688: @*/
2689: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2690: {
2692:   ((DM_Plex *)dm->data)->minradius = minradius;
2693:   return 0;
2694: }

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

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

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

2721:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2722:     DMPlexGetConeSize(dm, c, &numFaces);
2723:     DMPlexGetCone(dm, c, &faces);
2725:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2726:       PetscFVCellGeom *cg1;
2727:       PetscFVFaceGeom *fg;
2728:       const PetscInt  *fcells;
2729:       PetscInt         ncell, side;

2731:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2732:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2733:       if ((ghost >= 0) || boundary) continue;
2734:       DMPlexGetSupport(dm, faces[f], &fcells);
2735:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2736:       ncell = fcells[!side];    /* the neighbor */
2737:       DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2738:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2739:       for (d = 0; d < dim; ++d) dx[usedFaces * dim + d] = cg1->centroid[d] - cg->centroid[d];
2740:       gref[usedFaces++] = fg->grad[side]; /* Gradient reconstruction term will go here */
2741:     }
2743:     PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2744:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2745:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2746:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2747:       if ((ghost >= 0) || boundary) continue;
2748:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces * dim + d];
2749:       ++usedFaces;
2750:     }
2751:   }
2752:   PetscFree3(dx, grad, gref);
2753:   return 0;
2754: }

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

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

2779:     if (ghostLabel) DMLabelGetValue(ghostLabel, f, &ghost);
2780:     DMIsBoundaryPoint(dm, f, &boundary);
2781:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2782:     if ((ghost >= 0) || boundary || numChildren) continue;
2783:     DMPlexGetSupportSize(dm, f, &numCells);
2784:     if (numCells == 2) {
2785:       DMPlexGetSupport(dm, f, &fcells);
2786:       for (c = 0; c < 2; c++) {
2787:         PetscInt cell = fcells[c];

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

2806:     if (ghostLabel) DMLabelGetValue(ghostLabel, f, &ghost);
2807:     DMIsBoundaryPoint(dm, f, &boundary);
2808:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2809:     if ((ghost >= 0) || boundary || numChildren) continue;
2810:     DMPlexGetSupportSize(dm, f, &numCells);
2811:     if (numCells == 2) {
2812:       DMPlexGetSupport(dm, f, &fcells);
2813:       for (c = 0; c < 2; c++) {
2814:         PetscInt cell = fcells[c], off;

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

2831:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2832:     PetscSectionGetDof(neighSec, c, &numFaces);
2833:     PetscSectionGetOffset(neighSec, c, &off);

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

2840:     for (f = 0; f < numFaces; ++f) {
2841:       PetscFVCellGeom *cg1;
2842:       PetscFVFaceGeom *fg;
2843:       const PetscInt  *fcells;
2844:       PetscInt         ncell, side, nface;

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

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

2869:   Collective on dm

2871:   Input Parameters:
2872: + dm  - The DM
2873: . fvm - The PetscFV
2874: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()

2876:   Input/Output Parameter:
2877: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM(); on output
2878:                  the geometric factors for gradient calculation are inserted

2880:   Output Parameter:
2881: . dmGrad - The DM describing the layout of gradient data

2883:   Level: developer

2885: .seealso: `DMPlexGetFaceGeometryFVM()`, `DMPlexGetCellGeometryFVM()`
2886: @*/
2887: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2888: {
2889:   DM           dmFace, dmCell;
2890:   PetscScalar *fgeom, *cgeom;
2891:   PetscSection sectionGrad, parentSection;
2892:   PetscInt     dim, pdim, cStart, cEnd, cEndInterior, c;

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

2922: /*@
2923:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2925:   Collective on dm

2927:   Input Parameters:
2928: + dm  - The DM
2929: - fv  - The PetscFV

2931:   Output Parameters:
2932: + cellGeometry - The cell geometry
2933: . faceGeometry - The face geometry
2934: - gradDM       - The gradient matrices

2936:   Level: developer

2938: .seealso: `DMPlexComputeGeometryFVM()`
2939: @*/
2940: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2941: {
2942:   PetscObject cellgeomobj, facegeomobj;

2944:   PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2945:   if (!cellgeomobj) {
2946:     Vec cellgeomInt, facegeomInt;

2948:     DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2949:     PetscObjectCompose((PetscObject)dm, "DMPlex_cellgeom_fvm", (PetscObject)cellgeomInt);
2950:     PetscObjectCompose((PetscObject)dm, "DMPlex_facegeom_fvm", (PetscObject)facegeomInt);
2951:     VecDestroy(&cellgeomInt);
2952:     VecDestroy(&facegeomInt);
2953:     PetscObjectQuery((PetscObject)dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2954:   }
2955:   PetscObjectQuery((PetscObject)dm, "DMPlex_facegeom_fvm", &facegeomobj);
2956:   if (cellgeom) *cellgeom = (Vec)cellgeomobj;
2957:   if (facegeom) *facegeom = (Vec)facegeomobj;
2958:   if (gradDM) {
2959:     PetscObject gradobj;
2960:     PetscBool   computeGradients;

2962:     PetscFVGetComputeGradients(fv, &computeGradients);
2963:     if (!computeGradients) {
2964:       *gradDM = NULL;
2965:       return 0;
2966:     }
2967:     PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj);
2968:     if (!gradobj) {
2969:       DM dmGradInt;

2971:       DMPlexComputeGradientFVM(dm, fv, (Vec)facegeomobj, (Vec)cellgeomobj, &dmGradInt);
2972:       PetscObjectCompose((PetscObject)dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2973:       DMDestroy(&dmGradInt);
2974:       PetscObjectQuery((PetscObject)dm, "DMPlex_dmgrad_fvm", &gradobj);
2975:     }
2976:     *gradDM = (DM)gradobj;
2977:   }
2978:   return 0;
2979: }

2981: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work, PetscReal *resNeg, PetscReal *guess)
2982: {
2983:   PetscInt l, m;

2986:   if (dimC == dimR && dimR <= 3) {
2987:     /* invert Jacobian, multiply */
2988:     PetscScalar det, idet;

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

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

3035:     PetscCallBLAS("LAPACKgels", LAPACKgels_(&transpose, &m, &n, &one, J, &m, invJ, &n, work, &worksize, &info));

3038:     for (l = 0; l < dimR; l++) guess[l] += PetscRealPart(invJ[l]);
3039:   }
3040:   return 0;
3041: }

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

3051:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
3053:   DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
3054:   DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
3055:   cellCoords = &cellData[0];
3056:   cellCoeffs = &cellData[coordSize];
3057:   extJ       = &cellData[2 * coordSize];
3058:   resNeg     = &cellData[2 * coordSize + dimR];
3059:   invJ       = &J[dimR * dimC];
3060:   work       = &J[2 * dimR * dimC];
3061:   if (dimR == 2) {
3062:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

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

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

3075:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3076:     }
3077:   } else {
3078:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3079:   }
3080:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3081:   for (i = 0; i < dimR; i++) {
3082:     PetscReal *swap;

3084:     for (j = 0; j < (numV / 2); j++) {
3085:       for (k = 0; k < dimC; k++) {
3086:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3087:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3088:       }
3089:     }

3091:     if (i < dimR - 1) {
3092:       swap       = cellCoeffs;
3093:       cellCoeffs = cellCoords;
3094:       cellCoords = swap;
3095:     }
3096:   }
3097:   PetscArrayzero(refCoords, numPoints * dimR);
3098:   for (j = 0; j < numPoints; j++) {
3099:     for (i = 0; i < maxIts; i++) {
3100:       PetscReal *guess = &refCoords[dimR * j];

3102:       /* compute -residual and Jacobian */
3103:       for (k = 0; k < dimC; k++) resNeg[k] = realCoords[dimC * j + k];
3104:       for (k = 0; k < dimC * dimR; k++) J[k] = 0.;
3105:       for (k = 0; k < numV; k++) {
3106:         PetscReal extCoord = 1.;
3107:         for (l = 0; l < dimR; l++) {
3108:           PetscReal coord = guess[l];
3109:           PetscInt  dep   = (k & (1 << l)) >> l;

3111:           extCoord *= dep * coord + !dep;
3112:           extJ[l] = dep;

3114:           for (m = 0; m < dimR; m++) {
3115:             PetscReal coord = guess[m];
3116:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
3117:             PetscReal mult  = dep * coord + !dep;

3119:             extJ[l] *= mult;
3120:           }
3121:         }
3122:         for (l = 0; l < dimC; l++) {
3123:           PetscReal coeff = cellCoeffs[dimC * k + l];

3125:           resNeg[l] -= coeff * extCoord;
3126:           for (m = 0; m < dimR; m++) J[dimR * l + m] += coeff * extJ[m];
3127:         }
3128:       }
3129:       if (0 && PetscDefined(USE_DEBUG)) {
3130:         PetscReal maxAbs = 0.;

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

3136:       DMPlexCoordinatesToReference_NewtonUpdate(dimC, dimR, J, invJ, work, resNeg, guess);
3137:     }
3138:   }
3139:   DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
3140:   DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
3141:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
3142:   return 0;
3143: }

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

3152:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
3154:   DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
3155:   cellCoords = &cellData[0];
3156:   cellCoeffs = &cellData[coordSize];
3157:   if (dimR == 2) {
3158:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

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

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

3171:       for (j = 0; j < dimC; j++) cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
3172:     }
3173:   } else {
3174:     for (i = 0; i < coordSize; i++) cellCoords[i] = PetscRealPart(coordsScalar[i]);
3175:   }
3176:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
3177:   for (i = 0; i < dimR; i++) {
3178:     PetscReal *swap;

3180:     for (j = 0; j < (numV / 2); j++) {
3181:       for (k = 0; k < dimC; k++) {
3182:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
3183:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
3184:       }
3185:     }

3187:     if (i < dimR - 1) {
3188:       swap       = cellCoeffs;
3189:       cellCoeffs = cellCoords;
3190:       cellCoords = swap;
3191:     }
3192:   }
3193:   PetscArrayzero(realCoords, numPoints * dimC);
3194:   for (j = 0; j < numPoints; j++) {
3195:     const PetscReal *guess  = &refCoords[dimR * j];
3196:     PetscReal       *mapped = &realCoords[dimC * j];

3198:     for (k = 0; k < numV; k++) {
3199:       PetscReal extCoord = 1.;
3200:       for (l = 0; l < dimR; l++) {
3201:         PetscReal coord = guess[l];
3202:         PetscInt  dep   = (k & (1 << l)) >> l;

3204:         extCoord *= dep * coord + !dep;
3205:       }
3206:       for (l = 0; l < dimC; l++) {
3207:         PetscReal coeff = cellCoeffs[dimC * k + l];

3209:         mapped[l] += coeff * extCoord;
3210:       }
3211:     }
3212:   }
3213:   DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
3214:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
3215:   return 0;
3216: }

3218: /* TODO: TOBY please fix this for Nc > 1 */
3219: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3220: {
3221:   PetscInt     numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
3222:   PetscScalar *nodes = NULL;
3223:   PetscReal   *invV, *modes;
3224:   PetscReal   *B, *D, *resNeg;
3225:   PetscScalar *J, *invJ, *work;

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

3260:         for (l = 0; l < Nc; l++) maxAbs = PetscMax(maxAbs, PetscAbsReal(resNeg[l]));
3261:         PetscInfo(dm, "cell %" PetscInt_FMT ", point %" PetscInt_FMT ", iter %" PetscInt_FMT ": res %g\n", cell, j, i, (double)maxAbs);
3262:       }
3263:       DMPlexCoordinatesToReference_NewtonUpdate(Nc, dimR, J, invJ, work, resNeg, guess);
3264:     }
3265:   }
3266:   DMRestoreWorkArray(dm, 3 * Nc * dimR, MPIU_SCALAR, &J);
3267:   DMRestoreWorkArray(dm, pdim * Nc + pdim * Nc * dimR + Nc, MPIU_REAL, &B);
3268:   DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes);
3269:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3270:   return 0;
3271: }

3273: /* TODO: TOBY please fix this for Nc > 1 */
3274: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
3275: {
3276:   PetscInt     numComp, pdim, i, j, k, l, coordSize;
3277:   PetscScalar *nodes = NULL;
3278:   PetscReal   *invV, *modes;
3279:   PetscReal   *B;

3281:   PetscFEGetDimension(fe, &pdim);
3282:   PetscFEGetNumComponents(fe, &numComp);
3284:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3285:   /* convert nodes to values in the stable evaluation basis */
3286:   DMGetWorkArray(dm, pdim, MPIU_REAL, &modes);
3287:   invV = fe->invV;
3288:   for (i = 0; i < pdim; ++i) {
3289:     modes[i] = 0.;
3290:     for (j = 0; j < pdim; ++j) modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
3291:   }
3292:   DMGetWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B);
3293:   PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
3294:   for (i = 0; i < numPoints * Nc; i++) realCoords[i] = 0.;
3295:   for (j = 0; j < numPoints; j++) {
3296:     PetscReal *mapped = &realCoords[j * Nc];

3298:     for (k = 0; k < pdim; k++) {
3299:       for (l = 0; l < Nc; l++) mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
3300:     }
3301:   }
3302:   DMRestoreWorkArray(dm, numPoints * pdim * Nc, MPIU_REAL, &B);
3303:   DMRestoreWorkArray(dm, pdim, MPIU_REAL, &modes);
3304:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
3305:   return 0;
3306: }

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

3313:   Not collective

3315:   Input Parameters:
3316: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3317:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3318:                as a multilinear map for tensor-product elements
3319: . cell       - the cell whose map is used.
3320: . numPoints  - the number of points to locate
3321: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

3323:   Output Parameters:
3324: . refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

3326:   Level: intermediate

3328: .seealso: `DMPlexReferenceToCoordinates()`
3329: @*/
3330: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
3331: {
3332:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3333:   DM       coordDM = NULL;
3334:   Vec      coords;
3335:   PetscFE  fe = NULL;

3338:   DMGetDimension(dm, &dimR);
3339:   DMGetCoordinateDim(dm, &dimC);
3340:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return 0;
3341:   DMPlexGetDepth(dm, &depth);
3342:   DMGetCoordinatesLocal(dm, &coords);
3343:   DMGetCoordinateDM(dm, &coordDM);
3344:   if (coordDM) {
3345:     PetscInt coordFields;

3347:     DMGetNumFields(coordDM, &coordFields);
3348:     if (coordFields) {
3349:       PetscClassId id;
3350:       PetscObject  disc;

3352:       DMGetField(coordDM, 0, NULL, &disc);
3353:       PetscObjectGetClassId(disc, &id);
3354:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3355:     }
3356:   }
3357:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3359:   if (!fe) { /* implicit discretization: affine or multilinear */
3360:     PetscInt  coneSize;
3361:     PetscBool isSimplex, isTensor;

3363:     DMPlexGetConeSize(dm, cell, &coneSize);
3364:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3365:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3366:     if (isSimplex) {
3367:       PetscReal detJ, *v0, *J, *invJ;

3369:       DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3370:       J    = &v0[dimC];
3371:       invJ = &J[dimC * dimC];
3372:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3373:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3374:         const PetscReal x0[3] = {-1., -1., -1.};

3376:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3377:       }
3378:       DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3379:     } else if (isTensor) {
3380:       DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3381:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3382:   } else {
3383:     DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3384:   }
3385:   return 0;
3386: }

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

3391:   Not collective

3393:   Input Parameters:
3394: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3395:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3396:                as a multilinear map for tensor-product elements
3397: . cell       - the cell whose map is used.
3398: . numPoints  - the number of points to locate
3399: - refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

3401:   Output Parameters:
3402: . realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

3404:    Level: intermediate

3406: .seealso: `DMPlexCoordinatesToReference()`
3407: @*/
3408: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3409: {
3410:   PetscInt dimC, dimR, depth, cStart, cEnd, i;
3411:   DM       coordDM = NULL;
3412:   Vec      coords;
3413:   PetscFE  fe = NULL;

3416:   DMGetDimension(dm, &dimR);
3417:   DMGetCoordinateDim(dm, &dimC);
3418:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return 0;
3419:   DMPlexGetDepth(dm, &depth);
3420:   DMGetCoordinatesLocal(dm, &coords);
3421:   DMGetCoordinateDM(dm, &coordDM);
3422:   if (coordDM) {
3423:     PetscInt coordFields;

3425:     DMGetNumFields(coordDM, &coordFields);
3426:     if (coordFields) {
3427:       PetscClassId id;
3428:       PetscObject  disc;

3430:       DMGetField(coordDM, 0, NULL, &disc);
3431:       PetscObjectGetClassId(disc, &id);
3432:       if (id == PETSCFE_CLASSID) fe = (PetscFE)disc;
3433:     }
3434:   }
3435:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3437:   if (!fe) { /* implicit discretization: affine or multilinear */
3438:     PetscInt  coneSize;
3439:     PetscBool isSimplex, isTensor;

3441:     DMPlexGetConeSize(dm, cell, &coneSize);
3442:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3443:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3444:     if (isSimplex) {
3445:       PetscReal detJ, *v0, *J;

3447:       DMGetWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3448:       J = &v0[dimC];
3449:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, NULL, &detJ);
3450:       for (i = 0; i < numPoints; i++) { /* Apply the affine transformation for each point */
3451:         const PetscReal xi0[3] = {-1., -1., -1.};

3453:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3454:       }
3455:       DMRestoreWorkArray(dm, dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3456:     } else if (isTensor) {
3457:       DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3458:     } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unrecognized cone size %" PetscInt_FMT, coneSize);
3459:   } else {
3460:     DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3461:   }
3462:   return 0;
3463: }

3465: /*@C
3466:   DMPlexRemapGeometry - This function maps the original DM coordinates to new coordinates.

3468:   Not collective

3470:   Input Parameters:
3471: + dm      - The DM
3472: . time    - The time
3473: - func    - The function transforming current coordinates to new coordaintes

3475:    Calling sequence of func:
3476: $    func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3477: $         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3478: $         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3479: $         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);

3481: +  dim          - The spatial dimension
3482: .  Nf           - The number of input fields (here 1)
3483: .  NfAux        - The number of input auxiliary fields
3484: .  uOff         - The offset of the coordinates in u[] (here 0)
3485: .  uOff_x       - The offset of the coordinates in u_x[] (here 0)
3486: .  u            - The coordinate values at this point in space
3487: .  u_t          - The coordinate time derivative at this point in space (here NULL)
3488: .  u_x          - The coordinate derivatives at this point in space
3489: .  aOff         - The offset of each auxiliary field in u[]
3490: .  aOff_x       - The offset of each auxiliary field in u_x[]
3491: .  a            - The auxiliary field values at this point in space
3492: .  a_t          - The auxiliary field time derivative at this point in space (or NULL)
3493: .  a_x          - The auxiliary field derivatives at this point in space
3494: .  t            - The current time
3495: .  x            - The coordinates of this point (here not used)
3496: .  numConstants - The number of constants
3497: .  constants    - The value of each constant
3498: -  f            - The new coordinates at this point in space

3500:   Level: intermediate

3502: .seealso: `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMProjectFieldLocal()`, `DMProjectFieldLabelLocal()`
3503: @*/
3504: 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[]))
3505: {
3506:   DM      cdm;
3507:   DMField cf;
3508:   Vec     lCoords, tmpCoords;

3510:   DMGetCoordinateDM(dm, &cdm);
3511:   DMGetCoordinatesLocal(dm, &lCoords);
3512:   DMGetLocalVector(cdm, &tmpCoords);
3513:   VecCopy(lCoords, tmpCoords);
3514:   /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3515:   DMGetCoordinateField(dm, &cf);
3516:   cdm->coordinates[0].field = cf;
3517:   DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3518:   cdm->coordinates[0].field = NULL;
3519:   DMRestoreLocalVector(cdm, &tmpCoords);
3520:   DMSetCoordinatesLocal(dm, lCoords);
3521:   return 0;
3522: }

3524: /* Shear applies the transformation, assuming we fix z,
3525:   / 1  0  m_0 \
3526:   | 0  1  m_1 |
3527:   \ 0  0   1  /
3528: */
3529: 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[])
3530: {
3531:   const PetscInt Nc = uOff[1] - uOff[0];
3532:   const PetscInt ax = (PetscInt)PetscRealPart(constants[0]);
3533:   PetscInt       c;

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

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

3541:   Not collective

3543:   Input Parameters:
3544: + dm          - The DM
3545: . direction   - The shear coordinate direction, e.g. 0 is the x-axis
3546: - multipliers - The multiplier m for each direction which is not the shear direction

3548:   Level: intermediate

3550: .seealso: `DMPlexRemapGeometry()`
3551: @*/
3552: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3553: {
3554:   DM             cdm;
3555:   PetscDS        cds;
3556:   PetscObject    obj;
3557:   PetscClassId   id;
3558:   PetscScalar   *moduli;
3559:   const PetscInt dir = (PetscInt)direction;
3560:   PetscInt       dE, d, e;

3562:   DMGetCoordinateDM(dm, &cdm);
3563:   DMGetCoordinateDim(dm, &dE);
3564:   PetscMalloc1(dE + 1, &moduli);
3565:   moduli[0] = dir;
3566:   for (d = 0, e = 0; d < dE; ++d) moduli[d + 1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3567:   DMGetDS(cdm, &cds);
3568:   PetscDSGetDiscretization(cds, 0, &obj);
3569:   PetscObjectGetClassId(obj, &id);
3570:   if (id != PETSCFE_CLASSID) {
3571:     Vec          lCoords;
3572:     PetscSection cSection;
3573:     PetscScalar *coords;
3574:     PetscInt     vStart, vEnd, v;

3576:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3577:     DMGetCoordinateSection(dm, &cSection);
3578:     DMGetCoordinatesLocal(dm, &lCoords);
3579:     VecGetArray(lCoords, &coords);
3580:     for (v = vStart; v < vEnd; ++v) {
3581:       PetscReal ds;
3582:       PetscInt  off, c;

3584:       PetscSectionGetOffset(cSection, v, &off);
3585:       ds = PetscRealPart(coords[off + dir]);
3586:       for (c = 0; c < dE; ++c) coords[off + c] += moduli[c] * ds;
3587:     }
3588:     VecRestoreArray(lCoords, &coords);
3589:   } else {
3590:     PetscDSSetConstants(cds, dE + 1, moduli);
3591:     DMPlexRemapGeometry(dm, 0.0, f0_shear);
3592:   }
3593:   PetscFree(moduli);
3594:   return 0;
3595: }