Actual source code: ctetgenerate.c
petsc-master 2021-02-25
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;
15: #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0)
16: for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]);
17: #undef SWAP
18: return(0);
19: }
21: PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm)
22: {
23: MPI_Comm comm;
24: const PetscInt dim = 3;
25: PLC *in, *out;
26: DMUniversalLabel universal;
27: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, verbose = 0;
28: DMPlexInterpolatedFlag isInterpolated;
29: PetscMPIInt rank;
30: PetscErrorCode ierr;
33: PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);
34: PetscObjectGetComm((PetscObject)boundary,&comm);
35: MPI_Comm_rank(comm, &rank);
36: DMPlexIsInterpolatedCollective(boundary, &isInterpolated);
37: DMUniversalLabelCreate(boundary, &universal);
39: PLCCreate(&in);
40: PLCCreate(&out);
42: DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);
43: in->numberofpoints = vEnd - vStart;
44: if (in->numberofpoints > 0) {
45: PetscSection coordSection;
46: Vec coordinates;
47: const PetscScalar *array;
49: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
50: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
51: DMGetCoordinatesLocal(boundary, &coordinates);
52: DMGetCoordinateSection(boundary, &coordSection);
53: VecGetArrayRead(coordinates, &array);
54: for (v = vStart; v < vEnd; ++v) {
55: const PetscInt idx = v - vStart;
56: PetscInt off, d, m;
58: PetscSectionGetOffset(coordSection, v, &off);
59: for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
60: DMLabelGetValue(universal->label, v, &m);
61: in->pointmarkerlist[idx] = (int) m;
62: }
63: VecRestoreArrayRead(coordinates, &array);
64: }
66: DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);
67: in->numberofedges = eEnd - eStart;
68: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
69: PetscMalloc1(in->numberofedges*2, &in->edgelist);
70: PetscMalloc1(in->numberofedges, &in->edgemarkerlist);
71: for (e = eStart; e < eEnd; ++e) {
72: const PetscInt idx = e - eStart;
73: const PetscInt *cone;
74: PetscInt coneSize, val;
76: DMPlexGetConeSize(boundary, e, &coneSize);
77: DMPlexGetCone(boundary, e, &cone);
78: in->edgelist[idx*2] = cone[0] - vStart;
79: in->edgelist[idx*2 + 1] = cone[1] - vStart;
81: DMLabelGetValue(universal->label, e, &val);
82: in->edgemarkerlist[idx] = (int) val;
83: }
84: }
86: DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);
87: in->numberoffacets = fEnd - fStart;
88: if (in->numberoffacets > 0) {
89: PetscMalloc1(in->numberoffacets, &in->facetlist);
90: PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);
91: for (f = fStart; f < fEnd; ++f) {
92: const PetscInt idx = f - fStart;
93: PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1;
94: polygon *poly;
96: in->facetlist[idx].numberofpolygons = 1;
97: PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);
98: in->facetlist[idx].numberofholes = 0;
99: in->facetlist[idx].holelist = NULL;
101: DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
102: for (p = 0; p < numPoints*2; p += 2) {
103: const PetscInt point = points[p];
104: if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point;
105: }
107: poly = in->facetlist[idx].polygonlist;
108: poly->numberofvertices = numVertices;
109: PetscMalloc1(poly->numberofvertices, &poly->vertexlist);
110: for (v = 0; v < numVertices; ++v) {
111: const PetscInt vIdx = points[v] - vStart;
112: poly->vertexlist[v] = vIdx;
113: }
114: DMLabelGetValue(universal->label, f, &m);
115: in->facetmarkerlist[idx] = (int) m;
116: DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);
117: }
118: }
119: if (!rank) {
120: TetGenOpts t;
122: TetGenOptsInitialize(&t);
123: t.in = boundary; /* Should go away */
124: t.plc = 1;
125: t.quality = 1;
126: t.edgesout = 1;
127: t.zeroindex = 1;
128: t.quiet = 1;
129: t.verbose = verbose;
130: #if 0
131: #ifdef PETSC_HAVE_EGADS
132: /* Need to add in more TetGen code */
133: t.nobisect = 1; /* Add Y to preserve Surface Mesh for EGADS */
134: #endif
135: #endif
137: TetGenCheckOpts(&t);
138: TetGenTetrahedralize(&t, in, out);
139: }
140: {
141: const PetscInt numCorners = 4;
142: const PetscInt numCells = out->numberoftetrahedra;
143: const PetscInt numVertices = out->numberofpoints;
144: PetscReal *meshCoords = NULL;
145: PetscInt *cells = NULL;
147: if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
148: meshCoords = (PetscReal *) out->pointlist;
149: } else {
150: PetscInt i;
152: PetscMalloc1(dim * numVertices, &meshCoords);
153: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
154: }
155: if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
156: cells = (PetscInt *) out->tetrahedronlist;
157: } else {
158: PetscInt i;
160: PetscMalloc1(numCells * numCorners, &cells);
161: for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out->tetrahedronlist[i];
162: }
164: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
165: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);
166: if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {
167: PetscFree(meshCoords);
168: }
169: if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {
170: PetscFree(cells);
171: }
173: /* Set labels */
174: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);
175: for (v = 0; v < numVertices; ++v) {
176: if (out->pointmarkerlist[v]) {
177: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
178: }
179: }
180: if (interpolate) {
181: PetscInt e;
183: for (e = 0; e < out->numberofedges; e++) {
184: if (out->edgemarkerlist[e]) {
185: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
186: const PetscInt *edges;
187: PetscInt numEdges;
189: DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);
190: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
191: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
192: DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);
193: }
194: }
195: for (f = 0; f < out->numberoftrifaces; f++) {
196: if (out->trifacemarkerlist[f]) {
197: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
198: const PetscInt *faces;
199: PetscInt numFaces;
201: DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);
202: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
203: DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
204: DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);
205: }
206: }
207: }
209: #ifdef PETSC_HAVE_EGADS
210: {
211: DMLabel bodyLabel;
212: PetscContainer modelObj;
213: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
214: ego *bodies;
215: ego model, geom;
216: int Nb, oclass, mtype, *senses;
218: /* Get Attached EGADS Model from Original DMPlex */
219: PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);
220: if (modelObj) {
221: PetscContainerGetPointer(modelObj, (void **) &model);
222: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
223: /* Transfer EGADS Model to Volumetric Mesh */
224: PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);
226: /* Set Cell Labels */
227: DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);
228: DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);
229: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
230: DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);
232: for (c = cStart; c < cEnd; ++c) {
233: PetscReal centroid[3] = {0., 0., 0.};
234: PetscInt b;
236: /* Deterimine what body the cell's centroid is located in */
237: if (!interpolate) {
238: PetscSection coordSection;
239: Vec coordinates;
240: PetscScalar *coords = NULL;
241: PetscInt coordSize, s, d;
243: DMGetCoordinatesLocal(*dm, &coordinates);
244: DMGetCoordinateSection(*dm, &coordSection);
245: DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
246: for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
247: DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);
248: } else {
249: DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);
250: }
251: for (b = 0; b < Nb; ++b) {
252: 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: }
276: }
277: #endif
278: DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
279: }
281: DMUniversalLabelDestroy(&universal);
282: PLCDestroy(&in);
283: PLCDestroy(&out);
284: return(0);
285: }
287: PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined)
288: {
289: MPI_Comm comm;
290: const PetscInt dim = 3;
291: PLC *in, *out;
292: DMUniversalLabel universal;
293: PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, verbose = 0;
294: DMPlexInterpolatedFlag isInterpolated;
295: PetscMPIInt rank;
296: PetscErrorCode ierr;
299: PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);
300: PetscObjectGetComm((PetscObject)dm,&comm);
301: MPI_Comm_rank(comm, &rank);
302: DMPlexIsInterpolatedCollective(dm, &isInterpolated);
303: DMUniversalLabelCreate(dm, &universal);
305: PLCCreate(&in);
306: PLCCreate(&out);
308: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
309: in->numberofpoints = vEnd - vStart;
310: if (in->numberofpoints > 0) {
311: PetscSection coordSection;
312: Vec coordinates;
313: PetscScalar *array;
315: PetscMalloc1(in->numberofpoints*dim, &in->pointlist);
316: PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);
317: DMGetCoordinatesLocal(dm, &coordinates);
318: DMGetCoordinateSection(dm, &coordSection);
319: VecGetArray(coordinates, &array);
320: for (v = vStart; v < vEnd; ++v) {
321: const PetscInt idx = v - vStart;
322: PetscInt off, d, m;
324: PetscSectionGetOffset(coordSection, v, &off);
325: for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]);
326: DMLabelGetValue(universal->label, v, &m);
327: in->pointmarkerlist[idx] = (int) m;
328: }
329: VecRestoreArray(coordinates, &array);
330: }
332: DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);
333: in->numberofedges = eEnd - eStart;
334: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) {
335: PetscMalloc1(in->numberofedges * 2, &in->edgelist);
336: PetscMalloc1(in->numberofedges, &in->edgemarkerlist);
337: for (e = eStart; e < eEnd; ++e) {
338: const PetscInt idx = e - eStart;
339: const PetscInt *cone;
340: PetscInt coneSize, val;
342: DMPlexGetConeSize(dm, e, &coneSize);
343: DMPlexGetCone(dm, e, &cone);
344: in->edgelist[idx*2] = cone[0] - vStart;
345: in->edgelist[idx*2 + 1] = cone[1] - vStart;
347: DMLabelGetValue(universal->label, e, &val);
348: in->edgemarkerlist[idx] = (int) val;
349: }
350: }
352: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
353: in->numberoftrifaces = 0;
354: for (f = fStart; f < fEnd; ++f) {
355: PetscInt supportSize;
357: DMPlexGetSupportSize(dm, f, &supportSize);
358: if (supportSize == 1) ++in->numberoftrifaces;
359: }
360: if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) {
361: PetscInt tf = 0;
363: PetscMalloc1(in->numberoftrifaces*3, &in->trifacelist);
364: PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist);
365: for (f = fStart; f < fEnd; ++f) {
366: PetscInt *points = NULL;
367: PetscInt supportSize, numPoints, p, Nv = 0, val;
369: DMPlexGetSupportSize(dm, f, &supportSize);
370: if (supportSize != 1) continue;
371: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
372: for (p = 0; p < numPoints*2; p += 2) {
373: const PetscInt point = points[p];
374: if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf*3 + Nv++] = point - vStart;
375: }
376: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);
377: if (Nv != 3) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %D has %D vertices, not 3", f, Nv);
378: DMLabelGetValue(universal->label, f, &val);
379: in->trifacemarkerlist[tf] = (int) val;
380: ++tf;
381: }
382: }
384: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
385: in->numberofcorners = 4;
386: in->numberoftetrahedra = cEnd - cStart;
387: in->tetrahedronvolumelist = maxVolumes;
388: if (in->numberoftetrahedra > 0) {
389: PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);
390: for (c = cStart; c < cEnd; ++c) {
391: const PetscInt idx = c - cStart;
392: PetscInt *closure = NULL;
393: PetscInt closureSize;
395: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
396: if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize);
397: for (v = 0; v < 4; ++v) in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart;
398: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);
399: }
400: }
402: if (!rank) {
403: TetGenOpts t;
405: TetGenOptsInitialize(&t);
407: t.in = dm; /* Should go away */
408: t.refine = 1;
409: t.varvolume = 1;
410: t.quality = 1;
411: t.edgesout = 1;
412: t.zeroindex = 1;
413: t.quiet = 1;
414: t.verbose = verbose; /* Change this */
416: TetGenCheckOpts(&t);
417: TetGenTetrahedralize(&t, in, out);
418: }
420: in->tetrahedronvolumelist = NULL;
421: {
422: const PetscInt numCorners = 4;
423: const PetscInt numCells = out->numberoftetrahedra;
424: const PetscInt numVertices = out->numberofpoints;
425: PetscReal *meshCoords = NULL;
426: PetscInt *cells = NULL;
427: PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE;
429: if (sizeof (PetscReal) == sizeof (out->pointlist[0])) {
430: meshCoords = (PetscReal *) out->pointlist;
431: } else {
432: PetscInt i;
434: PetscMalloc1(dim * numVertices, &meshCoords);
435: for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i];
436: }
437: if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) {
438: cells = (PetscInt *) out->tetrahedronlist;
439: } else {
440: PetscInt i;
442: PetscMalloc1(numCells * numCorners, &cells);
443: for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt) out->tetrahedronlist[i];
444: }
446: DMPlexInvertCells_CTetgen(numCells, numCorners, cells);
447: DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);
448: if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {PetscFree(meshCoords);}
449: if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {PetscFree(cells);}
451: /* Set labels */
452: DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);
453: for (v = 0; v < numVertices; ++v) {
454: if (out->pointmarkerlist[v]) {
455: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);
456: }
457: }
458: if (interpolate) {
459: PetscInt e, f;
461: for (e = 0; e < out->numberofedges; e++) {
462: if (out->edgemarkerlist[e]) {
463: const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells};
464: const PetscInt *edges;
465: PetscInt numEdges;
467: DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);
468: if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges);
469: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);
470: DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);
471: }
472: }
473: for (f = 0; f < out->numberoftrifaces; f++) {
474: if (out->trifacemarkerlist[f]) {
475: const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells};
476: const PetscInt *faces;
477: PetscInt numFaces;
479: DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);
480: if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces);
481: DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);
482: DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);
483: }
484: }
485: }
487: #ifdef PETSC_HAVE_EGADS
488: {
489: DMLabel bodyLabel;
490: PetscContainer modelObj;
491: PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd;
492: ego *bodies;
493: ego model, geom;
494: int Nb, oclass, mtype, *senses;
496: /* Get Attached EGADS Model from Original DMPlex */
497: PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
498: if (modelObj) {
499: PetscContainerGetPointer(modelObj, (void **) &model);
500: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
501: /* Transfer EGADS Model to Volumetric Mesh */
502: PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);
504: /* Set Cell Labels */
505: DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);
506: DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);
507: DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);
508: DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);
510: for (c = cStart; c < cEnd; ++c) {
511: PetscReal centroid[3] = {0., 0., 0.};
512: PetscInt b;
514: /* Deterimine what body the cell's centroid is located in */
515: if (!interpolate) {
516: PetscSection coordSection;
517: Vec coordinates;
518: PetscScalar *coords = NULL;
519: PetscInt coordSize, s, d;
521: DMGetCoordinatesLocal(*dmRefined, &coordinates);
522: DMGetCoordinateSection(*dmRefined, &coordSection);
523: DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
524: for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d];
525: DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);
526: } else {
527: DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);
528: }
529: for (b = 0; b < Nb; ++b) {
530: 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: }
554: }
555: #endif
556: DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);
557: }
558: DMUniversalLabelDestroy(&universal);
559: PLCDestroy(&in);
560: PLCDestroy(&out);
561: return(0);
562: }