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