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