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;

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

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

 41:   PetscObjectGetComm((PetscObject)boundary,&comm);
 42:   MPI_Comm_rank(comm, &rank);
 43:   DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
 44:   DMUniversalLabelCreate(boundary, &universal);

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

 53:     in.pointlist       = new double[in.numberofpoints*dim];
 54:     in.pointmarkerlist = new int[in.numberofpoints];

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

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

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

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

 86:       DMLabelGetValue(universal->label, e, &val);
 87:       in.edgemarkerlist[idx] = (int) val;
 88:     }
 89:   }

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

100:       in.facetlist[idx].numberofpolygons = 1;
101:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
102:       in.facetlist[idx].numberofholes    = 0;
103:       in.facetlist[idx].holelist         = NULL;

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

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

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

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

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

156:       cells = new PetscInt[numCells * numCorners];
157:       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out.tetrahedronlist[i];
158:     }

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

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

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

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

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

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

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

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

234:       for (c = cStart; c < cEnd; ++c) {
235:         PetscReal centroid[3] = {0., 0., 0.};
236:         PetscInt  b;

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

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

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

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

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

301:   PetscObjectGetComm((PetscObject)dm,&comm);
302:   MPI_Comm_rank(comm, &rank);
303:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
304:   DMUniversalLabelCreate(dm, &universal);

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

313:     in.pointlist       = new double[in.numberofpoints*dim];
314:     in.pointmarkerlist = new int[in.numberofpoints];

316:     DMGetCoordinatesLocal(dm, &coordinates);
317:     DMGetCoordinateSection(dm, &coordSection);
318:     VecGetArray(coordinates, &array);
319:     for (v = vStart; v < vEnd; ++v) {
320:       const PetscInt idx = v - vStart;
321:       PetscInt       off, d, val;

323:       PetscSectionGetOffset(coordSection, v, &off);
324:       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
325:       DMLabelGetValue(universal->label, v, &val);
326:       in.pointmarkerlist[idx] = (int) val;
327:     }
328:     VecRestoreArray(coordinates, &array);
329:   }

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

341:       DMPlexGetConeSize(dm, e, &coneSize);
342:       DMPlexGetCone(dm, e, &cone);
343:       in.edgelist[idx*2]     = cone[0] - vStart;
344:       in.edgelist[idx*2 + 1] = cone[1] - vStart;

346:       DMLabelGetValue(universal->label, e, &val);
347:       in.edgemarkerlist[idx] = (int) val;
348:     }
349:   }

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

360:       in.facetlist[idx].numberofpolygons = 1;
361:       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
362:       in.facetlist[idx].numberofholes    = 0;
363:       in.facetlist[idx].holelist         = NULL;

365:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
366:       for (p = 0; p < numPoints*2; p += 2) {
367:         const PetscInt point = points[p];
368:         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
369:       }

371:       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
372:       poly->numberofvertices = numVertices;
373:       poly->vertexlist       = new int[poly->numberofvertices];
374:       for (v = 0; v < numVertices; ++v) {
375:         const PetscInt vIdx = points[v] - vStart;
376:         poly->vertexlist[v] = vIdx;
377:       }

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

382:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
383:     }
384:   }

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

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

404:   if (rank == 0) {
405:     char args[32];

407:     /* Take away 'Q' for verbose output */
408:     PetscStrcpy(args, "qezQra");
409:     ::tetrahedralize(args, &in, &out);
410:   }

412:   in.tetrahedronvolumelist = NULL;
413:   {
414:     const PetscInt   numCorners  = 4;
415:     const PetscInt   numCells    = out.numberoftetrahedra;
416:     const PetscInt   numVertices = out.numberofpoints;
417:     PetscReal        *meshCoords = NULL;
418:     PetscInt         *cells      = NULL;
419:     PetscBool        interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;

421:     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
422:       meshCoords = (PetscReal *) out.pointlist;
423:     } else {
424:       PetscInt i;

426:       meshCoords = new PetscReal[dim * numVertices];
427:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
428:     }
429:     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
430:       cells = (PetscInt *) out.tetrahedronlist;
431:     } else {
432:       PetscInt i;

434:       cells = new PetscInt[numCells * numCorners];
435:       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
436:     }

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

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

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

459:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
460:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
461:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
462:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
463:         }
464:       }
465:       for (f = 0; f < out.numberoftrifaces; ++f) {
466:         if (out.trifacemarkerlist[f]) {
467:           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
468:           const PetscInt *faces;
469:           PetscInt        numFaces;

471:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
472:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
473:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
474:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
475:         }
476:       }
477:     }

479:     PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
480:     if (modelObj) {
481: #ifdef PETSC_HAVE_EGADS
482:       DMLabel        bodyLabel;
483:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
484:       PetscBool      islite = PETSC_FALSE;
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);
496:       } else {
497:         PetscObjectQuery((PetscObject) dm, "EGADSLite Model", (PetscObject *) &modelObj);
498:         if (modelObj) {
499:           PetscContainerGetPointer(modelObj, (void **) &model);
500:           EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
501:           /* Transfer EGADS Model to Volumetric Mesh */
502:           PetscObjectCompose((PetscObject) *dmRefined, "EGADSLite Model", (PetscObject) modelObj);
503:           islite = PETSC_TRUE;
504:         }
505:       }
506:       if (!modelObj) goto skip_egads;

508:       /* Set Cell Labels */
509:       DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
510:       DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
511:       DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
512:       DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);

514:       for (c = cStart; c < cEnd; ++c) {
515:         PetscReal centroid[3] = {0., 0., 0.};
516:         PetscInt  b;

518:         /* Deterimine what body the cell's centroid is located in */
519:         if (!interpolate) {
520:           PetscSection   coordSection;
521:           Vec            coordinates;
522:           PetscScalar   *coords = NULL;
523:           PetscInt       coordSize, s, d;

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

541:           DMLabelSetValue(bodyLabel, c, cval);
542:           DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
543:           for (cl = 0; cl < Ncl; cl += 2) {
544:             const PetscInt p = closure[cl];

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