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