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