xref: /petsc/src/dm/impls/plex/generators/tetgen/tetgenerate.cxx (revision ccb4e88a40f0b86eaeca07ff64c64e4de2fae686) 
1 #include <petsc/private/dmpleximpl.h>   /*I      "petscdmplex.h"   I*/
2 
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
9 
10 #if defined(PETSC_HAVE_TETGEN_TETLIBRARY_NEEDED)
11 #define TETLIBRARY
12 #endif
13 #include <tetgen.h>
14 
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;
19 
20   PetscFunctionBegin;
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   PetscFunctionReturn(0);
25 }
26 
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;
39 
40   PetscFunctionBegin;
41   ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr);
42   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
43   ierr = DMPlexIsInterpolatedCollective(boundary, &isInterpolated);CHKERRQ(ierr);
44   ierr = DMUniversalLabelCreate(boundary, &universal);CHKERRQ(ierr);
45 
46   ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr);
47   in.numberofpoints = vEnd - vStart;
48   if (in.numberofpoints > 0) {
49     PetscSection       coordSection;
50     Vec                coordinates;
51     const PetscScalar *array;
52 
53     in.pointlist       = new double[in.numberofpoints*dim];
54     in.pointmarkerlist = new int[in.numberofpoints];
55 
56     ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr);
57     ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr);
58     ierr = VecGetArrayRead(coordinates, &array);CHKERRQ(ierr);
59     for (v = vStart; v < vEnd; ++v) {
60       const PetscInt idx = v - vStart;
61       PetscInt       off, d, val;
62 
63       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
64       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
65       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
66       in.pointmarkerlist[idx] = (int) val;
67     }
68     ierr = VecRestoreArrayRead(coordinates, &array);CHKERRQ(ierr);
69   }
70 
71   ierr = DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);CHKERRQ(ierr);
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;
80 
81       ierr = DMPlexGetConeSize(boundary, e, &coneSize);CHKERRQ(ierr);
82       ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr);
83       in.edgelist[idx*2]     = cone[0] - vStart;
84       in.edgelist[idx*2 + 1] = cone[1] - vStart;
85 
86       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
87       in.edgemarkerlist[idx] = (int) val;
88     }
89   }
90 
91   ierr = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr);
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;
99 
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;
104 
105       ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
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       }
110 
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       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
119       in.facetmarkerlist[idx] = (int) val;
120       ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
121     }
122   }
123   if (rank == 0) {
124     DM_Plex *mesh = (DM_Plex *) boundary->data;
125     char     args[32];
126 
127     /* Take away 'Q' for verbose output */
128 #ifdef PETSC_HAVE_EGADS
129     ierr = PetscStrcpy(args, "pqezQY");CHKERRQ(ierr);
130 #else
131     ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr);
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;
142 
143     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
144       meshCoords = (PetscReal *) out.pointlist;
145     } else {
146       PetscInt i;
147 
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;
155 
156       cells = new PetscInt[numCells * numCorners];
157       for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out.tetrahedronlist[i];
158     }
159 
160     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
161     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr);
162 
163     /* Set labels */
164     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);CHKERRQ(ierr);
165     for (v = 0; v < numVertices; ++v) {
166       if (out.pointmarkerlist[v]) {
167         ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
168       }
169     }
170     if (interpolate) {
171       PetscInt e;
172 
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;
178 
179           ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
180           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
181           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
182           ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
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;
190 
191           ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
192           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
193           ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
194           ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
195         }
196       }
197     }
198 
199     ierr = PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
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;
208 
209       /* Get Attached EGADS Model from Original DMPlex */
210       ierr = PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
211       if (modelObj) {
212         ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
213         ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
214         /* Transfer EGADS Model to Volumetric Mesh */
215         ierr = PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr);
216       } else {
217         ierr = PetscObjectQuery((PetscObject) boundary, "EGADSLite Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
218         if (modelObj) {
219           ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
220           ierr = EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
221           /* Transfer EGADS Model to Volumetric Mesh */
222           ierr = PetscObjectCompose((PetscObject) *dm, "EGADSLite Model", (PetscObject) modelObj);CHKERRQ(ierr);
223           islite = PETSC_TRUE;
224         }
225       }
226       if (!modelObj) goto skip_egads;
227 
228       /* Set Cell Labels */
229       ierr = DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr);
230       ierr = DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
231       ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
232       ierr = DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
233 
234       for (c = cStart; c < cEnd; ++c) {
235         PetscReal centroid[3] = {0., 0., 0.};
236         PetscInt  b;
237 
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;
244 
245           ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr);
246           ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
247           ierr = DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
248           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
249           ierr = DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
250         } else {
251           ierr = DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);CHKERRQ(ierr);
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;
260 
261           ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr);
262           ierr = DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
263           for (cl = 0; cl < Ncl; cl += 2) {
264             const PetscInt p = closure[cl];
265 
266             if (p >= eStart && p < eEnd) {
267               ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr);
268               if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
269             }
270             if (p >= fStart && p < fEnd) {
271               ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr);
272               if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
273             }
274           }
275           ierr = DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
276         }
277       }
278 skip_egads: ;
279 #endif
280     }
281     ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
282   }
283   ierr = DMUniversalLabelDestroy(&universal);CHKERRQ(ierr);
284   PetscFunctionReturn(0);
285 }
286 
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;
299 
300   PetscFunctionBegin;
301   ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
302   ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr);
303   ierr = DMPlexIsInterpolatedCollective(dm, &isInterpolated);CHKERRQ(ierr);
304   ierr = DMUniversalLabelCreate(dm, &universal);CHKERRQ(ierr);
305 
306   ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
307   in.numberofpoints = vEnd - vStart;
308   if (in.numberofpoints > 0) {
309     PetscSection coordSection;
310     Vec          coordinates;
311     PetscScalar *array;
312 
313     in.pointlist       = new double[in.numberofpoints*dim];
314     in.pointmarkerlist = new int[in.numberofpoints];
315 
316     ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr);
317     ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
318     ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr);
319     for (v = vStart; v < vEnd; ++v) {
320       const PetscInt idx = v - vStart;
321       PetscInt       off, d, val;
322 
323       ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
324       for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
325       ierr = DMLabelGetValue(universal->label, v, &val);CHKERRQ(ierr);
326       in.pointmarkerlist[idx] = (int) val;
327     }
328     ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr);
329   }
330 
331   ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr);
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;
340 
341       ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr);
342       ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr);
343       in.edgelist[idx*2]     = cone[0] - vStart;
344       in.edgelist[idx*2 + 1] = cone[1] - vStart;
345 
346       ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr);
347       in.edgemarkerlist[idx] = (int) val;
348     }
349   }
350 
351   ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
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;
359 
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;
364 
365       ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
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       }
370 
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       }
378 
379       ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr);
380       in.facetmarkerlist[idx] = (int) val;
381 
382       ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr);
383     }
384   }
385 
386   ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
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;
396 
397       ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
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       ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
401     }
402   }
403 
404   if (rank == 0) {
405     char args[32];
406 
407     /* Take away 'Q' for verbose output */
408     ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr);
409     ::tetrahedralize(args, &in, &out);
410   }
411 
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;
420 
421     if (sizeof (PetscReal) == sizeof (out.pointlist[0])) {
422       meshCoords = (PetscReal *) out.pointlist;
423     } else {
424       PetscInt i;
425 
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;
433 
434       cells = new PetscInt[numCells * numCorners];
435       for (i = 0; i < numCells * numCorners; ++i)cells[i] = (PetscInt) out.tetrahedronlist[i];
436     }
437 
438     ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr);
439     ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr);
440     if (sizeof (PetscReal) != sizeof (out.pointlist[0])) {delete [] meshCoords;}
441     if (sizeof (PetscInt) != sizeof (out.tetrahedronlist[0])) {delete [] cells;}
442 
443     /* Set labels */
444     ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);CHKERRQ(ierr);
445     for (v = 0; v < numVertices; ++v) {
446       if (out.pointmarkerlist[v]) {
447         ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);
448       }
449     }
450     if (interpolate) {
451       PetscInt e, f;
452 
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;
458 
459           ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
460           if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
461           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);
462           ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr);
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;
470 
471           ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
472           if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
473           ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);
474           ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr);
475         }
476       }
477     }
478 
479     ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
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;
488 
489       /* Get Attached EGADS Model from Original DMPlex */
490       ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
491       if (modelObj) {
492         ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
493         ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
494         /* Transfer EGADS Model to Volumetric Mesh */
495         ierr = PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr);
496       } else {
497         ierr = PetscObjectQuery((PetscObject) dm, "EGADSLite Model", (PetscObject *) &modelObj);CHKERRQ(ierr);
498         if (modelObj) {
499           ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr);
500           ierr = EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr);
501           /* Transfer EGADS Model to Volumetric Mesh */
502           ierr = PetscObjectCompose((PetscObject) *dmRefined, "EGADSLite Model", (PetscObject) modelObj);CHKERRQ(ierr);
503           islite = PETSC_TRUE;
504         }
505       }
506       if (!modelObj) goto skip_egads;
507 
508       /* Set Cell Labels */
509       ierr = DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr);
510       ierr = DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);CHKERRQ(ierr);
511       ierr = DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);CHKERRQ(ierr);
512       ierr = DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);CHKERRQ(ierr);
513 
514       for (c = cStart; c < cEnd; ++c) {
515         PetscReal centroid[3] = {0., 0., 0.};
516         PetscInt  b;
517 
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;
524 
525           ierr = DMGetCoordinatesLocal(*dmRefined, &coordinates);CHKERRQ(ierr);
526           ierr = DMGetCoordinateSection(*dmRefined, &coordSection);CHKERRQ(ierr);
527           ierr = DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
528           for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
529           ierr = DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr);
530         } else {
531           ierr = DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);CHKERRQ(ierr);
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;
540 
541           ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr);
542           ierr = DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
543           for (cl = 0; cl < Ncl; cl += 2) {
544             const PetscInt p = closure[cl];
545 
546             if (p >= eStart && p < eEnd) {
547               ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr);
548               if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
549             }
550             if (p >= fStart && p < fEnd) {
551               ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr);
552               if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);}
553             }
554           }
555           ierr = DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr);
556         }
557       }
558 skip_egads: ;
559 #endif
560     }
561     ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr);
562   }
563   PetscFunctionReturn(0);
564 }
565