Actual source code: ctetgenerate.c

  1: #include <petsc/private/dmpleximpl.h>

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

  7: #include <ctetgen.h>

  9: /* This is to fix the tetrahedron orientation from TetGen */
 10: static PetscErrorCode DMPlexInvertCells_CTetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
 11: {
 12:   PetscInt bound = numCells*numCorners, coff;

 14: #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
 15:   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
 16: #undef SWAP
 17:   return 0;
 18: }

 20: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
 21: {
 22:   MPI_Comm               comm;
 23:   const PetscInt         dim = 3;
 24:   PLC                   *in, *out;
 25:   DMUniversalLabel       universal;
 26:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, verbose = 0;
 27:   DMPlexInterpolatedFlag isInterpolated;
 28:   PetscMPIInt            rank;

 30:   PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
 31:   PetscObjectGetComm((PetscObject)boundary,&comm);
 32:   MPI_Comm_rank(comm, &rank);
 33:   DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
 34:   DMUniversalLabelCreate(boundary, &universal);

 36:   PLCCreate(&in);
 37:   PLCCreate(&out);

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

 46:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
 47:     PetscMalloc1(in->numberofpoints,     &in->pointmarkerlist);
 48:     DMGetCoordinatesLocal(boundary, &coordinates);
 49:     DMGetCoordinateSection(boundary, &coordSection);
 50:     VecGetArrayRead(coordinates, &array);
 51:     for (v = vStart; v < vEnd; ++v) {
 52:       const PetscInt idx = v - vStart;
 53:       PetscInt       off, d, m;

 55:       PetscSectionGetOffset(coordSection, v, &off);
 56:       for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
 57:       DMLabelGetValue(universal->label, v, &m);
 58:       in->pointmarkerlist[idx] = (int) m;
 59:     }
 60:     VecRestoreArrayRead(coordinates, &array);
 61:   }

 63:   DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
 64:   in->numberofedges = eEnd - eStart;
 65:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
 66:     PetscMalloc1(in->numberofedges*2, &in->edgelist);
 67:     PetscMalloc1(in->numberofedges,   &in->edgemarkerlist);
 68:     for (e = eStart; e < eEnd; ++e) {
 69:       const PetscInt  idx = e - eStart;
 70:       const PetscInt *cone;
 71:       PetscInt        coneSize, val;

 73:       DMPlexGetConeSize(boundary, e, &coneSize);
 74:       DMPlexGetCone(boundary, e, &cone);
 75:       in->edgelist[idx*2]     = cone[0] - vStart;
 76:       in->edgelist[idx*2 + 1] = cone[1] - vStart;

 78:       DMLabelGetValue(universal->label, e, &val);
 79:       in->edgemarkerlist[idx] = (int) val;
 80:     }
 81:   }

 83:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
 84:   in->numberoffacets = fEnd - fStart;
 85:   if (in->numberoffacets > 0) {
 86:     PetscMalloc1(in->numberoffacets, &in->facetlist);
 87:     PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
 88:     for (f = fStart; f < fEnd; ++f) {
 89:       const PetscInt idx    = f - fStart;
 90:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
 91:       polygon       *poly;

 93:       in->facetlist[idx].numberofpolygons = 1;
 94:       PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
 95:       in->facetlist[idx].numberofholes    = 0;
 96:       in->facetlist[idx].holelist         = NULL;

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

104:       poly                   = in->facetlist[idx].polygonlist;
105:       poly->numberofvertices = numVertices;
106:       PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
107:       for (v = 0; v < numVertices; ++v) {
108:         const PetscInt vIdx = points[v] - vStart;
109:         poly->vertexlist[v] = vIdx;
110:       }
111:       DMLabelGetValue(universal->label, f, &m);
112:       in->facetmarkerlist[idx] = (int) m;
113:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
114:     }
115:   }
116:   if (rank == 0) {
117:     TetGenOpts t;

119:     TetGenOptsInitialize(&t);
120:     t.in        = boundary; /* Should go away */
121:     t.plc       = 1;
122:     t.quality   = 1;
123:     t.edgesout  = 1;
124:     t.zeroindex = 1;
125:     t.quiet     = 1;
126:     t.verbose   = verbose;
127: #if 0
128: #ifdef PETSC_HAVE_EGADS
129:     /* Need to add in more TetGen code */
130:     t.nobisect  = 1; /* Add Y to preserve Surface Mesh for EGADS */
131: #endif
132: #endif

134:     TetGenCheckOpts(&t);
135:     TetGenTetrahedralize(&t, in, out);
136:   }
137:   {
138:     const PetscInt numCorners  = 4;
139:     const PetscInt numCells    = out->numberoftetrahedra;
140:     const PetscInt numVertices = out->numberofpoints;
141:     PetscReal      *meshCoords = NULL;
142:     PetscInt       *cells      = NULL;

144:     if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
145:       meshCoords = (PetscReal *) out->pointlist;
146:     } else {
147:       PetscInt i;

149:       PetscMalloc1(dim * numVertices, &meshCoords);
150:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
151:     }
152:     if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
153:       cells = (PetscInt *) out->tetrahedronlist;
154:     } else {
155:       PetscInt i;

157:       PetscMalloc1(numCells * numCorners, &cells);
158:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out->tetrahedronlist[i];
159:     }

161:     DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
162:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
163:     if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {
164:       PetscFree(meshCoords);
165:     }
166:     if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {
167:       PetscFree(cells);
168:     }

170:     /* Set labels */
171:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
172:     for (v = 0; v < numVertices; ++v) {
173:       if (out->pointmarkerlist[v]) {
174:         DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
175:       }
176:     }
177:     if (interpolate) {
178:       PetscInt e;

180:       for (e = 0; e < out->numberofedges; e++) {
181:         if (out->edgemarkerlist[e]) {
182:           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
183:           const PetscInt *edges;
184:           PetscInt        numEdges;

186:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
188:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
189:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
190:         }
191:       }
192:       for (f = 0; f < out->numberoftrifaces; f++) {
193:         if (out->trifacemarkerlist[f]) {
194:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
195:           const PetscInt *faces;
196:           PetscInt        numFaces;

198:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
200:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
201:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
202:         }
203:       }
204:     }

206: #ifdef PETSC_HAVE_EGADS
207:     {
208:       DMLabel        bodyLabel;
209:       PetscContainer modelObj;
210:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
211:       ego           *bodies;
212:       ego            model, geom;
213:       int            Nb, oclass, mtype, *senses;

215:       /* Get Attached EGADS Model from Original DMPlex */
216:       PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);
217:       if (modelObj) {
218:         PetscContainerGetPointer(modelObj, (void **) &model);
219:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
220:         /* Transfer EGADS Model to Volumetric Mesh */
221:         PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);

223:         /* Set Cell Labels */
224:         DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
225:         DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
226:         DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
227:         DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);

229:         for (c = cStart; c < cEnd; ++c) {
230:           PetscReal centroid[3] = {0., 0., 0.};
231:           PetscInt  b;

233:           /* Deterimine what body the cell's centroid is located in */
234:           if (!interpolate) {
235:             PetscSection   coordSection;
236:             Vec            coordinates;
237:             PetscScalar   *coords = NULL;
238:             PetscInt       coordSize, s, d;

240:             DMGetCoordinatesLocal(*dm, &coordinates);
241:             DMGetCoordinateSection(*dm, &coordSection);
242:             DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
243:             for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
244:             DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
245:           } else {
246:             DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
247:           }
248:           for (b = 0; b < Nb; ++b) {
249:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
250:           }
251:           if (b < Nb) {
252:             PetscInt   cval = b, eVal, fVal;
253:             PetscInt *closure = NULL, Ncl, cl;

255:             DMLabelSetValue(bodyLabel, c, cval);
256:             DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
257:             for (cl = 0; cl < Ncl; ++cl) {
258:               const PetscInt p = closure[cl*2];

260:               if (p >= eStart && p < eEnd) {
261:                 DMLabelGetValue(bodyLabel, p, &eVal);
262:                 if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
263:               }
264:               if (p >= fStart && p < fEnd) {
265:                 DMLabelGetValue(bodyLabel, p, &fVal);
266:                 if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
267:               }
268:             }
269:             DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
270:           }
271:         }
272:       }
273:     }
274: #endif
275:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
276:   }

278:   DMUniversalLabelDestroy(&universal);
279:   PLCDestroy(&in);
280:   PLCDestroy(&out);
281:   return 0;
282: }

284: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
285: {
286:   MPI_Comm               comm;
287:   const PetscInt         dim = 3;
288:   PLC                   *in, *out;
289:   DMUniversalLabel       universal;
290:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, verbose = 0;
291:   DMPlexInterpolatedFlag isInterpolated;
292:   PetscMPIInt            rank;

294:   PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
295:   PetscObjectGetComm((PetscObject)dm,&comm);
296:   MPI_Comm_rank(comm, &rank);
297:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
298:   DMUniversalLabelCreate(dm, &universal);

300:   PLCCreate(&in);
301:   PLCCreate(&out);

303:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
304:   in->numberofpoints = vEnd - vStart;
305:   if (in->numberofpoints > 0) {
306:     PetscSection coordSection;
307:     Vec          coordinates;
308:     PetscScalar *array;

310:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
311:     PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
312:     DMGetCoordinatesLocal(dm, &coordinates);
313:     DMGetCoordinateSection(dm, &coordSection);
314:     VecGetArray(coordinates, &array);
315:     for (v = vStart; v < vEnd; ++v) {
316:       const PetscInt idx = v - vStart;
317:       PetscInt       off, d, m;

319:       PetscSectionGetOffset(coordSection, v, &off);
320:       for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
321:       DMLabelGetValue(universal->label, v, &m);
322:       in->pointmarkerlist[idx] = (int) m;
323:     }
324:     VecRestoreArray(coordinates, &array);
325:   }

327:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
328:   in->numberofedges = eEnd - eStart;
329:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
330:     PetscMalloc1(in->numberofedges * 2, &in->edgelist);
331:     PetscMalloc1(in->numberofedges,     &in->edgemarkerlist);
332:     for (e = eStart; e < eEnd; ++e) {
333:       const PetscInt  idx = e - eStart;
334:       const PetscInt *cone;
335:       PetscInt        coneSize, val;

337:       DMPlexGetConeSize(dm, e, &coneSize);
338:       DMPlexGetCone(dm, e, &cone);
339:       in->edgelist[idx*2]     = cone[0] - vStart;
340:       in->edgelist[idx*2 + 1] = cone[1] - vStart;

342:       DMLabelGetValue(universal->label, e, &val);
343:       in->edgemarkerlist[idx] = (int) val;
344:     }
345:   }

347:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
348:   in->numberoftrifaces = 0;
349:   for (f = fStart; f < fEnd; ++f) {
350:     PetscInt supportSize;

352:     DMPlexGetSupportSize(dm, f, &supportSize);
353:     if (supportSize == 1) ++in->numberoftrifaces;
354:   }
355:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
356:     PetscInt tf = 0;

358:     PetscMalloc1(in->numberoftrifaces*3, &in->trifacelist);
359:     PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist);
360:     for (f = fStart; f < fEnd; ++f) {
361:       PetscInt *points = NULL;
362:       PetscInt supportSize, numPoints, p, Nv = 0, val;

364:       DMPlexGetSupportSize(dm, f, &supportSize);
365:       if (supportSize != 1) continue;
366:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
367:       for (p = 0; p < numPoints*2; p += 2) {
368:         const PetscInt point = points[p];
369:         if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf*3 + Nv++] = point - vStart;
370:       }
371:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
373:       DMLabelGetValue(universal->label, f, &val);
374:       in->trifacemarkerlist[tf] = (int) val;
375:       ++tf;
376:     }
377:   }

379:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
380:   in->numberofcorners       = 4;
381:   in->numberoftetrahedra    = cEnd - cStart;
382:   in->tetrahedronvolumelist = maxVolumes;
383:   if (in->numberoftetrahedra > 0) {
384:     PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
385:     for (c = cStart; c < cEnd; ++c) {
386:       const PetscInt idx     = c - cStart;
387:       PetscInt      *closure = NULL;
388:       PetscInt       closureSize;

390:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
392:       for (v = 0; v < 4; ++v) in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
393:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
394:     }
395:   }

397:   if (rank == 0) {
398:     TetGenOpts t;

400:     TetGenOptsInitialize(&t);

402:     t.in        = dm; /* Should go away */
403:     t.refine    = 1;
404:     t.varvolume = 1;
405:     t.quality   = 1;
406:     t.edgesout  = 1;
407:     t.zeroindex = 1;
408:     t.quiet     = 1;
409:     t.verbose   = verbose; /* Change this */

411:     TetGenCheckOpts(&t);
412:     TetGenTetrahedralize(&t, in, out);
413:   }

415:   in->tetrahedronvolumelist = NULL;
416:   {
417:     const PetscInt numCorners  = 4;
418:     const PetscInt numCells    = out->numberoftetrahedra;
419:     const PetscInt numVertices = out->numberofpoints;
420:     PetscReal      *meshCoords = NULL;
421:     PetscInt       *cells      = NULL;
422:     PetscBool      interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

424:     if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
425:       meshCoords = (PetscReal *) out->pointlist;
426:     } else {
427:       PetscInt i;

429:       PetscMalloc1(dim * numVertices, &meshCoords);
430:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
431:     }
432:     if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
433:       cells = (PetscInt *) out->tetrahedronlist;
434:     } else {
435:       PetscInt i;

437:       PetscMalloc1(numCells * numCorners, &cells);
438:       for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt) out->tetrahedronlist[i];
439:     }

441:     DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
442:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
443:     if (sizeof (PetscReal) != sizeof (out->pointlist[0])) PetscFree(meshCoords);
444:     if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) PetscFree(cells);

446:     /* Set labels */
447:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
448:     for (v = 0; v < numVertices; ++v) {
449:       if (out->pointmarkerlist[v]) {
450:         DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
451:       }
452:     }
453:     if (interpolate) {
454:       PetscInt e, f;

456:       for (e = 0; e < out->numberofedges; e++) {
457:         if (out->edgemarkerlist[e]) {
458:           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
459:           const PetscInt *edges;
460:           PetscInt        numEdges;

462:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
464:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
465:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
466:         }
467:       }
468:       for (f = 0; f < out->numberoftrifaces; f++) {
469:         if (out->trifacemarkerlist[f]) {
470:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
471:           const PetscInt *faces;
472:           PetscInt        numFaces;

474:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
476:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
477:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
478:         }
479:       }
480:     }

482: #ifdef PETSC_HAVE_EGADS
483:     {
484:       DMLabel        bodyLabel;
485:       PetscContainer modelObj;
486:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
487:       ego           *bodies;
488:       ego            model, geom;
489:       int            Nb, oclass, mtype, *senses;

491:       /* Get Attached EGADS Model from Original DMPlex */
492:       PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
493:       if (modelObj) {
494:         PetscContainerGetPointer(modelObj, (void **) &model);
495:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
496:         /* Transfer EGADS Model to Volumetric Mesh */
497:         PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);

499:         /* Set Cell Labels */
500:         DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
501:         DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
502:         DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
503:         DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);

505:         for (c = cStart; c < cEnd; ++c) {
506:           PetscReal centroid[3] = {0., 0., 0.};
507:           PetscInt  b;

509:           /* Deterimine what body the cell's centroid is located in */
510:           if (!interpolate) {
511:             PetscSection   coordSection;
512:             Vec            coordinates;
513:             PetscScalar   *coords = NULL;
514:             PetscInt       coordSize, s, d;

516:             DMGetCoordinatesLocal(*dmRefined, &coordinates);
517:             DMGetCoordinateSection(*dmRefined, &coordSection);
518:             DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
519:             for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
520:             DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
521:           } else {
522:             DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
523:           }
524:           for (b = 0; b < Nb; ++b) {
525:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
526:           }
527:           if (b < Nb) {
528:             PetscInt   cval = b, eVal, fVal;
529:             PetscInt *closure = NULL, Ncl, cl;

531:             DMLabelSetValue(bodyLabel, c, cval);
532:             DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
533:             for (cl = 0; cl < Ncl; cl += 2) {
534:               const PetscInt p = closure[cl];

536:               if (p >= eStart && p < eEnd) {
537:                 DMLabelGetValue(bodyLabel, p, &eVal);
538:                 if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
539:               }
540:               if (p >= fStart && p < fEnd) {
541:                 DMLabelGetValue(bodyLabel, p, &fVal);
542:                 if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
543:               }
544:             }
545:             DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
546:           }
547:         }
548:       }
549:     }
550: #endif
551:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
552:   }
553:   DMUniversalLabelDestroy(&universal);
554:   PLCDestroy(&in);
555:   PLCDestroy(&out);
556:   return 0;
557: }