Actual source code: plexgeometry.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/

  5: static PetscErrorCode DMPlexGetLineIntersection_2D_Internal(const PetscReal segmentA[], const PetscReal segmentB[], PetscReal intersection[], PetscBool *hasIntersection)
  6: {
  7:   const PetscReal p0_x  = segmentA[0*2+0];
  8:   const PetscReal p0_y  = segmentA[0*2+1];
  9:   const PetscReal p1_x  = segmentA[1*2+0];
 10:   const PetscReal p1_y  = segmentA[1*2+1];
 11:   const PetscReal p2_x  = segmentB[0*2+0];
 12:   const PetscReal p2_y  = segmentB[0*2+1];
 13:   const PetscReal p3_x  = segmentB[1*2+0];
 14:   const PetscReal p3_y  = segmentB[1*2+1];
 15:   const PetscReal s1_x  = p1_x - p0_x;
 16:   const PetscReal s1_y  = p1_y - p0_y;
 17:   const PetscReal s2_x  = p3_x - p2_x;
 18:   const PetscReal s2_y  = p3_y - p2_y;
 19:   const PetscReal denom = (-s2_x * s1_y + s1_x * s2_y);

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

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

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

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

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

 63: static PetscErrorCode DMPlexLocatePoint_General_2D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 64: {
 65:   PetscSection       coordSection;
 66:   Vec             coordsLocal;
 67:   PetscScalar    *coords = NULL;
 68:   const PetscInt  faces[8]  = {0, 1, 1, 2, 2, 3, 3, 0};
 69:   PetscReal       x         = PetscRealPart(point[0]);
 70:   PetscReal       y         = PetscRealPart(point[1]);
 71:   PetscInt        crossings = 0, f;
 72:   PetscErrorCode  ierr;

 75:   DMGetCoordinatesLocal(dm, &coordsLocal);
 76:   DMGetCoordinateSection(dm, &coordSection);
 77:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
 78:   for (f = 0; f < 4; ++f) {
 79:     PetscReal x_i   = PetscRealPart(coords[faces[2*f+0]*2+0]);
 80:     PetscReal y_i   = PetscRealPart(coords[faces[2*f+0]*2+1]);
 81:     PetscReal x_j   = PetscRealPart(coords[faces[2*f+1]*2+0]);
 82:     PetscReal y_j   = PetscRealPart(coords[faces[2*f+1]*2+1]);
 83:     PetscReal slope = (y_j - y_i) / (x_j - x_i);
 84:     PetscBool cond1 = (x_i <= x) && (x < x_j) ? PETSC_TRUE : PETSC_FALSE;
 85:     PetscBool cond2 = (x_j <= x) && (x < x_i) ? PETSC_TRUE : PETSC_FALSE;
 86:     PetscBool above = (y < slope * (x - x_i) + y_i) ? PETSC_TRUE : PETSC_FALSE;
 87:     if ((cond1 || cond2)  && above) ++crossings;
 88:   }
 89:   if (crossings % 2) *cell = c;
 90:   else *cell = -1;
 91:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
 92:   return(0);
 93: }

 97: static PetscErrorCode DMPlexLocatePoint_Simplex_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
 98: {
 99:   const PetscInt embedDim = 3;
100:   PetscReal      v0[3], J[9], invJ[9], detJ;
101:   PetscReal      x = PetscRealPart(point[0]);
102:   PetscReal      y = PetscRealPart(point[1]);
103:   PetscReal      z = PetscRealPart(point[2]);
104:   PetscReal      xi, eta, zeta;

108:   DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);
109:   xi   = invJ[0*embedDim+0]*(x - v0[0]) + invJ[0*embedDim+1]*(y - v0[1]) + invJ[0*embedDim+2]*(z - v0[2]);
110:   eta  = invJ[1*embedDim+0]*(x - v0[0]) + invJ[1*embedDim+1]*(y - v0[1]) + invJ[1*embedDim+2]*(z - v0[2]);
111:   zeta = invJ[2*embedDim+0]*(x - v0[0]) + invJ[2*embedDim+1]*(y - v0[1]) + invJ[2*embedDim+2]*(z - v0[2]);

113:   if ((xi >= 0.0) && (eta >= 0.0) && (zeta >= 0.0) && (xi + eta + zeta <= 2.0)) *cell = c;
114:   else *cell = -1;
115:   return(0);
116: }

120: static PetscErrorCode DMPlexLocatePoint_General_3D_Internal(DM dm, const PetscScalar point[], PetscInt c, PetscInt *cell)
121: {
122:   PetscSection   coordSection;
123:   Vec            coordsLocal;
124:   PetscScalar   *coords;
125:   const PetscInt faces[24] = {0, 3, 2, 1,  5, 4, 7, 6,  3, 0, 4, 5,
126:                               1, 2, 6, 7,  3, 5, 6, 2,  0, 1, 7, 4};
127:   PetscBool      found = PETSC_TRUE;
128:   PetscInt       f;

132:   DMGetCoordinatesLocal(dm, &coordsLocal);
133:   DMGetCoordinateSection(dm, &coordSection);
134:   DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
135:   for (f = 0; f < 6; ++f) {
136:     /* Check the point is under plane */
137:     /*   Get face normal */
138:     PetscReal v_i[3];
139:     PetscReal v_j[3];
140:     PetscReal normal[3];
141:     PetscReal pp[3];
142:     PetscReal dot;

144:     v_i[0]    = PetscRealPart(coords[faces[f*4+3]*3+0]-coords[faces[f*4+0]*3+0]);
145:     v_i[1]    = PetscRealPart(coords[faces[f*4+3]*3+1]-coords[faces[f*4+0]*3+1]);
146:     v_i[2]    = PetscRealPart(coords[faces[f*4+3]*3+2]-coords[faces[f*4+0]*3+2]);
147:     v_j[0]    = PetscRealPart(coords[faces[f*4+1]*3+0]-coords[faces[f*4+0]*3+0]);
148:     v_j[1]    = PetscRealPart(coords[faces[f*4+1]*3+1]-coords[faces[f*4+0]*3+1]);
149:     v_j[2]    = PetscRealPart(coords[faces[f*4+1]*3+2]-coords[faces[f*4+0]*3+2]);
150:     normal[0] = v_i[1]*v_j[2] - v_i[2]*v_j[1];
151:     normal[1] = v_i[2]*v_j[0] - v_i[0]*v_j[2];
152:     normal[2] = v_i[0]*v_j[1] - v_i[1]*v_j[0];
153:     pp[0]     = PetscRealPart(coords[faces[f*4+0]*3+0] - point[0]);
154:     pp[1]     = PetscRealPart(coords[faces[f*4+0]*3+1] - point[1]);
155:     pp[2]     = PetscRealPart(coords[faces[f*4+0]*3+2] - point[2]);
156:     dot       = normal[0]*pp[0] + normal[1]*pp[1] + normal[2]*pp[2];

158:     /* Check that projected point is in face (2D location problem) */
159:     if (dot < 0.0) {
160:       found = PETSC_FALSE;
161:       break;
162:     }
163:   }
164:   if (found) *cell = c;
165:   else *cell = -1;
166:   DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &coords);
167:   return(0);
168: }

172: static PetscErrorCode PetscGridHashInitialize_Internal(PetscGridHash box, PetscInt dim, const PetscScalar point[])
173: {
174:   PetscInt d;

177:   box->dim = dim;
178:   for (d = 0; d < dim; ++d) box->lower[d] = box->upper[d] = PetscRealPart(point[d]);
179:   return(0);
180: }

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

189:   PetscMalloc1(1, box);
190:   PetscGridHashInitialize_Internal(*box, dim, point);
191:   return(0);
192: }

196: PetscErrorCode PetscGridHashEnlarge(PetscGridHash box, const PetscScalar point[])
197: {
198:   PetscInt d;

201:   for (d = 0; d < box->dim; ++d) {
202:     box->lower[d] = PetscMin(box->lower[d], PetscRealPart(point[d]));
203:     box->upper[d] = PetscMax(box->upper[d], PetscRealPart(point[d]));
204:   }
205:   return(0);
206: }

210: PetscErrorCode PetscGridHashSetGrid(PetscGridHash box, const PetscInt n[], const PetscReal h[])
211: {
212:   PetscInt d;

215:   for (d = 0; d < box->dim; ++d) {
216:     box->extent[d] = box->upper[d] - box->lower[d];
217:     if (n[d] == PETSC_DETERMINE) {
218:       box->h[d] = h[d];
219:       box->n[d] = PetscCeilReal(box->extent[d]/h[d]);
220:     } else {
221:       box->n[d] = n[d];
222:       box->h[d] = box->extent[d]/n[d];
223:     }
224:   }
225:   return(0);
226: }

230: PetscErrorCode PetscGridHashGetEnclosingBox(PetscGridHash box, PetscInt numPoints, const PetscScalar points[], PetscInt dboxes[], PetscInt boxes[])
231: {
232:   const PetscReal *lower = box->lower;
233:   const PetscReal *upper = box->upper;
234:   const PetscReal *h     = box->h;
235:   const PetscInt  *n     = box->n;
236:   const PetscInt   dim   = box->dim;
237:   PetscInt         d, p;

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

244:       if (dbox == n[d] && PetscAbsReal(PetscRealPart(points[p*dim+d]) - upper[d]) < 1.0e-9) dbox = n[d]-1;
245:       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",
246:                                              p, PetscRealPart(points[p*dim+0]), dim > 1 ? PetscRealPart(points[p*dim+1]) : 0.0, dim > 2 ? PetscRealPart(points[p*dim+2]) : 0.0);
247:       dboxes[p*dim+d] = dbox;
248:     }
249:     if (boxes) for (d = 1, boxes[p] = dboxes[p*dim]; d < dim; ++d) boxes[p] += dboxes[p*dim+d]*n[d-1];
250:   }
251:   return(0);
252: }

256: PetscErrorCode PetscGridHashDestroy(PetscGridHash *box)
257: {

261:   if (*box) {
262:     PetscSectionDestroy(&(*box)->cellSection);
263:     ISDestroy(&(*box)->cells);
264:     DMLabelDestroy(&(*box)->cellsSparse);
265:   }
266:   PetscFree(*box);
267:   return(0);
268: }

272: PetscErrorCode DMPlexLocatePoint_Internal(DM dm, PetscInt dim, const PetscScalar point[], PetscInt cellStart, PetscInt *cell)
273: {
274:   PetscInt       coneSize;

278:   switch (dim) {
279:   case 2:
280:     DMPlexGetConeSize(dm, cellStart, &coneSize);
281:     switch (coneSize) {
282:     case 3:
283:       DMPlexLocatePoint_Simplex_2D_Internal(dm, point, cellStart, cell);
284:       break;
285:     case 4:
286:       DMPlexLocatePoint_General_2D_Internal(dm, point, cellStart, cell);
287:       break;
288:     default:
289:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
290:     }
291:     break;
292:   case 3:
293:     DMPlexGetConeSize(dm, cellStart, &coneSize);
294:     switch (coneSize) {
295:     case 4:
296:       DMPlexLocatePoint_Simplex_3D_Internal(dm, point, cellStart, cell);
297:       break;
298:     case 6:
299:       DMPlexLocatePoint_General_3D_Internal(dm, point, cellStart, cell);
300:       break;
301:     default:
302:       SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for cell with cone size %D", coneSize);
303:     }
304:     break;
305:   default:
306:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No point location for mesh dimension %D", dim);
307:   }
308:   return(0);
309: }

313: PetscErrorCode DMPlexComputeGridHash_Internal(DM dm, PetscGridHash *localBox)
314: {
315:   MPI_Comm           comm;
316:   PetscGridHash      lbox;
317:   Vec                coordinates;
318:   PetscSection       coordSection;
319:   Vec                coordsLocal;
320:   const PetscScalar *coords;
321:   PetscInt          *dboxes, *boxes;
322:   PetscInt           n[3] = {10, 10, 10};
323:   PetscInt           dim, N, cStart, cEnd, cMax, c, i;
324:   PetscErrorCode     ierr;

327:   PetscObjectGetComm((PetscObject) dm, &comm);
328:   DMGetCoordinatesLocal(dm, &coordinates);
329:   DMGetCoordinateDim(dm, &dim);
330:   VecGetLocalSize(coordinates, &N);
331:   VecGetArrayRead(coordinates, &coords);
332:   PetscGridHashCreate(comm, dim, coords, &lbox);
333:   for (i = 0; i < N; i += dim) {PetscGridHashEnlarge(lbox, &coords[i]);}
334:   VecRestoreArrayRead(coordinates, &coords);
335:   PetscGridHashSetGrid(lbox, n, NULL);
336: #if 0
337:   /* Could define a custom reduction to merge these */
338:   MPIU_Allreduce(lbox->lower, gbox->lower, 3, MPIU_REAL, MPI_MIN, comm);
339:   MPIU_Allreduce(lbox->upper, gbox->upper, 3, MPIU_REAL, MPI_MAX, comm);
340: #endif
341:   /* Is there a reason to snap the local bounding box to a division of the global box? */
342:   /* Should we compute all overlaps of local boxes? We could do this with a rendevouz scheme partitioning the global box */
343:   /* Create label */
344:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
345:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
346:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
347:   DMLabelCreate("cells", &lbox->cellsSparse);
348:   DMLabelCreateIndex(lbox->cellsSparse, cStart, cEnd);
349:   /* Compute boxes which overlap each cell: http://stackoverflow.com/questions/13790208/triangle-square-intersection-test-in-2d */
350:   DMGetCoordinatesLocal(dm, &coordsLocal);
351:   DMGetCoordinateSection(dm, &coordSection);
352:   PetscCalloc2(16 * dim, &dboxes, 16, &boxes);
353:   for (c = cStart; c < cEnd; ++c) {
354:     const PetscReal *h       = lbox->h;
355:     PetscScalar     *ccoords = NULL;
356:     PetscInt         csize   = 0;
357:     PetscScalar      point[3];
358:     PetscInt         dlim[6], d, e, i, j, k;

360:     /* Find boxes enclosing each vertex */
361:     DMPlexVecGetClosure(dm, coordSection, coordsLocal, c, &csize, &ccoords);
362:     PetscGridHashGetEnclosingBox(lbox, csize/dim, ccoords, dboxes, boxes);
363:     /* Mark cells containing the vertices */
364:     for (e = 0; e < csize/dim; ++e) {DMLabelSetValue(lbox->cellsSparse, c, boxes[e]);}
365:     /* Get grid of boxes containing these */
366:     for (d = 0;   d < dim; ++d) {dlim[d*2+0] = dlim[d*2+1] = dboxes[d];}
367:     for (d = dim; d < 3;   ++d) {dlim[d*2+0] = dlim[d*2+1] = 0;}
368:     for (e = 1; e < dim+1; ++e) {
369:       for (d = 0; d < dim; ++d) {
370:         dlim[d*2+0] = PetscMin(dlim[d*2+0], dboxes[e*dim+d]);
371:         dlim[d*2+1] = PetscMax(dlim[d*2+1], dboxes[e*dim+d]);
372:       }
373:     }
374:     /* Check for intersection of box with cell */
375:     for (k = dlim[2*2+0], point[2] = lbox->lower[2] + k*h[2]; k <= dlim[2*2+1]; ++k, point[2] += h[2]) {
376:       for (j = dlim[1*2+0], point[1] = lbox->lower[1] + j*h[1]; j <= dlim[1*2+1]; ++j, point[1] += h[1]) {
377:         for (i = dlim[0*2+0], point[0] = lbox->lower[0] + i*h[0]; i <= dlim[0*2+1]; ++i, point[0] += h[0]) {
378:           const PetscInt box = (k*lbox->n[1] + j)*lbox->n[0] + i;
379:           PetscScalar    cpoint[3];
380:           PetscInt       cell, edge, ii, jj, kk;

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

387:                 DMPlexLocatePoint_Internal(dm, dim, cpoint, c, &cell);
388:                 if (cell >= 0) {DMLabelSetValue(lbox->cellsSparse, c, box); ii = jj = kk = 2;}
389:               }
390:             }
391:           }
392:           /* Check whether cell edge intersects any edge of these subboxes TODO vectorize this */
393:           for (edge = 0; edge < dim+1; ++edge) {
394:             PetscReal segA[6], segB[6];

396:             for (d = 0; d < dim; ++d) {segA[d] = PetscRealPart(ccoords[edge*dim+d]); segA[dim+d] = PetscRealPart(ccoords[((edge+1)%(dim+1))*dim+d]);}
397:             for (kk = 0; kk < (dim > 2 ? 2 : 1); ++kk) {
398:               if (dim > 2) {segB[2]     = PetscRealPart(point[2]);
399:                             segB[dim+2] = PetscRealPart(point[2]) + kk*h[2];}
400:               for (jj = 0; jj < (dim > 1 ? 2 : 1); ++jj) {
401:                 if (dim > 1) {segB[1]     = PetscRealPart(point[1]);
402:                               segB[dim+1] = PetscRealPart(point[1]) + jj*h[1];}
403:                 for (ii = 0; ii < 2; ++ii) {
404:                   PetscBool intersects;

406:                   segB[0]     = PetscRealPart(point[0]);
407:                   segB[dim+0] = PetscRealPart(point[0]) + ii*h[0];
408:                   DMPlexGetLineIntersection_2D_Internal(segA, segB, NULL, &intersects);
409:                   if (intersects) {DMLabelSetValue(lbox->cellsSparse, c, box); edge = ii = jj = kk = dim+1;}
410:                 }
411:               }
412:             }
413:           }
414:         }
415:       }
416:     }
417:     DMPlexVecRestoreClosure(dm, coordSection, coordsLocal, c, NULL, &ccoords);
418:   }
419:   PetscFree2(dboxes, boxes);
420:   DMLabelConvertToSection(lbox->cellsSparse, &lbox->cellSection, &lbox->cells);
421:   DMLabelDestroy(&lbox->cellsSparse);
422:   *localBox = lbox;
423:   return(0);
424: }

428: PetscErrorCode DMLocatePoints_Plex(DM dm, Vec v, PetscSF cellSF)
429: {
430:   DM_Plex        *mesh = (DM_Plex *) dm->data;
431:   PetscBool       hash = mesh->useHashLocation;
432:   PetscInt        bs, numPoints, p, numFound, *found = NULL;
433:   PetscInt        dim, cStart, cEnd, cMax, numCells, c;
434:   const PetscInt *boxCells;
435:   PetscSFNode    *cells;
436:   PetscScalar    *a;
437:   PetscMPIInt     result;
438:   PetscErrorCode  ierr;

441:   DMGetCoordinateDim(dm, &dim);
442:   VecGetBlockSize(v, &bs);
443:   MPI_Comm_compare(PetscObjectComm((PetscObject)cellSF),PETSC_COMM_SELF,&result);
444:   if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)cellSF),PETSC_ERR_SUP, "Trying parallel point location: only local point location supported");
445:   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);
446:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
447:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
448:   if (cMax >= 0) cEnd = PetscMin(cEnd, cMax);
449:   VecGetLocalSize(v, &numPoints);
450:   VecGetArray(v, &a);
451:   numPoints /= bs;
452:   PetscMalloc1(numPoints, &cells);
453:   if (hash) {
454:     if (!mesh->lbox) {PetscInfo(dm, "Initializing grid hashing");DMPlexComputeGridHash_Internal(dm, &mesh->lbox);}
455:     /* Designate the local box for each point */
456:     /* Send points to correct process */
457:     /* Search cells that lie in each subbox */
458:     /*   Should we bin points before doing search? */
459:     ISGetIndices(mesh->lbox->cells, &boxCells);
460:   }
461:   for (p = 0, numFound = 0; p < numPoints; ++p) {
462:     const PetscScalar *point = &a[p*bs];
463:     PetscInt           dbin[3], bin, cell = -1, cellOffset;

465:     cells[p].rank  = -1;
466:     cells[p].index = -1;
467:     if (hash) {
468:       PetscGridHashGetEnclosingBox(mesh->lbox, 1, point, dbin, &bin);
469:       /* TODO Lay an interface over this so we can switch between Section (dense) and Label (sparse) */
470:       PetscSectionGetDof(mesh->lbox->cellSection, bin, &numCells);
471:       PetscSectionGetOffset(mesh->lbox->cellSection, bin, &cellOffset);
472:       for (c = cellOffset; c < cellOffset + numCells; ++c) {
473:         DMPlexLocatePoint_Internal(dm, dim, point, boxCells[c], &cell);
474:         if (cell >= 0) {
475:           cells[p].rank = 0;
476:           cells[p].index = cell;
477:           numFound++;
478:           break;
479:         }
480:       }
481:     } else {
482:       for (c = cStart; c < cEnd; ++c) {
483:         DMPlexLocatePoint_Internal(dm, dim, point, c, &cell);
484:         if (cell >= 0) {
485:           cells[p].rank = 0;
486:           cells[p].index = cell;
487:           numFound++;
488:           break;
489:         }
490:       }
491:     }
492:   }
493:   if (hash) {ISRestoreIndices(mesh->lbox->cells, &boxCells);}
494:   /* Check for highest numbered proc that claims a point (do we care?) */
495:   VecRestoreArray(v, &a);
496:   if (numFound < numPoints) {
497:     PetscMalloc1(numFound,&found);
498:     for (p = 0, numFound = 0; p < numPoints; p++) {
499:       if (cells[p].rank >= 0 && cells[p].index >= 0) {
500:         if (numFound < p) {
501:           cells[numFound] = cells[p];
502:         }
503:         found[numFound++] = p;
504:       }
505:     }
506:   }
507:   PetscSFSetGraph(cellSF, cEnd - cStart, numFound, found, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);
508:   return(0);
509: }

513: /*
514:   DMPlexComputeProjection2Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 2D
515: */
516: PetscErrorCode DMPlexComputeProjection2Dto1D_Internal(PetscScalar coords[], PetscReal R[])
517: {
518:   const PetscReal x = PetscRealPart(coords[2] - coords[0]);
519:   const PetscReal y = PetscRealPart(coords[3] - coords[1]);
520:   const PetscReal r = PetscSqrtReal(x*x + y*y), c = x/r, s = y/r;

523:   R[0] = c; R[1] = -s;
524:   R[2] = s; R[3] =  c;
525:   coords[0] = 0.0;
526:   coords[1] = r;
527:   return(0);
528: }

532: /*
533:   DMPlexComputeProjection3Dto1D_Internal - Rewrite coordinates to be the 1D projection of the 3D

535:   This uses the basis completion described by Frisvad,

537:   http://www.imm.dtu.dk/~jerf/papers/abstracts/onb.html
538:   DOI:10.1080/2165347X.2012.689606
539: */
540: PetscErrorCode DMPlexComputeProjection3Dto1D_Internal(PetscScalar coords[], PetscReal R[])
541: {
542:   PetscReal      x    = PetscRealPart(coords[3] - coords[0]);
543:   PetscReal      y    = PetscRealPart(coords[4] - coords[1]);
544:   PetscReal      z    = PetscRealPart(coords[5] - coords[2]);
545:   PetscReal      r    = PetscSqrtReal(x*x + y*y + z*z);
546:   PetscReal      rinv = 1. / r;

549:   x *= rinv; y *= rinv; z *= rinv;
550:   if (x > 0.) {
551:     PetscReal inv1pX   = 1./ (1. + x);

553:     R[0] = x; R[1] = -y;              R[2] = -z;
554:     R[3] = y; R[4] = 1. - y*y*inv1pX; R[5] =     -y*z*inv1pX;
555:     R[6] = z; R[7] =     -y*z*inv1pX; R[8] = 1. - z*z*inv1pX;
556:   }
557:   else {
558:     PetscReal inv1mX   = 1./ (1. - x);

560:     R[0] = x; R[1] = z;               R[2] = y;
561:     R[3] = y; R[4] =     -y*z*inv1mX; R[5] = 1. - y*y*inv1mX;
562:     R[6] = z; R[7] = 1. - z*z*inv1mX; R[8] =     -y*z*inv1mX;
563:   }
564:   coords[0] = 0.0;
565:   coords[1] = r;
566:   return(0);
567: }

571: /*
572:   DMPlexComputeProjection3Dto2D_Internal - Rewrite coordinates to be the 2D projection of the 3D
573: */
574: PetscErrorCode DMPlexComputeProjection3Dto2D_Internal(PetscInt coordSize, PetscScalar coords[], PetscReal R[])
575: {
576:   PetscReal      x1[3],  x2[3], n[3], norm;
577:   PetscReal      x1p[3], x2p[3], xnp[3];
578:   PetscReal      sqrtz, alpha;
579:   const PetscInt dim = 3;
580:   PetscInt       d, e, p;

583:   /* 0) Calculate normal vector */
584:   for (d = 0; d < dim; ++d) {
585:     x1[d] = PetscRealPart(coords[1*dim+d] - coords[0*dim+d]);
586:     x2[d] = PetscRealPart(coords[2*dim+d] - coords[0*dim+d]);
587:   }
588:   n[0] = x1[1]*x2[2] - x1[2]*x2[1];
589:   n[1] = x1[2]*x2[0] - x1[0]*x2[2];
590:   n[2] = x1[0]*x2[1] - x1[1]*x2[0];
591:   norm = PetscSqrtReal(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
592:   n[0] /= norm;
593:   n[1] /= norm;
594:   n[2] /= norm;
595:   /* 1) Take the normal vector and rotate until it is \hat z

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

599:     R = /  alpha nx nz  alpha ny nz -1/alpha \
600:         | -alpha ny     alpha nx        0    |
601:         \     nx            ny         nz    /

603:     will rotate the normal vector to \hat z
604:   */
605:   sqrtz = PetscSqrtReal(1.0 - n[2]*n[2]);
606:   /* Check for n = z */
607:   if (sqrtz < 1.0e-10) {
608:     const PetscInt s = PetscSign(n[2]);
609:     /* If nz < 0, rotate 180 degrees around x-axis */
610:     for (p = 3; p < coordSize/3; ++p) {
611:       coords[p*2+0] = PetscRealPart(coords[p*dim+0] - coords[0*dim+0]);
612:       coords[p*2+1] = (PetscRealPart(coords[p*dim+1] - coords[0*dim+1])) * s;
613:     }
614:     coords[0] = 0.0;
615:     coords[1] = 0.0;
616:     coords[2] = x1[0];
617:     coords[3] = x1[1] * s;
618:     coords[4] = x2[0];
619:     coords[5] = x2[1] * s;
620:     R[0] = 1.0;     R[1] = 0.0;     R[2] = 0.0;
621:     R[3] = 0.0;     R[4] = 1.0 * s; R[5] = 0.0;
622:     R[6] = 0.0;     R[7] = 0.0;     R[8] = 1.0 * s;
623:     return(0);
624:   }
625:   alpha = 1.0/sqrtz;
626:   R[0] =  alpha*n[0]*n[2]; R[1] = alpha*n[1]*n[2]; R[2] = -sqrtz;
627:   R[3] = -alpha*n[1];      R[4] = alpha*n[0];      R[5] = 0.0;
628:   R[6] =  n[0];            R[7] = n[1];            R[8] = n[2];
629:   for (d = 0; d < dim; ++d) {
630:     x1p[d] = 0.0;
631:     x2p[d] = 0.0;
632:     for (e = 0; e < dim; ++e) {
633:       x1p[d] += R[d*dim+e]*x1[e];
634:       x2p[d] += R[d*dim+e]*x2[e];
635:     }
636:   }
637:   if (PetscAbsReal(x1p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
638:   if (PetscAbsReal(x2p[2]) > 1.0e-9) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid rotation calculated");
639:   /* 2) Project to (x, y) */
640:   for (p = 3; p < coordSize/3; ++p) {
641:     for (d = 0; d < dim; ++d) {
642:       xnp[d] = 0.0;
643:       for (e = 0; e < dim; ++e) {
644:         xnp[d] += R[d*dim+e]*PetscRealPart(coords[p*dim+e] - coords[0*dim+e]);
645:       }
646:       if (d < dim-1) coords[p*2+d] = xnp[d];
647:     }
648:   }
649:   coords[0] = 0.0;
650:   coords[1] = 0.0;
651:   coords[2] = x1p[0];
652:   coords[3] = x1p[1];
653:   coords[4] = x2p[0];
654:   coords[5] = x2p[1];
655:   /* Output R^T which rotates \hat z to the input normal */
656:   for (d = 0; d < dim; ++d) {
657:     for (e = d+1; e < dim; ++e) {
658:       PetscReal tmp;

660:       tmp        = R[d*dim+e];
661:       R[d*dim+e] = R[e*dim+d];
662:       R[e*dim+d] = tmp;
663:     }
664:   }
665:   return(0);
666: }

670: PETSC_UNUSED
671: PETSC_STATIC_INLINE void Volume_Triangle_Internal(PetscReal *vol, PetscReal coords[])
672: {
673:   /* Signed volume is 1/2 the determinant

675:    |  1  1  1 |
676:    | x0 x1 x2 |
677:    | y0 y1 y2 |

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

681:    | x1 x2 |
682:    | y1 y2 |
683:   */
684:   const PetscReal x1 = coords[2] - coords[0], y1 = coords[3] - coords[1];
685:   const PetscReal x2 = coords[4] - coords[0], y2 = coords[5] - coords[1];
686:   PetscReal       M[4], detM;
687:   M[0] = x1; M[1] = x2;
688:   M[2] = y1; M[3] = y2;
689:   DMPlex_Det2D_Internal(&detM, M);
690:   *vol = 0.5*detM;
691:   (void)PetscLogFlops(5.0);
692: }

696: PETSC_STATIC_INLINE void Volume_Triangle_Origin_Internal(PetscReal *vol, PetscReal coords[])
697: {
698:   DMPlex_Det2D_Internal(vol, coords);
699:   *vol *= 0.5;
700: }

704: PETSC_UNUSED
705: PETSC_STATIC_INLINE void Volume_Tetrahedron_Internal(PetscReal *vol, PetscReal coords[])
706: {
707:   /* Signed volume is 1/6th of the determinant

709:    |  1  1  1  1 |
710:    | x0 x1 x2 x3 |
711:    | y0 y1 y2 y3 |
712:    | z0 z1 z2 z3 |

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

716:    | x1 x2 x3 |
717:    | y1 y2 y3 |
718:    | z1 z2 z3 |
719:   */
720:   const PetscReal x1 = coords[3] - coords[0], y1 = coords[4]  - coords[1], z1 = coords[5]  - coords[2];
721:   const PetscReal x2 = coords[6] - coords[0], y2 = coords[7]  - coords[1], z2 = coords[8]  - coords[2];
722:   const PetscReal x3 = coords[9] - coords[0], y3 = coords[10] - coords[1], z3 = coords[11] - coords[2];
723:   PetscReal       M[9], detM;
724:   M[0] = x1; M[1] = x2; M[2] = x3;
725:   M[3] = y1; M[4] = y2; M[5] = y3;
726:   M[6] = z1; M[7] = z2; M[8] = z3;
727:   DMPlex_Det3D_Internal(&detM, M);
728:   *vol = -0.16666666666666666666666*detM;
729:   (void)PetscLogFlops(10.0);
730: }

734: PETSC_STATIC_INLINE void Volume_Tetrahedron_Origin_Internal(PetscReal *vol, PetscReal coords[])
735: {
736:   DMPlex_Det3D_Internal(vol, coords);
737:   *vol *= -0.16666666666666666666666;
738: }

742: static PetscErrorCode DMPlexComputeLineGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
743: {
744:   PetscSection   coordSection;
745:   Vec            coordinates;
746:   PetscScalar   *coords = NULL;
747:   PetscInt       numCoords, d;

751:   DMGetCoordinatesLocal(dm, &coordinates);
752:   DMGetCoordinateSection(dm, &coordSection);
753:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
754:   if (invJ && !J) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "In order to compute invJ, J must not be NULL");
755:   *detJ = 0.0;
756:   if (numCoords == 6) {
757:     const PetscInt dim = 3;
758:     PetscReal      R[9], J0;

760:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
761:     DMPlexComputeProjection3Dto1D_Internal(coords, R);
762:     if (J)    {
763:       J0   = 0.5*PetscRealPart(coords[1]);
764:       J[0] = R[0]*J0; J[1] = R[1]; J[2] = R[2];
765:       J[3] = R[3]*J0; J[4] = R[4]; J[5] = R[5];
766:       J[6] = R[6]*J0; J[7] = R[7]; J[8] = R[8];
767:       DMPlex_Det3D_Internal(detJ, J);
768:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
769:     }
770:   } else if (numCoords == 4) {
771:     const PetscInt dim = 2;
772:     PetscReal      R[4], J0;

774:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
775:     DMPlexComputeProjection2Dto1D_Internal(coords, R);
776:     if (J)    {
777:       J0   = 0.5*PetscRealPart(coords[1]);
778:       J[0] = R[0]*J0; J[1] = R[1];
779:       J[2] = R[2]*J0; J[3] = R[3];
780:       DMPlex_Det2D_Internal(detJ, J);
781:       if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
782:     }
783:   } else if (numCoords == 2) {
784:     const PetscInt dim = 1;

786:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
787:     if (J)    {
788:       J[0]  = 0.5*(PetscRealPart(coords[1]) - PetscRealPart(coords[0]));
789:       *detJ = J[0];
790:       PetscLogFlops(2.0);
791:       if (invJ) {invJ[0] = 1.0/J[0]; PetscLogFlops(1.0);}
792:     }
793:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this segment is %D != 2", numCoords);
794:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
795:   return(0);
796: }

800: static PetscErrorCode DMPlexComputeTriangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
801: {
802:   PetscSection   coordSection;
803:   Vec            coordinates;
804:   PetscScalar   *coords = NULL;
805:   PetscInt       numCoords, d, f, g;

809:   DMGetCoordinatesLocal(dm, &coordinates);
810:   DMGetCoordinateSection(dm, &coordSection);
811:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
812:   *detJ = 0.0;
813:   if (numCoords == 9) {
814:     const PetscInt dim = 3;
815:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

817:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
818:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
819:     if (J)    {
820:       const PetscInt pdim = 2;

822:       for (d = 0; d < pdim; d++) {
823:         for (f = 0; f < pdim; f++) {
824:           J0[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
825:         }
826:       }
827:       PetscLogFlops(8.0);
828:       DMPlex_Det3D_Internal(detJ, J0);
829:       for (d = 0; d < dim; d++) {
830:         for (f = 0; f < dim; f++) {
831:           J[d*dim+f] = 0.0;
832:           for (g = 0; g < dim; g++) {
833:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
834:           }
835:         }
836:       }
837:       PetscLogFlops(18.0);
838:     }
839:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
840:   } else if (numCoords == 6) {
841:     const PetscInt dim = 2;

843:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
844:     if (J)    {
845:       for (d = 0; d < dim; d++) {
846:         for (f = 0; f < dim; f++) {
847:           J[d*dim+f] = 0.5*(PetscRealPart(coords[(f+1)*dim+d]) - PetscRealPart(coords[0*dim+d]));
848:         }
849:       }
850:       PetscLogFlops(8.0);
851:       DMPlex_Det2D_Internal(detJ, J);
852:     }
853:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
854:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this triangle is %D != 6 or 9", numCoords);
855:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
856:   return(0);
857: }

861: static PetscErrorCode DMPlexComputeRectangleGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
862: {
863:   PetscSection   coordSection;
864:   Vec            coordinates;
865:   PetscScalar   *coords = NULL;
866:   PetscInt       numCoords, d, f, g;

870:   DMGetCoordinatesLocal(dm, &coordinates);
871:   DMGetCoordinateSection(dm, &coordSection);
872:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
873:   *detJ = 0.0;
874:   if (numCoords == 12) {
875:     const PetscInt dim = 3;
876:     PetscReal      R[9], J0[9] = {1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0};

878:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
879:     DMPlexComputeProjection3Dto2D_Internal(numCoords, coords, R);
880:     if (J)    {
881:       const PetscInt pdim = 2;

883:       for (d = 0; d < pdim; d++) {
884:         J0[d*dim+0] = 0.5*(PetscRealPart(coords[1*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
885:         J0[d*dim+1] = 0.5*(PetscRealPart(coords[3*pdim+d]) - PetscRealPart(coords[0*pdim+d]));
886:       }
887:       PetscLogFlops(8.0);
888:       DMPlex_Det3D_Internal(detJ, J0);
889:       for (d = 0; d < dim; d++) {
890:         for (f = 0; f < dim; f++) {
891:           J[d*dim+f] = 0.0;
892:           for (g = 0; g < dim; g++) {
893:             J[d*dim+f] += R[d*dim+g]*J0[g*dim+f];
894:           }
895:         }
896:       }
897:       PetscLogFlops(18.0);
898:     }
899:     if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
900:   } else if ((numCoords == 8) || (numCoords == 16)) {
901:     const PetscInt dim = 2;

903:     if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
904:     if (J)    {
905:       for (d = 0; d < dim; d++) {
906:         J[d*dim+0] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
907:         J[d*dim+1] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
908:       }
909:       PetscLogFlops(8.0);
910:       DMPlex_Det2D_Internal(detJ, J);
911:     }
912:     if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
913:   } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The number of coordinates for this quadrilateral is %D != 8 or 12", numCoords);
914:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, &numCoords, &coords);
915:   return(0);
916: }

920: static PetscErrorCode DMPlexComputeTetrahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
921: {
922:   PetscSection   coordSection;
923:   Vec            coordinates;
924:   PetscScalar   *coords = NULL;
925:   const PetscInt dim = 3;
926:   PetscInt       d;

930:   DMGetCoordinatesLocal(dm, &coordinates);
931:   DMGetCoordinateSection(dm, &coordSection);
932:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
933:   *detJ = 0.0;
934:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
935:   if (J)    {
936:     for (d = 0; d < dim; d++) {
937:       /* I orient with outward face normals */
938:       J[d*dim+0] = 0.5*(PetscRealPart(coords[2*dim+d]) - PetscRealPart(coords[0*dim+d]));
939:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
940:       J[d*dim+2] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
941:     }
942:     PetscLogFlops(18.0);
943:     DMPlex_Det3D_Internal(detJ, J);
944:   }
945:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
946:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
947:   return(0);
948: }

952: static PetscErrorCode DMPlexComputeHexahedronGeometry_Internal(DM dm, PetscInt e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
953: {
954:   PetscSection   coordSection;
955:   Vec            coordinates;
956:   PetscScalar   *coords = NULL;
957:   const PetscInt dim = 3;
958:   PetscInt       d;

962:   DMGetCoordinatesLocal(dm, &coordinates);
963:   DMGetCoordinateSection(dm, &coordSection);
964:   DMPlexVecGetClosure(dm, coordSection, coordinates, e, NULL, &coords);
965:   *detJ = 0.0;
966:   if (v0)   {for (d = 0; d < dim; d++) v0[d] = PetscRealPart(coords[d]);}
967:   if (J)    {
968:     for (d = 0; d < dim; d++) {
969:       J[d*dim+0] = 0.5*(PetscRealPart(coords[3*dim+d]) - PetscRealPart(coords[0*dim+d]));
970:       J[d*dim+1] = 0.5*(PetscRealPart(coords[1*dim+d]) - PetscRealPart(coords[0*dim+d]));
971:       J[d*dim+2] = 0.5*(PetscRealPart(coords[4*dim+d]) - PetscRealPart(coords[0*dim+d]));
972:     }
973:     PetscLogFlops(18.0);
974:     DMPlex_Det3D_Internal(detJ, J);
975:   }
976:   if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
977:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, e, NULL, &coords);
978:   return(0);
979: }

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

986:   Collective on DM

988:   Input Arguments:
989: + dm   - the DM
990: - cell - the cell

992:   Output Arguments:
993: + v0   - the translation part of this affine transform
994: . J    - the Jacobian of the transform from the reference element
995: . invJ - the inverse of the Jacobian
996: - detJ - the Jacobian determinant

998:   Level: advanced

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

1004: .seealso: DMPlexComputeCellGeometryFEM(), DMGetCoordinateSection(), DMGetCoordinateVec()
1005: @*/
1006: PetscErrorCode DMPlexComputeCellGeometryAffineFEM(DM dm, PetscInt cell, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1007: {
1008:   PetscInt       depth, dim, coneSize;

1012:   DMPlexGetDepth(dm, &depth);
1013:   DMPlexGetConeSize(dm, cell, &coneSize);
1014:   if (depth == 1) {
1015:     DMGetDimension(dm, &dim);
1016:   } else {
1017:     DMLabel depth;

1019:     DMPlexGetDepthLabel(dm, &depth);
1020:     DMLabelGetValue(depth, cell, &dim);
1021:   }
1022:   switch (dim) {
1023:   case 1:
1024:     DMPlexComputeLineGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1025:     break;
1026:   case 2:
1027:     switch (coneSize) {
1028:     case 3:
1029:       DMPlexComputeTriangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1030:       break;
1031:     case 4:
1032:       DMPlexComputeRectangleGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1033:       break;
1034:     default:
1035:       SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1036:     }
1037:     break;
1038:   case 3:
1039:     switch (coneSize) {
1040:     case 4:
1041:       DMPlexComputeTetrahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1042:       break;
1043:     case 6: /* Faces */
1044:     case 8: /* Vertices */
1045:       DMPlexComputeHexahedronGeometry_Internal(dm, cell, v0, J, invJ, detJ);
1046:       break;
1047:     default:
1048:         SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported number of faces %D in cell %D for element geometry computation", coneSize, cell);
1049:     }
1050:       break;
1051:   default:
1052:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1053:   }
1054:   return(0);
1055: }

1059: static PetscErrorCode DMPlexComputeIsoparametricGeometry_Internal(DM dm, PetscFE fe, PetscInt point, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1060: {
1061:   PetscQuadrature  quad;
1062:   PetscSection     coordSection;
1063:   Vec              coordinates;
1064:   PetscScalar     *coords = NULL;
1065:   const PetscReal *quadPoints;
1066:   PetscReal       *basisDer;
1067:   PetscInt         dim, cdim, pdim, qdim, Nq, numCoords, d, q;
1068:   PetscErrorCode   ierr;

1071:   DMGetCoordinatesLocal(dm, &coordinates);
1072:   DMGetCoordinateSection(dm, &coordSection);
1073:   DMPlexVecGetClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1074:   DMGetDimension(dm, &dim);
1075:   DMGetCoordinateDim(dm, &cdim);
1076:   PetscFEGetQuadrature(fe, &quad);
1077:   PetscFEGetDimension(fe, &pdim);
1078:   PetscQuadratureGetData(quad, &qdim, &Nq, &quadPoints, NULL);
1079:   PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
1080:   *detJ = 0.0;
1081:   if (qdim != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Point dimension %d != quadrature dimension %d", dim, qdim);
1082:   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);
1083:   if (v0) {for (d = 0; d < cdim; d++) v0[d] = PetscRealPart(coords[d]);}
1084:   if (J) {
1085:     for (q = 0; q < Nq; ++q) {
1086:       PetscInt i, j, k, c, r;

1088:       /* J = dx_i/d\xi_j = sum[k=0,n-1] dN_k/d\xi_j * x_i(k) */
1089:       for (k = 0; k < pdim; ++k)
1090:         for (j = 0; j < dim; ++j)
1091:           for (i = 0; i < cdim; ++i)
1092:             J[(q*cdim + i)*dim + j] += basisDer[(q*pdim + k)*dim + j] * PetscRealPart(coords[k*cdim + i]);
1093:       PetscLogFlops(2.0*pdim*dim*cdim);
1094:       if (cdim > dim) {
1095:         for (c = dim; c < cdim; ++c)
1096:           for (r = 0; r < cdim; ++r)
1097:             J[r*cdim+c] = r == c ? 1.0 : 0.0;
1098:       }
1099:       switch (cdim) {
1100:       case 3:
1101:         DMPlex_Det3D_Internal(detJ, J);
1102:         if (invJ) {DMPlex_Invert3D_Internal(invJ, J, *detJ);}
1103:         break;
1104:       case 2:
1105:         DMPlex_Det2D_Internal(detJ, J);
1106:         if (invJ) {DMPlex_Invert2D_Internal(invJ, J, *detJ);}
1107:         break;
1108:       case 1:
1109:         *detJ = J[0];
1110:         if (invJ) invJ[0] = 1.0/J[0];
1111:       }
1112:     }
1113:   }
1114:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, point, &numCoords, &coords);
1115:   return(0);
1116: }

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

1123:   Collective on DM

1125:   Input Arguments:
1126: + dm   - the DM
1127: . cell - the cell
1128: - fe   - the finite element containing the quadrature

1130:   Output Arguments:
1131: + v0   - the translation part of this transform
1132: . J    - the Jacobian of the transform from the reference element at each quadrature point
1133: . invJ - the inverse of the Jacobian at each quadrature point
1134: - detJ - the Jacobian determinant at each quadrature point

1136:   Level: advanced

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

1142: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
1143: @*/
1144: PetscErrorCode DMPlexComputeCellGeometryFEM(DM dm, PetscInt cell, PetscFE fe, PetscReal *v0, PetscReal *J, PetscReal *invJ, PetscReal *detJ)
1145: {

1149:   if (!fe) {DMPlexComputeCellGeometryAffineFEM(dm, cell, v0, J, invJ, detJ);}
1150:   else     {DMPlexComputeIsoparametricGeometry_Internal(dm, fe, cell, v0, J, invJ, detJ);}
1151:   return(0);
1152: }

1156: static PetscErrorCode DMPlexComputeGeometryFVM_1D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1157: {
1158:   PetscSection   coordSection;
1159:   Vec            coordinates;
1160:   PetscScalar   *coords = NULL;
1161:   PetscScalar    tmp[2];
1162:   PetscInt       coordSize;

1166:   DMGetCoordinatesLocal(dm, &coordinates);
1167:   DMGetCoordinateSection(dm, &coordSection);
1168:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1169:   if (dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support 2D edges right now");
1170:   DMLocalizeCoordinate_Internal(dm, dim, coords, &coords[dim], tmp);
1171:   if (centroid) {
1172:     centroid[0] = 0.5*PetscRealPart(coords[0] + tmp[0]);
1173:     centroid[1] = 0.5*PetscRealPart(coords[1] + tmp[1]);
1174:   }
1175:   if (normal) {
1176:     PetscReal norm;

1178:     normal[0]  = -PetscRealPart(coords[1] - tmp[1]);
1179:     normal[1]  =  PetscRealPart(coords[0] - tmp[0]);
1180:     norm       = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1]);
1181:     normal[0] /= norm;
1182:     normal[1] /= norm;
1183:   }
1184:   if (vol) {
1185:     *vol = PetscSqrtReal(PetscSqr(PetscRealPart(coords[0] - tmp[0])) + PetscSqr(PetscRealPart(coords[1] - tmp[1])));
1186:   }
1187:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1188:   return(0);
1189: }

1193: /* Centroid_i = (\sum_n A_n Cn_i ) / A */
1194: static PetscErrorCode DMPlexComputeGeometryFVM_2D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1195: {
1196:   PetscSection   coordSection;
1197:   Vec            coordinates;
1198:   PetscScalar   *coords = NULL;
1199:   PetscReal      vsum = 0.0, csum[3] = {0.0, 0.0, 0.0}, vtmp, ctmp[4], v0[3], R[9];
1200:   PetscInt       tdim = 2, coordSize, numCorners, p, d, e;

1204:   DMGetCoordinatesLocal(dm, &coordinates);
1205:   DMPlexGetConeSize(dm, cell, &numCorners);
1206:   DMGetCoordinateSection(dm, &coordSection);
1207:   DMPlexVecGetClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1208:   DMGetCoordinateDim(dm, &dim);
1209:   if (dim > 2 && centroid) {
1210:     v0[0] = PetscRealPart(coords[0]);
1211:     v0[1] = PetscRealPart(coords[1]);
1212:     v0[2] = PetscRealPart(coords[2]);
1213:   }
1214:   if (normal) {
1215:     if (dim > 2) {
1216:       const PetscReal x0 = PetscRealPart(coords[dim+0] - coords[0]), x1 = PetscRealPart(coords[dim*2+0] - coords[0]);
1217:       const PetscReal y0 = PetscRealPart(coords[dim+1] - coords[1]), y1 = PetscRealPart(coords[dim*2+1] - coords[1]);
1218:       const PetscReal z0 = PetscRealPart(coords[dim+2] - coords[2]), z1 = PetscRealPart(coords[dim*2+2] - coords[2]);
1219:       PetscReal       norm;

1221:       normal[0] = y0*z1 - z0*y1;
1222:       normal[1] = z0*x1 - x0*z1;
1223:       normal[2] = x0*y1 - y0*x1;
1224:       norm = PetscSqrtReal(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
1225:       normal[0] /= norm;
1226:       normal[1] /= norm;
1227:       normal[2] /= norm;
1228:     } else {
1229:       for (d = 0; d < dim; ++d) normal[d] = 0.0;
1230:     }
1231:   }
1232:   if (dim == 3) {DMPlexComputeProjection3Dto2D_Internal(coordSize, coords, R);}
1233:   for (p = 0; p < numCorners; ++p) {
1234:     /* Need to do this copy to get types right */
1235:     for (d = 0; d < tdim; ++d) {
1236:       ctmp[d]      = PetscRealPart(coords[p*tdim+d]);
1237:       ctmp[tdim+d] = PetscRealPart(coords[((p+1)%numCorners)*tdim+d]);
1238:     }
1239:     Volume_Triangle_Origin_Internal(&vtmp, ctmp);
1240:     vsum += vtmp;
1241:     for (d = 0; d < tdim; ++d) {
1242:       csum[d] += (ctmp[d] + ctmp[tdim+d])*vtmp;
1243:     }
1244:   }
1245:   for (d = 0; d < tdim; ++d) {
1246:     csum[d] /= (tdim+1)*vsum;
1247:   }
1248:   DMPlexVecRestoreClosure(dm, coordSection, coordinates, cell, &coordSize, &coords);
1249:   if (vol) *vol = PetscAbsReal(vsum);
1250:   if (centroid) {
1251:     if (dim > 2) {
1252:       for (d = 0; d < dim; ++d) {
1253:         centroid[d] = v0[d];
1254:         for (e = 0; e < dim; ++e) {
1255:           centroid[d] += R[d*dim+e]*csum[e];
1256:         }
1257:       }
1258:     } else for (d = 0; d < dim; ++d) centroid[d] = csum[d];
1259:   }
1260:   return(0);
1261: }

1265: /* Centroid_i = (\sum_n V_n Cn_i ) / V */
1266: static PetscErrorCode DMPlexComputeGeometryFVM_3D_Internal(DM dm, PetscInt dim, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1267: {
1268:   PetscSection    coordSection;
1269:   Vec             coordinates;
1270:   PetscScalar    *coords = NULL;
1271:   PetscReal       vsum = 0.0, vtmp, coordsTmp[3*3];
1272:   const PetscInt *faces, *facesO;
1273:   PetscInt        numFaces, f, coordSize, numCorners, p, d;
1274:   PetscErrorCode  ierr;

1277:   if (PetscUnlikely(dim > 3)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"No support for dim %D > 3",dim);
1278:   DMGetCoordinatesLocal(dm, &coordinates);
1279:   DMGetCoordinateSection(dm, &coordSection);

1281:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] = 0.0;
1282:   DMPlexGetConeSize(dm, cell, &numFaces);
1283:   DMPlexGetCone(dm, cell, &faces);
1284:   DMPlexGetConeOrientation(dm, cell, &facesO);
1285:   for (f = 0; f < numFaces; ++f) {
1286:     DMPlexVecGetClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1287:     numCorners = coordSize/dim;
1288:     switch (numCorners) {
1289:     case 3:
1290:       for (d = 0; d < dim; ++d) {
1291:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1292:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1293:         coordsTmp[2*dim+d] = PetscRealPart(coords[2*dim+d]);
1294:       }
1295:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1296:       if (facesO[f] < 0) vtmp = -vtmp;
1297:       vsum += vtmp;
1298:       if (centroid) {           /* Centroid of OABC = (a+b+c)/4 */
1299:         for (d = 0; d < dim; ++d) {
1300:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1301:         }
1302:       }
1303:       break;
1304:     case 4:
1305:       /* DO FOR PYRAMID */
1306:       /* First tet */
1307:       for (d = 0; d < dim; ++d) {
1308:         coordsTmp[0*dim+d] = PetscRealPart(coords[0*dim+d]);
1309:         coordsTmp[1*dim+d] = PetscRealPart(coords[1*dim+d]);
1310:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1311:       }
1312:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1313:       if (facesO[f] < 0) vtmp = -vtmp;
1314:       vsum += vtmp;
1315:       if (centroid) {
1316:         for (d = 0; d < dim; ++d) {
1317:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1318:         }
1319:       }
1320:       /* Second tet */
1321:       for (d = 0; d < dim; ++d) {
1322:         coordsTmp[0*dim+d] = PetscRealPart(coords[1*dim+d]);
1323:         coordsTmp[1*dim+d] = PetscRealPart(coords[2*dim+d]);
1324:         coordsTmp[2*dim+d] = PetscRealPart(coords[3*dim+d]);
1325:       }
1326:       Volume_Tetrahedron_Origin_Internal(&vtmp, coordsTmp);
1327:       if (facesO[f] < 0) vtmp = -vtmp;
1328:       vsum += vtmp;
1329:       if (centroid) {
1330:         for (d = 0; d < dim; ++d) {
1331:           for (p = 0; p < 3; ++p) centroid[d] += coordsTmp[p*dim+d]*vtmp;
1332:         }
1333:       }
1334:       break;
1335:     default:
1336:       SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle faces with %D vertices", numCorners);
1337:     }
1338:     DMPlexVecRestoreClosure(dm, coordSection, coordinates, faces[f], &coordSize, &coords);
1339:   }
1340:   if (vol)     *vol = PetscAbsReal(vsum);
1341:   if (normal)   for (d = 0; d < dim; ++d) normal[d]    = 0.0;
1342:   if (centroid) for (d = 0; d < dim; ++d) centroid[d] /= (vsum*4);
1343:   return(0);
1344: }

1348: /*@C
1349:   DMPlexComputeCellGeometryFVM - Compute the volume for a given cell

1351:   Collective on DM

1353:   Input Arguments:
1354: + dm   - the DM
1355: - cell - the cell

1357:   Output Arguments:
1358: + volume   - the cell volume
1359: . centroid - the cell centroid
1360: - normal - the cell normal, if appropriate

1362:   Level: advanced

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

1368: .seealso: DMGetCoordinateSection(), DMGetCoordinateVec()
1369: @*/
1370: PetscErrorCode DMPlexComputeCellGeometryFVM(DM dm, PetscInt cell, PetscReal *vol, PetscReal centroid[], PetscReal normal[])
1371: {
1372:   PetscInt       depth, dim;

1376:   DMPlexGetDepth(dm, &depth);
1377:   DMGetDimension(dm, &dim);
1378:   if (depth != dim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mesh must be interpolated");
1379:   /* We need to keep a pointer to the depth label */
1380:   DMGetLabelValue(dm, "depth", cell, &depth);
1381:   /* Cone size is now the number of faces */
1382:   switch (depth) {
1383:   case 1:
1384:     DMPlexComputeGeometryFVM_1D_Internal(dm, dim, cell, vol, centroid, normal);
1385:     break;
1386:   case 2:
1387:     DMPlexComputeGeometryFVM_2D_Internal(dm, dim, cell, vol, centroid, normal);
1388:     break;
1389:   case 3:
1390:     DMPlexComputeGeometryFVM_3D_Internal(dm, dim, cell, vol, centroid, normal);
1391:     break;
1392:   default:
1393:     SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unsupported dimension %D for element geometry computation", dim);
1394:   }
1395:   return(0);
1396: }

1400: /* This should also take a PetscFE argument I think */
1401: PetscErrorCode DMPlexComputeGeometryFEM(DM dm, Vec *cellgeom)
1402: {
1403:   DM             dmCell;
1404:   Vec            coordinates;
1405:   PetscSection   coordSection, sectionCell;
1406:   PetscScalar   *cgeom;
1407:   PetscInt       cStart, cEnd, cMax, c;

1411:   DMClone(dm, &dmCell);
1412:   DMGetCoordinateSection(dm, &coordSection);
1413:   DMGetCoordinatesLocal(dm, &coordinates);
1414:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1415:   DMSetCoordinatesLocal(dmCell, coordinates);
1416:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
1417:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1418:   DMPlexGetHybridBounds(dm, &cMax, NULL, NULL, NULL);
1419:   cEnd = cMax < 0 ? cEnd : cMax;
1420:   PetscSectionSetChart(sectionCell, cStart, cEnd);
1421:   /* TODO This needs to be multiplied by Nq for non-affine */
1422:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFECellGeom))/sizeof(PetscScalar)));}
1423:   PetscSectionSetUp(sectionCell);
1424:   DMSetDefaultSection(dmCell, sectionCell);
1425:   PetscSectionDestroy(&sectionCell);
1426:   DMCreateLocalVector(dmCell, cellgeom);
1427:   VecGetArray(*cellgeom, &cgeom);
1428:   for (c = cStart; c < cEnd; ++c) {
1429:     PetscFECellGeom *cg;

1431:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1432:     PetscMemzero(cg, sizeof(*cg));
1433:     DMPlexComputeCellGeometryFEM(dmCell, c, NULL, cg->v0, cg->J, cg->invJ, &cg->detJ);
1434:     if (cg->detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cg->detJ, c);
1435:   }
1436:   VecRestoreArray(*cellgeom, &cgeom);
1437:   DMDestroy(&dmCell);
1438:   return(0);
1439: }

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

1446:   Input Parameter:
1447: . dm - The DM

1449:   Output Parameters:
1450: + cellgeom - A Vec of PetscFVCellGeom data
1451: . facegeom - A Vec of PetscFVFaceGeom data

1453:   Level: developer

1455: .seealso: PetscFVFaceGeom, PetscFVCellGeom, DMPlexComputeGeometryFEM()
1456: @*/
1457: PetscErrorCode DMPlexComputeGeometryFVM(DM dm, Vec *cellgeom, Vec *facegeom)
1458: {
1459:   DM             dmFace, dmCell;
1460:   DMLabel        ghostLabel;
1461:   PetscSection   sectionFace, sectionCell;
1462:   PetscSection   coordSection;
1463:   Vec            coordinates;
1464:   PetscScalar   *fgeom, *cgeom;
1465:   PetscReal      minradius, gminradius;
1466:   PetscInt       dim, cStart, cEnd, cEndInterior, c, fStart, fEnd, f;

1470:   DMGetDimension(dm, &dim);
1471:   DMGetCoordinateSection(dm, &coordSection);
1472:   DMGetCoordinatesLocal(dm, &coordinates);
1473:   /* Make cell centroids and volumes */
1474:   DMClone(dm, &dmCell);
1475:   DMSetCoordinateSection(dmCell, PETSC_DETERMINE, coordSection);
1476:   DMSetCoordinatesLocal(dmCell, coordinates);
1477:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionCell);
1478:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1479:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1480:   PetscSectionSetChart(sectionCell, cStart, cEnd);
1481:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionCell, c, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVCellGeom))/sizeof(PetscScalar)));}
1482:   PetscSectionSetUp(sectionCell);
1483:   DMSetDefaultSection(dmCell, sectionCell);
1484:   PetscSectionDestroy(&sectionCell);
1485:   DMCreateLocalVector(dmCell, cellgeom);
1486:   if (cEndInterior < 0) {
1487:     cEndInterior = cEnd;
1488:   }
1489:   VecGetArray(*cellgeom, &cgeom);
1490:   for (c = cStart; c < cEndInterior; ++c) {
1491:     PetscFVCellGeom *cg;

1493:     DMPlexPointLocalRef(dmCell, c, cgeom, &cg);
1494:     PetscMemzero(cg, sizeof(*cg));
1495:     DMPlexComputeCellGeometryFVM(dmCell, c, &cg->volume, cg->centroid, NULL);
1496:   }
1497:   /* Compute face normals and minimum cell radius */
1498:   DMClone(dm, &dmFace);
1499:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionFace);
1500:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1501:   PetscSectionSetChart(sectionFace, fStart, fEnd);
1502:   for (f = fStart; f < fEnd; ++f) {PetscSectionSetDof(sectionFace, f, (PetscInt) PetscCeilReal(((PetscReal) sizeof(PetscFVFaceGeom))/sizeof(PetscScalar)));}
1503:   PetscSectionSetUp(sectionFace);
1504:   DMSetDefaultSection(dmFace, sectionFace);
1505:   PetscSectionDestroy(&sectionFace);
1506:   DMCreateLocalVector(dmFace, facegeom);
1507:   VecGetArray(*facegeom, &fgeom);
1508:   DMGetLabel(dm, "ghost", &ghostLabel);
1509:   minradius = PETSC_MAX_REAL;
1510:   for (f = fStart; f < fEnd; ++f) {
1511:     PetscFVFaceGeom *fg;
1512:     PetscReal        area;
1513:     PetscInt         ghost = -1, d, numChildren;

1515:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1516:     DMPlexGetTreeChildren(dm,f,&numChildren,NULL);
1517:     if (ghost >= 0 || numChildren) continue;
1518:     DMPlexPointLocalRef(dmFace, f, fgeom, &fg);
1519:     DMPlexComputeCellGeometryFVM(dm, f, &area, fg->centroid, fg->normal);
1520:     for (d = 0; d < dim; ++d) fg->normal[d] *= area;
1521:     /* Flip face orientation if necessary to match ordering in support, and Update minimum radius */
1522:     {
1523:       PetscFVCellGeom *cL, *cR;
1524:       PetscInt         ncells;
1525:       const PetscInt  *cells;
1526:       PetscReal       *lcentroid, *rcentroid;
1527:       PetscReal        l[3], r[3], v[3];

1529:       DMPlexGetSupport(dm, f, &cells);
1530:       DMPlexGetSupportSize(dm, f, &ncells);
1531:       DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cL);
1532:       lcentroid = cells[0] >= cEndInterior ? fg->centroid : cL->centroid;
1533:       if (ncells > 1) {
1534:         DMPlexPointLocalRead(dmCell, cells[1], cgeom, &cR);
1535:         rcentroid = cells[1] >= cEndInterior ? fg->centroid : cR->centroid;
1536:       }
1537:       else {
1538:         rcentroid = fg->centroid;
1539:       }
1540:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, lcentroid, l);
1541:       DMLocalizeCoordinateReal_Internal(dm, dim, fg->centroid, rcentroid, r);
1542:       DMPlex_WaxpyD_Internal(dim, -1, l, r, v);
1543:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) < 0) {
1544:         for (d = 0; d < dim; ++d) fg->normal[d] = -fg->normal[d];
1545:       }
1546:       if (DMPlex_DotRealD_Internal(dim, fg->normal, v) <= 0) {
1547:         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]);
1548:         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]);
1549:         SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Direction for face %d could not be fixed", f);
1550:       }
1551:       if (cells[0] < cEndInterior) {
1552:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cL->centroid, v);
1553:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1554:       }
1555:       if (ncells > 1 && cells[1] < cEndInterior) {
1556:         DMPlex_WaxpyD_Internal(dim, -1, fg->centroid, cR->centroid, v);
1557:         minradius = PetscMin(minradius, DMPlex_NormD_Internal(dim, v));
1558:       }
1559:     }
1560:   }
1561:   MPIU_Allreduce(&minradius, &gminradius, 1, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
1562:   DMPlexSetMinRadius(dm, gminradius);
1563:   /* Compute centroids of ghost cells */
1564:   for (c = cEndInterior; c < cEnd; ++c) {
1565:     PetscFVFaceGeom *fg;
1566:     const PetscInt  *cone,    *support;
1567:     PetscInt         coneSize, supportSize, s;

1569:     DMPlexGetConeSize(dmCell, c, &coneSize);
1570:     if (coneSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Ghost cell %d has cone size %d != 1", c, coneSize);
1571:     DMPlexGetCone(dmCell, c, &cone);
1572:     DMPlexGetSupportSize(dmCell, cone[0], &supportSize);
1573:     if (supportSize != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %d has support size %d != 2", cone[0], supportSize);
1574:     DMPlexGetSupport(dmCell, cone[0], &support);
1575:     DMPlexPointLocalRef(dmFace, cone[0], fgeom, &fg);
1576:     for (s = 0; s < 2; ++s) {
1577:       /* Reflect ghost centroid across plane of face */
1578:       if (support[s] == c) {
1579:         PetscFVCellGeom       *ci;
1580:         PetscFVCellGeom       *cg;
1581:         PetscReal              c2f[3], a;

1583:         DMPlexPointLocalRead(dmCell, support[(s+1)%2], cgeom, &ci);
1584:         DMPlex_WaxpyD_Internal(dim, -1, ci->centroid, fg->centroid, c2f); /* cell to face centroid */
1585:         a    = DMPlex_DotRealD_Internal(dim, c2f, fg->normal)/DMPlex_DotRealD_Internal(dim, fg->normal, fg->normal);
1586:         DMPlexPointLocalRef(dmCell, support[s], cgeom, &cg);
1587:         DMPlex_WaxpyD_Internal(dim, 2*a, fg->normal, ci->centroid, cg->centroid);
1588:         cg->volume = ci->volume;
1589:       }
1590:     }
1591:   }
1592:   VecRestoreArray(*facegeom, &fgeom);
1593:   VecRestoreArray(*cellgeom, &cgeom);
1594:   DMDestroy(&dmCell);
1595:   DMDestroy(&dmFace);
1596:   return(0);
1597: }

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

1604:   Not collective

1606:   Input Argument:
1607: . dm - the DM

1609:   Output Argument:
1610: . minradius - the minium cell radius

1612:   Level: developer

1614: .seealso: DMGetCoordinates()
1615: @*/
1616: PetscErrorCode DMPlexGetMinRadius(DM dm, PetscReal *minradius)
1617: {
1621:   *minradius = ((DM_Plex*) dm->data)->minradius;
1622:   return(0);
1623: }

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

1630:   Logically collective

1632:   Input Arguments:
1633: + dm - the DM
1634: - minradius - the minium cell radius

1636:   Level: developer

1638: .seealso: DMSetCoordinates()
1639: @*/
1640: PetscErrorCode DMPlexSetMinRadius(DM dm, PetscReal minradius)
1641: {
1644:   ((DM_Plex*) dm->data)->minradius = minradius;
1645:   return(0);
1646: }

1650: static PetscErrorCode BuildGradientReconstruction_Internal(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1651: {
1652:   DMLabel        ghostLabel;
1653:   PetscScalar   *dx, *grad, **gref;
1654:   PetscInt       dim, cStart, cEnd, c, cEndInterior, maxNumFaces;

1658:   DMGetDimension(dm, &dim);
1659:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1660:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1661:   DMPlexGetMaxSizes(dm, &maxNumFaces, NULL);
1662:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
1663:   DMGetLabel(dm, "ghost", &ghostLabel);
1664:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
1665:   for (c = cStart; c < cEndInterior; c++) {
1666:     const PetscInt        *faces;
1667:     PetscInt               numFaces, usedFaces, f, d;
1668:     PetscFVCellGeom        *cg;
1669:     PetscBool              boundary;
1670:     PetscInt               ghost;

1672:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1673:     DMPlexGetConeSize(dm, c, &numFaces);
1674:     DMPlexGetCone(dm, c, &faces);
1675:     if (numFaces < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cell %D has only %D faces, not enough for gradient reconstruction", c, numFaces);
1676:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1677:       PetscFVCellGeom       *cg1;
1678:       PetscFVFaceGeom       *fg;
1679:       const PetscInt        *fcells;
1680:       PetscInt               ncell, side;

1682:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
1683:       DMIsBoundaryPoint(dm, faces[f], &boundary);
1684:       if ((ghost >= 0) || boundary) continue;
1685:       DMPlexGetSupport(dm, faces[f], &fcells);
1686:       side  = (c != fcells[0]); /* c is on left=0 or right=1 of face */
1687:       ncell = fcells[!side];    /* the neighbor */
1688:       DMPlexPointLocalRef(dmFace, faces[f], fgeom, &fg);
1689:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
1690:       for (d = 0; d < dim; ++d) dx[usedFaces*dim+d] = cg1->centroid[d] - cg->centroid[d];
1691:       gref[usedFaces++] = fg->grad[side];  /* Gradient reconstruction term will go here */
1692:     }
1693:     if (!usedFaces) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Mesh contains isolated cell (no neighbors). Is it intentional?");
1694:     PetscFVComputeGradient(fvm, usedFaces, dx, grad);
1695:     for (f = 0, usedFaces = 0; f < numFaces; ++f) {
1696:       DMLabelGetValue(ghostLabel, faces[f], &ghost);
1697:       DMIsBoundaryPoint(dm, faces[f], &boundary);
1698:       if ((ghost >= 0) || boundary) continue;
1699:       for (d = 0; d < dim; ++d) gref[usedFaces][d] = grad[usedFaces*dim+d];
1700:       ++usedFaces;
1701:     }
1702:   }
1703:   PetscFree3(dx, grad, gref);
1704:   return(0);
1705: }

1709: static PetscErrorCode BuildGradientReconstruction_Internal_Tree(DM dm, PetscFV fvm, DM dmFace, PetscScalar *fgeom, DM dmCell, PetscScalar *cgeom)
1710: {
1711:   DMLabel        ghostLabel;
1712:   PetscScalar   *dx, *grad, **gref;
1713:   PetscInt       dim, cStart, cEnd, c, cEndInterior, fStart, fEnd, f, nStart, nEnd, maxNumFaces = 0;
1714:   PetscSection   neighSec;
1715:   PetscInt     (*neighbors)[2];
1716:   PetscInt      *counter;

1720:   DMGetDimension(dm, &dim);
1721:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1722:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1723:   if (cEndInterior < 0) {
1724:     cEndInterior = cEnd;
1725:   }
1726:   PetscSectionCreate(PetscObjectComm((PetscObject)dm),&neighSec);
1727:   PetscSectionSetChart(neighSec,cStart,cEndInterior);
1728:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
1729:   DMGetLabel(dm, "ghost", &ghostLabel);
1730:   for (f = fStart; f < fEnd; f++) {
1731:     const PetscInt        *fcells;
1732:     PetscBool              boundary;
1733:     PetscInt               ghost = -1;
1734:     PetscInt               numChildren, numCells, c;

1736:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1737:     DMIsBoundaryPoint(dm, f, &boundary);
1738:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
1739:     if ((ghost >= 0) || boundary || numChildren) continue;
1740:     DMPlexGetSupportSize(dm, f, &numCells);
1741:     if (numCells == 2) {
1742:       DMPlexGetSupport(dm, f, &fcells);
1743:       for (c = 0; c < 2; c++) {
1744:         PetscInt cell = fcells[c];

1746:         if (cell >= cStart && cell < cEndInterior) {
1747:           PetscSectionAddDof(neighSec,cell,1);
1748:         }
1749:       }
1750:     }
1751:   }
1752:   PetscSectionSetUp(neighSec);
1753:   PetscSectionGetMaxDof(neighSec,&maxNumFaces);
1754:   PetscFVLeastSquaresSetMaxFaces(fvm, maxNumFaces);
1755:   nStart = 0;
1756:   PetscSectionGetStorageSize(neighSec,&nEnd);
1757:   PetscMalloc1((nEnd-nStart),&neighbors);
1758:   PetscCalloc1((cEndInterior-cStart),&counter);
1759:   for (f = fStart; f < fEnd; f++) {
1760:     const PetscInt        *fcells;
1761:     PetscBool              boundary;
1762:     PetscInt               ghost = -1;
1763:     PetscInt               numChildren, numCells, c;

1765:     if (ghostLabel) {DMLabelGetValue(ghostLabel, f, &ghost);}
1766:     DMIsBoundaryPoint(dm, f, &boundary);
1767:     DMPlexGetTreeChildren(dm, f, &numChildren, NULL);
1768:     if ((ghost >= 0) || boundary || numChildren) continue;
1769:     DMPlexGetSupportSize(dm, f, &numCells);
1770:     if (numCells == 2) {
1771:       DMPlexGetSupport(dm, f, &fcells);
1772:       for (c = 0; c < 2; c++) {
1773:         PetscInt cell = fcells[c], off;

1775:         if (cell >= cStart && cell < cEndInterior) {
1776:           PetscSectionGetOffset(neighSec,cell,&off);
1777:           off += counter[cell - cStart]++;
1778:           neighbors[off][0] = f;
1779:           neighbors[off][1] = fcells[1 - c];
1780:         }
1781:       }
1782:     }
1783:   }
1784:   PetscFree(counter);
1785:   PetscMalloc3(maxNumFaces*dim, &dx, maxNumFaces*dim, &grad, maxNumFaces, &gref);
1786:   for (c = cStart; c < cEndInterior; c++) {
1787:     PetscInt               numFaces, f, d, off, ghost = -1;
1788:     PetscFVCellGeom        *cg;

1790:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1791:     PetscSectionGetDof(neighSec, c, &numFaces);
1792:     PetscSectionGetOffset(neighSec, c, &off);
1793:     if (ghostLabel) {DMLabelGetValue(ghostLabel, c, &ghost);}
1794:     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);
1795:     for (f = 0; f < numFaces; ++f) {
1796:       PetscFVCellGeom       *cg1;
1797:       PetscFVFaceGeom       *fg;
1798:       const PetscInt        *fcells;
1799:       PetscInt               ncell, side, nface;

1801:       nface = neighbors[off + f][0];
1802:       ncell = neighbors[off + f][1];
1803:       DMPlexGetSupport(dm,nface,&fcells);
1804:       side  = (c != fcells[0]);
1805:       DMPlexPointLocalRef(dmFace, nface, fgeom, &fg);
1806:       DMPlexPointLocalRead(dmCell, ncell, cgeom, &cg1);
1807:       for (d = 0; d < dim; ++d) dx[f*dim+d] = cg1->centroid[d] - cg->centroid[d];
1808:       gref[f] = fg->grad[side];  /* Gradient reconstruction term will go here */
1809:     }
1810:     PetscFVComputeGradient(fvm, numFaces, dx, grad);
1811:     for (f = 0; f < numFaces; ++f) {
1812:       for (d = 0; d < dim; ++d) gref[f][d] = grad[f*dim+d];
1813:     }
1814:   }
1815:   PetscFree3(dx, grad, gref);
1816:   PetscSectionDestroy(&neighSec);
1817:   PetscFree(neighbors);
1818:   return(0);
1819: }

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

1826:   Collective on DM

1828:   Input Arguments:
1829: + dm  - The DM
1830: . fvm - The PetscFV
1831: . faceGeometry - The face geometry from DMPlexGetFaceGeometryFVM()
1832: - cellGeometry - The face geometry from DMPlexGetCellGeometryFVM()

1834:   Output Parameters:
1835: + faceGeometry - The geometric factors for gradient calculation are inserted
1836: - dmGrad - The DM describing the layout of gradient data

1838:   Level: developer

1840: .seealso: DMPlexGetFaceGeometryFVM(), DMPlexGetCellGeometryFVM()
1841: @*/
1842: PetscErrorCode DMPlexComputeGradientFVM(DM dm, PetscFV fvm, Vec faceGeometry, Vec cellGeometry, DM *dmGrad)
1843: {
1844:   DM             dmFace, dmCell;
1845:   PetscScalar   *fgeom, *cgeom;
1846:   PetscSection   sectionGrad, parentSection;
1847:   PetscInt       dim, pdim, cStart, cEnd, cEndInterior, c;

1851:   DMGetDimension(dm, &dim);
1852:   PetscFVGetNumComponents(fvm, &pdim);
1853:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1854:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
1855:   /* Construct the interpolant corresponding to each face from the least-square solution over the cell neighborhood */
1856:   VecGetDM(faceGeometry, &dmFace);
1857:   VecGetDM(cellGeometry, &dmCell);
1858:   VecGetArray(faceGeometry, &fgeom);
1859:   VecGetArray(cellGeometry, &cgeom);
1860:   DMPlexGetTree(dm,&parentSection,NULL,NULL,NULL,NULL);
1861:   if (!parentSection) {
1862:     BuildGradientReconstruction_Internal(dm, fvm, dmFace, fgeom, dmCell, cgeom);
1863:   }
1864:   else {
1865:     BuildGradientReconstruction_Internal_Tree(dm, fvm, dmFace, fgeom, dmCell, cgeom);
1866:   }
1867:   VecRestoreArray(faceGeometry, &fgeom);
1868:   VecRestoreArray(cellGeometry, &cgeom);
1869:   /* Create storage for gradients */
1870:   DMClone(dm, dmGrad);
1871:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &sectionGrad);
1872:   PetscSectionSetChart(sectionGrad, cStart, cEnd);
1873:   for (c = cStart; c < cEnd; ++c) {PetscSectionSetDof(sectionGrad, c, pdim*dim);}
1874:   PetscSectionSetUp(sectionGrad);
1875:   DMSetDefaultSection(*dmGrad, sectionGrad);
1876:   PetscSectionDestroy(&sectionGrad);
1877:   return(0);
1878: }