Actual source code: tetgenerate.cxx

petsc-master 2021-02-25
Report Typos and Errors
  1: #include <petsc/private/dmpleximpl.h>

  3: #ifdef PETSC_HAVE_EGADS
  4: #include <egads.h>
  5: #endif

  7: #if defined(PETSC_HAVE_TETGEN_TETLIBRARY_NEEDED)
  8: #define TETLIBRARY
  9: #endif
 10: #include <tetgen.h>

 12: /* This is to fix the tetrahedron orientation from TetGen */
 13: static PetscErrorCode DMPlexInvertCells_Tetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
 14: {
 15:   PetscInt bound = numCells*numCorners, coff;

 18: #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
 19:   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
 20: #undef SWAP
 21:   return(0);
 22: }

 24: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
 25: {
 26:   MPI_Comm               comm;
 27:   const PetscInt         dim = 3;
 28:   ::tetgenio             in;
 29:   ::tetgenio             out;
 30:   DMUniversalLabel       universal;
 31:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f;
 32:   DMPlexInterpolatedFlag isInterpolated;
 33:   PetscMPIInt            rank;
 34:   PetscErrorCode         ierr;

 37:   PetscObjectGetComm((PetscObject)boundary,&comm);
 38:   MPI_Comm_rank(comm, &rank);
 39:   DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
 40:   DMUniversalLabelCreate(boundary, &universal);

 42:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
 43:   in.numberofpoints = vEnd - vStart;
 44:   if (in.numberofpoints > 0) {
 45:     PetscSection       coordSection;
 46:     Vec                coordinates;
 47:     const PetscScalar *array;

 49:     in.pointlist       = new double[in.numberofpoints*dim];
 50:     in.pointmarkerlist = new int[in.numberofpoints];

 52:     DMGetCoordinatesLocal(boundary, &coordinates);
 53:     DMGetCoordinateSection(boundary, &coordSection);
 54:     VecGetArrayRead(coordinates, &array);
 55:     for (v = vStart; v < vEnd; ++v) {
 56:       const PetscInt idx = v - vStart;
 57:       PetscInt       off, d, val;

 59:       PetscSectionGetOffset(coordSection, v, &off);
 60:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
 61:       DMLabelGetValue(universal->label, v, &val);
 62:       in.pointmarkerlist[idx] = (int) val;
 63:     }
 64:     VecRestoreArrayRead(coordinates, &array);
 65:   }

 67:   DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
 68:   in.numberofedges = eEnd - eStart;
 69:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
 70:     in.edgelist       = new int[in.numberofedges * 2];
 71:     in.edgemarkerlist = new int[in.numberofedges];
 72:     for (e = eStart; e < eEnd; ++e) {
 73:       const PetscInt  idx = e - eStart;
 74:       const PetscInt *cone;
 75:       PetscInt        coneSize, val;

 77:       DMPlexGetConeSize(boundary, e, &coneSize);
 78:       DMPlexGetCone(boundary, e, &cone);
 79:       in.edgelist[idx*2]     = cone[0] - vStart;
 80:       in.edgelist[idx*2 + 1] = cone[1] - vStart;

 82:       DMLabelGetValue(universal->label, e, &val);
 83:       in.edgemarkerlist[idx] = (int) val;
 84:     }
 85:   }

 87:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
 88:   in.numberoffacets = fEnd - fStart;
 89:   if (in.numberoffacets > 0) {
 90:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
 91:     in.facetmarkerlist = new int[in.numberoffacets];
 92:     for (f = fStart; f < fEnd; ++f) {
 93:       const PetscInt idx    = f - fStart;
 94:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val = -1;

 96:       in.facetlist[idx].numberofpolygons = 1;
 97:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
 98:       in.facetlist[idx].numberofholes    = 0;
 99:       in.facetlist[idx].holelist         = NULL;

101:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
102:       for (p = 0; p < numPoints*2; p += 2) {
103:         const PetscInt point = points[p];
104:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
105:       }

107:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
108:       poly->numberofvertices = numVertices;
109:       poly->vertexlist       = new int[poly->numberofvertices];
110:       for (v = 0; v < numVertices; ++v) {
111:         const PetscInt vIdx = points[v] - vStart;
112:         poly->vertexlist[v] = vIdx;
113:       }
114:       DMLabelGetValue(universal->label, f, &val);
115:       in.facetmarkerlist[idx] = (int) val;
116:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
117:     }
118:   }
119:   if (!rank) {
120:     DM_Plex *mesh = (DM_Plex *) boundary->data;
121:     char     args[32];

123:     /* Take away 'Q' for verbose output */
124: #ifdef PETSC_HAVE_EGADS
125:     PetscStrcpy(args, "pqezQY");
126: #else
127:     PetscStrcpy(args, "pqezQ");
128: #endif
129:     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
130:     else                  {::tetrahedralize(args, &in, &out);}
131:   }
132:   {
133:     const PetscInt   numCorners  = 4;
134:     const PetscInt   numCells    = out.numberoftetrahedra;
135:     const PetscInt   numVertices = out.numberofpoints;
136:     PetscReal        *meshCoords = NULL;
137:     PetscInt         *cells      = NULL;

139:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
140:       meshCoords = (PetscReal *) out.pointlist;
141:     } else {
142:       PetscInt i;

144:       meshCoords = new PetscReal[dim * numVertices];
145:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
146:     }
147:     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
148:       cells = (PetscInt *) out.tetrahedronlist;
149:     } else {
150:       PetscInt i;

152:       cells = new PetscInt[numCells * numCorners];
153:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out.tetrahedronlist[i];
154:     }

156:     DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
157:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);

159:     /* Set labels */
160:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
161:     for (v = 0; v < numVertices; ++v) {
162:       if (out.pointmarkerlist[v]) {
163:         DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);
164:       }
165:     }
166:     if (interpolate) {
167:       PetscInt e;

169:       for (e = 0; e < out.numberofedges; e++) {
170:         if (out.edgemarkerlist[e]) {
171:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
172:           const PetscInt *edges;
173:           PetscInt        numEdges;

175:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
176:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
177:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
178:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
179:         }
180:       }
181:       for (f = 0; f < out.numberoftrifaces; f++) {
182:         if (out.trifacemarkerlist[f]) {
183:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
184:           const PetscInt *faces;
185:           PetscInt        numFaces;

187:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
188:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
189:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
190:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
191:         }
192:       }
193:     }

195: #ifdef PETSC_HAVE_EGADS
196:     {
197:       DMLabel        bodyLabel;
198:       PetscContainer modelObj;
199:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
200:       ego           *bodies;
201:       ego            model, geom;
202:       int            Nb, oclass, mtype, *senses;

204:       /* Get Attached EGADS Model from Original DMPlex */
205:       PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);
206:       PetscContainerGetPointer(modelObj, (void **) &model);
207:       EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
208:       /* Transfer EGADS Model to Volumetric Mesh */
209:       PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);

211:       /* Set Cell Labels */
212:       DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
213:       DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
214:       DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
215:       DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);

217:       for (c = cStart; c < cEnd; ++c) {
218:         PetscReal centroid[3] = {0., 0., 0.};
219:         PetscInt  b;

221:         /* Deterimine what body the cell's centroid is located in */
222:         if (!interpolate) {
223:           PetscSection   coordSection;
224:           Vec            coordinates;
225:           PetscScalar   *coords = NULL;
226:           PetscInt       coordSize, s, d;

228:           DMGetCoordinatesLocal(*dm, &coordinates);
229:           DMGetCoordinateSection(*dm, &coordSection);
230:           DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
231:           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
232:           DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
233:         } else {
234:           DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
235:         }
236:         for (b = 0; b < Nb; ++b) {
237:           if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
238:         }
239:         if (b < Nb) {
240:           PetscInt   cval = b, eVal, fVal;
241:           PetscInt *closure = NULL, Ncl, cl;

243:           DMLabelSetValue(bodyLabel, c, cval);
244:           DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
245:           for (cl = 0; cl < Ncl; cl += 2) {
246:             const PetscInt p = closure[cl];

248:             if (p >= eStart && p < eEnd) {
249:               DMLabelGetValue(bodyLabel, p, &eVal);
250:               if (eVal < 0) {DMLabelSetValue(bodyLabel, p, cval);}
251:             }
252:             if (p >= fStart && p < fEnd) {
253:               DMLabelGetValue(bodyLabel, p, &fVal);
254:               if (fVal < 0) {DMLabelSetValue(bodyLabel, p, cval);}
255:             }
256:           }
257:           DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
258:         }
259:       }
260:     }
261: #endif
262:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
263:   }
264:   return(0);
265: }

267: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
268: {
269:   MPI_Comm               comm;
270:   const PetscInt         dim = 3;
271:   ::tetgenio             in;
272:   ::tetgenio             out;
273:   DMUniversalLabel       universal;
274:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c;
275:   DMPlexInterpolatedFlag isInterpolated;
276:   PetscMPIInt            rank;
277:   PetscErrorCode         ierr;

280:   PetscObjectGetComm((PetscObject)dm,&comm);
281:   MPI_Comm_rank(comm, &rank);
282:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
283:   DMUniversalLabelCreate(dm, &universal);

285:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
286:   in.numberofpoints = vEnd - vStart;
287:   if (in.numberofpoints > 0) {
288:     PetscSection coordSection;
289:     Vec          coordinates;
290:     PetscScalar *array;

292:     in.pointlist       = new double[in.numberofpoints*dim];
293:     in.pointmarkerlist = new int[in.numberofpoints];

295:     DMGetCoordinatesLocal(dm, &coordinates);
296:     DMGetCoordinateSection(dm, &coordSection);
297:     VecGetArray(coordinates, &array);
298:     for (v = vStart; v < vEnd; ++v) {
299:       const PetscInt idx = v - vStart;
300:       PetscInt       off, d, val;

302:       PetscSectionGetOffset(coordSection, v, &off);
303:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
304:       DMLabelGetValue(universal->label, v, &val);
305:       in.pointmarkerlist[idx] = (int) val;
306:     }
307:     VecRestoreArray(coordinates, &array);
308:   }

310:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
311:   in.numberofedges = eEnd - eStart;
312:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
313:     in.edgelist       = new int[in.numberofedges * 2];
314:     in.edgemarkerlist = new int[in.numberofedges];
315:     for (e = eStart; e < eEnd; ++e) {
316:       const PetscInt  idx = e - eStart;
317:       const PetscInt *cone;
318:       PetscInt        coneSize, val;

320:       DMPlexGetConeSize(dm, e, &coneSize);
321:       DMPlexGetCone(dm, e, &cone);
322:       in.edgelist[idx*2]     = cone[0] - vStart;
323:       in.edgelist[idx*2 + 1] = cone[1] - vStart;

325:       DMLabelGetValue(universal->label, e, &val);
326:       in.edgemarkerlist[idx] = (int) val;
327:     }
328:   }

330:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
331:   in.numberoffacets = fEnd - fStart;
332:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
333:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
334:     in.facetmarkerlist = new int[in.numberoffacets];
335:     for (f = fStart; f < fEnd; ++f) {
336:       const PetscInt idx    = f - fStart;
337:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val;

339:       in.facetlist[idx].numberofpolygons = 1;
340:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
341:       in.facetlist[idx].numberofholes    = 0;
342:       in.facetlist[idx].holelist         = NULL;

344:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
345:       for (p = 0; p < numPoints*2; p += 2) {
346:         const PetscInt point = points[p];
347:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
348:       }

350:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
351:       poly->numberofvertices = numVertices;
352:       poly->vertexlist       = new int[poly->numberofvertices];
353:       for (v = 0; v < numVertices; ++v) {
354:         const PetscInt vIdx = points[v] - vStart;
355:         poly->vertexlist[v] = vIdx;
356:       }

358:       DMLabelGetValue(universal->label, f, &val);
359:       in.facetmarkerlist[idx] = (int) val;

361:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
362:     }
363:   }

365:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
366:   in.numberofcorners       = 4;
367:   in.numberoftetrahedra    = cEnd - cStart;
368:   in.tetrahedronvolumelist = (double *) maxVolumes;
369:   if (in.numberoftetrahedra > 0) {
370:     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
371:     for (c = cStart; c < cEnd; ++c) {
372:       const PetscInt idx     = c - cStart;
373:       PetscInt      *closure = NULL;
374:       PetscInt       closureSize;

376:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
377:       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
378:       for (v = 0; v < 4; ++v) in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
379:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
380:     }
381:   }

383:   if (!rank) {
384:     char args[32];

386:     /* Take away 'Q' for verbose output */
387:     PetscStrcpy(args, "qezQra");
388:     ::tetrahedralize(args, &in, &out);
389:   }

391:   in.tetrahedronvolumelist = NULL;
392:   {
393:     const PetscInt   numCorners  = 4;
394:     const PetscInt   numCells    = out.numberoftetrahedra;
395:     const PetscInt   numVertices = out.numberofpoints;
396:     PetscReal        *meshCoords = NULL;
397:     PetscInt         *cells      = NULL;
398:     PetscBool        interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

400:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
401:       meshCoords = (PetscReal *) out.pointlist;
402:     } else {
403:       PetscInt i;

405:       meshCoords = new PetscReal[dim * numVertices];
406:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
407:     }
408:     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
409:       cells = (PetscInt *) out.tetrahedronlist;
410:     } else {
411:       PetscInt i;

413:       cells = new PetscInt[numCells * numCorners];
414:       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
415:     }

417:     DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
418:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
419:     if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {delete [] meshCoords;}
420:     if (sizeof (PetscInt) != sizeof (out.tetrahedronlist[0])) {delete [] cells;}

422:     /* Set labels */
423:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
424:     for (v = 0; v < numVertices; ++v) {
425:       if (out.pointmarkerlist[v]) {
426:         DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);
427:       }
428:     }
429:     if (interpolate) {
430:       PetscInt e, f;

432:       for (e = 0; e < out.numberofedges; ++e) {
433:         if (out.edgemarkerlist[e]) {
434:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
435:           const PetscInt *edges;
436:           PetscInt        numEdges;

438:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
439:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
440:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
441:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
442:         }
443:       }
444:       for (f = 0; f < out.numberoftrifaces; ++f) {
445:         if (out.trifacemarkerlist[f]) {
446:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
447:           const PetscInt *faces;
448:           PetscInt        numFaces;

450:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
451:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
452:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
453:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
454:         }
455:       }
456:     }

458: #ifdef PETSC_HAVE_EGADS
459:     {
460:       DMLabel        bodyLabel;
461:       PetscContainer modelObj;
462:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
463:       ego           *bodies;
464:       ego            model, geom;
465:       int            Nb, oclass, mtype, *senses;

467:       /* Get Attached EGADS Model from Original DMPlex */
468:       PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
469:       PetscContainerGetPointer(modelObj, (void **) &model);
470:       EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
471:       /* Transfer EGADS Model to Volumetric Mesh */
472:       PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);

474:       /* Set Cell Labels */
475:       DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
476:       DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
477:       DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
478:       DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);

480:       for (c = cStart; c < cEnd; ++c) {
481:         PetscReal centroid[3] = {0., 0., 0.};
482:         PetscInt  b;

484:         /* Deterimine what body the cell's centroid is located in */
485:         if (!interpolate) {
486:           PetscSection   coordSection;
487:           Vec            coordinates;
488:           PetscScalar   *coords = NULL;
489:           PetscInt       coordSize, s, d;

491:           DMGetCoordinatesLocal(*dmRefined, &coordinates);
492:           DMGetCoordinateSection(*dmRefined, &coordSection);
493:           DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
494:           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
495:           DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
496:         } else {
497:           DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
498:         }
499:         for (b = 0; b < Nb; ++b) {
500:           if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
501:         }
502:         if (b < Nb) {
503:           PetscInt   cval = b, eVal, fVal;
504:           PetscInt *closure = NULL, Ncl, cl;

506:           DMLabelSetValue(bodyLabel, c, cval);
507:           DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
508:           for (cl = 0; cl < Ncl; cl += 2) {
509:             const PetscInt p = closure[cl];

511:             if (p >= eStart && p < eEnd) {
512:               DMLabelGetValue(bodyLabel, p, &eVal);
513:               if (eVal < 0) {DMLabelSetValue(bodyLabel, p, cval);}
514:             }
515:             if (p >= fStart && p < fEnd) {
516:               DMLabelGetValue(bodyLabel, p, &fVal);
517:               if (fVal < 0) {DMLabelSetValue(bodyLabel, p, cval);}
518:             }
519:           }
520:           DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
521:         }
522:       }
523:     }
524: #endif
525:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
526:   }
527:   return(0);
528: }