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) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
15: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
16: #undef SWAP
17: return 0;
18: }
20: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
21: {
22: MPI_Comm comm;
23: const PetscInt dim = 3;
24: PLC *in, *out;
25: DMUniversalLabel universal;
26: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, verbose = 0;
27: DMPlexInterpolatedFlag isInterpolated;
28: PetscMPIInt rank;
30: PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
31: PetscObjectGetComm((PetscObject)boundary,&comm);
32: MPI_Comm_rank(comm, &rank);
33: DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
34: DMUniversalLabelCreate(boundary, &universal);
36: PLCCreate(&in);
37: PLCCreate(&out);
39: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
40: in->numberofpoints = vEnd - vStart;
41: if (in->numberofpoints > 0) {
42: PetscSection coordSection;
43: Vec coordinates;
44: const PetscScalar *array;
46: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
47: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
48: DMGetCoordinatesLocal(boundary, &coordinates);
49: DMGetCoordinateSection(boundary, &coordSection);
50: VecGetArrayRead(coordinates, &array);
51: for (v = vStart; v < vEnd; ++v) {
52: const PetscInt idx = v - vStart;
53: PetscInt off, d, m;
55: PetscSectionGetOffset(coordSection, v, &off);
56: for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
57: DMLabelGetValue(universal->label, v, &m);
58: in->pointmarkerlist[idx] = (int) m;
59: }
60: VecRestoreArrayRead(coordinates, &array);
61: }
63: DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
64: in->numberofedges = eEnd - eStart;
65: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
66: PetscMalloc1(in->numberofedges*2, &in->edgelist);
67: PetscMalloc1(in->numberofedges, &in->edgemarkerlist);
68: for (e = eStart; e < eEnd; ++e) {
69: const PetscInt idx = e - eStart;
70: const PetscInt *cone;
71: PetscInt coneSize, val;
73: DMPlexGetConeSize(boundary, e, &coneSize);
74: DMPlexGetCone(boundary, e, &cone);
75: in->edgelist[idx*2] = cone[0] - vStart;
76: in->edgelist[idx*2 + 1] = cone[1] - vStart;
78: DMLabelGetValue(universal->label, e, &val);
79: in->edgemarkerlist[idx] = (int) val;
80: }
81: }
83: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
84: in->numberoffacets = fEnd - fStart;
85: if (in->numberoffacets > 0) {
86: PetscMalloc1(in->numberoffacets, &in->facetlist);
87: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
88: for (f = fStart; f < fEnd; ++f) {
89: const PetscInt idx = f - fStart;
90: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
91: polygon *poly;
93: in->facetlist[idx].numberofpolygons = 1;
94: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
95: in->facetlist[idx].numberofholes = 0;
96: in->facetlist[idx].holelist = NULL;
98: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
99: for (p = 0; p < numPoints*2; p += 2) {
100: const PetscInt point = points[p];
101: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
102: }
104: poly = in->facetlist[idx].polygonlist;
105: poly->numberofvertices = numVertices;
106: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
107: for (v = 0; v < numVertices; ++v) {
108: const PetscInt vIdx = points[v] - vStart;
109: poly->vertexlist[v] = vIdx;
110: }
111: DMLabelGetValue(universal->label, f, &m);
112: in->facetmarkerlist[idx] = (int) m;
113: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
114: }
115: }
116: if (rank == 0) {
117: TetGenOpts t;
119: TetGenOptsInitialize(&t);
120: t.in = boundary; /* Should go away */
121: t.plc = 1;
122: t.quality = 1;
123: t.edgesout = 1;
124: t.zeroindex = 1;
125: t.quiet = 1;
126: t.verbose = verbose;
127: #if 0
128: #ifdef PETSC_HAVE_EGADS
129: /* Need to add in more TetGen code */
130: t.nobisect = 1; /* Add Y to preserve Surface Mesh for EGADS */
131: #endif
132: #endif
134: TetGenCheckOpts(&t);
135: TetGenTetrahedralize(&t, in, out);
136: }
137: {
138: const PetscInt numCorners = 4;
139: const PetscInt numCells = out->numberoftetrahedra;
140: const PetscInt numVertices = out->numberofpoints;
141: PetscReal *meshCoords = NULL;
142: PetscInt *cells = NULL;
144: if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
145: meshCoords = (PetscReal *) out->pointlist;
146: } else {
147: PetscInt i;
149: PetscMalloc1(dim * numVertices, &meshCoords);
150: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
151: }
152: if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
153: cells = (PetscInt *) out->tetrahedronlist;
154: } else {
155: PetscInt i;
157: PetscMalloc1(numCells * numCorners, &cells);
158: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out->tetrahedronlist[i];
159: }
161: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
162: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
163: if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {
164: PetscFree(meshCoords);
165: }
166: if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {
167: PetscFree(cells);
168: }
170: /* Set labels */
171: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
172: for (v = 0; v < numVertices; ++v) {
173: if (out->pointmarkerlist[v]) {
174: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
175: }
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: /* Deterimine 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) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
244: DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
245: } else {
246: DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
247: }
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: DMLabelSetValue(bodyLabel, c, cval);
256: 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: DMLabelGetValue(bodyLabel, p, &eVal);
262: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
263: }
264: if (p >= fStart && p < fEnd) {
265: DMLabelGetValue(bodyLabel, p, &fVal);
266: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
267: }
268: }
269: DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);
270: }
271: }
272: }
273: }
274: #endif
275: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
276: }
278: DMUniversalLabelDestroy(&universal);
279: PLCDestroy(&in);
280: PLCDestroy(&out);
281: return 0;
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, verbose = 0;
291: DMPlexInterpolatedFlag isInterpolated;
292: PetscMPIInt rank;
294: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
295: PetscObjectGetComm((PetscObject)dm,&comm);
296: MPI_Comm_rank(comm, &rank);
297: DMPlexIsInterpolatedCollective(dm, &isInterpolated);
298: DMUniversalLabelCreate(dm, &universal);
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: PetscMalloc1(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: 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: 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: 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]) {
450: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
451: }
452: }
453: if (interpolate) {
454: PetscInt e, f;
456: for (e = 0; e < out->numberofedges; e++) {
457: if (out->edgemarkerlist[e]) {
458: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
459: const PetscInt *edges;
460: PetscInt numEdges;
462: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
464: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
465: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
466: }
467: }
468: for (f = 0; f < out->numberoftrifaces; f++) {
469: if (out->trifacemarkerlist[f]) {
470: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
471: const PetscInt *faces;
472: PetscInt numFaces;
474: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
476: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
477: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
478: }
479: }
480: }
482: #ifdef PETSC_HAVE_EGADS
483: {
484: DMLabel bodyLabel;
485: PetscContainer modelObj;
486: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
487: ego *bodies;
488: ego model, geom;
489: int Nb, oclass, mtype, *senses;
491: /* Get Attached EGADS Model from Original DMPlex */
492: PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
493: if (modelObj) {
494: PetscContainerGetPointer(modelObj, (void **) &model);
495: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
496: /* Transfer EGADS Model to Volumetric Mesh */
497: PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);
499: /* Set Cell Labels */
500: DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
501: DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
502: DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
503: DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);
505: for (c = cStart; c < cEnd; ++c) {
506: PetscReal centroid[3] = {0., 0., 0.};
507: PetscInt b;
509: /* Deterimine what body the cell's centroid is located in */
510: if (!interpolate) {
511: PetscSection coordSection;
512: Vec coordinates;
513: PetscScalar *coords = NULL;
514: PetscInt coordSize, s, d;
516: DMGetCoordinatesLocal(*dmRefined, &coordinates);
517: DMGetCoordinateSection(*dmRefined, &coordSection);
518: DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
519: for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
520: DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
521: } else {
522: DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
523: }
524: for (b = 0; b < Nb; ++b) {
525: if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break;
526: }
527: if (b < Nb) {
528: PetscInt cval = b, eVal, fVal;
529: PetscInt *closure = NULL, Ncl, cl;
531: DMLabelSetValue(bodyLabel, c, cval);
532: DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
533: for (cl = 0; cl < Ncl; cl += 2) {
534: const PetscInt p = closure[cl];
536: if (p >= eStart && p < eEnd) {
537: DMLabelGetValue(bodyLabel, p, &eVal);
538: if (eVal < 0) DMLabelSetValue(bodyLabel, p, cval);
539: }
540: if (p >= fStart && p < fEnd) {
541: DMLabelGetValue(bodyLabel, p, &fVal);
542: if (fVal < 0) DMLabelSetValue(bodyLabel, p, cval);
543: }
544: }
545: DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);
546: }
547: }
548: }
549: }
550: #endif
551: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
552: }
553: DMUniversalLabelDestroy(&universal);
554: PLCDestroy(&in);
555: PLCDestroy(&out);
556: return 0;
557: }