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