xref: /petsc/src/dm/impls/plex/generators/tetgen/tetgenerate.cxx (revision 0db4d2e0165a0ea245cec0549c2de0bb7b39e2c0) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
3 #ifdef PETSC_HAVE_EGADS
4 #include <egads.h>
5 #endif
6 
7 #if defined(PETSC_HAVE_TETGEN_TETLIBRARY_NEEDED)
8 #define TETLIBRARY
9 #endif
10 #include <tetgen.h>
11 
12 /* This is to fix the tetrahedron orientation from TetGen */
13 static PetscErrorCode DMPlexInvertCells_Tetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
14 {
15   PetscInt bound = numCells*numCorners, coff;
16 
17   PetscFunctionBegin;
18 #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
19   for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
20 #undef SWAP
21   PetscFunctionReturn(0);
22 }
23 
24 PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
25 {
26   MPI_Comm               comm;
27   const PetscInt         dim = 3;
28   ::tetgenio             in;
29   ::tetgenio             out;
30   DMUniversalLabel       universal;
31   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f;
32   DMPlexInterpolatedFlag isInterpolated;
33   PetscMPIInt            rank;
34   PetscErrorCode         ierr;
35 
36   PetscFunctionBegin;
37   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
38   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
39   ierr = DMPlexIsInterpolatedCollective(boundary, &isInterpolated);CHKERRQ(ierr);
40   ierr = DMUniversalLabelCreate(boundary, &universal);CHKERRQ(ierr);
41 
42   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
43   in.numberofpoints = vEnd - vStart;
44   if (in.numberofpoints > 0) {
45     PetscSection       coordSection;
46     Vec                coordinates;
47     const PetscScalar *array;
48 
49     in.pointlist       = new double[in.numberofpoints*dim];
50     in.pointmarkerlist = new int[in.numberofpoints];
51 
52     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
53     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
54     ierr = VecGetArrayRead(coordinates, &array);CHKERRQ(ierr);
55     for (v = vStart; v < vEnd; ++v) {
56       const PetscInt idx = v - vStart;
57       PetscInt       off, d, val;
58 
59       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
60       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
61       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
62       in.pointmarkerlist[idx] = (int) val;
63     }
64     ierr = VecRestoreArrayRead(coordinates, &array);CHKERRQ(ierr);
65   }
66 
67   ierr = DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);CHKERRQ(ierr);
68   in.numberofedges = eEnd - eStart;
69   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
70     in.edgelist       = new int[in.numberofedges * 2];
71     in.edgemarkerlist = new int[in.numberofedges];
72     for (e = eStart; e < eEnd; ++e) {
73       const PetscInt  idx = e - eStart;
74       const PetscInt *cone;
75       PetscInt        coneSize, val;
76 
77       ierr = DMPlexGetConeSize(boundary, e, &coneSize);CHKERRQ(ierr);
78       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
79       in.edgelist[idx*2]     = cone[0] - vStart;
80       in.edgelist[idx*2 + 1] = cone[1] - vStart;
81 
82       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
83       in.edgemarkerlist[idx] = (int) val;
84     }
85   }
86 
87   ierr = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
88   in.numberoffacets = fEnd - fStart;
89   if (in.numberoffacets > 0) {
90     in.facetlist       = new tetgenio::facet[in.numberoffacets];
91     in.facetmarkerlist = new int[in.numberoffacets];
92     for (f = fStart; f < fEnd; ++f) {
93       const PetscInt idx    = f - fStart;
94       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val = -1;
95 
96       in.facetlist[idx].numberofpolygons = 1;
97       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
98       in.facetlist[idx].numberofholes    = 0;
99       in.facetlist[idx].holelist         = NULL;
100 
101       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
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       }
106 
107       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
108       poly->numberofvertices = numVertices;
109       poly->vertexlist       = new int[poly->numberofvertices];
110       for (v = 0; v < numVertices; ++v) {
111         const PetscInt vIdx = points[v] - vStart;
112         poly->vertexlist[v] = vIdx;
113       }
114       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
115       in.facetmarkerlist[idx] = (int) val;
116       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
117     }
118   }
119   if (!rank) {
120     DM_Plex *mesh = (DM_Plex *) boundary->data;
121     char     args[32];
122 
123     /* Take away 'Q' for verbose output */
124 #ifdef PETSC_HAVE_EGADS
125     ierr = PetscStrcpy(args, "pqezQY");CHKERRQ(ierr);
126 #else
127     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
128 #endif
129     if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);}
130     else                  {::tetrahedralize(args, &in, &out);}
131   }
132   {
133     const PetscInt   numCorners  = 4;
134     const PetscInt   numCells    = out.numberoftetrahedra;
135     const PetscInt   numVertices = out.numberofpoints;
136     PetscReal        *meshCoords = NULL;
137     PetscInt         *cells      = NULL;
138 
139     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
140       meshCoords = (PetscReal *) out.pointlist;
141     } else {
142       PetscInt i;
143 
144       meshCoords = new PetscReal[dim * numVertices];
145       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
146     }
147     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
148       cells = (PetscInt *) out.tetrahedronlist;
149     } else {
150       PetscInt i;
151 
152       cells = new PetscInt[numCells * numCorners];
153       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out.tetrahedronlist[i];
154     }
155 
156     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
157     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
158 
159     /* Set labels */
160     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);CHKERRQ(ierr);
161     for (v = 0; v < numVertices; ++v) {
162       if (out.pointmarkerlist[v]) {
163         ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
164       }
165     }
166     if (interpolate) {
167       PetscInt e;
168 
169       for (e = 0; e < out.numberofedges; e++) {
170         if (out.edgemarkerlist[e]) {
171           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
172           const PetscInt *edges;
173           PetscInt        numEdges;
174 
175           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
176           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
177           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
178           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
179         }
180       }
181       for (f = 0; f < out.numberoftrifaces; f++) {
182         if (out.trifacemarkerlist[f]) {
183           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
184           const PetscInt *faces;
185           PetscInt        numFaces;
186 
187           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
188           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
189           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
190           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
191         }
192       }
193     }
194 
195 #ifdef PETSC_HAVE_EGADS
196     {
197       DMLabel        bodyLabel;
198       PetscContainer modelObj;
199       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
200       ego           *bodies;
201       ego            model, geom;
202       int            Nb, oclass, mtype, *senses;
203 
204       /* Get Attached EGADS Model from Original DMPlex */
205       ierr = PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
206       ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
207       ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
208       /* Transfer EGADS Model to Volumetric Mesh */
209       ierr = PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr);
210 
211       /* Set Cell Labels */
212       ierr = DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr);
213       ierr = DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
214       ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
215       ierr = DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
216 
217       for (c = cStart; c < cEnd; ++c) {
218         PetscReal centroid[3] = {0., 0., 0.};
219         PetscInt  b;
220 
221         /* Deterimine what body the cell's centroid is located in */
222         if (!interpolate) {
223           PetscSection   coordSection;
224           Vec            coordinates;
225           PetscScalar   *coords = NULL;
226           PetscInt       coordSize, s, d;
227 
228           ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr);
229           ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
230           ierr = DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
231           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
232           ierr = DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
233         } else {
234           ierr = DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);CHKERRQ(ierr);
235         }
236         for (b = 0; b < Nb; ++b) {
237           if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
238         }
239         if (b < Nb) {
240           PetscInt   cval = b, eVal, fVal;
241           PetscInt *closure = NULL, Ncl, cl;
242 
243           ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr);
244           ierr = DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
245           for (cl = 0; cl < Ncl; cl += 2) {
246             const PetscInt p = closure[cl];
247 
248             if (p >= eStart && p < eEnd) {
249               ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr);
250               if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
251             }
252             if (p >= fStart && p < fEnd) {
253               ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr);
254               if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
255             }
256           }
257           ierr = DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
258         }
259       }
260     }
261 #endif
262     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
263   }
264   PetscFunctionReturn(0);
265 }
266 
267 PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
268 {
269   MPI_Comm               comm;
270   const PetscInt         dim = 3;
271   ::tetgenio             in;
272   ::tetgenio             out;
273   DMUniversalLabel       universal;
274   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c;
275   DMPlexInterpolatedFlag isInterpolated;
276   PetscMPIInt            rank;
277   PetscErrorCode         ierr;
278 
279   PetscFunctionBegin;
280   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
281   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
282   ierr = DMPlexIsInterpolatedCollective(dm, &isInterpolated);CHKERRQ(ierr);
283   ierr = DMUniversalLabelCreate(dm, &universal);CHKERRQ(ierr);
284 
285   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
286   in.numberofpoints = vEnd - vStart;
287   if (in.numberofpoints > 0) {
288     PetscSection coordSection;
289     Vec          coordinates;
290     PetscScalar *array;
291 
292     in.pointlist       = new double[in.numberofpoints*dim];
293     in.pointmarkerlist = new int[in.numberofpoints];
294 
295     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
296     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
297     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
298     for (v = vStart; v < vEnd; ++v) {
299       const PetscInt idx = v - vStart;
300       PetscInt       off, d, val;
301 
302       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
303       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
304       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
305       in.pointmarkerlist[idx] = (int) val;
306     }
307     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
308   }
309 
310   ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
311   in.numberofedges = eEnd - eStart;
312   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
313     in.edgelist       = new int[in.numberofedges * 2];
314     in.edgemarkerlist = new int[in.numberofedges];
315     for (e = eStart; e < eEnd; ++e) {
316       const PetscInt  idx = e - eStart;
317       const PetscInt *cone;
318       PetscInt        coneSize, val;
319 
320       ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
321       ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
322       in.edgelist[idx*2]     = cone[0] - vStart;
323       in.edgelist[idx*2 + 1] = cone[1] - vStart;
324 
325       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
326       in.edgemarkerlist[idx] = (int) val;
327     }
328   }
329 
330   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
331   in.numberoffacets = fEnd - fStart;
332   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
333     in.facetlist       = new tetgenio::facet[in.numberoffacets];
334     in.facetmarkerlist = new int[in.numberoffacets];
335     for (f = fStart; f < fEnd; ++f) {
336       const PetscInt idx    = f - fStart;
337       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val;
338 
339       in.facetlist[idx].numberofpolygons = 1;
340       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
341       in.facetlist[idx].numberofholes    = 0;
342       in.facetlist[idx].holelist         = NULL;
343 
344       ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
345       for (p = 0; p < numPoints*2; p += 2) {
346         const PetscInt point = points[p];
347         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
348       }
349 
350       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
351       poly->numberofvertices = numVertices;
352       poly->vertexlist       = new int[poly->numberofvertices];
353       for (v = 0; v < numVertices; ++v) {
354         const PetscInt vIdx = points[v] - vStart;
355         poly->vertexlist[v] = vIdx;
356       }
357 
358       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
359       in.facetmarkerlist[idx] = (int) val;
360 
361       ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
362     }
363   }
364 
365   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
366   in.numberofcorners       = 4;
367   in.numberoftetrahedra    = cEnd - cStart;
368   in.tetrahedronvolumelist = (double *) maxVolumes;
369   if (in.numberoftetrahedra > 0) {
370     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
371     for (c = cStart; c < cEnd; ++c) {
372       const PetscInt idx     = c - cStart;
373       PetscInt      *closure = NULL;
374       PetscInt       closureSize;
375 
376       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
377       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
378       for (v = 0; v < 4; ++v) in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
379       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
380     }
381   }
382 
383   if (!rank) {
384     char args[32];
385 
386     /* Take away 'Q' for verbose output */
387     ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr);
388     ::tetrahedralize(args, &in, &out);
389   }
390 
391   in.tetrahedronvolumelist = NULL;
392   {
393     const PetscInt   numCorners  = 4;
394     const PetscInt   numCells    = out.numberoftetrahedra;
395     const PetscInt   numVertices = out.numberofpoints;
396     PetscReal        *meshCoords = NULL;
397     PetscInt         *cells      = NULL;
398     PetscBool        interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
399 
400     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
401       meshCoords = (PetscReal *) out.pointlist;
402     } else {
403       PetscInt i;
404 
405       meshCoords = new PetscReal[dim * numVertices];
406       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
407     }
408     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
409       cells = (PetscInt *) out.tetrahedronlist;
410     } else {
411       PetscInt i;
412 
413       cells = new PetscInt[numCells * numCorners];
414       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
415     }
416 
417     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
418     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
419     if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {delete [] meshCoords;}
420     if (sizeof (PetscInt) != sizeof (out.tetrahedronlist[0])) {delete [] cells;}
421 
422     /* Set labels */
423     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);CHKERRQ(ierr);
424     for (v = 0; v < numVertices; ++v) {
425       if (out.pointmarkerlist[v]) {
426         ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
427       }
428     }
429     if (interpolate) {
430       PetscInt e, f;
431 
432       for (e = 0; e < out.numberofedges; ++e) {
433         if (out.edgemarkerlist[e]) {
434           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
435           const PetscInt *edges;
436           PetscInt        numEdges;
437 
438           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
439           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
440           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
441           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
442         }
443       }
444       for (f = 0; f < out.numberoftrifaces; ++f) {
445         if (out.trifacemarkerlist[f]) {
446           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
447           const PetscInt *faces;
448           PetscInt        numFaces;
449 
450           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
451           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
452           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
453           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
454         }
455       }
456     }
457 
458 #ifdef PETSC_HAVE_EGADS
459     {
460       DMLabel        bodyLabel;
461       PetscContainer modelObj;
462       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
463       ego           *bodies;
464       ego            model, geom;
465       int            Nb, oclass, mtype, *senses;
466 
467       /* Get Attached EGADS Model from Original DMPlex */
468       ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
469       ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
470       ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
471       /* Transfer EGADS Model to Volumetric Mesh */
472       ierr = PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr);
473 
474       /* Set Cell Labels */
475       ierr = DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr);
476       ierr = DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);CHKERRQ(ierr);
477       ierr = DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);CHKERRQ(ierr);
478       ierr = DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);CHKERRQ(ierr);
479 
480       for (c = cStart; c < cEnd; ++c) {
481         PetscReal centroid[3] = {0., 0., 0.};
482         PetscInt  b;
483 
484         /* Deterimine what body the cell's centroid is located in */
485         if (!interpolate) {
486           PetscSection   coordSection;
487           Vec            coordinates;
488           PetscScalar   *coords = NULL;
489           PetscInt       coordSize, s, d;
490 
491           ierr = DMGetCoordinatesLocal(*dmRefined, &coordinates);CHKERRQ(ierr);
492           ierr = DMGetCoordinateSection(*dmRefined, &coordSection);CHKERRQ(ierr);
493           ierr = DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
494           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
495           ierr = DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
496         } else {
497           ierr = DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);CHKERRQ(ierr);
498         }
499         for (b = 0; b < Nb; ++b) {
500           if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
501         }
502         if (b < Nb) {
503           PetscInt   cval = b, eVal, fVal;
504           PetscInt *closure = NULL, Ncl, cl;
505 
506           ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr);
507           ierr = DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
508           for (cl = 0; cl < Ncl; cl += 2) {
509             const PetscInt p = closure[cl];
510 
511             if (p >= eStart && p < eEnd) {
512               ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr);
513               if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
514             }
515             if (p >= fStart && p < fEnd) {
516               ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr);
517               if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
518             }
519           }
520           ierr = DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
521         }
522       }
523     }
524 #endif
525     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
526   }
527   PetscFunctionReturn(0);
528 }
529