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;
54: PetscBool flg;
55: char opts[64];
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));
63: PetscCall(PetscOptionsGetString(((PetscObject)boundary)->options, ((PetscObject)boundary)->prefix, "-dm_plex_generate_tetgen_opts", opts, sizeof(opts), &flg));
64: if (flg) PetscCall(DMPlexTetgenSetOptions(boundary, opts));
66: PetscCall(DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd));
67: in.numberofpoints = vEnd - vStart;
68: if (in.numberofpoints > 0) {
69: PetscSection coordSection;
70: Vec coordinates;
71: const PetscScalar *array;
73: in.pointlist = new double[in.numberofpoints * dim];
74: in.pointmarkerlist = new int[in.numberofpoints];
76: PetscCall(PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints));
77: PetscCall(DMGetCoordinatesLocal(boundary, &coordinates));
78: PetscCall(DMGetCoordinateSection(boundary, &coordSection));
79: PetscCall(VecGetArrayRead(coordinates, &array));
80: for (v = vStart; v < vEnd; ++v) {
81: const PetscInt idx = v - vStart;
82: PetscInt off, d, val;
84: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
85: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
86: PetscCall(DMLabelGetValue(universal->label, v, &val));
87: if (val != defVal) in.pointmarkerlist[idx] = (int)val;
88: }
89: PetscCall(VecRestoreArrayRead(coordinates, &array));
90: }
92: PetscCall(DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd));
93: in.numberofedges = eEnd - eStart;
94: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
95: in.edgelist = new int[in.numberofedges * 2];
96: in.edgemarkerlist = new int[in.numberofedges];
97: for (e = eStart; e < eEnd; ++e) {
98: const PetscInt idx = e - eStart;
99: const PetscInt *cone;
100: PetscInt coneSize, val;
102: PetscCall(DMPlexGetConeSize(boundary, e, &coneSize));
103: PetscCall(DMPlexGetCone(boundary, e, &cone));
104: in.edgelist[idx * 2] = cone[0] - vStart;
105: in.edgelist[idx * 2 + 1] = cone[1] - vStart;
107: PetscCall(DMLabelGetValue(universal->label, e, &val));
108: if (val != defVal) in.edgemarkerlist[idx] = (int)val;
109: }
110: }
112: PetscCall(DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd));
113: in.numberoffacets = fEnd - fStart;
114: if (in.numberoffacets > 0) {
115: in.facetlist = new tetgenio::facet[in.numberoffacets];
116: in.facetmarkerlist = new int[in.numberoffacets];
117: for (f = fStart; f < fEnd; ++f) {
118: const PetscInt idx = f - fStart;
119: PetscInt *points = nullptr, numPoints, p, numVertices = 0, v, val = -1;
121: in.facetlist[idx].numberofpolygons = 1;
122: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
123: in.facetlist[idx].numberofholes = 0;
124: in.facetlist[idx].holelist = nullptr;
126: PetscCall(DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
127: for (p = 0; p < numPoints * 2; p += 2) {
128: const PetscInt point = points[p];
129: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
130: }
132: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
133: poly->numberofvertices = numVertices;
134: poly->vertexlist = new int[poly->numberofvertices];
135: for (v = 0; v < numVertices; ++v) {
136: const PetscInt vIdx = points[v] - vStart;
137: poly->vertexlist[v] = vIdx;
138: }
139: PetscCall(DMLabelGetValue(universal->label, f, &val));
140: if (val != defVal) in.facetmarkerlist[idx] = (int)val;
141: PetscCall(DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
142: }
143: }
144: if (rank == 0) {
145: DM_Plex *mesh = (DM_Plex *)boundary->data;
146: char args[32];
148: /* Take away 'Q' for verbose output */
149: #ifdef PETSC_HAVE_EGADS
150: if (mesh->tetgenRadiusEdgeBound > 0.) {
151: PetscCall(PetscSNPrintf(args, sizeof(args), "pYq%.2f/%.0fezQY", (double)mesh->tetgenRadiusEdgeBound, (double)mesh->tetgenDihedralBound));
152: } else {
153: PetscCall(PetscStrncpy(args, "pYqezQY", sizeof(args)));
154: }
155: #else
156: if (mesh->tetgenRadiusEdgeBound > 0.) {
157: PetscCall(PetscSNPrintf(args, sizeof(args), "pq%.2f/%.0fezQ", (double)mesh->tetgenRadiusEdgeBound, (double)mesh->tetgenDihedralBound));
158: } else {
159: PetscCall(PetscStrncpy(args, "pqezQ", sizeof(args)));
160: }
161: #endif
162: if (mesh->tetgenOpts) {
163: ::tetrahedralize(mesh->tetgenOpts, &in, &out);
164: } else {
165: ::tetrahedralize(args, &in, &out);
166: }
167: }
168: {
169: const PetscInt numCorners = 4;
170: const PetscInt numCells = out.numberoftetrahedra;
171: const PetscInt numVertices = out.numberofpoints;
172: PetscReal *meshCoords = nullptr;
173: PetscInt *cells = nullptr;
175: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
176: meshCoords = (PetscReal *)out.pointlist;
177: } else {
178: PetscInt i;
180: meshCoords = new PetscReal[dim * numVertices];
181: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out.pointlist[i];
182: }
183: if (sizeof(PetscInt) == sizeof(out.tetrahedronlist[0])) {
184: cells = (PetscInt *)out.tetrahedronlist;
185: } else {
186: PetscInt i;
188: cells = new PetscInt[numCells * numCorners];
189: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.tetrahedronlist[i];
190: }
192: PetscCall(DMPlexInvertCells_Tetgen(numCells, numCorners, cells));
193: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
195: /* Set labels */
196: PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm));
197: for (v = 0; v < numVertices; ++v) {
198: if (out.pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out.pointmarkerlist[v]));
199: }
200: if (interpolate) {
201: PetscInt e;
203: for (e = 0; e < out.numberofedges; e++) {
204: if (out.edgemarkerlist[e]) {
205: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
206: const PetscInt *edges;
207: PetscInt numEdges;
209: PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
210: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
211: PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]));
212: PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
213: }
214: }
215: for (f = 0; f < out.numberoftrifaces; f++) {
216: if (out.trifacemarkerlist[f]) {
217: const PetscInt vertices[3] = {out.trifacelist[f * 3 + 0] + numCells, out.trifacelist[f * 3 + 1] + numCells, out.trifacelist[f * 3 + 2] + numCells};
218: const PetscInt *faces;
219: PetscInt numFaces;
221: PetscCall(DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces));
222: PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
223: PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]));
224: PetscCall(DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces));
225: }
226: }
227: }
229: PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj));
230: if (!modelObj) PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADSlite Model", (PetscObject *)&modelObj));
232: if (modelObj) {
233: #ifdef PETSC_HAVE_EGADS
234: DMLabel bodyLabel;
235: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
236: PetscBool islite = PETSC_FALSE;
237: ego *bodies;
238: ego model, geom;
239: int Nb, oclass, mtype, *senses;
241: PetscCall(DMPlexCopyEGADSInfo_Internal(boundary, *dm));
243: // Get Attached EGADS Model from Original DMPlex
244: PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj));
245: if (modelObj) {
246: PetscCall(PetscContainerGetPointer(modelObj, &model));
247: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, nullptr, &Nb, &bodies, &senses));
248: } else {
249: PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADSlite Model", (PetscObject *)&modelObj));
250: if (modelObj) {
251: PetscCall(PetscContainerGetPointer(modelObj, &model));
252: PetscCall(EGlite_getTopology(model, &geom, &oclass, &mtype, nullptr, &Nb, &bodies, &senses));
253: islite = PETSC_TRUE;
254: }
255: }
256: if (!modelObj) goto skip_egads;
258: /* Set Cell Labels */
259: PetscCall(DMGetLabel(*dm, "EGADS Body ID", &bodyLabel));
260: PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
261: PetscCall(DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd));
262: PetscCall(DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd));
264: for (c = cStart; c < cEnd; ++c) {
265: PetscReal centroid[3] = {0., 0., 0.};
266: PetscInt b;
268: /* Determine what body the cell's centroid is located in */
269: if (!interpolate) {
270: PetscSection coordSection;
271: Vec coordinates;
272: PetscScalar *coords = nullptr;
273: PetscInt coordSize, s, d;
275: PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
276: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
277: PetscCall(DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
278: for (s = 0; s < coordSize; ++s)
279: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
280: PetscCall(DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
281: } else PetscCall(DMPlexComputeCellGeometryFVM(*dm, c, nullptr, centroid, nullptr));
282: for (b = 0; b < Nb; ++b) {
283: if (islite) {
284: if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
285: } else {
286: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
287: }
288: }
289: if (b < Nb) {
290: PetscInt cval = b, eVal, fVal;
291: PetscInt *closure = nullptr, Ncl, cl;
293: PetscCall(DMLabelSetValue(bodyLabel, c, cval));
294: PetscCall(DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
295: for (cl = 0; cl < Ncl; cl += 2) {
296: const PetscInt p = closure[cl];
298: if (p >= eStart && p < eEnd) {
299: PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
300: if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
301: }
302: if (p >= fStart && p < fEnd) {
303: PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
304: if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
305: }
306: }
307: PetscCall(DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
308: }
309: }
310: skip_egads:;
311: #endif
312: }
313: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
314: }
315: PetscCall(DMUniversalLabelDestroy(&universal));
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
320: {
321: MPI_Comm comm;
322: const PetscInt dim = 3;
323: ::tetgenio in;
324: ::tetgenio out;
325: PetscContainer modelObj;
326: DMUniversalLabel universal;
327: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal;
328: DMPlexInterpolatedFlag isInterpolated;
329: PetscMPIInt rank;
331: PetscFunctionBegin;
332: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
333: PetscCallMPI(MPI_Comm_rank(comm, &rank));
334: PetscCall(DMPlexIsInterpolatedCollective(dm, &isInterpolated));
335: PetscCall(DMUniversalLabelCreate(dm, &universal));
336: PetscCall(DMLabelGetDefaultValue(universal->label, &defVal));
338: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
339: in.numberofpoints = vEnd - vStart;
340: if (in.numberofpoints > 0) {
341: PetscSection coordSection;
342: Vec coordinates;
343: PetscScalar *array;
345: in.pointlist = new double[in.numberofpoints * dim];
346: in.pointmarkerlist = new int[in.numberofpoints];
348: PetscCall(PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints));
349: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
350: PetscCall(DMGetCoordinateSection(dm, &coordSection));
351: PetscCall(VecGetArray(coordinates, &array));
352: for (v = vStart; v < vEnd; ++v) {
353: const PetscInt idx = v - vStart;
354: PetscInt off, d, val;
356: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
357: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
358: PetscCall(DMLabelGetValue(universal->label, v, &val));
359: if (val != defVal) in.pointmarkerlist[idx] = (int)val;
360: }
361: PetscCall(VecRestoreArray(coordinates, &array));
362: }
364: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
365: in.numberofedges = eEnd - eStart;
366: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
367: in.edgelist = new int[in.numberofedges * 2];
368: in.edgemarkerlist = new int[in.numberofedges];
369: for (e = eStart; e < eEnd; ++e) {
370: const PetscInt idx = e - eStart;
371: const PetscInt *cone;
372: PetscInt coneSize, val;
374: PetscCall(DMPlexGetConeSize(dm, e, &coneSize));
375: PetscCall(DMPlexGetCone(dm, e, &cone));
376: in.edgelist[idx * 2] = cone[0] - vStart;
377: in.edgelist[idx * 2 + 1] = cone[1] - vStart;
379: PetscCall(DMLabelGetValue(universal->label, e, &val));
380: if (val != defVal) in.edgemarkerlist[idx] = (int)val;
381: }
382: }
384: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
385: in.numberoffacets = fEnd - fStart;
386: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
387: in.facetlist = new tetgenio::facet[in.numberoffacets];
388: in.facetmarkerlist = new int[in.numberoffacets];
389: for (f = fStart; f < fEnd; ++f) {
390: const PetscInt idx = f - fStart;
391: PetscInt *points = nullptr, numPoints, p, numVertices = 0, v, val;
393: in.facetlist[idx].numberofpolygons = 1;
394: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
395: in.facetlist[idx].numberofholes = 0;
396: in.facetlist[idx].holelist = nullptr;
398: PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
399: for (p = 0; p < numPoints * 2; p += 2) {
400: const PetscInt point = points[p];
401: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
402: }
404: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
405: poly->numberofvertices = numVertices;
406: poly->vertexlist = new int[poly->numberofvertices];
407: for (v = 0; v < numVertices; ++v) {
408: const PetscInt vIdx = points[v] - vStart;
409: poly->vertexlist[v] = vIdx;
410: }
412: PetscCall(DMLabelGetValue(universal->label, f, &val));
413: if (val != defVal) in.facetmarkerlist[idx] = (int)val;
415: PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
416: }
417: }
419: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
420: in.numberofcorners = 4;
421: in.numberoftetrahedra = cEnd - cStart;
422: in.tetrahedronvolumelist = (double *)maxVolumes;
423: if (in.numberoftetrahedra > 0) {
424: in.tetrahedronlist = new int[in.numberoftetrahedra * in.numberofcorners];
425: for (c = cStart; c < cEnd; ++c) {
426: const PetscInt idx = c - cStart;
427: PetscInt *closure = nullptr;
428: PetscInt closureSize;
430: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
431: PetscCheck(!(closureSize != 5) || !(closureSize != 15), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %" PetscInt_FMT " vertices in closure", closureSize);
432: for (v = 0; v < 4; ++v) in.tetrahedronlist[idx * in.numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
433: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
434: }
435: }
437: if (rank == 0) {
438: char args[32];
440: /* Take away 'Q' for verbose output */
441: PetscCall(PetscStrncpy(args, "qezQra", sizeof(args)));
442: ::tetrahedralize(args, &in, &out);
443: }
445: in.tetrahedronvolumelist = nullptr;
446: {
447: const PetscInt numCorners = 4;
448: const PetscInt numCells = out.numberoftetrahedra;
449: const PetscInt numVertices = out.numberofpoints;
450: PetscReal *meshCoords = nullptr;
451: PetscInt *cells = nullptr;
452: PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
454: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
455: meshCoords = (PetscReal *)out.pointlist;
456: } else {
457: PetscInt i;
459: meshCoords = new PetscReal[dim * numVertices];
460: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out.pointlist[i];
461: }
462: if (sizeof(PetscInt) == sizeof(out.tetrahedronlist[0])) {
463: cells = (PetscInt *)out.tetrahedronlist;
464: } else {
465: PetscInt i;
467: cells = new PetscInt[numCells * numCorners];
468: for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out.tetrahedronlist[i];
469: }
471: PetscCall(DMPlexInvertCells_Tetgen(numCells, numCorners, cells));
472: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
473: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) delete[] meshCoords;
474: if (sizeof(PetscInt) != sizeof(out.tetrahedronlist[0])) delete[] cells;
476: /* Set labels */
477: PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined));
478: for (v = 0; v < numVertices; ++v) {
479: if (out.pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out.pointmarkerlist[v]));
480: }
481: if (interpolate) {
482: PetscInt e, f;
484: for (e = 0; e < out.numberofedges; ++e) {
485: if (out.edgemarkerlist[e]) {
486: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
487: const PetscInt *edges;
488: PetscInt numEdges;
490: PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
491: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
492: PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]));
493: PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
494: }
495: }
496: for (f = 0; f < out.numberoftrifaces; ++f) {
497: if (out.trifacemarkerlist[f]) {
498: const PetscInt vertices[3] = {out.trifacelist[f * 3 + 0] + numCells, out.trifacelist[f * 3 + 1] + numCells, out.trifacelist[f * 3 + 2] + numCells};
499: const PetscInt *faces;
500: PetscInt numFaces;
502: PetscCall(DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces));
503: PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
504: PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]));
505: PetscCall(DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces));
506: }
507: }
508: }
510: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
511: if (modelObj) {
512: #ifdef PETSC_HAVE_EGADS
513: DMLabel bodyLabel;
514: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
515: PetscBool islite = PETSC_FALSE;
516: ego *bodies;
517: ego model, geom;
518: int Nb, oclass, mtype, *senses;
520: PetscCall(DMPlexCopyEGADSInfo_Internal(dm, *dmRefined));
522: /* Get Attached EGADS Model from Original DMPlex */
523: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
524: if (modelObj) {
525: PetscCall(PetscContainerGetPointer(modelObj, &model));
526: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, nullptr, &Nb, &bodies, &senses));
527: } else {
528: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", (PetscObject *)&modelObj));
529: if (modelObj) {
530: PetscCall(PetscContainerGetPointer(modelObj, &model));
531: PetscCall(EGlite_getTopology(model, &geom, &oclass, &mtype, nullptr, &Nb, &bodies, &senses));
532: islite = PETSC_TRUE;
533: }
534: }
535: if (!modelObj) goto skip_egads;
537: /* Set Cell Labels */
538: PetscCall(DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel));
539: PetscCall(DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd));
540: PetscCall(DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd));
541: PetscCall(DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd));
543: for (c = cStart; c < cEnd; ++c) {
544: PetscReal centroid[3] = {0., 0., 0.};
545: PetscInt b;
547: /* Determine what body the cell's centroid is located in */
548: if (!interpolate) {
549: PetscSection coordSection;
550: Vec coordinates;
551: PetscScalar *coords = nullptr;
552: PetscInt coordSize, s, d;
554: PetscCall(DMGetCoordinatesLocal(*dmRefined, &coordinates));
555: PetscCall(DMGetCoordinateSection(*dmRefined, &coordSection));
556: PetscCall(DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
557: for (s = 0; s < coordSize; ++s)
558: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
559: PetscCall(DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
560: } else PetscCall(DMPlexComputeCellGeometryFVM(*dmRefined, c, nullptr, centroid, nullptr));
561: for (b = 0; b < Nb; ++b) {
562: if (islite) {
563: if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
564: } else {
565: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
566: }
567: }
568: if (b < Nb) {
569: PetscInt cval = b, eVal, fVal;
570: PetscInt *closure = nullptr, Ncl, cl;
572: PetscCall(DMLabelSetValue(bodyLabel, c, cval));
573: PetscCall(DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
574: for (cl = 0; cl < Ncl; cl += 2) {
575: const PetscInt p = closure[cl];
577: if (p >= eStart && p < eEnd) {
578: PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
579: if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
580: }
581: if (p >= fStart && p < fEnd) {
582: PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
583: if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
584: }
585: }
586: PetscCall(DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
587: }
588: }
589: skip_egads:;
590: #endif
591: }
592: PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
593: }
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }