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 */ 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 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, (void **)&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, (void **)&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 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, (void **)&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, (void **)&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