Actual source code: tetgenerate.cxx

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

  3: #ifdef PETSC_HAVE_EGADS
  4:   #include <egads.h>
  5: /* Need to make EGADSLite header compatible */
  6: extern "C" int EGlite_getTopology(const ego, ego *, int *, int *, double *, int *, ego **, int **);
  7: extern "C" int EGlite_inTopology(const ego, const double *);
  8: #endif

 10: #if defined(PETSC_HAVE_TETGEN_TETLIBRARY_NEEDED)
 11:   #define TETLIBRARY
 12: #endif
 13: #include <tetgen.h>

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

 20: #define SWAP(a, b) \
 21:   do { \
 22:     PetscInt tmp = (a); \
 23:     (a)          = (b); \
 24:     (b)          = tmp; \
 25:   } while (0)
 26:   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
 27: #undef SWAP
 28:   return 0;
 29: }

 31: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
 32: {
 33:   MPI_Comm               comm;
 34:   const PetscInt         dim = 3;
 35:   ::tetgenio             in;
 36:   ::tetgenio             out;
 37:   PetscContainer         modelObj;
 38:   DMUniversalLabel       universal;
 39:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, defVal;
 40:   DMPlexInterpolatedFlag isInterpolated;
 41:   PetscMPIInt            rank;

 43:   PetscObjectGetComm((PetscObject)boundary, &comm);
 44:   MPI_Comm_rank(comm, &rank);
 45:   DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
 46:   DMUniversalLabelCreate(boundary, &universal);
 47:   DMLabelGetDefaultValue(universal->label, &defVal);

 49:   DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
 50:   in.numberofpoints = vEnd - vStart;
 51:   if (in.numberofpoints > 0) {
 52:     PetscSection       coordSection;
 53:     Vec                coordinates;
 54:     const PetscScalar *array;

 56:     in.pointlist       = new double[in.numberofpoints * dim];
 57:     in.pointmarkerlist = new int[in.numberofpoints];

 59:     PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints);
 60:     DMGetCoordinatesLocal(boundary, &coordinates);
 61:     DMGetCoordinateSection(boundary, &coordSection);
 62:     VecGetArrayRead(coordinates, &array);
 63:     for (v = vStart; v < vEnd; ++v) {
 64:       const PetscInt idx = v - vStart;
 65:       PetscInt       off, d, val;

 67:       PetscSectionGetOffset(coordSection, v, &off);
 68:       for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
 69:       DMLabelGetValue(universal->label, v, &val);
 70:       if (val != defVal) in.pointmarkerlist[idx] = (int)val;
 71:     }
 72:     VecRestoreArrayRead(coordinates, &array);
 73:   }

 75:   DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
 76:   in.numberofedges = eEnd - eStart;
 77:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
 78:     in.edgelist       = new int[in.numberofedges * 2];
 79:     in.edgemarkerlist = new int[in.numberofedges];
 80:     for (e = eStart; e < eEnd; ++e) {
 81:       const PetscInt  idx = e - eStart;
 82:       const PetscInt *cone;
 83:       PetscInt        coneSize, val;

 85:       DMPlexGetConeSize(boundary, e, &coneSize);
 86:       DMPlexGetCone(boundary, e, &cone);
 87:       in.edgelist[idx * 2]     = cone[0] - vStart;
 88:       in.edgelist[idx * 2 + 1] = cone[1] - vStart;

 90:       DMLabelGetValue(universal->label, e, &val);
 91:       if (val != defVal) in.edgemarkerlist[idx] = (int)val;
 92:     }
 93:   }

 95:   DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
 96:   in.numberoffacets = fEnd - fStart;
 97:   if (in.numberoffacets > 0) {
 98:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
 99:     in.facetmarkerlist = new int[in.numberoffacets];
100:     for (f = fStart; f < fEnd; ++f) {
101:       const PetscInt idx    = f - fStart;
102:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val = -1;

104:       in.facetlist[idx].numberofpolygons = 1;
105:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
106:       in.facetlist[idx].numberofholes    = 0;
107:       in.facetlist[idx].holelist         = NULL;

109:       DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
110:       for (p = 0; p < numPoints * 2; p += 2) {
111:         const PetscInt point = points[p];
112:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
113:       }

115:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
116:       poly->numberofvertices  = numVertices;
117:       poly->vertexlist        = new int[poly->numberofvertices];
118:       for (v = 0; v < numVertices; ++v) {
119:         const PetscInt vIdx = points[v] - vStart;
120:         poly->vertexlist[v] = vIdx;
121:       }
122:       DMLabelGetValue(universal->label, f, &val);
123:       if (val != defVal) in.facetmarkerlist[idx] = (int)val;
124:       DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
125:     }
126:   }
127:   if (rank == 0) {
128:     DM_Plex *mesh = (DM_Plex *)boundary->data;
129:     char     args[32];

131:     /* Take away 'Q' for verbose output */
132: #ifdef PETSC_HAVE_EGADS
133:     PetscStrcpy(args, "pqezQY");
134: #else
135:     PetscStrcpy(args, "pqezQ");
136: #endif
137:     if (mesh->tetgenOpts) {
138:       ::tetrahedralize(mesh->tetgenOpts, &in, &out);
139:     } else {
140:       ::tetrahedralize(args, &in, &out);
141:     }
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:       meshCoords = new PetscReal[dim * numVertices];
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:       cells = new PetscInt[numCells * numCorners];
164:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.tetrahedronlist[i];
165:     }

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

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

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

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

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

204:     PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj);
205:     if (modelObj) {
206: #ifdef PETSC_HAVE_EGADS
207:       DMLabel   bodyLabel;
208:       PetscInt  cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
209:       PetscBool islite = PETSC_FALSE;
210:       ego      *bodies;
211:       ego       model, geom;
212:       int       Nb, oclass, mtype, *senses;

214:       /* Get Attached EGADS Model from Original DMPlex */
215:       PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj);
216:       if (modelObj) {
217:         PetscContainerGetPointer(modelObj, (void **)&model);
218:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
219:         /* Transfer EGADS Model to Volumetric Mesh */
220:         PetscObjectCompose((PetscObject)*dm, "EGADS Model", (PetscObject)modelObj);
221:       } else {
222:         PetscObjectQuery((PetscObject)boundary, "EGADSLite Model", (PetscObject *)&modelObj);
223:         if (modelObj) {
224:           PetscContainerGetPointer(modelObj, (void **)&model);
225:           EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
226:           /* Transfer EGADS Model to Volumetric Mesh */
227:           PetscObjectCompose((PetscObject)*dm, "EGADSLite Model", (PetscObject)modelObj);
228:           islite = PETSC_TRUE;
229:         }
230:       }
231:       if (!modelObj) goto skip_egads;

233:       /* Set Cell Labels */
234:       DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
235:       DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
236:       DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
237:       DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);

239:       for (c = cStart; c < cEnd; ++c) {
240:         PetscReal centroid[3] = {0., 0., 0.};
241:         PetscInt  b;

243:         /* Determine what body the cell's centroid is located in */
244:         if (!interpolate) {
245:           PetscSection coordSection;
246:           Vec          coordinates;
247:           PetscScalar *coords = NULL;
248:           PetscInt     coordSize, s, d;

250:           DMGetCoordinatesLocal(*dm, &coordinates);
251:           DMGetCoordinateSection(*dm, &coordSection);
252:           DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
253:           for (s = 0; s < coordSize; ++s)
254:             for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
255:           DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
256:         } else DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
257:         for (b = 0; b < Nb; ++b) {
258:           if (islite) {
259:             if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
260:           } else {
261:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
262:           }
263:         }
264:         if (b < Nb) {
265:           PetscInt  cval    = b, eVal, fVal;
266:           PetscInt *closure = NULL, Ncl, cl;

268:           DMLabelSetValue(bodyLabel, c, cval);
269:           DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
270:           for (cl = 0; cl < Ncl; cl += 2) {
271:             const PetscInt p = closure[cl];

273:             if (p >= eStart && p < eEnd) {
274:               DMLabelGetValue(bodyLabel, p, &eVal);
275:               if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
276:             }
277:             if (p >= fStart && p < fEnd) {
278:               DMLabelGetValue(bodyLabel, p, &fVal);
279:               if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
280:             }
281:           }
282:           DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
283:         }
284:       }
285:     skip_egads:;
286: #endif
287:     }
288:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
289:   }
290:   DMUniversalLabelDestroy(&universal);
291:   return 0;
292: }

294: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
295: {
296:   MPI_Comm               comm;
297:   const PetscInt         dim = 3;
298:   ::tetgenio             in;
299:   ::tetgenio             out;
300:   PetscContainer         modelObj;
301:   DMUniversalLabel       universal;
302:   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal;
303:   DMPlexInterpolatedFlag isInterpolated;
304:   PetscMPIInt            rank;

306:   PetscObjectGetComm((PetscObject)dm, &comm);
307:   MPI_Comm_rank(comm, &rank);
308:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
309:   DMUniversalLabelCreate(dm, &universal);
310:   DMLabelGetDefaultValue(universal->label, &defVal);

312:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
313:   in.numberofpoints = vEnd - vStart;
314:   if (in.numberofpoints > 0) {
315:     PetscSection coordSection;
316:     Vec          coordinates;
317:     PetscScalar *array;

319:     in.pointlist       = new double[in.numberofpoints * dim];
320:     in.pointmarkerlist = new int[in.numberofpoints];

322:     PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints);
323:     DMGetCoordinatesLocal(dm, &coordinates);
324:     DMGetCoordinateSection(dm, &coordSection);
325:     VecGetArray(coordinates, &array);
326:     for (v = vStart; v < vEnd; ++v) {
327:       const PetscInt idx = v - vStart;
328:       PetscInt       off, d, val;

330:       PetscSectionGetOffset(coordSection, v, &off);
331:       for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
332:       DMLabelGetValue(universal->label, v, &val);
333:       if (val != defVal) in.pointmarkerlist[idx] = (int)val;
334:     }
335:     VecRestoreArray(coordinates, &array);
336:   }

338:   DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
339:   in.numberofedges = eEnd - eStart;
340:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
341:     in.edgelist       = new int[in.numberofedges * 2];
342:     in.edgemarkerlist = new int[in.numberofedges];
343:     for (e = eStart; e < eEnd; ++e) {
344:       const PetscInt  idx = e - eStart;
345:       const PetscInt *cone;
346:       PetscInt        coneSize, val;

348:       DMPlexGetConeSize(dm, e, &coneSize);
349:       DMPlexGetCone(dm, e, &cone);
350:       in.edgelist[idx * 2]     = cone[0] - vStart;
351:       in.edgelist[idx * 2 + 1] = cone[1] - vStart;

353:       DMLabelGetValue(universal->label, e, &val);
354:       if (val != defVal) in.edgemarkerlist[idx] = (int)val;
355:     }
356:   }

358:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
359:   in.numberoffacets = fEnd - fStart;
360:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
361:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
362:     in.facetmarkerlist = new int[in.numberoffacets];
363:     for (f = fStart; f < fEnd; ++f) {
364:       const PetscInt idx    = f - fStart;
365:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val;

367:       in.facetlist[idx].numberofpolygons = 1;
368:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
369:       in.facetlist[idx].numberofholes    = 0;
370:       in.facetlist[idx].holelist         = NULL;

372:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
373:       for (p = 0; p < numPoints * 2; p += 2) {
374:         const PetscInt point = points[p];
375:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
376:       }

378:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
379:       poly->numberofvertices  = numVertices;
380:       poly->vertexlist        = new int[poly->numberofvertices];
381:       for (v = 0; v < numVertices; ++v) {
382:         const PetscInt vIdx = points[v] - vStart;
383:         poly->vertexlist[v] = vIdx;
384:       }

386:       DMLabelGetValue(universal->label, f, &val);
387:       if (val != defVal) in.facetmarkerlist[idx] = (int)val;

389:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
390:     }
391:   }

393:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
394:   in.numberofcorners       = 4;
395:   in.numberoftetrahedra    = cEnd - cStart;
396:   in.tetrahedronvolumelist = (double *)maxVolumes;
397:   if (in.numberoftetrahedra > 0) {
398:     in.tetrahedronlist = new int[in.numberoftetrahedra * in.numberofcorners];
399:     for (c = cStart; c < cEnd; ++c) {
400:       const PetscInt idx     = c - cStart;
401:       PetscInt      *closure = NULL;
402:       PetscInt       closureSize;

404:       DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
406:       for (v = 0; v < 4; ++v) in.tetrahedronlist[idx * in.numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
407:       DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
408:     }
409:   }

411:   if (rank == 0) {
412:     char args[32];

414:     /* Take away 'Q' for verbose output */
415:     PetscStrcpy(args, "qezQra");
416:     ::tetrahedralize(args, &in, &out);
417:   }

419:   in.tetrahedronvolumelist = NULL;
420:   {
421:     const PetscInt numCorners  = 4;
422:     const PetscInt numCells    = out.numberoftetrahedra;
423:     const PetscInt numVertices = out.numberofpoints;
424:     PetscReal     *meshCoords  = NULL;
425:     PetscInt      *cells       = NULL;
426:     PetscBool      interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

428:     if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
429:       meshCoords = (PetscReal *)out.pointlist;
430:     } else {
431:       PetscInt i;

433:       meshCoords = new PetscReal[dim * numVertices];
434:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out.pointlist[i];
435:     }
436:     if (sizeof(PetscInt) == sizeof(out.tetrahedronlist[0])) {
437:       cells = (PetscInt *)out.tetrahedronlist;
438:     } else {
439:       PetscInt i;

441:       cells = new PetscInt[numCells * numCorners];
442:       for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out.tetrahedronlist[i];
443:     }

445:     DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
446:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
447:     if (sizeof(PetscReal) != sizeof(out.pointlist[0])) delete[] meshCoords;
448:     if (sizeof(PetscInt) != sizeof(out.tetrahedronlist[0])) delete[] cells;

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

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

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

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

484:     PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj);
485:     if (modelObj) {
486: #ifdef PETSC_HAVE_EGADS
487:       DMLabel   bodyLabel;
488:       PetscInt  cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
489:       PetscBool islite = PETSC_FALSE;
490:       ego      *bodies;
491:       ego       model, geom;
492:       int       Nb, oclass, mtype, *senses;

494:       /* Get Attached EGADS Model from Original DMPlex */
495:       PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj);
496:       if (modelObj) {
497:         PetscContainerGetPointer(modelObj, (void **)&model);
498:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
499:         /* Transfer EGADS Model to Volumetric Mesh */
500:         PetscObjectCompose((PetscObject)*dmRefined, "EGADS Model", (PetscObject)modelObj);
501:       } else {
502:         PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj);
503:         if (modelObj) {
504:           PetscContainerGetPointer(modelObj, (void **)&model);
505:           EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
506:           /* Transfer EGADS Model to Volumetric Mesh */
507:           PetscObjectCompose((PetscObject)*dmRefined, "EGADSLite Model", (PetscObject)modelObj);
508:           islite = PETSC_TRUE;
509:         }
510:       }
511:       if (!modelObj) goto skip_egads;

513:       /* Set Cell Labels */
514:       DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
515:       DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
516:       DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
517:       DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);

519:       for (c = cStart; c < cEnd; ++c) {
520:         PetscReal centroid[3] = {0., 0., 0.};
521:         PetscInt  b;

523:         /* Determine what body the cell's centroid is located in */
524:         if (!interpolate) {
525:           PetscSection coordSection;
526:           Vec          coordinates;
527:           PetscScalar *coords = NULL;
528:           PetscInt     coordSize, s, d;

530:           DMGetCoordinatesLocal(*dmRefined, &coordinates);
531:           DMGetCoordinateSection(*dmRefined, &coordSection);
532:           DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
533:           for (s = 0; s < coordSize; ++s)
534:             for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
535:           DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
536:         } else DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
537:         for (b = 0; b < Nb; ++b) {
538:           if (islite) {
539:             if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
540:           } else {
541:             if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
542:           }
543:         }
544:         if (b < Nb) {
545:           PetscInt  cval    = b, eVal, fVal;
546:           PetscInt *closure = NULL, Ncl, cl;

548:           DMLabelSetValue(bodyLabel, c, cval);
549:           DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
550:           for (cl = 0; cl < Ncl; cl += 2) {
551:             const PetscInt p = closure[cl];

553:             if (p >= eStart && p < eEnd) {
554:               DMLabelGetValue(bodyLabel, p, &eVal);
555:               if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
556:             }
557:             if (p >= fStart && p < fEnd) {
558:               DMLabelGetValue(bodyLabel, p, &fVal);
559:               if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
560:             }
561:           }
562:           DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
563:         }
564:       }
565:     skip_egads:;
566: #endif
567:     }
568:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
569:   }
570:   return 0;
571: }