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 if (bodyID >= Nb) SETERRQ2(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 if (Nv > 16) SETERRQ2(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 if ((params[0] < range[0]) || (params[0] > range[1])) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D had bad interpolation", p); 87 if (Np > 1 && ((params[1] < range[2]) || (params[1] > range[3]))) SETERRQ1(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 - 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. 109 110 Level: intermediate 111 112 .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform() 113 @*/ 114 PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[]) 115 { 116 PetscInt dE, d; 117 PetscErrorCode ierr; 118 119 PetscFunctionBeginHot; 120 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 121 #ifdef PETSC_HAVE_EGADS 122 { 123 DM_Plex *plex = (DM_Plex *) dm->data; 124 DMLabel bodyLabel, faceLabel, edgeLabel; 125 PetscInt bodyID, faceID, edgeID; 126 PetscContainer modelObj; 127 ego model; 128 PetscBool islite = PETSC_FALSE; 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 SETERRQ1(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 if (Nf != 1) SETERRQ3(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 if (numEdges + numQuads != newVertices) SETERRQ2(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 if (Nv != 2) SETERRQ2(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 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); 476 } 477 } 478 if (nc != 2*Ner) SETERRQ2(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 if (nc > maxCorners) SETERRQ2(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: SETERRQ2(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 if (cOff != numCells) SETERRQ2(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 = PetscInfo2(dm, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);CHKERRQ(ierr); 568 ierr = PetscInfo2(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 if (Nf > 1) SETERRQ2(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 if (numEdges != 1) SETERRQ3(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 if (numEdges != 1) SETERRQ3(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: SETERRQ1(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 if (!BVfound) SETERRQ1(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 if (!BEfound) SETERRQ1(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 if (!BEGfound) SETERRQ1(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 if (!EMfound) SETERRQ1(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 if (!BFfound) SETERRQ1(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 if (Nl > 1) SETERRQ1(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 if (!EMfound) SETERRQ3(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 = PetscInfo1(dm, " Total Number of Unique Cells = %D \n", numCells);CHKERRQ(ierr); 947 ierr = PetscInfo1(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 if (!BVfound) SETERRQ1(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 if (!BEfound) SETERRQ1(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 if (!BFfound) SETERRQ1(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 if (!BEGfound) SETERRQ1(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 if (Nl > 1) SETERRQ2(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 if (!EMfound) SETERRQ3(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 if (!PISfound) SETERRQ1(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 if (!PTLfound) SETERRQ1(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 if (!PILfound) SETERRQ1(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 if (!PBLfound) SETERRQ1(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 if (!TFLfound) SETERRQ1(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 if (!TBLfound) SETERRQ1(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 if (bodyID >= Nb) SETERRQ2(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