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: #include <tetgen.h>
15: /* This is to fix the tetrahedron orientation from TetGen */
16: static PetscErrorCode DMPlexInvertCells_Tetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[])
17: {
18: PetscInt bound = numCells * numCorners, coff;
20: #define SWAP(a, b) \
21: do { \
22: PetscInt tmp = (a); \
23: (a) = (b); \
24: (b) = tmp; \
25: } while (0)
26: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
27: #undef SWAP
28: return 0;
29: }
31: PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm)
32: {
33: MPI_Comm comm;
34: const PetscInt dim = 3;
35: ::tetgenio in;
36: ::tetgenio out;
37: PetscContainer modelObj;
38: DMUniversalLabel universal;
39: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, defVal;
40: DMPlexInterpolatedFlag isInterpolated;
41: PetscMPIInt rank;
43: PetscObjectGetComm((PetscObject)boundary, &comm);
44: MPI_Comm_rank(comm, &rank);
45: DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
46: DMUniversalLabelCreate(boundary, &universal);
47: DMLabelGetDefaultValue(universal->label, &defVal);
49: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
50: in.numberofpoints = vEnd - vStart;
51: if (in.numberofpoints > 0) {
52: PetscSection coordSection;
53: Vec coordinates;
54: const PetscScalar *array;
56: in.pointlist = new double[in.numberofpoints * dim];
57: in.pointmarkerlist = new int[in.numberofpoints];
59: PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints);
60: DMGetCoordinatesLocal(boundary, &coordinates);
61: DMGetCoordinateSection(boundary, &coordSection);
62: VecGetArrayRead(coordinates, &array);
63: for (v = vStart; v < vEnd; ++v) {
64: const PetscInt idx = v - vStart;
65: PetscInt off, d, val;
67: PetscSectionGetOffset(coordSection, v, &off);
68: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
69: DMLabelGetValue(universal->label, v, &val);
70: if (val != defVal) in.pointmarkerlist[idx] = (int)val;
71: }
72: VecRestoreArrayRead(coordinates, &array);
73: }
75: DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
76: in.numberofedges = eEnd - eStart;
77: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
78: in.edgelist = new int[in.numberofedges * 2];
79: in.edgemarkerlist = new int[in.numberofedges];
80: for (e = eStart; e < eEnd; ++e) {
81: const PetscInt idx = e - eStart;
82: const PetscInt *cone;
83: PetscInt coneSize, val;
85: DMPlexGetConeSize(boundary, e, &coneSize);
86: DMPlexGetCone(boundary, e, &cone);
87: in.edgelist[idx * 2] = cone[0] - vStart;
88: in.edgelist[idx * 2 + 1] = cone[1] - vStart;
90: DMLabelGetValue(universal->label, e, &val);
91: if (val != defVal) in.edgemarkerlist[idx] = (int)val;
92: }
93: }
95: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
96: in.numberoffacets = fEnd - fStart;
97: if (in.numberoffacets > 0) {
98: in.facetlist = new tetgenio::facet[in.numberoffacets];
99: in.facetmarkerlist = new int[in.numberoffacets];
100: for (f = fStart; f < fEnd; ++f) {
101: const PetscInt idx = f - fStart;
102: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, val = -1;
104: in.facetlist[idx].numberofpolygons = 1;
105: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
106: in.facetlist[idx].numberofholes = 0;
107: in.facetlist[idx].holelist = NULL;
109: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
110: for (p = 0; p < numPoints * 2; p += 2) {
111: const PetscInt point = points[p];
112: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
113: }
115: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
116: poly->numberofvertices = numVertices;
117: poly->vertexlist = new int[poly->numberofvertices];
118: for (v = 0; v < numVertices; ++v) {
119: const PetscInt vIdx = points[v] - vStart;
120: poly->vertexlist[v] = vIdx;
121: }
122: DMLabelGetValue(universal->label, f, &val);
123: if (val != defVal) in.facetmarkerlist[idx] = (int)val;
124: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
125: }
126: }
127: if (rank == 0) {
128: DM_Plex *mesh = (DM_Plex *)boundary->data;
129: char args[32];
131: /* Take away 'Q' for verbose output */
132: #ifdef PETSC_HAVE_EGADS
133: PetscStrcpy(args, "pqezQY");
134: #else
135: PetscStrcpy(args, "pqezQ");
136: #endif
137: if (mesh->tetgenOpts) {
138: ::tetrahedralize(mesh->tetgenOpts, &in, &out);
139: } else {
140: ::tetrahedralize(args, &in, &out);
141: }
142: }
143: {
144: const PetscInt numCorners = 4;
145: const PetscInt numCells = out.numberoftetrahedra;
146: const PetscInt numVertices = out.numberofpoints;
147: PetscReal *meshCoords = NULL;
148: PetscInt *cells = NULL;
150: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
151: meshCoords = (PetscReal *)out.pointlist;
152: } else {
153: PetscInt i;
155: meshCoords = new PetscReal[dim * numVertices];
156: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out.pointlist[i];
157: }
158: if (sizeof(PetscInt) == sizeof(out.tetrahedronlist[0])) {
159: cells = (PetscInt *)out.tetrahedronlist;
160: } else {
161: PetscInt i;
163: cells = new PetscInt[numCells * numCorners];
164: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out.tetrahedronlist[i];
165: }
167: DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
168: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
170: /* Set labels */
171: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
172: for (v = 0; v < numVertices; ++v) {
173: if (out.pointmarkerlist[v]) DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out.pointmarkerlist[v]);
174: }
175: if (interpolate) {
176: PetscInt e;
178: for (e = 0; e < out.numberofedges; e++) {
179: if (out.edgemarkerlist[e]) {
180: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
181: const PetscInt *edges;
182: PetscInt numEdges;
184: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
186: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
187: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
188: }
189: }
190: for (f = 0; f < out.numberoftrifaces; f++) {
191: if (out.trifacemarkerlist[f]) {
192: const PetscInt vertices[3] = {out.trifacelist[f * 3 + 0] + numCells, out.trifacelist[f * 3 + 1] + numCells, out.trifacelist[f * 3 + 2] + numCells};
193: const PetscInt *faces;
194: PetscInt numFaces;
196: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
198: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
199: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
200: }
201: }
202: }
204: PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj);
205: if (modelObj) {
206: #ifdef PETSC_HAVE_EGADS
207: DMLabel bodyLabel;
208: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
209: PetscBool islite = PETSC_FALSE;
210: ego *bodies;
211: ego model, geom;
212: int Nb, oclass, mtype, *senses;
214: /* Get Attached EGADS Model from Original DMPlex */
215: PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj);
216: if (modelObj) {
217: PetscContainerGetPointer(modelObj, (void **)&model);
218: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
219: /* Transfer EGADS Model to Volumetric Mesh */
220: PetscObjectCompose((PetscObject)*dm, "EGADS Model", (PetscObject)modelObj);
221: } else {
222: PetscObjectQuery((PetscObject)boundary, "EGADSLite Model", (PetscObject *)&modelObj);
223: if (modelObj) {
224: PetscContainerGetPointer(modelObj, (void **)&model);
225: EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
226: /* Transfer EGADS Model to Volumetric Mesh */
227: PetscObjectCompose((PetscObject)*dm, "EGADSLite Model", (PetscObject)modelObj);
228: islite = PETSC_TRUE;
229: }
230: }
231: if (!modelObj) goto skip_egads;
233: /* Set Cell Labels */
234: DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
235: DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
236: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
237: DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);
239: for (c = cStart; c < cEnd; ++c) {
240: PetscReal centroid[3] = {0., 0., 0.};
241: PetscInt b;
243: /* Determine what body the cell's centroid is located in */
244: if (!interpolate) {
245: PetscSection coordSection;
246: Vec coordinates;
247: PetscScalar *coords = NULL;
248: PetscInt coordSize, s, d;
250: DMGetCoordinatesLocal(*dm, &coordinates);
251: DMGetCoordinateSection(*dm, &coordSection);
252: DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
253: for (s = 0; s < coordSize; ++s)
254: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
255: DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
256: } else DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
257: for (b = 0; b < Nb; ++b) {
258: if (islite) {
259: if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
260: } else {
261: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
262: }
263: }
264: if (b < Nb) {
265: PetscInt cval = b, eVal, fVal;
266: PetscInt *closure = NULL, Ncl, cl;
268: DMLabelSetValue(bodyLabel, c, cval);
269: DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
270: for (cl = 0; cl < Ncl; cl += 2) {
271: const PetscInt p = closure[cl];
273: if (p >= eStart && p < eEnd) {
274: DMLabelGetValue(bodyLabel, p, &eVal);
275: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
276: }
277: if (p >= fStart && p < fEnd) {
278: DMLabelGetValue(bodyLabel, p, &fVal);
279: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
280: }
281: }
282: DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
283: }
284: }
285: skip_egads:;
286: #endif
287: }
288: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
289: }
290: DMUniversalLabelDestroy(&universal);
291: return 0;
292: }
294: PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined)
295: {
296: MPI_Comm comm;
297: const PetscInt dim = 3;
298: ::tetgenio in;
299: ::tetgenio out;
300: PetscContainer modelObj;
301: DMUniversalLabel universal;
302: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal;
303: DMPlexInterpolatedFlag isInterpolated;
304: PetscMPIInt rank;
306: PetscObjectGetComm((PetscObject)dm, &comm);
307: MPI_Comm_rank(comm, &rank);
308: DMPlexIsInterpolatedCollective(dm, &isInterpolated);
309: DMUniversalLabelCreate(dm, &universal);
310: DMLabelGetDefaultValue(universal->label, &defVal);
312: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
313: in.numberofpoints = vEnd - vStart;
314: if (in.numberofpoints > 0) {
315: PetscSection coordSection;
316: Vec coordinates;
317: PetscScalar *array;
319: in.pointlist = new double[in.numberofpoints * dim];
320: in.pointmarkerlist = new int[in.numberofpoints];
322: PetscArrayzero(in.pointmarkerlist, (size_t)in.numberofpoints);
323: DMGetCoordinatesLocal(dm, &coordinates);
324: DMGetCoordinateSection(dm, &coordSection);
325: VecGetArray(coordinates, &array);
326: for (v = vStart; v < vEnd; ++v) {
327: const PetscInt idx = v - vStart;
328: PetscInt off, d, val;
330: PetscSectionGetOffset(coordSection, v, &off);
331: for (d = 0; d < dim; ++d) in.pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
332: DMLabelGetValue(universal->label, v, &val);
333: if (val != defVal) in.pointmarkerlist[idx] = (int)val;
334: }
335: VecRestoreArray(coordinates, &array);
336: }
338: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
339: in.numberofedges = eEnd - eStart;
340: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberofedges > 0) {
341: in.edgelist = new int[in.numberofedges * 2];
342: in.edgemarkerlist = new int[in.numberofedges];
343: for (e = eStart; e < eEnd; ++e) {
344: const PetscInt idx = e - eStart;
345: const PetscInt *cone;
346: PetscInt coneSize, val;
348: DMPlexGetConeSize(dm, e, &coneSize);
349: DMPlexGetCone(dm, e, &cone);
350: in.edgelist[idx * 2] = cone[0] - vStart;
351: in.edgelist[idx * 2 + 1] = cone[1] - vStart;
353: DMLabelGetValue(universal->label, e, &val);
354: if (val != defVal) in.edgemarkerlist[idx] = (int)val;
355: }
356: }
358: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
359: in.numberoffacets = fEnd - fStart;
360: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in.numberoffacets > 0) {
361: in.facetlist = new tetgenio::facet[in.numberoffacets];
362: in.facetmarkerlist = new int[in.numberoffacets];
363: for (f = fStart; f < fEnd; ++f) {
364: const PetscInt idx = f - fStart;
365: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, val;
367: in.facetlist[idx].numberofpolygons = 1;
368: in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons];
369: in.facetlist[idx].numberofholes = 0;
370: in.facetlist[idx].holelist = NULL;
372: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
373: for (p = 0; p < numPoints * 2; p += 2) {
374: const PetscInt point = points[p];
375: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
376: }
378: tetgenio::polygon *poly = in.facetlist[idx].polygonlist;
379: poly->numberofvertices = numVertices;
380: poly->vertexlist = new int[poly->numberofvertices];
381: for (v = 0; v < numVertices; ++v) {
382: const PetscInt vIdx = points[v] - vStart;
383: poly->vertexlist[v] = vIdx;
384: }
386: DMLabelGetValue(universal->label, f, &val);
387: if (val != defVal) in.facetmarkerlist[idx] = (int)val;
389: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
390: }
391: }
393: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
394: in.numberofcorners = 4;
395: in.numberoftetrahedra = cEnd - cStart;
396: in.tetrahedronvolumelist = (double *)maxVolumes;
397: if (in.numberoftetrahedra > 0) {
398: in.tetrahedronlist = new int[in.numberoftetrahedra * in.numberofcorners];
399: for (c = cStart; c < cEnd; ++c) {
400: const PetscInt idx = c - cStart;
401: PetscInt *closure = NULL;
402: PetscInt closureSize;
404: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
406: for (v = 0; v < 4; ++v) in.tetrahedronlist[idx * in.numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
407: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
408: }
409: }
411: if (rank == 0) {
412: char args[32];
414: /* Take away 'Q' for verbose output */
415: PetscStrcpy(args, "qezQra");
416: ::tetrahedralize(args, &in, &out);
417: }
419: in.tetrahedronvolumelist = NULL;
420: {
421: const PetscInt numCorners = 4;
422: const PetscInt numCells = out.numberoftetrahedra;
423: const PetscInt numVertices = out.numberofpoints;
424: PetscReal *meshCoords = NULL;
425: PetscInt *cells = NULL;
426: PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
428: if (sizeof(PetscReal) == sizeof(out.pointlist[0])) {
429: meshCoords = (PetscReal *)out.pointlist;
430: } else {
431: PetscInt i;
433: meshCoords = new PetscReal[dim * numVertices];
434: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal)out.pointlist[i];
435: }
436: if (sizeof(PetscInt) == sizeof(out.tetrahedronlist[0])) {
437: cells = (PetscInt *)out.tetrahedronlist;
438: } else {
439: PetscInt i;
441: cells = new PetscInt[numCells * numCorners];
442: for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out.tetrahedronlist[i];
443: }
445: DMPlexInvertCells_Tetgen(numCells, numCorners, cells);
446: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
447: if (sizeof(PetscReal) != sizeof(out.pointlist[0])) delete[] meshCoords;
448: if (sizeof(PetscInt) != sizeof(out.tetrahedronlist[0])) delete[] cells;
450: /* Set labels */
451: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
452: for (v = 0; v < numVertices; ++v) {
453: if (out.pointmarkerlist[v]) DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out.pointmarkerlist[v]);
454: }
455: if (interpolate) {
456: PetscInt e, f;
458: for (e = 0; e < out.numberofedges; ++e) {
459: if (out.edgemarkerlist[e]) {
460: const PetscInt vertices[2] = {out.edgelist[e * 2 + 0] + numCells, out.edgelist[e * 2 + 1] + numCells};
461: const PetscInt *edges;
462: PetscInt numEdges;
464: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
466: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out.edgemarkerlist[e]);
467: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
468: }
469: }
470: for (f = 0; f < out.numberoftrifaces; ++f) {
471: if (out.trifacemarkerlist[f]) {
472: const PetscInt vertices[3] = {out.trifacelist[f * 3 + 0] + numCells, out.trifacelist[f * 3 + 1] + numCells, out.trifacelist[f * 3 + 2] + numCells};
473: const PetscInt *faces;
474: PetscInt numFaces;
476: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
478: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out.trifacemarkerlist[f]);
479: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
480: }
481: }
482: }
484: PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj);
485: if (modelObj) {
486: #ifdef PETSC_HAVE_EGADS
487: DMLabel bodyLabel;
488: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
489: PetscBool islite = PETSC_FALSE;
490: ego *bodies;
491: ego model, geom;
492: int Nb, oclass, mtype, *senses;
494: /* Get Attached EGADS Model from Original DMPlex */
495: PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj);
496: if (modelObj) {
497: PetscContainerGetPointer(modelObj, (void **)&model);
498: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
499: /* Transfer EGADS Model to Volumetric Mesh */
500: PetscObjectCompose((PetscObject)*dmRefined, "EGADS Model", (PetscObject)modelObj);
501: } else {
502: PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj);
503: if (modelObj) {
504: PetscContainerGetPointer(modelObj, (void **)&model);
505: EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
506: /* Transfer EGADS Model to Volumetric Mesh */
507: PetscObjectCompose((PetscObject)*dmRefined, "EGADSLite Model", (PetscObject)modelObj);
508: islite = PETSC_TRUE;
509: }
510: }
511: if (!modelObj) goto skip_egads;
513: /* Set Cell Labels */
514: DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
515: DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
516: DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
517: DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);
519: for (c = cStart; c < cEnd; ++c) {
520: PetscReal centroid[3] = {0., 0., 0.};
521: PetscInt b;
523: /* Determine what body the cell's centroid is located in */
524: if (!interpolate) {
525: PetscSection coordSection;
526: Vec coordinates;
527: PetscScalar *coords = NULL;
528: PetscInt coordSize, s, d;
530: DMGetCoordinatesLocal(*dmRefined, &coordinates);
531: DMGetCoordinateSection(*dmRefined, &coordSection);
532: DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
533: for (s = 0; s < coordSize; ++s)
534: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
535: DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
536: } else DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
537: for (b = 0; b < Nb; ++b) {
538: if (islite) {
539: if (EGlite_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
540: } else {
541: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
542: }
543: }
544: if (b < Nb) {
545: PetscInt cval = b, eVal, fVal;
546: PetscInt *closure = NULL, Ncl, cl;
548: DMLabelSetValue(bodyLabel, c, cval);
549: DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
550: for (cl = 0; cl < Ncl; cl += 2) {
551: const PetscInt p = closure[cl];
553: if (p >= eStart && p < eEnd) {
554: DMLabelGetValue(bodyLabel, p, &eVal);
555: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
556: }
557: if (p >= fStart && p < fEnd) {
558: DMLabelGetValue(bodyLabel, p, &fVal);
559: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
560: }
561: }
562: DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
563: }
564: }
565: skip_egads:;
566: #endif
567: }
568: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
569: }
570: return 0;
571: }