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