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