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