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