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