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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

176:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
178:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
179:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
180:         }
181:       }
182:       for (f = 0; f < out.numberoftrifaces; f++) {
183:         if (out.trifacemarkerlist[f]) {
184:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
185:           const PetscInt *faces;
186:           PetscInt        numFaces;

188:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
190:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
191:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
192:         }
193:       }
194:     }

196:     PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);
197:     if (modelObj) {
198: #ifdef PETSC_HAVE_EGADS
199:       DMLabel        bodyLabel;
200:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
201:       PetscBool      islite = PETSC_FALSE;
202:       ego           *bodies;
203:       ego            model, geom;
204:       int            Nb, oclass, mtype, *senses;

206:       /* Get Attached EGADS Model from Original DMPlex */
207:       PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);
208:       if (modelObj) {
209:         PetscContainerGetPointer(modelObj, (void **) &model);
210:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
211:         /* Transfer EGADS Model to Volumetric Mesh */
212:         PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);
213:       } else {
214:         PetscObjectQuery((PetscObject) boundary, "EGADSLite Model", (PetscObject *) &modelObj);
215:         if (modelObj) {
216:           PetscContainerGetPointer(modelObj, (void **) &model);
217:           EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
218:           /* Transfer EGADS Model to Volumetric Mesh */
219:           PetscObjectCompose((PetscObject) *dm, "EGADSLite Model", (PetscObject) modelObj);
220:           islite = PETSC_TRUE;
221:         }
222:       }
223:       if (!modelObj) goto skip_egads;

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

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

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

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

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

263:             if (p >= eStart && p < eEnd) {
264:               DMLabelGetValue(bodyLabel, p, &eVal);
265:               if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
266:             }
267:             if (p >= fStart && p < fEnd) {
268:               DMLabelGetValue(bodyLabel, p, &fVal);
269:               if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
270:             }
271:           }
272:           DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
273:         }
274:       }
275: skip_egads: ;
276: #endif
277:     }
278:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
279:   }
280:   DMUniversalLabelDestroy(&universal);
281:   return 0;
282: }

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

296:   PetscObjectGetComm((PetscObject)dm,&comm);
297:   MPI_Comm_rank(comm, &rank);
298:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
299:   DMUniversalLabelCreate(dm, &universal);

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

308:     in.pointlist       = new double[in.numberofpoints*dim];
309:     in.pointmarkerlist = new int[in.numberofpoints];

311:     DMGetCoordinatesLocal(dm, &coordinates);
312:     DMGetCoordinateSection(dm, &coordSection);
313:     VecGetArray(coordinates, &array);
314:     for (v = vStart; v < vEnd; ++v) {
315:       const PetscInt idx = v - vStart;
316:       PetscInt       off, d, val;

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

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

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

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

346:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
347:   in.numberoffacets = fEnd - fStart;
348:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
349:     in.facetlist       = new tetgenio::facet[in.numberoffacets];
350:     in.facetmarkerlist = new int[in.numberoffacets];
351:     for (f = fStart; f < fEnd; ++f) {
352:       const PetscInt idx    = f - fStart;
353:       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val;

355:       in.facetlist[idx].numberofpolygons = 1;
356:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
357:       in.facetlist[idx].numberofholes    = 0;
358:       in.facetlist[idx].holelist         = NULL;

360:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
361:       for (p = 0; p < numPoints*2; p += 2) {
362:         const PetscInt point = points[p];
363:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
364:       }

366:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
367:       poly->numberofvertices = numVertices;
368:       poly->vertexlist       = new int[poly->numberofvertices];
369:       for (v = 0; v < numVertices; ++v) {
370:         const PetscInt vIdx = points[v] - vStart;
371:         poly->vertexlist[v] = vIdx;
372:       }

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

377:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
378:     }
379:   }

381:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
382:   in.numberofcorners       = 4;
383:   in.numberoftetrahedra    = cEnd - cStart;
384:   in.tetrahedronvolumelist = (double *) maxVolumes;
385:   if (in.numberoftetrahedra > 0) {
386:     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
387:     for (c = cStart; c < cEnd; ++c) {
388:       const PetscInt idx     = c - cStart;
389:       PetscInt      *closure = NULL;
390:       PetscInt       closureSize;

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

399:   if (rank == 0) {
400:     char args[32];

402:     /* Take away 'Q' for verbose output */
403:     PetscStrcpy(args, "qezQra");
404:     ::tetrahedralize(args, &in, &out);
405:   }

407:   in.tetrahedronvolumelist = NULL;
408:   {
409:     const PetscInt   numCorners  = 4;
410:     const PetscInt   numCells    = out.numberoftetrahedra;
411:     const PetscInt   numVertices = out.numberofpoints;
412:     PetscReal        *meshCoords = NULL;
413:     PetscInt         *cells      = NULL;
414:     PetscBool        interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

416:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
417:       meshCoords = (PetscReal *) out.pointlist;
418:     } else {
419:       PetscInt i;

421:       meshCoords = new PetscReal[dim * numVertices];
422:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
423:     }
424:     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
425:       cells = (PetscInt *) out.tetrahedronlist;
426:     } else {
427:       PetscInt i;

429:       cells = new PetscInt[numCells * numCorners];
430:       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
431:     }

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

438:     /* Set labels */
439:     DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
440:     for (v = 0; v < numVertices; ++v) {
441:       if (out.pointmarkerlist[v]) {
442:         DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);
443:       }
444:     }
445:     if (interpolate) {
446:       PetscInt e, f;

448:       for (e = 0; e < out.numberofedges; ++e) {
449:         if (out.edgemarkerlist[e]) {
450:           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
451:           const PetscInt *edges;
452:           PetscInt        numEdges;

454:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
456:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
457:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
458:         }
459:       }
460:       for (f = 0; f < out.numberoftrifaces; ++f) {
461:         if (out.trifacemarkerlist[f]) {
462:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
463:           const PetscInt *faces;
464:           PetscInt        numFaces;

466:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
468:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
469:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
470:         }
471:       }
472:     }

474:     PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
475:     if (modelObj) {
476: #ifdef PETSC_HAVE_EGADS
477:       DMLabel        bodyLabel;
478:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
479:       PetscBool      islite = PETSC_FALSE;
480:       ego           *bodies;
481:       ego            model, geom;
482:       int            Nb, oclass, mtype, *senses;

484:       /* Get Attached EGADS Model from Original DMPlex */
485:       PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
486:       if (modelObj) {
487:         PetscContainerGetPointer(modelObj, (void **) &model);
488:         EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
489:         /* Transfer EGADS Model to Volumetric Mesh */
490:         PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);
491:       } else {
492:         PetscObjectQuery((PetscObject) dm, "EGADSLite Model", (PetscObject *) &modelObj);
493:         if (modelObj) {
494:           PetscContainerGetPointer(modelObj, (void **) &model);
495:           EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
496:           /* Transfer EGADS Model to Volumetric Mesh */
497:           PetscObjectCompose((PetscObject) *dmRefined, "EGADSLite Model", (PetscObject) modelObj);
498:           islite = PETSC_TRUE;
499:         }
500:       }
501:       if (!modelObj) goto skip_egads;

503:       /* Set Cell Labels */
504:       DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
505:       DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
506:       DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
507:       DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);

509:       for (c = cStart; c < cEnd; ++c) {
510:         PetscReal centroid[3] = {0., 0., 0.};
511:         PetscInt  b;

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

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

536:           DMLabelSetValue(bodyLabel, c, cval);
537:           DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
538:           for (cl = 0; cl < Ncl; cl += 2) {
539:             const PetscInt p = closure[cl];

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