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