Actual source code: ctetgenerate.c
1: #include <petsc/private/dmpleximpl.h>
3: #ifdef PETSC_HAVE_EGADS
4: #include <egads.h>
5: #endif
7: #include <ctetgen.h>
9: /* This is to fix the tetrahedron orientation from TetGen */
10: static PetscErrorCode DMPlexInvertCells_CTetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
11: {
12: PetscInt bound = numCells * numCorners, coff;
14: PetscFunctionBegin;
15: #define SWAP(a, b) \
16: do { \
17: PetscInt tmp = (a); \
18: (a) = (b); \
19: (b) = tmp; \
20: } while (0)
21: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
22: #undef SWAP
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
27: {
28: MPI_Comm comm;
29: const PetscInt dim = 3;
30: PLC *in, *out;
31: DMUniversalLabel universal;
32: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, defVal, verbose = 0;
33: DMPlexInterpolatedFlag isInterpolated;
34: PetscMPIInt rank;
36: PetscFunctionBegin;
37: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)boundary)->prefix, "-ctetgen_verbose", &verbose, NULL));
38: PetscCall(PetscObjectGetComm((PetscObject)boundary, &comm));
39: PetscCallMPI(MPI_Comm_rank(comm, &rank));
40: PetscCall(DMPlexIsInterpolatedCollective(boundary, &isInterpolated));
41: PetscCall(DMUniversalLabelCreate(boundary, &universal));
42: PetscCall(DMLabelGetDefaultValue(universal->label, &defVal));
44: PetscCall(PLCCreate(&in));
45: PetscCall(PLCCreate(&out));
47: PetscCall(DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd));
48: in->numberofpoints = vEnd - vStart;
49: if (in->numberofpoints > 0) {
50: PetscSection coordSection;
51: Vec coordinates;
52: const PetscScalar *array;
54: PetscCall(PetscMalloc1(in->numberofpoints * dim, &in->pointlist));
55: PetscCall(PetscCalloc1(in->numberofpoints, &in->pointmarkerlist));
56: PetscCall(DMGetCoordinatesLocal(boundary, &coordinates));
57: PetscCall(DMGetCoordinateSection(boundary, &coordSection));
58: PetscCall(VecGetArrayRead(coordinates, &array));
59: for (v = vStart; v < vEnd; ++v) {
60: const PetscInt idx = v - vStart;
61: PetscInt off, d, m;
63: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
64: for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
65: PetscCall(DMLabelGetValue(universal->label, v, &m));
66: if (m != defVal) in->pointmarkerlist[idx] = (int)m;
67: }
68: PetscCall(VecRestoreArrayRead(coordinates, &array));
69: }
71: PetscCall(DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd));
72: in->numberofedges = eEnd - eStart;
73: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
74: PetscCall(PetscMalloc1(in->numberofedges * 2, &in->edgelist));
75: PetscCall(PetscMalloc1(in->numberofedges, &in->edgemarkerlist));
76: for (e = eStart; e < eEnd; ++e) {
77: const PetscInt idx = e - eStart;
78: const PetscInt *cone;
79: PetscInt coneSize, val;
81: PetscCall(DMPlexGetConeSize(boundary, e, &coneSize));
82: PetscCall(DMPlexGetCone(boundary, e, &cone));
83: in->edgelist[idx * 2] = cone[0] - vStart;
84: in->edgelist[idx * 2 + 1] = cone[1] - vStart;
86: PetscCall(DMLabelGetValue(universal->label, e, &val));
87: if (val != defVal) in->edgemarkerlist[idx] = (int)val;
88: }
89: }
91: PetscCall(DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd));
92: in->numberoffacets = fEnd - fStart;
93: if (in->numberoffacets > 0) {
94: PetscCall(PetscMalloc1(in->numberoffacets, &in->facetlist));
95: PetscCall(PetscMalloc1(in->numberoffacets, &in->facetmarkerlist));
96: for (f = fStart; f < fEnd; ++f) {
97: const PetscInt idx = f - fStart;
98: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
99: polygon *poly;
101: in->facetlist[idx].numberofpolygons = 1;
102: PetscCall(PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist));
103: in->facetlist[idx].numberofholes = 0;
104: in->facetlist[idx].holelist = NULL;
106: PetscCall(DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
107: for (p = 0; p < numPoints * 2; p += 2) {
108: const PetscInt point = points[p];
109: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
110: }
112: poly = in->facetlist[idx].polygonlist;
113: poly->numberofvertices = numVertices;
114: PetscCall(PetscMalloc1(poly->numberofvertices, &poly->vertexlist));
115: for (v = 0; v < numVertices; ++v) {
116: const PetscInt vIdx = points[v] - vStart;
117: poly->vertexlist[v] = vIdx;
118: }
119: PetscCall(DMLabelGetValue(universal->label, f, &m));
120: if (m != defVal) in->facetmarkerlist[idx] = (int)m;
121: PetscCall(DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points));
122: }
123: }
124: if (rank == 0) {
125: TetGenOpts t;
126: PetscReal minratio, mindihedral;
128: PetscCall(TetGenOptsInitialize(&t));
129: t.in = boundary; /* Should go away */
130: t.plc = 1;
131: t.quality = 1;
132: t.edgesout = 1;
133: t.zeroindex = 1;
134: t.quiet = 1;
135: t.verbose = verbose;
136: PetscCall(DMPlexTetgenGetRadiusEdgeBound(boundary, &minratio));
137: PetscCall(DMPlexTetgenGetDihedralBound(boundary, &mindihedral));
138: if (minratio > 0.) t.minratio = minratio;
139: if (mindihedral > 0.) t.mindihedral = mindihedral;
140: #if 0
141: #ifdef PETSC_HAVE_EGADS
142: /* Need to add in more TetGen code */
143: t.nobisect = 1; /* Add Y to preserve Surface Mesh for EGADS */
144: #endif
145: #endif
147: PetscCall(TetGenCheckOpts(&t));
148: PetscCall(TetGenTetrahedralize(&t, in, out));
149: }
150: {
151: const PetscInt numCorners = 4;
152: const PetscInt numCells = out->numberoftetrahedra;
153: const PetscInt numVertices = out->numberofpoints;
154: PetscReal *meshCoords = NULL;
155: PetscInt *cells = NULL;
157: if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
158: meshCoords = (PetscReal *)out->pointlist;
159: } else {
160: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
161: for (PetscInt i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out->pointlist[i];
162: }
163: if (sizeof(PetscInt) == sizeof(out->tetrahedronlist[0])) {
164: cells = (PetscInt *)out->tetrahedronlist;
165: } else {
166: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
167: for (PetscInt i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out->tetrahedronlist[i];
168: }
170: PetscCall(DMPlexInvertCells_CTetgen(numCells, numCorners, cells));
171: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm));
172: if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscCall(PetscFree(meshCoords));
173: if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscCall(PetscFree(cells));
175: /* Set labels */
176: PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm));
177: for (v = 0; v < numVertices; ++v) {
178: if (out->pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]));
179: }
180: if (interpolate) {
181: for (PetscInt e = 0; e < out->numberofedges; e++) {
182: if (out->edgemarkerlist[e]) {
183: const PetscInt vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
184: const PetscInt *edges;
185: PetscInt numEdges;
187: PetscCall(DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges));
188: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
189: PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]));
190: PetscCall(DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges));
191: }
192: }
193: for (f = 0; f < out->numberoftrifaces; f++) {
194: if (out->trifacemarkerlist[f]) {
195: const PetscInt vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
196: const PetscInt *faces;
197: PetscInt numFaces;
199: PetscCall(DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces));
200: PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
201: PetscCall(DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]));
202: PetscCall(DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces));
203: }
204: }
205: }
207: #ifdef PETSC_HAVE_EGADS
208: {
209: DMLabel bodyLabel;
210: PetscContainer modelObj;
211: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
212: ego *bodies;
213: ego model, geom;
214: int Nb, oclass, mtype, *senses;
216: /* Get Attached EGADS Model from Original DMPlex */
217: PetscCall(PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj));
218: if (modelObj) {
219: PetscCall(PetscContainerGetPointer(modelObj, &model));
220: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
221: /* Transfer EGADS Model to Volumetric Mesh */
222: PetscCall(PetscObjectCompose((PetscObject)*dm, "EGADS Model", (PetscObject)modelObj));
224: /* Set Cell Labels */
225: PetscCall(DMGetLabel(*dm, "EGADS Body ID", &bodyLabel));
226: PetscCall(DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd));
227: PetscCall(DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd));
228: PetscCall(DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd));
230: for (c = cStart; c < cEnd; ++c) {
231: PetscReal centroid[3] = {0., 0., 0.};
232: PetscInt b;
234: /* Determine what body the cell's centroid is located in */
235: if (!interpolate) {
236: PetscSection coordSection;
237: Vec coordinates;
238: PetscScalar *coords = NULL;
239: PetscInt coordSize, s, d;
241: PetscCall(DMGetCoordinatesLocal(*dm, &coordinates));
242: PetscCall(DMGetCoordinateSection(*dm, &coordSection));
243: PetscCall(DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
244: for (s = 0; s < coordSize; ++s)
245: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
246: PetscCall(DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords));
247: } else PetscCall(DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL));
248: for (b = 0; b < Nb; ++b) {
249: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
250: }
251: if (b < Nb) {
252: PetscInt cval = b, eVal, fVal;
253: PetscInt *closure = NULL, Ncl, cl;
255: PetscCall(DMLabelSetValue(bodyLabel, c, cval));
256: PetscCall(DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
257: for (cl = 0; cl < Ncl; ++cl) {
258: const PetscInt p = closure[cl * 2];
260: if (p >= eStart && p < eEnd) {
261: PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
262: if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
263: }
264: if (p >= fStart && p < fEnd) {
265: PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
266: if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
267: }
268: }
269: PetscCall(DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure));
270: }
271: }
272: }
273: }
274: #endif
275: PetscCall(DMPlexSetRefinementUniform(*dm, PETSC_FALSE));
276: }
278: PetscCall(DMUniversalLabelDestroy(&universal));
279: PetscCall(PLCDestroy(&in));
280: PetscCall(PLCDestroy(&out));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
285: {
286: MPI_Comm comm;
287: const PetscInt dim = 3;
288: PLC *in, *out;
289: DMUniversalLabel universal;
290: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal, verbose = 0;
291: DMPlexInterpolatedFlag isInterpolated;
292: PetscMPIInt rank;
294: PetscFunctionBegin;
295: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)dm)->prefix, "-ctetgen_verbose", &verbose, NULL));
296: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
297: PetscCallMPI(MPI_Comm_rank(comm, &rank));
298: PetscCall(DMPlexIsInterpolatedCollective(dm, &isInterpolated));
299: PetscCall(DMUniversalLabelCreate(dm, &universal));
300: PetscCall(DMLabelGetDefaultValue(universal->label, &defVal));
302: PetscCall(PLCCreate(&in));
303: PetscCall(PLCCreate(&out));
305: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
306: in->numberofpoints = vEnd - vStart;
307: if (in->numberofpoints > 0) {
308: PetscSection coordSection;
309: Vec coordinates;
310: PetscScalar *array;
312: PetscCall(PetscMalloc1(in->numberofpoints * dim, &in->pointlist));
313: PetscCall(PetscCalloc1(in->numberofpoints, &in->pointmarkerlist));
314: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
315: PetscCall(DMGetCoordinateSection(dm, &coordSection));
316: PetscCall(VecGetArray(coordinates, &array));
317: for (v = vStart; v < vEnd; ++v) {
318: const PetscInt idx = v - vStart;
319: PetscInt off, d, m;
321: PetscCall(PetscSectionGetOffset(coordSection, v, &off));
322: for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
323: PetscCall(DMLabelGetValue(universal->label, v, &m));
324: if (m != defVal) in->pointmarkerlist[idx] = (int)m;
325: }
326: PetscCall(VecRestoreArray(coordinates, &array));
327: }
329: PetscCall(DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd));
330: in->numberofedges = eEnd - eStart;
331: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
332: PetscCall(PetscMalloc1(in->numberofedges * 2, &in->edgelist));
333: PetscCall(PetscMalloc1(in->numberofedges, &in->edgemarkerlist));
334: for (e = eStart; e < eEnd; ++e) {
335: const PetscInt idx = e - eStart;
336: const PetscInt *cone;
337: PetscInt coneSize, val;
339: PetscCall(DMPlexGetConeSize(dm, e, &coneSize));
340: PetscCall(DMPlexGetCone(dm, e, &cone));
341: in->edgelist[idx * 2] = cone[0] - vStart;
342: in->edgelist[idx * 2 + 1] = cone[1] - vStart;
344: PetscCall(DMLabelGetValue(universal->label, e, &val));
345: if (val != defVal) in->edgemarkerlist[idx] = (int)val;
346: }
347: }
349: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
350: in->numberoftrifaces = 0;
351: for (f = fStart; f < fEnd; ++f) {
352: PetscInt supportSize;
354: PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
355: if (supportSize == 1) ++in->numberoftrifaces;
356: }
357: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
358: PetscInt tf = 0;
360: PetscCall(PetscMalloc1(in->numberoftrifaces * 3, &in->trifacelist));
361: PetscCall(PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist));
362: for (f = fStart; f < fEnd; ++f) {
363: PetscInt *points = NULL;
364: PetscInt supportSize, numPoints, p, Nv = 0, val;
366: PetscCall(DMPlexGetSupportSize(dm, f, &supportSize));
367: if (supportSize != 1) continue;
368: PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
369: for (p = 0; p < numPoints * 2; p += 2) {
370: const PetscInt point = points[p];
371: if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf * 3 + Nv++] = point - vStart;
372: }
373: PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points));
374: PetscCheck(Nv == 3, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %" PetscInt_FMT " has %" PetscInt_FMT " vertices, not 3", f, Nv);
375: PetscCall(DMLabelGetValue(universal->label, f, &val));
376: if (val != defVal) in->trifacemarkerlist[tf] = (int)val;
377: ++tf;
378: }
379: }
381: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
382: in->numberofcorners = 4;
383: in->numberoftetrahedra = cEnd - cStart;
384: in->tetrahedronvolumelist = maxVolumes;
385: if (in->numberoftetrahedra > 0) {
386: PetscCall(PetscMalloc1(in->numberoftetrahedra * in->numberofcorners, &in->tetrahedronlist));
387: for (c = cStart; c < cEnd; ++c) {
388: const PetscInt idx = c - cStart;
389: PetscInt *closure = NULL;
390: PetscInt closureSize;
392: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
393: PetscCheck((closureSize == 5) || (closureSize == 15), comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %" PetscInt_FMT " vertices in closure", closureSize);
394: for (v = 0; v < 4; ++v) in->tetrahedronlist[idx * in->numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
395: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure));
396: }
397: }
399: if (rank == 0) {
400: TetGenOpts t;
402: PetscCall(TetGenOptsInitialize(&t));
404: t.in = dm; /* Should go away */
405: t.refine = 1;
406: t.varvolume = 1;
407: t.quality = 1;
408: t.edgesout = 1;
409: t.zeroindex = 1;
410: t.quiet = 1;
411: t.verbose = verbose; /* Change this */
413: PetscCall(TetGenCheckOpts(&t));
414: PetscCall(TetGenTetrahedralize(&t, in, out));
415: }
417: in->tetrahedronvolumelist = NULL;
418: {
419: const PetscInt numCorners = 4;
420: const PetscInt numCells = out->numberoftetrahedra;
421: const PetscInt numVertices = out->numberofpoints;
422: PetscReal *meshCoords = NULL;
423: PetscInt *cells = NULL;
424: PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
426: if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
427: meshCoords = (PetscReal *)out->pointlist;
428: } else {
429: PetscCall(PetscMalloc1(dim * numVertices, &meshCoords));
430: for (PetscInt i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out->pointlist[i];
431: }
432: if (sizeof(PetscInt) == sizeof(out->tetrahedronlist[0])) {
433: cells = (PetscInt *)out->tetrahedronlist;
434: } else {
435: PetscCall(PetscMalloc1(numCells * numCorners, &cells));
436: for (PetscInt i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out->tetrahedronlist[i];
437: }
439: PetscCall(DMPlexInvertCells_CTetgen(numCells, numCorners, cells));
440: PetscCall(DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined));
441: if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscCall(PetscFree(meshCoords));
442: if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscCall(PetscFree(cells));
444: /* Set labels */
445: PetscCall(DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined));
446: for (v = 0; v < numVertices; ++v) {
447: if (out->pointmarkerlist[v]) PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]));
448: }
449: if (interpolate) {
450: for (PetscInt e = 0; e < out->numberofedges; e++) {
451: if (out->edgemarkerlist[e]) {
452: const PetscInt vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
453: const PetscInt *edges;
454: PetscInt numEdges;
456: PetscCall(DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges));
457: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %" PetscInt_FMT, numEdges);
458: PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]));
459: PetscCall(DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges));
460: }
461: }
462: for (PetscInt f = 0; f < out->numberoftrifaces; f++) {
463: if (out->trifacemarkerlist[f]) {
464: const PetscInt vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
465: const PetscInt *faces;
466: PetscInt numFaces;
468: PetscCall(DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces));
469: PetscCheck(numFaces == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %" PetscInt_FMT, numFaces);
470: PetscCall(DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]));
471: PetscCall(DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces));
472: }
473: }
474: }
476: #ifdef PETSC_HAVE_EGADS
477: {
478: DMLabel bodyLabel;
479: PetscContainer modelObj;
480: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
481: ego *bodies;
482: ego model, geom;
483: int Nb, oclass, mtype, *senses;
485: /* Get Attached EGADS Model from Original DMPlex */
486: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
487: if (modelObj) {
488: PetscCall(PetscContainerGetPointer(modelObj, &model));
489: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
490: /* Transfer EGADS Model to Volumetric Mesh */
491: PetscCall(PetscObjectCompose((PetscObject)*dmRefined, "EGADS Model", (PetscObject)modelObj));
493: /* Set Cell Labels */
494: PetscCall(DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel));
495: PetscCall(DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd));
496: PetscCall(DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd));
497: PetscCall(DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd));
499: for (c = cStart; c < cEnd; ++c) {
500: PetscReal centroid[3] = {0., 0., 0.};
501: PetscInt b;
503: /* Determine what body the cell's centroid is located in */
504: if (!interpolate) {
505: PetscSection coordSection;
506: Vec coordinates;
507: PetscScalar *coords = NULL;
508: PetscInt coordSize, s, d;
510: PetscCall(DMGetCoordinatesLocal(*dmRefined, &coordinates));
511: PetscCall(DMGetCoordinateSection(*dmRefined, &coordSection));
512: PetscCall(DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
513: for (s = 0; s < coordSize; ++s)
514: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
515: PetscCall(DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords));
516: } else PetscCall(DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL));
517: for (b = 0; b < Nb; ++b) {
518: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
519: }
520: if (b < Nb) {
521: PetscInt cval = b, eVal, fVal;
522: PetscInt *closure = NULL, Ncl, cl;
524: PetscCall(DMLabelSetValue(bodyLabel, c, cval));
525: PetscCall(DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
526: for (cl = 0; cl < Ncl; cl += 2) {
527: const PetscInt p = closure[cl];
529: if (p >= eStart && p < eEnd) {
530: PetscCall(DMLabelGetValue(bodyLabel, p, &eVal));
531: if (eVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
532: }
533: if (p >= fStart && p < fEnd) {
534: PetscCall(DMLabelGetValue(bodyLabel, p, &fVal));
535: if (fVal < 0) PetscCall(DMLabelSetValue(bodyLabel, p, cval));
536: }
537: }
538: PetscCall(DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure));
539: }
540: }
541: }
542: }
543: #endif
544: PetscCall(DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE));
545: }
546: PetscCall(DMUniversalLabelDestroy(&universal));
547: PetscCall(PLCDestroy(&in));
548: PetscCall(PLCDestroy(&out));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }