Actual source code: plexgeometry.c

petsc-master 2020-02-18
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 defined(PETSC_USE_DEBUG)
 50:   /* check coordinate section is consistent with DM dimension */
 51:   {
 52:     PetscSection      cs;
 53:     PetscInt          ndof;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

311: /*
312:   PetscGridHashSetGrid - Divide the grid into boxes

314:   Not collective

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

321:   Level: developer

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

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

343: /*
344:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

346:   Not collective

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

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

357:   Level: developer

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

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

375:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
376:       if (dbox == -1   && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
377:       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",
378:                                              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);
379:       dboxes[p*dim+d] = dbox;
380:     }
381:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
382:   }
383:   return(0);
384: }

386: /*
387:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

389:  Not collective

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

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

401:   Level: developer

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

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

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

432: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
433: {

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

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

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

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

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

493: /*
494:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex

496:   Collective on dm

498:   Input Parameter:
499: . dm - The Plex

501:   Output Parameter:
502: . localBox - The grid hash object

504:   Level: developer

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

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

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

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

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

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

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

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

642:   PetscTime(&t0);
643:   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.");
644:   DMGetCoordinateDim(dm, &dim);
645:   VecGetBlockSize(v, &bs);
646:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
647:   if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
648:   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);
649:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
650:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
651:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
652:   VecGetLocalSize(v, &numPoints);
653:   VecGetArray(v, &a);
654:   numPoints /= bs;
655:   {
656:     const PetscSFNode *sf_cells;

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

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

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

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

722:     if (hash) {
723:       PetscBool found_box;

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

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

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

812:   Not collective

814:   Input Parameter:
815: . coords - The coordinates of a segment

817:   Output Parameters:
818: + coords - The new y-coordinate, and 0 for x
819: - R - The rotation which accomplishes the projection

821:   Level: developer

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

832:   R[0] = c; R[1] = -s;
833:   R[2] = s; R[3] =  c;
834:   coords[0] = 0.0;
835:   coords[1] = r;
836:   return(0);
837: }

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

842:   Not collective

844:   Input Parameter:
845: . coords - The coordinates of a segment

847:   Output Parameters:
848: + coords - The new y-coordinate, and 0 for x and z
849: - R - The rotation which accomplishes the projection

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

853:   Level: developer

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

866:   x *= rinv; y *= rinv; z *= rinv;
867:   if (x > 0.) {
868:     PetscReal inv1pX   = 1./ (1. + x);

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

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

886: /*@
887:   DMPlexComputeProjection3Dto2D - Rewrite coordinates to be the 2D projection of the 3D coordinates

889:   Not collective

891:   Input Parameter:
892: . coords - The coordinates of a segment

894:   Output Parameters:
895: + coords - The new y- and z-coordinates, and 0 for x
896: - R - The rotation which accomplishes the projection

898:   Level: developer

900: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
901: @*/
902: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
903: {
904:   PetscReal      x1[3],  x2[3], n[3], norm;
905:   PetscReal      x1p[3], x2p[3], xnp[3];
906:   PetscReal      sqrtz, alpha;
907:   const PetscInt dim = 3;
908:   PetscInt       d, e, p;

911:   /* 0) Calculate normal vector */
912:   for (d = 0; d < dim; ++d) {
913:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
914:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
915:   }
916:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
917:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
918:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
919:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
920:   n[0] /= norm;
921:   n[1] /= norm;
922:   n[2] /= norm;
923:   /* 1) Take the normal vector and rotate until it is \hat z

925:     Let the normal vector be <nx, ny, nz> and alpha = 1/sqrt(1 - nz^2), then

927:     R = /  alpha nx nz  alpha ny nz -1/alpha \
928:         | -alpha ny     alpha nx        0    |
929:         \     nx            ny         nz    /

931:     will rotate the normal vector to \hat z
932:   */
933:   sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
934:   /* Check for n = z */
935:   if (sqrtz < 1.0e-10) {
936:     const PetscInt s = PetscSign(n[2]);
937:     /* If nz < 0, rotate 180 degrees around x-axis */
938:     for (p = 3; p < coordSize/3; ++p) {
939:       coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
940:       coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
941:     }
942:     coords[0] = 0.0;
943:     coords[1] = 0.0;
944:     coords[2] = x1[0];
945:     coords[3] = x1[1] * s;
946:     coords[4] = x2[0];
947:     coords[5] = x2[1] * s;
948:     R[0] = 1.0;     R[1] = 0.0;     R[2] = 0.0;
949:     R[3] = 0.0;     R[4] = 1.0 * s; R[5] = 0.0;
950:     R[6] = 0.0;     R[7] = 0.0;     R[8] = 1.0 * s;
951:     return(0);
952:   }
953:   alpha = 1.0/sqrtz;
954:   R[0] =  alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
955:   R[3] = -alpha*n[1];      R[4] = alpha*n[0];      R[5] = 0.0;
956:   R[6] =  n[0];            R[7] = n[1];            R[8] = n[2];
957:   for (d = 0; d < dim; ++d) {
958:     x1p[d] = 0.0;
959:     x2p[d] = 0.0;
960:     for (e = 0; e < dim; ++e) {
961:       x1p[d] += R[d*dim+e]*x1[e];
962:       x2p[d] += R[d*dim+e]*x2[e];
963:     }
964:   }
965:   if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
966:   if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
967:   /* 2) Project to (x, y) */
968:   for (p = 3; p < coordSize/3; ++p) {
969:     for (d = 0; d < dim; ++d) {
970:       xnp[d] = 0.0;
971:       for (e = 0; e < dim; ++e) {
972:         xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
973:       }
974:       if (d < dim-1) coords[p*2+d] = xnp[d];
975:     }
976:   }
977:   coords[0] = 0.0;
978:   coords[1] = 0.0;
979:   coords[2] = x1p[0];
980:   coords[3] = x1p[1];
981:   coords[4] = x2p[0];
982:   coords[5] = x2p[1];
983:   /* Output R^T which rotates \hat z to the input normal */
984:   for (d = 0; d < dim; ++d) {
985:     for (e = d+1; e < dim; ++e) {
986:       PetscReal tmp;

988:       tmp        = R[d*dim+e];
989:       R[d*dim+e] = R[e*dim+d];
990:       R[e*dim+d] = tmp;
991:     }
992:   }
993:   return(0);
994: }

996: PETSC_UNUSED
997: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
998: {
999:   /* Signed volume is 1/2 the determinant

1001:    |  1  1  1 |
1002:    | x0 x1 x2 |
1003:    | y0 y1 y2 |

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

1007:    | x1 x2 |
1008:    | y1 y2 |
1009:   */
1010:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1011:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1012:   PetscReal       M[4], detM;
1013:   M[0] = x1; M[1] = x2;
1014:   M[2] = y1; M[3] = y2;
1015:   DMPlex_Det2D_Internal(&detM, M);
1016:   *vol = 0.5*detM;
1017:   (void)PetscLogFlops(5.0);
1018: }

1020: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
1021: {
1022:   DMPlex_Det2D_Internal(vol, coords);
1023:   *vol *= 0.5;
1024: }

1026: PETSC_UNUSED
1027: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1028: {
1029:   /* Signed volume is 1/6th of the determinant

1031:    |  1  1  1  1 |
1032:    | x0 x1 x2 x3 |
1033:    | y0 y1 y2 y3 |
1034:    | z0 z1 z2 z3 |

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

1038:    | x1 x2 x3 |
1039:    | y1 y2 y3 |
1040:    | z1 z2 z3 |
1041:   */
1042:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
1043:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
1044:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1045:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1046:   PetscReal       M[9], detM;
1047:   M[0] = x1; M[1] = x2; M[2] = x3;
1048:   M[3] = y1; M[4] = y2; M[5] = y3;
1049:   M[6] = z1; M[7] = z2; M[8] = z3;
1050:   DMPlex_Det3D_Internal(&detM, M);
1051:   *vol = -onesixth*detM;
1052:   (void)PetscLogFlops(10.0);
1053: }

1055: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1056: {
1057:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1058:   DMPlex_Det3D_Internal(vol, coords);
1059:   *vol *= -onesixth;
1060: }

1062: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1063: {
1064:   PetscSection   coordSection;
1065:   Vec            coordinates;
1066:   const PetscScalar *coords;
1067:   PetscInt       dim, d, off;

1071:   DMGetCoordinatesLocal(dm, &coordinates);
1072:   DMGetCoordinateSection(dm, &coordSection);
1073:   PetscSectionGetDof(coordSection,e,&dim);
1074:   if (!dim) return(0);
1075:   PetscSectionGetOffset(coordSection,e,&off);
1076:   VecGetArrayRead(coordinates,&coords);
1077:   if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1078:   VecRestoreArrayRead(coordinates,&coords);
1079:   *detJ = 1.;
1080:   if (J) {
1081:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1082:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1083:     if (invJ) {
1084:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1085:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1086:     }
1087:   }
1088:   return(0);
1089: }

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

1100:   DMGetCoordinatesLocal(dm, &coordinates);
1101:   DMGetCoordinateSection(dm, &coordSection);
1102:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1103:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1104:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1105:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1106:   if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1107:   *detJ = 0.0;
1108:   if (numCoords == 6) {
1109:     const PetscInt dim = 3;
1110:     PetscReal      R[9], J0;

1112:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1113:     DMPlexComputeProjection3Dto1D(coords, R);
1114:     if (J)    {
1115:       J0   = 0.5*PetscRealPart(coords[1]);
1116:       J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1117:       J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1118:       J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1119:       DMPlex_Det3D_Internal(detJ, J);
1120:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1121:     }
1122:   } else if (numCoords == 4) {
1123:     const PetscInt dim = 2;
1124:     PetscReal      R[4], J0;

1126:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1127:     DMPlexComputeProjection2Dto1D(coords, R);
1128:     if (J)    {
1129:       J0   = 0.5*PetscRealPart(coords[1]);
1130:       J[0] = R[0]*J0; J[1] = R[1];
1131:       J[2] = R[2]*J0; J[3] = R[3];
1132:       DMPlex_Det2D_Internal(detJ, J);
1133:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1134:     }
1135:   } else if (numCoords == 2) {
1136:     const PetscInt dim = 1;

1138:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1139:     if (J)    {
1140:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1141:       *detJ = J[0];
1142:       PetscLogFlops(2.0);
1143:       if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1144:     }
1145:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1146:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1147:   return(0);
1148: }

1150: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1151: {
1152:   PetscSection   coordSection;
1153:   Vec            coordinates;
1154:   PetscScalar   *coords = NULL;
1155:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1159:   DMGetCoordinatesLocal(dm, &coordinates);
1160:   DMGetCoordinateSection(dm, &coordSection);
1161:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1162:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1163:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1164:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1165:   *detJ = 0.0;
1166:   if (numCoords == 9) {
1167:     const PetscInt dim = 3;
1168:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1170:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1171:     DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1172:     if (J)    {
1173:       const PetscInt pdim = 2;

1175:       for (d = 0; d < pdim; d++) {
1176:         for (f = 0; f < pdim; f++) {
1177:           J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1178:         }
1179:       }
1180:       PetscLogFlops(8.0);
1181:       DMPlex_Det3D_Internal(detJ, J0);
1182:       for (d = 0; d < dim; d++) {
1183:         for (f = 0; f < dim; f++) {
1184:           J[d*dim+f] = 0.0;
1185:           for (g = 0; g < dim; g++) {
1186:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1187:           }
1188:         }
1189:       }
1190:       PetscLogFlops(18.0);
1191:     }
1192:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1193:   } else if (numCoords == 6) {
1194:     const PetscInt dim = 2;

1196:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1197:     if (J)    {
1198:       for (d = 0; d < dim; d++) {
1199:         for (f = 0; f < dim; f++) {
1200:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1201:         }
1202:       }
1203:       PetscLogFlops(8.0);
1204:       DMPlex_Det2D_Internal(detJ, J);
1205:     }
1206:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1207:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1208:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1209:   return(0);
1210: }

1212: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1213: {
1214:   PetscSection   coordSection;
1215:   Vec            coordinates;
1216:   PetscScalar   *coords = NULL;
1217:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1221:   DMGetCoordinatesLocal(dm, &coordinates);
1222:   DMGetCoordinateSection(dm, &coordSection);
1223:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1224:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1225:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1226:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1227:   if (!Nq) {
1228:     *detJ = 0.0;
1229:     if (numCoords == 12) {
1230:       const PetscInt dim = 3;
1231:       PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1233:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1234:       DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1235:       if (J)    {
1236:         const PetscInt pdim = 2;

1238:         for (d = 0; d < pdim; d++) {
1239:           J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1240:           J0[d*dim+1] = 0.5*(PetscRealPart(coords[2*pdim+d]) - PetscRealPart(coords[1*pdim+d]));
1241:         }
1242:         PetscLogFlops(8.0);
1243:         DMPlex_Det3D_Internal(detJ, J0);
1244:         for (d = 0; d < dim; d++) {
1245:           for (f = 0; f < dim; f++) {
1246:             J[d*dim+f] = 0.0;
1247:             for (g = 0; g < dim; g++) {
1248:               J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1249:             }
1250:           }
1251:         }
1252:         PetscLogFlops(18.0);
1253:       }
1254:       if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1255:     } else if (numCoords == 8) {
1256:       const PetscInt dim = 2;

1258:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1259:       if (J)    {
1260:         for (d = 0; d < dim; d++) {
1261:           J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1262:           J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1263:         }
1264:         PetscLogFlops(8.0);
1265:         DMPlex_Det2D_Internal(detJ, J);
1266:       }
1267:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1268:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1269:   } else {
1270:     const PetscInt Nv = 4;
1271:     const PetscInt dimR = 2;
1272:     const PetscInt zToPlex[4] = {0, 1, 3, 2};
1273:     PetscReal zOrder[12];
1274:     PetscReal zCoeff[12];
1275:     PetscInt  i, j, k, l, dim;

1277:     if (numCoords == 12) {
1278:       dim = 3;
1279:     } else if (numCoords == 8) {
1280:       dim = 2;
1281:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1282:     for (i = 0; i < Nv; i++) {
1283:       PetscInt zi = zToPlex[i];

1285:       for (j = 0; j < dim; j++) {
1286:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1287:       }
1288:     }
1289:     for (j = 0; j < dim; j++) {
1290:       zCoeff[dim * 0 + j] = 0.25 * (  zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1291:       zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1292:       zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1293:       zCoeff[dim * 3 + j] = 0.25 * (  zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1294:     }
1295:     for (i = 0; i < Nq; i++) {
1296:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1298:       if (v) {
1299:         PetscReal extPoint[4];

1301:         extPoint[0] = 1.;
1302:         extPoint[1] = xi;
1303:         extPoint[2] = eta;
1304:         extPoint[3] = xi * eta;
1305:         for (j = 0; j < dim; j++) {
1306:           PetscReal val = 0.;

1308:           for (k = 0; k < Nv; k++) {
1309:             val += extPoint[k] * zCoeff[dim * k + j];
1310:           }
1311:           v[i * dim + j] = val;
1312:         }
1313:       }
1314:       if (J) {
1315:         PetscReal extJ[8];

1317:         extJ[0] = 0.;
1318:         extJ[1] = 0.;
1319:         extJ[2] = 1.;
1320:         extJ[3] = 0.;
1321:         extJ[4] = 0.;
1322:         extJ[5] = 1.;
1323:         extJ[6] = eta;
1324:         extJ[7] = xi;
1325:         for (j = 0; j < dim; j++) {
1326:           for (k = 0; k < dimR; k++) {
1327:             PetscReal val = 0.;

1329:             for (l = 0; l < Nv; l++) {
1330:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1331:             }
1332:             J[i * dim * dim + dim * j + k] = val;
1333:           }
1334:         }
1335:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1336:           PetscReal x, y, z;
1337:           PetscReal *iJ = &J[i * dim * dim];
1338:           PetscReal norm;

1340:           x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1341:           y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1342:           z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1343:           norm = PetscSqrtReal(x * x + y * y + z * z);
1344:           iJ[2] = x / norm;
1345:           iJ[5] = y / norm;
1346:           iJ[8] = z / norm;
1347:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1348:           if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1349:         } else {
1350:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1351:           if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1352:         }
1353:       }
1354:     }
1355:   }
1356:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1357:   return(0);
1358: }

1360: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1361: {
1362:   PetscSection   coordSection;
1363:   Vec            coordinates;
1364:   PetscScalar   *coords = NULL;
1365:   const PetscInt dim = 3;
1366:   PetscInt       d;

1370:   DMGetCoordinatesLocal(dm, &coordinates);
1371:   DMGetCoordinateSection(dm, &coordSection);
1372:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1373:   *detJ = 0.0;
1374:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1375:   if (J)    {
1376:     for (d = 0; d < dim; d++) {
1377:       /* I orient with outward face normals */
1378:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1379:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1380:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1381:     }
1382:     PetscLogFlops(18.0);
1383:     DMPlex_Det3D_Internal(detJ, J);
1384:   }
1385:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1386:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1387:   return(0);
1388: }

1390: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1391: {
1392:   PetscSection   coordSection;
1393:   Vec            coordinates;
1394:   PetscScalar   *coords = NULL;
1395:   const PetscInt dim = 3;
1396:   PetscInt       d;

1400:   DMGetCoordinatesLocal(dm, &coordinates);
1401:   DMGetCoordinateSection(dm, &coordSection);
1402:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1403:   if (!Nq) {
1404:     *detJ = 0.0;
1405:     if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1406:     if (J)    {
1407:       for (d = 0; d < dim; d++) {
1408:         J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1409:         J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1410:         J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1411:       }
1412:       PetscLogFlops(18.0);
1413:       DMPlex_Det3D_Internal(detJ, J);
1414:     }
1415:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1416:   } else {
1417:     const PetscInt Nv = 8;
1418:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1419:     const PetscInt dim = 3;
1420:     const PetscInt dimR = 3;
1421:     PetscReal zOrder[24];
1422:     PetscReal zCoeff[24];
1423:     PetscInt  i, j, k, l;

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

1428:       for (j = 0; j < dim; j++) {
1429:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1430:       }
1431:     }
1432:     for (j = 0; j < dim; j++) {
1433:       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]);
1434:       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]);
1435:       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]);
1436:       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]);
1437:       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]);
1438:       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]);
1439:       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]);
1440:       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]);
1441:     }
1442:     for (i = 0; i < Nq; i++) {
1443:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

1445:       if (v) {
1446:         PetscReal extPoint[8];

1448:         extPoint[0] = 1.;
1449:         extPoint[1] = xi;
1450:         extPoint[2] = eta;
1451:         extPoint[3] = xi * eta;
1452:         extPoint[4] = theta;
1453:         extPoint[5] = theta * xi;
1454:         extPoint[6] = theta * eta;
1455:         extPoint[7] = theta * eta * xi;
1456:         for (j = 0; j < dim; j++) {
1457:           PetscReal val = 0.;

1459:           for (k = 0; k < Nv; k++) {
1460:             val += extPoint[k] * zCoeff[dim * k + j];
1461:           }
1462:           v[i * dim + j] = val;
1463:         }
1464:       }
1465:       if (J) {
1466:         PetscReal extJ[24];

1468:         extJ[0]  = 0.         ; extJ[1]  = 0.        ; extJ[2]  = 0.      ;
1469:         extJ[3]  = 1.         ; extJ[4]  = 0.        ; extJ[5]  = 0.      ;
1470:         extJ[6]  = 0.         ; extJ[7]  = 1.        ; extJ[8]  = 0.      ;
1471:         extJ[9]  = eta        ; extJ[10] = xi        ; extJ[11] = 0.      ;
1472:         extJ[12] = 0.         ; extJ[13] = 0.        ; extJ[14] = 1.      ;
1473:         extJ[15] = theta      ; extJ[16] = 0.        ; extJ[17] = xi      ;
1474:         extJ[18] = 0.         ; extJ[19] = theta     ; extJ[20] = eta     ;
1475:         extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;

1477:         for (j = 0; j < dim; j++) {
1478:           for (k = 0; k < dimR; k++) {
1479:             PetscReal val = 0.;

1481:             for (l = 0; l < Nv; l++) {
1482:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1483:             }
1484:             J[i * dim * dim + dim * j + k] = val;
1485:           }
1486:         }
1487:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1488:         if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1489:       }
1490:     }
1491:   }
1492:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1493:   return(0);
1494: }

1496: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1497: {
1498:   DMPolytopeType  ct;
1499:   PetscInt        depth, dim, coordDim, coneSize, i;
1500:   PetscInt        Nq = 0;
1501:   const PetscReal *points = NULL;
1502:   DMLabel         depthLabel;
1503:   PetscReal       xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1504:   PetscBool       isAffine = PETSC_TRUE;
1505:   PetscErrorCode  ierr;

1508:   DMPlexGetDepth(dm, &depth);
1509:   DMPlexGetConeSize(dm, cell, &coneSize);
1510:   DMPlexGetDepthLabel(dm, &depthLabel);
1511:   DMLabelGetValue(depthLabel, cell, &dim);
1512:   if (depth == 1 && dim == 1) {
1513:     DMGetDimension(dm, &dim);
1514:   }
1515:   DMGetCoordinateDim(dm, &coordDim);
1516:   if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1517:   if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1518:   DMPlexGetCellType(dm, cell, &ct);
1519:   switch (ct) {
1520:     case DM_POLYTOPE_POINT:
1521:     DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1522:     isAffine = PETSC_FALSE;
1523:     break;
1524:     case DM_POLYTOPE_SEGMENT:
1525:     if (Nq) {DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1526:     else    {DMPlexComputeLineGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1527:     break;
1528:     case DM_POLYTOPE_TRIANGLE:
1529:     if (Nq) {DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1530:     else    {DMPlexComputeTriangleGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1531:     break;
1532:     case DM_POLYTOPE_QUADRILATERAL:
1533:     DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1534:     isAffine = PETSC_FALSE;
1535:     break;
1536:     case DM_POLYTOPE_TETRAHEDRON:
1537:     if (Nq) {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);}
1538:     else    {DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v,  J,  invJ,  detJ);}
1539:     break;
1540:     case DM_POLYTOPE_HEXAHEDRON:
1541:     DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1542:     isAffine = PETSC_FALSE;
1543:     break;
1544:     default: SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No element geometry for cell %D with type %s", cell, DMPolytopeTypes[ct]);
1545:   }
1546:   if (isAffine && Nq) {
1547:     if (v) {
1548:       for (i = 0; i < Nq; i++) {
1549:         CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1550:       }
1551:     }
1552:     if (detJ) {
1553:       for (i = 0; i < Nq; i++) {
1554:         detJ[i] = detJ0;
1555:       }
1556:     }
1557:     if (J) {
1558:       PetscInt k;

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

1563:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1564:           J[k] = J0[j];
1565:         }
1566:       }
1567:     }
1568:     if (invJ) {
1569:       PetscInt k;
1570:       switch (coordDim) {
1571:       case 0:
1572:         break;
1573:       case 1:
1574:         invJ[0] = 1./J0[0];
1575:         break;
1576:       case 2:
1577:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1578:         break;
1579:       case 3:
1580:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1581:         break;
1582:       }
1583:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1584:         PetscInt j;

1586:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1587:           invJ[k] = invJ[j];
1588:         }
1589:       }
1590:     }
1591:   }
1592:   return(0);
1593: }

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

1598:   Collective on dm

1600:   Input Arguments:
1601: + dm   - the DM
1602: - cell - the cell

1604:   Output Arguments:
1605: + v0   - the translation part of this affine transform
1606: . J    - the Jacobian of the transform from the reference element
1607: . invJ - the inverse of the Jacobian
1608: - detJ - the Jacobian determinant

1610:   Level: advanced

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

1616: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1617: @*/
1618: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1619: {

1623:   DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1624:   return(0);
1625: }

1627: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1628: {
1629:   PetscQuadrature   feQuad;
1630:   PetscSection      coordSection;
1631:   Vec               coordinates;
1632:   PetscScalar      *coords = NULL;
1633:   const PetscReal  *quadPoints;
1634:   PetscTabulation T;
1635:   PetscInt          dim, cdim, pdim, qdim, Nq, numCoords, q;
1636:   PetscErrorCode    ierr;

1639:   DMGetCoordinatesLocal(dm, &coordinates);
1640:   DMGetCoordinateSection(dm, &coordSection);
1641:   DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1642:   DMGetDimension(dm, &dim);
1643:   DMGetCoordinateDim(dm, &cdim);
1644:   if (!quad) { /* use the first point of the first functional of the dual space */
1645:     PetscDualSpace dsp;

1647:     PetscFEGetDualSpace(fe, &dsp);
1648:     PetscDualSpaceGetFunctional(dsp, 0, &quad);
1649:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1650:     Nq = 1;
1651:   } else {
1652:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1653:   }
1654:   PetscFEGetDimension(fe, &pdim);
1655:   PetscFEGetQuadrature(fe, &feQuad);
1656:   if (feQuad == quad) {
1657:     PetscFEGetCellTabulation(fe, &T);
1658:     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);
1659:   } else {
1660:     PetscFECreateTabulation(fe, 1, Nq, quadPoints, J ? 1 : 0, &T);
1661:   }
1662:   if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1663:   {
1664:     const PetscReal *basis    = T->T[0];
1665:     const PetscReal *basisDer = T->T[1];
1666:     PetscReal        detJt;

1668:     if (v) {
1669:       PetscArrayzero(v, Nq*cdim);
1670:       for (q = 0; q < Nq; ++q) {
1671:         PetscInt i, k;

1673:         for (k = 0; k < pdim; ++k)
1674:           for (i = 0; i < cdim; ++i)
1675:             v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1676:         PetscLogFlops(2.0*pdim*cdim);
1677:       }
1678:     }
1679:     if (J) {
1680:       PetscArrayzero(J, Nq*cdim*cdim);
1681:       for (q = 0; q < Nq; ++q) {
1682:         PetscInt i, j, k, c, r;

1684:         /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1685:         for (k = 0; k < pdim; ++k)
1686:           for (j = 0; j < dim; ++j)
1687:             for (i = 0; i < cdim; ++i)
1688:               J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1689:         PetscLogFlops(2.0*pdim*dim*cdim);
1690:         if (cdim > dim) {
1691:           for (c = dim; c < cdim; ++c)
1692:             for (r = 0; r < cdim; ++r)
1693:               J[r*cdim+c] = r == c ? 1.0 : 0.0;
1694:         }
1695:         if (!detJ && !invJ) continue;
1696:         detJt = 0.;
1697:         switch (cdim) {
1698:         case 3:
1699:           DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1700:           if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1701:           break;
1702:         case 2:
1703:           DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1704:           if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1705:           break;
1706:         case 1:
1707:           detJt = J[q*cdim*dim];
1708:           if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1709:         }
1710:         if (detJ) detJ[q] = detJt;
1711:       }
1712:     } else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1713:   }
1714:   if (feQuad != quad) {PetscTabulationDestroy(&T);}
1715:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1716:   return(0);
1717: }

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

1722:   Collective on dm

1724:   Input Arguments:
1725: + dm   - the DM
1726: . cell - the cell
1727: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If quad == NULL, geometry will be
1728:          evaluated at the first vertex of the reference element

1730:   Output Arguments:
1731: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1732: . J    - the Jacobian of the transform from the reference element at each quadrature point
1733: . invJ - the inverse of the Jacobian at each quadrature point
1734: - detJ - the Jacobian determinant at each quadrature point

1736:   Level: advanced

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

1742: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1743: @*/
1744: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1745: {
1746:   DM             cdm;
1747:   PetscFE        fe = NULL;

1752:   DMGetCoordinateDM(dm, &cdm);
1753:   if (cdm) {
1754:     PetscClassId id;
1755:     PetscInt     numFields;
1756:     PetscDS      prob;
1757:     PetscObject  disc;

1759:     DMGetNumFields(cdm, &numFields);
1760:     if (numFields) {
1761:       DMGetDS(cdm, &prob);
1762:       PetscDSGetDiscretization(prob,0,&disc);
1763:       PetscObjectGetClassId(disc,&id);
1764:       if (id == PETSCFE_CLASSID) {
1765:         fe = (PetscFE) disc;
1766:       }
1767:     }
1768:   }
1769:   if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1770:   else     {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1771:   return(0);
1772: }

1774: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1775: {
1776:   PetscSection   coordSection;
1777:   Vec            coordinates;
1778:   PetscScalar   *coords = NULL;
1779:   PetscScalar    tmp[2];
1780:   PetscInt       coordSize, d;

1784:   DMGetCoordinatesLocal(dm, &coordinates);
1785:   DMGetCoordinateSection(dm, &coordSection);
1786:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1787:   DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1788:   if (centroid) {
1789:     for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1790:   }
1791:   if (normal) {
1792:     PetscReal norm;

1794:     if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1795:     normal[0]  = -PetscRealPart(coords[1] - tmp[1]);
1796:     normal[1]  =  PetscRealPart(coords[0] - tmp[0]);
1797:     norm       = DMPlex_NormD_Internal(dim, normal);
1798:     for (d = 0; d < dim; ++d) normal[d] /= norm;
1799:   }
1800:   if (vol) {
1801:     *vol = 0.0;
1802:     for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1803:     *vol = PetscSqrtReal(*vol);
1804:   }
1805:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1806:   return(0);
1807: }

1809: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1810: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1811: {
1812:   DMLabel        depth;
1813:   PetscSection   coordSection;
1814:   Vec            coordinates;
1815:   PetscScalar   *coords = NULL;
1816:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1817:   PetscBool      isHybrid = PETSC_FALSE;
1818:   PetscInt       fv[4] = {0, 1, 2, 3};
1819:   PetscInt       pEndInterior[4], cdepth, tdim = 2, coordSize, numCorners, p, d, e;

1823:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1824:   DMPlexGetDepthLabel(dm, &depth);
1825:   DMLabelGetValue(depth, cell, &cdepth);
1826:   DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1827:   if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;
1828:   DMGetCoordinatesLocal(dm, &coordinates);
1829:   DMPlexGetConeSize(dm, cell, &numCorners);
1830:   DMGetCoordinateSection(dm, &coordSection);
1831:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1832:   DMGetCoordinateDim(dm, &dim);
1833:   /* Side faces for hybrid cells are are stored as tensor products */
1834:   if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}

1836:   if (dim > 2 && centroid) {
1837:     v0[0] = PetscRealPart(coords[0]);
1838:     v0[1] = PetscRealPart(coords[1]);
1839:     v0[2] = PetscRealPart(coords[2]);
1840:   }
1841:   if (normal) {
1842:     if (dim > 2) {
1843:       const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1844:       const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1845:       const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1846:       PetscReal       norm;

1848:       normal[0] = y0*z1 - z0*y1;
1849:       normal[1] = z0*x1 - x0*z1;
1850:       normal[2] = x0*y1 - y0*x1;
1851:       norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1852:       normal[0] /= norm;
1853:       normal[1] /= norm;
1854:       normal[2] /= norm;
1855:     } else {
1856:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
1857:     }
1858:   }
1859:   if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1860:   for (p = 0; p < numCorners; ++p) {
1861:     const PetscInt pi  = p < 4 ? fv[p] : p;
1862:     const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1863:     /* Need to do this copy to get types right */
1864:     for (d = 0; d < tdim; ++d) {
1865:       ctmp[d]      = PetscRealPart(coords[pi*tdim+d]);
1866:       ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1867:     }
1868:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1869:     vsum += vtmp;
1870:     for (d = 0; d < tdim; ++d) {
1871:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1872:     }
1873:   }
1874:   for (d = 0; d < tdim; ++d) {
1875:     csum[d] /= (tdim+1)*vsum;
1876:   }
1877:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1878:   if (vol) *vol = PetscAbsReal(vsum);
1879:   if (centroid) {
1880:     if (dim > 2) {
1881:       for (d = 0; d < dim; ++d) {
1882:         centroid[d] = v0[d];
1883:         for (e = 0; e < dim; ++e) {
1884:           centroid[d] += R[d*dim+e]*csum[e];
1885:         }
1886:       }
1887:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1888:   }
1889:   return(0);
1890: }

1892: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1893: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1894: {
1895:   DMLabel         depth;
1896:   PetscSection    coordSection;
1897:   Vec             coordinates;
1898:   PetscScalar    *coords = NULL;
1899:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
1900:   const PetscInt *faces, *facesO;
1901:   PetscBool       isHybrid = PETSC_FALSE;
1902:   PetscInt        pEndInterior[4], cdepth, numFaces, f, coordSize, p, d;
1903:   PetscErrorCode  ierr;

1906:   if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1907:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1908:   DMPlexGetDepthLabel(dm, &depth);
1909:   DMLabelGetValue(depth, cell, &cdepth);
1910:   DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1911:   if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;

1913:   DMGetCoordinatesLocal(dm, &coordinates);
1914:   DMGetCoordinateSection(dm, &coordSection);

1916:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1917:   DMPlexGetConeSize(dm, cell, &numFaces);
1918:   DMPlexGetCone(dm, cell, &faces);
1919:   DMPlexGetConeOrientation(dm, cell, &facesO);
1920:   for (f = 0; f < numFaces; ++f) {
1921:     PetscBool      flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */
1922:     DMPolytopeType ct;

1924:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1925:     DMPlexGetCellType(dm, faces[f], &ct);
1926:     switch (ct) {
1927:     case DM_POLYTOPE_TRIANGLE:
1928:       for (d = 0; d < dim; ++d) {
1929:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1930:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1931:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1932:       }
1933:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1934:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1935:       vsum += vtmp;
1936:       if (centroid) {           /* Centroid of OABC = (a+b+c)/4 */
1937:         for (d = 0; d < dim; ++d) {
1938:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1939:         }
1940:       }
1941:       break;
1942:     case DM_POLYTOPE_QUADRILATERAL:
1943:     {
1944:       PetscInt fv[4] = {0, 1, 2, 3};

1946:       /* Side faces for hybrid cells are are stored as tensor products */
1947:       if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1948:       /* DO FOR PYRAMID */
1949:       /* First tet */
1950:       for (d = 0; d < dim; ++d) {
1951:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1952:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1953:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1954:       }
1955:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1956:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1957:       vsum += vtmp;
1958:       if (centroid) {
1959:         for (d = 0; d < dim; ++d) {
1960:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1961:         }
1962:       }
1963:       /* Second tet */
1964:       for (d = 0; d < dim; ++d) {
1965:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1966:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1967:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1968:       }
1969:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1970:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1971:       vsum += vtmp;
1972:       if (centroid) {
1973:         for (d = 0; d < dim; ++d) {
1974:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1975:         }
1976:       }
1977:       break;
1978:     }
1979:     default:
1980:       SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle face %D of type %s", faces[f], DMPolytopeTypes[ct]);
1981:     }
1982:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1983:   }
1984:   if (vol)     *vol = PetscAbsReal(vsum);
1985:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
1986:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1987:   return(0);
1988: }

1990: /*@C
1991:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

1993:   Collective on dm

1995:   Input Arguments:
1996: + dm   - the DM
1997: - cell - the cell

1999:   Output Arguments:
2000: + volume   - the cell volume
2001: . centroid - the cell centroid
2002: - normal - the cell normal, if appropriate

2004:   Level: advanced

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

2010: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
2011: @*/
2012: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2013: {
2014:   PetscInt       depth, dim;

2018:   DMPlexGetDepth(dm, &depth);
2019:   DMGetDimension(dm, &dim);
2020:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2021:   DMPlexGetPointDepth(dm, cell, &depth);
2022:   switch (depth) {
2023:   case 1:
2024:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2025:     break;
2026:   case 2:
2027:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2028:     break;
2029:   case 3:
2030:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2031:     break;
2032:   default:
2033:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2034:   }
2035:   return(0);
2036: }

2038: /*@
2039:   DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh

2041:   Collective on dm

2043:   Input Parameter:
2044: . dm - The DMPlex

2046:   Output Parameter:
2047: . cellgeom - A vector with the cell geometry data for each cell

2049:   Level: beginner

2051: @*/
2052: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2053: {
2054:   DM             dmCell;
2055:   Vec            coordinates;
2056:   PetscSection   coordSection, sectionCell;
2057:   PetscScalar   *cgeom;
2058:   PetscInt       cStart, cEnd, cMax, c;

2062:   DMClone(dm, &dmCell);
2063:   DMGetCoordinateSection(dm, &coordSection);
2064:   DMGetCoordinatesLocal(dm, &coordinates);
2065:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2066:   DMSetCoordinatesLocal(dmCell, coordinates);
2067:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2068:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2069:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2070:   cEnd = cMax < 0 ? cEnd : cMax;
2071:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2072:   /* TODO This needs to be multiplied by Nq for non-affine */
2073:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2074:   PetscSectionSetUp(sectionCell);
2075:   DMSetLocalSection(dmCell, sectionCell);
2076:   PetscSectionDestroy(&sectionCell);
2077:   DMCreateLocalVector(dmCell, cellgeom);
2078:   VecGetArray(*cellgeom, &cgeom);
2079:   for (c = cStart; c < cEnd; ++c) {
2080:     PetscFEGeom *cg;

2082:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2083:     PetscArrayzero(cg, 1);
2084:     DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2085:     if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2086:   }
2087:   return(0);
2088: }

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

2093:   Input Parameter:
2094: . dm - The DM

2096:   Output Parameters:
2097: + cellgeom - A Vec of PetscFVCellGeom data
2098: - facegeom - A Vec of PetscFVFaceGeom data

2100:   Level: developer

2102: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2103: @*/
2104: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2105: {
2106:   DM             dmFace, dmCell;
2107:   DMLabel        ghostLabel;
2108:   PetscSection   sectionFace, sectionCell;
2109:   PetscSection   coordSection;
2110:   Vec            coordinates;
2111:   PetscScalar   *fgeom, *cgeom;
2112:   PetscReal      minradius, gminradius;
2113:   PetscInt       dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

2117:   DMGetDimension(dm, &dim);
2118:   DMGetCoordinateSection(dm, &coordSection);
2119:   DMGetCoordinatesLocal(dm, &coordinates);
2120:   /* Make cell centroids and volumes */
2121:   DMClone(dm, &dmCell);
2122:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2123:   DMSetCoordinatesLocal(dmCell, coordinates);
2124:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2125:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2126:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2127:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2128:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2129:   PetscSectionSetUp(sectionCell);
2130:   DMSetLocalSection(dmCell, sectionCell);
2131:   PetscSectionDestroy(&sectionCell);
2132:   DMCreateLocalVector(dmCell, cellgeom);
2133:   if (cEndInterior < 0) cEndInterior = cEnd;
2134:   VecGetArray(*cellgeom, &cgeom);
2135:   for (c = cStart; c < cEndInterior; ++c) {
2136:     PetscFVCellGeom *cg;

2138:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2139:     PetscArrayzero(cg, 1);
2140:     DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2141:   }
2142:   /* Compute face normals and minimum cell radius */
2143:   DMClone(dm, &dmFace);
2144:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionFace);
2145:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2146:   PetscSectionSetChart(sectionFace, fStart, fEnd);
2147:   for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2148:   PetscSectionSetUp(sectionFace);
2149:   DMSetLocalSection(dmFace, sectionFace);
2150:   PetscSectionDestroy(&sectionFace);
2151:   DMCreateLocalVector(dmFace, facegeom);
2152:   VecGetArray(*facegeom, &fgeom);
2153:   DMGetLabel(dm, "ghost", &ghostLabel);
2154:   minradius = PETSC_MAX_REAL;
2155:   for (f = fStart; f < fEnd; ++f) {
2156:     PetscFVFaceGeom *fg;
2157:     PetscReal        area;
2158:     PetscInt         ghost = -1, d, numChildren;

2160:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2161:     DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2162:     if (ghost >= 0 || numChildren) continue;
2163:     DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2164:     DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2165:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2166:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2167:     {
2168:       PetscFVCellGeom *cL, *cR;
2169:       PetscInt         ncells;
2170:       const PetscInt  *cells;
2171:       PetscReal       *lcentroid, *rcentroid;
2172:       PetscReal        l[3], r[3], v[3];

2174:       DMPlexGetSupport(dm, f, &cells);
2175:       DMPlexGetSupportSize(dm, f, &ncells);
2176:       DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2177:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2178:       if (ncells > 1) {
2179:         DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2180:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2181:       }
2182:       else {
2183:         rcentroid = fg->centroid;
2184:       }
2185:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2186:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2187:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2188:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2189:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2190:       }
2191:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2192:         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]);
2193:         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]);
2194:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2195:       }
2196:       if (cells[0] < cEndInterior) {
2197:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2198:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2199:       }
2200:       if (ncells > 1 && cells[1] < cEndInterior) {
2201:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2202:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2203:       }
2204:     }
2205:   }
2206:   MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2207:   DMPlexSetMinRadius(dm, gminradius);
2208:   /* Compute centroids of ghost cells */
2209:   for (c = cEndInterior; c < cEnd; ++c) {
2210:     PetscFVFaceGeom *fg;
2211:     const PetscInt  *cone,    *support;
2212:     PetscInt         coneSize, supportSize, s;

2214:     DMPlexGetConeSize(dmCell, c, &coneSize);
2215:     if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2216:     DMPlexGetCone(dmCell, c, &cone);
2217:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2218:     if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2219:     DMPlexGetSupport(dmCell, cone[0], &support);
2220:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2221:     for (s = 0; s < 2; ++s) {
2222:       /* Reflect ghost centroid across plane of face */
2223:       if (support[s] == c) {
2224:         PetscFVCellGeom       *ci;
2225:         PetscFVCellGeom       *cg;
2226:         PetscReal              c2f[3], a;

2228:         DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2229:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2230:         a    = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2231:         DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2232:         DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2233:         cg->volume = ci->volume;
2234:       }
2235:     }
2236:   }
2237:   VecRestoreArray(*facegeom, &fgeom);
2238:   VecRestoreArray(*cellgeom, &cgeom);
2239:   DMDestroy(&dmCell);
2240:   DMDestroy(&dmFace);
2241:   return(0);
2242: }

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

2247:   Not collective

2249:   Input Argument:
2250: . dm - the DM

2252:   Output Argument:
2253: . minradius - the minium cell radius

2255:   Level: developer

2257: .seealso: DMGetCoordinates()
2258: @*/
2259: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2260: {
2264:   *minradius = ((DM_Plex*) dm->data)->minradius;
2265:   return(0);
2266: }

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

2271:   Logically collective

2273:   Input Arguments:
2274: + dm - the DM
2275: - minradius - the minium cell radius

2277:   Level: developer

2279: .seealso: DMSetCoordinates()
2280: @*/
2281: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2282: {
2285:   ((DM_Plex*) dm->data)->minradius = minradius;
2286:   return(0);
2287: }

2289: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2290: {
2291:   DMLabel        ghostLabel;
2292:   PetscScalar   *dx, *grad, **gref;
2293:   PetscInt       dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

2297:   DMGetDimension(dm, &dim);
2298:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2299:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2300:   DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2301:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2302:   DMGetLabel(dm, "ghost", &ghostLabel);
2303:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2304:   for (c = cStart; c < cEndInterior; c++) {
2305:     const PetscInt        *faces;
2306:     PetscInt               numFaces, usedFaces, f, d;
2307:     PetscFVCellGeom        *cg;
2308:     PetscBool              boundary;
2309:     PetscInt               ghost;

2311:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2312:     DMPlexGetConeSize(dm, c, &numFaces);
2313:     DMPlexGetCone(dm, c, &faces);
2314:     if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2315:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2316:       PetscFVCellGeom       *cg1;
2317:       PetscFVFaceGeom       *fg;
2318:       const PetscInt        *fcells;
2319:       PetscInt               ncell, side;

2321:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2322:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2323:       if ((ghost >= 0) || boundary) continue;
2324:       DMPlexGetSupport(dm, faces[f], &fcells);
2325:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2326:       ncell = fcells[!side];    /* the neighbor */
2327:       DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2328:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2329:       for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2330:       gref[usedFaces++] = fg->grad[side];  /* Gradient reconstruction term will go here */
2331:     }
2332:     if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2333:     PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2334:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2335:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2336:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2337:       if ((ghost >= 0) || boundary) continue;
2338:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2339:       ++usedFaces;
2340:     }
2341:   }
2342:   PetscFree3(dx, grad, gref);
2343:   return(0);
2344: }

2346: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2347: {
2348:   DMLabel        ghostLabel;
2349:   PetscScalar   *dx, *grad, **gref;
2350:   PetscInt       dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2351:   PetscSection   neighSec;
2352:   PetscInt     (*neighbors)[2];
2353:   PetscInt      *counter;

2357:   DMGetDimension(dm, &dim);
2358:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2359:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2360:   if (cEndInterior < 0) cEndInterior = cEnd;
2361:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2362:   PetscSectionSetChart(neighSec,cStart,cEndInterior);
2363:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2364:   DMGetLabel(dm, "ghost", &ghostLabel);
2365:   for (f = fStart; f < fEnd; f++) {
2366:     const PetscInt        *fcells;
2367:     PetscBool              boundary;
2368:     PetscInt               ghost = -1;
2369:     PetscInt               numChildren, numCells, c;

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

2381:         if (cell >= cStart && cell < cEndInterior) {
2382:           PetscSectionAddDof(neighSec,cell,1);
2383:         }
2384:       }
2385:     }
2386:   }
2387:   PetscSectionSetUp(neighSec);
2388:   PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2389:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2390:   nStart = 0;
2391:   PetscSectionGetStorageSize(neighSec,&nEnd);
2392:   PetscMalloc1((nEnd-nStart),&neighbors);
2393:   PetscCalloc1((cEndInterior-cStart),&counter);
2394:   for (f = fStart; f < fEnd; f++) {
2395:     const PetscInt        *fcells;
2396:     PetscBool              boundary;
2397:     PetscInt               ghost = -1;
2398:     PetscInt               numChildren, numCells, c;

2400:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2401:     DMIsBoundaryPoint(dm, f, &boundary);
2402:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2403:     if ((ghost >= 0) || boundary || numChildren) continue;
2404:     DMPlexGetSupportSize(dm, f, &numCells);
2405:     if (numCells == 2) {
2406:       DMPlexGetSupport(dm, f, &fcells);
2407:       for (c = 0; c < 2; c++) {
2408:         PetscInt cell = fcells[c], off;

2410:         if (cell >= cStart && cell < cEndInterior) {
2411:           PetscSectionGetOffset(neighSec,cell,&off);
2412:           off += counter[cell - cStart]++;
2413:           neighbors[off][0] = f;
2414:           neighbors[off][1] = fcells[1 - c];
2415:         }
2416:       }
2417:     }
2418:   }
2419:   PetscFree(counter);
2420:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2421:   for (c = cStart; c < cEndInterior; c++) {
2422:     PetscInt               numFaces, f, d, off, ghost = -1;
2423:     PetscFVCellGeom        *cg;

2425:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2426:     PetscSectionGetDof(neighSec, c, &numFaces);
2427:     PetscSectionGetOffset(neighSec, c, &off);
2428:     if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2429:     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);
2430:     for (f = 0; f < numFaces; ++f) {
2431:       PetscFVCellGeom       *cg1;
2432:       PetscFVFaceGeom       *fg;
2433:       const PetscInt        *fcells;
2434:       PetscInt               ncell, side, nface;

2436:       nface = neighbors[off + f][0];
2437:       ncell = neighbors[off + f][1];
2438:       DMPlexGetSupport(dm,nface,&fcells);
2439:       side  = (c != fcells[0]);
2440:       DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2441:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2442:       for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2443:       gref[f] = fg->grad[side];  /* Gradient reconstruction term will go here */
2444:     }
2445:     PetscFVComputeGradient(fvm, numFaces, dx, grad);
2446:     for (f = 0; f < numFaces; ++f) {
2447:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2448:     }
2449:   }
2450:   PetscFree3(dx, grad, gref);
2451:   PetscSectionDestroy(&neighSec);
2452:   PetscFree(neighbors);
2453:   return(0);
2454: }

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

2459:   Collective on dm

2461:   Input Arguments:
2462: + dm  - The DM
2463: . fvm - The PetscFV
2464: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2465: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()

2467:   Output Parameters:
2468: + faceGeometry - The geometric factors for gradient calculation are inserted
2469: - dmGrad - The DM describing the layout of gradient data

2471:   Level: developer

2473: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2474: @*/
2475: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2476: {
2477:   DM             dmFace, dmCell;
2478:   PetscScalar   *fgeom, *cgeom;
2479:   PetscSection   sectionGrad, parentSection;
2480:   PetscInt       dim, pdim, cStart, cEnd, cEndInterior, c;

2484:   DMGetDimension(dm, &dim);
2485:   PetscFVGetNumComponents(fvm, &pdim);
2486:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2487:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2488:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2489:   VecGetDM(faceGeometry, &dmFace);
2490:   VecGetDM(cellGeometry, &dmCell);
2491:   VecGetArray(faceGeometry, &fgeom);
2492:   VecGetArray(cellGeometry, &cgeom);
2493:   DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2494:   if (!parentSection) {
2495:     BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2496:   } else {
2497:     BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2498:   }
2499:   VecRestoreArray(faceGeometry, &fgeom);
2500:   VecRestoreArray(cellGeometry, &cgeom);
2501:   /* Create storage for gradients */
2502:   DMClone(dm, dmGrad);
2503:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionGrad);
2504:   PetscSectionSetChart(sectionGrad, cStart, cEnd);
2505:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2506:   PetscSectionSetUp(sectionGrad);
2507:   DMSetLocalSection(*dmGrad, sectionGrad);
2508:   PetscSectionDestroy(&sectionGrad);
2509:   return(0);
2510: }

2512: /*@
2513:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2515:   Collective on dm

2517:   Input Arguments:
2518: + dm  - The DM
2519: - fvm - The PetscFV

2521:   Output Parameters:
2522: + cellGeometry - The cell geometry
2523: . faceGeometry - The face geometry
2524: - dmGrad       - The gradient matrices

2526:   Level: developer

2528: .seealso: DMPlexComputeGeometryFVM()
2529: @*/
2530: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2531: {
2532:   PetscObject    cellgeomobj, facegeomobj;

2536:   PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2537:   if (!cellgeomobj) {
2538:     Vec cellgeomInt, facegeomInt;

2540:     DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2541:     PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2542:     PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2543:     VecDestroy(&cellgeomInt);
2544:     VecDestroy(&facegeomInt);
2545:     PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2546:   }
2547:   PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2548:   if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2549:   if (facegeom) *facegeom = (Vec) facegeomobj;
2550:   if (gradDM) {
2551:     PetscObject gradobj;
2552:     PetscBool   computeGradients;

2554:     PetscFVGetComputeGradients(fv,&computeGradients);
2555:     if (!computeGradients) {
2556:       *gradDM = NULL;
2557:       return(0);
2558:     }
2559:     PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2560:     if (!gradobj) {
2561:       DM dmGradInt;

2563:       DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2564:       PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2565:       DMDestroy(&dmGradInt);
2566:       PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2567:     }
2568:     *gradDM = (DM) gradobj;
2569:   }
2570:   return(0);
2571: }

2573: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work,  PetscReal *resNeg, PetscReal *guess)
2574: {
2575:   PetscInt l, m;

2578:   if (dimC == dimR && dimR <= 3) {
2579:     /* invert Jacobian, multiply */
2580:     PetscScalar det, idet;

2582:     switch (dimR) {
2583:     case 1:
2584:       invJ[0] = 1./ J[0];
2585:       break;
2586:     case 2:
2587:       det = J[0] * J[3] - J[1] * J[2];
2588:       idet = 1./det;
2589:       invJ[0] =  J[3] * idet;
2590:       invJ[1] = -J[1] * idet;
2591:       invJ[2] = -J[2] * idet;
2592:       invJ[3] =  J[0] * idet;
2593:       break;
2594:     case 3:
2595:       {
2596:         invJ[0] = J[4] * J[8] - J[5] * J[7];
2597:         invJ[1] = J[2] * J[7] - J[1] * J[8];
2598:         invJ[2] = J[1] * J[5] - J[2] * J[4];
2599:         det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2600:         idet = 1./det;
2601:         invJ[0] *= idet;
2602:         invJ[1] *= idet;
2603:         invJ[2] *= idet;
2604:         invJ[3]  = idet * (J[5] * J[6] - J[3] * J[8]);
2605:         invJ[4]  = idet * (J[0] * J[8] - J[2] * J[6]);
2606:         invJ[5]  = idet * (J[2] * J[3] - J[0] * J[5]);
2607:         invJ[6]  = idet * (J[3] * J[7] - J[4] * J[6]);
2608:         invJ[7]  = idet * (J[1] * J[6] - J[0] * J[7]);
2609:         invJ[8]  = idet * (J[0] * J[4] - J[1] * J[3]);
2610:       }
2611:       break;
2612:     }
2613:     for (l = 0; l < dimR; l++) {
2614:       for (m = 0; m < dimC; m++) {
2615:         guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2616:       }
2617:     }
2618:   } else {
2619: #if defined(PETSC_USE_COMPLEX)
2620:     char transpose = 'C';
2621: #else
2622:     char transpose = 'T';
2623: #endif
2624:     PetscBLASInt m = dimR;
2625:     PetscBLASInt n = dimC;
2626:     PetscBLASInt one = 1;
2627:     PetscBLASInt worksize = dimR * dimC, info;

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

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

2634:     for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2635:   }
2636:   return(0);
2637: }

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

2649:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2650:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2651:   DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2652:   DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2653:   cellCoords = &cellData[0];
2654:   cellCoeffs = &cellData[coordSize];
2655:   extJ       = &cellData[2 * coordSize];
2656:   resNeg     = &cellData[2 * coordSize + dimR];
2657:   invJ       = &J[dimR * dimC];
2658:   work       = &J[2 * dimR * dimC];
2659:   if (dimR == 2) {
2660:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2665:       for (j = 0; j < dimC; j++) {
2666:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2667:       }
2668:     }
2669:   } else if (dimR == 3) {
2670:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2675:       for (j = 0; j < dimC; j++) {
2676:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2677:       }
2678:     }
2679:   } else {
2680:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2681:   }
2682:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2683:   for (i = 0; i < dimR; i++) {
2684:     PetscReal *swap;

2686:     for (j = 0; j < (numV / 2); j++) {
2687:       for (k = 0; k < dimC; k++) {
2688:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2689:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2690:       }
2691:     }

2693:     if (i < dimR - 1) {
2694:       swap = cellCoeffs;
2695:       cellCoeffs = cellCoords;
2696:       cellCoords = swap;
2697:     }
2698:   }
2699:   PetscArrayzero(refCoords,numPoints * dimR);
2700:   for (j = 0; j < numPoints; j++) {
2701:     for (i = 0; i < maxIts; i++) {
2702:       PetscReal *guess = &refCoords[dimR * j];

2704:       /* compute -residual and Jacobian */
2705:       for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2706:       for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2707:       for (k = 0; k < numV; k++) {
2708:         PetscReal extCoord = 1.;
2709:         for (l = 0; l < dimR; l++) {
2710:           PetscReal coord = guess[l];
2711:           PetscInt  dep   = (k & (1 << l)) >> l;

2713:           extCoord *= dep * coord + !dep;
2714:           extJ[l] = dep;

2716:           for (m = 0; m < dimR; m++) {
2717:             PetscReal coord = guess[m];
2718:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
2719:             PetscReal mult  = dep * coord + !dep;

2721:             extJ[l] *= mult;
2722:           }
2723:         }
2724:         for (l = 0; l < dimC; l++) {
2725:           PetscReal coeff = cellCoeffs[dimC * k + l];

2727:           resNeg[l] -= coeff * extCoord;
2728:           for (m = 0; m < dimR; m++) {
2729:             J[dimR * l + m] += coeff * extJ[m];
2730:           }
2731:         }
2732:       }
2733: #if 0 && defined(PETSC_USE_DEBUG)
2734:       {
2735:         PetscReal maxAbs = 0.;

2737:         for (l = 0; l < dimC; l++) {
2738:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2739:         }
2740:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2741:       }
2742: #endif

2744:       DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2745:     }
2746:   }
2747:   DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2748:   DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2749:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2750:   return(0);
2751: }

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

2762:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2763:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2764:   DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2765:   cellCoords = &cellData[0];
2766:   cellCoeffs = &cellData[coordSize];
2767:   if (dimR == 2) {
2768:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2773:       for (j = 0; j < dimC; j++) {
2774:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2775:       }
2776:     }
2777:   } else if (dimR == 3) {
2778:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2783:       for (j = 0; j < dimC; j++) {
2784:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2785:       }
2786:     }
2787:   } else {
2788:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2789:   }
2790:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2791:   for (i = 0; i < dimR; i++) {
2792:     PetscReal *swap;

2794:     for (j = 0; j < (numV / 2); j++) {
2795:       for (k = 0; k < dimC; k++) {
2796:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2797:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2798:       }
2799:     }

2801:     if (i < dimR - 1) {
2802:       swap = cellCoeffs;
2803:       cellCoeffs = cellCoords;
2804:       cellCoords = swap;
2805:     }
2806:   }
2807:   PetscArrayzero(realCoords,numPoints * dimC);
2808:   for (j = 0; j < numPoints; j++) {
2809:     const PetscReal *guess  = &refCoords[dimR * j];
2810:     PetscReal       *mapped = &realCoords[dimC * j];

2812:     for (k = 0; k < numV; k++) {
2813:       PetscReal extCoord = 1.;
2814:       for (l = 0; l < dimR; l++) {
2815:         PetscReal coord = guess[l];
2816:         PetscInt  dep   = (k & (1 << l)) >> l;

2818:         extCoord *= dep * coord + !dep;
2819:       }
2820:       for (l = 0; l < dimC; l++) {
2821:         PetscReal coeff = cellCoeffs[dimC * k + l];

2823:         mapped[l] += coeff * extCoord;
2824:       }
2825:     }
2826:   }
2827:   DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2828:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2829:   return(0);
2830: }

2832: /* TODO: TOBY please fix this for Nc > 1 */
2833: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2834: {
2835:   PetscInt       numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2836:   PetscScalar    *nodes = NULL;
2837:   PetscReal      *invV, *modes;
2838:   PetscReal      *B, *D, *resNeg;
2839:   PetscScalar    *J, *invJ, *work;

2843:   PetscFEGetDimension(fe, &pdim);
2844:   PetscFEGetNumComponents(fe, &numComp);
2845:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2846:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2847:   /* convert nodes to values in the stable evaluation basis */
2848:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2849:   invV = fe->invV;
2850:   for (i = 0; i < pdim; ++i) {
2851:     modes[i] = 0.;
2852:     for (j = 0; j < pdim; ++j) {
2853:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2854:     }
2855:   }
2856:   DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2857:   D      = &B[pdim*Nc];
2858:   resNeg = &D[pdim*Nc * dimR];
2859:   DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2860:   invJ = &J[Nc * dimR];
2861:   work = &invJ[Nc * dimR];
2862:   for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2863:   for (j = 0; j < numPoints; j++) {
2864:       for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2865:       PetscReal *guess = &refCoords[j * dimR];
2866:       PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2867:       for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2868:       for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2869:       for (k = 0; k < pdim; k++) {
2870:         for (l = 0; l < Nc; l++) {
2871:           resNeg[l] -= modes[k] * B[k * Nc + l];
2872:           for (m = 0; m < dimR; m++) {
2873:             J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2874:           }
2875:         }
2876:       }
2877: #if 0 && defined(PETSC_USE_DEBUG)
2878:       {
2879:         PetscReal maxAbs = 0.;

2881:         for (l = 0; l < Nc; l++) {
2882:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2883:         }
2884:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2885:       }
2886: #endif
2887:       DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2888:     }
2889:   }
2890:   DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2891:   DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2892:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2893:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2894:   return(0);
2895: }

2897: /* TODO: TOBY please fix this for Nc > 1 */
2898: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2899: {
2900:   PetscInt       numComp, pdim, i, j, k, l, coordSize;
2901:   PetscScalar    *nodes = NULL;
2902:   PetscReal      *invV, *modes;
2903:   PetscReal      *B;

2907:   PetscFEGetDimension(fe, &pdim);
2908:   PetscFEGetNumComponents(fe, &numComp);
2909:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2910:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2911:   /* convert nodes to values in the stable evaluation basis */
2912:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2913:   invV = fe->invV;
2914:   for (i = 0; i < pdim; ++i) {
2915:     modes[i] = 0.;
2916:     for (j = 0; j < pdim; ++j) {
2917:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2918:     }
2919:   }
2920:   DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2921:   PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2922:   for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2923:   for (j = 0; j < numPoints; j++) {
2924:     PetscReal *mapped = &realCoords[j * Nc];

2926:     for (k = 0; k < pdim; k++) {
2927:       for (l = 0; l < Nc; l++) {
2928:         mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2929:       }
2930:     }
2931:   }
2932:   DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2933:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2934:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2935:   return(0);
2936: }

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

2943:   Not collective

2945:   Input Parameters:
2946: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2947:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2948:                as a multilinear map for tensor-product elements
2949: . cell       - the cell whose map is used.
2950: . numPoints  - the number of points to locate
2951: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

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

2956:   Level: intermediate

2958: .seealso: DMPlexReferenceToCoordinates()
2959: @*/
2960: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2961: {
2962:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
2963:   DM             coordDM = NULL;
2964:   Vec            coords;
2965:   PetscFE        fe = NULL;

2970:   DMGetDimension(dm,&dimR);
2971:   DMGetCoordinateDim(dm,&dimC);
2972:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2973:   DMPlexGetDepth(dm,&depth);
2974:   DMGetCoordinatesLocal(dm,&coords);
2975:   DMGetCoordinateDM(dm,&coordDM);
2976:   if (coordDM) {
2977:     PetscInt coordFields;

2979:     DMGetNumFields(coordDM,&coordFields);
2980:     if (coordFields) {
2981:       PetscClassId id;
2982:       PetscObject  disc;

2984:       DMGetField(coordDM,0,NULL,&disc);
2985:       PetscObjectGetClassId(disc,&id);
2986:       if (id == PETSCFE_CLASSID) {
2987:         fe = (PetscFE) disc;
2988:       }
2989:     }
2990:   }
2991:   DMPlexGetInteriorCellStratum(dm,&cStart,&cEnd);
2992:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
2993:   if (!fe) { /* implicit discretization: affine or multilinear */
2994:     PetscInt  coneSize;
2995:     PetscBool isSimplex, isTensor;

2997:     DMPlexGetConeSize(dm,cell,&coneSize);
2998:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
2999:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3000:     if (isSimplex) {
3001:       PetscReal detJ, *v0, *J, *invJ;

3003:       DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3004:       J    = &v0[dimC];
3005:       invJ = &J[dimC * dimC];
3006:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3007:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3008:         const PetscReal x0[3] = {-1.,-1.,-1.};

3010:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3011:       }
3012:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3013:     } else if (isTensor) {
3014:       DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3015:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3016:   } else {
3017:     DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3018:   }
3019:   return(0);
3020: }

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

3025:   Not collective

3027:   Input Parameters:
3028: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3029:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3030:                as a multilinear map for tensor-product elements
3031: . cell       - the cell whose map is used.
3032: . numPoints  - the number of points to locate
3033: - refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

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

3038:    Level: intermediate

3040: .seealso: DMPlexCoordinatesToReference()
3041: @*/
3042: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3043: {
3044:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
3045:   DM             coordDM = NULL;
3046:   Vec            coords;
3047:   PetscFE        fe = NULL;

3052:   DMGetDimension(dm,&dimR);
3053:   DMGetCoordinateDim(dm,&dimC);
3054:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3055:   DMPlexGetDepth(dm,&depth);
3056:   DMGetCoordinatesLocal(dm,&coords);
3057:   DMGetCoordinateDM(dm,&coordDM);
3058:   if (coordDM) {
3059:     PetscInt coordFields;

3061:     DMGetNumFields(coordDM,&coordFields);
3062:     if (coordFields) {
3063:       PetscClassId id;
3064:       PetscObject  disc;

3066:       DMGetField(coordDM,0,NULL,&disc);
3067:       PetscObjectGetClassId(disc,&id);
3068:       if (id == PETSCFE_CLASSID) {
3069:         fe = (PetscFE) disc;
3070:       }
3071:     }
3072:   }
3073:   DMPlexGetInteriorCellStratum(dm,&cStart,&cEnd);
3074:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3075:   if (!fe) { /* implicit discretization: affine or multilinear */
3076:     PetscInt  coneSize;
3077:     PetscBool isSimplex, isTensor;

3079:     DMPlexGetConeSize(dm,cell,&coneSize);
3080:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3081:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3082:     if (isSimplex) {
3083:       PetscReal detJ, *v0, *J;

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

3091:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3092:       }
3093:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3094:     } else if (isTensor) {
3095:       DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3096:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3097:   } else {
3098:     DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3099:   }
3100:   return(0);
3101: }

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

3106:   Not collective

3108:   Input Parameters:
3109: + dm      - The DM
3110: . time    - The time
3111: - func    - The function transforming current coordinates to new coordaintes

3113:    Calling sequence of func:
3114: $    func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3115: $         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3116: $         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3117: $         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);

3119: +  dim          - The spatial dimension
3120: .  Nf           - The number of input fields (here 1)
3121: .  NfAux        - The number of input auxiliary fields
3122: .  uOff         - The offset of the coordinates in u[] (here 0)
3123: .  uOff_x       - The offset of the coordinates in u_x[] (here 0)
3124: .  u            - The coordinate values at this point in space
3125: .  u_t          - The coordinate time derivative at this point in space (here NULL)
3126: .  u_x          - The coordinate derivatives at this point in space
3127: .  aOff         - The offset of each auxiliary field in u[]
3128: .  aOff_x       - The offset of each auxiliary field in u_x[]
3129: .  a            - The auxiliary field values at this point in space
3130: .  a_t          - The auxiliary field time derivative at this point in space (or NULL)
3131: .  a_x          - The auxiliary field derivatives at this point in space
3132: .  t            - The current time
3133: .  x            - The coordinates of this point (here not used)
3134: .  numConstants - The number of constants
3135: .  constants    - The value of each constant
3136: -  f            - The new coordinates at this point in space

3138:   Level: intermediate

3140: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3141: @*/
3142: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3143:                                    void (*func)(PetscInt, PetscInt, PetscInt,
3144:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3145:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3146:                                                 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3147: {
3148:   DM             cdm;
3149:   DMField        cf;
3150:   Vec            lCoords, tmpCoords;

3154:   DMGetCoordinateDM(dm, &cdm);
3155:   DMGetCoordinatesLocal(dm, &lCoords);
3156:   DMGetLocalVector(cdm, &tmpCoords);
3157:   VecCopy(lCoords, tmpCoords);
3158:   /* We have to do the coordinate field manually right now since the coordinate DM will not have its own */
3159:   DMGetCoordinateField(dm, &cf);
3160:   cdm->coordinateField = cf;
3161:   DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3162:   cdm->coordinateField = NULL;
3163:   DMRestoreLocalVector(cdm, &tmpCoords);
3164:   return(0);
3165: }

3167: /* Shear applies the transformation, assuming we fix z,
3168:   / 1  0  m_0 \
3169:   | 0  1  m_1 |
3170:   \ 0  0   1  /
3171: */
3172: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3173:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3174:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3175:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3176: {
3177:   const PetscInt Nc = uOff[1]-uOff[0];
3178:   const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3179:   PetscInt       c;

3181:   for (c = 0; c < Nc; ++c) {
3182:     coords[c] = u[c] + constants[c+1]*u[ax];
3183:   }
3184: }

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

3189:   Not collective

3191:   Input Parameters:
3192: + dm          - The DM
3193: . direction   - The shear coordinate direction, e.g. 0 is the x-axis
3194: - multipliers - The multiplier m for each direction which is not the shear direction

3196:   Level: intermediate

3198: .seealso: DMPlexRemapGeometry()
3199: @*/
3200: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3201: {
3202:   DM             cdm;
3203:   PetscDS        cds;
3204:   PetscObject    obj;
3205:   PetscClassId   id;
3206:   PetscScalar   *moduli;
3207:   const PetscInt dir = (PetscInt) direction;
3208:   PetscInt       dE, d, e;

3212:   DMGetCoordinateDM(dm, &cdm);
3213:   DMGetCoordinateDim(dm, &dE);
3214:   PetscMalloc1(dE+1, &moduli);
3215:   moduli[0] = dir;
3216:   for (d = 0, e = 0; d < dE; ++d) moduli[d] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3217:   DMGetDS(cdm, &cds);
3218:   PetscDSGetDiscretization(cds, 0, &obj);
3219:   PetscObjectGetClassId(obj, &id);
3220:   if (id != PETSCFE_CLASSID) {
3221:     Vec           lCoords;
3222:     PetscSection  cSection;
3223:     PetscScalar  *coords;
3224:     PetscInt      vStart, vEnd, v;

3226:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3227:     DMGetCoordinateSection(dm, &cSection);
3228:     DMGetCoordinatesLocal(dm, &lCoords);
3229:     VecGetArray(lCoords, &coords);
3230:     for (v = vStart; v < vEnd; ++v) {
3231:       PetscReal ds;
3232:       PetscInt  off, c;

3234:       PetscSectionGetOffset(cSection, v, &off);
3235:       ds   = PetscRealPart(coords[off+dir]);
3236:       for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3237:     }
3238:     VecRestoreArray(lCoords, &coords);
3239:   } else {
3240:     PetscDSSetConstants(cds, dE+1, moduli);
3241:     DMPlexRemapGeometry(dm, 0.0, f0_shear);
3242:   }
3243:   PetscFree(moduli);
3244:   return(0);
3245: }