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 CHKERRQ(DMGetCoordinateDM(dm, &cdm)); 40 CHKERRQ(DMGetCoordinateDim(dm, &dE)); 41 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinatesLocal)); 42 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 43 PetscCheckFalse(bodyID >= Nb,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb); 44 body = bodies[bodyID]; 45 46 if (edgeID >= 0) {CHKERRQ(EG_objectBodyTopo(body, EDGE, edgeID, &obj)); Np = 1;} 47 else if (faceID >= 0) {CHKERRQ(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 CHKERRQ(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); 55 Nv /= dE; 56 if (Nv == 1) { 57 CHKERRQ(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); 58 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 59 PetscFunctionReturn(0); 60 } 61 PetscCheckFalse(Nv > 16,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p); 62 63 /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */ 64 CHKERRQ(EG_getRange(obj, range, &peri)); 65 for (v = 0; v < Nv; ++v) { 66 CHKERRQ(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 CHKERRQ(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 PetscCheckFalse((params[0] < range[0]) || (params[0] > range[1]),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D had bad interpolation", p); 86 PetscCheckFalse(Np > 1 && ((params[1] < range[2]) || (params[1] > range[3])),PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D had bad interpolation", p); 87 /* Put coordinates for new vertex in result[] */ 88 CHKERRQ(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 CHKERRQ(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 129 CHKERRQ(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 130 CHKERRQ(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 CHKERRQ(PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj)); 136 if (!modelObj) { 137 CHKERRQ(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 CHKERRQ(PetscContainerGetPointer(modelObj, (void **) &model)); 145 CHKERRQ(DMLabelGetValue(bodyLabel, p, &bodyID)); 146 CHKERRQ(DMLabelGetValue(faceLabel, p, &faceID)); 147 CHKERRQ(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) CHKERRQ(DMPlexSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords)); 154 else CHKERRQ(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 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 172 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs)); 180 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh)); 181 EG_free(objs); 182 183 CHKERRQ(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs)); 184 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf)); 185 EG_free(objs); 186 187 CHKERRQ(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); 188 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl)); 189 190 CHKERRQ(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs)); 191 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne)); 192 EG_free(objs); 193 194 CHKERRQ(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs)); 195 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id)); 203 204 /* Get EDGE info which associated with the current LOOP */ 205 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e)); 217 218 CHKERRQ(EG_getRange(edge, range, &peri)); 219 CHKERRQ(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 CHKERRQ(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 223 if (mtype == DEGENERATE) { 224 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id)); 225 } else { 226 params[0] = range[0]; 227 CHKERRQ(EG_evaluate(edge, params, result)); 228 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2])); 229 params[0] = range[1]; 230 CHKERRQ(EG_evaluate(edge, params, result)); 231 CHKERRQ(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 CHKERRQ(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); 240 id = EG_indexBodyTopo(body, vertex); 241 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id)); 242 CHKERRQ(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 CHKERRMPI(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 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 291 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 310 if (mtype == DEGENERATE) continue; 311 id = EG_indexBodyTopo(body, edge); 312 CHKERRQ(PetscHMapIFind(edgeMap, id-1, &iter, &found)); 313 if (!found) CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 385 if (mtype == DEGENERATE) continue; 386 ++Ner; 387 eid = EG_indexBodyTopo(body, edge); 388 CHKERRQ(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 CHKERRQ(PetscHMapISet(edgeMap, eid-1, numEdges++)); 395 CHKERRQ(EG_getRange(edge, range, periodic)); 396 params[0] = 0.5*(range[0] + range[1]); 397 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); 410 face = fobjs[0]; 411 fid = EG_indexBodyTopo(body, face); 412 PetscCheckFalse(Nf != 1,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid); 413 CHKERRQ(EG_getRange(face, range, periodic)); 414 params[0] = 0.5*(range[0] + range[1]); 415 params[1] = 0.5*(range[2] + range[3]); 416 CHKERRQ(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 PetscCheckFalse(numEdges + numQuads != newVertices,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D 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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); 447 if (mtype == DEGENERATE) continue; 448 else ++Ner; 449 PetscCheckFalse(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 CHKERRQ(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 PetscCheckFalse(nc != 2*Ner,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner); 475 if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;} 476 PetscCheckFalse(nc > maxCorners,PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D 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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t)); 549 for (c = 0; c < numCorners; ++c) { 550 if (c > 0) CHKERRQ(PetscPrintf(PETSC_COMM_SELF, ", ")); 551 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c])); 552 } 553 CHKERRQ(PetscPrintf(PETSC_COMM_SELF, ")\n")); 554 } 555 } 556 } 557 EG_free(lobjs); 558 } 559 } 560 PetscCheckFalse(cOff != numCells,PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells); 561 CHKERRQ(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 562 CHKERRQ(PetscFree3(coords, cells, cone)); 563 CHKERRQ(PetscInfo(dm, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells)); 564 CHKERRQ(PetscInfo(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices)); 565 /* Embed EGADS model in DM */ 566 { 567 PetscContainer modelObj, contextObj; 568 569 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 570 CHKERRQ(PetscContainerSetPointer(modelObj, model)); 571 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj)); 572 CHKERRQ(PetscContainerDestroy(&modelObj)); 573 574 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 575 CHKERRQ(PetscContainerSetPointer(contextObj, context)); 576 CHKERRQ(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 577 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj)); 578 CHKERRQ(PetscContainerDestroy(&contextObj)); 579 } 580 /* Label points */ 581 CHKERRQ(DMCreateLabel(dm, "EGADS Body ID")); 582 CHKERRQ(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 583 CHKERRQ(DMCreateLabel(dm, "EGADS Face ID")); 584 CHKERRQ(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 585 CHKERRQ(DMCreateLabel(dm, "EGADS Edge ID")); 586 CHKERRQ(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 587 CHKERRQ(DMCreateLabel(dm, "EGADS Vertex ID")); 588 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); 602 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(DMLabelSetValue(edgeLabel, numCells + id-1, eid)); 621 points[v*2] = numCells + id-1; 622 } 623 CHKERRQ(PetscHMapIGet(edgeMap, eid-1, &edgeNum)); 624 points[1] = numCells + numVertices - newVertices + edgeNum; 625 626 CHKERRQ(DMLabelSetValue(edgeLabel, points[1], eid)); 627 support[0] = points[0]; 628 support[1] = points[1]; 629 CHKERRQ(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); 630 PetscCheckFalse(numEdges != 1,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges); 631 CHKERRQ(DMLabelSetValue(edgeLabel, edges[0], eid)); 632 CHKERRQ(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges)); 633 support[0] = points[1]; 634 support[1] = points[2]; 635 CHKERRQ(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); 636 PetscCheckFalse(numEdges != 1,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges); 637 CHKERRQ(DMLabelSetValue(edgeLabel, edges[0], eid)); 638 CHKERRQ(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 CHKERRQ(DMLabelSetValue(bodyLabel, cOff+t, b)); 648 CHKERRQ(DMLabelSetValue(faceLabel, cOff+t, fid)); 649 } 650 cOff += Nt; 651 } 652 EG_free(lobjs); 653 } 654 CHKERRQ(PetscHMapIDestroy(&edgeMap)); 655 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 656 for (c = cStart; c < cEnd; ++c) { 657 PetscInt *closure = NULL; 658 PetscInt clSize, cl, bval, fval; 659 660 CHKERRQ(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); 661 CHKERRQ(DMLabelGetValue(bodyLabel, c, &bval)); 662 CHKERRQ(DMLabelGetValue(faceLabel, c, &fval)); 663 for (cl = 0; cl < clSize*2; cl += 2) { 664 CHKERRQ(DMLabelSetValue(bodyLabel, closure[cl], bval)); 665 CHKERRQ(DMLabelSetValue(faceLabel, closure[cl], fval)); 666 } 667 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 693 if (!rank) { 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 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 702 703 CHKERRQ(PetscHMapICreate(&edgeMap)); 704 CHKERRQ(PetscHMapICreate(&bodyIndexMap)); 705 CHKERRQ(PetscHMapICreate(&bodyVertexMap)); 706 CHKERRQ(PetscHMapICreate(&bodyEdgeMap)); 707 CHKERRQ(PetscHMapICreate(&bodyEdgeGlobalMap)); 708 CHKERRQ(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 CHKERRQ(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound)); 717 CHKERRQ(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 718 CHKERRQ(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 719 CHKERRQ(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 720 CHKERRQ(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 721 722 if (!BIfound) CHKERRQ(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices)); 723 if (!BVfound) CHKERRQ(PetscHMapISet(bodyVertexMap, b, numVertices)); 724 if (!BEfound) CHKERRQ(PetscHMapISet(bodyEdgeMap, b, numEdges)); 725 if (!BEGfound) CHKERRQ(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr)); 726 if (!BFfound) CHKERRQ(PetscHMapISet(bodyFaceMap, b, numFaces)); 727 728 CHKERRQ(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 729 CHKERRQ(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); 730 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); 743 eid = EG_indexBodyTopo(body, edge); 744 745 CHKERRQ(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound)); 746 if (mtype == DEGENERATE) { 747 if (!EMfound) CHKERRQ(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1)); 748 } 749 else { 750 ++Netemp; 751 if (!EMfound) CHKERRQ(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp)); 752 } 753 } 754 EG_free(eobjs); 755 756 // Determine Number of Cells 757 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs)); 763 for (int e = 0; e < Ne; ++e) { 764 ego edge = eobjs[e]; 765 766 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); 798 799 CHKERRQ(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 800 PetscCheckFalse(!BVfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); 801 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); 819 820 CHKERRQ(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 821 PetscCheckFalse(!BEfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); 822 CHKERRQ(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); 823 824 CHKERRQ(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 825 PetscCheckFalse(!BEGfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); 826 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 841 PetscCheckFalse(!EMfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1); 842 CHKERRQ(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 843 844 CHKERRQ(EG_getRange(edge, range, &periodic)); 845 avgt[0] = (range[0] + range[1]) / 2.; 846 847 CHKERRQ(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 CHKERRQ(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); 856 857 CHKERRQ(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 858 PetscCheckFalse(!BFfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); 859 CHKERRQ(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 CHKERRQ(EG_getRange(face, range, &peri)); 868 869 avgUV[0] = (range[0] + range[1]) / 2.; 870 avgUV[1] = (range[2] + range[3]) / 2.; 871 CHKERRQ(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 CHKERRQ(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 CHKERRQ(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses)); 895 896 PetscCheckFalse(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 911 PetscCheckFalse(!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 CHKERRQ(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 913 914 midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1; 915 916 CHKERRQ(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 CHKERRQ(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 940 CHKERRQ(PetscFree2(coords, cells)); 941 CHKERRQ(PetscInfo(dm, " Total Number of Unique Cells = %D \n", numCells)); 942 CHKERRQ(PetscInfo(dm, " Total Number of Unique Vertices = %D \n", numVertices)); 943 944 // Embed EGADS model in DM 945 { 946 PetscContainer modelObj, contextObj; 947 948 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 949 CHKERRQ(PetscContainerSetPointer(modelObj, model)); 950 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj)); 951 CHKERRQ(PetscContainerDestroy(&modelObj)); 952 953 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 954 CHKERRQ(PetscContainerSetPointer(contextObj, context)); 955 CHKERRQ(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 956 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj)); 957 CHKERRQ(PetscContainerDestroy(&contextObj)); 958 } 959 // Label points 960 PetscInt nStart, nEnd; 961 962 CHKERRQ(DMCreateLabel(dm, "EGADS Body ID")); 963 CHKERRQ(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 964 CHKERRQ(DMCreateLabel(dm, "EGADS Face ID")); 965 CHKERRQ(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 966 CHKERRQ(DMCreateLabel(dm, "EGADS Edge ID")); 967 CHKERRQ(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 968 CHKERRQ(DMCreateLabel(dm, "EGADS Vertex ID")); 969 CHKERRQ(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 970 971 CHKERRQ(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 CHKERRQ(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); 982 PetscCheckFalse(!BVfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); 983 CHKERRQ(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart)); 984 985 CHKERRQ(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); 986 PetscCheckFalse(!BEfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); 987 CHKERRQ(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); 988 989 CHKERRQ(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); 990 PetscCheckFalse(!BFfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); 991 CHKERRQ(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart)); 992 993 CHKERRQ(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); 994 PetscCheckFalse(!BEGfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); 995 CHKERRQ(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart)); 996 997 CHKERRQ(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 CHKERRQ(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 PetscCheckFalse(Nl > 1,PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf); 1011 1012 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); 1024 PetscCheckFalse(!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 CHKERRQ(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); 1026 1027 CHKERRQ(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 CHKERRQ(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b)); 1034 CHKERRQ(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID)); 1035 } 1036 EG_free(nobjs); 1037 1038 CHKERRQ(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b)); 1039 CHKERRQ(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid)); 1040 1041 // Define Cell faces 1042 for (int jj = 0; jj < 2; ++jj){ 1043 CHKERRQ(DMLabelSetValue(bodyLabel, cellCntr, b)); 1044 CHKERRQ(DMLabelSetValue(faceLabel, cellCntr, fID)); 1045 CHKERRQ(DMPlexGetCone(dm, cellCntr, &cone)); 1046 1047 CHKERRQ(DMLabelSetValue(bodyLabel, cone[0], b)); 1048 CHKERRQ(DMLabelSetValue(faceLabel, cone[0], fID)); 1049 1050 CHKERRQ(DMLabelSetValue(bodyLabel, cone[1], b)); 1051 CHKERRQ(DMLabelSetValue(edgeLabel, cone[1], eid)); 1052 1053 CHKERRQ(DMLabelSetValue(bodyLabel, cone[2], b)); 1054 CHKERRQ(DMLabelSetValue(faceLabel, cone[2], fID)); 1055 1056 cellCntr = cellCntr + 1; 1057 } 1058 } 1059 } 1060 EG_free(lobjs); 1061 1062 CHKERRQ(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b)); 1063 CHKERRQ(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID)); 1064 } 1065 EG_free(fobjs); 1066 } 1067 1068 CHKERRQ(PetscHMapIDestroy(&edgeMap)); 1069 CHKERRQ(PetscHMapIDestroy(&bodyIndexMap)); 1070 CHKERRQ(PetscHMapIDestroy(&bodyVertexMap)); 1071 CHKERRQ(PetscHMapIDestroy(&bodyEdgeMap)); 1072 CHKERRQ(PetscHMapIDestroy(&bodyEdgeGlobalMap)); 1073 CHKERRQ(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 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1099 if (!rank) { 1100 // --------------------------------------------------------------------------------------------------- 1101 // Generate Petsc Plex from EGADSlite created Tessellation of geometry 1102 // --------------------------------------------------------------------------------------------------- 1103 1104 // Caluculate cell and vertex sizes 1105 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); 1106 1107 CHKERRQ(PetscHMapICreate(&pointIndexStartMap)); 1108 CHKERRQ(PetscHMapICreate(&triIndexStartMap)); 1109 CHKERRQ(PetscHMapICreate(&pTypeLabelMap)); 1110 CHKERRQ(PetscHMapICreate(&pIndexLabelMap)); 1111 CHKERRQ(PetscHMapICreate(&pBodyIndexLabelMap)); 1112 CHKERRQ(PetscHMapICreate(&triFaceIDLabelMap)); 1113 CHKERRQ(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 CHKERRQ(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); 1127 CHKERRQ(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound)); 1128 1129 if (!PISfound) CHKERRQ(PetscHMapISet(pointIndexStartMap, b, totalNumPoints)); 1130 if (!TISfound) CHKERRQ(PetscHMapISet(triIndexStartMap, b, totalNumTris)); 1131 1132 /* Calculate Tessellation parameters based on Bounding Box */ 1133 /* Get Bounding Box Dimensions of the BODY */ 1134 CHKERRQ(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 CHKERRQ(EG_makeTessBody(body, params, &tessArray[b])); 1145 1146 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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) 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 CHKERRQ(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 CHKERRQ(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); 1197 PetscCheckFalse(!PISfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b); 1198 CHKERRQ(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart)); 1199 1200 CHKERRQ(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 CHKERRQ(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 CHKERRQ(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 CHKERRQ(PetscHMapIFind(pTypeLabelMap, global-1+pointIndexStart, &PTLiter, &PTLfound)); 1230 CHKERRQ(PetscHMapIFind(pIndexLabelMap, global-1+pointIndexStart, &PILiter, &PILfound)); 1231 CHKERRQ(PetscHMapIFind(pBodyIndexLabelMap, global-1+pointIndexStart, &PBLiter, &PBLfound)); 1232 1233 if (!PTLfound) CHKERRQ(PetscHMapISet(pTypeLabelMap, global-1+pointIndexStart, ptype[p])); 1234 if (!PILfound) CHKERRQ(PetscHMapISet(pIndexLabelMap, global-1+pointIndexStart, pindex[p])); 1235 if (!PBLfound) CHKERRQ(PetscHMapISet(pBodyIndexLabelMap, global-1+pointIndexStart, b)); 1236 1237 if (ptype[p] < 0) CHKERRQ(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 CHKERRQ(EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 0], &global)); 1246 cells[(counter*3) +0] = global-1+pointIndexStart; 1247 1248 CHKERRQ(EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 1], &globalA)); 1249 cells[(counter*3) +1] = globalA-1+pointIndexStart; 1250 1251 CHKERRQ(EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 2], &globalB)); 1252 cells[(counter*3) +2] = globalB-1+pointIndexStart; 1253 1254 CHKERRQ(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound)); 1255 CHKERRQ(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound)); 1256 1257 if (!TFLfound) CHKERRQ(PetscHMapISet(triFaceIDLabelMap, counter, fID)); 1258 if (!TBLfound) CHKERRQ(PetscHMapISet(triBodyIDLabelMap, counter, b)); 1259 1260 counter += 1; 1261 } 1262 } 1263 EG_free(fobjs); 1264 } 1265 } 1266 1267 //Build DMPlex 1268 CHKERRQ(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); 1269 CHKERRQ(PetscFree2(coords, cells)); 1270 1271 // Embed EGADS model in DM 1272 { 1273 PetscContainer modelObj, contextObj; 1274 1275 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); 1276 CHKERRQ(PetscContainerSetPointer(modelObj, model)); 1277 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj)); 1278 CHKERRQ(PetscContainerDestroy(&modelObj)); 1279 1280 CHKERRQ(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); 1281 CHKERRQ(PetscContainerSetPointer(contextObj, context)); 1282 CHKERRQ(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); 1283 CHKERRQ(PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj)); 1284 CHKERRQ(PetscContainerDestroy(&contextObj)); 1285 } 1286 1287 // Label Points 1288 CHKERRQ(DMCreateLabel(dm, "EGADS Body ID")); 1289 CHKERRQ(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 1290 CHKERRQ(DMCreateLabel(dm, "EGADS Face ID")); 1291 CHKERRQ(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 1292 CHKERRQ(DMCreateLabel(dm, "EGADS Edge ID")); 1293 CHKERRQ(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 1294 CHKERRQ(DMCreateLabel(dm, "EGADS Vertex ID")); 1295 CHKERRQ(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 CHKERRQ(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd)); 1301 CHKERRQ(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd)); 1302 CHKERRQ(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 CHKERRQ(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound)); 1312 PetscCheckFalse(!PTLfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n); 1313 CHKERRQ(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal)); 1314 1315 CHKERRQ(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound)); 1316 PetscCheckFalse(!PILfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n); 1317 CHKERRQ(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal)); 1318 1319 CHKERRQ(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound)); 1320 PetscCheckFalse(!PBLfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n); 1321 CHKERRQ(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal)); 1322 1323 CHKERRQ(DMLabelSetValue(bodyLabel, n, pBodyVal)); 1324 if (pTypeVal == 0) CHKERRQ(DMLabelSetValue(vertexLabel, n, pIndexVal)); 1325 if (pTypeVal > 0) CHKERRQ(DMLabelSetValue(edgeLabel, n, pIndexVal)); 1326 if (pTypeVal < 0) CHKERRQ(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 CHKERRQ(DMPlexGetCone(dm, e, &cone)); 1334 CHKERRQ(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end? 1335 CHKERRQ(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0)); 1336 CHKERRQ(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1)); 1337 CHKERRQ(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0)); 1338 CHKERRQ(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1)); 1339 CHKERRQ(DMLabelGetValue(faceLabel, cone[0], &faceID_0)); 1340 CHKERRQ(DMLabelGetValue(faceLabel, cone[1], &faceID_1)); 1341 1342 CHKERRQ(DMLabelSetValue(bodyLabel, e, bodyID_0)); 1343 1344 if (edgeID_0 == edgeID_1) CHKERRQ(DMLabelSetValue(edgeLabel, e, edgeID_0)); 1345 else if (vertexID_0 > 0 && edgeID_1 > 0) CHKERRQ(DMLabelSetValue(edgeLabel, e, edgeID_1)); 1346 else if (vertexID_1 > 0 && edgeID_0 > 0) CHKERRQ(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 CHKERRQ(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound)); 1359 PetscCheckFalse(!TFLfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f); 1360 CHKERRQ(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal)); 1361 1362 CHKERRQ(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound)); 1363 PetscCheckFalse(!TBLfound,PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f); 1364 CHKERRQ(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal)); 1365 1366 CHKERRQ(DMLabelSetValue(bodyLabel, f, triBodyVal)); 1367 CHKERRQ(DMLabelSetValue(faceLabel, f, triFaceVal)); 1368 1369 /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */ 1370 CHKERRQ(DMPlexGetCone(dm, f, &cone)); 1371 1372 for (int jj = 0; jj < 3; ++jj) { 1373 CHKERRQ(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0)); 1374 1375 if (edgeID_0 < 0) { 1376 CHKERRQ(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal)); 1377 CHKERRQ(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 CHKERRQ(PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj)); 1423 if (!modelObj) PetscFunctionReturn(0); 1424 CHKERRQ(DMGetCoordinateDim(dm, &cdim)); 1425 CHKERRQ(DMGetCoordinateDM(dm, &cdm)); 1426 CHKERRQ(DMGetCoordinatesLocal(dm, &coordinates)); 1427 CHKERRQ(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 1428 CHKERRQ(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 1429 CHKERRQ(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); 1430 CHKERRQ(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); 1431 1432 CHKERRQ(PetscContainerGetPointer(modelObj, (void **) &model)); 1433 CHKERRQ(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); 1434 1435 CHKERRQ(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 1436 CHKERRQ(VecGetArrayWrite(coordinates, &coords)); 1437 for (v = vStart; v < vEnd; ++v) { 1438 PetscScalar *vcoords; 1439 1440 CHKERRQ(DMLabelGetValue(bodyLabel, v, &bodyID)); 1441 CHKERRQ(DMLabelGetValue(faceLabel, v, &faceID)); 1442 CHKERRQ(DMLabelGetValue(edgeLabel, v, &edgeID)); 1443 CHKERRQ(DMLabelGetValue(vertexLabel, v, &vertexID)); 1444 1445 PetscCheckFalse(bodyID >= Nb,PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb); 1446 body = bodies[bodyID]; 1447 1448 CHKERRQ(DMPlexPointLocalRef(cdm, v, coords, (void *) &vcoords)); 1449 if (edgeID > 0) { 1450 /* Snap to EDGE at nearest location */ 1451 double params[1]; 1452 CHKERRQ(EG_objectBodyTopo(body, EDGE, edgeID, &edge)); 1453 CHKERRQ(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 CHKERRQ(EG_objectBodyTopo(body, FACE, faceID, &face)); 1459 CHKERRQ(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 CHKERRQ(VecRestoreArrayWrite(coordinates, &coords)); 1464 /* Clear out global coordinates */ 1465 CHKERRQ(VecDestroy(&dm->coordinates)); 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 CHKERRQ(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL)); 1497 CHKERRQ(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL)); 1498 CHKERRQ(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL)); 1499 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 1500 #if defined(PETSC_HAVE_EGADS) 1501 if (rank == 0) { 1502 1503 CHKERRQ(EG_open(&context)); 1504 CHKERRQ(EG_loadModel(context, 0, filename, &model)); 1505 if (printModel) CHKERRQ(DMPlexEGADSPrintModel_Internal(model)); 1506 1507 } 1508 if (tessModel) CHKERRQ(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm)); 1509 else if (newModel) CHKERRQ(DMPlexCreateEGADS(comm, context, model, dm)); 1510 else CHKERRQ(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