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