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: #define SWAP(a, b) \
15: do { \
16: PetscInt tmp = (a); \
17: (a) = (b); \
18: (b) = tmp; \
19: } while (0)
20: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff], cells[coff + 1]);
21: #undef SWAP
22: return 0;
23: }
25: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
26: {
27: MPI_Comm comm;
28: const PetscInt dim = 3;
29: PLC *in, *out;
30: DMUniversalLabel universal;
31: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, defVal, verbose = 0;
32: DMPlexInterpolatedFlag isInterpolated;
33: PetscMPIInt rank;
35: PetscOptionsGetInt(NULL, ((PetscObject)boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
36: PetscObjectGetComm((PetscObject)boundary, &comm);
37: MPI_Comm_rank(comm, &rank);
38: DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
39: DMUniversalLabelCreate(boundary, &universal);
40: DMLabelGetDefaultValue(universal->label, &defVal);
42: PLCCreate(&in);
43: PLCCreate(&out);
45: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
46: in->numberofpoints = vEnd - vStart;
47: if (in->numberofpoints > 0) {
48: PetscSection coordSection;
49: Vec coordinates;
50: const PetscScalar *array;
52: PetscMalloc1(in->numberofpoints * dim, &in->pointlist);
53: PetscCalloc1(in->numberofpoints, &in->pointmarkerlist);
54: DMGetCoordinatesLocal(boundary, &coordinates);
55: DMGetCoordinateSection(boundary, &coordSection);
56: VecGetArrayRead(coordinates, &array);
57: for (v = vStart; v < vEnd; ++v) {
58: const PetscInt idx = v - vStart;
59: PetscInt off, d, m;
61: PetscSectionGetOffset(coordSection, v, &off);
62: for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
63: DMLabelGetValue(universal->label, v, &m);
64: if (m != defVal) in->pointmarkerlist[idx] = (int)m;
65: }
66: VecRestoreArrayRead(coordinates, &array);
67: }
69: DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
70: in->numberofedges = eEnd - eStart;
71: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
72: PetscMalloc1(in->numberofedges * 2, &in->edgelist);
73: PetscMalloc1(in->numberofedges, &in->edgemarkerlist);
74: for (e = eStart; e < eEnd; ++e) {
75: const PetscInt idx = e - eStart;
76: const PetscInt *cone;
77: PetscInt coneSize, val;
79: DMPlexGetConeSize(boundary, e, &coneSize);
80: DMPlexGetCone(boundary, e, &cone);
81: in->edgelist[idx * 2] = cone[0] - vStart;
82: in->edgelist[idx * 2 + 1] = cone[1] - vStart;
84: DMLabelGetValue(universal->label, e, &val);
85: if (val != defVal) in->edgemarkerlist[idx] = (int)val;
86: }
87: }
89: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
90: in->numberoffacets = fEnd - fStart;
91: if (in->numberoffacets > 0) {
92: PetscMalloc1(in->numberoffacets, &in->facetlist);
93: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
94: for (f = fStart; f < fEnd; ++f) {
95: const PetscInt idx = f - fStart;
96: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
97: polygon *poly;
99: in->facetlist[idx].numberofpolygons = 1;
100: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
101: in->facetlist[idx].numberofholes = 0;
102: in->facetlist[idx].holelist = NULL;
104: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
105: for (p = 0; p < numPoints * 2; p += 2) {
106: const PetscInt point = points[p];
107: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
108: }
110: poly = in->facetlist[idx].polygonlist;
111: poly->numberofvertices = numVertices;
112: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
113: for (v = 0; v < numVertices; ++v) {
114: const PetscInt vIdx = points[v] - vStart;
115: poly->vertexlist[v] = vIdx;
116: }
117: DMLabelGetValue(universal->label, f, &m);
118: if (m != defVal) in->facetmarkerlist[idx] = (int)m;
119: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
120: }
121: }
122: if (rank == 0) {
123: TetGenOpts t;
125: TetGenOptsInitialize(&t);
126: t.in = boundary; /* Should go away */
127: t.plc = 1;
128: t.quality = 1;
129: t.edgesout = 1;
130: t.zeroindex = 1;
131: t.quiet = 1;
132: t.verbose = verbose;
133: #if 0
134: #ifdef PETSC_HAVE_EGADS
135: /* Need to add in more TetGen code */
136: t.nobisect = 1; /* Add Y to preserve Surface Mesh for EGADS */
137: #endif
138: #endif
140: TetGenCheckOpts(&t);
141: TetGenTetrahedralize(&t, in, out);
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: PetscMalloc1(dim * numVertices, &meshCoords);
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: PetscMalloc1(numCells * numCorners, &cells);
164: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt)out->tetrahedronlist[i];
165: }
167: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
168: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
169: if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscFree(meshCoords);
170: if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscFree(cells);
172: /* Set labels */
173: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
174: for (v = 0; v < numVertices; ++v) {
175: if (out->pointmarkerlist[v]) DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]);
176: }
177: if (interpolate) {
178: PetscInt e;
180: for (e = 0; e < out->numberofedges; e++) {
181: if (out->edgemarkerlist[e]) {
182: const PetscInt vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
183: const PetscInt *edges;
184: PetscInt numEdges;
186: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
188: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
189: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
190: }
191: }
192: for (f = 0; f < out->numberoftrifaces; f++) {
193: if (out->trifacemarkerlist[f]) {
194: const PetscInt vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
195: const PetscInt *faces;
196: PetscInt numFaces;
198: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
200: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
201: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
202: }
203: }
204: }
206: #ifdef PETSC_HAVE_EGADS
207: {
208: DMLabel bodyLabel;
209: PetscContainer modelObj;
210: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
211: ego *bodies;
212: ego model, geom;
213: int Nb, oclass, mtype, *senses;
215: /* Get Attached EGADS Model from Original DMPlex */
216: PetscObjectQuery((PetscObject)boundary, "EGADS Model", (PetscObject *)&modelObj);
217: if (modelObj) {
218: PetscContainerGetPointer(modelObj, (void **)&model);
219: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
220: /* Transfer EGADS Model to Volumetric Mesh */
221: PetscObjectCompose((PetscObject)*dm, "EGADS Model", (PetscObject)modelObj);
223: /* Set Cell Labels */
224: DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
225: DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
226: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
227: DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);
229: for (c = cStart; c < cEnd; ++c) {
230: PetscReal centroid[3] = {0., 0., 0.};
231: PetscInt b;
233: /* Determine what body the cell's centroid is located in */
234: if (!interpolate) {
235: PetscSection coordSection;
236: Vec coordinates;
237: PetscScalar *coords = NULL;
238: PetscInt coordSize, s, d;
240: DMGetCoordinatesLocal(*dm, &coordinates);
241: DMGetCoordinateSection(*dm, &coordSection);
242: DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
243: for (s = 0; s < coordSize; ++s)
244: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
245: DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
246: } else DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
247: for (b = 0; b < Nb; ++b) {
248: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
249: }
250: if (b < Nb) {
251: PetscInt cval = b, eVal, fVal;
252: PetscInt *closure = NULL, Ncl, cl;
254: DMLabelSetValue(bodyLabel, c, cval);
255: DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
256: for (cl = 0; cl < Ncl; ++cl) {
257: const PetscInt p = closure[cl * 2];
259: if (p >= eStart && p < eEnd) {
260: DMLabelGetValue(bodyLabel, p, &eVal);
261: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
262: }
263: if (p >= fStart && p < fEnd) {
264: DMLabelGetValue(bodyLabel, p, &fVal);
265: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
266: }
267: }
268: DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
269: }
270: }
271: }
272: }
273: #endif
274: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
275: }
277: DMUniversalLabelDestroy(&universal);
278: PLCDestroy(&in);
279: PLCDestroy(&out);
280: return 0;
281: }
283: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
284: {
285: MPI_Comm comm;
286: const PetscInt dim = 3;
287: PLC *in, *out;
288: DMUniversalLabel universal;
289: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, defVal, verbose = 0;
290: DMPlexInterpolatedFlag isInterpolated;
291: PetscMPIInt rank;
293: PetscOptionsGetInt(NULL, ((PetscObject)dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
294: PetscObjectGetComm((PetscObject)dm, &comm);
295: MPI_Comm_rank(comm, &rank);
296: DMPlexIsInterpolatedCollective(dm, &isInterpolated);
297: DMUniversalLabelCreate(dm, &universal);
298: DMLabelGetDefaultValue(universal->label, &defVal);
300: PLCCreate(&in);
301: PLCCreate(&out);
303: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
304: in->numberofpoints = vEnd - vStart;
305: if (in->numberofpoints > 0) {
306: PetscSection coordSection;
307: Vec coordinates;
308: PetscScalar *array;
310: PetscMalloc1(in->numberofpoints * dim, &in->pointlist);
311: PetscCalloc1(in->numberofpoints, &in->pointmarkerlist);
312: DMGetCoordinatesLocal(dm, &coordinates);
313: DMGetCoordinateSection(dm, &coordSection);
314: VecGetArray(coordinates, &array);
315: for (v = vStart; v < vEnd; ++v) {
316: const PetscInt idx = v - vStart;
317: PetscInt off, d, m;
319: PetscSectionGetOffset(coordSection, v, &off);
320: for (d = 0; d < dim; ++d) in->pointlist[idx * dim + d] = PetscRealPart(array[off + d]);
321: DMLabelGetValue(universal->label, v, &m);
322: if (m != defVal) in->pointmarkerlist[idx] = (int)m;
323: }
324: VecRestoreArray(coordinates, &array);
325: }
327: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
328: in->numberofedges = eEnd - eStart;
329: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
330: PetscMalloc1(in->numberofedges * 2, &in->edgelist);
331: PetscMalloc1(in->numberofedges, &in->edgemarkerlist);
332: for (e = eStart; e < eEnd; ++e) {
333: const PetscInt idx = e - eStart;
334: const PetscInt *cone;
335: PetscInt coneSize, val;
337: DMPlexGetConeSize(dm, e, &coneSize);
338: DMPlexGetCone(dm, e, &cone);
339: in->edgelist[idx * 2] = cone[0] - vStart;
340: in->edgelist[idx * 2 + 1] = cone[1] - vStart;
342: DMLabelGetValue(universal->label, e, &val);
343: if (val != defVal) in->edgemarkerlist[idx] = (int)val;
344: }
345: }
347: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
348: in->numberoftrifaces = 0;
349: for (f = fStart; f < fEnd; ++f) {
350: PetscInt supportSize;
352: DMPlexGetSupportSize(dm, f, &supportSize);
353: if (supportSize == 1) ++in->numberoftrifaces;
354: }
355: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
356: PetscInt tf = 0;
358: PetscMalloc1(in->numberoftrifaces * 3, &in->trifacelist);
359: PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist);
360: for (f = fStart; f < fEnd; ++f) {
361: PetscInt *points = NULL;
362: PetscInt supportSize, numPoints, p, Nv = 0, val;
364: DMPlexGetSupportSize(dm, f, &supportSize);
365: if (supportSize != 1) continue;
366: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
367: for (p = 0; p < numPoints * 2; p += 2) {
368: const PetscInt point = points[p];
369: if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf * 3 + Nv++] = point - vStart;
370: }
371: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
373: DMLabelGetValue(universal->label, f, &val);
374: if (val != defVal) in->trifacemarkerlist[tf] = (int)val;
375: ++tf;
376: }
377: }
379: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
380: in->numberofcorners = 4;
381: in->numberoftetrahedra = cEnd - cStart;
382: in->tetrahedronvolumelist = maxVolumes;
383: if (in->numberoftetrahedra > 0) {
384: PetscMalloc1(in->numberoftetrahedra * in->numberofcorners, &in->tetrahedronlist);
385: for (c = cStart; c < cEnd; ++c) {
386: const PetscInt idx = c - cStart;
387: PetscInt *closure = NULL;
388: PetscInt closureSize;
390: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
392: for (v = 0; v < 4; ++v) in->tetrahedronlist[idx * in->numberofcorners + v] = closure[(v + closureSize - 4) * 2] - vStart;
393: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
394: }
395: }
397: if (rank == 0) {
398: TetGenOpts t;
400: TetGenOptsInitialize(&t);
402: t.in = dm; /* Should go away */
403: t.refine = 1;
404: t.varvolume = 1;
405: t.quality = 1;
406: t.edgesout = 1;
407: t.zeroindex = 1;
408: t.quiet = 1;
409: t.verbose = verbose; /* Change this */
411: TetGenCheckOpts(&t);
412: TetGenTetrahedralize(&t, in, out);
413: }
415: in->tetrahedronvolumelist = NULL;
416: {
417: const PetscInt numCorners = 4;
418: const PetscInt numCells = out->numberoftetrahedra;
419: const PetscInt numVertices = out->numberofpoints;
420: PetscReal *meshCoords = NULL;
421: PetscInt *cells = NULL;
422: PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
424: if (sizeof(PetscReal) == sizeof(out->pointlist[0])) {
425: meshCoords = (PetscReal *)out->pointlist;
426: } else {
427: PetscInt i;
429: PetscMalloc1(dim * numVertices, &meshCoords);
430: for (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: PetscInt i;
437: PetscMalloc1(numCells * numCorners, &cells);
438: for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt)out->tetrahedronlist[i];
439: }
441: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
442: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
443: if (sizeof(PetscReal) != sizeof(out->pointlist[0])) PetscFree(meshCoords);
444: if (sizeof(PetscInt) != sizeof(out->tetrahedronlist[0])) PetscFree(cells);
446: /* Set labels */
447: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
448: for (v = 0; v < numVertices; ++v) {
449: if (out->pointmarkerlist[v]) DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v + numCells, out->pointmarkerlist[v]);
450: }
451: if (interpolate) {
452: PetscInt e, f;
454: for (e = 0; e < out->numberofedges; e++) {
455: if (out->edgemarkerlist[e]) {
456: const PetscInt vertices[2] = {out->edgelist[e * 2 + 0] + numCells, out->edgelist[e * 2 + 1] + numCells};
457: const PetscInt *edges;
458: PetscInt numEdges;
460: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
462: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
463: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
464: }
465: }
466: for (f = 0; f < out->numberoftrifaces; f++) {
467: if (out->trifacemarkerlist[f]) {
468: const PetscInt vertices[3] = {out->trifacelist[f * 3 + 0] + numCells, out->trifacelist[f * 3 + 1] + numCells, out->trifacelist[f * 3 + 2] + numCells};
469: const PetscInt *faces;
470: PetscInt numFaces;
472: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
474: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
475: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
476: }
477: }
478: }
480: #ifdef PETSC_HAVE_EGADS
481: {
482: DMLabel bodyLabel;
483: PetscContainer modelObj;
484: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
485: ego *bodies;
486: ego model, geom;
487: int Nb, oclass, mtype, *senses;
489: /* Get Attached EGADS Model from Original DMPlex */
490: PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj);
491: if (modelObj) {
492: PetscContainerGetPointer(modelObj, (void **)&model);
493: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
494: /* Transfer EGADS Model to Volumetric Mesh */
495: PetscObjectCompose((PetscObject)*dmRefined, "EGADS Model", (PetscObject)modelObj);
497: /* Set Cell Labels */
498: DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
499: DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
500: DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
501: DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);
503: for (c = cStart; c < cEnd; ++c) {
504: PetscReal centroid[3] = {0., 0., 0.};
505: PetscInt b;
507: /* Determine what body the cell's centroid is located in */
508: if (!interpolate) {
509: PetscSection coordSection;
510: Vec coordinates;
511: PetscScalar *coords = NULL;
512: PetscInt coordSize, s, d;
514: DMGetCoordinatesLocal(*dmRefined, &coordinates);
515: DMGetCoordinateSection(*dmRefined, &coordSection);
516: DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
517: for (s = 0; s < coordSize; ++s)
518: for (d = 0; d < dim; ++d) centroid[d] += coords[s * dim + d];
519: DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
520: } else DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
521: for (b = 0; b < Nb; ++b) {
522: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
523: }
524: if (b < Nb) {
525: PetscInt cval = b, eVal, fVal;
526: PetscInt *closure = NULL, Ncl, cl;
528: DMLabelSetValue(bodyLabel, c, cval);
529: DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
530: for (cl = 0; cl < Ncl; cl += 2) {
531: const PetscInt p = closure[cl];
533: if (p >= eStart && p < eEnd) {
534: DMLabelGetValue(bodyLabel, p, &eVal);
535: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
536: }
537: if (p >= fStart && p < fEnd) {
538: DMLabelGetValue(bodyLabel, p, &fVal);
539: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
540: }
541: }
542: DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
543: }
544: }
545: }
546: }
547: #endif
548: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
549: }
550: DMUniversalLabelDestroy(&universal);
551: PLCDestroy(&in);
552: PLCDestroy(&out);
553: return 0;
554: }