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