Actual source code: plexgeometry.c

petsc-master 2019-11-16
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: .seealso: DMPlexCreate(), DMGetCoordinates()
 32: @*/
 33: PetscErrorCode DMPlexFindVertices(DM dm, PetscInt npoints, const PetscReal coord[], PetscReal eps, PetscInt dagPoints[])
 34: {
 35:   PetscInt          c, dim, i, j, ndof, o, p, vStart, vEnd;
 36:   PetscSection      cs;
 37:   Vec               allCoordsVec;
 38:   const PetscScalar *allCoords;
 39:   PetscReal         norm;
 40:   PetscErrorCode    ierr;

 43:   if (eps < 0) eps = PETSC_SQRT_MACHINE_EPSILON;
 44:   DMGetDimension(dm, &dim);
 45:   DMGetCoordinateSection(dm, &cs);
 46:   DMGetCoordinatesLocal(dm, &allCoordsVec);
 47:   VecGetArrayRead(allCoordsVec, &allCoords);
 48:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 49:   for (i=0,j=0; i < npoints; i++,j+=dim) {
 50:     dagPoints[i] = -1;
 51:     for (p = vStart; p < vEnd; p++) {
 52:       PetscSectionGetOffset(cs, p, &o);
 53:       PetscSectionGetDof(cs, p, &ndof);
 54:       if (PetscUnlikely(ndof != dim)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point %D: ndof = %D != %D = dim", p, ndof, dim);
 55:       norm = 0.0;
 56:       for (c = 0; c < dim; c++) {
 57:         norm += PetscSqr(coord[j+c] - PetscRealPart(allCoords[o+c]));
 58:       }
 59:       norm = PetscSqrtReal(norm);
 60:       if (norm <= eps) {
 61:         dagPoints[i] = p;
 62:         break;
 63:       }
 64:     }
 65:   }
 66:   VecRestoreArrayRead(allCoordsVec, &allCoords);
 67:   return(0);
 68: }

 70: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
 71: {
 72:   const PetscReal p0_x  = segmentA[0*2+0];
 73:   const PetscReal p0_y  = segmentA[0*2+1];
 74:   const PetscReal p1_x  = segmentA[1*2+0];
 75:   const PetscReal p1_y  = segmentA[1*2+1];
 76:   const PetscReal p2_x  = segmentB[0*2+0];
 77:   const PetscReal p2_y  = segmentB[0*2+1];
 78:   const PetscReal p3_x  = segmentB[1*2+0];
 79:   const PetscReal p3_y  = segmentB[1*2+1];
 80:   const PetscReal s1_x  = p1_x - p0_x;
 81:   const PetscReal s1_y  = p1_y - p0_y;
 82:   const PetscReal s2_x  = p3_x - p2_x;
 83:   const PetscReal s2_y  = p3_y - p2_y;
 84:   const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

 87:   *hasIntersection = PETSC_FALSE;
 88:   /* Non-parallel lines */
 89:   if (denom != 0.0) {
 90:     const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
 91:     const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;

 93:     if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
 94:       *hasIntersection = PETSC_TRUE;
 95:       if (intersection) {
 96:         intersection[0] = p0_x + (t * s1_x);
 97:         intersection[1] = p0_y + (t * s1_y);
 98:       }
 99:     }
100:   }
101:   return(0);
102: }

104: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
105: {
106:   const PetscInt  embedDim = 2;
107:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
108:   PetscReal       x        = PetscRealPart(point[0]);
109:   PetscReal       y        = PetscRealPart(point[1]);
110:   PetscReal       v0[2], J[4], invJ[4], detJ;
111:   PetscReal       xi, eta;
112:   PetscErrorCode  ierr;

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

119:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
120:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
121:   return(0);
122: }

124: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
125: {
126:   const PetscInt  embedDim = 2;
127:   PetscReal       x        = PetscRealPart(point[0]);
128:   PetscReal       y        = PetscRealPart(point[1]);
129:   PetscReal       v0[2], J[4], invJ[4], detJ;
130:   PetscReal       xi, eta, r;
131:   PetscErrorCode  ierr;

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

138:   xi  = PetscMax(xi,  0.0);
139:   eta = PetscMax(eta, 0.0);
140:   if (xi + eta > 2.0) {
141:     r    = (xi + eta)/2.0;
142:     xi  /= r;
143:     eta /= r;
144:   }
145:   cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
146:   cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
147:   return(0);
148: }

150: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
151: {
152:   PetscSection       coordSection;
153:   Vec             coordsLocal;
154:   PetscScalar    *coords = NULL;
155:   const PetscInt  faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
156:   PetscReal       x         = PetscRealPart(point[0]);
157:   PetscReal       y         = PetscRealPart(point[1]);
158:   PetscInt        crossings = 0, f;
159:   PetscErrorCode  ierr;

162:   DMGetCoordinatesLocal(dm, &coordsLocal);
163:   DMGetCoordinateSection(dm, &coordSection);
164:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
165:   for (f = 0; f < 4; ++f) {
166:     PetscReal x_i   = PetscRealPart(coords[faces[2*f+0]*2+0]);
167:     PetscReal y_i   = PetscRealPart(coords[faces[2*f+0]*2+1]);
168:     PetscReal x_j   = PetscRealPart(coords[faces[2*f+1]*2+0]);
169:     PetscReal y_j   = PetscRealPart(coords[faces[2*f+1]*2+1]);
170:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
171:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
172:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
173:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
174:     if ((cond1 || cond2)  && above) ++crossings;
175:   }
176:   if (crossings % 2) *cell = c;
177:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
178:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
179:   return(0);
180: }

182: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
183: {
184:   const PetscInt embedDim = 3;
185:   PetscReal      v0[3], J[9], invJ[9], detJ;
186:   PetscReal      x = PetscRealPart(point[0]);
187:   PetscReal      y = PetscRealPart(point[1]);
188:   PetscReal      z = PetscRealPart(point[2]);
189:   PetscReal      xi, eta, zeta;

193:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
194:   xi   = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
195:   eta  = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
196:   zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);

198:   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
199:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
200:   return(0);
201: }

203: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
204: {
205:   PetscSection   coordSection;
206:   Vec            coordsLocal;
207:   PetscScalar   *coords = NULL;
208:   const PetscInt faces[24] = {0, 3, 2, 1,  5, 4, 7, 6,  3, 0, 4, 5,
209:                               1, 2, 6, 7,  3, 5, 6, 2,  0, 1, 7, 4};
210:   PetscBool      found = PETSC_TRUE;
211:   PetscInt       f;

215:   DMGetCoordinatesLocal(dm, &coordsLocal);
216:   DMGetCoordinateSection(dm, &coordSection);
217:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
218:   for (f = 0; f < 6; ++f) {
219:     /* Check the point is under plane */
220:     /*   Get face normal */
221:     PetscReal v_i[3];
222:     PetscReal v_j[3];
223:     PetscReal normal[3];
224:     PetscReal pp[3];
225:     PetscReal dot;

227:     v_i[0]    = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
228:     v_i[1]    = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
229:     v_i[2]    = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
230:     v_j[0]    = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
231:     v_j[1]    = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
232:     v_j[2]    = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
233:     normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
234:     normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
235:     normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
236:     pp[0]     = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
237:     pp[1]     = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
238:     pp[2]     = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
239:     dot       = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];

241:     /* Check that projected point is in face (2D location problem) */
242:     if (dot < 0.0) {
243:       found = PETSC_FALSE;
244:       break;
245:     }
246:   }
247:   if (found) *cell = c;
248:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
249:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
250:   return(0);
251: }

253: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
254: {
255:   PetscInt d;

258:   box->dim = dim;
259:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
260:   return(0);
261: }

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

268:   PetscMalloc1(1, box);
269:   PetscGridHashInitialize_Internal(*box, dim, point);
270:   return(0);
271: }

273: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
274: {
275:   PetscInt d;

278:   for (d = 0; d < box->dim; ++d) {
279:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
280:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
281:   }
282:   return(0);
283: }

285: /*
286:   PetscGridHashSetGrid - Divide the grid into boxes

288:   Not collective

290:   Input Parameters:
291: + box - The grid hash object
292: . n   - The number of boxes in each dimension, or PETSC_DETERMINE
293: - h   - The box size in each dimension, only used if n[d] == PETSC_DETERMINE

295:   Level: developer

297: .seealso: PetscGridHashCreate()
298: */
299: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
300: {
301:   PetscInt d;

304:   for (d = 0; d < box->dim; ++d) {
305:     box->extent[d] = box->upper[d] - box->lower[d];
306:     if (n[d] == PETSC_DETERMINE) {
307:       box->h[d] = h[d];
308:       box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
309:     } else {
310:       box->n[d] = n[d];
311:       box->h[d] = box->extent[d]/n[d];
312:     }
313:   }
314:   return(0);
315: }

317: /*
318:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

320:   Not collective

322:   Input Parameters:
323: + box       - The grid hash object
324: . numPoints - The number of input points
325: - points    - The input point coordinates

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

331:   Level: developer

333: .seealso: PetscGridHashCreate()
334: */
335: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
336: {
337:   const PetscReal *lower = box->lower;
338:   const PetscReal *upper = box->upper;
339:   const PetscReal *h     = box->h;
340:   const PetscInt  *n     = box->n;
341:   const PetscInt   dim   = box->dim;
342:   PetscInt         d, p;

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

349:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
350:       if (dbox == -1   && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
351:       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",
352:                                              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);
353:       dboxes[p*dim+d] = dbox;
354:     }
355:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
356:   }
357:   return(0);
358: }

360: /*
361:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

363:  Not collective

365:   Input Parameters:
366: + box       - The grid hash object
367: . numPoints - The number of input points
368: - points    - The input point coordinates

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

375:   Level: developer

377: .seealso: PetscGridHashGetEnclosingBox()
378: */
379: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
380: {
381:   const PetscReal *lower = box->lower;
382:   const PetscReal *upper = box->upper;
383:   const PetscReal *h     = box->h;
384:   const PetscInt  *n     = box->n;
385:   const PetscInt   dim   = box->dim;
386:   PetscInt         d, p;

389:   *found = PETSC_FALSE;
390:   for (p = 0; p < numPoints; ++p) {
391:     for (d = 0; d < dim; ++d) {
392:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);

394:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
395:       if (dbox < 0 || dbox >= n[d]) {
396:         return(0);
397:       }
398:       dboxes[p*dim+d] = dbox;
399:     }
400:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
401:   }
402:   *found = PETSC_TRUE;
403:   return(0);
404: }

406: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
407: {

411:   if (*box) {
412:     PetscSectionDestroy(&(*box)->cellSection);
413:     ISDestroy(&(*box)->cells);
414:     DMLabelDestroy(&(*box)->cellsSparse);
415:   }
416:   PetscFree(*box);
417:   return(0);
418: }

420: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
421: {
422:   PetscInt       coneSize;

426:   switch (dim) {
427:   case 2:
428:     DMPlexGetConeSize(dm, cellStart, &coneSize);
429:     switch (coneSize) {
430:     case 3:
431:       DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
432:       break;
433:     case 4:
434:       DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
435:       break;
436:     default:
437:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
438:     }
439:     break;
440:   case 3:
441:     DMPlexGetConeSize(dm, cellStart, &coneSize);
442:     switch (coneSize) {
443:     case 4:
444:       DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
445:       break;
446:     case 6:
447:       DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
448:       break;
449:     default:
450:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
451:     }
452:     break;
453:   default:
454:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
455:   }
456:   return(0);
457: }

459: /*
460:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
461: */
462: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
463: {
464:   PetscInt       coneSize;

468:   switch (dim) {
469:   case 2:
470:     DMPlexGetConeSize(dm, cell, &coneSize);
471:     switch (coneSize) {
472:     case 3:
473:       DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);
474:       break;
475: #if 0
476:     case 4:
477:       DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);
478:       break;
479: #endif
480:     default:
481:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
482:     }
483:     break;
484: #if 0
485:   case 3:
486:     DMPlexGetConeSize(dm, cell, &coneSize);
487:     switch (coneSize) {
488:     case 4:
489:       DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);
490:       break;
491:     case 6:
492:       DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);
493:       break;
494:     default:
495:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
496:     }
497:     break;
498: #endif
499:   default:
500:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim);
501:   }
502:   return(0);
503: }

505: /*
506:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex

508:   Collective on dm

510:   Input Parameter:
511: . dm - The Plex

513:   Output Parameter:
514: . localBox - The grid hash object

516:   Level: developer

518: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
519: */
520: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
521: {
522:   MPI_Comm           comm;
523:   PetscGridHash      lbox;
524:   Vec                coordinates;
525:   PetscSection       coordSection;
526:   Vec                coordsLocal;
527:   const PetscScalar *coords;
528:   PetscInt          *dboxes, *boxes;
529:   PetscInt           n[3] = {10, 10, 10};
530:   PetscInt           dim, N, cStart, cEnd, cMax, c, i;
531:   PetscErrorCode     ierr;

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

571:     /* Find boxes enclosing each vertex */
572:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
573:     PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
574:     /* Mark cells containing the vertices */
575:     for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
576:     /* Get grid of boxes containing these */
577:     for (d = 0;   d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
578:     for (d = dim; d < 3;   ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
579:     for (e = 1; e < dim+1; ++e) {
580:       for (d = 0; d < dim; ++d) {
581:         dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
582:         dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
583:       }
584:     }
585:     /* Check for intersection of box with cell */
586:     for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
587:       for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
588:         for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
589:           const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
590:           PetscScalar    cpoint[3];
591:           PetscInt       cell, edge, ii, jj, kk;

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

598:                 DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
599:                 if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
600:               }
601:             }
602:           }
603:           /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
604:           for (edge = 0; edge < dim+1; ++edge) {
605:             PetscReal segA[6], segB[6];

607:             if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
608:             for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
609:             for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
610:               if (dim > 2) {segB[2]     = PetscRealPart(point[2]);
611:                             segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
612:               for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
613:                 if (dim > 1) {segB[1]     = PetscRealPart(point[1]);
614:                               segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
615:                 for (ii = 0; ii < 2; ++ii) {
616:                   PetscBool intersects;

618:                   segB[0]     = PetscRealPart(point[0]);
619:                   segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
620:                   DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
621:                   if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
622:                 }
623:               }
624:             }
625:           }
626:         }
627:       }
628:     }
629:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
630:   }
631:   PetscFree2(dboxes, boxes);
632:   DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
633:   DMLabelDestroy(&lbox->cellsSparse);
634:   *localBox = lbox;
635:   return(0);
636: }

638: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
639: {
640:   DM_Plex        *mesh = (DM_Plex *) dm->data;
641:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
642:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
643:   PetscInt        dim, cStart, cEnd, cMax, numCells, c, d;
644:   const PetscInt *boxCells;
645:   PetscSFNode    *cells;
646:   PetscScalar    *a;
647:   PetscMPIInt     result;
648:   PetscLogDouble  t0,t1;
649:   PetscReal       gmin[3],gmax[3];
650:   PetscInt        terminating_query_type[] = { 0, 0, 0 };
651:   PetscErrorCode  ierr;

654:   PetscTime(&t0);
655:   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.");
656:   DMGetCoordinateDim(dm, &dim);
657:   VecGetBlockSize(v, &bs);
658:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
659:   if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
660:   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);
661:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
662:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
663:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
664:   VecGetLocalSize(v, &numPoints);
665:   VecGetArray(v, &a);
666:   numPoints /= bs;
667:   {
668:     const PetscSFNode *sf_cells;

670:     PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
671:     if (sf_cells) {
672:       PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
673:       cells = (PetscSFNode*)sf_cells;
674:       reuse = PETSC_TRUE;
675:     } else {
676:       PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
677:       PetscMalloc1(numPoints, &cells);
678:       /* initialize cells if created */
679:       for (p=0; p<numPoints; p++) {
680:         cells[p].rank  = 0;
681:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
682:       }
683:     }
684:   }
685:   /* define domain bounding box */
686:   {
687:     Vec coorglobal;

689:     DMGetCoordinates(dm,&coorglobal);
690:     VecStrideMaxAll(coorglobal,NULL,gmax);
691:     VecStrideMinAll(coorglobal,NULL,gmin);
692:   }
693:   if (hash) {
694:     if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
695:     /* Designate the local box for each point */
696:     /* Send points to correct process */
697:     /* Search cells that lie in each subbox */
698:     /*   Should we bin points before doing search? */
699:     ISGetIndices(mesh->lbox->cells, &boxCells);
700:   }
701:   for (p = 0, numFound = 0; p < numPoints; ++p) {
702:     const PetscScalar *point = &a[p*bs];
703:     PetscInt           dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
704:     PetscBool          point_outside_domain = PETSC_FALSE;

706:     /* check bounding box of domain */
707:     for (d=0; d<dim; d++) {
708:       if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
709:       if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
710:     }
711:     if (point_outside_domain) {
712:       cells[p].rank = 0;
713:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
714:       terminating_query_type[0]++;
715:       continue;
716:     }

718:     /* check initial values in cells[].index - abort early if found */
719:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
720:       c = cells[p].index;
721:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
722:       DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
723:       if (cell >= 0) {
724:         cells[p].rank = 0;
725:         cells[p].index = cell;
726:         numFound++;
727:       }
728:     }
729:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
730:       terminating_query_type[1]++;
731:       continue;
732:     }

734:     if (hash) {
735:       PetscBool found_box;

737:       /* allow for case that point is outside box - abort early */
738:       PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
739:       if (found_box) {
740:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
741:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
742:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
743:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
744:           DMPlexLocatePoint_Internal(dm, dim, point, boxCells[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:     } else {
755:       for (c = cStart; c < cEnd; ++c) {
756:         DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
757:         if (cell >= 0) {
758:           cells[p].rank = 0;
759:           cells[p].index = cell;
760:           numFound++;
761:           terminating_query_type[2]++;
762:           break;
763:         }
764:       }
765:     }
766:   }
767:   if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
768:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
769:     for (p = 0; p < numPoints; p++) {
770:       const PetscScalar *point = &a[p*bs];
771:       PetscReal          cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
772:       PetscInt           dbin[3] = {-1,-1,-1}, bin, cellOffset, d;

774:       if (cells[p].index < 0) {
775:         ++numFound;
776:         PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
777:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
778:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
779:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
780:           DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
781:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
782:           dist = DMPlex_NormD_Internal(dim, diff);
783:           if (dist < distMax) {
784:             for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
785:             cells[p].rank  = 0;
786:             cells[p].index = boxCells[c];
787:             distMax = dist;
788:             break;
789:           }
790:         }
791:       }
792:     }
793:   }
794:   /* This code is only be relevant when interfaced to parallel point location */
795:   /* Check for highest numbered proc that claims a point (do we care?) */
796:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
797:     PetscMalloc1(numFound,&found);
798:     for (p = 0, numFound = 0; p < numPoints; p++) {
799:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
800:         if (numFound < p) {
801:           cells[numFound] = cells[p];
802:         }
803:         found[numFound++] = p;
804:       }
805:     }
806:   }
807:   VecRestoreArray(v, &a);
808:   if (!reuse) {
809:     PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
810:   }
811:   PetscTime(&t1);
812:   if (hash) {
813:     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]);
814:   } else {
815:     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]);
816:   }
817:   PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
818:   return(0);
819: }

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

824:   Not collective

826:   Input Parameter:
827: . coords - The coordinates of a segment

829:   Output Parameters:
830: + coords - The new y-coordinate, and 0 for x
831: - R - The rotation which accomplishes the projection

833:   Level: developer

835: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
836: @*/
837: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
838: {
839:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
840:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
841:   const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;

844:   R[0] = c; R[1] = -s;
845:   R[2] = s; R[3] =  c;
846:   coords[0] = 0.0;
847:   coords[1] = r;
848:   return(0);
849: }

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

854:   Not collective

856:   Input Parameter:
857: . coords - The coordinates of a segment

859:   Output Parameters:
860: + coords - The new y-coordinate, and 0 for x and z
861: - R - The rotation which accomplishes the projection

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

865:   Level: developer

867: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
868: @*/
869: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
870: {
871:   PetscReal      x    = PetscRealPart(coords[3] - coords[0]);
872:   PetscReal      y    = PetscRealPart(coords[4] - coords[1]);
873:   PetscReal      z    = PetscRealPart(coords[5] - coords[2]);
874:   PetscReal      r    = PetscSqrtReal(x*x + y*y + z*z);
875:   PetscReal      rinv = 1. / r;

878:   x *= rinv; y *= rinv; z *= rinv;
879:   if (x > 0.) {
880:     PetscReal inv1pX   = 1./ (1. + x);

882:     R[0] = x; R[1] = -y;              R[2] = -z;
883:     R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] =     -y*z*inv1pX;
884:     R[6] = z; R[7] =     -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
885:   }
886:   else {
887:     PetscReal inv1mX   = 1./ (1. - x);

889:     R[0] = x; R[1] = z;               R[2] = y;
890:     R[3] = y; R[4] =     -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
891:     R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] =     -y*z*inv1mX;
892:   }
893:   coords[0] = 0.0;
894:   coords[1] = r;
895:   return(0);
896: }

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

901:   Not collective

903:   Input Parameter:
904: . coords - The coordinates of a segment

906:   Output Parameters:
907: + coords - The new y- and z-coordinates, and 0 for x
908: - R - The rotation which accomplishes the projection

910:   Level: developer

912: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
913: @*/
914: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
915: {
916:   PetscReal      x1[3],  x2[3], n[3], norm;
917:   PetscReal      x1p[3], x2p[3], xnp[3];
918:   PetscReal      sqrtz, alpha;
919:   const PetscInt dim = 3;
920:   PetscInt       d, e, p;

923:   /* 0) Calculate normal vector */
924:   for (d = 0; d < dim; ++d) {
925:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
926:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
927:   }
928:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
929:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
930:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
931:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
932:   n[0] /= norm;
933:   n[1] /= norm;
934:   n[2] /= norm;
935:   /* 1) Take the normal vector and rotate until it is \hat z

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

939:     R = /  alpha nx nz  alpha ny nz -1/alpha \
940:         | -alpha ny     alpha nx        0    |
941:         \     nx            ny         nz    /

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

1000:       tmp        = R[d*dim+e];
1001:       R[d*dim+e] = R[e*dim+d];
1002:       R[e*dim+d] = tmp;
1003:     }
1004:   }
1005:   return(0);
1006: }

1008: PETSC_UNUSED
1009: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
1010: {
1011:   /* Signed volume is 1/2 the determinant

1013:    |  1  1  1 |
1014:    | x0 x1 x2 |
1015:    | y0 y1 y2 |

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

1019:    | x1 x2 |
1020:    | y1 y2 |
1021:   */
1022:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
1023:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
1024:   PetscReal       M[4], detM;
1025:   M[0] = x1; M[1] = x2;
1026:   M[2] = y1; M[3] = y2;
1027:   DMPlex_Det2D_Internal(&detM, M);
1028:   *vol = 0.5*detM;
1029:   (void)PetscLogFlops(5.0);
1030: }

1032: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
1033: {
1034:   DMPlex_Det2D_Internal(vol, coords);
1035:   *vol *= 0.5;
1036: }

1038: PETSC_UNUSED
1039: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
1040: {
1041:   /* Signed volume is 1/6th of the determinant

1043:    |  1  1  1  1 |
1044:    | x0 x1 x2 x3 |
1045:    | y0 y1 y2 y3 |
1046:    | z0 z1 z2 z3 |

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

1050:    | x1 x2 x3 |
1051:    | y1 y2 y3 |
1052:    | z1 z2 z3 |
1053:   */
1054:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
1055:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
1056:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
1057:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1058:   PetscReal       M[9], detM;
1059:   M[0] = x1; M[1] = x2; M[2] = x3;
1060:   M[3] = y1; M[4] = y2; M[5] = y3;
1061:   M[6] = z1; M[7] = z2; M[8] = z3;
1062:   DMPlex_Det3D_Internal(&detM, M);
1063:   *vol = -onesixth*detM;
1064:   (void)PetscLogFlops(10.0);
1065: }

1067: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1068: {
1069:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1070:   DMPlex_Det3D_Internal(vol, coords);
1071:   *vol *= -onesixth;
1072: }

1074: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1075: {
1076:   PetscSection   coordSection;
1077:   Vec            coordinates;
1078:   const PetscScalar *coords;
1079:   PetscInt       dim, d, off;

1083:   DMGetCoordinatesLocal(dm, &coordinates);
1084:   DMGetCoordinateSection(dm, &coordSection);
1085:   PetscSectionGetDof(coordSection,e,&dim);
1086:   if (!dim) return(0);
1087:   PetscSectionGetOffset(coordSection,e,&off);
1088:   VecGetArrayRead(coordinates,&coords);
1089:   if (v0) {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[off + d]);}
1090:   VecRestoreArrayRead(coordinates,&coords);
1091:   *detJ = 1.;
1092:   if (J) {
1093:     for (d = 0; d < dim * dim; d++) J[d] = 0.;
1094:     for (d = 0; d < dim; d++) J[d * dim + d] = 1.;
1095:     if (invJ) {
1096:       for (d = 0; d < dim * dim; d++) invJ[d] = 0.;
1097:       for (d = 0; d < dim; d++) invJ[d * dim + d] = 1.;
1098:     }
1099:   }
1100:   return(0);
1101: }

1103: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1104: {
1105:   PetscSection   coordSection;
1106:   Vec            coordinates;
1107:   PetscScalar   *coords = NULL;
1108:   PetscInt       numCoords, d, pStart, pEnd, numSelfCoords = 0;

1112:   DMGetCoordinatesLocal(dm, &coordinates);
1113:   DMGetCoordinateSection(dm, &coordSection);
1114:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1115:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1116:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1117:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1118:   if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
1119:   *detJ = 0.0;
1120:   if (numCoords == 6) {
1121:     const PetscInt dim = 3;
1122:     PetscReal      R[9], J0;

1124:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1125:     DMPlexComputeProjection3Dto1D(coords, R);
1126:     if (J)    {
1127:       J0   = 0.5*PetscRealPart(coords[1]);
1128:       J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
1129:       J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
1130:       J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
1131:       DMPlex_Det3D_Internal(detJ, J);
1132:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1133:     }
1134:   } else if (numCoords == 4) {
1135:     const PetscInt dim = 2;
1136:     PetscReal      R[4], J0;

1138:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1139:     DMPlexComputeProjection2Dto1D(coords, R);
1140:     if (J)    {
1141:       J0   = 0.5*PetscRealPart(coords[1]);
1142:       J[0] = R[0]*J0; J[1] = R[1];
1143:       J[2] = R[2]*J0; J[3] = R[3];
1144:       DMPlex_Det2D_Internal(detJ, J);
1145:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1146:     }
1147:   } else if (numCoords == 2) {
1148:     const PetscInt dim = 1;

1150:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1151:     if (J)    {
1152:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
1153:       *detJ = J[0];
1154:       PetscLogFlops(2.0);
1155:       if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
1156:     }
1157:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
1158:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1159:   return(0);
1160: }

1162: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1163: {
1164:   PetscSection   coordSection;
1165:   Vec            coordinates;
1166:   PetscScalar   *coords = NULL;
1167:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1171:   DMGetCoordinatesLocal(dm, &coordinates);
1172:   DMGetCoordinateSection(dm, &coordSection);
1173:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1174:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1175:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1176:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1177:   *detJ = 0.0;
1178:   if (numCoords == 9) {
1179:     const PetscInt dim = 3;
1180:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1182:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1183:     DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1184:     if (J)    {
1185:       const PetscInt pdim = 2;

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

1208:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1209:     if (J)    {
1210:       for (d = 0; d < dim; d++) {
1211:         for (f = 0; f < dim; f++) {
1212:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
1213:         }
1214:       }
1215:       PetscLogFlops(8.0);
1216:       DMPlex_Det2D_Internal(detJ, J);
1217:     }
1218:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1219:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
1220:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1221:   return(0);
1222: }

1224: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1225: {
1226:   PetscSection   coordSection;
1227:   Vec            coordinates;
1228:   PetscScalar   *coords = NULL;
1229:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1233:   DMGetCoordinatesLocal(dm, &coordinates);
1234:   DMGetCoordinateSection(dm, &coordSection);
1235:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1236:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1237:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1238:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1239:   if (!Nq) {
1240:     *detJ = 0.0;
1241:     if (numCoords == 12) {
1242:       const PetscInt dim = 3;
1243:       PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1245:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1246:       DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1247:       if (J)    {
1248:         const PetscInt pdim = 2;

1250:         for (d = 0; d < pdim; d++) {
1251:           J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
1252:           J0[d*dim+1] = 0.5*(PetscRealPart(coords[2*pdim+d]) - PetscRealPart(coords[1*pdim+d]));
1253:         }
1254:         PetscLogFlops(8.0);
1255:         DMPlex_Det3D_Internal(detJ, J0);
1256:         for (d = 0; d < dim; d++) {
1257:           for (f = 0; f < dim; f++) {
1258:             J[d*dim+f] = 0.0;
1259:             for (g = 0; g < dim; g++) {
1260:               J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
1261:             }
1262:           }
1263:         }
1264:         PetscLogFlops(18.0);
1265:       }
1266:       if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1267:     } else if (numCoords == 8) {
1268:       const PetscInt dim = 2;

1270:       if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1271:       if (J)    {
1272:         for (d = 0; d < dim; d++) {
1273:           J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1274:           J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1275:         }
1276:         PetscLogFlops(8.0);
1277:         DMPlex_Det2D_Internal(detJ, J);
1278:       }
1279:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1280:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1281:   } else {
1282:     const PetscInt Nv = 4;
1283:     const PetscInt dimR = 2;
1284:     const PetscInt zToPlex[4] = {0, 1, 3, 2};
1285:     PetscReal zOrder[12];
1286:     PetscReal zCoeff[12];
1287:     PetscInt  i, j, k, l, dim;

1289:     if (numCoords == 12) {
1290:       dim = 3;
1291:     } else if (numCoords == 8) {
1292:       dim = 2;
1293:     } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
1294:     for (i = 0; i < Nv; i++) {
1295:       PetscInt zi = zToPlex[i];

1297:       for (j = 0; j < dim; j++) {
1298:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1299:       }
1300:     }
1301:     for (j = 0; j < dim; j++) {
1302:       zCoeff[dim * 0 + j] = 0.25 * (  zOrder[dim * 0 + j] + zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1303:       zCoeff[dim * 1 + j] = 0.25 * (- zOrder[dim * 0 + j] + zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1304:       zCoeff[dim * 2 + j] = 0.25 * (- zOrder[dim * 0 + j] - zOrder[dim * 1 + j] + zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1305:       zCoeff[dim * 3 + j] = 0.25 * (  zOrder[dim * 0 + j] - zOrder[dim * 1 + j] - zOrder[dim * 2 + j] + zOrder[dim * 3 + j]);
1306:     }
1307:     for (i = 0; i < Nq; i++) {
1308:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1];

1310:       if (v) {
1311:         PetscReal extPoint[4];

1313:         extPoint[0] = 1.;
1314:         extPoint[1] = xi;
1315:         extPoint[2] = eta;
1316:         extPoint[3] = xi * eta;
1317:         for (j = 0; j < dim; j++) {
1318:           PetscReal val = 0.;

1320:           for (k = 0; k < Nv; k++) {
1321:             val += extPoint[k] * zCoeff[dim * k + j];
1322:           }
1323:           v[i * dim + j] = val;
1324:         }
1325:       }
1326:       if (J) {
1327:         PetscReal extJ[8];

1329:         extJ[0] = 0.;
1330:         extJ[1] = 0.;
1331:         extJ[2] = 1.;
1332:         extJ[3] = 0.;
1333:         extJ[4] = 0.;
1334:         extJ[5] = 1.;
1335:         extJ[6] = eta;
1336:         extJ[7] = xi;
1337:         for (j = 0; j < dim; j++) {
1338:           for (k = 0; k < dimR; k++) {
1339:             PetscReal val = 0.;

1341:             for (l = 0; l < Nv; l++) {
1342:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1343:             }
1344:             J[i * dim * dim + dim * j + k] = val;
1345:           }
1346:         }
1347:         if (dim == 3) { /* put the cross product in the third component of the Jacobian */
1348:           PetscReal x, y, z;
1349:           PetscReal *iJ = &J[i * dim * dim];
1350:           PetscReal norm;

1352:           x = iJ[1 * dim + 0] * iJ[2 * dim + 1] - iJ[1 * dim + 1] * iJ[2 * dim + 0];
1353:           y = iJ[0 * dim + 1] * iJ[2 * dim + 0] - iJ[0 * dim + 0] * iJ[2 * dim + 1];
1354:           z = iJ[0 * dim + 0] * iJ[1 * dim + 1] - iJ[0 * dim + 1] * iJ[1 * dim + 0];
1355:           norm = PetscSqrtReal(x * x + y * y + z * z);
1356:           iJ[2] = x / norm;
1357:           iJ[5] = y / norm;
1358:           iJ[8] = z / norm;
1359:           DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1360:           if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1361:         } else {
1362:           DMPlex_Det2D_Internal(&detJ[i], &J[i * dim * dim]);
1363:           if (invJ) {DMPlex_Invert2D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1364:         }
1365:       }
1366:     }
1367:   }
1368:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1369:   return(0);
1370: }

1372: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1373: {
1374:   PetscSection   coordSection;
1375:   Vec            coordinates;
1376:   PetscScalar   *coords = NULL;
1377:   const PetscInt dim = 3;
1378:   PetscInt       d;

1382:   DMGetCoordinatesLocal(dm, &coordinates);
1383:   DMGetCoordinateSection(dm, &coordSection);
1384:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1385:   *detJ = 0.0;
1386:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1387:   if (J)    {
1388:     for (d = 0; d < dim; d++) {
1389:       /* I orient with outward face normals */
1390:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
1391:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1392:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1393:     }
1394:     PetscLogFlops(18.0);
1395:     DMPlex_Det3D_Internal(detJ, J);
1396:   }
1397:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1398:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1399:   return(0);
1400: }

1402: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1403: {
1404:   PetscSection   coordSection;
1405:   Vec            coordinates;
1406:   PetscScalar   *coords = NULL;
1407:   const PetscInt dim = 3;
1408:   PetscInt       d;

1412:   DMGetCoordinatesLocal(dm, &coordinates);
1413:   DMGetCoordinateSection(dm, &coordSection);
1414:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
1415:   if (!Nq) {
1416:     *detJ = 0.0;
1417:     if (v)   {for (d = 0; d < dim; d++) v[d] = PetscRealPart(coords[d]);}
1418:     if (J)    {
1419:       for (d = 0; d < dim; d++) {
1420:         J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
1421:         J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
1422:         J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
1423:       }
1424:       PetscLogFlops(18.0);
1425:       DMPlex_Det3D_Internal(detJ, J);
1426:     }
1427:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1428:   } else {
1429:     const PetscInt Nv = 8;
1430:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};
1431:     const PetscInt dim = 3;
1432:     const PetscInt dimR = 3;
1433:     PetscReal zOrder[24];
1434:     PetscReal zCoeff[24];
1435:     PetscInt  i, j, k, l;

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

1440:       for (j = 0; j < dim; j++) {
1441:         zOrder[dim * i + j] = PetscRealPart(coords[dim * zi + j]);
1442:       }
1443:     }
1444:     for (j = 0; j < dim; j++) {
1445:       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]);
1446:       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]);
1447:       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]);
1448:       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]);
1449:       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]);
1450:       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]);
1451:       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]);
1452:       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]);
1453:     }
1454:     for (i = 0; i < Nq; i++) {
1455:       PetscReal xi = points[dimR * i], eta = points[dimR * i + 1], theta = points[dimR * i + 2];

1457:       if (v) {
1458:         PetscReal extPoint[8];

1460:         extPoint[0] = 1.;
1461:         extPoint[1] = xi;
1462:         extPoint[2] = eta;
1463:         extPoint[3] = xi * eta;
1464:         extPoint[4] = theta;
1465:         extPoint[5] = theta * xi;
1466:         extPoint[6] = theta * eta;
1467:         extPoint[7] = theta * eta * xi;
1468:         for (j = 0; j < dim; j++) {
1469:           PetscReal val = 0.;

1471:           for (k = 0; k < Nv; k++) {
1472:             val += extPoint[k] * zCoeff[dim * k + j];
1473:           }
1474:           v[i * dim + j] = val;
1475:         }
1476:       }
1477:       if (J) {
1478:         PetscReal extJ[24];

1480:         extJ[0]  = 0.         ; extJ[1]  = 0.        ; extJ[2]  = 0.      ;
1481:         extJ[3]  = 1.         ; extJ[4]  = 0.        ; extJ[5]  = 0.      ;
1482:         extJ[6]  = 0.         ; extJ[7]  = 1.        ; extJ[8]  = 0.      ;
1483:         extJ[9]  = eta        ; extJ[10] = xi        ; extJ[11] = 0.      ;
1484:         extJ[12] = 0.         ; extJ[13] = 0.        ; extJ[14] = 1.      ;
1485:         extJ[15] = theta      ; extJ[16] = 0.        ; extJ[17] = xi      ;
1486:         extJ[18] = 0.         ; extJ[19] = theta     ; extJ[20] = eta     ;
1487:         extJ[21] = theta * eta; extJ[22] = theta * xi; extJ[23] = eta * xi;

1489:         for (j = 0; j < dim; j++) {
1490:           for (k = 0; k < dimR; k++) {
1491:             PetscReal val = 0.;

1493:             for (l = 0; l < Nv; l++) {
1494:               val += zCoeff[dim * l + j] * extJ[dimR * l + k];
1495:             }
1496:             J[i * dim * dim + dim * j + k] = val;
1497:           }
1498:         }
1499:         DMPlex_Det3D_Internal(&detJ[i], &J[i * dim * dim]);
1500:         if (invJ) {DMPlex_Invert3D_Internal(&invJ[i * dim * dim], &J[i * dim * dim], detJ[i]);}
1501:       }
1502:     }
1503:   }
1504:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
1505:   return(0);
1506: }

1508: static PetscErrorCode DMPlexComputeCellGeometryFEM_Implicit(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1509: {
1510:   PetscInt        depth, dim, coordDim, coneSize, i;
1511:   PetscInt        Nq = 0;
1512:   const PetscReal *points = NULL;
1513:   DMLabel         depthLabel;
1514:   PetscReal       xi0[3] = {-1.,-1.,-1.}, v0[3], J0[9], detJ0;
1515:   PetscBool       isAffine = PETSC_TRUE;
1516:   PetscErrorCode  ierr;

1519:   DMPlexGetDepth(dm, &depth);
1520:   DMPlexGetConeSize(dm, cell, &coneSize);
1521:   DMPlexGetDepthLabel(dm, &depthLabel);
1522:   DMLabelGetValue(depthLabel, cell, &dim);
1523:   if (depth == 1 && dim == 1) {
1524:     DMGetDimension(dm, &dim);
1525:   }
1526:   DMGetCoordinateDim(dm, &coordDim);
1527:   if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1528:   if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1529:   switch (dim) {
1530:   case 0:
1531:     DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1532:     isAffine = PETSC_FALSE;
1533:     break;
1534:   case 1:
1535:     if (Nq) {
1536:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1537:     } else {
1538:       DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);
1539:     }
1540:     break;
1541:   case 2:
1542:     switch (coneSize) {
1543:     case 3:
1544:       if (Nq) {
1545:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1546:       } else {
1547:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);
1548:       }
1549:       break;
1550:     case 4:
1551:       DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1552:       isAffine = PETSC_FALSE;
1553:       break;
1554:     default:
1555:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1556:     }
1557:     break;
1558:   case 3:
1559:     switch (coneSize) {
1560:     case 4:
1561:       if (Nq) {
1562:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1563:       } else {
1564:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);
1565:       }
1566:       break;
1567:     case 6: /* Faces */
1568:     case 8: /* Vertices */
1569:       DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1570:       isAffine = PETSC_FALSE;
1571:       break;
1572:     default:
1573:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1574:     }
1575:     break;
1576:   default:
1577:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1578:   }
1579:   if (isAffine && Nq) {
1580:     if (v) {
1581:       for (i = 0; i < Nq; i++) {
1582:         CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1583:       }
1584:     }
1585:     if (detJ) {
1586:       for (i = 0; i < Nq; i++) {
1587:         detJ[i] = detJ0;
1588:       }
1589:     }
1590:     if (J) {
1591:       PetscInt k;

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

1596:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1597:           J[k] = J0[j];
1598:         }
1599:       }
1600:     }
1601:     if (invJ) {
1602:       PetscInt k;
1603:       switch (coordDim) {
1604:       case 0:
1605:         break;
1606:       case 1:
1607:         invJ[0] = 1./J0[0];
1608:         break;
1609:       case 2:
1610:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1611:         break;
1612:       case 3:
1613:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1614:         break;
1615:       }
1616:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1617:         PetscInt j;

1619:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1620:           invJ[k] = invJ[j];
1621:         }
1622:       }
1623:     }
1624:   }
1625:   return(0);
1626: }

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

1631:   Collective on dm

1633:   Input Arguments:
1634: + dm   - the DM
1635: - cell - the cell

1637:   Output Arguments:
1638: + v0   - the translation part of this affine transform
1639: . J    - the Jacobian of the transform from the reference element
1640: . invJ - the inverse of the Jacobian
1641: - detJ - the Jacobian determinant

1643:   Level: advanced

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

1649: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1650: @*/
1651: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1652: {

1656:   DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1657:   return(0);
1658: }

1660: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1661: {
1662:   PetscQuadrature  feQuad;
1663:   PetscSection     coordSection;
1664:   Vec              coordinates;
1665:   PetscScalar     *coords = NULL;
1666:   const PetscReal *quadPoints;
1667:   PetscReal       *basisDer, *basis, detJt;
1668:   PetscInt         dim, cdim, pdim, qdim, Nq, numCoords, q;
1669:   PetscErrorCode   ierr;

1672:   DMGetCoordinatesLocal(dm, &coordinates);
1673:   DMGetCoordinateSection(dm, &coordSection);
1674:   DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1675:   DMGetDimension(dm, &dim);
1676:   DMGetCoordinateDim(dm, &cdim);
1677:   if (!quad) { /* use the first point of the first functional of the dual space */
1678:     PetscDualSpace dsp;

1680:     PetscFEGetDualSpace(fe, &dsp);
1681:     PetscDualSpaceGetFunctional(dsp, 0, &quad);
1682:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1683:     Nq = 1;
1684:   } else {
1685:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1686:   }
1687:   PetscFEGetDimension(fe, &pdim);
1688:   PetscFEGetQuadrature(fe, &feQuad);
1689:   if (feQuad == quad) {
1690:     PetscFEGetDefaultTabulation(fe, &basis, J ? &basisDer : NULL, NULL);
1691:     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);
1692:   } else {
1693:     PetscFEGetTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1694:   }
1695:   if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1696:   if (v) {
1697:     PetscArrayzero(v, Nq*cdim);
1698:     for (q = 0; q < Nq; ++q) {
1699:       PetscInt i, k;

1701:       for (k = 0; k < pdim; ++k)
1702:         for (i = 0; i < cdim; ++i)
1703:           v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1704:       PetscLogFlops(2.0*pdim*cdim);
1705:     }
1706:   }
1707:   if (J) {
1708:     PetscArrayzero(J, Nq*cdim*cdim);
1709:     for (q = 0; q < Nq; ++q) {
1710:       PetscInt i, j, k, c, r;

1712:       /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1713:       for (k = 0; k < pdim; ++k)
1714:         for (j = 0; j < dim; ++j)
1715:           for (i = 0; i < cdim; ++i)
1716:             J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1717:       PetscLogFlops(2.0*pdim*dim*cdim);
1718:       if (cdim > dim) {
1719:         for (c = dim; c < cdim; ++c)
1720:           for (r = 0; r < cdim; ++r)
1721:             J[r*cdim+c] = r == c ? 1.0 : 0.0;
1722:       }
1723:       if (!detJ && !invJ) continue;
1724:       detJt = 0.;
1725:       switch (cdim) {
1726:       case 3:
1727:         DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1728:         if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1729:         break;
1730:       case 2:
1731:         DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1732:         if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1733:         break;
1734:       case 1:
1735:         detJt = J[q*cdim*dim];
1736:         if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1737:       }
1738:       if (detJ) detJ[q] = detJt;
1739:     }
1740:   }
1741:   else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1742:   if (feQuad != quad) {
1743:     PetscFERestoreTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1744:   }
1745:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1746:   return(0);
1747: }

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

1752:   Collective on dm

1754:   Input Arguments:
1755: + dm   - the DM
1756: . cell - the cell
1757: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If quad == NULL, geometry will be
1758:          evaluated at the first vertex of the reference element

1760:   Output Arguments:
1761: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1762: . J    - the Jacobian of the transform from the reference element at each quadrature point
1763: . invJ - the inverse of the Jacobian at each quadrature point
1764: - detJ - the Jacobian determinant at each quadrature point

1766:   Level: advanced

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

1772: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1773: @*/
1774: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1775: {
1776:   DM             cdm;
1777:   PetscFE        fe = NULL;

1782:   DMGetCoordinateDM(dm, &cdm);
1783:   if (cdm) {
1784:     PetscClassId id;
1785:     PetscInt     numFields;
1786:     PetscDS      prob;
1787:     PetscObject  disc;

1789:     DMGetNumFields(cdm, &numFields);
1790:     if (numFields) {
1791:       DMGetDS(cdm, &prob);
1792:       PetscDSGetDiscretization(prob,0,&disc);
1793:       PetscObjectGetClassId(disc,&id);
1794:       if (id == PETSCFE_CLASSID) {
1795:         fe = (PetscFE) disc;
1796:       }
1797:     }
1798:   }
1799:   if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1800:   else     {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1801:   return(0);
1802: }

1804: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1805: {
1806:   PetscSection   coordSection;
1807:   Vec            coordinates;
1808:   PetscScalar   *coords = NULL;
1809:   PetscScalar    tmp[2];
1810:   PetscInt       coordSize, d;

1814:   DMGetCoordinatesLocal(dm, &coordinates);
1815:   DMGetCoordinateSection(dm, &coordSection);
1816:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1817:   DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1818:   if (centroid) {
1819:     for (d = 0; d < dim; ++d) centroid[d] = 0.5*PetscRealPart(coords[d] + tmp[d]);
1820:   }
1821:   if (normal) {
1822:     PetscReal norm;

1824:     if (dim != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support 2D edges right now");
1825:     normal[0]  = -PetscRealPart(coords[1] - tmp[1]);
1826:     normal[1]  =  PetscRealPart(coords[0] - tmp[0]);
1827:     norm       = DMPlex_NormD_Internal(dim, normal);
1828:     for (d = 0; d < dim; ++d) normal[d] /= norm;
1829:   }
1830:   if (vol) {
1831:     *vol = 0.0;
1832:     for (d = 0; d < dim; ++d) *vol += PetscSqr(PetscRealPart(coords[d] - tmp[d]));
1833:     *vol = PetscSqrtReal(*vol);
1834:   }
1835:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1836:   return(0);
1837: }

1839: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1840: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1841: {
1842:   DMLabel        depth;
1843:   PetscSection   coordSection;
1844:   Vec            coordinates;
1845:   PetscScalar   *coords = NULL;
1846:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1847:   PetscBool      isHybrid = PETSC_FALSE;
1848:   PetscInt       fv[4] = {0, 1, 2, 3};
1849:   PetscInt       pEndInterior[4], cdepth, tdim = 2, coordSize, numCorners, p, d, e;

1853:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1854:   DMPlexGetDepthLabel(dm, &depth);
1855:   DMLabelGetValue(depth, cell, &cdepth);
1856:   DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1857:   if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;
1858:   DMGetCoordinatesLocal(dm, &coordinates);
1859:   DMPlexGetConeSize(dm, cell, &numCorners);
1860:   DMGetCoordinateSection(dm, &coordSection);
1861:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1862:   DMGetCoordinateDim(dm, &dim);
1863:   /* Side faces for hybrid cells are are stored as tensor products */
1864:   if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}

1866:   if (dim > 2 && centroid) {
1867:     v0[0] = PetscRealPart(coords[0]);
1868:     v0[1] = PetscRealPart(coords[1]);
1869:     v0[2] = PetscRealPart(coords[2]);
1870:   }
1871:   if (normal) {
1872:     if (dim > 2) {
1873:       const PetscReal x0 = PetscRealPart(coords[dim*fv[1]+0] - coords[0]), x1 = PetscRealPart(coords[dim*fv[2]+0] - coords[0]);
1874:       const PetscReal y0 = PetscRealPart(coords[dim*fv[1]+1] - coords[1]), y1 = PetscRealPart(coords[dim*fv[2]+1] - coords[1]);
1875:       const PetscReal z0 = PetscRealPart(coords[dim*fv[1]+2] - coords[2]), z1 = PetscRealPart(coords[dim*fv[2]+2] - coords[2]);
1876:       PetscReal       norm;

1878:       normal[0] = y0*z1 - z0*y1;
1879:       normal[1] = z0*x1 - x0*z1;
1880:       normal[2] = x0*y1 - y0*x1;
1881:       norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1882:       normal[0] /= norm;
1883:       normal[1] /= norm;
1884:       normal[2] /= norm;
1885:     } else {
1886:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
1887:     }
1888:   }
1889:   if (dim == 3) {DMPlexComputeProjection3Dto2D(coordSize, coords, R);}
1890:   for (p = 0; p < numCorners; ++p) {
1891:     const PetscInt pi  = p < 4 ? fv[p] : p;
1892:     const PetscInt pin = p < 3 ? fv[(p+1)%numCorners] : (p+1)%numCorners;
1893:     /* Need to do this copy to get types right */
1894:     for (d = 0; d < tdim; ++d) {
1895:       ctmp[d]      = PetscRealPart(coords[pi*tdim+d]);
1896:       ctmp[tdim+d] = PetscRealPart(coords[pin*tdim+d]);
1897:     }
1898:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1899:     vsum += vtmp;
1900:     for (d = 0; d < tdim; ++d) {
1901:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1902:     }
1903:   }
1904:   for (d = 0; d < tdim; ++d) {
1905:     csum[d] /= (tdim+1)*vsum;
1906:   }
1907:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1908:   if (vol) *vol = PetscAbsReal(vsum);
1909:   if (centroid) {
1910:     if (dim > 2) {
1911:       for (d = 0; d < dim; ++d) {
1912:         centroid[d] = v0[d];
1913:         for (e = 0; e < dim; ++e) {
1914:           centroid[d] += R[d*dim+e]*csum[e];
1915:         }
1916:       }
1917:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1918:   }
1919:   return(0);
1920: }

1922: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1923: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1924: {
1925:   DMLabel         depth;
1926:   PetscSection    coordSection;
1927:   Vec             coordinates;
1928:   PetscScalar    *coords = NULL;
1929:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
1930:   const PetscInt *faces, *facesO;
1931:   PetscBool       isHybrid = PETSC_FALSE;
1932:   PetscInt        pEndInterior[4], cdepth, numFaces, f, coordSize, numCorners, p, d;
1933:   PetscErrorCode  ierr;

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

1943:   DMGetCoordinatesLocal(dm, &coordinates);
1944:   DMGetCoordinateSection(dm, &coordSection);

1946:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1947:   DMPlexGetConeSize(dm, cell, &numFaces);
1948:   DMPlexGetCone(dm, cell, &faces);
1949:   DMPlexGetConeOrientation(dm, cell, &facesO);
1950:   for (f = 0; f < numFaces; ++f) {
1951:     PetscBool flip = isHybrid && f == 0 ? PETSC_TRUE : PETSC_FALSE; /* The first hybrid face is reversed */

1953:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1954:     numCorners = coordSize/dim;
1955:     switch (numCorners) {
1956:     case 3:
1957:       for (d = 0; d < dim; ++d) {
1958:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1959:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1960:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1961:       }
1962:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1963:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1964:       vsum += vtmp;
1965:       if (centroid) {           /* Centroid of OABC = (a+b+c)/4 */
1966:         for (d = 0; d < dim; ++d) {
1967:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1968:         }
1969:       }
1970:       break;
1971:     case 4:
1972:     {
1973:       PetscInt fv[4] = {0, 1, 2, 3};

1975:       /* Side faces for hybrid cells are are stored as tensor products */
1976:       if (isHybrid && f > 1) {fv[2] = 3; fv[3] = 2;}
1977:       /* DO FOR PYRAMID */
1978:       /* First tet */
1979:       for (d = 0; d < dim; ++d) {
1980:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[0]*dim+d]);
1981:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1982:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1983:       }
1984:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1985:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
1986:       vsum += vtmp;
1987:       if (centroid) {
1988:         for (d = 0; d < dim; ++d) {
1989:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1990:         }
1991:       }
1992:       /* Second tet */
1993:       for (d = 0; d < dim; ++d) {
1994:         coordsTmp[0*dim+d] = PetscRealPart(coords[fv[1]*dim+d]);
1995:         coordsTmp[1*dim+d] = PetscRealPart(coords[fv[2]*dim+d]);
1996:         coordsTmp[2*dim+d] = PetscRealPart(coords[fv[3]*dim+d]);
1997:       }
1998:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1999:       if (facesO[f] < 0 || flip) vtmp = -vtmp;
2000:       vsum += vtmp;
2001:       if (centroid) {
2002:         for (d = 0; d < dim; ++d) {
2003:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
2004:         }
2005:       }
2006:       break;
2007:     }
2008:     default:
2009:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
2010:     }
2011:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
2012:   }
2013:   if (vol)     *vol = PetscAbsReal(vsum);
2014:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
2015:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
2016:   return(0);
2017: }

2019: /*@C
2020:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

2022:   Collective on dm

2024:   Input Arguments:
2025: + dm   - the DM
2026: - cell - the cell

2028:   Output Arguments:
2029: + volume   - the cell volume
2030: . centroid - the cell centroid
2031: - normal - the cell normal, if appropriate

2033:   Level: advanced

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

2039: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
2040: @*/
2041: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
2042: {
2043:   PetscInt       depth, dim;

2047:   DMPlexGetDepth(dm, &depth);
2048:   DMGetDimension(dm, &dim);
2049:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
2050:   /* We need to keep a pointer to the depth label */
2051:   DMGetLabelValue(dm, "depth", cell, &depth);
2052:   /* Cone size is now the number of faces */
2053:   switch (depth) {
2054:   case 1:
2055:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
2056:     break;
2057:   case 2:
2058:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
2059:     break;
2060:   case 3:
2061:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
2062:     break;
2063:   default:
2064:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2065:   }
2066:   return(0);
2067: }

2069: /*@
2070:   DMPlexComputeGeometryFEM - Precompute cell geometry for the entire mesh

2072:   Collective on dm

2074:   Input Parameter:
2075: . dm - The DMPlex

2077:   Output Parameter:
2078: . cellgeom - A vector with the cell geometry data for each cell

2080:   Level: beginner

2082: @*/
2083: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
2084: {
2085:   DM             dmCell;
2086:   Vec            coordinates;
2087:   PetscSection   coordSection, sectionCell;
2088:   PetscScalar   *cgeom;
2089:   PetscInt       cStart, cEnd, cMax, c;

2093:   DMClone(dm, &dmCell);
2094:   DMGetCoordinateSection(dm, &coordSection);
2095:   DMGetCoordinatesLocal(dm, &coordinates);
2096:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2097:   DMSetCoordinatesLocal(dmCell, coordinates);
2098:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2099:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2100:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
2101:   cEnd = cMax < 0 ? cEnd : cMax;
2102:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2103:   /* TODO This needs to be multiplied by Nq for non-affine */
2104:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFEGeom))/sizeof(PetscScalar)));}
2105:   PetscSectionSetUp(sectionCell);
2106:   DMSetLocalSection(dmCell, sectionCell);
2107:   PetscSectionDestroy(&sectionCell);
2108:   DMCreateLocalVector(dmCell, cellgeom);
2109:   VecGetArray(*cellgeom, &cgeom);
2110:   for (c = cStart; c < cEnd; ++c) {
2111:     PetscFEGeom *cg;

2113:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2114:     PetscArrayzero(cg, 1);
2115:     DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v, cg->J, cg->invJ, cg->detJ);
2116:     if (*cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D", (double) *cg->detJ, c);
2117:   }
2118:   return(0);
2119: }

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

2124:   Input Parameter:
2125: . dm - The DM

2127:   Output Parameters:
2128: + cellgeom - A Vec of PetscFVCellGeom data
2129: - facegeom - A Vec of PetscFVFaceGeom data

2131:   Level: developer

2133: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
2134: @*/
2135: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
2136: {
2137:   DM             dmFace, dmCell;
2138:   DMLabel        ghostLabel;
2139:   PetscSection   sectionFace, sectionCell;
2140:   PetscSection   coordSection;
2141:   Vec            coordinates;
2142:   PetscScalar   *fgeom, *cgeom;
2143:   PetscReal      minradius, gminradius;
2144:   PetscInt       dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

2148:   DMGetDimension(dm, &dim);
2149:   DMGetCoordinateSection(dm, &coordSection);
2150:   DMGetCoordinatesLocal(dm, &coordinates);
2151:   /* Make cell centroids and volumes */
2152:   DMClone(dm, &dmCell);
2153:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
2154:   DMSetCoordinatesLocal(dmCell, coordinates);
2155:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
2156:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2157:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2158:   PetscSectionSetChart(sectionCell, cStart, cEnd);
2159:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
2160:   PetscSectionSetUp(sectionCell);
2161:   DMSetLocalSection(dmCell, sectionCell);
2162:   PetscSectionDestroy(&sectionCell);
2163:   DMCreateLocalVector(dmCell, cellgeom);
2164:   if (cEndInterior < 0) cEndInterior = cEnd;
2165:   VecGetArray(*cellgeom, &cgeom);
2166:   for (c = cStart; c < cEndInterior; ++c) {
2167:     PetscFVCellGeom *cg;

2169:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
2170:     PetscArrayzero(cg, 1);
2171:     DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
2172:   }
2173:   /* Compute face normals and minimum cell radius */
2174:   DMClone(dm, &dmFace);
2175:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionFace);
2176:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2177:   PetscSectionSetChart(sectionFace, fStart, fEnd);
2178:   for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
2179:   PetscSectionSetUp(sectionFace);
2180:   DMSetLocalSection(dmFace, sectionFace);
2181:   PetscSectionDestroy(&sectionFace);
2182:   DMCreateLocalVector(dmFace, facegeom);
2183:   VecGetArray(*facegeom, &fgeom);
2184:   DMGetLabel(dm, "ghost", &ghostLabel);
2185:   minradius = PETSC_MAX_REAL;
2186:   for (f = fStart; f < fEnd; ++f) {
2187:     PetscFVFaceGeom *fg;
2188:     PetscReal        area;
2189:     PetscInt         ghost = -1, d, numChildren;

2191:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2192:     DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
2193:     if (ghost >= 0 || numChildren) continue;
2194:     DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
2195:     DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
2196:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
2197:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
2198:     {
2199:       PetscFVCellGeom *cL, *cR;
2200:       PetscInt         ncells;
2201:       const PetscInt  *cells;
2202:       PetscReal       *lcentroid, *rcentroid;
2203:       PetscReal        l[3], r[3], v[3];

2205:       DMPlexGetSupport(dm, f, &cells);
2206:       DMPlexGetSupportSize(dm, f, &ncells);
2207:       DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
2208:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
2209:       if (ncells > 1) {
2210:         DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
2211:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
2212:       }
2213:       else {
2214:         rcentroid = fg->centroid;
2215:       }
2216:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
2217:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
2218:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
2219:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
2220:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
2221:       }
2222:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
2223:         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]);
2224:         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]);
2225:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
2226:       }
2227:       if (cells[0] < cEndInterior) {
2228:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
2229:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2230:       }
2231:       if (ncells > 1 && cells[1] < cEndInterior) {
2232:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
2233:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
2234:       }
2235:     }
2236:   }
2237:   MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
2238:   DMPlexSetMinRadius(dm, gminradius);
2239:   /* Compute centroids of ghost cells */
2240:   for (c = cEndInterior; c < cEnd; ++c) {
2241:     PetscFVFaceGeom *fg;
2242:     const PetscInt  *cone,    *support;
2243:     PetscInt         coneSize, supportSize, s;

2245:     DMPlexGetConeSize(dmCell, c, &coneSize);
2246:     if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
2247:     DMPlexGetCone(dmCell, c, &cone);
2248:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
2249:     if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
2250:     DMPlexGetSupport(dmCell, cone[0], &support);
2251:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
2252:     for (s = 0; s < 2; ++s) {
2253:       /* Reflect ghost centroid across plane of face */
2254:       if (support[s] == c) {
2255:         PetscFVCellGeom       *ci;
2256:         PetscFVCellGeom       *cg;
2257:         PetscReal              c2f[3], a;

2259:         DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
2260:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
2261:         a    = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
2262:         DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
2263:         DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
2264:         cg->volume = ci->volume;
2265:       }
2266:     }
2267:   }
2268:   VecRestoreArray(*facegeom, &fgeom);
2269:   VecRestoreArray(*cellgeom, &cgeom);
2270:   DMDestroy(&dmCell);
2271:   DMDestroy(&dmFace);
2272:   return(0);
2273: }

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

2278:   Not collective

2280:   Input Argument:
2281: . dm - the DM

2283:   Output Argument:
2284: . minradius - the minium cell radius

2286:   Level: developer

2288: .seealso: DMGetCoordinates()
2289: @*/
2290: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2291: {
2295:   *minradius = ((DM_Plex*) dm->data)->minradius;
2296:   return(0);
2297: }

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

2302:   Logically collective

2304:   Input Arguments:
2305: + dm - the DM
2306: - minradius - the minium cell radius

2308:   Level: developer

2310: .seealso: DMSetCoordinates()
2311: @*/
2312: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2313: {
2316:   ((DM_Plex*) dm->data)->minradius = minradius;
2317:   return(0);
2318: }

2320: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2321: {
2322:   DMLabel        ghostLabel;
2323:   PetscScalar   *dx, *grad, **gref;
2324:   PetscInt       dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

2328:   DMGetDimension(dm, &dim);
2329:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2330:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2331:   DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
2332:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2333:   DMGetLabel(dm, "ghost", &ghostLabel);
2334:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2335:   for (c = cStart; c < cEndInterior; c++) {
2336:     const PetscInt        *faces;
2337:     PetscInt               numFaces, usedFaces, f, d;
2338:     PetscFVCellGeom        *cg;
2339:     PetscBool              boundary;
2340:     PetscInt               ghost;

2342:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2343:     DMPlexGetConeSize(dm, c, &numFaces);
2344:     DMPlexGetCone(dm, c, &faces);
2345:     if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
2346:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2347:       PetscFVCellGeom       *cg1;
2348:       PetscFVFaceGeom       *fg;
2349:       const PetscInt        *fcells;
2350:       PetscInt               ncell, side;

2352:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2353:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2354:       if ((ghost >= 0) || boundary) continue;
2355:       DMPlexGetSupport(dm, faces[f], &fcells);
2356:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
2357:       ncell = fcells[!side];    /* the neighbor */
2358:       DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
2359:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2360:       for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
2361:       gref[usedFaces++] = fg->grad[side];  /* Gradient reconstruction term will go here */
2362:     }
2363:     if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
2364:     PetscFVComputeGradient(fvm, usedFaces, dx, grad);
2365:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
2366:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
2367:       DMIsBoundaryPoint(dm, faces[f], &boundary);
2368:       if ((ghost >= 0) || boundary) continue;
2369:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
2370:       ++usedFaces;
2371:     }
2372:   }
2373:   PetscFree3(dx, grad, gref);
2374:   return(0);
2375: }

2377: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
2378: {
2379:   DMLabel        ghostLabel;
2380:   PetscScalar   *dx, *grad, **gref;
2381:   PetscInt       dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
2382:   PetscSection   neighSec;
2383:   PetscInt     (*neighbors)[2];
2384:   PetscInt      *counter;

2388:   DMGetDimension(dm, &dim);
2389:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2390:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2391:   if (cEndInterior < 0) cEndInterior = cEnd;
2392:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
2393:   PetscSectionSetChart(neighSec,cStart,cEndInterior);
2394:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
2395:   DMGetLabel(dm, "ghost", &ghostLabel);
2396:   for (f = fStart; f < fEnd; f++) {
2397:     const PetscInt        *fcells;
2398:     PetscBool              boundary;
2399:     PetscInt               ghost = -1;
2400:     PetscInt               numChildren, numCells, c;

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

2412:         if (cell >= cStart && cell < cEndInterior) {
2413:           PetscSectionAddDof(neighSec,cell,1);
2414:         }
2415:       }
2416:     }
2417:   }
2418:   PetscSectionSetUp(neighSec);
2419:   PetscSectionGetMaxDof(neighSec,&maxNumFaces);
2420:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
2421:   nStart = 0;
2422:   PetscSectionGetStorageSize(neighSec,&nEnd);
2423:   PetscMalloc1((nEnd-nStart),&neighbors);
2424:   PetscCalloc1((cEndInterior-cStart),&counter);
2425:   for (f = fStart; f < fEnd; f++) {
2426:     const PetscInt        *fcells;
2427:     PetscBool              boundary;
2428:     PetscInt               ghost = -1;
2429:     PetscInt               numChildren, numCells, c;

2431:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
2432:     DMIsBoundaryPoint(dm, f, &boundary);
2433:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
2434:     if ((ghost >= 0) || boundary || numChildren) continue;
2435:     DMPlexGetSupportSize(dm, f, &numCells);
2436:     if (numCells == 2) {
2437:       DMPlexGetSupport(dm, f, &fcells);
2438:       for (c = 0; c < 2; c++) {
2439:         PetscInt cell = fcells[c], off;

2441:         if (cell >= cStart && cell < cEndInterior) {
2442:           PetscSectionGetOffset(neighSec,cell,&off);
2443:           off += counter[cell - cStart]++;
2444:           neighbors[off][0] = f;
2445:           neighbors[off][1] = fcells[1 - c];
2446:         }
2447:       }
2448:     }
2449:   }
2450:   PetscFree(counter);
2451:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
2452:   for (c = cStart; c < cEndInterior; c++) {
2453:     PetscInt               numFaces, f, d, off, ghost = -1;
2454:     PetscFVCellGeom        *cg;

2456:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2457:     PetscSectionGetDof(neighSec, c, &numFaces);
2458:     PetscSectionGetOffset(neighSec, c, &off);
2459:     if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2460:     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);
2461:     for (f = 0; f < numFaces; ++f) {
2462:       PetscFVCellGeom       *cg1;
2463:       PetscFVFaceGeom       *fg;
2464:       const PetscInt        *fcells;
2465:       PetscInt               ncell, side, nface;

2467:       nface = neighbors[off + f][0];
2468:       ncell = neighbors[off + f][1];
2469:       DMPlexGetSupport(dm,nface,&fcells);
2470:       side  = (c != fcells[0]);
2471:       DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
2472:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
2473:       for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
2474:       gref[f] = fg->grad[side];  /* Gradient reconstruction term will go here */
2475:     }
2476:     PetscFVComputeGradient(fvm, numFaces, dx, grad);
2477:     for (f = 0; f < numFaces; ++f) {
2478:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
2479:     }
2480:   }
2481:   PetscFree3(dx, grad, gref);
2482:   PetscSectionDestroy(&neighSec);
2483:   PetscFree(neighbors);
2484:   return(0);
2485: }

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

2490:   Collective on dm

2492:   Input Arguments:
2493: + dm  - The DM
2494: . fvm - The PetscFV
2495: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2496: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()

2498:   Output Parameters:
2499: + faceGeometry - The geometric factors for gradient calculation are inserted
2500: - dmGrad - The DM describing the layout of gradient data

2502:   Level: developer

2504: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
2505: @*/
2506: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
2507: {
2508:   DM             dmFace, dmCell;
2509:   PetscScalar   *fgeom, *cgeom;
2510:   PetscSection   sectionGrad, parentSection;
2511:   PetscInt       dim, pdim, cStart, cEnd, cEndInterior, c;

2515:   DMGetDimension(dm, &dim);
2516:   PetscFVGetNumComponents(fvm, &pdim);
2517:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
2518:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
2519:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
2520:   VecGetDM(faceGeometry, &dmFace);
2521:   VecGetDM(cellGeometry, &dmCell);
2522:   VecGetArray(faceGeometry, &fgeom);
2523:   VecGetArray(cellGeometry, &cgeom);
2524:   DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
2525:   if (!parentSection) {
2526:     BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2527:   } else {
2528:     BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
2529:   }
2530:   VecRestoreArray(faceGeometry, &fgeom);
2531:   VecRestoreArray(cellGeometry, &cgeom);
2532:   /* Create storage for gradients */
2533:   DMClone(dm, dmGrad);
2534:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionGrad);
2535:   PetscSectionSetChart(sectionGrad, cStart, cEnd);
2536:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
2537:   PetscSectionSetUp(sectionGrad);
2538:   DMSetLocalSection(*dmGrad, sectionGrad);
2539:   PetscSectionDestroy(&sectionGrad);
2540:   return(0);
2541: }

2543: /*@
2544:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2546:   Collective on dm

2548:   Input Arguments:
2549: + dm  - The DM
2550: - fvm - The PetscFV

2552:   Output Parameters:
2553: + cellGeometry - The cell geometry
2554: . faceGeometry - The face geometry
2555: - dmGrad       - The gradient matrices

2557:   Level: developer

2559: .seealso: DMPlexComputeGeometryFVM()
2560: @*/
2561: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2562: {
2563:   PetscObject    cellgeomobj, facegeomobj;

2567:   PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2568:   if (!cellgeomobj) {
2569:     Vec cellgeomInt, facegeomInt;

2571:     DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2572:     PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2573:     PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2574:     VecDestroy(&cellgeomInt);
2575:     VecDestroy(&facegeomInt);
2576:     PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2577:   }
2578:   PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2579:   if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2580:   if (facegeom) *facegeom = (Vec) facegeomobj;
2581:   if (gradDM) {
2582:     PetscObject gradobj;
2583:     PetscBool   computeGradients;

2585:     PetscFVGetComputeGradients(fv,&computeGradients);
2586:     if (!computeGradients) {
2587:       *gradDM = NULL;
2588:       return(0);
2589:     }
2590:     PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2591:     if (!gradobj) {
2592:       DM dmGradInt;

2594:       DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2595:       PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2596:       DMDestroy(&dmGradInt);
2597:       PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2598:     }
2599:     *gradDM = (DM) gradobj;
2600:   }
2601:   return(0);
2602: }

2604: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work,  PetscReal *resNeg, PetscReal *guess)
2605: {
2606:   PetscInt l, m;

2609:   if (dimC == dimR && dimR <= 3) {
2610:     /* invert Jacobian, multiply */
2611:     PetscScalar det, idet;

2613:     switch (dimR) {
2614:     case 1:
2615:       invJ[0] = 1./ J[0];
2616:       break;
2617:     case 2:
2618:       det = J[0] * J[3] - J[1] * J[2];
2619:       idet = 1./det;
2620:       invJ[0] =  J[3] * idet;
2621:       invJ[1] = -J[1] * idet;
2622:       invJ[2] = -J[2] * idet;
2623:       invJ[3] =  J[0] * idet;
2624:       break;
2625:     case 3:
2626:       {
2627:         invJ[0] = J[4] * J[8] - J[5] * J[7];
2628:         invJ[1] = J[2] * J[7] - J[1] * J[8];
2629:         invJ[2] = J[1] * J[5] - J[2] * J[4];
2630:         det = invJ[0] * J[0] + invJ[1] * J[3] + invJ[2] * J[6];
2631:         idet = 1./det;
2632:         invJ[0] *= idet;
2633:         invJ[1] *= idet;
2634:         invJ[2] *= idet;
2635:         invJ[3]  = idet * (J[5] * J[6] - J[3] * J[8]);
2636:         invJ[4]  = idet * (J[0] * J[8] - J[2] * J[6]);
2637:         invJ[5]  = idet * (J[2] * J[3] - J[0] * J[5]);
2638:         invJ[6]  = idet * (J[3] * J[7] - J[4] * J[6]);
2639:         invJ[7]  = idet * (J[1] * J[6] - J[0] * J[7]);
2640:         invJ[8]  = idet * (J[0] * J[4] - J[1] * J[3]);
2641:       }
2642:       break;
2643:     }
2644:     for (l = 0; l < dimR; l++) {
2645:       for (m = 0; m < dimC; m++) {
2646:         guess[l] += PetscRealPart(invJ[l * dimC + m]) * resNeg[m];
2647:       }
2648:     }
2649:   } else {
2650: #if defined(PETSC_USE_COMPLEX)
2651:     char transpose = 'C';
2652: #else
2653:     char transpose = 'T';
2654: #endif
2655:     PetscBLASInt m = dimR;
2656:     PetscBLASInt n = dimC;
2657:     PetscBLASInt one = 1;
2658:     PetscBLASInt worksize = dimR * dimC, info;

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

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

2665:     for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2666:   }
2667:   return(0);
2668: }

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

2680:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2681:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2682:   DMGetWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2683:   DMGetWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2684:   cellCoords = &cellData[0];
2685:   cellCoeffs = &cellData[coordSize];
2686:   extJ       = &cellData[2 * coordSize];
2687:   resNeg     = &cellData[2 * coordSize + dimR];
2688:   invJ       = &J[dimR * dimC];
2689:   work       = &J[2 * dimR * dimC];
2690:   if (dimR == 2) {
2691:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2696:       for (j = 0; j < dimC; j++) {
2697:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2698:       }
2699:     }
2700:   } else if (dimR == 3) {
2701:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2706:       for (j = 0; j < dimC; j++) {
2707:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2708:       }
2709:     }
2710:   } else {
2711:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2712:   }
2713:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2714:   for (i = 0; i < dimR; i++) {
2715:     PetscReal *swap;

2717:     for (j = 0; j < (numV / 2); j++) {
2718:       for (k = 0; k < dimC; k++) {
2719:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2720:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2721:       }
2722:     }

2724:     if (i < dimR - 1) {
2725:       swap = cellCoeffs;
2726:       cellCoeffs = cellCoords;
2727:       cellCoords = swap;
2728:     }
2729:   }
2730:   PetscArrayzero(refCoords,numPoints * dimR);
2731:   for (j = 0; j < numPoints; j++) {
2732:     for (i = 0; i < maxIts; i++) {
2733:       PetscReal *guess = &refCoords[dimR * j];

2735:       /* compute -residual and Jacobian */
2736:       for (k = 0; k < dimC; k++) {resNeg[k] = realCoords[dimC * j + k];}
2737:       for (k = 0; k < dimC * dimR; k++) {J[k] = 0.;}
2738:       for (k = 0; k < numV; k++) {
2739:         PetscReal extCoord = 1.;
2740:         for (l = 0; l < dimR; l++) {
2741:           PetscReal coord = guess[l];
2742:           PetscInt  dep   = (k & (1 << l)) >> l;

2744:           extCoord *= dep * coord + !dep;
2745:           extJ[l] = dep;

2747:           for (m = 0; m < dimR; m++) {
2748:             PetscReal coord = guess[m];
2749:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
2750:             PetscReal mult  = dep * coord + !dep;

2752:             extJ[l] *= mult;
2753:           }
2754:         }
2755:         for (l = 0; l < dimC; l++) {
2756:           PetscReal coeff = cellCoeffs[dimC * k + l];

2758:           resNeg[l] -= coeff * extCoord;
2759:           for (m = 0; m < dimR; m++) {
2760:             J[dimR * l + m] += coeff * extJ[m];
2761:           }
2762:         }
2763:       }
2764: #if 0 && defined(PETSC_USE_DEBUG)
2765:       {
2766:         PetscReal maxAbs = 0.;

2768:         for (l = 0; l < dimC; l++) {
2769:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2770:         }
2771:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2772:       }
2773: #endif

2775:       DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2776:     }
2777:   }
2778:   DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2779:   DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2780:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2781:   return(0);
2782: }

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

2793:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2794:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2795:   DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2796:   cellCoords = &cellData[0];
2797:   cellCoeffs = &cellData[coordSize];
2798:   if (dimR == 2) {
2799:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2804:       for (j = 0; j < dimC; j++) {
2805:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2806:       }
2807:     }
2808:   } else if (dimR == 3) {
2809:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2814:       for (j = 0; j < dimC; j++) {
2815:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2816:       }
2817:     }
2818:   } else {
2819:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2820:   }
2821:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2822:   for (i = 0; i < dimR; i++) {
2823:     PetscReal *swap;

2825:     for (j = 0; j < (numV / 2); j++) {
2826:       for (k = 0; k < dimC; k++) {
2827:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2828:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2829:       }
2830:     }

2832:     if (i < dimR - 1) {
2833:       swap = cellCoeffs;
2834:       cellCoeffs = cellCoords;
2835:       cellCoords = swap;
2836:     }
2837:   }
2838:   PetscArrayzero(realCoords,numPoints * dimC);
2839:   for (j = 0; j < numPoints; j++) {
2840:     const PetscReal *guess  = &refCoords[dimR * j];
2841:     PetscReal       *mapped = &realCoords[dimC * j];

2843:     for (k = 0; k < numV; k++) {
2844:       PetscReal extCoord = 1.;
2845:       for (l = 0; l < dimR; l++) {
2846:         PetscReal coord = guess[l];
2847:         PetscInt  dep   = (k & (1 << l)) >> l;

2849:         extCoord *= dep * coord + !dep;
2850:       }
2851:       for (l = 0; l < dimC; l++) {
2852:         PetscReal coeff = cellCoeffs[dimC * k + l];

2854:         mapped[l] += coeff * extCoord;
2855:       }
2856:     }
2857:   }
2858:   DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2859:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2860:   return(0);
2861: }

2863: /* TODO: TOBY please fix this for Nc > 1 */
2864: static PetscErrorCode DMPlexCoordinatesToReference_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2865: {
2866:   PetscInt       numComp, pdim, i, j, k, l, m, maxIter = 7, coordSize;
2867:   PetscScalar    *nodes = NULL;
2868:   PetscReal      *invV, *modes;
2869:   PetscReal      *B, *D, *resNeg;
2870:   PetscScalar    *J, *invJ, *work;

2874:   PetscFEGetDimension(fe, &pdim);
2875:   PetscFEGetNumComponents(fe, &numComp);
2876:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2877:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2878:   /* convert nodes to values in the stable evaluation basis */
2879:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2880:   invV = fe->invV;
2881:   for (i = 0; i < pdim; ++i) {
2882:     modes[i] = 0.;
2883:     for (j = 0; j < pdim; ++j) {
2884:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2885:     }
2886:   }
2887:   DMGetWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2888:   D      = &B[pdim*Nc];
2889:   resNeg = &D[pdim*Nc * dimR];
2890:   DMGetWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2891:   invJ = &J[Nc * dimR];
2892:   work = &invJ[Nc * dimR];
2893:   for (i = 0; i < numPoints * dimR; i++) {refCoords[i] = 0.;}
2894:   for (j = 0; j < numPoints; j++) {
2895:       for (i = 0; i < maxIter; i++) { /* we could batch this so that we're not making big B and D arrays all the time */
2896:       PetscReal *guess = &refCoords[j * dimR];
2897:       PetscSpaceEvaluate(fe->basisSpace, 1, guess, B, D, NULL);
2898:       for (k = 0; k < Nc; k++) {resNeg[k] = realCoords[j * Nc + k];}
2899:       for (k = 0; k < Nc * dimR; k++) {J[k] = 0.;}
2900:       for (k = 0; k < pdim; k++) {
2901:         for (l = 0; l < Nc; l++) {
2902:           resNeg[l] -= modes[k] * B[k * Nc + l];
2903:           for (m = 0; m < dimR; m++) {
2904:             J[l * dimR + m] += modes[k] * D[(k * Nc + l) * dimR + m];
2905:           }
2906:         }
2907:       }
2908: #if 0 && defined(PETSC_USE_DEBUG)
2909:       {
2910:         PetscReal maxAbs = 0.;

2912:         for (l = 0; l < Nc; l++) {
2913:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2914:         }
2915:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,(double) maxAbs);
2916:       }
2917: #endif
2918:       DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2919:     }
2920:   }
2921:   DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2922:   DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2923:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2924:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2925:   return(0);
2926: }

2928: /* TODO: TOBY please fix this for Nc > 1 */
2929: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2930: {
2931:   PetscInt       numComp, pdim, i, j, k, l, coordSize;
2932:   PetscScalar    *nodes = NULL;
2933:   PetscReal      *invV, *modes;
2934:   PetscReal      *B;

2938:   PetscFEGetDimension(fe, &pdim);
2939:   PetscFEGetNumComponents(fe, &numComp);
2940:   if (numComp != Nc) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"coordinate discretization must have as many components (%D) as embedding dimension (!= %D)",numComp,Nc);
2941:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2942:   /* convert nodes to values in the stable evaluation basis */
2943:   DMGetWorkArray(dm,pdim,MPIU_REAL,&modes);
2944:   invV = fe->invV;
2945:   for (i = 0; i < pdim; ++i) {
2946:     modes[i] = 0.;
2947:     for (j = 0; j < pdim; ++j) {
2948:       modes[i] += invV[i * pdim + j] * PetscRealPart(nodes[j]);
2949:     }
2950:   }
2951:   DMGetWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2952:   PetscSpaceEvaluate(fe->basisSpace, numPoints, refCoords, B, NULL, NULL);
2953:   for (i = 0; i < numPoints * Nc; i++) {realCoords[i] = 0.;}
2954:   for (j = 0; j < numPoints; j++) {
2955:     PetscReal *mapped = &realCoords[j * Nc];

2957:     for (k = 0; k < pdim; k++) {
2958:       for (l = 0; l < Nc; l++) {
2959:         mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2960:       }
2961:     }
2962:   }
2963:   DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2964:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2965:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2966:   return(0);
2967: }

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

2974:   Not collective

2976:   Input Parameters:
2977: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2978:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2979:                as a multilinear map for tensor-product elements
2980: . cell       - the cell whose map is used.
2981: . numPoints  - the number of points to locate
2982: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

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

2987:   Level: intermediate

2989: .seealso: DMPlexReferenceToCoordinates()
2990: @*/
2991: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2992: {
2993:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
2994:   DM             coordDM = NULL;
2995:   Vec            coords;
2996:   PetscFE        fe = NULL;

3001:   DMGetDimension(dm,&dimR);
3002:   DMGetCoordinateDim(dm,&dimC);
3003:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3004:   DMPlexGetDepth(dm,&depth);
3005:   DMGetCoordinatesLocal(dm,&coords);
3006:   DMGetCoordinateDM(dm,&coordDM);
3007:   if (coordDM) {
3008:     PetscInt coordFields;

3010:     DMGetNumFields(coordDM,&coordFields);
3011:     if (coordFields) {
3012:       PetscClassId id;
3013:       PetscObject  disc;

3015:       DMGetField(coordDM,0,NULL,&disc);
3016:       PetscObjectGetClassId(disc,&id);
3017:       if (id == PETSCFE_CLASSID) {
3018:         fe = (PetscFE) disc;
3019:       }
3020:     }
3021:   }
3022:   DMPlexGetInteriorCellStratum(dm,&cStart,&cEnd);
3023:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3024:   if (!fe) { /* implicit discretization: affine or multilinear */
3025:     PetscInt  coneSize;
3026:     PetscBool isSimplex, isTensor;

3028:     DMPlexGetConeSize(dm,cell,&coneSize);
3029:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3030:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3031:     if (isSimplex) {
3032:       PetscReal detJ, *v0, *J, *invJ;

3034:       DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3035:       J    = &v0[dimC];
3036:       invJ = &J[dimC * dimC];
3037:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
3038:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
3039:         const PetscReal x0[3] = {-1.,-1.,-1.};

3041:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
3042:       }
3043:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3044:     } else if (isTensor) {
3045:       DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3046:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3047:   } else {
3048:     DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
3049:   }
3050:   return(0);
3051: }

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

3056:   Not collective

3058:   Input Parameters:
3059: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3060:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3061:                as a multilinear map for tensor-product elements
3062: . cell       - the cell whose map is used.
3063: . numPoints  - the number of points to locate
3064: - refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

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

3069:    Level: intermediate

3071: .seealso: DMPlexCoordinatesToReference()
3072: @*/
3073: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3074: {
3075:   PetscInt       dimC, dimR, depth, cStart, cEnd, i;
3076:   DM             coordDM = NULL;
3077:   Vec            coords;
3078:   PetscFE        fe = NULL;

3083:   DMGetDimension(dm,&dimR);
3084:   DMGetCoordinateDim(dm,&dimC);
3085:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3086:   DMPlexGetDepth(dm,&depth);
3087:   DMGetCoordinatesLocal(dm,&coords);
3088:   DMGetCoordinateDM(dm,&coordDM);
3089:   if (coordDM) {
3090:     PetscInt coordFields;

3092:     DMGetNumFields(coordDM,&coordFields);
3093:     if (coordFields) {
3094:       PetscClassId id;
3095:       PetscObject  disc;

3097:       DMGetField(coordDM,0,NULL,&disc);
3098:       PetscObjectGetClassId(disc,&id);
3099:       if (id == PETSCFE_CLASSID) {
3100:         fe = (PetscFE) disc;
3101:       }
3102:     }
3103:   }
3104:   DMPlexGetInteriorCellStratum(dm,&cStart,&cEnd);
3105:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3106:   if (!fe) { /* implicit discretization: affine or multilinear */
3107:     PetscInt  coneSize;
3108:     PetscBool isSimplex, isTensor;

3110:     DMPlexGetConeSize(dm,cell,&coneSize);
3111:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3112:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3113:     if (isSimplex) {
3114:       PetscReal detJ, *v0, *J;

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

3122:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3123:       }
3124:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3125:     } else if (isTensor) {
3126:       DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3127:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3128:   } else {
3129:     DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3130:   }
3131:   return(0);
3132: }

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

3137:   Not collective

3139:   Input Parameters:
3140: + dm      - The DM
3141: . time    - The time
3142: - func    - The function transforming current coordinates to new coordaintes

3144:    Calling sequence of func:
3145: $    func(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3146: $         const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3147: $         const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3148: $         PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f[]);

3150: +  dim          - The spatial dimension
3151: .  Nf           - The number of input fields (here 1)
3152: .  NfAux        - The number of input auxiliary fields
3153: .  uOff         - The offset of the coordinates in u[] (here 0)
3154: .  uOff_x       - The offset of the coordinates in u_x[] (here 0)
3155: .  u            - The coordinate values at this point in space
3156: .  u_t          - The coordinate time derivative at this point in space (here NULL)
3157: .  u_x          - The coordinate derivatives at this point in space
3158: .  aOff         - The offset of each auxiliary field in u[]
3159: .  aOff_x       - The offset of each auxiliary field in u_x[]
3160: .  a            - The auxiliary field values at this point in space
3161: .  a_t          - The auxiliary field time derivative at this point in space (or NULL)
3162: .  a_x          - The auxiliary field derivatives at this point in space
3163: .  t            - The current time
3164: .  x            - The coordinates of this point (here not used)
3165: .  numConstants - The number of constants
3166: .  constants    - The value of each constant
3167: -  f            - The new coordinates at this point in space

3169:   Level: intermediate

3171: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMProjectFieldLocal(), DMProjectFieldLabelLocal()
3172: @*/
3173: PetscErrorCode DMPlexRemapGeometry(DM dm, PetscReal time,
3174:                                    void (*func)(PetscInt, PetscInt, PetscInt,
3175:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3176:                                                 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
3177:                                                 PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]))
3178: {
3179:   DM             cdm;
3180:   Vec            lCoords, tmpCoords;

3184:   DMGetCoordinateDM(dm, &cdm);
3185:   DMGetCoordinatesLocal(dm, &lCoords);
3186:   DMGetLocalVector(cdm, &tmpCoords);
3187:   VecCopy(lCoords, tmpCoords);
3188:   DMProjectFieldLocal(cdm, time, tmpCoords, &func, INSERT_VALUES, lCoords);
3189:   DMRestoreLocalVector(cdm, &tmpCoords);
3190:   return(0);
3191: }

3193: /* Shear applies the transformation, assuming we fix z,
3194:   / 1  0  m_0 \
3195:   | 0  1  m_1 |
3196:   \ 0  0   1  /
3197: */
3198: static void f0_shear(PetscInt dim, PetscInt Nf, PetscInt NfAux,
3199:                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
3200:                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
3201:                      PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar coords[])
3202: {
3203:   const PetscInt Nc = uOff[1]-uOff[0];
3204:   const PetscInt ax = (PetscInt) PetscRealPart(constants[0]);
3205:   PetscInt       c;

3207:   for (c = 0; c < Nc; ++c) {
3208:     coords[c] = u[c] + constants[c+1]*u[ax];
3209:   }
3210: }

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

3215:   Not collective

3217:   Input Parameters:
3218: + dm          - The DM
3219: . direction   - The shear coordinate direction, e.g. 0 is the x-axis
3220: - multipliers - The multiplier m for each direction which is not the shear direction

3222:   Level: intermediate

3224: .seealso: DMPlexRemapGeometry()
3225: @*/
3226: PetscErrorCode DMPlexShearGeometry(DM dm, DMDirection direction, PetscReal multipliers[])
3227: {
3228:   DM             cdm;
3229:   PetscDS        cds;
3230:   PetscObject    obj;
3231:   PetscClassId   id;
3232:   PetscScalar   *moduli;
3233:   const PetscInt dir = (PetscInt) direction;
3234:   PetscInt       dE, d, e;

3238:   DMGetCoordinateDM(dm, &cdm);
3239:   DMGetCoordinateDim(dm, &dE);
3240:   PetscMalloc1(dE+1, &moduli);
3241:   moduli[0] = dir;
3242:   for (d = 0, e = 0; d < dE; ++d) moduli[d] = d == dir ? 0.0 : (multipliers ? multipliers[e++] : 1.0);
3243:   DMGetDS(cdm, &cds);
3244:   PetscDSGetDiscretization(cds, 0, &obj);
3245:   PetscObjectGetClassId(obj, &id);
3246:   if (id != PETSCFE_CLASSID) {
3247:     Vec           lCoords;
3248:     PetscSection  cSection;
3249:     PetscScalar  *coords;
3250:     PetscInt      vStart, vEnd, v;

3252:     DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
3253:     DMGetCoordinateSection(dm, &cSection);
3254:     DMGetCoordinatesLocal(dm, &lCoords);
3255:     VecGetArray(lCoords, &coords);
3256:     for (v = vStart; v < vEnd; ++v) {
3257:       PetscReal ds;
3258:       PetscInt  off, c;

3260:       PetscSectionGetOffset(cSection, v, &off);
3261:       ds   = PetscRealPart(coords[off+dir]);
3262:       for (c = 0; c < dE; ++c) coords[off+c] += moduli[c]*ds;
3263:     }
3264:     VecRestoreArray(lCoords, &coords);
3265:   } else {
3266:     PetscDSSetConstants(cds, dE+1, moduli);
3267:     DMPlexRemapGeometry(dm, 0.0, f0_shear);
3268:   }
3269:   PetscFree(moduli);
3270:   return(0);
3271: }