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