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) \
 15:   do { \
 16:     PetscInt tmp = (a); \
 17:     (a)          = (b); \
 18:     (b)          = tmp; \
 19:   } while (0)
 20:   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
 21: #undef SWAP
 22:   return 0;
 23: }

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

 35:   PetscOptionsGetInt(NULL, ((PetscObject)boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
 36:   PetscObjectGetComm((PetscObject)boundary, &comm);
 37:   MPI_Comm_rank(comm, &rank);
 38:   DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
 39:   DMUniversalLabelCreate(boundary, &universal);
 40:   DMLabelGetDefaultValue(universal->label, &defVal);

 42:   PLCCreate(&in);
 43:   PLCCreate(&out);

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

 52:     PetscMalloc1(in->numberofpoints * dim, &in->pointlist);
 53:     PetscCalloc1(in->numberofpoints, &in->pointmarkerlist);
 54:     DMGetCoordinatesLocal(boundary, &coordinates);
 55:     DMGetCoordinateSection(boundary, &coordSection);
 56:     VecGetArrayRead(coordinates, &array);
 57:     for (v = vStart; v < vEnd; ++v) {
 58:       const PetscInt idx = v - vStart;
 59:       PetscInt       off, d, m;

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

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

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

 84:       DMLabelGetValue(universal->label, e, &val);
 85:       if (val != defVal) in->edgemarkerlist[idx] = (int)val;
 86:     }
 87:   }

 89:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
 90:   in->numberoffacets = fEnd - fStart;
 91:   if (in->numberoffacets > 0) {
 92:     PetscMalloc1(in->numberoffacets, &in->facetlist);
 93:     PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
 94:     for (f = fStart; f < fEnd; ++f) {
 95:       const PetscInt idx    = f - fStart;
 96:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
 97:       polygon       *poly;

 99:       in->facetlist[idx].numberofpolygons = 1;
100:       PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
101:       in->facetlist[idx].numberofholes = 0;
102:       in->facetlist[idx].holelist      = NULL;

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

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

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

140:     TetGenCheckOpts(&t);
141:     TetGenTetrahedralize(&t, in, out);
142:   }
143:   {
144:     const PetscInt numCorners  = 4;
145:     const PetscInt numCells    = out->numberoftetrahedra;
146:     const PetscInt numVertices = out->numberofpoints;
147:     PetscReal     *meshCoords  = NULL;
148:     PetscInt      *cells       = NULL;

150:     if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
151:       meshCoords = (PetscReal *)out->pointlist;
152:     } else {
153:       PetscInt i;

155:       PetscMalloc1(dim * numVertices, &meshCoords);
156:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out->pointlist[i];
157:     }
158:     if (sizeof(PetscInt) == sizeof(out->tetrahedronlist[0])) {
159:       cells = (PetscInt *)out->tetrahedronlist;
160:     } else {
161:       PetscInt i;

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

167:     DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
168:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
169:     if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscFree(meshCoords);
170:     if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscFree(cells);

172:     /* Set labels */
173:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
174:     for (v = 0; v < numVertices; ++v) {
175:       if (out->pointmarkerlist[v]) DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]);
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:           /* Determine 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)
244:               for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
245:             DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
246:           } else DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
247:           for (b = 0; b < Nb; ++b) {
248:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
249:           }
250:           if (b < Nb) {
251:             PetscInt  cval    = b, eVal, fVal;
252:             PetscInt *closure = NULL, Ncl, cl;

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

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

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

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

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

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:     PetscCalloc1(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:       if (m != defVal) 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:       if (val != defVal) 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:       if (val != defVal) 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]) DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]);
450:     }
451:     if (interpolate) {
452:       PetscInt e, f;

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

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

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

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

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

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

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

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

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

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

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