Actual source code: plexgeometry.c

petsc-master 2020-08-09
Report Typos and Errors
  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: . npoints - The number of sought points
 14: . coords - The array of coordinates of the sought points
 15: - eps - The tolerance or PETSC_DEFAULT

 17:   Output Parameters:
 18: . dagPoints - The array of found DAG points, or -1 if not found

 20:   Level: intermediate

 22:   Notes:
 23:   The length of the array coords must be npoints * dim where dim is the spatial dimension returned by DMGetDimension().

 25:   The output array dagPoints is NOT newly allocated; the user must pass an array of length npoints.

 27:   Each rank does the search independently; a nonnegative value is returned only if this rank's local DMPlex portion contains the point.

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

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

 33: .seealso: DMPlexCreate(), DMGetCoordinates()
 34: @*/
 35: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
 36: {
 37:   PetscInt          c, dim, i, j, o, p, vStart, vEnd;
 38:   Vec               allCoordsVec;
 39:   const PetscScalar *allCoords;
 40:   PetscReal         norm;
 41:   PetscErrorCode    ierr;

 44:   if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
 45:   DMGetDimension(dm, &dim);
 46:   DMGetCoordinatesLocal(dm, &allCoordsVec);
 47:   VecGetArrayRead(allCoordsVec, &allCoords);
 48:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 49:   if (PetscDefined(USE_DEBUG)) {
 50:     /* check coordinate section is consistent with DM dimension */
 51:     PetscSection      cs;
 52:     PetscInt          ndof;

 54:     DMGetCoordinateSection(dm, &cs);
 55:     for (p = vStart; p < vEnd; p++) {
 56:       PetscSectionGetDof(cs, p, &ndof);
 57:       if (PetscUnlikely(ndof != dim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = dim", p, ndof, dim);
 58:     }
 59:   }
 60:   if (eps == 0.0) {
 61:     for (i=0,j=0; i < npoints; i++,j+=dim) {
 62:       dagPoints[i] = -1;
 63:       for (p = vStart,o=0; p < vEnd; p++,o+=dim) {
 64:         for (c = 0; c < dim; c++) {
 65:           if (coord[j+c] != PetscRealPart(allCoords[o+c])) break;
 66:         }
 67:         if (c == dim) {
 68:           dagPoints[i] = p;
 69:           break;
 70:         }
 71:       }
 72:     }
 73:     VecRestoreArrayRead(allCoordsVec, &allCoords);
 74:     return(0);
 75:   }
 76:   for (i=0,j=0; i < npoints; i++,j+=dim) {
 77:     dagPoints[i] = -1;
 78:     for (p = vStart,o=0; p < vEnd; p++,o+=dim) {
 79:       norm = 0.0;
 80:       for (c = 0; c < dim; c++) {
 81:         norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
 82:       }
 83:       norm = PetscSqrtReal(norm);
 84:       if (norm <= eps) {
 85:         dagPoints[i] = p;
 86:         break;
 87:       }
 88:     }
 89:   }
 90:   VecRestoreArrayRead(allCoordsVec, &allCoords);
 91:   return(0);
 92: }

 94: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
 95: {
 96:   const PetscReal p0_x  = segmentA[0*2+0];
 97:   const PetscReal p0_y  = segmentA[0*2+1];
 98:   const PetscReal p1_x  = segmentA[1*2+0];
 99:   const PetscReal p1_y  = segmentA[1*2+1];
100:   const PetscReal p2_x  = segmentB[0*2+0];
101:   const PetscReal p2_y  = segmentB[0*2+1];
102:   const PetscReal p3_x  = segmentB[1*2+0];
103:   const PetscReal p3_y  = segmentB[1*2+1];
104:   const PetscReal s1_x  = p1_x - p0_x;
105:   const PetscReal s1_y  = p1_y - p0_y;
106:   const PetscReal s2_x  = p3_x - p2_x;
107:   const PetscReal s2_y  = p3_y - p2_y;
108:   const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

111:   *hasIntersection = PETSC_FALSE;
112:   /* Non-parallel lines */
113:   if (denom != 0.0) {
114:     const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
115:     const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;

117:     if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
118:       *hasIntersection = PETSC_TRUE;
119:       if (intersection) {
120:         intersection[0] = p0_x + (t * s1_x);
121:         intersection[1] = p0_y + (t * s1_y);
122:       }
123:     }
124:   }
125:   return(0);
126: }

128: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
129: {
130:   const PetscInt  embedDim = 2;
131:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
132:   PetscReal       x        = PetscRealPart(point[0]);
133:   PetscReal       y        = PetscRealPart(point[1]);
134:   PetscReal       v0[2], J[4], invJ[4], detJ;
135:   PetscReal       xi, eta;
136:   PetscErrorCode  ierr;

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

143:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
144:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
145:   return(0);
146: }

148: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
149: {
150:   const PetscInt  embedDim = 2;
151:   PetscReal       x        = PetscRealPart(point[0]);
152:   PetscReal       y        = PetscRealPart(point[1]);
153:   PetscReal       v0[2], J[4], invJ[4], detJ;
154:   PetscReal       xi, eta, r;
155:   PetscErrorCode  ierr;

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

162:   xi  = PetscMax(xi,  0.0);
163:   eta = PetscMax(eta, 0.0);
164:   if (xi + eta > 2.0) {
165:     r    = (xi + eta)/2.0;
166:     xi  /= r;
167:     eta /= r;
168:   }
169:   cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
170:   cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
171:   return(0);
172: }

174: static PetscErrorCode DMPlexLocatePoint_Quad_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
175: {
176:   PetscSection       coordSection;
177:   Vec             coordsLocal;
178:   PetscScalar    *coords = NULL;
179:   const PetscInt  faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
180:   PetscReal       x         = PetscRealPart(point[0]);
181:   PetscReal       y         = PetscRealPart(point[1]);
182:   PetscInt        crossings = 0, f;
183:   PetscErrorCode  ierr;

186:   DMGetCoordinatesLocal(dm, &coordsLocal);
187:   DMGetCoordinateSection(dm, &coordSection);
188:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
189:   for (f = 0; f < 4; ++f) {
190:     PetscReal x_i   = PetscRealPart(coords[faces[2*f+0]*2+0]);
191:     PetscReal y_i   = PetscRealPart(coords[faces[2*f+0]*2+1]);
192:     PetscReal x_j   = PetscRealPart(coords[faces[2*f+1]*2+0]);
193:     PetscReal y_j   = PetscRealPart(coords[faces[2*f+1]*2+1]);
194:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
195:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
196:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
197:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
198:     if ((cond1 || cond2)  && above) ++crossings;
199:   }
200:   if (crossings % 2) *cell = c;
201:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
202:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
203:   return(0);
204: }

206: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
207: {
208:   const PetscInt embedDim = 3;
209:   PetscReal      v0[3], J[9], invJ[9], detJ;
210:   PetscReal      x = PetscRealPart(point[0]);
211:   PetscReal      y = PetscRealPart(point[1]);
212:   PetscReal      z = PetscRealPart(point[2]);
213:   PetscReal      xi, eta, zeta;

217:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
218:   xi   = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
219:   eta  = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
220:   zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);

222:   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
223:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
224:   return(0);
225: }

227: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
228: {
229:   PetscSection   coordSection;
230:   Vec            coordsLocal;
231:   PetscScalar   *coords = NULL;
232:   const PetscInt faces[24] = {0, 3, 2, 1,  5, 4, 7, 6,  3, 0, 4, 5,
233:                               1, 2, 6, 7,  3, 5, 6, 2,  0, 1, 7, 4};
234:   PetscBool      found = PETSC_TRUE;
235:   PetscInt       f;

239:   DMGetCoordinatesLocal(dm, &coordsLocal);
240:   DMGetCoordinateSection(dm, &coordSection);
241:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
242:   for (f = 0; f < 6; ++f) {
243:     /* Check the point is under plane */
244:     /*   Get face normal */
245:     PetscReal v_i[3];
246:     PetscReal v_j[3];
247:     PetscReal normal[3];
248:     PetscReal pp[3];
249:     PetscReal dot;

251:     v_i[0]    = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
252:     v_i[1]    = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
253:     v_i[2]    = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
254:     v_j[0]    = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
255:     v_j[1]    = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
256:     v_j[2]    = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
257:     normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
258:     normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
259:     normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
260:     pp[0]     = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
261:     pp[1]     = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
262:     pp[2]     = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
263:     dot       = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];

265:     /* Check that projected point is in face (2D location problem) */
266:     if (dot < 0.0) {
267:       found = PETSC_FALSE;
268:       break;
269:     }
270:   }
271:   if (found) *cell = c;
272:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
273:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
274:   return(0);
275: }

277: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
278: {
279:   PetscInt d;

282:   box->dim = dim;
283:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
284:   return(0);
285: }

287: PetscErrorCode PetscGridHashCreate(MPI_Comm comm, PetscInt dim, const PetscScalar point[], PetscGridHash *box)
288: {

292:   PetscMalloc1(1, box);
293:   PetscGridHashInitialize_Internal(*box, dim, point);
294:   return(0);
295: }

297: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
298: {
299:   PetscInt d;

302:   for (d = 0; d < box->dim; ++d) {
303:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
304:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
305:   }
306:   return(0);
307: }

309: /*
310:   PetscGridHashSetGrid - Divide the grid into boxes

312:   Not collective

314:   Input Parameters:
315: + box - The grid hash object
316: . n   - The number of boxes in each dimension, or PETSC_DETERMINE
317: - h   - The box size in each dimension, only used if n[d] == PETSC_DETERMINE

319:   Level: developer

321: .seealso: PetscGridHashCreate()
322: */
323: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
324: {
325:   PetscInt d;

328:   for (d = 0; d < box->dim; ++d) {
329:     box->extent[d] = box->upper[d] - box->lower[d];
330:     if (n[d] == PETSC_DETERMINE) {
331:       box->h[d] = h[d];
332:       box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
333:     } else {
334:       box->n[d] = n[d];
335:       box->h[d] = box->extent[d]/n[d];
336:     }
337:   }
338:   return(0);
339: }

341: /*
342:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

344:   Not collective

346:   Input Parameters:
347: + box       - The grid hash object
348: . numPoints - The number of input points
349: - points    - The input point coordinates

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

355:   Level: developer

357: .seealso: PetscGridHashCreate()
358: */
359: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
360: {
361:   const PetscReal *lower = box->lower;
362:   const PetscReal *upper = box->upper;
363:   const PetscReal *h     = box->h;
364:   const PetscInt  *n     = box->n;
365:   const PetscInt   dim   = box->dim;
366:   PetscInt         d, p;

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

373:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
374:       if (dbox == -1   && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
375:       if (dbox < 0 || dbox >= n[d]) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input point %d (%g, %g, %g) is outside of our bounding box",
376:                                              p, (double) PetscRealPart(points[p*dim+0]), dim > 1 ? (double) PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? (double) PetscRealPart(points[p*dim+2]) : 0.0);
377:       dboxes[p*dim+d] = dbox;
378:     }
379:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
380:   }
381:   return(0);
382: }

384: /*
385:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

387:  Not collective

389:   Input Parameters:
390: + box       - The grid hash object
391: . numPoints - The number of input points
392: - points    - The input point coordinates

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

399:   Level: developer

401: .seealso: PetscGridHashGetEnclosingBox()
402: */
403: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
404: {
405:   const PetscReal *lower = box->lower;
406:   const PetscReal *upper = box->upper;
407:   const PetscReal *h     = box->h;
408:   const PetscInt  *n     = box->n;
409:   const PetscInt   dim   = box->dim;
410:   PetscInt         d, p;

413:   *found = PETSC_FALSE;
414:   for (p = 0; p < numPoints; ++p) {
415:     for (d = 0; d < dim; ++d) {
416:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);

418:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
419:       if (dbox < 0 || dbox >= n[d]) {
420:         return(0);
421:       }
422:       dboxes[p*dim+d] = dbox;
423:     }
424:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
425:   }
426:   *found = PETSC_TRUE;
427:   return(0);
428: }

430: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
431: {

435:   if (*box) {
436:     PetscSectionDestroy(&(*box)->cellSection);
437:     ISDestroy(&(*box)->cells);
438:     DMLabelDestroy(&(*box)->cellsSparse);
439:   }
440:   PetscFree(*box);
441:   return(0);
442: }

444: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
445: {
446:   DMPolytopeType ct;

450:   DMPlexGetCellType(dm, cellStart, &ct);
451:   switch (ct) {
452:     case DM_POLYTOPE_TRIANGLE:
453:     DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);break;
454:     case DM_POLYTOPE_QUADRILATERAL:
455:     DMPlexLocatePoint_Quad_2D_Internal(dm, point, cellStart, cell);break;
456:     case DM_POLYTOPE_TETRAHEDRON:
457:     DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);break;
458:     case DM_POLYTOPE_HEXAHEDRON:
459:     DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);break;
460:     default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell %D with type %s", cellStart, DMPolytopeTypes[ct]);
461:   }
462:   return(0);
463: }

465: /*
466:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
467: */
468: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
469: {
470:   DMPolytopeType ct;

474:   DMPlexGetCellType(dm, cell, &ct);
475:   switch (ct) {
476:     case DM_POLYTOPE_TRIANGLE:
477:     DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);break;
478: #if 0
479:     case DM_POLYTOPE_QUADRILATERAL:
480:     DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);break;
481:     case DM_POLYTOPE_TETRAHEDRON:
482:     DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);break;
483:     case DM_POLYTOPE_HEXAHEDRON:
484:     DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);break;
485: #endif
486:     default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell %D with type %s", cell, DMPolytopeTypes[ct]);
487:   }
488:   return(0);
489: }

491: /*
492:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex

494:   Collective on dm

496:   Input Parameter:
497: . dm - The Plex

499:   Output Parameter:
500: . localBox - The grid hash object

502:   Level: developer

504: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
505: */
506: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
507: {
508:   MPI_Comm           comm;
509:   PetscGridHash      lbox;
510:   Vec                coordinates;
511:   PetscSection       coordSection;
512:   Vec                coordsLocal;
513:   const PetscScalar *coords;
514:   PetscInt          *dboxes, *boxes;
515:   PetscInt           n[3] = {10, 10, 10};
516:   PetscInt           dim, N, cStart, cEnd, c, i;
517:   PetscErrorCode     ierr;

520:   PetscObjectGetComm((PetscObject) dm, &comm);
521:   DMGetCoordinatesLocal(dm, &coordinates);
522:   DMGetCoordinateDim(dm, &dim);
523:   if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
524:   VecGetLocalSize(coordinates, &N);
525:   VecGetArrayRead(coordinates, &coords);
526:   PetscGridHashCreate(comm, dim, coords, &lbox);
527:   for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
528:   VecRestoreArrayRead(coordinates, &coords);
529:   PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
530:   n[1] = n[0];
531:   n[2] = n[0];
532:   PetscGridHashSetGrid(lbox, n, NULL);
533: #if 0
534:   /* Could define a custom reduction to merge these */
535:   MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
536:   MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
537: #endif
538:   /* Is there a reason to snap the local bounding box to a division of the global box? */
539:   /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
540:   /* Create label */
541:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
542:   DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
543:   DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
544:   /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
545:   DMGetCoordinatesLocal(dm, &coordsLocal);
546:   DMGetCoordinateSection(dm, &coordSection);
547:   PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
548:   for (c = cStart; c < cEnd; ++c) {
549:     const PetscReal *h       = lbox->h;
550:     PetscScalar     *ccoords = NULL;
551:     PetscInt         csize   = 0;
552:     PetscScalar      point[3];
553:     PetscInt         dlim[6], d, e, i, j, k;

555:     /* Find boxes enclosing each vertex */
556:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
557:     PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
558:     /* Mark cells containing the vertices */
559:     for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
560:     /* Get grid of boxes containing these */
561:     for (d = 0;   d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
562:     for (d = dim; d < 3;   ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
563:     for (e = 1; e < dim+1; ++e) {
564:       for (d = 0; d < dim; ++d) {
565:         dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
566:         dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
567:       }
568:     }
569:     /* Check for intersection of box with cell */
570:     for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
571:       for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
572:         for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
573:           const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
574:           PetscScalar    cpoint[3];
575:           PetscInt       cell, edge, ii, jj, kk;

577:           /* Check whether cell contains any vertex of these subboxes TODO vectorize this */
578:           for (kk = 0, cpoint[2] = point[2]; kk < (dim > 2 ? 2 : 1); ++kk, cpoint[2] += h[2]) {
579:             for (jj = 0, cpoint[1] = point[1]; jj < (dim > 1 ? 2 : 1); ++jj, cpoint[1] += h[1]) {
580:               for (ii = 0, cpoint[0] = point[0]; ii < 2; ++ii, cpoint[0] += h[0]) {

582:                 DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
583:                 if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
584:               }
585:             }
586:           }
587:           /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
588:           for (edge = 0; edge < dim+1; ++edge) {
589:             PetscReal segA[6], segB[6];

591:             if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
592:             for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
593:             for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
594:               if (dim > 2) {segB[2]     = PetscRealPart(point[2]);
595:                             segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
596:               for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
597:                 if (dim > 1) {segB[1]     = PetscRealPart(point[1]);
598:                               segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
599:                 for (ii = 0; ii < 2; ++ii) {
600:                   PetscBool intersects;

602:                   segB[0]     = PetscRealPart(point[0]);
603:                   segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
604:                   DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
605:                   if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
606:                 }
607:               }
608:             }
609:           }
610:         }
611:       }
612:     }
613:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
614:   }
615:   PetscFree2(dboxes, boxes);
616:   DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
617:   DMLabelDestroy(&lbox->cellsSparse);
618:   *localBox = lbox;
619:   return(0);
620: }

622: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
623: {
624:   DM_Plex        *mesh = (DM_Plex *) dm->data;
625:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
626:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
627:   PetscInt        dim, cStart, cEnd, numCells, c, d;
628:   const PetscInt *boxCells;
629:   PetscSFNode    *cells;
630:   PetscScalar    *a;
631:   PetscMPIInt     result;
632:   PetscLogDouble  t0,t1;
633:   PetscReal       gmin[3],gmax[3];
634:   PetscInt        terminating_query_type[] = { 0, 0, 0 };
635:   PetscErrorCode  ierr;

638:   PetscLogEventBegin(DMPLEX_LocatePoints,0,0,0,0);
639:   PetscTime(&t0);
640:   if (ltype == DM_POINTLOCATION_NEAREST && !hash) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Nearest point location only supported with grid hashing. Use -dm_plex_hash_location to enable it.");
641:   DMGetCoordinateDim(dm, &dim);
642:   VecGetBlockSize(v, &bs);
643:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
644:   if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
645:   if (bs != dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Block size for point vector %D must be the mesh coordinate dimension %D", bs, dim);
646:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
647:   VecGetLocalSize(v, &numPoints);
648:   VecGetArray(v, &a);
649:   numPoints /= bs;
650:   {
651:     const PetscSFNode *sf_cells;

653:     PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
654:     if (sf_cells) {
655:       PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
656:       cells = (PetscSFNode*)sf_cells;
657:       reuse = PETSC_TRUE;
658:     } else {
659:       PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
660:       PetscMalloc1(numPoints, &cells);
661:       /* initialize cells if created */
662:       for (p=0; p<numPoints; p++) {
663:         cells[p].rank  = 0;
664:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
665:       }
666:     }
667:   }
668:   /* define domain bounding box */
669:   {
670:     Vec coorglobal;

672:     DMGetCoordinates(dm,&coorglobal);
673:     VecStrideMaxAll(coorglobal,NULL,gmax);
674:     VecStrideMinAll(coorglobal,NULL,gmin);
675:   }
676:   if (hash) {
677:     if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
678:     /* Designate the local box for each point */
679:     /* Send points to correct process */
680:     /* Search cells that lie in each subbox */
681:     /*   Should we bin points before doing search? */
682:     ISGetIndices(mesh->lbox->cells, &boxCells);
683:   }
684:   for (p = 0, numFound = 0; p < numPoints; ++p) {
685:     const PetscScalar *point = &a[p*bs];
686:     PetscInt           dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
687:     PetscBool          point_outside_domain = PETSC_FALSE;

689:     /* check bounding box of domain */
690:     for (d=0; d<dim; d++) {
691:       if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
692:       if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
693:     }
694:     if (point_outside_domain) {
695:       cells[p].rank = 0;
696:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
697:       terminating_query_type[0]++;
698:       continue;
699:     }

701:     /* check initial values in cells[].index - abort early if found */
702:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
703:       c = cells[p].index;
704:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
705:       DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
706:       if (cell >= 0) {
707:         cells[p].rank = 0;
708:         cells[p].index = cell;
709:         numFound++;
710:       }
711:     }
712:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
713:       terminating_query_type[1]++;
714:       continue;
715:     }

717:     if (hash) {
718:       PetscBool found_box;

720:       /* allow for case that point is outside box - abort early */
721:       PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
722:       if (found_box) {
723:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
724:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
725:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
726:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
727:           DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
728:           if (cell >= 0) {
729:             cells[p].rank = 0;
730:             cells[p].index = cell;
731:             numFound++;
732:             terminating_query_type[2]++;
733:             break;
734:           }
735:         }
736:       }
737:     } else {
738:       for (c = cStart; c < cEnd; ++c) {
739:         DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
740:         if (cell >= 0) {
741:           cells[p].rank = 0;
742:           cells[p].index = cell;
743:           numFound++;
744:           terminating_query_type[2]++;
745:           break;
746:         }
747:       }
748:     }
749:   }
750:   if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
751:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
752:     for (p = 0; p < numPoints; p++) {
753:       const PetscScalar *point = &a[p*bs];
754:       PetscReal          cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
755:       PetscInt           dbin[3] = {-1,-1,-1}, bin, cellOffset, d;

757:       if (cells[p].index < 0) {
758:         ++numFound;
759:         PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
760:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
761:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
762:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
763:           DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
764:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
765:           dist = DMPlex_NormD_Internal(dim, diff);
766:           if (dist < distMax) {
767:             for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
768:             cells[p].rank  = 0;
769:             cells[p].index = boxCells[c];
770:             distMax = dist;
771:             break;
772:           }
773:         }
774:       }
775:     }
776:   }
777:   /* This code is only be relevant when interfaced to parallel point location */
778:   /* Check for highest numbered proc that claims a point (do we care?) */
779:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
780:     PetscMalloc1(numFound,&found);
781:     for (p = 0, numFound = 0; p < numPoints; p++) {
782:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
783:         if (numFound < p) {
784:           cells[numFound] = cells[p];
785:         }
786:         found[numFound++] = p;
787:       }
788:     }
789:   }
790:   VecRestoreArray(v, &a);
791:   if (!reuse) {
792:     PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
793:   }
794:   PetscTime(&t1);
795:   if (hash) {
796:     PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside initial cell] : %D [hash]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
797:   } else {
798:     PetscInfo3(dm,"[DMLocatePoints_Plex] terminating_query_type : %D [outside domain] : %D [inside initial cell] : %D [brute-force]\n",terminating_query_type[0],terminating_query_type[1],terminating_query_type[2]);
799:   }
800:   PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
801:   PetscLogEventEnd(DMPLEX_LocatePoints,0,0,0,0);
802:   return(0);
803: }

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

808:   Not collective

810:   Input Parameter:
811: . coords - The coordinates of a segment

813:   Output Parameters:
814: + coords - The new y-coordinate, and 0 for x
815: - R - The rotation which accomplishes the projection

817:   Level: developer

819: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
820: @*/
821: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
822: {
823:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
824:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
825:   const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;

828:   R[0] = c; R[1] = -s;
829:   R[2] = s; R[3] =  c;
830:   coords[0] = 0.0;
831:   coords[1] = r;
832:   return(0);
833: }

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

838:   Not collective

840:   Input Parameter:
841: . coords - The coordinates of a segment

843:   Output Parameters:
844: + coords - The new y-coordinate, and 0 for x and z
845: - R - The rotation which accomplishes the projection

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

849:   Level: developer

851: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
852: @*/
853: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
854: {
855:   PetscReal      x    = PetscRealPart(coords[3] - coords[0]);
856:   PetscReal      y    = PetscRealPart(coords[4] - coords[1]);
857:   PetscReal      z    = PetscRealPart(coords[5] - coords[2]);
858:   PetscReal      r    = PetscSqrtReal(x*x + y*y + z*z);
859:   PetscReal      rinv = 1. / r;

862:   x *= rinv; y *= rinv; z *= rinv;
863:   if (x > 0.) {
864:     PetscReal inv1pX   = 1./ (1. + x);

866:     R[0] = x; R[1] = -y;              R[2] = -z;
867:     R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] =     -y*z*inv1pX;
868:     R[6] = z; R[7] =     -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
869:   }
870:   else {
871:     PetscReal inv1mX   = 1./ (1. - x);

873:     R[0] = x; R[1] = z;               R[2] = y;
874:     R[3] = y; R[4] =     -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
875:     R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] =     -y*z*inv1mX;
876:   }
877:   coords[0] = 0.0;
878:   coords[1] = r;
879:   return(0);
880: }

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

886:   Not collective

888:   Input Parameter:
889: + coordSize - Length of coordinate array (3x number of points); must be at least 9 (3 points)
890: - coords - The interlaced coordinates of each coplanar 3D point

892:   Output Parameters:
893: + coords - The first 2*coordSize/3 entries contain interlaced 2D points, with the rest undefined
894: - 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.

896:   Level: developer

898: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
899: @*/
900: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
901: {
902:   PetscReal x1[3], x2[3], n[3], c[3], norm;
903:   const PetscInt dim = 3;
904:   PetscInt       d, p;

907:   /* 0) Calculate normal vector */
908:   for (d = 0; d < dim; ++d) {
909:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
910:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
911:   }
912:   // n = x1 \otimes x2
913:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
914:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
915:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
916:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
917:   for (d = 0; d < dim; d++) n[d] /= norm;
918:   norm = PetscSqrtReal(x1[0] * x1[0] + x1[1] * x1[1] + x1[2] * x1[2]);
919:   for (d = 0; d < dim; d++) x1[d] /= norm;
920:   // x2 = n \otimes x1
921:   x2[0] = n[1] * x1[2] - n[2] * x1[1];
922:   x2[1] = n[2] * x1[0] - n[0] * x1[2];
923:   x2[2] = n[0] * x1[1] - n[1] * x1[0];
924:   for (d=0; d<dim; d++) {
925:     R[d * dim + 0] = x1[d];
926:     R[d * dim + 1] = x2[d];
927:     R[d * dim + 2] = n[d];
928:     c[d] = PetscRealPart(coords[0*dim + d]);
929:   }
930:   for (p=0; p<coordSize/dim; p++) {
931:     PetscReal y[3];
932:     for (d=0; d<dim; d++) y[d] = PetscRealPart(coords[p*dim + d]) - c[d];
933:     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];
934:   }
935:   return(0);
936: }

938: PETSC_UNUSED
939: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
940: {
941:   /* Signed volume is 1/2 the determinant

943:    |  1  1  1 |
944:    | x0 x1 x2 |
945:    | y0 y1 y2 |

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

949:    | x1 x2 |
950:    | y1 y2 |
951:   */
952:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
953:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
954:   PetscReal       M[4], detM;
955:   M[0] = x1; M[1] = x2;
956:   M[2] = y1; M[3] = y2;
957:   DMPlex_Det2D_Internal(&detM, M);
958:   *vol = 0.5*detM;
959:   (void)PetscLogFlops(5.0);
960: }

962: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
963: {
964:   DMPlex_Det2D_Internal(vol, coords);
965:   *vol *= 0.5;
966: }

968: PETSC_UNUSED
969: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
970: {
971:   /* Signed volume is 1/6th of the determinant

973:    |  1  1  1  1 |
974:    | x0 x1 x2 x3 |
975:    | y0 y1 y2 y3 |
976:    | z0 z1 z2 z3 |

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

980:    | x1 x2 x3 |
981:    | y1 y2 y3 |
982:    | z1 z2 z3 |
983:   */
984:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
985:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
986:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
987:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
988:   PetscReal       M[9], detM;
989:   M[0] = x1; M[1] = x2; M[2] = x3;
990:   M[3] = y1; M[4] = y2; M[5] = y3;
991:   M[6] = z1; M[7] = z2; M[8] = z3;
992:   DMPlex_Det3D_Internal(&detM, M);
993:   *vol = -onesixth*detM;
994:   (void)PetscLogFlops(10.0);
995: }

997: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
998: {
999:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1000:   DMPlex_Det3D_Internal(vol, coords);
1001:   *vol *= -onesixth;
1002: }

1004: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1005: {
1006:   PetscSection   coordSection;
1007:   Vec            coordinates;
1008:   const PetscScalar *coords;
1009:   PetscInt       dim, d, off;

1013:   DMGetCoordinatesLocal(dm, &coordinates);
1014:   DMGetCoordinateSection(dm, &coordSection);
1015:   PetscSectionGetDof(coordSection,e,&dim);
1016:   if (!dim) return(0);
1017:   PetscSectionGetOffset(coordSection,e,&off);
1018:   VecGetArrayRead(coordinates,&coords);
1019:   if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1020:   VecRestoreArrayRead(coordinates,&coords);
1021:   *detJ = 1.;
1022:   if (J) {
1023:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1024:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1025:     if (invJ) {
1026:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1027:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1028:     }
1029:   }
1030:   return(0);
1031: }

1033: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1034: {
1035:   PetscSection   coordSection;
1036:   Vec            coordinates;
1037:   PetscScalar   *coords = NULL;
1038:   PetscInt       numCoords, d, pStart, pEnd, numSelfCoords = 0;

1042:   DMGetCoordinatesLocal(dm, &coordinates);
1043:   DMGetCoordinateSection(dm, &coordSection);
1044:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1045:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1046:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1047:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1048:   if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1049:   *detJ = 0.0;
1050:   if (numCoords == 6) {
1051:     const PetscInt dim = 3;
1052:     PetscReal      R[9], J0;

1054:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1055:     DMPlexComputeProjection3Dto1D(coords, R);
1056:     if (J)    {
1057:       J0   = 0.5*PetscRealPart(coords[1]);
1058:       J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1059:       J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1060:       J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1061:       DMPlex_Det3D_Internal(detJ, J);
1062:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1063:     }
1064:   } else if (numCoords == 4) {
1065:     const PetscInt dim = 2;
1066:     PetscReal      R[4], J0;

1068:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1069:     DMPlexComputeProjection2Dto1D(coords, R);
1070:     if (J)    {
1071:       J0   = 0.5*PetscRealPart(coords[1]);
1072:       J[0] = R[0]*J0; J[1] = R[1];
1073:       J[2] = R[2]*J0; J[3] = R[3];
1074:       DMPlex_Det2D_Internal(detJ, J);
1075:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1076:     }
1077:   } else if (numCoords == 2) {
1078:     const PetscInt dim = 1;

1080:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1081:     if (J)    {
1082:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1083:       *detJ = J[0];
1084:       PetscLogFlops(2.0);
1085:       if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1086:     }
1087:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1088:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1089:   return(0);
1090: }

1092: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1093: {
1094:   PetscSection   coordSection;
1095:   Vec            coordinates;
1096:   PetscScalar   *coords = NULL;
1097:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1101:   DMGetCoordinatesLocal(dm, &coordinates);
1102:   DMGetCoordinateSection(dm, &coordSection);
1103:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1104:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1105:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1106:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1107:   *detJ = 0.0;
1108:   if (numCoords == 9) {
1109:     const PetscInt dim = 3;
1110:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1112:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1113:     DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1114:     if (J)    {
1115:       const PetscInt pdim = 2;

1117:       for (d = 0; d < pdim; d++) {
1118:         for (f = 0; f < pdim; f++) {
1119:           J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1120:         }
1121:       }
1122:       PetscLogFlops(8.0);
1123:       DMPlex_Det3D_Internal(detJ, J0);
1124:       for (d = 0; d < dim; d++) {
1125:         for (f = 0; f < dim; f++) {
1126:           J[d*dim+f] = 0.0;
1127:           for (g = 0; g < dim; g++) {
1128:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1129:           }
1130:         }
1131:       }
1132:       PetscLogFlops(18.0);
1133:     }
1134:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1135:   } else if (numCoords == 6) {
1136:     const PetscInt dim = 2;

1138:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1139:     if (J)    {
1140:       for (d = 0; d < dim; d++) {
1141:         for (f = 0; f < dim; f++) {
1142:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1143:         }
1144:       }
1145:       PetscLogFlops(8.0);
1146:       DMPlex_Det2D_Internal(detJ, J);
1147:     }
1148:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1149:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1150:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1151:   return(0);
1152: }

1154: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscBool isTensor, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1155: {
1156:   PetscSection   coordSection;
1157:   Vec            coordinates;
1158:   PetscScalar   *coords = NULL;
1159:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1163:   DMGetCoordinatesLocal(dm, &coordinates);
1164:   DMGetCoordinateSection(dm, &coordSection);
1165:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1166:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1167:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1168:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1169:   if (!Nq) {
1170:     PetscInt vorder[4] = {0, 1, 2, 3};

1172:     if (isTensor) {vorder[2] = 3; vorder[3] = 2;}
1173:     *detJ = 0.0;
1174:     if (numCoords == 12) {
1175:       const PetscInt dim = 3;
1176:       PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1178:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1179:       DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1180:       if (J)    {
1181:         const PetscInt pdim = 2;

1183:         for (d = 0; d < pdim; d++) {
1184:           J0[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*pdim+d]) - PetscRealPart(coords[vorder[0]*pdim+d]));
1185:           J0[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[2]*pdim+d]) - PetscRealPart(coords[vorder[1]*pdim+d]));
1186:         }
1187:         PetscLogFlops(8.0);
1188:         DMPlex_Det3D_Internal(detJ, J0);
1189:         for (d = 0; d < dim; d++) {
1190:           for (f = 0; f < dim; f++) {
1191:             J[d*dim+f] = 0.0;
1192:             for (g = 0; g < dim; g++) {
1193:               J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1194:             }
1195:           }
1196:         }
1197:         PetscLogFlops(18.0);
1198:       }
1199:       if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1200:     } else if (numCoords == 8) {
1201:       const PetscInt dim = 2;

1203:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1204:       if (J)    {
1205:         for (d = 0; d < dim; d++) {
1206:           J[d*dim+0] = 0.5*(PetscRealPart(coords[vorder[1]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1207:           J[d*dim+1] = 0.5*(PetscRealPart(coords[vorder[3]*dim+d]) - PetscRealPart(coords[vorder[0]*dim+d]));
1208:         }
1209:         PetscLogFlops(8.0);
1210:         DMPlex_Det2D_Internal(detJ, J);
1211:       }
1212:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1213:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1214:   } else {
1215:     const PetscInt Nv = 4;
1216:     const PetscInt dimR = 2;
1217:     PetscInt  zToPlex[4] = {0, 1, 3, 2};
1218:     PetscReal zOrder[12];
1219:     PetscReal zCoeff[12];
1220:     PetscInt  i, j, k, l, dim;

1222:     if (isTensor) {zToPlex[2] = 2; zToPlex[3] = 3;}
1223:     if (numCoords == 12) {
1224:       dim = 3;
1225:     } else if (numCoords == 8) {
1226:       dim = 2;
1227:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1228:     for (i = 0; i < Nv; i++) {
1229:       PetscInt zi = zToPlex[i];

1231:       for (j = 0; j < dim; j++) {
1232:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1233:       }
1234:     }
1235:     for (j = 0; j < dim; j++) {
1236:       zCoeff[dim * 0 + j] = 0.25 * (  zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1237:       zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1238:       zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1239:       zCoeff[dim * 3 + j] = 0.25 * (  zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1240:     }
1241:     for (i = 0; i < Nq; i++) {
1242:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1244:       if (v) {
1245:         PetscReal extPoint[4];

1247:         extPoint[0] = 1.;
1248:         extPoint[1] = xi;
1249:         extPoint[2] = eta;
1250:         extPoint[3] = xi * eta;
1251:         for (j = 0; j < dim; j++) {
1252:           PetscReal val = 0.;

1254:           for (k = 0; k < Nv; k++) {
1255:             val += extPoint[k] * zCoeff[dim * k + j];
1256:           }
1257:           v[i * dim + j] = val;
1258:         }
1259:       }
1260:       if (J) {
1261:         PetscReal extJ[8];

1263:         extJ[0] = 0.;
1264:         extJ[1] = 0.;
1265:         extJ[2] = 1.;
1266:         extJ[3] = 0.;
1267:         extJ[4] = 0.;
1268:         extJ[5] = 1.;
1269:         extJ[6] = eta;
1270:         extJ[7] = xi;
1271:         for (j = 0; j < dim; j++) {
1272:           for (k = 0; k < dimR; k++) {
1273:             PetscReal val = 0.;

1275:             for (l = 0; l < Nv; l++) {
1276:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1277:             }
1278:             J[i * dim * dim + dim * j + k] = val;
1279:           }
1280:         }
1281:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1282:           PetscReal x, y, z;
1283:           PetscReal *iJ = &J[i * dim * dim];
1284:           PetscReal norm;

1286:           x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1287:           y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1288:           z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1289:           norm = PetscSqrtReal(x * x + y * y + z * z);
1290:           iJ[2] = x / norm;
1291:           iJ[5] = y / norm;
1292:           iJ[8] = z / norm;
1293:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1294:           if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1295:         } else {
1296:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1297:           if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1298:         }
1299:       }
1300:     }
1301:   }
1302:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1303:   return(0);
1304: }

1306: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1307: {
1308:   PetscSection   coordSection;
1309:   Vec            coordinates;
1310:   PetscScalar   *coords = NULL;
1311:   const PetscInt dim = 3;
1312:   PetscInt       d;

1316:   DMGetCoordinatesLocal(dm, &coordinates);
1317:   DMGetCoordinateSection(dm, &coordSection);
1318:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1319:   *detJ = 0.0;
1320:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1321:   if (J)    {
1322:     for (d = 0; d < dim; d++) {
1323:       /* I orient with outward face normals */
1324:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1325:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1326:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1327:     }
1328:     PetscLogFlops(18.0);
1329:     DMPlex_Det3D_Internal(detJ, J);
1330:   }
1331:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1332:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1333:   return(0);
1334: }

1336: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1337: {
1338:   PetscSection   coordSection;
1339:   Vec            coordinates;
1340:   PetscScalar   *coords = NULL;
1341:   const PetscInt dim = 3;
1342:   PetscInt       d;

1346:   DMGetCoordinatesLocal(dm, &coordinates);
1347:   DMGetCoordinateSection(dm, &coordSection);
1348:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1349:   if (!Nq) {
1350:     *detJ = 0.0;
1351:     if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1352:     if (J)    {
1353:       for (d = 0; d < dim; d++) {
1354:         J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1355:         J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1356:         J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1357:       }
1358:       PetscLogFlops(18.0);
1359:       DMPlex_Det3D_Internal(detJ, J);
1360:     }
1361:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1362:   } else {
1363:     const PetscInt Nv = 8;
1364:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1365:     const PetscInt dim = 3;
1366:     const PetscInt dimR = 3;
1367:     PetscReal zOrder[24];
1368:     PetscReal zCoeff[24];
1369:     PetscInt  i, j, k, l;

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

1374:       for (j = 0; j < dim; j++) {
1375:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1376:       }
1377:     }
1378:     for (j = 0; j < dim; j++) {
1379:       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]);
1380:       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]);
1381:       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]);
1382:       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]);
1383:       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]);
1384:       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]);
1385:       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]);
1386:       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]);
1387:     }
1388:     for (i = 0; i < Nq; i++) {
1389:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

1391:       if (v) {
1392:         PetscReal extPoint[8];

1394:         extPoint[0] = 1.;
1395:         extPoint[1] = xi;
1396:         extPoint[2] = eta;
1397:         extPoint[3] = xi * eta;
1398:         extPoint[4] = theta;
1399:         extPoint[5] = theta * xi;
1400:         extPoint[6] = theta * eta;
1401:         extPoint[7] = theta * eta * xi;
1402:         for (j = 0; j < dim; j++) {
1403:           PetscReal val = 0.;

1405:           for (k = 0; k < Nv; k++) {
1406:             val += extPoint[k] * zCoeff[dim * k + j];
1407:           }
1408:           v[i * dim + j] = val;
1409:         }
1410:       }
1411:       if (J) {
1412:         PetscReal extJ[24];

1414:         extJ[0]  = 0.         ; extJ[1]  = 0.        ; extJ[2]  = 0.      ;
1415:         extJ[3]  = 1.         ; extJ[4]  = 0.        ; extJ[5]  = 0.      ;
1416:         extJ[6]  = 0.         ; extJ[7]  = 1.        ; extJ[8]  = 0.      ;
1417:         extJ[9]  = eta        ; extJ[10] = xi        ; extJ[11] = 0.      ;
1418:         extJ[12] = 0.         ; extJ[13] = 0.        ; extJ[14] = 1.      ;
1419:         extJ[15] = theta      ; extJ[16] = 0.        ; extJ[17] = xi      ;
1420:         extJ[18] = 0.         ; extJ[19] = theta     ; extJ[20] = eta     ;
1421:         extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;

1423:         for (j = 0; j < dim; j++) {
1424:           for (k = 0; k < dimR; k++) {
1425:             PetscReal val = 0.;

1427:             for (l = 0; l < Nv; l++) {
1428:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1429:             }
1430:             J[i * dim * dim + dim * j + k] = val;
1431:           }
1432:         }
1433:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1434:         if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1435:       }
1436:     }
1437:   }
1438:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1439:   return(0);
1440: }

1442: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1443: {
1444:   DMPolytopeType  ct;
1445:   PetscInt        depth, dim, coordDim, coneSize, i;
1446:   PetscInt        Nq = 0;
1447:   const PetscReal *points = NULL;
1448:   DMLabel         depthLabel;
1449:   PetscReal       xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1450:   PetscBool       isAffine = PETSC_TRUE;
1451:   PetscErrorCode  ierr;

1454:   DMPlexGetDepth(dm, &depth);
1455:   DMPlexGetConeSize(dm, cell, &coneSize);
1456:   DMPlexGetDepthLabel(dm, &depthLabel);
1457:   DMLabelGetValue(depthLabel, cell, &dim);
1458:   if (depth == 1 && dim == 1) {
1459:     DMGetDimension(dm, &dim);
1460:   }
1461:   DMGetCoordinateDim(dm, &coordDim);
1462:   if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1463:   if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1464:   DMPlexGetCellType(dm, cell, &ct);
1465:   switch (ct) {
1466:     case DM_POLYTOPE_POINT:
1467:     DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1468:     isAffine = PETSC_FALSE;
1469:     break;
1470:     case DM_POLYTOPE_SEGMENT:
1471:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1472:     if (Nq) {DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1473:     else    {DMPlexComputeLineGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1474:     break;
1475:     case DM_POLYTOPE_TRIANGLE:
1476:     if (Nq) {DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1477:     else    {DMPlexComputeTriangleGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1478:     break;
1479:     case DM_POLYTOPE_QUADRILATERAL:
1480:     DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_FALSE, Nq, points, v, J, invJ, detJ);
1481:     isAffine = PETSC_FALSE;
1482:     break;
1483:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1484:     DMPlexComputeRectangleGeometry_Internal(dm, cell, PETSC_TRUE, Nq, points, v, J, invJ, detJ);
1485:     isAffine = PETSC_FALSE;
1486:     break;
1487:     case DM_POLYTOPE_TETRAHEDRON:
1488:     if (Nq) {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1489:     else    {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1490:     break;
1491:     case DM_POLYTOPE_HEXAHEDRON:
1492:     DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1493:     isAffine = PETSC_FALSE;
1494:     break;
1495:     default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]);
1496:   }
1497:   if (isAffine && Nq) {
1498:     if (v) {
1499:       for (i = 0; i < Nq; i++) {
1500:         CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1501:       }
1502:     }
1503:     if (detJ) {
1504:       for (i = 0; i < Nq; i++) {
1505:         detJ[i] = detJ0;
1506:       }
1507:     }
1508:     if (J) {
1509:       PetscInt k;

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

1514:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1515:           J[k] = J0[j];
1516:         }
1517:       }
1518:     }
1519:     if (invJ) {
1520:       PetscInt k;
1521:       switch (coordDim) {
1522:       case 0:
1523:         break;
1524:       case 1:
1525:         invJ[0] = 1./J0[0];
1526:         break;
1527:       case 2:
1528:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1529:         break;
1530:       case 3:
1531:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1532:         break;
1533:       }
1534:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1535:         PetscInt j;

1537:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1538:           invJ[k] = invJ[j];
1539:         }
1540:       }
1541:     }
1542:   }
1543:   return(0);
1544: }

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

1549:   Collective on dm

1551:   Input Arguments:
1552: + dm   - the DM
1553: - cell - the cell

1555:   Output Arguments:
1556: + v0   - the translation part of this affine transform
1557: . J    - the Jacobian of the transform from the reference element
1558: . invJ - the inverse of the Jacobian
1559: - detJ - the Jacobian determinant

1561:   Level: advanced

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

1567: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1568: @*/
1569: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1570: {

1574:   DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1575:   return(0);
1576: }

1578: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1579: {
1580:   PetscQuadrature   feQuad;
1581:   PetscSection      coordSection;
1582:   Vec               coordinates;
1583:   PetscScalar      *coords = NULL;
1584:   const PetscReal  *quadPoints;
1585:   PetscTabulation T;
1586:   PetscInt          dim, cdim, pdim, qdim, Nq, numCoords, q;
1587:   PetscErrorCode    ierr;

1590:   DMGetCoordinatesLocal(dm, &coordinates);
1591:   DMGetCoordinateSection(dm, &coordSection);
1592:   DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1593:   DMGetDimension(dm, &dim);
1594:   DMGetCoordinateDim(dm, &cdim);
1595:   if (!quad) { /* use the first point of the first functional of the dual space */
1596:     PetscDualSpace dsp;

1598:     PetscFEGetDualSpace(fe, &dsp);
1599:     PetscDualSpaceGetFunctional(dsp, 0, &quad);
1600:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1601:     Nq = 1;
1602:   } else {
1603:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1604:   }
1605:   PetscFEGetDimension(fe, &pdim);
1606:   PetscFEGetQuadrature(fe, &feQuad);
1607:   if (feQuad == quad) {
1608:     PetscFEGetCellTabulation(fe, &T);
1609:     if (numCoords != pdim*cdim) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "There are %d coordinates for point %d != %d*%d", numCoords, point, pdim, cdim);
1610:   } else {
1611:     PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
1612:   }
1613:   if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1614:   {
1615:     const PetscReal *basis    = T->T[0];
1616:     const PetscReal *basisDer = T->T[1];
1617:     PetscReal        detJt;

1619:     if (v) {
1620:       PetscArrayzero(v, Nq*cdim);
1621:       for (q = 0; q < Nq; ++q) {
1622:         PetscInt i, k;

1624:         for (k = 0; k < pdim; ++k) {
1625:           const PetscInt vertex = k/cdim;
1626:           for (i = 0; i < cdim; ++i) {
1627:             v[q*cdim + i] += basis[(q*pdim + k)*cdim + i] * PetscRealPart(coords[vertex*cdim + i]);
1628:           }
1629:         }
1630:         PetscLogFlops(2.0*pdim*cdim);
1631:       }
1632:     }
1633:     if (J) {
1634:       PetscArrayzero(J, Nq*cdim*cdim);
1635:       for (q = 0; q < Nq; ++q) {
1636:         PetscInt i, j, k, c, r;

1638:         /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1639:         for (k = 0; k < pdim; ++k) {
1640:           const PetscInt vertex = k/cdim;
1641:           for (j = 0; j < dim; ++j) {
1642:             for (i = 0; i < cdim; ++i) {
1643:               J[(q*cdim + i)*cdim + j] += basisDer[((q*pdim + k)*cdim + i)*dim + j] * PetscRealPart(coords[vertex*cdim + i]);
1644:             }
1645:           }
1646:         }
1647:         PetscLogFlops(2.0*pdim*dim*cdim);
1648:         if (cdim > dim) {
1649:           for (c = dim; c < cdim; ++c)
1650:             for (r = 0; r < cdim; ++r)
1651:               J[r*cdim+c] = r == c ? 1.0 : 0.0;
1652:         }
1653:         if (!detJ && !invJ) continue;
1654:         detJt = 0.;
1655:         switch (cdim) {
1656:         case 3:
1657:           DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1658:           if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1659:           break;
1660:         case 2:
1661:           DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1662:           if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1663:           break;
1664:         case 1:
1665:           detJt = J[q*cdim*dim];
1666:           if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1667:         }
1668:         if (detJ) detJ[q] = detJt;
1669:       }
1670:     } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1671:   }
1672:   if (feQuad != quad) {PetscTabulationDestroy(&T);}
1673:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1674:   return(0);
1675: }

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

1680:   Collective on dm

1682:   Input Arguments:
1683: + dm   - the DM
1684: . cell - the cell
1685: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If quad == NULL, geometry will be
1686:          evaluated at the first vertex of the reference element

1688:   Output Arguments:
1689: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1690: . J    - the Jacobian of the transform from the reference element at each quadrature point
1691: . invJ - the inverse of the Jacobian at each quadrature point
1692: - detJ - the Jacobian determinant at each quadrature point

1694:   Level: advanced

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

1700: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1701: @*/
1702: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1703: {
1704:   DM             cdm;
1705:   PetscFE        fe = NULL;

1710:   DMGetCoordinateDM(dm, &cdm);
1711:   if (cdm) {
1712:     PetscClassId id;
1713:     PetscInt     numFields;
1714:     PetscDS      prob;
1715:     PetscObject  disc;

1717:     DMGetNumFields(cdm, &numFields);
1718:     if (numFields) {
1719:       DMGetDS(cdm, &prob);
1720:       PetscDSGetDiscretization(prob,0,&disc);
1721:       PetscObjectGetClassId(disc,&id);
1722:       if (id == PETSCFE_CLASSID) {
1723:         fe = (PetscFE) disc;
1724:       }
1725:     }
1726:   }
1727:   if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1728:   else     {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1729:   return(0);
1730: }

1732: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1733: {
1734:   PetscSection   coordSection;
1735:   Vec            coordinates;
1736:   PetscScalar   *coords = NULL;
1737:   PetscScalar    tmp[2];
1738:   PetscInt       coordSize, d;

1742:   DMGetCoordinatesLocal(dm, &coordinates);
1743:   DMGetCoordinateSection(dm, &coordSection);
1744:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1745:   DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1746:   if (centroid) {
1747:     for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1748:   }
1749:   if (normal) {
1750:     PetscReal norm;

1752:     if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1753:     normal[0]  = -PetscRealPart(coords[1] - tmp[1]);
1754:     normal[1]  =  PetscRealPart(coords[0] - tmp[0]);
1755:     norm       = DMPlex_NormD_Internal(dim, normal);
1756:     for (d = 0; d < dim; ++d) normal[d] /= norm;
1757:   }
1758:   if (vol) {
1759:     *vol = 0.0;
1760:     for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1761:     *vol = PetscSqrtReal(*vol);
1762:   }
1763:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1764:   return(0);
1765: }

1767: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1768: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1769: {
1770:   DMPolytopeType ct;
1771:   PetscSection   coordSection;
1772:   Vec            coordinates;
1773:   PetscScalar   *coords = NULL;
1774:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1775:   PetscBool      isHybrid = PETSC_FALSE;
1776:   PetscInt       fv[4] = {0, 1, 2, 3};
1777:   PetscInt       tdim = 2, coordSize, numCorners, p, d, e;

1781:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1782:   DMPlexGetCellType(dm, cell, &ct);
1783:   switch (ct) {
1784:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1785:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1786:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1787:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1788:       isHybrid = PETSC_TRUE;
1789:     default: break;
1790:   }
1791:   DMGetCoordinatesLocal(dm, &coordinates);
1792:   DMPlexGetConeSize(dm, cell, &numCorners);
1793:   DMGetCoordinateSection(dm, &coordSection);
1794:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1795:   DMGetCoordinateDim(dm, &dim);
1796:   /* Side faces for hybrid cells are are stored as tensor products */
1797:   if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}

1799:   if (dim > 2 && centroid) {
1800:     v0[0] = PetscRealPart(coords[0]);
1801:     v0[1] = PetscRealPart(coords[1]);
1802:     v0[2] = PetscRealPart(coords[2]);
1803:   }
1804:   if (normal) {
1805:     if (dim > 2) {
1806:       const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1807:       const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1808:       const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1809:       PetscReal       norm;

1811:       normal[0] = y0*z1 - z0*y1;
1812:       normal[1] = z0*x1 - x0*z1;
1813:       normal[2] = x0*y1 - y0*x1;
1814:       norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1815:       normal[0] /= norm;
1816:       normal[1] /= norm;
1817:       normal[2] /= norm;
1818:     } else {
1819:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
1820:     }
1821:   }
1822:   if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1823:   for (p = 0; p < numCorners; ++p) {
1824:     const PetscInt pi  = p < 4 ? fv[p] : p;
1825:     const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1826:     /* Need to do this copy to get types right */
1827:     for (d = 0; d < tdim; ++d) {
1828:       ctmp[d]      = PetscRealPart(coords[pi*tdim+d]);
1829:       ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1830:     }
1831:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1832:     vsum += vtmp;
1833:     for (d = 0; d < tdim; ++d) {
1834:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1835:     }
1836:   }
1837:   for (d = 0; d < tdim; ++d) {
1838:     csum[d] /= (tdim+1)*vsum;
1839:   }
1840:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1841:   if (vol) *vol = PetscAbsReal(vsum);
1842:   if (centroid) {
1843:     if (dim > 2) {
1844:       for (d = 0; d < dim; ++d) {
1845:         centroid[d] = v0[d];
1846:         for (e = 0; e < dim; ++e) {
1847:           centroid[d] += R[d*dim+e]*csum[e];
1848:         }
1849:       }
1850:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1851:   }
1852:   return(0);
1853: }

1855: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1856: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1857: {
1858:   DMPolytopeType  ct;
1859:   PetscSection    coordSection;
1860:   Vec             coordinates;
1861:   PetscScalar    *coords = NULL;
1862:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
1863:   const PetscInt *faces, *facesO;
1864:   PetscBool       isHybrid = PETSC_FALSE;
1865:   PetscInt        numFaces, f, coordSize, p, d;
1866:   PetscErrorCode  ierr;

1869:   if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1870:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1871:   DMPlexGetCellType(dm, cell, &ct);
1872:   switch (ct) {
1873:     case DM_POLYTOPE_POINT_PRISM_TENSOR:
1874:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1875:     case DM_POLYTOPE_TRI_PRISM_TENSOR:
1876:     case DM_POLYTOPE_QUAD_PRISM_TENSOR:
1877:       isHybrid = PETSC_TRUE;
1878:     default: break;
1879:   }

1881:   DMGetCoordinatesLocal(dm, &coordinates);
1882:   DMGetCoordinateSection(dm, &coordSection);

1884:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1885:   DMPlexGetConeSize(dm, cell, &numFaces);
1886:   DMPlexGetCone(dm, cell, &faces);
1887:   DMPlexGetConeOrientation(dm, cell, &facesO);
1888:   for (f = 0; f < numFaces; ++f) {
1889:     PetscBool      flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
1890:     DMPolytopeType ct;

1892:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1893:     DMPlexGetCellType(dm, faces[f], &ct);
1894:     switch (ct) {
1895:     case DM_POLYTOPE_TRIANGLE:
1896:       for (d = 0; d < dim; ++d) {
1897:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1898:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1899:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1900:       }
1901:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1902:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1903:       vsum += vtmp;
1904:       if (centroid) {           /* Centroid of OABC = (a+b+c)/4 */
1905:         for (d = 0; d < dim; ++d) {
1906:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1907:         }
1908:       }
1909:       break;
1910:     case DM_POLYTOPE_QUADRILATERAL:
1911:     case DM_POLYTOPE_SEG_PRISM_TENSOR:
1912:     {
1913:       PetscInt fv[4] = {0, 1, 2, 3};

1915:       /* Side faces for hybrid cells are are stored as tensor products */
1916:       if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1917:       /* DO FOR PYRAMID */
1918:       /* First tet */
1919:       for (d = 0; d < dim; ++d) {
1920:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1921:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1922:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1923:       }
1924:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1925:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1926:       vsum += vtmp;
1927:       if (centroid) {
1928:         for (d = 0; d < dim; ++d) {
1929:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1930:         }
1931:       }
1932:       /* Second tet */
1933:       for (d = 0; d < dim; ++d) {
1934:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1935:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1936:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1937:       }
1938:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1939:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1940:       vsum += vtmp;
1941:       if (centroid) {
1942:         for (d = 0; d < dim; ++d) {
1943:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1944:         }
1945:       }
1946:       break;
1947:     }
1948:     default:
1949:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]);
1950:     }
1951:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1952:   }
1953:   if (vol)     *vol = PetscAbsReal(vsum);
1954:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
1955:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1956:   return(0);
1957: }

1959: /*@C
1960:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

1962:   Collective on dm

1964:   Input Arguments:
1965: + dm   - the DM
1966: - cell - the cell

1968:   Output Arguments:
1969: + volume   - the cell volume
1970: . centroid - the cell centroid
1971: - normal - the cell normal, if appropriate

1973:   Level: advanced

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

1979: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1980: @*/
1981: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1982: {
1983:   PetscInt       depth, dim;

1987:   DMPlexGetDepth(dm, &depth);
1988:   DMGetDimension(dm, &dim);
1989:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1990:   DMPlexGetPointDepth(dm, cell, &depth);
1991:   switch (depth) {
1992:   case 1:
1993:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1994:     break;
1995:   case 2:
1996:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1997:     break;
1998:   case 3:
1999:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2000:     break;
2001:   default:
2002:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2003:   }
2004:   return(0);
2005: }

2007: /*@
2008:   DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh

2010:   Collective on dm

2012:   Input Parameter:
2013: . dm - The DMPlex

2015:   Output Parameter:
2016: . cellgeom - A vector with the cell geometry data for each cell

2018:   Level: beginner

2020: @*/
2021: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2022: {
2023:   DM             dmCell;
2024:   Vec            coordinates;
2025:   PetscSection   coordSection, sectionCell;
2026:   PetscScalar   *cgeom;
2027:   PetscInt       cStart, cEnd, c;

2031:   DMClone(dm, &dmCell);
2032:   DMGetCoordinateSection(dm, &coordSection);
2033:   DMGetCoordinatesLocal(dm, &coordinates);
2034:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2035:   DMSetCoordinatesLocal(dmCell, coordinates);
2036:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2037:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2038:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2039:   /* TODO This needs to be multiplied by Nq for non-affine */
2040:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2041:   PetscSectionSetUp(sectionCell);
2042:   DMSetLocalSection(dmCell, sectionCell);
2043:   PetscSectionDestroy(&sectionCell);
2044:   DMCreateLocalVector(dmCell, cellgeom);
2045:   VecGetArray(*cellgeom, &cgeom);
2046:   for (c = cStart; c < cEnd; ++c) {
2047:     PetscFEGeom *cg;

2049:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2050:     PetscArrayzero(cg, 1);
2051:     DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2052:     if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2053:   }
2054:   return(0);
2055: }

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

2060:   Input Parameter:
2061: . dm - The DM

2063:   Output Parameters:
2064: + cellgeom - A Vec of PetscFVCellGeom data
2065: - facegeom - A Vec of PetscFVFaceGeom data

2067:   Level: developer

2069: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2070: @*/
2071: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2072: {
2073:   DM             dmFace, dmCell;
2074:   DMLabel        ghostLabel;
2075:   PetscSection   sectionFace, sectionCell;
2076:   PetscSection   coordSection;
2077:   Vec            coordinates;
2078:   PetscScalar   *fgeom, *cgeom;
2079:   PetscReal      minradius, gminradius;
2080:   PetscInt       dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

2084:   DMGetDimension(dm, &dim);
2085:   DMGetCoordinateSection(dm, &coordSection);
2086:   DMGetCoordinatesLocal(dm, &coordinates);
2087:   /* Make cell centroids and volumes */
2088:   DMClone(dm, &dmCell);
2089:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2090:   DMSetCoordinatesLocal(dmCell, coordinates);
2091:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2092:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2093:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2094:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2095:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2096:   PetscSectionSetUp(sectionCell);
2097:   DMSetLocalSection(dmCell, sectionCell);
2098:   PetscSectionDestroy(&sectionCell);
2099:   DMCreateLocalVector(dmCell, cellgeom);
2100:   if (cEndInterior < 0) cEndInterior = cEnd;
2101:   VecGetArray(*cellgeom, &cgeom);
2102:   for (c = cStart; c < cEndInterior; ++c) {
2103:     PetscFVCellGeom *cg;

2105:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2106:     PetscArrayzero(cg, 1);
2107:     DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2108:   }
2109:   /* Compute face normals and minimum cell radius */
2110:   DMClone(dm, &dmFace);
2111:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionFace);
2112:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2113:   PetscSectionSetChart(sectionFace, fStart, fEnd);
2114:   for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2115:   PetscSectionSetUp(sectionFace);
2116:   DMSetLocalSection(dmFace, sectionFace);
2117:   PetscSectionDestroy(&sectionFace);
2118:   DMCreateLocalVector(dmFace, facegeom);
2119:   VecGetArray(*facegeom, &fgeom);
2120:   DMGetLabel(dm, "ghost", &ghostLabel);
2121:   minradius = PETSC_MAX_REAL;
2122:   for (f = fStart; f < fEnd; ++f) {
2123:     PetscFVFaceGeom *fg;
2124:     PetscReal        area;
2125:     const PetscInt  *cells;
2126:     PetscInt         ncells, ghost = -1, d, numChildren;

2128:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2129:     DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2130:     DMPlexGetSupport(dm, f, &cells);
2131:     DMPlexGetSupportSize(dm, f, &ncells);
2132:     /* It is possible to get a face with no support when using partition overlap */
2133:     if (!ncells || ghost >= 0 || numChildren) continue;
2134:     DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2135:     DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2136:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2137:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2138:     {
2139:       PetscFVCellGeom *cL, *cR;
2140:       PetscReal       *lcentroid, *rcentroid;
2141:       PetscReal        l[3], r[3], v[3];

2143:       DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2144:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2145:       if (ncells > 1) {
2146:         DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2147:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2148:       }
2149:       else {
2150:         rcentroid = fg->centroid;
2151:       }
2152:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2153:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2154:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2155:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2156:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2157:       }
2158:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2159:         if (dim == 2) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g) v (%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) v[0], (double) v[1]);
2160:         if (dim == 3) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed, normal (%g,%g,%g) v (%g,%g,%g)", f, (double) fg->normal[0], (double) fg->normal[1], (double) fg->normal[2], (double) v[0], (double) v[1], (double) v[2]);
2161:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2162:       }
2163:       if (cells[0] < cEndInterior) {
2164:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2165:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2166:       }
2167:       if (ncells > 1 && cells[1] < cEndInterior) {
2168:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2169:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2170:       }
2171:     }
2172:   }
2173:   MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2174:   DMPlexSetMinRadius(dm, gminradius);
2175:   /* Compute centroids of ghost cells */
2176:   for (c = cEndInterior; c < cEnd; ++c) {
2177:     PetscFVFaceGeom *fg;
2178:     const PetscInt  *cone,    *support;
2179:     PetscInt         coneSize, supportSize, s;

2181:     DMPlexGetConeSize(dmCell, c, &coneSize);
2182:     if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2183:     DMPlexGetCone(dmCell, c, &cone);
2184:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2185:     if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2186:     DMPlexGetSupport(dmCell, cone[0], &support);
2187:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2188:     for (s = 0; s < 2; ++s) {
2189:       /* Reflect ghost centroid across plane of face */
2190:       if (support[s] == c) {
2191:         PetscFVCellGeom       *ci;
2192:         PetscFVCellGeom       *cg;
2193:         PetscReal              c2f[3], a;

2195:         DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2196:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2197:         a    = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2198:         DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2199:         DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2200:         cg->volume = ci->volume;
2201:       }
2202:     }
2203:   }
2204:   VecRestoreArray(*facegeom, &fgeom);
2205:   VecRestoreArray(*cellgeom, &cgeom);
2206:   DMDestroy(&dmCell);
2207:   DMDestroy(&dmFace);
2208:   return(0);
2209: }

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

2214:   Not collective

2216:   Input Argument:
2217: . dm - the DM

2219:   Output Argument:
2220: . minradius - the minium cell radius

2222:   Level: developer

2224: .seealso: DMGetCoordinates()
2225: @*/
2226: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2227: {
2231:   *minradius = ((DM_Plex*) dm->data)->minradius;
2232:   return(0);
2233: }

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

2238:   Logically collective

2240:   Input Arguments:
2241: + dm - the DM
2242: - minradius - the minium cell radius

2244:   Level: developer

2246: .seealso: DMSetCoordinates()
2247: @*/
2248: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2249: {
2252:   ((DM_Plex*) dm->data)->minradius = minradius;
2253:   return(0);
2254: }

2256: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2257: {
2258:   DMLabel        ghostLabel;
2259:   PetscScalar   *dx, *grad, **gref;
2260:   PetscInt       dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

2264:   DMGetDimension(dm, &dim);
2265:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2266:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2267:   DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2268:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2269:   DMGetLabel(dm, "ghost", &ghostLabel);
2270:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2271:   for (c = cStart; c < cEndInterior; c++) {
2272:     const PetscInt        *faces;
2273:     PetscInt               numFaces, usedFaces, f, d;
2274:     PetscFVCellGeom        *cg;
2275:     PetscBool              boundary;
2276:     PetscInt               ghost;

2278:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2279:     DMPlexGetConeSize(dm, c, &numFaces);
2280:     DMPlexGetCone(dm, c, &faces);
2281:     if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2282:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2283:       PetscFVCellGeom       *cg1;
2284:       PetscFVFaceGeom       *fg;
2285:       const PetscInt        *fcells;
2286:       PetscInt               ncell, side;

2288:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2289:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2290:       if ((ghost >= 0) || boundary) continue;
2291:       DMPlexGetSupport(dm, faces[f], &fcells);
2292:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2293:       ncell = fcells[!side];    /* the neighbor */
2294:       DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2295:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2296:       for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2297:       gref[usedFaces++] = fg->grad[side];  /* Gradient reconstruction term will go here */
2298:     }
2299:     if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2300:     PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2301:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2302:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2303:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2304:       if ((ghost >= 0) || boundary) continue;
2305:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2306:       ++usedFaces;
2307:     }
2308:   }
2309:   PetscFree3(dx, grad, gref);
2310:   return(0);
2311: }

2313: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2314: {
2315:   DMLabel        ghostLabel;
2316:   PetscScalar   *dx, *grad, **gref;
2317:   PetscInt       dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2318:   PetscSection   neighSec;
2319:   PetscInt     (*neighbors)[2];
2320:   PetscInt      *counter;

2324:   DMGetDimension(dm, &dim);
2325:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2326:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2327:   if (cEndInterior < 0) cEndInterior = cEnd;
2328:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2329:   PetscSectionSetChart(neighSec,cStart,cEndInterior);
2330:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2331:   DMGetLabel(dm, "ghost", &ghostLabel);
2332:   for (f = fStart; f < fEnd; f++) {
2333:     const PetscInt        *fcells;
2334:     PetscBool              boundary;
2335:     PetscInt               ghost = -1;
2336:     PetscInt               numChildren, numCells, c;

2338:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2339:     DMIsBoundaryPoint(dm, f, &boundary);
2340:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2341:     if ((ghost >= 0) || boundary || numChildren) continue;
2342:     DMPlexGetSupportSize(dm, f, &numCells);
2343:     if (numCells == 2) {
2344:       DMPlexGetSupport(dm, f, &fcells);
2345:       for (c = 0; c < 2; c++) {
2346:         PetscInt cell = fcells[c];

2348:         if (cell >= cStart && cell < cEndInterior) {
2349:           PetscSectionAddDof(neighSec,cell,1);
2350:         }
2351:       }
2352:     }
2353:   }
2354:   PetscSectionSetUp(neighSec);
2355:   PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2356:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2357:   nStart = 0;
2358:   PetscSectionGetStorageSize(neighSec,&nEnd);
2359:   PetscMalloc1((nEnd-nStart),&neighbors);
2360:   PetscCalloc1((cEndInterior-cStart),&counter);
2361:   for (f = fStart; f < fEnd; f++) {
2362:     const PetscInt        *fcells;
2363:     PetscBool              boundary;
2364:     PetscInt               ghost = -1;
2365:     PetscInt               numChildren, numCells, c;

2367:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2368:     DMIsBoundaryPoint(dm, f, &boundary);
2369:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2370:     if ((ghost >= 0) || boundary || numChildren) continue;
2371:     DMPlexGetSupportSize(dm, f, &numCells);
2372:     if (numCells == 2) {
2373:       DMPlexGetSupport(dm, f, &fcells);
2374:       for (c = 0; c < 2; c++) {
2375:         PetscInt cell = fcells[c], off;

2377:         if (cell >= cStart && cell < cEndInterior) {
2378:           PetscSectionGetOffset(neighSec,cell,&off);
2379:           off += counter[cell - cStart]++;
2380:           neighbors[off][0] = f;
2381:           neighbors[off][1] = fcells[1 - c];
2382:         }
2383:       }
2384:     }
2385:   }
2386:   PetscFree(counter);
2387:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2388:   for (c = cStart; c < cEndInterior; c++) {
2389:     PetscInt               numFaces, f, d, off, ghost = -1;
2390:     PetscFVCellGeom        *cg;

2392:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2393:     PetscSectionGetDof(neighSec, c, &numFaces);
2394:     PetscSectionGetOffset(neighSec, c, &off);
2395:     if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2396:     if (ghost < 0 && numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2397:     for (f = 0; f < numFaces; ++f) {
2398:       PetscFVCellGeom       *cg1;
2399:       PetscFVFaceGeom       *fg;
2400:       const PetscInt        *fcells;
2401:       PetscInt               ncell, side, nface;

2403:       nface = neighbors[off + f][0];
2404:       ncell = neighbors[off + f][1];
2405:       DMPlexGetSupport(dm,nface,&fcells);
2406:       side  = (c != fcells[0]);
2407:       DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2408:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2409:       for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2410:       gref[f] = fg->grad[side];  /* Gradient reconstruction term will go here */
2411:     }
2412:     PetscFVComputeGradient(fvm, numFaces, dx, grad);
2413:     for (f = 0; f < numFaces; ++f) {
2414:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2415:     }
2416:   }
2417:   PetscFree3(dx, grad, gref);
2418:   PetscSectionDestroy(&neighSec);
2419:   PetscFree(neighbors);
2420:   return(0);
2421: }

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

2426:   Collective on dm

2428:   Input Arguments:
2429: + dm  - The DM
2430: . fvm - The PetscFV
2431: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2432: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()

2434:   Output Parameters:
2435: + faceGeometry - The geometric factors for gradient calculation are inserted
2436: - dmGrad - The DM describing the layout of gradient data

2438:   Level: developer

2440: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2441: @*/
2442: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2443: {
2444:   DM             dmFace, dmCell;
2445:   PetscScalar   *fgeom, *cgeom;
2446:   PetscSection   sectionGrad, parentSection;
2447:   PetscInt       dim, pdim, cStart, cEnd, cEndInterior, c;

2451:   DMGetDimension(dm, &dim);
2452:   PetscFVGetNumComponents(fvm, &pdim);
2453:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2454:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2455:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2456:   VecGetDM(faceGeometry, &dmFace);
2457:   VecGetDM(cellGeometry, &dmCell);
2458:   VecGetArray(faceGeometry, &fgeom);
2459:   VecGetArray(cellGeometry, &cgeom);
2460:   DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2461:   if (!parentSection) {
2462:     BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2463:   } else {
2464:     BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2465:   }
2466:   VecRestoreArray(faceGeometry, &fgeom);
2467:   VecRestoreArray(cellGeometry, &cgeom);
2468:   /* Create storage for gradients */
2469:   DMClone(dm, dmGrad);
2470:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionGrad);
2471:   PetscSectionSetChart(sectionGrad, cStart, cEnd);
2472:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2473:   PetscSectionSetUp(sectionGrad);
2474:   DMSetLocalSection(*dmGrad, sectionGrad);
2475:   PetscSectionDestroy(&sectionGrad);
2476:   return(0);
2477: }

2479: /*@
2480:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2482:   Collective on dm

2484:   Input Arguments:
2485: + dm  - The DM
2486: - fvm - The PetscFV

2488:   Output Parameters:
2489: + cellGeometry - The cell geometry
2490: . faceGeometry - The face geometry
2491: - dmGrad       - The gradient matrices

2493:   Level: developer

2495: .seealso: DMPlexComputeGeometryFVM()
2496: @*/
2497: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2498: {
2499:   PetscObject    cellgeomobj, facegeomobj;

2503:   PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2504:   if (!cellgeomobj) {
2505:     Vec cellgeomInt, facegeomInt;

2507:     DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2508:     PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2509:     PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2510:     VecDestroy(&cellgeomInt);
2511:     VecDestroy(&facegeomInt);
2512:     PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2513:   }
2514:   PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2515:   if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2516:   if (facegeom) *facegeom = (Vec) facegeomobj;
2517:   if (gradDM) {
2518:     PetscObject gradobj;
2519:     PetscBool   computeGradients;

2521:     PetscFVGetComputeGradients(fv,&computeGradients);
2522:     if (!computeGradients) {
2523:       *gradDM = NULL;
2524:       return(0);
2525:     }
2526:     PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2527:     if (!gradobj) {
2528:       DM dmGradInt;

2530:       DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2531:       PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2532:       DMDestroy(&dmGradInt);
2533:       PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2534:     }
2535:     *gradDM = (DM) gradobj;
2536:   }
2537:   return(0);
2538: }

2540: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work,  PetscReal *resNeg, PetscReal *guess)
2541: {
2542:   PetscInt l, m;

2545:   if (dimC == dimR && dimR <= 3) {
2546:     /* invert Jacobian, multiply */
2547:     PetscScalar det, idet;

2549:     switch (dimR) {
2550:     case 1:
2551:       invJ[0] = 1./ J[0];
2552:       break;
2553:     case 2:
2554:       det = J[0] * J[3] - J[1] * J[2];
2555:       idet = 1./det;
2556:       invJ[0] =  J[3] * idet;
2557:       invJ[1] = -J[1] * idet;
2558:       invJ[2] = -J[2] * idet;
2559:       invJ[3] =  J[0] * idet;
2560:       break;
2561:     case 3:
2562:       {
2563:         invJ[0] = J[4] * J[8] - J[5] * J[7];
2564:         invJ[1] = J[2] * J[7] - J[1] * J[8];
2565:         invJ[2] = J[1] * J[5] - J[2] * J[4];
2566:         det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2567:         idet = 1./det;
2568:         invJ[0] *= idet;
2569:         invJ[1] *= idet;
2570:         invJ[2] *= idet;
2571:         invJ[3]  = idet * (J[5] * J[6] - J[3] * J[8]);
2572:         invJ[4]  = idet * (J[0] * J[8] - J[2] * J[6]);
2573:         invJ[5]  = idet * (J[2] * J[3] - J[0] * J[5]);
2574:         invJ[6]  = idet * (J[3] * J[7] - J[4] * J[6]);
2575:         invJ[7]  = idet * (J[1] * J[6] - J[0] * J[7]);
2576:         invJ[8]  = idet * (J[0] * J[4] - J[1] * J[3]);
2577:       }
2578:       break;
2579:     }
2580:     for (l = 0; l < dimR; l++) {
2581:       for (m = 0; m < dimC; m++) {
2582:         guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2583:       }
2584:     }
2585:   } else {
2586: #if defined(PETSC_USE_COMPLEX)
2587:     char transpose = 'C';
2588: #else
2589:     char transpose = 'T';
2590: #endif
2591:     PetscBLASInt m = dimR;
2592:     PetscBLASInt n = dimC;
2593:     PetscBLASInt one = 1;
2594:     PetscBLASInt worksize = dimR * dimC, info;

2596:     for (l = 0; l < dimC; l++) {invJ[l] = resNeg[l];}

2598:     PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&m,&n,&one,J,&m,invJ,&n,work,&worksize, &info));
2599:     if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS");

2601:     for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2602:   }
2603:   return(0);
2604: }

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

2616:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2617:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2618:   DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2619:   DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2620:   cellCoords = &cellData[0];
2621:   cellCoeffs = &cellData[coordSize];
2622:   extJ       = &cellData[2 * coordSize];
2623:   resNeg     = &cellData[2 * coordSize + dimR];
2624:   invJ       = &J[dimR * dimC];
2625:   work       = &J[2 * dimR * dimC];
2626:   if (dimR == 2) {
2627:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2632:       for (j = 0; j < dimC; j++) {
2633:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2634:       }
2635:     }
2636:   } else if (dimR == 3) {
2637:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2642:       for (j = 0; j < dimC; j++) {
2643:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2644:       }
2645:     }
2646:   } else {
2647:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2648:   }
2649:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2650:   for (i = 0; i < dimR; i++) {
2651:     PetscReal *swap;

2653:     for (j = 0; j < (numV / 2); j++) {
2654:       for (k = 0; k < dimC; k++) {
2655:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2656:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2657:       }
2658:     }

2660:     if (i < dimR - 1) {
2661:       swap = cellCoeffs;
2662:       cellCoeffs = cellCoords;
2663:       cellCoords = swap;
2664:     }
2665:   }
2666:   PetscArrayzero(refCoords,numPoints * dimR);
2667:   for (j = 0; j < numPoints; j++) {
2668:     for (i = 0; i < maxIts; i++) {
2669:       PetscReal *guess = &refCoords[dimR * j];

2671:       /* compute -residual and Jacobian */
2672:       for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2673:       for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2674:       for (k = 0; k < numV; k++) {
2675:         PetscReal extCoord = 1.;
2676:         for (l = 0; l < dimR; l++) {
2677:           PetscReal coord = guess[l];
2678:           PetscInt  dep   = (k & (1 << l)) >> l;

2680:           extCoord *= dep * coord + !dep;
2681:           extJ[l] = dep;

2683:           for (m = 0; m < dimR; m++) {
2684:             PetscReal coord = guess[m];
2685:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
2686:             PetscReal mult  = dep * coord + !dep;

2688:             extJ[l] *= mult;
2689:           }
2690:         }
2691:         for (l = 0; l < dimC; l++) {
2692:           PetscReal coeff = cellCoeffs[dimC * k + l];

2694:           resNeg[l] -= coeff * extCoord;
2695:           for (m = 0; m < dimR; m++) {
2696:             J[dimR * l + m] += coeff * extJ[m];
2697:           }
2698:         }
2699:       }
2700:       if (0 && PetscDefined(USE_DEBUG)) {
2701:         PetscReal maxAbs = 0.;

2703:         for (l = 0; l < dimC; l++) {
2704:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2705:         }
2706:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2707:       }

2709:       DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2710:     }
2711:   }
2712:   DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2713:   DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2714:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2715:   return(0);
2716: }

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

2727:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2728:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2729:   DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2730:   cellCoords = &cellData[0];
2731:   cellCoeffs = &cellData[coordSize];
2732:   if (dimR == 2) {
2733:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2738:       for (j = 0; j < dimC; j++) {
2739:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2740:       }
2741:     }
2742:   } else if (dimR == 3) {
2743:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2748:       for (j = 0; j < dimC; j++) {
2749:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2750:       }
2751:     }
2752:   } else {
2753:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2754:   }
2755:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2756:   for (i = 0; i < dimR; i++) {
2757:     PetscReal *swap;

2759:     for (j = 0; j < (numV / 2); j++) {
2760:       for (k = 0; k < dimC; k++) {
2761:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2762:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2763:       }
2764:     }

2766:     if (i < dimR - 1) {
2767:       swap = cellCoeffs;
2768:       cellCoeffs = cellCoords;
2769:       cellCoords = swap;
2770:     }
2771:   }
2772:   PetscArrayzero(realCoords,numPoints * dimC);
2773:   for (j = 0; j < numPoints; j++) {
2774:     const PetscReal *guess  = &refCoords[dimR * j];
2775:     PetscReal       *mapped = &realCoords[dimC * j];

2777:     for (k = 0; k < numV; k++) {
2778:       PetscReal extCoord = 1.;
2779:       for (l = 0; l < dimR; l++) {
2780:         PetscReal coord = guess[l];
2781:         PetscInt  dep   = (k & (1 << l)) >> l;

2783:         extCoord *= dep * coord + !dep;
2784:       }
2785:       for (l = 0; l < dimC; l++) {
2786:         PetscReal coeff = cellCoeffs[dimC * k + l];

2788:         mapped[l] += coeff * extCoord;
2789:       }
2790:     }
2791:   }
2792:   DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2793:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2794:   return(0);
2795: }

2797: /* TODO: TOBY please fix this for Nc > 1 */
2798: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2799: {
2800:   PetscInt       numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2801:   PetscScalar    *nodes = NULL;
2802:   PetscReal      *invV, *modes;
2803:   PetscReal      *B, *D, *resNeg;
2804:   PetscScalar    *J, *invJ, *work;

2808:   PetscFEGetDimension(fe, &pdim);
2809:   PetscFEGetNumComponents(fe, &numComp);
2810:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2811:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2812:   /* convert nodes to values in the stable evaluation basis */
2813:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2814:   invV = fe->invV;
2815:   for (i = 0; i < pdim; ++i) {
2816:     modes[i] = 0.;
2817:     for (j = 0; j < pdim; ++j) {
2818:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2819:     }
2820:   }
2821:   DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2822:   D      = &B[pdim*Nc];
2823:   resNeg = &D[pdim*Nc * dimR];
2824:   DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2825:   invJ = &J[Nc * dimR];
2826:   work = &invJ[Nc * dimR];
2827:   for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2828:   for (j = 0; j < numPoints; j++) {
2829:       for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2830:       PetscReal *guess = &refCoords[j * dimR];
2831:       PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2832:       for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2833:       for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2834:       for (k = 0; k < pdim; k++) {
2835:         for (l = 0; l < Nc; l++) {
2836:           resNeg[l] -= modes[k] * B[k * Nc + l];
2837:           for (m = 0; m < dimR; m++) {
2838:             J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2839:           }
2840:         }
2841:       }
2842:       if (0 && PetscDefined(USE_DEBUG)) {
2843:         PetscReal maxAbs = 0.;

2845:         for (l = 0; l < Nc; l++) {
2846:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2847:         }
2848:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2849:       }
2850:       DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2851:     }
2852:   }
2853:   DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2854:   DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2855:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2856:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2857:   return(0);
2858: }

2860: /* TODO: TOBY please fix this for Nc > 1 */
2861: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2862: {
2863:   PetscInt       numComp, pdim, i, j, k, l, coordSize;
2864:   PetscScalar    *nodes = NULL;
2865:   PetscReal      *invV, *modes;
2866:   PetscReal      *B;

2870:   PetscFEGetDimension(fe, &pdim);
2871:   PetscFEGetNumComponents(fe, &numComp);
2872:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2873:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2874:   /* convert nodes to values in the stable evaluation basis */
2875:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2876:   invV = fe->invV;
2877:   for (i = 0; i < pdim; ++i) {
2878:     modes[i] = 0.;
2879:     for (j = 0; j < pdim; ++j) {
2880:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2881:     }
2882:   }
2883:   DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2884:   PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2885:   for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2886:   for (j = 0; j < numPoints; j++) {
2887:     PetscReal *mapped = &realCoords[j * Nc];

2889:     for (k = 0; k < pdim; k++) {
2890:       for (l = 0; l < Nc; l++) {
2891:         mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2892:       }
2893:     }
2894:   }
2895:   DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2896:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2897:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2898:   return(0);
2899: }

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

2906:   Not collective

2908:   Input Parameters:
2909: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2910:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2911:                as a multilinear map for tensor-product elements
2912: . cell       - the cell whose map is used.
2913: . numPoints  - the number of points to locate
2914: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

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

2919:   Level: intermediate

2921: .seealso: DMPlexReferenceToCoordinates()
2922: @*/
2923: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2924: {
2925:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
2926:   DM             coordDM = NULL;
2927:   Vec            coords;
2928:   PetscFE        fe = NULL;

2933:   DMGetDimension(dm,&dimR);
2934:   DMGetCoordinateDim(dm,&dimC);
2935:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2936:   DMPlexGetDepth(dm,&depth);
2937:   DMGetCoordinatesLocal(dm,&coords);
2938:   DMGetCoordinateDM(dm,&coordDM);
2939:   if (coordDM) {
2940:     PetscInt coordFields;

2942:     DMGetNumFields(coordDM,&coordFields);
2943:     if (coordFields) {
2944:       PetscClassId id;
2945:       PetscObject  disc;

2947:       DMGetField(coordDM,0,NULL,&disc);
2948:       PetscObjectGetClassId(disc,&id);
2949:       if (id == PETSCFE_CLASSID) {
2950:         fe = (PetscFE) disc;
2951:       }
2952:     }
2953:   }
2954:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
2955:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
2956:   if (!fe) { /* implicit discretization: affine or multilinear */
2957:     PetscInt  coneSize;
2958:     PetscBool isSimplex, isTensor;

2960:     DMPlexGetConeSize(dm,cell,&coneSize);
2961:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
2962:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
2963:     if (isSimplex) {
2964:       PetscReal detJ, *v0, *J, *invJ;

2966:       DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2967:       J    = &v0[dimC];
2968:       invJ = &J[dimC * dimC];
2969:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
2970:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
2971:         const PetscReal x0[3] = {-1.,-1.,-1.};

2973:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
2974:       }
2975:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2976:     } else if (isTensor) {
2977:       DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2978:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
2979:   } else {
2980:     DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2981:   }
2982:   return(0);
2983: }

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

2988:   Not collective

2990:   Input Parameters:
2991: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2992:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2993:                as a multilinear map for tensor-product elements
2994: . cell       - the cell whose map is used.
2995: . numPoints  - the number of points to locate
2996: - refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

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

3001:    Level: intermediate

3003: .seealso: DMPlexCoordinatesToReference()
3004: @*/
3005: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3006: {
3007:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
3008:   DM             coordDM = NULL;
3009:   Vec            coords;
3010:   PetscFE        fe = NULL;

3015:   DMGetDimension(dm,&dimR);
3016:   DMGetCoordinateDim(dm,&dimC);
3017:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3018:   DMPlexGetDepth(dm,&depth);
3019:   DMGetCoordinatesLocal(dm,&coords);
3020:   DMGetCoordinateDM(dm,&coordDM);
3021:   if (coordDM) {
3022:     PetscInt coordFields;

3024:     DMGetNumFields(coordDM,&coordFields);
3025:     if (coordFields) {
3026:       PetscClassId id;
3027:       PetscObject  disc;

3029:       DMGetField(coordDM,0,NULL,&disc);
3030:       PetscObjectGetClassId(disc,&id);
3031:       if (id == PETSCFE_CLASSID) {
3032:         fe = (PetscFE) disc;
3033:       }
3034:     }
3035:   }
3036:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
3037:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3038:   if (!fe) { /* implicit discretization: affine or multilinear */
3039:     PetscInt  coneSize;
3040:     PetscBool isSimplex, isTensor;

3042:     DMPlexGetConeSize(dm,cell,&coneSize);
3043:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3044:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3045:     if (isSimplex) {
3046:       PetscReal detJ, *v0, *J;

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

3054:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3055:       }
3056:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3057:     } else if (isTensor) {
3058:       DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3059:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3060:   } else {
3061:     DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3062:   }
3063:   return(0);
3064: }

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

3069:   Not collective

3071:   Input Parameters:
3072: + dm      - The DM
3073: . time    - The time
3074: - func    - The function transforming current coordinates to new coordaintes

3076:    Calling sequence of func:
3077: $    func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3078: $         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3079: $         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3080: $         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);

3082: +  dim          - The spatial dimension
3083: .  Nf           - The number of input fields (here 1)
3084: .  NfAux        - The number of input auxiliary fields
3085: .  uOff         - The offset of the coordinates in u[] (here 0)
3086: .  uOff_x       - The offset of the coordinates in u_x[] (here 0)
3087: .  u            - The coordinate values at this point in space
3088: .  u_t          - The coordinate time derivative at this point in space (here NULL)
3089: .  u_x          - The coordinate derivatives at this point in space
3090: .  aOff         - The offset of each auxiliary field in u[]
3091: .  aOff_x       - The offset of each auxiliary field in u_x[]
3092: .  a            - The auxiliary field values at this point in space
3093: .  a_t          - The auxiliary field time derivative at this point in space (or NULL)
3094: .  a_x          - The auxiliary field derivatives at this point in space
3095: .  t            - The current time
3096: .  x            - The coordinates of this point (here not used)
3097: .  numConstants - The number of constants
3098: .  constants    - The value of each constant
3099: -  f            - The new coordinates at this point in space

3101:   Level: intermediate

3103: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3104: @*/
3105: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3106:                                    void (*func)(PetscInt, PetscInt, PetscInt,
3107:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3108:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3109:                                                 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3110: {
3111:   DM             cdm;
3112:   DMField        cf;
3113:   Vec            lCoords, tmpCoords;

3117:   DMGetCoordinateDM(dm, &cdm);
3118:   DMGetCoordinatesLocal(dm, &lCoords);
3119:   DMGetLocalVector(cdm, &tmpCoords);
3120:   VecCopy(lCoords, tmpCoords);
3121:   /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3122:   DMGetCoordinateField(dm, &cf);
3123:   cdm->coordinateField = cf;
3124:   DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3125:   cdm->coordinateField = NULL;
3126:   DMRestoreLocalVector(cdm, &tmpCoords);
3127:   DMSetCoordinatesLocal(dm, lCoords);
3128:   return(0);
3129: }

3131: /* Shear applies the transformation, assuming we fix z,
3132:   / 1  0  m_0 \
3133:   | 0  1  m_1 |
3134:   \ 0  0   1  /
3135: */
3136: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3137:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3138:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3139:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3140: {
3141:   const PetscInt Nc = uOff[1]-uOff[0];
3142:   const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3143:   PetscInt       c;

3145:   for (c = 0; c < Nc; ++c) {
3146:     coords[c] = u[c] + constants[c+1]*u[ax];
3147:   }
3148: }

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

3153:   Not collective

3155:   Input Parameters:
3156: + dm          - The DM
3157: . direction   - The shear coordinate direction, e.g. 0 is the x-axis
3158: - multipliers - The multiplier m for each direction which is not the shear direction

3160:   Level: intermediate

3162: .seealso: DMPlexRemapGeometry()
3163: @*/
3164: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3165: {
3166:   DM             cdm;
3167:   PetscDS        cds;
3168:   PetscObject    obj;
3169:   PetscClassId   id;
3170:   PetscScalar   *moduli;
3171:   const PetscInt dir = (PetscInt) direction;
3172:   PetscInt       dE, d, e;

3176:   DMGetCoordinateDM(dm, &cdm);
3177:   DMGetCoordinateDim(dm, &dE);
3178:   PetscMalloc1(dE+1, &moduli);
3179:   moduli[0] = dir;
3180:   for (d = 0, e = 0; d < dE; ++d) moduli[d+1] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3181:   DMGetDS(cdm, &cds);
3182:   PetscDSGetDiscretization(cds, 0, &obj);
3183:   PetscObjectGetClassId(obj, &id);
3184:   if (id != PETSCFE_CLASSID) {
3185:     Vec           lCoords;
3186:     PetscSection  cSection;
3187:     PetscScalar  *coords;
3188:     PetscInt      vStart, vEnd, v;

3190:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3191:     DMGetCoordinateSection(dm, &cSection);
3192:     DMGetCoordinatesLocal(dm, &lCoords);
3193:     VecGetArray(lCoords, &coords);
3194:     for (v = vStart; v < vEnd; ++v) {
3195:       PetscReal ds;
3196:       PetscInt  off, c;

3198:       PetscSectionGetOffset(cSection, v, &off);
3199:       ds   = PetscRealPart(coords[off+dir]);
3200:       for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3201:     }
3202:     VecRestoreArray(lCoords, &coords);
3203:   } else {
3204:     PetscDSSetConstants(cds, dE+1, moduli);
3205:     DMPlexRemapGeometry(dm, 0.0, f0_shear);
3206:   }
3207:   PetscFree(moduli);
3208:   return(0);
3209: }