1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petsc/private/hashmapi.h> 3 4 #ifdef PETSC_HAVE_EGADS 5 #include <egads.h> 6 #endif 7 8 /* We need to understand how to natively parse STEP files. There seems to be only one open source implementation of 9 the STEP parser contained in the OpenCASCADE package. It is enough to make a strong man weep: 10 11 https://github.com/tpaviot/oce/tree/master/src/STEPControl 12 13 The STEP, and inner EXPRESS, formats are ISO standards, so they are documented 14 15 https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files 16 http://stepmod.sourceforge.net/express_model_spec/ 17 18 but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants. 19 */ 20 21 #ifdef PETSC_HAVE_EGADS 22 PETSC_INTERN PetscErrorCode DMPlexSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]); 23 PETSC_INTERN PetscErrorCode DMPlexSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]); 24 25 PetscErrorCode DMPlexSnapToGeomModel_EGADS_Internal(DM dm, PetscInt p, ego model, PetscInt bodyID, PetscInt faceID, PetscInt edgeID, const PetscScalar mcoords[], PetscScalar gcoords[]) 26 { 27 DM cdm; 28 ego *bodies; 29 ego geom, body, obj; 30 /* result has to hold derivatives, along with the value */ 31 double params[3], result[18], paramsV[16 * 3], resultV[16 * 3], range[4]; 32 int Nb, oclass, mtype, *senses, peri; 33 Vec coordinatesLocal; 34 PetscScalar *coords = NULL; 35 PetscInt Nv, v, Np = 0, pm; 36 PetscInt dE, d; 37 38 PetscFunctionBeginHot; 39 PetscCall(DMGetCoordinateDM(dm, &cdm)); 40 PetscCall(DMGetCoordinateDim(dm, &dE)); 41 PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 42 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 43 PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb); 44 body = bodies[bodyID]; 45 46 if (edgeID >= 0) { 47 PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &obj)); 48 Np = 1; 49 } else if (faceID >= 0) { 50 PetscCall(EG_objectBodyTopo(body, FACE, faceID, &obj)); 51 Np = 2; 52 } else { 53 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 54 PetscFunctionReturn(0); 55 } 56 57 /* Calculate parameters (t or u,v) for vertices */ 58 PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); 59 Nv /= dE; 60 if (Nv == 1) { 61 PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); 62 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 63 PetscFunctionReturn(0); 64 } 65 PetscCheck(Nv <= 16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %" PetscInt_FMT " coordinates associated to point %" PetscInt_FMT, Nv, p); 66 67 /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */ 68 PetscCall(EG_getRange(obj, range, &peri)); 69 for (v = 0; v < Nv; ++v) { 70 PetscCall(EG_invEvaluate(obj, &coords[v * dE], ¶msV[v * 3], &resultV[v * 3])); 71 #if 1 72 if (peri > 0) { 73 if (paramsV[v * 3 + 0] + 1.e-4 < range[0]) { 74 paramsV[v * 3 + 0] += 2. * PETSC_PI; 75 } else if (paramsV[v * 3 + 0] - 1.e-4 > range[1]) { 76 paramsV[v * 3 + 0] -= 2. * PETSC_PI; 77 } 78 } 79 if (peri > 1) { 80 if (paramsV[v * 3 + 1] + 1.e-4 < range[2]) { 81 paramsV[v * 3 + 1] += 2. * PETSC_PI; 82 } else if (paramsV[v * 3 + 1] - 1.e-4 > range[3]) { 83 paramsV[v * 3 + 1] -= 2. * PETSC_PI; 84 } 85 } 86 #endif 87 } 88 PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); 89 /* Calculate parameters (t or u,v) for new vertex at edge midpoint */ 90 for (pm = 0; pm < Np; ++pm) { 91 params[pm] = 0.; 92 for (v = 0; v < Nv; ++v) params[pm] += paramsV[v * 3 + pm]; 93 params[pm] /= Nv; 94 } 95 PetscCheck(!(params[0] < range[0]) && !(params[0] > range[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p); 96 PetscCheck(Np <= 1 || (params[1] >= range[2] && params[1] <= range[3]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p); 97 /* Put coordinates for new vertex in result[] */ 98 PetscCall(EG_evaluate(obj, params, result)); 99 for (d = 0; d < dE; ++d) gcoords[d] = result[d]; 100 PetscFunctionReturn(0); 101 } 102 #endif 103 104 /*@ 105 DMPlexSnapToGeomModel - Given a coordinate point 'mcoords' on the mesh point 'p', return the closest coordinate point 'gcoords' on the geometry model associated with that point. 106 107 Not collective 108 109 Input Parameters: 110 + dm - The `DMPLEX` object 111 . p - The mesh point 112 . dE - The coordinate dimension 113 - mcoords - A coordinate point lying on the mesh point 114 115 Output Parameter: 116 . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point 117 118 Level: intermediate 119 120 Note: 121 Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS. The coordinate dimension may be different from the coordinate dimension of the dm, for example if the transformation is extrusion. 122 123 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMPlexSetRefinementUniform()` 124 @*/ 125 PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[]) 126 { 127 PetscInt d; 128 129 PetscFunctionBeginHot; 130 #ifdef PETSC_HAVE_EGADS 131 { 132 DM_Plex *plex = (DM_Plex *)dm->data; 133 DMLabel bodyLabel, faceLabel, edgeLabel; 134 PetscInt bodyID, faceID, edgeID; 135 PetscContainer modelObj; 136 ego model; 137 PetscBool islite = PETSC_FALSE; 138 139 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 140 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 141 PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 142 if (!bodyLabel || !faceLabel || !edgeLabel || plex->ignoreModel) { 143 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 144 PetscFunctionReturn(0); 145 } 146 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj)); 147 if (!modelObj) { 148 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj)); 149 islite = PETSC_TRUE; 150 } 151 if (!modelObj) { 152 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 153 PetscFunctionReturn(0); 154 } 155 PetscCall(PetscContainerGetPointer(modelObj, (void **)&model)); 156 PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID)); 157 PetscCall(DMLabelGetValue(faceLabel, p, &faceID)); 158 PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID)); 159 /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */ 160 if (bodyID < 0) { 161 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 162 PetscFunctionReturn(0); 163 } 164 if (islite) PetscCall(DMPlexSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords)); 165 else PetscCall(DMPlexSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords)); 166 } 167 #else 168 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 169 #endif 170 PetscFunctionReturn(0); 171 } 172 173 #if defined(PETSC_HAVE_EGADS) 174 static PetscErrorCode DMPlexEGADSPrintModel_Internal(ego model) 175 { 176 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 177 int oclass, mtype, *senses; 178 int Nb, b; 179 180 PetscFunctionBeginUser; 181 /* test bodyTopo functions */ 182 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 183 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb)); 184 185 for (b = 0; b < Nb; ++b) { 186 ego body = bodies[b]; 187 int id, Nsh, Nf, Nl, l, Ne, e, Nv, v; 188 189 /* Output Basic Model Topology */ 190 PetscCall(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs)); 191 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh)); 192 EG_free(objs); 193 194 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs)); 195 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf)); 196 EG_free(objs); 197 198 PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 199 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl)); 200 201 PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs)); 202 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne)); 203 EG_free(objs); 204 205 PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs)); 206 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv)); 207 EG_free(objs); 208 209 for (l = 0; l < Nl; ++l) { 210 ego loop = lobjs[l]; 211 212 id = EG_indexBodyTopo(body, loop); 213 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id)); 214 215 /* Get EDGE info which associated with the current LOOP */ 216 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); 217 218 for (e = 0; e < Ne; ++e) { 219 ego edge = objs[e]; 220 double range[4] = {0., 0., 0., 0.}; 221 double point[3] = {0., 0., 0.}; 222 double params[3] = {0., 0., 0.}; 223 double result[18]; 224 int peri; 225 226 id = EG_indexBodyTopo(body, edge); 227 PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e)); 228 229 PetscCall(EG_getRange(edge, range, &peri)); 230 PetscCall(PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3])); 231 232 /* Get NODE info which associated with the current EDGE */ 233 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 234 if (mtype == DEGENERATE) { 235 PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id)); 236 } else { 237 params[0] = range[0]; 238 PetscCall(EG_evaluate(edge, params, result)); 239 PetscCall(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2])); 240 params[0] = range[1]; 241 PetscCall(EG_evaluate(edge, params, result)); 242 PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2])); 243 } 244 245 for (v = 0; v < Nv; ++v) { 246 ego vertex = nobjs[v]; 247 double limits[4]; 248 int dummy; 249 250 PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); 251 id = EG_indexBodyTopo(body, vertex); 252 PetscCall(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id)); 253 PetscCall(PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2])); 254 255 point[0] = point[0] + limits[0]; 256 point[1] = point[1] + limits[1]; 257 point[2] = point[2] + limits[2]; 258 } 259 } 260 } 261 EG_free(lobjs); 262 } 263 PetscFunctionReturn(0); 264 } 265 266 static PetscErrorCode DMPlexEGADSDestroy_Private(void *context) 267 { 268 if (context) EG_close((ego)context); 269 return 0; 270 } 271 272 static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm) 273 { 274 DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; 275 PetscInt cStart, cEnd, c; 276 /* EGADSLite variables */ 277 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 278 int oclass, mtype, nbodies, *senses; 279 int b; 280 /* PETSc variables */ 281 DM dm; 282 PetscHMapI edgeMap = NULL; 283 PetscInt dim = -1, cdim = -1, numCorners = 0, maxCorners = 0, numVertices = 0, newVertices = 0, numEdges = 0, numCells = 0, newCells = 0, numQuads = 0, cOff = 0, fOff = 0; 284 PetscInt *cells = NULL, *cone = NULL; 285 PetscReal *coords = NULL; 286 PetscMPIInt rank; 287 288 PetscFunctionBegin; 289 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 290 if (rank == 0) { 291 const PetscInt debug = 0; 292 293 /* --------------------------------------------------------------------------------------------------- 294 Generate Petsc Plex 295 Get all Nodes in model, record coordinates in a correctly formatted array 296 Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array 297 We need to uniformly refine the initial geometry to guarantee a valid mesh 298 */ 299 300 /* Calculate cell and vertex sizes */ 301 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 302 PetscCall(PetscHMapICreate(&edgeMap)); 303 numEdges = 0; 304 for (b = 0; b < nbodies; ++b) { 305 ego body = bodies[b]; 306 int id, Nl, l, Nv, v; 307 308 PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 309 for (l = 0; l < Nl; ++l) { 310 ego loop = lobjs[l]; 311 int Ner = 0, Ne, e, Nc; 312 313 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); 314 for (e = 0; e < Ne; ++e) { 315 ego edge = objs[e]; 316 int Nv, id; 317 PetscHashIter iter; 318 PetscBool found; 319 320 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 321 if (mtype == DEGENERATE) continue; 322 id = EG_indexBodyTopo(body, edge); 323 PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found)); 324 if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++)); 325 ++Ner; 326 } 327 if (Ner == 2) { 328 Nc = 2; 329 } else if (Ner == 3) { 330 Nc = 4; 331 } else if (Ner == 4) { 332 Nc = 8; 333 ++numQuads; 334 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner); 335 numCells += Nc; 336 newCells += Nc - 1; 337 maxCorners = PetscMax(Ner * 2 + 1, maxCorners); 338 } 339 PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); 340 for (v = 0; v < Nv; ++v) { 341 ego vertex = nobjs[v]; 342 343 id = EG_indexBodyTopo(body, vertex); 344 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 345 numVertices = PetscMax(id, numVertices); 346 } 347 EG_free(lobjs); 348 EG_free(nobjs); 349 } 350 PetscCall(PetscHMapIGetSize(edgeMap, &numEdges)); 351 newVertices = numEdges + numQuads; 352 numVertices += newVertices; 353 354 dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 355 cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 356 numCorners = 3; /* Split cells into triangles */ 357 PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone)); 358 359 /* Get vertex coordinates */ 360 for (b = 0; b < nbodies; ++b) { 361 ego body = bodies[b]; 362 int id, Nv, v; 363 364 PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); 365 for (v = 0; v < Nv; ++v) { 366 ego vertex = nobjs[v]; 367 double limits[4]; 368 int dummy; 369 370 PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); 371 id = EG_indexBodyTopo(body, vertex); 372 coords[(id - 1) * cdim + 0] = limits[0]; 373 coords[(id - 1) * cdim + 1] = limits[1]; 374 coords[(id - 1) * cdim + 2] = limits[2]; 375 } 376 EG_free(nobjs); 377 } 378 PetscCall(PetscHMapIClear(edgeMap)); 379 fOff = numVertices - newVertices + numEdges; 380 numEdges = 0; 381 numQuads = 0; 382 for (b = 0; b < nbodies; ++b) { 383 ego body = bodies[b]; 384 int Nl, l; 385 386 PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 387 for (l = 0; l < Nl; ++l) { 388 ego loop = lobjs[l]; 389 int lid, Ner = 0, Ne, e; 390 391 lid = EG_indexBodyTopo(body, loop); 392 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); 393 for (e = 0; e < Ne; ++e) { 394 ego edge = objs[e]; 395 int eid, Nv; 396 PetscHashIter iter; 397 PetscBool found; 398 399 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 400 if (mtype == DEGENERATE) continue; 401 ++Ner; 402 eid = EG_indexBodyTopo(body, edge); 403 PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found)); 404 if (!found) { 405 PetscInt v = numVertices - newVertices + numEdges; 406 double range[4], params[3] = {0., 0., 0.}, result[18]; 407 int periodic[2]; 408 409 PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++)); 410 PetscCall(EG_getRange(edge, range, periodic)); 411 params[0] = 0.5 * (range[0] + range[1]); 412 PetscCall(EG_evaluate(edge, params, result)); 413 coords[v * cdim + 0] = result[0]; 414 coords[v * cdim + 1] = result[1]; 415 coords[v * cdim + 2] = result[2]; 416 } 417 } 418 if (Ner == 4) { 419 PetscInt v = fOff + numQuads++; 420 ego *fobjs, face; 421 double range[4], params[3] = {0., 0., 0.}, result[18]; 422 int Nf, fid, periodic[2]; 423 424 PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); 425 face = fobjs[0]; 426 fid = EG_indexBodyTopo(body, face); 427 PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid); 428 PetscCall(EG_getRange(face, range, periodic)); 429 params[0] = 0.5 * (range[0] + range[1]); 430 params[1] = 0.5 * (range[2] + range[3]); 431 PetscCall(EG_evaluate(face, params, result)); 432 coords[v * cdim + 0] = result[0]; 433 coords[v * cdim + 1] = result[1]; 434 coords[v * cdim + 2] = result[2]; 435 } 436 } 437 } 438 PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices); 439 440 /* Get cell vertices by traversing loops */ 441 numQuads = 0; 442 cOff = 0; 443 for (b = 0; b < nbodies; ++b) { 444 ego body = bodies[b]; 445 int id, Nl, l; 446 447 PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 448 for (l = 0; l < Nl; ++l) { 449 ego loop = lobjs[l]; 450 int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t; 451 452 lid = EG_indexBodyTopo(body, loop); 453 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); 454 455 for (e = 0; e < Ne; ++e) { 456 ego edge = objs[e]; 457 int points[3]; 458 int eid, Nv, v, tmp; 459 460 eid = EG_indexBodyTopo(body, edge); 461 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 462 if (mtype == DEGENERATE) continue; 463 else ++Ner; 464 PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv); 465 466 for (v = 0; v < Nv; ++v) { 467 ego vertex = nobjs[v]; 468 469 id = EG_indexBodyTopo(body, vertex); 470 points[v * 2] = id - 1; 471 } 472 { 473 PetscInt edgeNum; 474 475 PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum)); 476 points[1] = numVertices - newVertices + edgeNum; 477 } 478 /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */ 479 if (!nc) { 480 for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v]; 481 } else { 482 if (cone[nc - 1] == points[0]) { 483 cone[nc++] = points[1]; 484 if (cone[0] != points[2]) cone[nc++] = points[2]; 485 } else if (cone[nc - 1] == points[2]) { 486 cone[nc++] = points[1]; 487 if (cone[0] != points[0]) cone[nc++] = points[0]; 488 } else if (cone[nc - 3] == points[0]) { 489 tmp = cone[nc - 3]; 490 cone[nc - 3] = cone[nc - 1]; 491 cone[nc - 1] = tmp; 492 cone[nc++] = points[1]; 493 if (cone[0] != points[2]) cone[nc++] = points[2]; 494 } else if (cone[nc - 3] == points[2]) { 495 tmp = cone[nc - 3]; 496 cone[nc - 3] = cone[nc - 1]; 497 cone[nc - 1] = tmp; 498 cone[nc++] = points[1]; 499 if (cone[0] != points[0]) cone[nc++] = points[0]; 500 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); 501 } 502 } 503 PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner); 504 if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++; 505 PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners); 506 /* Triangulate the loop */ 507 switch (Ner) { 508 case 2: /* Bi-Segment -> 2 triangles */ 509 Nt = 2; 510 cells[cOff * numCorners + 0] = cone[0]; 511 cells[cOff * numCorners + 1] = cone[1]; 512 cells[cOff * numCorners + 2] = cone[2]; 513 ++cOff; 514 cells[cOff * numCorners + 0] = cone[0]; 515 cells[cOff * numCorners + 1] = cone[2]; 516 cells[cOff * numCorners + 2] = cone[3]; 517 ++cOff; 518 break; 519 case 3: /* Triangle -> 4 triangles */ 520 Nt = 4; 521 cells[cOff * numCorners + 0] = cone[0]; 522 cells[cOff * numCorners + 1] = cone[1]; 523 cells[cOff * numCorners + 2] = cone[5]; 524 ++cOff; 525 cells[cOff * numCorners + 0] = cone[1]; 526 cells[cOff * numCorners + 1] = cone[2]; 527 cells[cOff * numCorners + 2] = cone[3]; 528 ++cOff; 529 cells[cOff * numCorners + 0] = cone[5]; 530 cells[cOff * numCorners + 1] = cone[3]; 531 cells[cOff * numCorners + 2] = cone[4]; 532 ++cOff; 533 cells[cOff * numCorners + 0] = cone[1]; 534 cells[cOff * numCorners + 1] = cone[3]; 535 cells[cOff * numCorners + 2] = cone[5]; 536 ++cOff; 537 break; 538 case 4: /* Quad -> 8 triangles */ 539 Nt = 8; 540 cells[cOff * numCorners + 0] = cone[0]; 541 cells[cOff * numCorners + 1] = cone[1]; 542 cells[cOff * numCorners + 2] = cone[7]; 543 ++cOff; 544 cells[cOff * numCorners + 0] = cone[1]; 545 cells[cOff * numCorners + 1] = cone[2]; 546 cells[cOff * numCorners + 2] = cone[3]; 547 ++cOff; 548 cells[cOff * numCorners + 0] = cone[3]; 549 cells[cOff * numCorners + 1] = cone[4]; 550 cells[cOff * numCorners + 2] = cone[5]; 551 ++cOff; 552 cells[cOff * numCorners + 0] = cone[5]; 553 cells[cOff * numCorners + 1] = cone[6]; 554 cells[cOff * numCorners + 2] = cone[7]; 555 ++cOff; 556 cells[cOff * numCorners + 0] = cone[8]; 557 cells[cOff * numCorners + 1] = cone[1]; 558 cells[cOff * numCorners + 2] = cone[3]; 559 ++cOff; 560 cells[cOff * numCorners + 0] = cone[8]; 561 cells[cOff * numCorners + 1] = cone[3]; 562 cells[cOff * numCorners + 2] = cone[5]; 563 ++cOff; 564 cells[cOff * numCorners + 0] = cone[8]; 565 cells[cOff * numCorners + 1] = cone[5]; 566 cells[cOff * numCorners + 2] = cone[7]; 567 ++cOff; 568 cells[cOff * numCorners + 0] = cone[8]; 569 cells[cOff * numCorners + 1] = cone[7]; 570 cells[cOff * numCorners + 2] = cone[1]; 571 ++cOff; 572 break; 573 default: 574 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner); 575 } 576 if (debug) { 577 for (t = 0; t < Nt; ++t) { 578 PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t)); 579 for (c = 0; c < numCorners; ++c) { 580 if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", ")); 581 PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c])); 582 } 583 PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n")); 584 } 585 } 586 } 587 EG_free(lobjs); 588 } 589 } 590 PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells); 591 PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 592 PetscCall(PetscFree3(coords, cells, cone)); 593 PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells)); 594 PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices)); 595 /* Embed EGADS model in DM */ 596 { 597 PetscContainer modelObj, contextObj; 598 599 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 600 PetscCall(PetscContainerSetPointer(modelObj, model)); 601 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); 602 PetscCall(PetscContainerDestroy(&modelObj)); 603 604 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 605 PetscCall(PetscContainerSetPointer(contextObj, context)); 606 PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 607 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); 608 PetscCall(PetscContainerDestroy(&contextObj)); 609 } 610 /* Label points */ 611 PetscCall(DMCreateLabel(dm, "EGADS Body ID")); 612 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 613 PetscCall(DMCreateLabel(dm, "EGADS Face ID")); 614 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 615 PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); 616 PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 617 PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); 618 PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 619 cOff = 0; 620 for (b = 0; b < nbodies; ++b) { 621 ego body = bodies[b]; 622 int id, Nl, l; 623 624 PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 625 for (l = 0; l < Nl; ++l) { 626 ego loop = lobjs[l]; 627 ego *fobjs; 628 int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t; 629 630 lid = EG_indexBodyTopo(body, loop); 631 PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); 632 PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf); 633 fid = EG_indexBodyTopo(body, fobjs[0]); 634 EG_free(fobjs); 635 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); 636 for (e = 0; e < Ne; ++e) { 637 ego edge = objs[e]; 638 int eid, Nv, v; 639 PetscInt points[3], support[2], numEdges, edgeNum; 640 const PetscInt *edges; 641 642 eid = EG_indexBodyTopo(body, edge); 643 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 644 if (mtype == DEGENERATE) continue; 645 else ++Ner; 646 for (v = 0; v < Nv; ++v) { 647 ego vertex = nobjs[v]; 648 649 id = EG_indexBodyTopo(body, vertex); 650 PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid)); 651 points[v * 2] = numCells + id - 1; 652 } 653 PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum)); 654 points[1] = numCells + numVertices - newVertices + edgeNum; 655 656 PetscCall(DMLabelSetValue(edgeLabel, points[1], eid)); 657 support[0] = points[0]; 658 support[1] = points[1]; 659 PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); 660 PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges); 661 PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid)); 662 PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges)); 663 support[0] = points[1]; 664 support[1] = points[2]; 665 PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); 666 PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges); 667 PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid)); 668 PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges)); 669 } 670 switch (Ner) { 671 case 2: 672 Nt = 2; 673 break; 674 case 3: 675 Nt = 4; 676 break; 677 case 4: 678 Nt = 8; 679 break; 680 default: 681 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner); 682 } 683 for (t = 0; t < Nt; ++t) { 684 PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b)); 685 PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid)); 686 } 687 cOff += Nt; 688 } 689 EG_free(lobjs); 690 } 691 PetscCall(PetscHMapIDestroy(&edgeMap)); 692 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 693 for (c = cStart; c < cEnd; ++c) { 694 PetscInt *closure = NULL; 695 PetscInt clSize, cl, bval, fval; 696 697 PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 698 PetscCall(DMLabelGetValue(bodyLabel, c, &bval)); 699 PetscCall(DMLabelGetValue(faceLabel, c, &fval)); 700 for (cl = 0; cl < clSize * 2; cl += 2) { 701 PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval)); 702 PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval)); 703 } 704 PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 705 } 706 *newdm = dm; 707 PetscFunctionReturn(0); 708 } 709 710 static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm) 711 { 712 DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; 713 // EGADS/EGADSLite variables 714 ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs; 715 ego topRef, prev, next; 716 int oclass, mtype, nbodies, *senses, *lSenses, *eSenses; 717 int b; 718 // PETSc variables 719 DM dm; 720 PetscHMapI edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL; 721 PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0; 722 PetscInt cellCntr = 0, numPoints = 0; 723 PetscInt *cells = NULL; 724 const PetscInt *cone = NULL; 725 PetscReal *coords = NULL; 726 PetscMPIInt rank; 727 728 PetscFunctionBeginUser; 729 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 730 if (rank == 0) { 731 // --------------------------------------------------------------------------------------------------- 732 // Generate Petsc Plex 733 // Get all Nodes in model, record coordinates in a correctly formatted array 734 // Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array 735 // We need to uniformly refine the initial geometry to guarantee a valid mesh 736 737 // Calculate cell and vertex sizes 738 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 739 740 PetscCall(PetscHMapICreate(&edgeMap)); 741 PetscCall(PetscHMapICreate(&bodyIndexMap)); 742 PetscCall(PetscHMapICreate(&bodyVertexMap)); 743 PetscCall(PetscHMapICreate(&bodyEdgeMap)); 744 PetscCall(PetscHMapICreate(&bodyEdgeGlobalMap)); 745 PetscCall(PetscHMapICreate(&bodyFaceMap)); 746 747 for (b = 0; b < nbodies; ++b) { 748 ego body = bodies[b]; 749 int Nf, Ne, Nv; 750 PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter; 751 PetscBool BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound; 752 753 PetscCall(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound)); 754 PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 755 PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 756 PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 757 PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 758 759 if (!BIfound) PetscCall(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices)); 760 if (!BVfound) PetscCall(PetscHMapISet(bodyVertexMap, b, numVertices)); 761 if (!BEfound) PetscCall(PetscHMapISet(bodyEdgeMap, b, numEdges)); 762 if (!BEGfound) PetscCall(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr)); 763 if (!BFfound) PetscCall(PetscHMapISet(bodyFaceMap, b, numFaces)); 764 765 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 766 PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); 767 PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); 768 EG_free(fobjs); 769 EG_free(eobjs); 770 EG_free(nobjs); 771 772 // Remove DEGENERATE EDGES from Edge count 773 PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); 774 int Netemp = 0; 775 for (int e = 0; e < Ne; ++e) { 776 ego edge = eobjs[e]; 777 int eid; 778 779 PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 780 eid = EG_indexBodyTopo(body, edge); 781 782 PetscCall(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound)); 783 if (mtype == DEGENERATE) { 784 if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1)); 785 } else { 786 ++Netemp; 787 if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp)); 788 } 789 } 790 EG_free(eobjs); 791 792 // Determine Number of Cells 793 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 794 for (int f = 0; f < Nf; ++f) { 795 ego face = fobjs[f]; 796 int edgeTemp = 0; 797 798 PetscCall(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs)); 799 for (int e = 0; e < Ne; ++e) { 800 ego edge = eobjs[e]; 801 802 PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 803 if (mtype != DEGENERATE) ++edgeTemp; 804 } 805 numCells += (2 * edgeTemp); 806 EG_free(eobjs); 807 } 808 EG_free(fobjs); 809 810 numFaces += Nf; 811 numEdges += Netemp; 812 numVertices += Nv; 813 edgeCntr += Ne; 814 } 815 816 // Set up basic DMPlex parameters 817 dim = 2; // Assumes 3D Models :: Need to handle 2D models in the future 818 cdim = 3; // Assumes 3D Models :: Need to update to handle 2D models in future 819 numCorners = 3; // Split Faces into triangles 820 numPoints = numVertices + numEdges + numFaces; // total number of coordinate points 821 822 PetscCall(PetscMalloc2(numPoints * cdim, &coords, numCells * numCorners, &cells)); 823 824 // Get Vertex Coordinates and Set up Cells 825 for (b = 0; b < nbodies; ++b) { 826 ego body = bodies[b]; 827 int Nf, Ne, Nv; 828 PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart; 829 PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter; 830 PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound; 831 832 // Vertices on Current Body 833 PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); 834 835 PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 836 PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); 837 PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart)); 838 839 for (int v = 0; v < Nv; ++v) { 840 ego vertex = nobjs[v]; 841 double limits[4]; 842 int id, dummy; 843 844 PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); 845 id = EG_indexBodyTopo(body, vertex); 846 847 coords[(bodyVertexIndexStart + id - 1) * cdim + 0] = limits[0]; 848 coords[(bodyVertexIndexStart + id - 1) * cdim + 1] = limits[1]; 849 coords[(bodyVertexIndexStart + id - 1) * cdim + 2] = limits[2]; 850 } 851 EG_free(nobjs); 852 853 // Edge Midpoint Vertices on Current Body 854 PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); 855 856 PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 857 PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); 858 PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); 859 860 PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 861 PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); 862 PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart)); 863 864 for (int e = 0; e < Ne; ++e) { 865 ego edge = eobjs[e]; 866 double range[2], avgt[1], cntrPnt[9]; 867 int eid, eOffset; 868 int periodic; 869 870 PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 871 if (mtype == DEGENERATE) continue; 872 873 eid = EG_indexBodyTopo(body, edge); 874 875 // get relative offset from globalEdgeID Vector 876 PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 877 PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1); 878 PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 879 880 PetscCall(EG_getRange(edge, range, &periodic)); 881 avgt[0] = (range[0] + range[1]) / 2.; 882 883 PetscCall(EG_evaluate(edge, avgt, cntrPnt)); 884 coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 0] = cntrPnt[0]; 885 coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 1] = cntrPnt[1]; 886 coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 2] = cntrPnt[2]; 887 } 888 EG_free(eobjs); 889 890 // Face Midpoint Vertices on Current Body 891 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 892 893 PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 894 PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); 895 PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart)); 896 897 for (int f = 0; f < Nf; ++f) { 898 ego face = fobjs[f]; 899 double range[4], avgUV[2], cntrPnt[18]; 900 int peri, id; 901 902 id = EG_indexBodyTopo(body, face); 903 PetscCall(EG_getRange(face, range, &peri)); 904 905 avgUV[0] = (range[0] + range[1]) / 2.; 906 avgUV[1] = (range[2] + range[3]) / 2.; 907 PetscCall(EG_evaluate(face, avgUV, cntrPnt)); 908 909 coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 0] = cntrPnt[0]; 910 coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 1] = cntrPnt[1]; 911 coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 2] = cntrPnt[2]; 912 } 913 EG_free(fobjs); 914 915 // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity 916 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 917 for (int f = 0; f < Nf; ++f) { 918 ego face = fobjs[f]; 919 int fID, midFaceID, midPntID, startID, endID, Nl; 920 921 fID = EG_indexBodyTopo(body, face); 922 midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1; 923 // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges. 924 // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are: 925 // 1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option. 926 // This will use a default EGADS tessellation as an initial surface mesh. 927 // 2) Create the initial surface mesh via a 2D mesher :: Currently not available (?future?) 928 // May I suggest the XXXX as a starting point? 929 930 PetscCall(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses)); 931 932 PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face has %d Loops. Can only handle Faces with 1 Loop. Please use --dm_plex_egads_with_tess = 1 Option", Nl); 933 for (int l = 0; l < Nl; ++l) { 934 ego loop = lobjs[l]; 935 936 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses)); 937 for (int e = 0; e < Ne; ++e) { 938 ego edge = eobjs[e]; 939 int eid, eOffset; 940 941 PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 942 eid = EG_indexBodyTopo(body, edge); 943 if (mtype == DEGENERATE) continue; 944 945 // get relative offset from globalEdgeID Vector 946 PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 947 PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1); 948 PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 949 950 midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1; 951 952 PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 953 954 if (eSenses[e] > 0) { 955 startID = EG_indexBodyTopo(body, nobjs[0]); 956 endID = EG_indexBodyTopo(body, nobjs[1]); 957 } else { 958 startID = EG_indexBodyTopo(body, nobjs[1]); 959 endID = EG_indexBodyTopo(body, nobjs[0]); 960 } 961 962 // Define 2 Cells per Edge with correct orientation 963 cells[cellCntr * numCorners + 0] = midFaceID; 964 cells[cellCntr * numCorners + 1] = bodyVertexIndexStart + startID - 1; 965 cells[cellCntr * numCorners + 2] = midPntID; 966 967 cells[cellCntr * numCorners + 3] = midFaceID; 968 cells[cellCntr * numCorners + 4] = midPntID; 969 cells[cellCntr * numCorners + 5] = bodyVertexIndexStart + endID - 1; 970 971 cellCntr = cellCntr + 2; 972 } 973 } 974 } 975 EG_free(fobjs); 976 } 977 } 978 979 // Generate DMPlex 980 PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 981 PetscCall(PetscFree2(coords, cells)); 982 PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " \n", numCells)); 983 PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " \n", numVertices)); 984 985 // Embed EGADS model in DM 986 { 987 PetscContainer modelObj, contextObj; 988 989 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 990 PetscCall(PetscContainerSetPointer(modelObj, model)); 991 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); 992 PetscCall(PetscContainerDestroy(&modelObj)); 993 994 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 995 PetscCall(PetscContainerSetPointer(contextObj, context)); 996 PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 997 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); 998 PetscCall(PetscContainerDestroy(&contextObj)); 999 } 1000 // Label points 1001 PetscInt nStart, nEnd; 1002 1003 PetscCall(DMCreateLabel(dm, "EGADS Body ID")); 1004 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 1005 PetscCall(DMCreateLabel(dm, "EGADS Face ID")); 1006 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 1007 PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); 1008 PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 1009 PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); 1010 PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 1011 1012 PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd)); 1013 1014 cellCntr = 0; 1015 for (b = 0; b < nbodies; ++b) { 1016 ego body = bodies[b]; 1017 int Nv, Ne, Nf; 1018 PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart; 1019 PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter; 1020 PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound; 1021 1022 PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 1023 PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); 1024 PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart)); 1025 1026 PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 1027 PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); 1028 PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); 1029 1030 PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 1031 PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); 1032 PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart)); 1033 1034 PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 1035 PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); 1036 PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart)); 1037 1038 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 1039 for (int f = 0; f < Nf; ++f) { 1040 ego face = fobjs[f]; 1041 int fID, Nl; 1042 1043 fID = EG_indexBodyTopo(body, face); 1044 1045 PetscCall(EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs)); 1046 for (int l = 0; l < Nl; ++l) { 1047 ego loop = lobjs[l]; 1048 int lid; 1049 1050 lid = EG_indexBodyTopo(body, loop); 1051 PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf); 1052 1053 PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses)); 1054 for (int e = 0; e < Ne; ++e) { 1055 ego edge = eobjs[e]; 1056 int eid, eOffset; 1057 1058 // Skip DEGENERATE Edges 1059 PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 1060 if (mtype == DEGENERATE) continue; 1061 eid = EG_indexBodyTopo(body, edge); 1062 1063 // get relative offset from globalEdgeID Vector 1064 PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 1065 PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1); 1066 PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 1067 1068 PetscCall(EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs)); 1069 for (int v = 0; v < Nv; ++v) { 1070 ego vertex = nobjs[v]; 1071 int vID; 1072 1073 vID = EG_indexBodyTopo(body, vertex); 1074 PetscCall(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b)); 1075 PetscCall(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID)); 1076 } 1077 EG_free(nobjs); 1078 1079 PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b)); 1080 PetscCall(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid)); 1081 1082 // Define Cell faces 1083 for (int jj = 0; jj < 2; ++jj) { 1084 PetscCall(DMLabelSetValue(bodyLabel, cellCntr, b)); 1085 PetscCall(DMLabelSetValue(faceLabel, cellCntr, fID)); 1086 PetscCall(DMPlexGetCone(dm, cellCntr, &cone)); 1087 1088 PetscCall(DMLabelSetValue(bodyLabel, cone[0], b)); 1089 PetscCall(DMLabelSetValue(faceLabel, cone[0], fID)); 1090 1091 PetscCall(DMLabelSetValue(bodyLabel, cone[1], b)); 1092 PetscCall(DMLabelSetValue(edgeLabel, cone[1], eid)); 1093 1094 PetscCall(DMLabelSetValue(bodyLabel, cone[2], b)); 1095 PetscCall(DMLabelSetValue(faceLabel, cone[2], fID)); 1096 1097 cellCntr = cellCntr + 1; 1098 } 1099 } 1100 } 1101 EG_free(lobjs); 1102 1103 PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b)); 1104 PetscCall(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID)); 1105 } 1106 EG_free(fobjs); 1107 } 1108 1109 PetscCall(PetscHMapIDestroy(&edgeMap)); 1110 PetscCall(PetscHMapIDestroy(&bodyIndexMap)); 1111 PetscCall(PetscHMapIDestroy(&bodyVertexMap)); 1112 PetscCall(PetscHMapIDestroy(&bodyEdgeMap)); 1113 PetscCall(PetscHMapIDestroy(&bodyEdgeGlobalMap)); 1114 PetscCall(PetscHMapIDestroy(&bodyFaceMap)); 1115 1116 *newdm = dm; 1117 PetscFunctionReturn(0); 1118 } 1119 1120 static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm) 1121 { 1122 DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; 1123 /* EGADSLite variables */ 1124 ego geom, *bodies, *fobjs; 1125 int b, oclass, mtype, nbodies, *senses; 1126 int totalNumTris = 0, totalNumPoints = 0; 1127 double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize; 1128 /* PETSc variables */ 1129 DM dm; 1130 PetscHMapI pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL; 1131 PetscHMapI pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL; 1132 PetscInt dim = -1, cdim = -1, numCorners = 0, counter = 0; 1133 PetscInt *cells = NULL; 1134 const PetscInt *cone = NULL; 1135 PetscReal *coords = NULL; 1136 PetscMPIInt rank; 1137 1138 PetscFunctionBeginUser; 1139 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1140 if (rank == 0) { 1141 // --------------------------------------------------------------------------------------------------- 1142 // Generate Petsc Plex from EGADSlite created Tessellation of geometry 1143 // --------------------------------------------------------------------------------------------------- 1144 1145 // Calculate cell and vertex sizes 1146 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 1147 1148 PetscCall(PetscHMapICreate(&pointIndexStartMap)); 1149 PetscCall(PetscHMapICreate(&triIndexStartMap)); 1150 PetscCall(PetscHMapICreate(&pTypeLabelMap)); 1151 PetscCall(PetscHMapICreate(&pIndexLabelMap)); 1152 PetscCall(PetscHMapICreate(&pBodyIndexLabelMap)); 1153 PetscCall(PetscHMapICreate(&triFaceIDLabelMap)); 1154 PetscCall(PetscHMapICreate(&triBodyIDLabelMap)); 1155 1156 /* Create Tessellation of Bodies */ 1157 ego tessArray[nbodies]; 1158 1159 for (b = 0; b < nbodies; ++b) { 1160 ego body = bodies[b]; 1161 double params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation 1162 int Nf, bodyNumPoints = 0, bodyNumTris = 0; 1163 PetscHashIter PISiter, TISiter; 1164 PetscBool PISfound, TISfound; 1165 1166 /* Store Start Index for each Body's Point and Tris */ 1167 PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); 1168 PetscCall(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound)); 1169 1170 if (!PISfound) PetscCall(PetscHMapISet(pointIndexStartMap, b, totalNumPoints)); 1171 if (!TISfound) PetscCall(PetscHMapISet(triIndexStartMap, b, totalNumTris)); 1172 1173 /* Calculate Tessellation parameters based on Bounding Box */ 1174 /* Get Bounding Box Dimensions of the BODY */ 1175 PetscCall(EG_getBoundingBox(body, boundBox)); 1176 tessSize = boundBox[3] - boundBox[0]; 1177 if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1]; 1178 if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2]; 1179 1180 // TODO :: May want to give users tessellation parameter options // 1181 params[0] = 0.0250 * tessSize; 1182 params[1] = 0.0075 * tessSize; 1183 params[2] = 15.0; 1184 1185 PetscCall(EG_makeTessBody(body, params, &tessArray[b])); 1186 1187 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 1188 1189 for (int f = 0; f < Nf; ++f) { 1190 ego face = fobjs[f]; 1191 int len, fID, ntris; 1192 const int *ptype, *pindex, *ptris, *ptric; 1193 const double *pxyz, *puv; 1194 1195 // Get Face ID // 1196 fID = EG_indexBodyTopo(body, face); 1197 1198 // Checkout the Surface Tessellation // 1199 PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric)); 1200 1201 // Determine total number of triangle cells in the tessellation // 1202 bodyNumTris += (int)ntris; 1203 1204 // Check out the point index and coordinate // 1205 for (int p = 0; p < len; ++p) { 1206 int global; 1207 1208 PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global)); 1209 1210 // Determine the total number of points in the tessellation // 1211 bodyNumPoints = PetscMax(bodyNumPoints, global); 1212 } 1213 } 1214 EG_free(fobjs); 1215 1216 totalNumPoints += bodyNumPoints; 1217 totalNumTris += bodyNumTris; 1218 } 1219 //} - Original End of (rank == 0) 1220 1221 dim = 2; 1222 cdim = 3; 1223 numCorners = 3; 1224 //PetscInt counter = 0; 1225 1226 /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA */ 1227 /* Fill in below and use to define DMLabels after DMPlex creation */ 1228 PetscCall(PetscMalloc2(totalNumPoints * cdim, &coords, totalNumTris * numCorners, &cells)); 1229 1230 for (b = 0; b < nbodies; ++b) { 1231 ego body = bodies[b]; 1232 int Nf; 1233 PetscInt pointIndexStart; 1234 PetscHashIter PISiter; 1235 PetscBool PISfound; 1236 1237 PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); 1238 PetscCheck(PISfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b); 1239 PetscCall(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart)); 1240 1241 PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 1242 1243 for (int f = 0; f < Nf; ++f) { 1244 /* Get Face Object */ 1245 ego face = fobjs[f]; 1246 int len, fID, ntris; 1247 const int *ptype, *pindex, *ptris, *ptric; 1248 const double *pxyz, *puv; 1249 1250 /* Get Face ID */ 1251 fID = EG_indexBodyTopo(body, face); 1252 1253 /* Checkout the Surface Tessellation */ 1254 PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric)); 1255 1256 /* Check out the point index and coordinate */ 1257 for (int p = 0; p < len; ++p) { 1258 int global; 1259 PetscHashIter PTLiter, PILiter, PBLiter; 1260 PetscBool PTLfound, PILfound, PBLfound; 1261 1262 PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global)); 1263 1264 /* Set the coordinates array for DAG */ 1265 coords[((global - 1 + pointIndexStart) * 3) + 0] = pxyz[(p * 3) + 0]; 1266 coords[((global - 1 + pointIndexStart) * 3) + 1] = pxyz[(p * 3) + 1]; 1267 coords[((global - 1 + pointIndexStart) * 3) + 2] = pxyz[(p * 3) + 2]; 1268 1269 /* Store Geometry Label Information for DMLabel assignment later */ 1270 PetscCall(PetscHMapIFind(pTypeLabelMap, global - 1 + pointIndexStart, &PTLiter, &PTLfound)); 1271 PetscCall(PetscHMapIFind(pIndexLabelMap, global - 1 + pointIndexStart, &PILiter, &PILfound)); 1272 PetscCall(PetscHMapIFind(pBodyIndexLabelMap, global - 1 + pointIndexStart, &PBLiter, &PBLfound)); 1273 1274 if (!PTLfound) PetscCall(PetscHMapISet(pTypeLabelMap, global - 1 + pointIndexStart, ptype[p])); 1275 if (!PILfound) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, pindex[p])); 1276 if (!PBLfound) PetscCall(PetscHMapISet(pBodyIndexLabelMap, global - 1 + pointIndexStart, b)); 1277 1278 if (ptype[p] < 0) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, fID)); 1279 } 1280 1281 for (int t = 0; t < (int)ntris; ++t) { 1282 int global, globalA, globalB; 1283 PetscHashIter TFLiter, TBLiter; 1284 PetscBool TFLfound, TBLfound; 1285 1286 PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 0], &global)); 1287 cells[(counter * 3) + 0] = global - 1 + pointIndexStart; 1288 1289 PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 1], &globalA)); 1290 cells[(counter * 3) + 1] = globalA - 1 + pointIndexStart; 1291 1292 PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 2], &globalB)); 1293 cells[(counter * 3) + 2] = globalB - 1 + pointIndexStart; 1294 1295 PetscCall(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound)); 1296 PetscCall(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound)); 1297 1298 if (!TFLfound) PetscCall(PetscHMapISet(triFaceIDLabelMap, counter, fID)); 1299 if (!TBLfound) PetscCall(PetscHMapISet(triBodyIDLabelMap, counter, b)); 1300 1301 counter += 1; 1302 } 1303 } 1304 EG_free(fobjs); 1305 } 1306 } 1307 1308 //Build DMPlex 1309 PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 1310 PetscCall(PetscFree2(coords, cells)); 1311 1312 // Embed EGADS model in DM 1313 { 1314 PetscContainer modelObj, contextObj; 1315 1316 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 1317 PetscCall(PetscContainerSetPointer(modelObj, model)); 1318 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); 1319 PetscCall(PetscContainerDestroy(&modelObj)); 1320 1321 PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 1322 PetscCall(PetscContainerSetPointer(contextObj, context)); 1323 PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 1324 PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); 1325 PetscCall(PetscContainerDestroy(&contextObj)); 1326 } 1327 1328 // Label Points 1329 PetscCall(DMCreateLabel(dm, "EGADS Body ID")); 1330 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 1331 PetscCall(DMCreateLabel(dm, "EGADS Face ID")); 1332 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 1333 PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); 1334 PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 1335 PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); 1336 PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 1337 1338 /* Get Number of DAG Nodes at each level */ 1339 int fStart, fEnd, eStart, eEnd, nStart, nEnd; 1340 1341 PetscCall(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd)); 1342 PetscCall(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd)); 1343 PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd)); 1344 1345 /* Set DMLabels for NODES */ 1346 for (int n = nStart; n < nEnd; ++n) { 1347 int pTypeVal, pIndexVal, pBodyVal; 1348 PetscHashIter PTLiter, PILiter, PBLiter; 1349 PetscBool PTLfound, PILfound, PBLfound; 1350 1351 //Converted to Hash Tables 1352 PetscCall(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound)); 1353 PetscCheck(PTLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n); 1354 PetscCall(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal)); 1355 1356 PetscCall(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound)); 1357 PetscCheck(PILfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n); 1358 PetscCall(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal)); 1359 1360 PetscCall(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound)); 1361 PetscCheck(PBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n); 1362 PetscCall(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal)); 1363 1364 PetscCall(DMLabelSetValue(bodyLabel, n, pBodyVal)); 1365 if (pTypeVal == 0) PetscCall(DMLabelSetValue(vertexLabel, n, pIndexVal)); 1366 if (pTypeVal > 0) PetscCall(DMLabelSetValue(edgeLabel, n, pIndexVal)); 1367 if (pTypeVal < 0) PetscCall(DMLabelSetValue(faceLabel, n, pIndexVal)); 1368 } 1369 1370 /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */ 1371 for (int e = eStart; e < eEnd; ++e) { 1372 int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1; 1373 1374 PetscCall(DMPlexGetCone(dm, e, &cone)); 1375 PetscCall(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end? 1376 PetscCall(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0)); 1377 PetscCall(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1)); 1378 PetscCall(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0)); 1379 PetscCall(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1)); 1380 PetscCall(DMLabelGetValue(faceLabel, cone[0], &faceID_0)); 1381 PetscCall(DMLabelGetValue(faceLabel, cone[1], &faceID_1)); 1382 1383 PetscCall(DMLabelSetValue(bodyLabel, e, bodyID_0)); 1384 1385 if (edgeID_0 == edgeID_1) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0)); 1386 else if (vertexID_0 > 0 && edgeID_1 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_1)); 1387 else if (vertexID_1 > 0 && edgeID_0 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0)); 1388 else { /* Do Nothing */ } 1389 } 1390 1391 /* Set DMLabels for Cells */ 1392 for (int f = fStart; f < fEnd; ++f) { 1393 int edgeID_0; 1394 PetscInt triBodyVal, triFaceVal; 1395 PetscHashIter TFLiter, TBLiter; 1396 PetscBool TFLfound, TBLfound; 1397 1398 // Convert to Hash Table 1399 PetscCall(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound)); 1400 PetscCheck(TFLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f); 1401 PetscCall(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal)); 1402 1403 PetscCall(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound)); 1404 PetscCheck(TBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f); 1405 PetscCall(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal)); 1406 1407 PetscCall(DMLabelSetValue(bodyLabel, f, triBodyVal)); 1408 PetscCall(DMLabelSetValue(faceLabel, f, triFaceVal)); 1409 1410 /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */ 1411 PetscCall(DMPlexGetCone(dm, f, &cone)); 1412 1413 for (int jj = 0; jj < 3; ++jj) { 1414 PetscCall(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0)); 1415 1416 if (edgeID_0 < 0) { 1417 PetscCall(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal)); 1418 PetscCall(DMLabelSetValue(faceLabel, cone[jj], triFaceVal)); 1419 } 1420 } 1421 } 1422 1423 *newdm = dm; 1424 PetscFunctionReturn(0); 1425 } 1426 #endif 1427 1428 /*@ 1429 DMPlexInflateToGeomModel - Snaps the vertex coordinates of a `DMPLEX` object representing the mesh to its geometry if some vertices depart from the model. This usually happens with non-conforming refinement. 1430 1431 Collective on dm 1432 1433 Input Parameter: 1434 . dm - The uninflated `DM` object representing the mesh 1435 1436 Output Parameter: 1437 . dm - The inflated `DM` object representing the mesh 1438 1439 Level: intermediate 1440 1441 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()` 1442 @*/ 1443 PetscErrorCode DMPlexInflateToGeomModel(DM dm) 1444 { 1445 #if defined(PETSC_HAVE_EGADS) 1446 /* EGADS Variables */ 1447 ego model, geom, body, face, edge; 1448 ego *bodies; 1449 int Nb, oclass, mtype, *senses; 1450 double result[3]; 1451 /* PETSc Variables */ 1452 DM cdm; 1453 PetscContainer modelObj; 1454 DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; 1455 Vec coordinates; 1456 PetscScalar *coords; 1457 PetscInt bodyID, faceID, edgeID, vertexID; 1458 PetscInt cdim, d, vStart, vEnd, v; 1459 #endif 1460 1461 PetscFunctionBegin; 1462 #if defined(PETSC_HAVE_EGADS) 1463 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj)); 1464 if (!modelObj) PetscFunctionReturn(0); 1465 PetscCall(DMGetCoordinateDim(dm, &cdim)); 1466 PetscCall(DMGetCoordinateDM(dm, &cdm)); 1467 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 1468 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 1469 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 1470 PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 1471 PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 1472 1473 PetscCall(PetscContainerGetPointer(modelObj, (void **)&model)); 1474 PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 1475 1476 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1477 PetscCall(VecGetArrayWrite(coordinates, &coords)); 1478 for (v = vStart; v < vEnd; ++v) { 1479 PetscScalar *vcoords; 1480 1481 PetscCall(DMLabelGetValue(bodyLabel, v, &bodyID)); 1482 PetscCall(DMLabelGetValue(faceLabel, v, &faceID)); 1483 PetscCall(DMLabelGetValue(edgeLabel, v, &edgeID)); 1484 PetscCall(DMLabelGetValue(vertexLabel, v, &vertexID)); 1485 1486 PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb); 1487 body = bodies[bodyID]; 1488 1489 PetscCall(DMPlexPointLocalRef(cdm, v, coords, (void *)&vcoords)); 1490 if (edgeID > 0) { 1491 /* Snap to EDGE at nearest location */ 1492 double params[1]; 1493 PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &edge)); 1494 PetscCall(EG_invEvaluate(edge, vcoords, params, result)); // Get (x,y,z) of nearest point on EDGE 1495 for (d = 0; d < cdim; ++d) vcoords[d] = result[d]; 1496 } else if (faceID > 0) { 1497 /* Snap to FACE at nearest location */ 1498 double params[2]; 1499 PetscCall(EG_objectBodyTopo(body, FACE, faceID, &face)); 1500 PetscCall(EG_invEvaluate(face, vcoords, params, result)); // Get (x,y,z) of nearest point on FACE 1501 for (d = 0; d < cdim; ++d) vcoords[d] = result[d]; 1502 } 1503 } 1504 PetscCall(VecRestoreArrayWrite(coordinates, &coords)); 1505 /* Clear out global coordinates */ 1506 PetscCall(VecDestroy(&dm->coordinates[0].x)); 1507 #endif 1508 PetscFunctionReturn(0); 1509 } 1510 1511 /*@C 1512 DMPlexCreateEGADSFromFile - Create a `DMPLEX` mesh from an EGADS, IGES, or STEP file. 1513 1514 Collective 1515 1516 Input Parameters: 1517 + comm - The MPI communicator 1518 - filename - The name of the EGADS, IGES, or STEP file 1519 1520 Output Parameter: 1521 . dm - The `DM` object representing the mesh 1522 1523 Level: beginner 1524 1525 .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSLiteFromFile()` 1526 @*/ 1527 PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm) 1528 { 1529 PetscMPIInt rank; 1530 #if defined(PETSC_HAVE_EGADS) 1531 ego context = NULL, model = NULL; 1532 #endif 1533 PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE; 1534 1535 PetscFunctionBegin; 1536 PetscValidCharPointer(filename, 2); 1537 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL)); 1538 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL)); 1539 PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL)); 1540 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 1541 #if defined(PETSC_HAVE_EGADS) 1542 if (rank == 0) { 1543 PetscCall(EG_open(&context)); 1544 PetscCall(EG_loadModel(context, 0, filename, &model)); 1545 if (printModel) PetscCall(DMPlexEGADSPrintModel_Internal(model)); 1546 } 1547 if (tessModel) PetscCall(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm)); 1548 else if (newModel) PetscCall(DMPlexCreateEGADS(comm, context, model, dm)); 1549 else PetscCall(DMPlexCreateEGADS_Internal(comm, context, model, dm)); 1550 PetscFunctionReturn(0); 1551 #else 1552 SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads"); 1553 #endif 1554 } 1555