xref: /petsc/src/dm/impls/plex/generators/tetgen/tetgenerate.cxx (revision 362febeeeb69b91ebadcb4b2dc0a22cb6dfc4097) 
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   PetscContainer         modelObj;
31   DMUniversalLabel       universal;
32   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f;
33   DMPlexInterpolatedFlag isInterpolated;
34   PetscMPIInt            rank;
35   PetscErrorCode         ierr;
36 
37   PetscFunctionBegin;
38   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
39   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
40   ierr = DMPlexIsInterpolatedCollective(boundary, &isInterpolated);CHKERRQ(ierr);
41   ierr = DMUniversalLabelCreate(boundary, &universal);CHKERRQ(ierr);
42 
43   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
44   in.numberofpoints = vEnd - vStart;
45   if (in.numberofpoints > 0) {
46     PetscSection       coordSection;
47     Vec                coordinates;
48     const PetscScalar *array;
49 
50     in.pointlist       = new double[in.numberofpoints*dim];
51     in.pointmarkerlist = new int[in.numberofpoints];
52 
53     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
54     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
55     ierr = VecGetArrayRead(coordinates, &array);CHKERRQ(ierr);
56     for (v = vStart; v < vEnd; ++v) {
57       const PetscInt idx = v - vStart;
58       PetscInt       off, d, val;
59 
60       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
61       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
62       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
63       in.pointmarkerlist[idx] = (int) val;
64     }
65     ierr = VecRestoreArrayRead(coordinates, &array);CHKERRQ(ierr);
66   }
67 
68   ierr = DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);CHKERRQ(ierr);
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;
77 
78       ierr = DMPlexGetConeSize(boundary, e, &coneSize);CHKERRQ(ierr);
79       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
80       in.edgelist[idx*2]     = cone[0] - vStart;
81       in.edgelist[idx*2 + 1] = cone[1] - vStart;
82 
83       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
84       in.edgemarkerlist[idx] = (int) val;
85     }
86   }
87 
88   ierr = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
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;
96 
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;
101 
102       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
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       }
107 
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       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
116       in.facetmarkerlist[idx] = (int) val;
117       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
118     }
119   }
120   if (!rank) {
121     DM_Plex *mesh = (DM_Plex *) boundary->data;
122     char     args[32];
123 
124     /* Take away 'Q' for verbose output */
125 #ifdef PETSC_HAVE_EGADS
126     ierr = PetscStrcpy(args, "pqezQY");CHKERRQ(ierr);
127 #else
128     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
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;
139 
140     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
141       meshCoords = (PetscReal *) out.pointlist;
142     } else {
143       PetscInt i;
144 
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;
152 
153       cells = new PetscInt[numCells * numCorners];
154       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out.tetrahedronlist[i];
155     }
156 
157     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
158     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
159 
160     /* Set labels */
161     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);CHKERRQ(ierr);
162     for (v = 0; v < numVertices; ++v) {
163       if (out.pointmarkerlist[v]) {
164         ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
165       }
166     }
167     if (interpolate) {
168       PetscInt e;
169 
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;
175 
176           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
177           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
178           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
179           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
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;
187 
188           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
189           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
190           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
191           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
192         }
193       }
194     }
195 
196     ierr = PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
197     if (modelObj) {
198 #ifdef PETSC_HAVE_EGADS
199       DMLabel        bodyLabel;
200       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
201       ego           *bodies;
202       ego            model, geom;
203       int            Nb, oclass, mtype, *senses;
204 
205       /* Get Attached EGADS Model from Original DMPlex */
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 #endif
261     }
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   PetscContainer         modelObj;
274   DMUniversalLabel       universal;
275   PetscInt               vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c;
276   DMPlexInterpolatedFlag isInterpolated;
277   PetscMPIInt            rank;
278   PetscErrorCode         ierr;
279 
280   PetscFunctionBegin;
281   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
282   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
283   ierr = DMPlexIsInterpolatedCollective(dm, &isInterpolated);CHKERRQ(ierr);
284   ierr = DMUniversalLabelCreate(dm, &universal);CHKERRQ(ierr);
285 
286   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
287   in.numberofpoints = vEnd - vStart;
288   if (in.numberofpoints > 0) {
289     PetscSection coordSection;
290     Vec          coordinates;
291     PetscScalar *array;
292 
293     in.pointlist       = new double[in.numberofpoints*dim];
294     in.pointmarkerlist = new int[in.numberofpoints];
295 
296     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
297     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
298     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
299     for (v = vStart; v < vEnd; ++v) {
300       const PetscInt idx = v - vStart;
301       PetscInt       off, d, val;
302 
303       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
304       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
305       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
306       in.pointmarkerlist[idx] = (int) val;
307     }
308     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
309   }
310 
311   ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
312   in.numberofedges = eEnd - eStart;
313   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
314     in.edgelist       = new int[in.numberofedges * 2];
315     in.edgemarkerlist = new int[in.numberofedges];
316     for (e = eStart; e < eEnd; ++e) {
317       const PetscInt  idx = e - eStart;
318       const PetscInt *cone;
319       PetscInt        coneSize, val;
320 
321       ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
322       ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
323       in.edgelist[idx*2]     = cone[0] - vStart;
324       in.edgelist[idx*2 + 1] = cone[1] - vStart;
325 
326       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
327       in.edgemarkerlist[idx] = (int) val;
328     }
329   }
330 
331   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
332   in.numberoffacets = fEnd - fStart;
333   if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
334     in.facetlist       = new tetgenio::facet[in.numberoffacets];
335     in.facetmarkerlist = new int[in.numberoffacets];
336     for (f = fStart; f < fEnd; ++f) {
337       const PetscInt idx    = f - fStart;
338       PetscInt      *points = NULL, numPoints, p, numVertices = 0, v, val;
339 
340       in.facetlist[idx].numberofpolygons = 1;
341       in.facetlist[idx].polygonlist      = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
342       in.facetlist[idx].numberofholes    = 0;
343       in.facetlist[idx].holelist         = NULL;
344 
345       ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
346       for (p = 0; p < numPoints*2; p += 2) {
347         const PetscInt point = points[p];
348         if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
349       }
350 
351       tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
352       poly->numberofvertices = numVertices;
353       poly->vertexlist       = new int[poly->numberofvertices];
354       for (v = 0; v < numVertices; ++v) {
355         const PetscInt vIdx = points[v] - vStart;
356         poly->vertexlist[v] = vIdx;
357       }
358 
359       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
360       in.facetmarkerlist[idx] = (int) val;
361 
362       ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
363     }
364   }
365 
366   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
367   in.numberofcorners       = 4;
368   in.numberoftetrahedra    = cEnd - cStart;
369   in.tetrahedronvolumelist = (double *) maxVolumes;
370   if (in.numberoftetrahedra > 0) {
371     in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners];
372     for (c = cStart; c < cEnd; ++c) {
373       const PetscInt idx     = c - cStart;
374       PetscInt      *closure = NULL;
375       PetscInt       closureSize;
376 
377       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
378       if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
379       for (v = 0; v < 4; ++v) in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
380       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
381     }
382   }
383 
384   if (!rank) {
385     char args[32];
386 
387     /* Take away 'Q' for verbose output */
388     ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr);
389     ::tetrahedralize(args, &in, &out);
390   }
391 
392   in.tetrahedronvolumelist = NULL;
393   {
394     const PetscInt   numCorners  = 4;
395     const PetscInt   numCells    = out.numberoftetrahedra;
396     const PetscInt   numVertices = out.numberofpoints;
397     PetscReal        *meshCoords = NULL;
398     PetscInt         *cells      = NULL;
399     PetscBool        interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
400 
401     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
402       meshCoords = (PetscReal *) out.pointlist;
403     } else {
404       PetscInt i;
405 
406       meshCoords = new PetscReal[dim * numVertices];
407       for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out.pointlist[i];
408     }
409     if (sizeof (PetscInt) == sizeof (out.tetrahedronlist[0])) {
410       cells = (PetscInt *) out.tetrahedronlist;
411     } else {
412       PetscInt i;
413 
414       cells = new PetscInt[numCells * numCorners];
415       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
416     }
417 
418     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
419     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
420     if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {delete [] meshCoords;}
421     if (sizeof (PetscInt) != sizeof (out.tetrahedronlist[0])) {delete [] cells;}
422 
423     /* Set labels */
424     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);CHKERRQ(ierr);
425     for (v = 0; v < numVertices; ++v) {
426       if (out.pointmarkerlist[v]) {
427         ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
428       }
429     }
430     if (interpolate) {
431       PetscInt e, f;
432 
433       for (e = 0; e < out.numberofedges; ++e) {
434         if (out.edgemarkerlist[e]) {
435           const PetscInt  vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells};
436           const PetscInt *edges;
437           PetscInt        numEdges;
438 
439           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
440           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
441           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
442           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
443         }
444       }
445       for (f = 0; f < out.numberoftrifaces; ++f) {
446         if (out.trifacemarkerlist[f]) {
447           const PetscInt  vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells};
448           const PetscInt *faces;
449           PetscInt        numFaces;
450 
451           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
452           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
453           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
454           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
455         }
456       }
457     }
458 
459     ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
460     if (modelObj) {
461 #ifdef PETSC_HAVE_EGADS
462       DMLabel        bodyLabel;
463       PetscInt       cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
464       ego           *bodies;
465       ego            model, geom;
466       int            Nb, oclass, mtype, *senses;
467 
468       /* Get Attached EGADS Model from Original DMPlex */
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 #endif
524     }
525     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
526   }
527   PetscFunctionReturn(0);
528 }
529