Actual source code: plexgeometry.c

petsc-master 2019-06-15
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: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
  7: {
  8:   const PetscReal p0_x  = segmentA[0*2+0];
  9:   const PetscReal p0_y  = segmentA[0*2+1];
 10:   const PetscReal p1_x  = segmentA[1*2+0];
 11:   const PetscReal p1_y  = segmentA[1*2+1];
 12:   const PetscReal p2_x  = segmentB[0*2+0];
 13:   const PetscReal p2_y  = segmentB[0*2+1];
 14:   const PetscReal p3_x  = segmentB[1*2+0];
 15:   const PetscReal p3_y  = segmentB[1*2+1];
 16:   const PetscReal s1_x  = p1_x - p0_x;
 17:   const PetscReal s1_y  = p1_y - p0_y;
 18:   const PetscReal s2_x  = p3_x - p2_x;
 19:   const PetscReal s2_y  = p3_y - p2_y;
 20:   const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

 23:   *hasIntersection = PETSC_FALSE;
 24:   /* Non-parallel lines */
 25:   if (denom != 0.0) {
 26:     const PetscReal s = (-s1_y * (p0_x - p2_x) + s1_x * (p0_y - p2_y)) / denom;
 27:     const PetscReal t = ( s2_x * (p0_y - p2_y) - s2_y * (p0_x - p2_x)) / denom;

 29:     if (s >= 0 && s <= 1 && t >= 0 && t <= 1) {
 30:       *hasIntersection = PETSC_TRUE;
 31:       if (intersection) {
 32:         intersection[0] = p0_x + (t * s1_x);
 33:         intersection[1] = p0_y + (t * s1_y);
 34:       }
 35:     }
 36:   }
 37:   return(0);
 38: }

 40: static PetscErrorCode DMPlexLocatePoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 41: {
 42:   const PetscInt  embedDim = 2;
 43:   const PetscReal eps      = PETSC_SQRT_MACHINE_EPSILON;
 44:   PetscReal       x        = PetscRealPart(point[0]);
 45:   PetscReal       y        = PetscRealPart(point[1]);
 46:   PetscReal       v0[2], J[4], invJ[4], detJ;
 47:   PetscReal       xi, eta;
 48:   PetscErrorCode  ierr;

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

 55:   if ((xi >= -eps) && (eta >= -eps) && (xi + eta <= 2.0+eps)) *cell = c;
 56:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
 57:   return(0);
 58: }

 60: static PetscErrorCode DMPlexClosestPoint_Simplex_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscReal cpoint[])
 61: {
 62:   const PetscInt  embedDim = 2;
 63:   PetscReal       x        = PetscRealPart(point[0]);
 64:   PetscReal       y        = PetscRealPart(point[1]);
 65:   PetscReal       v0[2], J[4], invJ[4], detJ;
 66:   PetscReal       xi, eta, r;
 67:   PetscErrorCode  ierr;

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

 74:   xi  = PetscMax(xi,  0.0);
 75:   eta = PetscMax(eta, 0.0);
 76:   if (xi + eta > 2.0) {
 77:     r    = (xi + eta)/2.0;
 78:     xi  /= r;
 79:     eta /= r;
 80:   }
 81:   cpoint[0] = J[0*embedDim+0]*xi + J[0*embedDim+1]*eta + v0[0];
 82:   cpoint[1] = J[1*embedDim+0]*xi + J[1*embedDim+1]*eta + v0[1];
 83:   return(0);
 84: }

 86: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 87: {
 88:   PetscSection       coordSection;
 89:   Vec             coordsLocal;
 90:   PetscScalar    *coords = NULL;
 91:   const PetscInt  faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
 92:   PetscReal       x         = PetscRealPart(point[0]);
 93:   PetscReal       y         = PetscRealPart(point[1]);
 94:   PetscInt        crossings = 0, f;
 95:   PetscErrorCode  ierr;

 98:   DMGetCoordinatesLocal(dm, &coordsLocal);
 99:   DMGetCoordinateSection(dm, &coordSection);
100:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
101:   for (f = 0; f < 4; ++f) {
102:     PetscReal x_i   = PetscRealPart(coords[faces[2*f+0]*2+0]);
103:     PetscReal y_i   = PetscRealPart(coords[faces[2*f+0]*2+1]);
104:     PetscReal x_j   = PetscRealPart(coords[faces[2*f+1]*2+0]);
105:     PetscReal y_j   = PetscRealPart(coords[faces[2*f+1]*2+1]);
106:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
107:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
108:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
109:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
110:     if ((cond1 || cond2)  && above) ++crossings;
111:   }
112:   if (crossings % 2) *cell = c;
113:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
114:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
115:   return(0);
116: }

118: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
119: {
120:   const PetscInt embedDim = 3;
121:   PetscReal      v0[3], J[9], invJ[9], detJ;
122:   PetscReal      x = PetscRealPart(point[0]);
123:   PetscReal      y = PetscRealPart(point[1]);
124:   PetscReal      z = PetscRealPart(point[2]);
125:   PetscReal      xi, eta, zeta;

129:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
130:   xi   = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
131:   eta  = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
132:   zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);

134:   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
135:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
136:   return(0);
137: }

139: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
140: {
141:   PetscSection   coordSection;
142:   Vec            coordsLocal;
143:   PetscScalar   *coords = NULL;
144:   const PetscInt faces[24] = {0, 3, 2, 1,  5, 4, 7, 6,  3, 0, 4, 5,
145:                               1, 2, 6, 7,  3, 5, 6, 2,  0, 1, 7, 4};
146:   PetscBool      found = PETSC_TRUE;
147:   PetscInt       f;

151:   DMGetCoordinatesLocal(dm, &coordsLocal);
152:   DMGetCoordinateSection(dm, &coordSection);
153:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
154:   for (f = 0; f < 6; ++f) {
155:     /* Check the point is under plane */
156:     /*   Get face normal */
157:     PetscReal v_i[3];
158:     PetscReal v_j[3];
159:     PetscReal normal[3];
160:     PetscReal pp[3];
161:     PetscReal dot;

163:     v_i[0]    = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
164:     v_i[1]    = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
165:     v_i[2]    = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
166:     v_j[0]    = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
167:     v_j[1]    = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
168:     v_j[2]    = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
169:     normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
170:     normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
171:     normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
172:     pp[0]     = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
173:     pp[1]     = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
174:     pp[2]     = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
175:     dot       = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];

177:     /* Check that projected point is in face (2D location problem) */
178:     if (dot < 0.0) {
179:       found = PETSC_FALSE;
180:       break;
181:     }
182:   }
183:   if (found) *cell = c;
184:   else *cell = DMLOCATEPOINT_POINT_NOT_FOUND;
185:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
186:   return(0);
187: }

189: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
190: {
191:   PetscInt d;

194:   box->dim = dim;
195:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
196:   return(0);
197: }

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

204:   PetscMalloc1(1, box);
205:   PetscGridHashInitialize_Internal(*box, dim, point);
206:   return(0);
207: }

209: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
210: {
211:   PetscInt d;

214:   for (d = 0; d < box->dim; ++d) {
215:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
216:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
217:   }
218:   return(0);
219: }

221: /*
222:   PetscGridHashSetGrid - Divide the grid into boxes

224:   Not collective

226:   Input Parameters:
227: + box - The grid hash object
228: . n   - The number of boxes in each dimension, or PETSC_DETERMINE
229: - h   - The box size in each dimension, only used if n[d] == PETSC_DETERMINE

231:   Level: developer

233: .seealso: PetscGridHashCreate()
234: */
235: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
236: {
237:   PetscInt d;

240:   for (d = 0; d < box->dim; ++d) {
241:     box->extent[d] = box->upper[d] - box->lower[d];
242:     if (n[d] == PETSC_DETERMINE) {
243:       box->h[d] = h[d];
244:       box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
245:     } else {
246:       box->n[d] = n[d];
247:       box->h[d] = box->extent[d]/n[d];
248:     }
249:   }
250:   return(0);
251: }

253: /*
254:   PetscGridHashGetEnclosingBox - Find the grid boxes containing each input point

256:   Not collective

258:   Input Parameters:
259: + box       - The grid hash object
260: . numPoints - The number of input points
261: - points    - The input point coordinates

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

267:   Level: developer

269: .seealso: PetscGridHashCreate()
270: */
271: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
272: {
273:   const PetscReal *lower = box->lower;
274:   const PetscReal *upper = box->upper;
275:   const PetscReal *h     = box->h;
276:   const PetscInt  *n     = box->n;
277:   const PetscInt   dim   = box->dim;
278:   PetscInt         d, p;

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

285:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
286:       if (dbox == -1   && PetscAbsReal(PetscRealPart(points[p*dim+d]) - lower[d]) < 1.0e-9) dbox = 0;
287:       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",
288:                                              p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0);
289:       dboxes[p*dim+d] = dbox;
290:     }
291:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
292:   }
293:   return(0);
294: }

296: /*
297:  PetscGridHashGetEnclosingBoxQuery - Find the grid boxes containing each input point

299:  Not collective

301:   Input Parameters:
302: + box       - The grid hash object
303: . numPoints - The number of input points
304: - points    - The input point coordinates

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

311:   Level: developer

313: .seealso: PetscGridHashGetEnclosingBox()
314: */
315: PetscErrorCode PetscGridHashGetEnclosingBoxQuery(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[],PetscBool *found)
316: {
317:   const PetscReal *lower = box->lower;
318:   const PetscReal *upper = box->upper;
319:   const PetscReal *h     = box->h;
320:   const PetscInt  *n     = box->n;
321:   const PetscInt   dim   = box->dim;
322:   PetscInt         d, p;

325:   *found = PETSC_FALSE;
326:   for (p = 0; p < numPoints; ++p) {
327:     for (d = 0; d < dim; ++d) {
328:       PetscInt dbox = PetscFloorReal((PetscRealPart(points[p*dim+d]) - lower[d])/h[d]);

330:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
331:       if (dbox < 0 || dbox >= n[d]) {
332:         return(0);
333:       }
334:       dboxes[p*dim+d] = dbox;
335:     }
336:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
337:   }
338:   *found = PETSC_TRUE;
339:   return(0);
340: }

342: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
343: {

347:   if (*box) {
348:     PetscSectionDestroy(&(*box)->cellSection);
349:     ISDestroy(&(*box)->cells);
350:     DMLabelDestroy(&(*box)->cellsSparse);
351:   }
352:   PetscFree(*box);
353:   return(0);
354: }

356: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
357: {
358:   PetscInt       coneSize;

362:   switch (dim) {
363:   case 2:
364:     DMPlexGetConeSize(dm, cellStart, &coneSize);
365:     switch (coneSize) {
366:     case 3:
367:       DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
368:       break;
369:     case 4:
370:       DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
371:       break;
372:     default:
373:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
374:     }
375:     break;
376:   case 3:
377:     DMPlexGetConeSize(dm, cellStart, &coneSize);
378:     switch (coneSize) {
379:     case 4:
380:       DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
381:       break;
382:     case 6:
383:       DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
384:       break;
385:     default:
386:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
387:     }
388:     break;
389:   default:
390:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
391:   }
392:   return(0);
393: }

395: /*
396:   DMPlexClosestPoint_Internal - Returns the closest point in the cell to the given point
397: */
398: PetscErrorCode DMPlexClosestPoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cell, PetscReal cpoint[])
399: {
400:   PetscInt       coneSize;

404:   switch (dim) {
405:   case 2:
406:     DMPlexGetConeSize(dm, cell, &coneSize);
407:     switch (coneSize) {
408:     case 3:
409:       DMPlexClosestPoint_Simplex_2D_Internal(dm, point, cell, cpoint);
410:       break;
411: #if 0
412:     case 4:
413:       DMPlexClosestPoint_General_2D_Internal(dm, point, cell, cpoint);
414:       break;
415: #endif
416:     default:
417:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
418:     }
419:     break;
420: #if 0
421:   case 3:
422:     DMPlexGetConeSize(dm, cell, &coneSize);
423:     switch (coneSize) {
424:     case 4:
425:       DMPlexClosestPoint_Simplex_3D_Internal(dm, point, cell, cpoint);
426:       break;
427:     case 6:
428:       DMPlexClosestPoint_General_3D_Internal(dm, point, cell, cpoint);
429:       break;
430:     default:
431:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for cell with cone size %D", coneSize);
432:     }
433:     break;
434: #endif
435:   default:
436:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No closest point location for mesh dimension %D", dim);
437:   }
438:   return(0);
439: }

441: /*
442:   DMPlexComputeGridHash_Internal - Create a grid hash structure covering the Plex

444:   Collective on dm

446:   Input Parameter:
447: . dm - The Plex

449:   Output Parameter:
450: . localBox - The grid hash object

452:   Level: developer

454: .seealso: PetscGridHashCreate(), PetscGridHashGetEnclosingBox()
455: */
456: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
457: {
458:   MPI_Comm           comm;
459:   PetscGridHash      lbox;
460:   Vec                coordinates;
461:   PetscSection       coordSection;
462:   Vec                coordsLocal;
463:   const PetscScalar *coords;
464:   PetscInt          *dboxes, *boxes;
465:   PetscInt           n[3] = {10, 10, 10};
466:   PetscInt           dim, N, cStart, cEnd, cMax, c, i;
467:   PetscErrorCode     ierr;

470:   PetscObjectGetComm((PetscObject) dm, &comm);
471:   DMGetCoordinatesLocal(dm, &coordinates);
472:   DMGetCoordinateDim(dm, &dim);
473:   if (dim != 2) SETERRQ(comm, PETSC_ERR_SUP, "I have only coded this for 2D");
474:   VecGetLocalSize(coordinates, &N);
475:   VecGetArrayRead(coordinates, &coords);
476:   PetscGridHashCreate(comm, dim, coords, &lbox);
477:   for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
478:   VecRestoreArrayRead(coordinates, &coords);
479:   PetscOptionsGetInt(NULL,NULL,"-dm_plex_hash_box_nijk",&n[0],NULL);
480:   n[1] = n[0];
481:   n[2] = n[0];
482:   PetscGridHashSetGrid(lbox, n, NULL);
483: #if 0
484:   /* Could define a custom reduction to merge these */
485:   MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
486:   MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
487: #endif
488:   /* Is there a reason to snap the local bounding box to a division of the global box? */
489:   /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
490:   /* Create label */
491:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
492:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
493:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
494:   DMLabelCreate(PETSC_COMM_SELF, "cells", &lbox->cellsSparse);
495:   DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
496:   /* Compute boxes which overlap each cell: https://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
497:   DMGetCoordinatesLocal(dm, &coordsLocal);
498:   DMGetCoordinateSection(dm, &coordSection);
499:   PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
500:   for (c = cStart; c < cEnd; ++c) {
501:     const PetscReal *h       = lbox->h;
502:     PetscScalar     *ccoords = NULL;
503:     PetscInt         csize   = 0;
504:     PetscScalar      point[3];
505:     PetscInt         dlim[6], d, e, i, j, k;

507:     /* Find boxes enclosing each vertex */
508:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
509:     PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
510:     /* Mark cells containing the vertices */
511:     for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
512:     /* Get grid of boxes containing these */
513:     for (d = 0;   d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
514:     for (d = dim; d < 3;   ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
515:     for (e = 1; e < dim+1; ++e) {
516:       for (d = 0; d < dim; ++d) {
517:         dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
518:         dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
519:       }
520:     }
521:     /* Check for intersection of box with cell */
522:     for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
523:       for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
524:         for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
525:           const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
526:           PetscScalar    cpoint[3];
527:           PetscInt       cell, edge, ii, jj, kk;

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

534:                 DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
535:                 if (cell >= 0) { DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
536:               }
537:             }
538:           }
539:           /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
540:           for (edge = 0; edge < dim+1; ++edge) {
541:             PetscReal segA[6], segB[6];

543:             if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Unexpected dim %d > 3",dim);
544:             for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
545:             for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
546:               if (dim > 2) {segB[2]     = PetscRealPart(point[2]);
547:                             segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
548:               for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
549:                 if (dim > 1) {segB[1]     = PetscRealPart(point[1]);
550:                               segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
551:                 for (ii = 0; ii < 2; ++ii) {
552:                   PetscBool intersects;

554:                   segB[0]     = PetscRealPart(point[0]);
555:                   segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
556:                   DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
557:                   if (intersects) { DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
558:                 }
559:               }
560:             }
561:           }
562:         }
563:       }
564:     }
565:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
566:   }
567:   PetscFree2(dboxes, boxes);
568:   DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
569:   DMLabelDestroy(&lbox->cellsSparse);
570:   *localBox = lbox;
571:   return(0);
572: }

574: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, DMPointLocationType ltype, PetscSF cellSF)
575: {
576:   DM_Plex        *mesh = (DM_Plex *) dm->data;
577:   PetscBool       hash = mesh->useHashLocation, reuse = PETSC_FALSE;
578:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
579:   PetscInt        dim, cStart, cEnd, cMax, numCells, c, d;
580:   const PetscInt *boxCells;
581:   PetscSFNode    *cells;
582:   PetscScalar    *a;
583:   PetscMPIInt     result;
584:   PetscLogDouble  t0,t1;
585:   PetscReal       gmin[3],gmax[3];
586:   PetscInt        terminating_query_type[] = { 0, 0, 0 };
587:   PetscErrorCode  ierr;

590:   PetscTime(&t0);
591:   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.");
592:   DMGetCoordinateDim(dm, &dim);
593:   VecGetBlockSize(v, &bs);
594:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
595:   if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
596:   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);
597:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
598:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
599:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
600:   VecGetLocalSize(v, &numPoints);
601:   VecGetArray(v, &a);
602:   numPoints /= bs;
603:   {
604:     const PetscSFNode *sf_cells;

606:     PetscSFGetGraph(cellSF,NULL,NULL,NULL,&sf_cells);
607:     if (sf_cells) {
608:       PetscInfo(dm,"[DMLocatePoints_Plex] Re-using existing StarForest node list\n");
609:       cells = (PetscSFNode*)sf_cells;
610:       reuse = PETSC_TRUE;
611:     } else {
612:       PetscInfo(dm,"[DMLocatePoints_Plex] Creating and initializing new StarForest node list\n");
613:       PetscMalloc1(numPoints, &cells);
614:       /* initialize cells if created */
615:       for (p=0; p<numPoints; p++) {
616:         cells[p].rank  = 0;
617:         cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
618:       }
619:     }
620:   }
621:   /* define domain bounding box */
622:   {
623:     Vec coorglobal;

625:     DMGetCoordinates(dm,&coorglobal);
626:     VecStrideMaxAll(coorglobal,NULL,gmax);
627:     VecStrideMinAll(coorglobal,NULL,gmin);
628:   }
629:   if (hash) {
630:     if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
631:     /* Designate the local box for each point */
632:     /* Send points to correct process */
633:     /* Search cells that lie in each subbox */
634:     /*   Should we bin points before doing search? */
635:     ISGetIndices(mesh->lbox->cells, &boxCells);
636:   }
637:   for (p = 0, numFound = 0; p < numPoints; ++p) {
638:     const PetscScalar *point = &a[p*bs];
639:     PetscInt           dbin[3] = {-1,-1,-1}, bin, cell = -1, cellOffset;
640:     PetscBool          point_outside_domain = PETSC_FALSE;

642:     /* check bounding box of domain */
643:     for (d=0; d<dim; d++) {
644:       if (PetscRealPart(point[d]) < gmin[d]) { point_outside_domain = PETSC_TRUE; break; }
645:       if (PetscRealPart(point[d]) > gmax[d]) { point_outside_domain = PETSC_TRUE; break; }
646:     }
647:     if (point_outside_domain) {
648:       cells[p].rank = 0;
649:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
650:       terminating_query_type[0]++;
651:       continue;
652:     }

654:     /* check initial values in cells[].index - abort early if found */
655:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
656:       c = cells[p].index;
657:       cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND;
658:       DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
659:       if (cell >= 0) {
660:         cells[p].rank = 0;
661:         cells[p].index = cell;
662:         numFound++;
663:       }
664:     }
665:     if (cells[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
666:       terminating_query_type[1]++;
667:       continue;
668:     }

670:     if (hash) {
671:       PetscBool found_box;

673:       /* allow for case that point is outside box - abort early */
674:       PetscGridHashGetEnclosingBoxQuery(mesh->lbox, 1, point, dbin, &bin,&found_box);
675:       if (found_box) {
676:         /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
677:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
678:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
679:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
680:           DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
681:           if (cell >= 0) {
682:             cells[p].rank = 0;
683:             cells[p].index = cell;
684:             numFound++;
685:             terminating_query_type[2]++;
686:             break;
687:           }
688:         }
689:       }
690:     } else {
691:       for (c = cStart; c < cEnd; ++c) {
692:         DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
693:         if (cell >= 0) {
694:           cells[p].rank = 0;
695:           cells[p].index = cell;
696:           numFound++;
697:           terminating_query_type[2]++;
698:           break;
699:         }
700:       }
701:     }
702:   }
703:   if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
704:   if (ltype == DM_POINTLOCATION_NEAREST && hash && numFound < numPoints) {
705:     for (p = 0; p < numPoints; p++) {
706:       const PetscScalar *point = &a[p*bs];
707:       PetscReal          cpoint[3], diff[3], dist, distMax = PETSC_MAX_REAL;
708:       PetscInt           dbin[3] = {-1,-1,-1}, bin, cellOffset, d;

710:       if (cells[p].index < 0) {
711:         ++numFound;
712:         PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
713:         PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
714:         PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
715:         for (c = cellOffset; c < cellOffset + numCells; ++c) {
716:           DMPlexClosestPoint_Internal(dm, dim, point, boxCells[c], cpoint);
717:           for (d = 0; d < dim; ++d) diff[d] = cpoint[d] - PetscRealPart(point[d]);
718:           dist = DMPlex_NormD_Internal(dim, diff);
719:           if (dist < distMax) {
720:             for (d = 0; d < dim; ++d) a[p*bs+d] = cpoint[d];
721:             cells[p].rank  = 0;
722:             cells[p].index = boxCells[c];
723:             distMax = dist;
724:             break;
725:           }
726:         }
727:       }
728:     }
729:   }
730:   /* This code is only be relevant when interfaced to parallel point location */
731:   /* Check for highest numbered proc that claims a point (do we care?) */
732:   if (ltype == DM_POINTLOCATION_REMOVE && numFound < numPoints) {
733:     PetscMalloc1(numFound,&found);
734:     for (p = 0, numFound = 0; p < numPoints; p++) {
735:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
736:         if (numFound < p) {
737:           cells[numFound] = cells[p];
738:         }
739:         found[numFound++] = p;
740:       }
741:     }
742:   }
743:   VecRestoreArray(v, &a);
744:   if (!reuse) {
745:     PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
746:   }
747:   PetscTime(&t1);
748:   if (hash) {
749:     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]);
750:   } else {
751:     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]);
752:   }
753:   PetscInfo3(dm,"[DMLocatePoints_Plex] npoints %D : time(rank0) %1.2e (sec): points/sec %1.4e\n",numPoints,t1-t0,(double)((double)numPoints/(t1-t0)));
754:   return(0);
755: }

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

760:   Not collective

762:   Input Parameter:
763: . coords - The coordinates of a segment

765:   Output Parameters:
766: + coords - The new y-coordinate, and 0 for x
767: - R - The rotation which accomplishes the projection

769:   Level: developer

771: .seealso: DMPlexComputeProjection3Dto1D(), DMPlexComputeProjection3Dto2D()
772: @*/
773: PetscErrorCode DMPlexComputeProjection2Dto1D(PetscScalar coords[], PetscReal R[])
774: {
775:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
776:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
777:   const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;

780:   R[0] = c; R[1] = -s;
781:   R[2] = s; R[3] =  c;
782:   coords[0] = 0.0;
783:   coords[1] = r;
784:   return(0);
785: }

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

790:   Not collective

792:   Input Parameter:
793: . coords - The coordinates of a segment

795:   Output Parameters:
796: + coords - The new y-coordinate, and 0 for x and z
797: - R - The rotation which accomplishes the projection

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

801:   Level: developer

803: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto2D()
804: @*/
805: PetscErrorCode DMPlexComputeProjection3Dto1D(PetscScalar coords[], PetscReal R[])
806: {
807:   PetscReal      x    = PetscRealPart(coords[3] - coords[0]);
808:   PetscReal      y    = PetscRealPart(coords[4] - coords[1]);
809:   PetscReal      z    = PetscRealPart(coords[5] - coords[2]);
810:   PetscReal      r    = PetscSqrtReal(x*x + y*y + z*z);
811:   PetscReal      rinv = 1. / r;

814:   x *= rinv; y *= rinv; z *= rinv;
815:   if (x > 0.) {
816:     PetscReal inv1pX   = 1./ (1. + x);

818:     R[0] = x; R[1] = -y;              R[2] = -z;
819:     R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] =     -y*z*inv1pX;
820:     R[6] = z; R[7] =     -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
821:   }
822:   else {
823:     PetscReal inv1mX   = 1./ (1. - x);

825:     R[0] = x; R[1] = z;               R[2] = y;
826:     R[3] = y; R[4] =     -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
827:     R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] =     -y*z*inv1mX;
828:   }
829:   coords[0] = 0.0;
830:   coords[1] = r;
831:   return(0);
832: }

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

837:   Not collective

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

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

846:   Level: developer

848: .seealso: DMPlexComputeProjection2Dto1D(), DMPlexComputeProjection3Dto1D()
849: @*/
850: PetscErrorCode DMPlexComputeProjection3Dto2D(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
851: {
852:   PetscReal      x1[3],  x2[3], n[3], norm;
853:   PetscReal      x1p[3], x2p[3], xnp[3];
854:   PetscReal      sqrtz, alpha;
855:   const PetscInt dim = 3;
856:   PetscInt       d, e, p;

859:   /* 0) Calculate normal vector */
860:   for (d = 0; d < dim; ++d) {
861:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
862:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
863:   }
864:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
865:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
866:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
867:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
868:   n[0] /= norm;
869:   n[1] /= norm;
870:   n[2] /= norm;
871:   /* 1) Take the normal vector and rotate until it is \hat z

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

875:     R = /  alpha nx nz  alpha ny nz -1/alpha \
876:         | -alpha ny     alpha nx        0    |
877:         \     nx            ny         nz    /

879:     will rotate the normal vector to \hat z
880:   */
881:   sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
882:   /* Check for n = z */
883:   if (sqrtz < 1.0e-10) {
884:     const PetscInt s = PetscSign(n[2]);
885:     /* If nz < 0, rotate 180 degrees around x-axis */
886:     for (p = 3; p < coordSize/3; ++p) {
887:       coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
888:       coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
889:     }
890:     coords[0] = 0.0;
891:     coords[1] = 0.0;
892:     coords[2] = x1[0];
893:     coords[3] = x1[1] * s;
894:     coords[4] = x2[0];
895:     coords[5] = x2[1] * s;
896:     R[0] = 1.0;     R[1] = 0.0;     R[2] = 0.0;
897:     R[3] = 0.0;     R[4] = 1.0 * s; R[5] = 0.0;
898:     R[6] = 0.0;     R[7] = 0.0;     R[8] = 1.0 * s;
899:     return(0);
900:   }
901:   alpha = 1.0/sqrtz;
902:   R[0] =  alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
903:   R[3] = -alpha*n[1];      R[4] = alpha*n[0];      R[5] = 0.0;
904:   R[6] =  n[0];            R[7] = n[1];            R[8] = n[2];
905:   for (d = 0; d < dim; ++d) {
906:     x1p[d] = 0.0;
907:     x2p[d] = 0.0;
908:     for (e = 0; e < dim; ++e) {
909:       x1p[d] += R[d*dim+e]*x1[e];
910:       x2p[d] += R[d*dim+e]*x2[e];
911:     }
912:   }
913:   if (PetscAbsReal(x1p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
914:   if (PetscAbsReal(x2p[2]) > 10. * PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
915:   /* 2) Project to (x, y) */
916:   for (p = 3; p < coordSize/3; ++p) {
917:     for (d = 0; d < dim; ++d) {
918:       xnp[d] = 0.0;
919:       for (e = 0; e < dim; ++e) {
920:         xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
921:       }
922:       if (d < dim-1) coords[p*2+d] = xnp[d];
923:     }
924:   }
925:   coords[0] = 0.0;
926:   coords[1] = 0.0;
927:   coords[2] = x1p[0];
928:   coords[3] = x1p[1];
929:   coords[4] = x2p[0];
930:   coords[5] = x2p[1];
931:   /* Output R^T which rotates \hat z to the input normal */
932:   for (d = 0; d < dim; ++d) {
933:     for (e = d+1; e < dim; ++e) {
934:       PetscReal tmp;

936:       tmp        = R[d*dim+e];
937:       R[d*dim+e] = R[e*dim+d];
938:       R[e*dim+d] = tmp;
939:     }
940:   }
941:   return(0);
942: }

944: PETSC_UNUSED
945: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
946: {
947:   /* Signed volume is 1/2 the determinant

949:    |  1  1  1 |
950:    | x0 x1 x2 |
951:    | y0 y1 y2 |

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

955:    | x1 x2 |
956:    | y1 y2 |
957:   */
958:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
959:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
960:   PetscReal       M[4], detM;
961:   M[0] = x1; M[1] = x2;
962:   M[2] = y1; M[3] = y2;
963:   DMPlex_Det2D_Internal(&detM, M);
964:   *vol = 0.5*detM;
965:   (void)PetscLogFlops(5.0);
966: }

968: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
969: {
970:   DMPlex_Det2D_Internal(vol, coords);
971:   *vol *= 0.5;
972: }

974: PETSC_UNUSED
975: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
976: {
977:   /* Signed volume is 1/6th of the determinant

979:    |  1  1  1  1 |
980:    | x0 x1 x2 x3 |
981:    | y0 y1 y2 y3 |
982:    | z0 z1 z2 z3 |

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

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

1003: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
1004: {
1005:   const PetscReal onesixth = ((PetscReal)1./(PetscReal)6.);
1006:   DMPlex_Det3D_Internal(vol, coords);
1007:   *vol *= -onesixth;
1008: }

1010: static PetscErrorCode DMPlexComputePointGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1011: {
1012:   PetscSection   coordSection;
1013:   Vec            coordinates;
1014:   const PetscScalar *coords;
1015:   PetscInt       dim, d, off;

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

1039: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1040: {
1041:   PetscSection   coordSection;
1042:   Vec            coordinates;
1043:   PetscScalar   *coords = NULL;
1044:   PetscInt       numCoords, d, pStart, pEnd, numSelfCoords = 0;

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

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

1074:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1075:     DMPlexComputeProjection2Dto1D(coords, R);
1076:     if (J)    {
1077:       J0   = 0.5*PetscRealPart(coords[1]);
1078:       J[0] = R[0]*J0; J[1] = R[1];
1079:       J[2] = R[2]*J0; J[3] = R[3];
1080:       DMPlex_Det2D_Internal(detJ, J);
1081:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1082:     }
1083:   } else if (numCoords == 2) {
1084:     const PetscInt dim = 1;

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

1098: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1099: {
1100:   PetscSection   coordSection;
1101:   Vec            coordinates;
1102:   PetscScalar   *coords = NULL;
1103:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

1107:   DMGetCoordinatesLocal(dm, &coordinates);
1108:   DMGetCoordinateSection(dm, &coordSection);
1109:   PetscSectionGetChart(coordSection,&pStart,&pEnd);
1110:   if (e >= pStart && e < pEnd) {PetscSectionGetDof(coordSection,e,&numSelfCoords);}
1111:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
1112:   numCoords = numSelfCoords ? numSelfCoords : numCoords;
1113:   *detJ = 0.0;
1114:   if (numCoords == 9) {
1115:     const PetscInt dim = 3;
1116:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

1118:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
1119:     DMPlexComputeProjection3Dto2D(numCoords, coords, R);
1120:     if (J)    {
1121:       const PetscInt pdim = 2;

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

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

1160: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscInt Nq, const PetscReal points[], PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1161: {
1162:   PetscSection   coordSection;
1163:   Vec            coordinates;
1164:   PetscScalar   *coords = NULL;
1165:   PetscInt       numCoords, numSelfCoords = 0, d, f, g, pStart, pEnd;

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

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

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

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

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

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

1246:       if (v) {
1247:         PetscReal extPoint[4];

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

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

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

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

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

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

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

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

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

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

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

1393:       if (v) {
1394:         PetscReal extPoint[8];

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

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

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

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

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

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

1455:   DMPlexGetDepth(dm, &depth);
1456:   DMPlexGetConeSize(dm, cell, &coneSize);
1457:   DMPlexGetDepthLabel(dm, &depthLabel);
1458:   DMLabelGetValue(depthLabel, cell, &dim);
1459:   if (depth == 1 && dim == 1) {
1460:     DMGetDimension(dm, &dim);
1461:   }
1462:   DMGetCoordinateDim(dm, &coordDim);
1463:   if (coordDim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported coordinate dimension %D > 3", coordDim);
1464:   if (quad) {PetscQuadratureGetData(quad, NULL, NULL, &Nq, &points, NULL);}
1465:   switch (dim) {
1466:   case 0:
1467:     DMPlexComputePointGeometry_Internal(dm, cell, v, J, invJ, detJ);
1468:     isAffine = PETSC_FALSE;
1469:     break;
1470:   case 1:
1471:     if (Nq) {
1472:       DMPlexComputeLineGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1473:     } else {
1474:       DMPlexComputeLineGeometry_Internal(dm, cell, v, J, invJ, detJ);
1475:     }
1476:     break;
1477:   case 2:
1478:     switch (coneSize) {
1479:     case 3:
1480:       if (Nq) {
1481:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1482:       } else {
1483:         DMPlexComputeTriangleGeometry_Internal(dm, cell, v, J, invJ, detJ);
1484:       }
1485:       break;
1486:     case 4:
1487:       DMPlexComputeRectangleGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1488:       isAffine = PETSC_FALSE;
1489:       break;
1490:     default:
1491:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1492:     }
1493:     break;
1494:   case 3:
1495:     switch (coneSize) {
1496:     case 4:
1497:       if (Nq) {
1498:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J0, NULL, &detJ0);
1499:       } else {
1500:         DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v, J, invJ, detJ);
1501:       }
1502:       break;
1503:     case 6: /* Faces */
1504:     case 8: /* Vertices */
1505:       DMPlexComputeHexahedronGeometry_Internal(dm, cell, Nq, points, v, J, invJ, detJ);
1506:       isAffine = PETSC_FALSE;
1507:       break;
1508:     default:
1509:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1510:     }
1511:     break;
1512:   default:
1513:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1514:   }
1515:   if (isAffine && Nq) {
1516:     if (v) {
1517:       for (i = 0; i < Nq; i++) {
1518:         CoordinatesRefToReal(coordDim, dim, xi0, v0, J0, &points[dim * i], &v[coordDim * i]);
1519:       }
1520:     }
1521:     if (detJ) {
1522:       for (i = 0; i < Nq; i++) {
1523:         detJ[i] = detJ0;
1524:       }
1525:     }
1526:     if (J) {
1527:       PetscInt k;

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

1532:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1533:           J[k] = J0[j];
1534:         }
1535:       }
1536:     }
1537:     if (invJ) {
1538:       PetscInt k;
1539:       switch (coordDim) {
1540:       case 0:
1541:         break;
1542:       case 1:
1543:         invJ[0] = 1./J0[0];
1544:         break;
1545:       case 2:
1546:         DMPlex_Invert2D_Internal(invJ, J0, detJ0);
1547:         break;
1548:       case 3:
1549:         DMPlex_Invert3D_Internal(invJ, J0, detJ0);
1550:         break;
1551:       }
1552:       for (i = 1, k = coordDim * coordDim; i < Nq; i++) {
1553:         PetscInt j;

1555:         for (j = 0; j < coordDim * coordDim; j++, k++) {
1556:           invJ[k] = invJ[j];
1557:         }
1558:       }
1559:     }
1560:   }
1561:   return(0);
1562: }

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

1567:   Collective on dm

1569:   Input Arguments:
1570: + dm   - the DM
1571: - cell - the cell

1573:   Output Arguments:
1574: + v0   - the translation part of this affine transform
1575: . J    - the Jacobian of the transform from the reference element
1576: . invJ - the inverse of the Jacobian
1577: - detJ - the Jacobian determinant

1579:   Level: advanced

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

1585: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinates()
1586: @*/
1587: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1588: {

1592:   DMPlexComputeCellGeometryFEM_Implicit(dm,cell,NULL,v0,J,invJ,detJ);
1593:   return(0);
1594: }

1596: static PetscErrorCode DMPlexComputeCellGeometryFEM_FE(DM dm, PetscFE fe, PetscInt point, PetscQuadrature quad, PetscReal v[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1597: {
1598:   PetscQuadrature  feQuad;
1599:   PetscSection     coordSection;
1600:   Vec              coordinates;
1601:   PetscScalar     *coords = NULL;
1602:   const PetscReal *quadPoints;
1603:   PetscReal       *basisDer, *basis, detJt;
1604:   PetscInt         dim, cdim, pdim, qdim, Nq, numCoords, q;
1605:   PetscErrorCode   ierr;

1608:   DMGetCoordinatesLocal(dm, &coordinates);
1609:   DMGetCoordinateSection(dm, &coordSection);
1610:   DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1611:   DMGetDimension(dm, &dim);
1612:   DMGetCoordinateDim(dm, &cdim);
1613:   if (!quad) { /* use the first point of the first functional of the dual space */
1614:     PetscDualSpace dsp;

1616:     PetscFEGetDualSpace(fe, &dsp);
1617:     PetscDualSpaceGetFunctional(dsp, 0, &quad);
1618:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1619:     Nq = 1;
1620:   } else {
1621:     PetscQuadratureGetData(quad, &qdim, NULL, &Nq, &quadPoints, NULL);
1622:   }
1623:   PetscFEGetDimension(fe, &pdim);
1624:   PetscFEGetQuadrature(fe, &feQuad);
1625:   if (feQuad == quad) {
1626:     PetscFEGetDefaultTabulation(fe, &basis, J ? &basisDer : NULL, NULL);
1627:     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);
1628:   } else {
1629:     PetscFEGetTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1630:   }
1631:   if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1632:   if (v) {
1633:     PetscMemzero(v, Nq*cdim*sizeof(PetscReal));
1634:     for (q = 0; q < Nq; ++q) {
1635:       PetscInt i, k;

1637:       for (k = 0; k < pdim; ++k)
1638:         for (i = 0; i < cdim; ++i)
1639:           v[q*cdim + i] += basis[q*pdim + k] * PetscRealPart(coords[k*cdim + i]);
1640:       PetscLogFlops(2.0*pdim*cdim);
1641:     }
1642:   }
1643:   if (J) {
1644:     PetscMemzero(J, Nq*cdim*cdim*sizeof(PetscReal));
1645:     for (q = 0; q < Nq; ++q) {
1646:       PetscInt i, j, k, c, r;

1648:       /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1649:       for (k = 0; k < pdim; ++k)
1650:         for (j = 0; j < dim; ++j)
1651:           for (i = 0; i < cdim; ++i)
1652:             J[(q*cdim + i)*cdim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1653:       PetscLogFlops(2.0*pdim*dim*cdim);
1654:       if (cdim > dim) {
1655:         for (c = dim; c < cdim; ++c)
1656:           for (r = 0; r < cdim; ++r)
1657:             J[r*cdim+c] = r == c ? 1.0 : 0.0;
1658:       }
1659:       if (!detJ && !invJ) continue;
1660:       detJt = 0.;
1661:       switch (cdim) {
1662:       case 3:
1663:         DMPlex_Det3D_Internal(&detJt, &J[q*cdim*dim]);
1664:         if (invJ) {DMPlex_Invert3D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1665:         break;
1666:       case 2:
1667:         DMPlex_Det2D_Internal(&detJt, &J[q*cdim*dim]);
1668:         if (invJ) {DMPlex_Invert2D_Internal(&invJ[q*cdim*dim], &J[q*cdim*dim], detJt);}
1669:         break;
1670:       case 1:
1671:         detJt = J[q*cdim*dim];
1672:         if (invJ) invJ[q*cdim*dim] = 1.0/detJt;
1673:       }
1674:       if (detJ) detJ[q] = detJt;
1675:     }
1676:   }
1677:   else if (detJ || invJ) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Need J to compute invJ or detJ");
1678:   if (feQuad != quad) {
1679:     PetscFERestoreTabulation(fe, Nq, quadPoints, &basis, J ? &basisDer : NULL, NULL);
1680:   }
1681:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1682:   return(0);
1683: }

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

1688:   Collective on dm

1690:   Input Arguments:
1691: + dm   - the DM
1692: . cell - the cell
1693: - quad - the quadrature containing the points in the reference element where the geometry will be evaluated.  If quad == NULL, geometry will be
1694:          evaluated at the first vertex of the reference element

1696:   Output Arguments:
1697: + v    - the image of the transformed quadrature points, otherwise the image of the first vertex in the closure of the reference element
1698: . J    - the Jacobian of the transform from the reference element at each quadrature point
1699: . invJ - the inverse of the Jacobian at each quadrature point
1700: - detJ - the Jacobian determinant at each quadrature point

1702:   Level: advanced

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

1708: .seealso: DMGetCoordinateSection(), DMGetCoordinates()
1709: @*/
1710: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal *v, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1711: {
1712:   DM             cdm;
1713:   PetscFE        fe = NULL;

1718:   DMGetCoordinateDM(dm, &cdm);
1719:   if (cdm) {
1720:     PetscClassId id;
1721:     PetscInt     numFields;
1722:     PetscDS      prob;
1723:     PetscObject  disc;

1725:     DMGetNumFields(cdm, &numFields);
1726:     if (numFields) {
1727:       DMGetDS(cdm, &prob);
1728:       PetscDSGetDiscretization(prob,0,&disc);
1729:       PetscObjectGetClassId(disc,&id);
1730:       if (id == PETSCFE_CLASSID) {
1731:         fe = (PetscFE) disc;
1732:       }
1733:     }
1734:   }
1735:   if (!fe) {DMPlexComputeCellGeometryFEM_Implicit(dm, cell, quad, v, J, invJ, detJ);}
1736:   else     {DMPlexComputeCellGeometryFEM_FE(dm, fe, cell, quad, v, J, invJ, detJ);}
1737:   return(0);
1738: }

1740: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1741: {
1742:   PetscSection   coordSection;
1743:   Vec            coordinates;
1744:   PetscScalar   *coords = NULL;
1745:   PetscScalar    tmp[2];
1746:   PetscInt       coordSize;

1750:   DMGetCoordinatesLocal(dm, &coordinates);
1751:   DMGetCoordinateSection(dm, &coordSection);
1752:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1753:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
1754:   DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1755:   if (centroid) {
1756:     centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]);
1757:     centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]);
1758:   }
1759:   if (normal) {
1760:     PetscReal norm;

1762:     normal[0]  = -PetscRealPart(coords[1] - tmp[1]);
1763:     normal[1]  =  PetscRealPart(coords[0] - tmp[0]);
1764:     norm       = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
1765:     normal[0] /= norm;
1766:     normal[1] /= norm;
1767:   }
1768:   if (vol) {
1769:     *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1])));
1770:   }
1771:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1772:   return(0);
1773: }

1775: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1776: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1777: {
1778:   DMLabel        depth;
1779:   PetscSection   coordSection;
1780:   Vec            coordinates;
1781:   PetscScalar   *coords = NULL;
1782:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1783:   PetscBool      isHybrid = PETSC_FALSE;
1784:   PetscInt       fv[4] = {0, 1, 2, 3};
1785:   PetscInt       cdepth, pEndInterior[4], tdim = 2, coordSize, numCorners, p, d, e;

1789:   /* Must check for hybrid cells because prisms have a different orientation scheme */
1790:   DMPlexGetDepthLabel(dm, &depth);
1791:   DMLabelGetValue(depth, cell, &cdepth);
1792:   DMPlexGetHybridBounds(dm, &pEndInterior[3], &pEndInterior[2], &pEndInterior[1], &pEndInterior[0]);
1793:   if ((pEndInterior[cdepth]) >=0 && (cell >= pEndInterior[cdepth])) isHybrid = PETSC_TRUE;

1795:   DMGetCoordinatesLocal(dm, &coordinates);
1796:   DMPlexGetConeSize(dm, cell, &numCorners);
1797:   DMGetCoordinateSection(dm, &coordSection);
1798:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1799:   DMGetCoordinateDim(dm, &dim);
1800:   /* Side faces for hybrid cells are are stored as tensor products */
1801:   if (isHybrid && numCorners == 4) {fv[2] = 3; fv[3] = 2;}

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

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

1859: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1860: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1861: {
1862:   DMLabel         depth;
1863:   PetscSection    coordSection;
1864:   Vec             coordinates;
1865:   PetscScalar    *coords = NULL;
1866:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
1867:   const PetscInt *faces, *facesO;
1868:   PetscBool       isHybrid = PETSC_FALSE;
1869:   PetscInt        pEndInterior[4], cdepth, numFaces, f, coordSize, numCorners, p, d;
1870:   PetscErrorCode  ierr;

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

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

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

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

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

1956: /*@C
1957:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

1959:   Collective on dm

1961:   Input Arguments:
1962: + dm   - the DM
1963: - cell - the cell

1965:   Output Arguments:
1966: + volume   - the cell volume
1967: . centroid - the cell centroid
1968: - normal - the cell normal, if appropriate

1970:   Level: advanced

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

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

1984:   DMPlexGetDepth(dm, &depth);
1985:   DMGetDimension(dm, &dim);
1986:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1987:   /* We need to keep a pointer to the depth label */
1988:   DMGetLabelValue(dm, "depth", cell, &depth);
1989:   /* Cone size is now the number of faces */
1990:   switch (depth) {
1991:   case 1:
1992:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1993:     break;
1994:   case 2:
1995:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1996:     break;
1997:   case 3:
1998:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1999:     break;
2000:   default:
2001:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D (depth %D) for element geometry computation", dim, depth);
2002:   }
2003:   return(0);
2004: }

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

2009:   Collective on dm

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

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

2017:   Level: beginner

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

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

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

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

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

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

2068:   Level: developer

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

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

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

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

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

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

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

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

2217:   Not collective

2219:   Input Argument:
2220: . dm - the DM

2222:   Output Argument:
2223: . minradius - the minium cell radius

2225:   Level: developer

2227: .seealso: DMGetCoordinates()
2228: @*/
2229: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
2230: {
2234:   *minradius = ((DM_Plex*) dm->data)->minradius;
2235:   return(0);
2236: }

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

2241:   Logically collective

2243:   Input Arguments:
2244: + dm - the DM
2245: - minradius - the minium cell radius

2247:   Level: developer

2249: .seealso: DMSetCoordinates()
2250: @*/
2251: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
2252: {
2255:   ((DM_Plex*) dm->data)->minradius = minradius;
2256:   return(0);
2257: }

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

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

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

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

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

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

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

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

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

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

2397:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
2398:     PetscSectionGetDof(neighSec, c, &numFaces);
2399:     PetscSectionGetOffset(neighSec, c, &off);
2400:     if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
2401:     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);
2402:     for (f = 0; f < numFaces; ++f) {
2403:       PetscFVCellGeom       *cg1;
2404:       PetscFVFaceGeom       *fg;
2405:       const PetscInt        *fcells;
2406:       PetscInt               ncell, side, nface;

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

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

2431:   Collective on dm

2433:   Input Arguments:
2434: + dm  - The DM
2435: . fvm - The PetscFV
2436: . faceGeometry - The face geometry from DMPlexComputeFaceGeometryFVM()
2437: - cellGeometry - The face geometry from DMPlexComputeCellGeometryFVM()

2439:   Output Parameters:
2440: + faceGeometry - The geometric factors for gradient calculation are inserted
2441: - dmGrad - The DM describing the layout of gradient data

2443:   Level: developer

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

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

2484: /*@
2485:   DMPlexGetDataFVM - Retrieve precomputed cell geometry

2487:   Collective on dm

2489:   Input Arguments:
2490: + dm  - The DM
2491: - fvm - The PetscFV

2493:   Output Parameters:
2494: + cellGeometry - The cell geometry
2495: . faceGeometry - The face geometry
2496: - dmGrad       - The gradient matrices

2498:   Level: developer

2500: .seealso: DMPlexComputeGeometryFVM()
2501: @*/
2502: PetscErrorCode DMPlexGetDataFVM(DM dm, PetscFV fv, Vec *cellgeom, Vec *facegeom, DM *gradDM)
2503: {
2504:   PetscObject    cellgeomobj, facegeomobj;

2508:   PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2509:   if (!cellgeomobj) {
2510:     Vec cellgeomInt, facegeomInt;

2512:     DMPlexComputeGeometryFVM(dm, &cellgeomInt, &facegeomInt);
2513:     PetscObjectCompose((PetscObject) dm, "DMPlex_cellgeom_fvm",(PetscObject)cellgeomInt);
2514:     PetscObjectCompose((PetscObject) dm, "DMPlex_facegeom_fvm",(PetscObject)facegeomInt);
2515:     VecDestroy(&cellgeomInt);
2516:     VecDestroy(&facegeomInt);
2517:     PetscObjectQuery((PetscObject) dm, "DMPlex_cellgeom_fvm", &cellgeomobj);
2518:   }
2519:   PetscObjectQuery((PetscObject) dm, "DMPlex_facegeom_fvm", &facegeomobj);
2520:   if (cellgeom) *cellgeom = (Vec) cellgeomobj;
2521:   if (facegeom) *facegeom = (Vec) facegeomobj;
2522:   if (gradDM) {
2523:     PetscObject gradobj;
2524:     PetscBool   computeGradients;

2526:     PetscFVGetComputeGradients(fv,&computeGradients);
2527:     if (!computeGradients) {
2528:       *gradDM = NULL;
2529:       return(0);
2530:     }
2531:     PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2532:     if (!gradobj) {
2533:       DM dmGradInt;

2535:       DMPlexComputeGradientFVM(dm,fv,(Vec) facegeomobj,(Vec) cellgeomobj,&dmGradInt);
2536:       PetscObjectCompose((PetscObject) dm, "DMPlex_dmgrad_fvm", (PetscObject)dmGradInt);
2537:       DMDestroy(&dmGradInt);
2538:       PetscObjectQuery((PetscObject) dm, "DMPlex_dmgrad_fvm", &gradobj);
2539:     }
2540:     *gradDM = (DM) gradobj;
2541:   }
2542:   return(0);
2543: }

2545: static PetscErrorCode DMPlexCoordinatesToReference_NewtonUpdate(PetscInt dimC, PetscInt dimR, PetscScalar *J, PetscScalar *invJ, PetscScalar *work,  PetscReal *resNeg, PetscReal *guess)
2546: {
2547:   PetscInt l, m;

2550:   if (dimC == dimR && dimR <= 3) {
2551:     /* invert Jacobian, multiply */
2552:     PetscScalar det, idet;

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

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

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

2606:     for (l = 0; l < dimR; l++) {guess[l] += PetscRealPart(invJ[l]);}
2607:   }
2608:   return(0);
2609: }

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

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

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

2637:       for (j = 0; j < dimC; j++) {
2638:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2639:       }
2640:     }
2641:   } else if (dimR == 3) {
2642:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

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

2658:     for (j = 0; j < (numV / 2); j++) {
2659:       for (k = 0; k < dimC; k++) {
2660:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2661:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2662:       }
2663:     }

2665:     if (i < dimR - 1) {
2666:       swap = cellCoeffs;
2667:       cellCoeffs = cellCoords;
2668:       cellCoords = swap;
2669:     }
2670:   }
2671:   PetscMemzero(refCoords,numPoints * dimR * sizeof (PetscReal));
2672:   for (j = 0; j < numPoints; j++) {
2673:     for (i = 0; i < maxIts; i++) {
2674:       PetscReal *guess = &refCoords[dimR * j];

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

2685:           extCoord *= dep * coord + !dep;
2686:           extJ[l] = dep;

2688:           for (m = 0; m < dimR; m++) {
2689:             PetscReal coord = guess[m];
2690:             PetscInt  dep   = ((k & (1 << m)) >> m) && (m != l);
2691:             PetscReal mult  = dep * coord + !dep;

2693:             extJ[l] *= mult;
2694:           }
2695:         }
2696:         for (l = 0; l < dimC; l++) {
2697:           PetscReal coeff = cellCoeffs[dimC * k + l];

2699:           resNeg[l] -= coeff * extCoord;
2700:           for (m = 0; m < dimR; m++) {
2701:             J[dimR * l + m] += coeff * extJ[m];
2702:           }
2703:         }
2704:       }
2705: #if 0 && defined(PETSC_USE_DEBUG)
2706:       {
2707:         PetscReal maxAbs = 0.;

2709:         for (l = 0; l < dimC; l++) {
2710:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2711:         }
2712:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);
2713:       }
2714: #endif

2716:       DMPlexCoordinatesToReference_NewtonUpdate(dimC,dimR,J,invJ,work,resNeg,guess);
2717:     }
2718:   }
2719:   DMRestoreWorkArray(dm, 3 * dimR * dimC, MPIU_SCALAR, &J);
2720:   DMRestoreWorkArray(dm, 2 * coordSize + dimR + dimC, MPIU_REAL, &cellData);
2721:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2722:   return(0);
2723: }

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

2734:   DMPlexVecGetClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2735:   if (coordSize < dimC * numV) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Expecting at least %D coordinates, got %D",dimC * (1 << dimR), coordSize);
2736:   DMGetWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2737:   cellCoords = &cellData[0];
2738:   cellCoeffs = &cellData[coordSize];
2739:   if (dimR == 2) {
2740:     const PetscInt zToPlex[4] = {0, 1, 3, 2};

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

2745:       for (j = 0; j < dimC; j++) {
2746:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2747:       }
2748:     }
2749:   } else if (dimR == 3) {
2750:     const PetscInt zToPlex[8] = {0, 3, 1, 2, 4, 5, 7, 6};

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

2755:       for (j = 0; j < dimC; j++) {
2756:         cellCoords[dimC * i + j] = PetscRealPart(coordsScalar[dimC * plexI + j]);
2757:       }
2758:     }
2759:   } else {
2760:     for (i = 0; i < coordSize; i++) {cellCoords[i] = PetscRealPart(coordsScalar[i]);}
2761:   }
2762:   /* Perform the shuffling transform that converts values at the corners of [-1,1]^d to coefficients */
2763:   for (i = 0; i < dimR; i++) {
2764:     PetscReal *swap;

2766:     for (j = 0; j < (numV / 2); j++) {
2767:       for (k = 0; k < dimC; k++) {
2768:         cellCoeffs[dimC * j + k]                = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] + cellCoords[dimC * 2 * j + k]);
2769:         cellCoeffs[dimC * (j + (numV / 2)) + k] = 0.5 * (cellCoords[dimC * (2 * j + 1) + k] - cellCoords[dimC * 2 * j + k]);
2770:       }
2771:     }

2773:     if (i < dimR - 1) {
2774:       swap = cellCoeffs;
2775:       cellCoeffs = cellCoords;
2776:       cellCoords = swap;
2777:     }
2778:   }
2779:   PetscMemzero(realCoords,numPoints * dimC * sizeof (PetscReal));
2780:   for (j = 0; j < numPoints; j++) {
2781:     const PetscReal *guess  = &refCoords[dimR * j];
2782:     PetscReal       *mapped = &realCoords[dimC * j];

2784:     for (k = 0; k < numV; k++) {
2785:       PetscReal extCoord = 1.;
2786:       for (l = 0; l < dimR; l++) {
2787:         PetscReal coord = guess[l];
2788:         PetscInt  dep   = (k & (1 << l)) >> l;

2790:         extCoord *= dep * coord + !dep;
2791:       }
2792:       for (l = 0; l < dimC; l++) {
2793:         PetscReal coeff = cellCoeffs[dimC * k + l];

2795:         mapped[l] += coeff * extCoord;
2796:       }
2797:     }
2798:   }
2799:   DMRestoreWorkArray(dm, 2 * coordSize, MPIU_REAL, &cellData);
2800:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &coordsScalar);
2801:   return(0);
2802: }

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

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

2853:         for (l = 0; l < Nc; l++) {
2854:           maxAbs = PetscMax(maxAbs,PetscAbsReal(resNeg[l]));
2855:         }
2856:         PetscInfo4(dm,"cell %D, point %D, iter %D: res %g\n",cell,j,i,maxAbs);
2857:       }
2858: #endif
2859:       DMPlexCoordinatesToReference_NewtonUpdate(Nc,dimR,J,invJ,work,resNeg,guess);
2860:     }
2861:   }
2862:   DMRestoreWorkArray(dm,3 * Nc * dimR,MPIU_SCALAR,&J);
2863:   DMRestoreWorkArray(dm,pdim * Nc + pdim * Nc * dimR + Nc,MPIU_REAL,&B);
2864:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2865:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2866:   return(0);
2867: }

2869: /* TODO: TOBY please fix this for Nc > 1 */
2870: static PetscErrorCode DMPlexReferenceToCoordinates_FE(DM dm, PetscFE fe, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[], Vec coords, PetscInt Nc, PetscInt dimR)
2871: {
2872:   PetscInt       numComp, pdim, i, j, k, l, coordSize;
2873:   PetscScalar    *nodes = NULL;
2874:   PetscReal      *invV, *modes;
2875:   PetscReal      *B;

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

2898:     for (k = 0; k < pdim; k++) {
2899:       for (l = 0; l < Nc; l++) {
2900:         mapped[l] += modes[k] * B[(j * pdim + k) * Nc + l];
2901:       }
2902:     }
2903:   }
2904:   DMRestoreWorkArray(dm,numPoints * pdim * Nc,MPIU_REAL,&B);
2905:   DMRestoreWorkArray(dm,pdim,MPIU_REAL,&modes);
2906:   DMPlexVecRestoreClosure(dm, NULL, coords, cell, &coordSize, &nodes);
2907:   return(0);
2908: }

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

2915:   Not collective

2917:   Input Parameters:
2918: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
2919:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
2920:                as a multilinear map for tensor-product elements
2921: . cell       - the cell whose map is used.
2922: . numPoints  - the number of points to locate
2923: - realCoords - (numPoints x coordinate dimension) array of coordinates (see DMGetCoordinateDim())

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

2928:   Level: intermediate

2930: .seealso: DMPlexReferenceToCoordinates()
2931: @*/
2932: PetscErrorCode DMPlexCoordinatesToReference(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal realCoords[], PetscReal refCoords[])
2933: {
2934:   PetscInt       dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
2935:   DM             coordDM = NULL;
2936:   Vec            coords;
2937:   PetscFE        fe = NULL;

2942:   DMGetDimension(dm,&dimR);
2943:   DMGetCoordinateDim(dm,&dimC);
2944:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
2945:   DMPlexGetDepth(dm,&depth);
2946:   DMGetCoordinatesLocal(dm,&coords);
2947:   DMGetCoordinateDM(dm,&coordDM);
2948:   if (coordDM) {
2949:     PetscInt coordFields;

2951:     DMGetNumFields(coordDM,&coordFields);
2952:     if (coordFields) {
2953:       PetscClassId id;
2954:       PetscObject  disc;

2956:       DMGetField(coordDM,0,NULL,&disc);
2957:       PetscObjectGetClassId(disc,&id);
2958:       if (id == PETSCFE_CLASSID) {
2959:         fe = (PetscFE) disc;
2960:       }
2961:     }
2962:   }
2963:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
2964:   DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
2965:   cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
2966:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
2967:   if (!fe) { /* implicit discretization: affine or multilinear */
2968:     PetscInt  coneSize;
2969:     PetscBool isSimplex, isTensor;

2971:     DMPlexGetConeSize(dm,cell,&coneSize);
2972:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
2973:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
2974:     if (isSimplex) {
2975:       PetscReal detJ, *v0, *J, *invJ;

2977:       DMGetWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2978:       J    = &v0[dimC];
2979:       invJ = &J[dimC * dimC];
2980:       DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, &detJ);
2981:       for (i = 0; i < numPoints; i++) { /* Apply the inverse affine transformation for each point */
2982:         const PetscReal x0[3] = {-1.,-1.,-1.};

2984:         CoordinatesRealToRef(dimC, dimR, x0, v0, invJ, &realCoords[dimC * i], &refCoords[dimR * i]);
2985:       }
2986:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
2987:     } else if (isTensor) {
2988:       DMPlexCoordinatesToReference_Tensor(coordDM, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2989:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
2990:   } else {
2991:     DMPlexCoordinatesToReference_FE(coordDM, fe, cell, numPoints, realCoords, refCoords, coords, dimC, dimR);
2992:   }
2993:   return(0);
2994: }

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

2999:   Not collective

3001:   Input Parameters:
3002: + dm         - The mesh, with coordinate maps defined either by a PetscDS for the coordinate DM (see DMGetCoordinateDM()) or
3003:                implicitly by the coordinates of the corner vertices of the cell: as an affine map for simplicial elements, or
3004:                as a multilinear map for tensor-product elements
3005: . cell       - the cell whose map is used.
3006: . numPoints  - the number of points to locate
3007: + refCoords  - (numPoints x dimension) array of reference coordinates (see DMGetDimension())

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

3012:    Level: intermediate
3013:    
3014: .seealso: DMPlexCoordinatesToReference()
3015: @*/
3016: PetscErrorCode DMPlexReferenceToCoordinates(DM dm, PetscInt cell, PetscInt numPoints, const PetscReal refCoords[], PetscReal realCoords[])
3017: {
3018:   PetscInt       dimC, dimR, depth, cStart, cEnd, cEndInterior, i;
3019:   DM             coordDM = NULL;
3020:   Vec            coords;
3021:   PetscFE        fe = NULL;

3026:   DMGetDimension(dm,&dimR);
3027:   DMGetCoordinateDim(dm,&dimC);
3028:   if (dimR <= 0 || dimC <= 0 || numPoints <= 0) return(0);
3029:   DMPlexGetDepth(dm,&depth);
3030:   DMGetCoordinatesLocal(dm,&coords);
3031:   DMGetCoordinateDM(dm,&coordDM);
3032:   if (coordDM) {
3033:     PetscInt coordFields;

3035:     DMGetNumFields(coordDM,&coordFields);
3036:     if (coordFields) {
3037:       PetscClassId id;
3038:       PetscObject  disc;

3040:       DMGetField(coordDM,0,NULL,&disc);
3041:       PetscObjectGetClassId(disc,&id);
3042:       if (id == PETSCFE_CLASSID) {
3043:         fe = (PetscFE) disc;
3044:       }
3045:     }
3046:   }
3047:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
3048:   DMPlexGetHybridBounds(dm,&cEndInterior,NULL,NULL,NULL);
3049:   cEnd = cEndInterior > 0 ? cEndInterior : cEnd;
3050:   if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"point %D not in cell range [%D,%D)",cell,cStart,cEnd);
3051:   if (!fe) { /* implicit discretization: affine or multilinear */
3052:     PetscInt  coneSize;
3053:     PetscBool isSimplex, isTensor;

3055:     DMPlexGetConeSize(dm,cell,&coneSize);
3056:     isSimplex = (coneSize == (dimR + 1)) ? PETSC_TRUE : PETSC_FALSE;
3057:     isTensor  = (coneSize == ((depth == 1) ? (1 << dimR) : (2 * dimR))) ? PETSC_TRUE : PETSC_FALSE;
3058:     if (isSimplex) {
3059:       PetscReal detJ, *v0, *J;

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

3067:         CoordinatesRefToReal(dimC, dimR, xi0, v0, J, &refCoords[dimR * i], &realCoords[dimC * i]);
3068:       }
3069:       DMRestoreWorkArray(dm,dimC + 2 * dimC * dimC, MPIU_REAL, &v0);
3070:     } else if (isTensor) {
3071:       DMPlexReferenceToCoordinates_Tensor(coordDM, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3072:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unrecognized cone size %D",coneSize);
3073:   } else {
3074:     DMPlexReferenceToCoordinates_FE(coordDM, fe, cell, numPoints, refCoords, realCoords, coords, dimC, dimR);
3075:   }
3076:   return(0);
3077: }