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