Actual source code: ctetgenerate.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  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;

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

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

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

 39:   PLCCreate(&in);
 40:   PLCCreate(&out);

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

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

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

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

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

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

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

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

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

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

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

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

147:     if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
148:       meshCoords = (PetscReal *) out->pointlist;
149:     } else {
150:       PetscInt i;

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

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

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

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

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

189:           DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
190:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
191:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
192:           DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
193:         }
194:       }
195:       for (f = 0; f < out->numberoftrifaces; f++) {
196:         if (out->trifacemarkerlist[f]) {
197:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
198:           const PetscInt *faces;
199:           PetscInt        numFaces;

201:           DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
202:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
203:           DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
204:           DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
205:         }
206:       }
207:     }

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

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

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

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

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

243:             DMGetCoordinatesLocal(*dm, &coordinates);
244:             DMGetCoordinateSection(*dm, &coordSection);
245:             DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
246:             for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
247:             DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
248:           } else {
249:             DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
250:           }
251:           for (b = 0; b < Nb; ++b) {
252:             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:       }
276:     }
277: #endif
278:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
279:   }

281:   DMUniversalLabelDestroy(&universal);
282:   PLCDestroy(&in);
283:   PLCDestroy(&out);
284:   return(0);
285: }

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

299:   PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
300:   PetscObjectGetComm((PetscObject)dm,&comm);
301:   MPI_Comm_rank(comm, &rank);
302:   DMPlexIsInterpolatedCollective(dm, &isInterpolated);
303:   DMUniversalLabelCreate(dm, &universal);

305:   PLCCreate(&in);
306:   PLCCreate(&out);

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

315:     PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
316:     PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
317:     DMGetCoordinatesLocal(dm, &coordinates);
318:     DMGetCoordinateSection(dm, &coordSection);
319:     VecGetArray(coordinates, &array);
320:     for (v = vStart; v < vEnd; ++v) {
321:       const PetscInt idx = v - vStart;
322:       PetscInt       off, d, m;

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

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

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

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

352:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
353:   in->numberoftrifaces = 0;
354:   for (f = fStart; f < fEnd; ++f) {
355:     PetscInt supportSize;

357:     DMPlexGetSupportSize(dm, f, &supportSize);
358:     if (supportSize == 1) ++in->numberoftrifaces;
359:   }
360:   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
361:     PetscInt tf = 0;

363:     PetscMalloc1(in->numberoftrifaces*3, &in->trifacelist);
364:     PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist);
365:     for (f = fStart; f < fEnd; ++f) {
366:       PetscInt *points = NULL;
367:       PetscInt supportSize, numPoints, p, Nv = 0, val;

369:       DMPlexGetSupportSize(dm, f, &supportSize);
370:       if (supportSize != 1) continue;
371:       DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
372:       for (p = 0; p < numPoints*2; p += 2) {
373:         const PetscInt point = points[p];
374:         if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf*3 + Nv++] = point - vStart;
375:       }
376:       DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
377:       if (Nv != 3) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %D has %D vertices, not 3", f, Nv);
378:       DMLabelGetValue(universal->label, f, &val);
379:       in->trifacemarkerlist[tf] = (int) val;
380:       ++tf;
381:     }
382:   }

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

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

402:   if (!rank) {
403:     TetGenOpts t;

405:     TetGenOptsInitialize(&t);

407:     t.in        = dm; /* Should go away */
408:     t.refine    = 1;
409:     t.varvolume = 1;
410:     t.quality   = 1;
411:     t.edgesout  = 1;
412:     t.zeroindex = 1;
413:     t.quiet     = 1;
414:     t.verbose   = verbose; /* Change this */

416:     TetGenCheckOpts(&t);
417:     TetGenTetrahedralize(&t, in, out);
418:   }

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

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

434:       PetscMalloc1(dim * numVertices, &meshCoords);
435:       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
436:     }
437:     if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
438:       cells = (PetscInt *) out->tetrahedronlist;
439:     } else {
440:       PetscInt i;

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

446:     DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
447:     DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
448:     if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {PetscFree(meshCoords);}
449:     if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {PetscFree(cells);}

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

461:       for (e = 0; e < out->numberofedges; e++) {
462:         if (out->edgemarkerlist[e]) {
463:           const PetscInt  vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
464:           const PetscInt *edges;
465:           PetscInt        numEdges;

467:           DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
468:           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
469:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
470:           DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
471:         }
472:       }
473:       for (f = 0; f < out->numberoftrifaces; f++) {
474:         if (out->trifacemarkerlist[f]) {
475:           const PetscInt  vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
476:           const PetscInt *faces;
477:           PetscInt        numFaces;

479:           DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
480:           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
481:           DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
482:           DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
483:         }
484:       }
485:     }

487: #ifdef PETSC_HAVE_EGADS
488:     {
489:       DMLabel        bodyLabel;
490:       PetscContainer modelObj;
491:       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
492:       ego           *bodies;
493:       ego            model, geom;
494:       int            Nb, oclass, mtype, *senses;

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

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

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

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

521:             DMGetCoordinatesLocal(*dmRefined, &coordinates);
522:             DMGetCoordinateSection(*dmRefined, &coordSection);
523:             DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
524:             for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
525:             DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
526:           } else {
527:             DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
528:           }
529:           for (b = 0; b < Nb; ++b) {
530:             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:       }
554:     }
555: #endif
556:     DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
557:   }
558:   DMUniversalLabelDestroy(&universal);
559:   PLCDestroy(&in);
560:   PLCDestroy(&out);
561:   return(0);
562: }