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 22 /*@ 23 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. 24 25 Not collective 26 27 Input Parameters: 28 + dm - The DMPlex object 29 . p - The mesh point 30 - mcoords - A coordinate point lying on the mesh point 31 32 Output Parameter: 33 . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point 34 35 Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS. 36 37 Level: intermediate 38 39 .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform() 40 @*/ 41 PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[]) 42 { 43 #ifdef PETSC_HAVE_EGADS 44 DM cdm; 45 DMLabel bodyLabel, faceLabel, edgeLabel; 46 PetscContainer modelObj; 47 PetscInt bodyID, faceID, edgeID; 48 ego *bodies; 49 ego model, geom, body, obj; 50 /* result has to hold derviatives, along with the value */ 51 double params[3], result[18], paramsV[16*3], resultV[16*3]; 52 int Nb, oclass, mtype, *senses; 53 Vec coordinatesLocal; 54 PetscScalar *coords = NULL; 55 PetscInt Nv, v, Np = 0, pm; 56 #endif 57 PetscInt dE, d; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 ierr = DMGetCoordinateDim(dm, &dE);CHKERRQ(ierr); 62 #ifdef PETSC_HAVE_EGADS 63 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 64 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 65 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 66 if (!bodyLabel || !faceLabel || !edgeLabel) { 67 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 68 PetscFunctionReturn(0); 69 } 70 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 71 ierr = DMGetCoordinatesLocal(dm, &coordinatesLocal);CHKERRQ(ierr); 72 ierr = DMLabelGetValue(bodyLabel, p, &bodyID);CHKERRQ(ierr); 73 ierr = DMLabelGetValue(faceLabel, p, &faceID);CHKERRQ(ierr); 74 ierr = DMLabelGetValue(edgeLabel, p, &edgeID);CHKERRQ(ierr); 75 ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr); 76 ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr); 77 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 78 if (bodyID < 0 || bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb); 79 body = bodies[bodyID]; 80 81 if (edgeID >= 0) {ierr = EG_objectBodyTopo(body, EDGE, edgeID, &obj);CHKERRQ(ierr); Np = 1;} 82 else if (faceID >= 0) {ierr = EG_objectBodyTopo(body, FACE, faceID, &obj);CHKERRQ(ierr); Np = 2;} 83 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D is not in edge or face label for EGADS", p); 84 /* Calculate parameters (t or u,v) for vertices */ 85 ierr = DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 86 Nv /= dE; 87 if (Nv == 1) { 88 ierr = DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 89 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 90 PetscFunctionReturn(0); 91 } 92 if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p); 93 for (v = 0; v < Nv; ++v) {ierr = EG_invEvaluate(obj, &coords[v*dE], ¶msV[v*3], &resultV[v*3]);CHKERRQ(ierr);} 94 ierr = DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);CHKERRQ(ierr); 95 /* Calculate parameters (t or u,v) for new vertex at edge midpoint */ 96 for (pm = 0; pm < Np; ++pm) { 97 params[pm] = 0.; 98 for (v = 0; v < Nv; ++v) {params[pm] += paramsV[v*3+pm];} 99 params[pm] /= Nv; 100 } 101 /* TODO Check 102 double range[4]; // [umin, umax, vmin, vmax] 103 int peri; 104 ierr = EG_getRange(face, range, &peri);CHKERRQ(ierr); 105 if ((paramsNew[0] < range[0]) || (paramsNew[0] > range[1]) || (paramsNew[1] < range[2]) || (paramsNew[1] > range[3])) SETERRQ(); 106 */ 107 /* Put coordinates for new vertex in result[] */ 108 ierr = EG_evaluate(obj, params, result);CHKERRQ(ierr); 109 for (d = 0; d < dE; ++d) gcoords[d] = result[d]; 110 #else 111 for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; 112 #endif 113 PetscFunctionReturn(0); 114 } 115 116 #if defined(PETSC_HAVE_EGADS) 117 static PetscErrorCode DMPlexEGADSPrintModel(ego model) 118 { 119 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 120 int oclass, mtype, *senses; 121 int Nb, b; 122 PetscErrorCode ierr; 123 124 PetscFunctionBeginUser; 125 /* test bodyTopo functions */ 126 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 127 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb);CHKERRQ(ierr); 128 129 for (b = 0; b < Nb; ++b) { 130 ego body = bodies[b]; 131 int id, Nsh, Nf, Nl, l, Ne, e, Nv, v; 132 133 /* Output Basic Model Topology */ 134 ierr = EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);CHKERRQ(ierr); 135 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh);CHKERRQ(ierr); 136 EG_free(objs); 137 138 ierr = EG_getBodyTopos(body, NULL, FACE, &Nf, &objs);CHKERRQ(ierr); 139 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf);CHKERRQ(ierr); 140 EG_free(objs); 141 142 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 143 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl);CHKERRQ(ierr); 144 145 ierr = EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs);CHKERRQ(ierr); 146 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne);CHKERRQ(ierr); 147 EG_free(objs); 148 149 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &objs);CHKERRQ(ierr); 150 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv);CHKERRQ(ierr); 151 EG_free(objs); 152 153 for (l = 0; l < Nl; ++l) { 154 ego loop = lobjs[l]; 155 156 id = EG_indexBodyTopo(body, loop); 157 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id);CHKERRQ(ierr); 158 159 /* Get EDGE info which associated with the current LOOP */ 160 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 161 162 for (e = 0; e < Ne; ++e) { 163 ego edge = objs[e]; 164 double range[4] = {0., 0., 0., 0.}; 165 double point[3] = {0., 0., 0.}; 166 int peri; 167 168 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 169 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d\n", id);CHKERRQ(ierr); 170 171 ierr = EG_getRange(objs[e], range, &peri);CHKERRQ(ierr); 172 ierr = PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]); 173 174 /* Get NODE info which associated with the current EDGE */ 175 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);CHKERRQ(ierr); 176 177 for (v = 0; v < Nv; ++v) { 178 ego vertex = nobjs[v]; 179 double limits[4]; 180 int dummy; 181 182 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 183 id = EG_indexBodyTopo(body, vertex); 184 ierr = PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);CHKERRQ(ierr); 185 ierr = PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]); 186 187 point[0] = point[0] + limits[0]; 188 point[1] = point[1] + limits[1]; 189 point[2] = point[2] + limits[2]; 190 } 191 } 192 } 193 EG_free(lobjs); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 static PetscErrorCode DMPlexEGADSDestroy_Private(void *context) 199 { 200 if (context) EG_close((ego) context); 201 return 0; 202 } 203 204 static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm) 205 { 206 DMLabel bodyLabel, faceLabel, edgeLabel; 207 PetscInt cStart, cEnd, c; 208 /* EGADSLite variables */ 209 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 210 int oclass, mtype, nbodies, *senses; 211 int b; 212 /* PETSc variables */ 213 DM dm; 214 PetscHMapI edgeMap = NULL; 215 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; 216 PetscInt *cells = NULL, *cone = NULL; 217 PetscReal *coords = NULL; 218 PetscMPIInt rank; 219 PetscErrorCode ierr; 220 221 PetscFunctionBegin; 222 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 223 if (!rank) { 224 /* --------------------------------------------------------------------------------------------------- 225 Generate Petsc Plex 226 Get all Nodes in model, record coordinates in a correctly formatted array 227 Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array 228 We need to uniformly refine the initial geometry to guarantee a valid mesh 229 */ 230 231 /* Calculate cell and vertex sizes */ 232 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);CHKERRQ(ierr); 233 ierr = PetscHMapICreate(&edgeMap);CHKERRQ(ierr); 234 numEdges = 0; 235 for (b = 0; b < nbodies; ++b) { 236 ego body = bodies[b]; 237 int id, Nl, l, Nv, v; 238 239 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 240 for (l = 0; l < Nl; ++l) { 241 ego loop = lobjs[l]; 242 int Ner = 0, Ne, e, Nc; 243 244 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 245 for (e = 0; e < Ne; ++e) { 246 ego edge = objs[e]; 247 int Nv, id; 248 PetscHashIter iter; 249 PetscBool found; 250 251 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 252 if (mtype == DEGENERATE) continue; 253 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 254 ierr = PetscHMapIFind(edgeMap, id-1, &iter, &found);CHKERRQ(ierr); 255 if (!found) {ierr = PetscHMapISet(edgeMap, id-1, numEdges++);CHKERRQ(ierr);} 256 ++Ner; 257 } 258 if (Ner == 2) {Nc = 2;} 259 else if (Ner == 3) {Nc = 4;} 260 else if (Ner == 4) {Nc = 8; ++numQuads;} 261 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner); 262 numCells += Nc; 263 newCells += Nc-1; 264 maxCorners = PetscMax(Ner*2+1, maxCorners); 265 } 266 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 267 for (v = 0; v < Nv; ++v) { 268 ego vertex = nobjs[v]; 269 270 id = EG_indexBodyTopo(body, vertex); 271 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 272 numVertices = PetscMax(id, numVertices); 273 } 274 EG_free(lobjs); 275 EG_free(nobjs); 276 } 277 ierr = PetscHMapIGetSize(edgeMap, &numEdges);CHKERRQ(ierr); 278 newVertices = numEdges + numQuads; 279 numVertices += newVertices; 280 ierr = PetscPrintf(PETSC_COMM_SELF, "\nPLEX Input Array Checkouts\n");CHKERRQ(ierr); 281 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);CHKERRQ(ierr); 282 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Vertices = %D (%D) \n", numVertices, newVertices);CHKERRQ(ierr); 283 284 dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 285 cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 286 numCorners = 3; /* Split cells into triangles */ 287 ierr = PetscMalloc3(numVertices*cdim, &coords, numCells*numCorners, &cells, maxCorners, &cone);CHKERRQ(ierr); 288 289 /* Get vertex coordinates */ 290 for (b = 0; b < nbodies; ++b) { 291 ego body = bodies[b]; 292 int id, Nv, v; 293 294 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 295 for (v = 0; v < Nv; ++v) { 296 ego vertex = nobjs[v]; 297 double limits[4]; 298 int dummy; 299 300 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 301 id = EG_indexBodyTopo(body, vertex);CHKERRQ(ierr); 302 coords[(id-1)*cdim+0] = limits[0]; 303 coords[(id-1)*cdim+1] = limits[1]; 304 coords[(id-1)*cdim+2] = limits[2]; 305 ierr = PetscPrintf(PETSC_COMM_SELF, " Node ID = %d (%d)\n", id, id-1); 306 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf) \n \n", coords[(id-1)*cdim+0], coords[(id-1)*cdim+1],coords[(id-1)*cdim+2]); 307 } 308 EG_free(nobjs); 309 } 310 ierr = PetscHMapIClear(edgeMap);CHKERRQ(ierr); 311 fOff = numVertices - newVertices + numEdges; 312 numEdges = 0; 313 numQuads = 0; 314 for (b = 0; b < nbodies; ++b) { 315 ego body = bodies[b]; 316 int Nl, l; 317 318 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 319 for (l = 0; l < Nl; ++l) { 320 ego loop = lobjs[l]; 321 int lid, Ner = 0, Ne, e; 322 323 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 324 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 325 for (e = 0; e < Ne; ++e) { 326 ego edge = objs[e]; 327 int eid, Nv; 328 PetscHashIter iter; 329 PetscBool found; 330 331 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 332 if (mtype == DEGENERATE) continue; 333 ++Ner; 334 eid = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 335 ierr = PetscHMapIFind(edgeMap, eid-1, &iter, &found);CHKERRQ(ierr); 336 if (!found) { 337 PetscInt v = numVertices - newVertices + numEdges; 338 double range[4], params[3] = {0., 0., 0.}, result[18]; 339 int periodic[2]; 340 341 ierr = PetscHMapISet(edgeMap, eid-1, numEdges++);CHKERRQ(ierr); 342 ierr = EG_getRange(edge, range, periodic);CHKERRQ(ierr); 343 params[0] = 0.5*(range[0] + range[1]); 344 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 345 coords[v*cdim+0] = result[0]; 346 coords[v*cdim+1] = result[1]; 347 coords[v*cdim+2] = result[2]; 348 ierr = PetscPrintf(PETSC_COMM_SELF, " Edge ID = %d (%D) \n", eid, v); 349 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf)\n", coords[v*cdim+0], coords[v*cdim+1],coords[v*cdim+2]); 350 params[0] = range[0]; 351 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 352 ierr = PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]); 353 params[0] = range[1]; 354 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 355 ierr = PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n\n", result[0], result[1], result[2]); 356 } 357 } 358 if (Ner == 4) { 359 PetscInt v = fOff + numQuads++; 360 ego *fobjs; 361 double range[4], params[3] = {0., 0., 0.}, result[18]; 362 int Nf, face, periodic[2]; 363 364 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 365 if (Nf != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1", lid-1, Nf); 366 face = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 367 ierr = EG_getRange(fobjs[0], range, periodic);CHKERRQ(ierr); 368 params[0] = 0.5*(range[0] + range[1]); 369 params[1] = 0.5*(range[2] + range[3]); 370 ierr = EG_evaluate(fobjs[0], params, result);CHKERRQ(ierr); 371 coords[v*cdim+0] = result[0]; 372 coords[v*cdim+1] = result[1]; 373 coords[v*cdim+2] = result[2]; 374 ierr = PetscPrintf(PETSC_COMM_SELF, " Face ID = %d (%D) \n", face-1, v); 375 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf) \n \n", coords[v*cdim+0], coords[v*cdim+1],coords[v*cdim+2]); 376 } 377 } 378 } 379 if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices); 380 CHKMEMQ; 381 382 /* Get cell vertices by traversing loops */ 383 numQuads = 0; 384 cOff = 0; 385 for (b = 0; b < nbodies; ++b) { 386 ego body = bodies[b]; 387 int id, Nl, l; 388 389 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 390 for (l = 0; l < Nl; ++l) { 391 ego loop = lobjs[l]; 392 int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t; 393 394 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 395 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d \n", lid);CHKERRQ(ierr); 396 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 397 398 for (e = 0; e < Ne; ++e) { 399 ego edge = objs[e]; 400 int points[3]; 401 int eid, Nv, v, tmp; 402 403 eid = EG_indexBodyTopo(body, edge); 404 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d \n", eid);CHKERRQ(ierr); 405 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 406 if (mtype == DEGENERATE) {ierr = PetscPrintf(PETSC_COMM_SELF, " EGDE %d is DEGENERATE \n", eid);CHKERRQ(ierr); continue;} 407 else {++Ner;} 408 if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv); 409 410 for (v = 0; v < Nv; ++v) { 411 ego vertex = nobjs[v]; 412 413 id = EG_indexBodyTopo(body, vertex); 414 points[v*2] = id-1; 415 } 416 { 417 PetscInt edgeNum; 418 419 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 420 points[1] = numVertices - newVertices + edgeNum; 421 } 422 /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */ 423 if (!nc) { 424 for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v]; 425 } else { 426 if (cone[nc-1] == points[0]) {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];} 427 else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];} 428 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];} 429 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];} 430 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); 431 } 432 } 433 if (nc != 2*Ner) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner); 434 if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;} 435 if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners); 436 /* Triangulate the loop */ 437 switch (Ner) { 438 case 2: /* Bi-Segment -> 2 triangles */ 439 Nt = 2; 440 cells[cOff*numCorners+0] = cone[0]; 441 cells[cOff*numCorners+1] = cone[1]; 442 cells[cOff*numCorners+2] = cone[2]; 443 ++cOff; 444 cells[cOff*numCorners+0] = cone[0]; 445 cells[cOff*numCorners+1] = cone[2]; 446 cells[cOff*numCorners+2] = cone[3]; 447 ++cOff; 448 break; 449 case 3: /* Triangle -> 4 triangles */ 450 Nt = 4; 451 cells[cOff*numCorners+0] = cone[0]; 452 cells[cOff*numCorners+1] = cone[1]; 453 cells[cOff*numCorners+2] = cone[5]; 454 ++cOff; 455 cells[cOff*numCorners+0] = cone[1]; 456 cells[cOff*numCorners+1] = cone[2]; 457 cells[cOff*numCorners+2] = cone[3]; 458 ++cOff; 459 cells[cOff*numCorners+0] = cone[5]; 460 cells[cOff*numCorners+1] = cone[3]; 461 cells[cOff*numCorners+2] = cone[4]; 462 ++cOff; 463 cells[cOff*numCorners+0] = cone[1]; 464 cells[cOff*numCorners+1] = cone[3]; 465 cells[cOff*numCorners+2] = cone[5]; 466 ++cOff; 467 break; 468 case 4: /* Quad -> 8 triangles */ 469 Nt = 8; 470 cells[cOff*numCorners+0] = cone[0]; 471 cells[cOff*numCorners+1] = cone[1]; 472 cells[cOff*numCorners+2] = cone[7]; 473 ++cOff; 474 cells[cOff*numCorners+0] = cone[1]; 475 cells[cOff*numCorners+1] = cone[2]; 476 cells[cOff*numCorners+2] = cone[3]; 477 ++cOff; 478 cells[cOff*numCorners+0] = cone[3]; 479 cells[cOff*numCorners+1] = cone[4]; 480 cells[cOff*numCorners+2] = cone[5]; 481 ++cOff; 482 cells[cOff*numCorners+0] = cone[5]; 483 cells[cOff*numCorners+1] = cone[6]; 484 cells[cOff*numCorners+2] = cone[7]; 485 ++cOff; 486 cells[cOff*numCorners+0] = cone[8]; 487 cells[cOff*numCorners+1] = cone[1]; 488 cells[cOff*numCorners+2] = cone[3]; 489 ++cOff; 490 cells[cOff*numCorners+0] = cone[8]; 491 cells[cOff*numCorners+1] = cone[3]; 492 cells[cOff*numCorners+2] = cone[5]; 493 ++cOff; 494 cells[cOff*numCorners+0] = cone[8]; 495 cells[cOff*numCorners+1] = cone[5]; 496 cells[cOff*numCorners+2] = cone[7]; 497 ++cOff; 498 cells[cOff*numCorners+0] = cone[8]; 499 cells[cOff*numCorners+1] = cone[7]; 500 cells[cOff*numCorners+2] = cone[1]; 501 ++cOff; 502 break; 503 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner); 504 } 505 for (t = 0; t < Nt; ++t) { 506 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t);CHKERRQ(ierr); 507 for (c = 0; c < numCorners; ++c) { 508 if (c > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 509 ierr = PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);CHKERRQ(ierr); 510 } 511 ierr = PetscPrintf(PETSC_COMM_SELF, ")\n");CHKERRQ(ierr); 512 } 513 } 514 EG_free(lobjs); 515 } 516 } 517 if (cOff != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells); 518 ierr = DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);CHKERRQ(ierr); 519 ierr = PetscFree3(coords, cells, cone);CHKERRQ(ierr); 520 /* Embed EGADS model in DM */ 521 { 522 PetscContainer modelObj, contextObj; 523 524 ierr = PetscContainerCreate(PETSC_COMM_SELF, &modelObj);CHKERRQ(ierr); 525 ierr = PetscContainerSetPointer(modelObj, model);CHKERRQ(ierr); 526 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 527 ierr = PetscContainerDestroy(&modelObj);CHKERRQ(ierr); 528 529 ierr = PetscContainerCreate(PETSC_COMM_SELF, &contextObj);CHKERRQ(ierr); 530 ierr = PetscContainerSetPointer(contextObj, context);CHKERRQ(ierr); 531 ierr = PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);CHKERRQ(ierr); 532 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);CHKERRQ(ierr); 533 ierr = PetscContainerDestroy(&contextObj);CHKERRQ(ierr); 534 } 535 /* Label points */ 536 ierr = DMCreateLabel(dm, "EGADS Body ID");CHKERRQ(ierr); 537 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 538 ierr = DMCreateLabel(dm, "EGADS Face ID");CHKERRQ(ierr); 539 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 540 ierr = DMCreateLabel(dm, "EGADS Edge ID");CHKERRQ(ierr); 541 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 542 cOff = 0; 543 for (b = 0; b < nbodies; ++b) { 544 ego body = bodies[b]; 545 int id, Nl, l; 546 547 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 548 for (l = 0; l < Nl; ++l) { 549 ego loop = lobjs[l]; 550 ego *fobjs; 551 int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t; 552 553 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 554 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 555 if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);CHKERRQ(ierr); 556 fid = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 557 EG_free(fobjs); 558 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 559 for (e = 0; e < Ne; ++e) { 560 ego edge = objs[e]; 561 int eid, Nv, v; 562 PetscInt points[3], support[2], numEdges, edgeNum; 563 const PetscInt *edges; 564 565 eid = EG_indexBodyTopo(body, edge); 566 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 567 if (mtype == DEGENERATE) continue; 568 else ++Ner; 569 for (v = 0; v < Nv; ++v) { 570 ego vertex = nobjs[v]; 571 572 id = EG_indexBodyTopo(body, vertex); 573 ierr = DMLabelSetValue(edgeLabel, numCells + id-1, eid);CHKERRQ(ierr); 574 points[v*2] = numCells + id-1; 575 } 576 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 577 points[1] = numCells + numVertices - newVertices + edgeNum; 578 579 ierr = DMLabelSetValue(edgeLabel, points[1], eid);CHKERRQ(ierr); 580 support[0] = points[0]; 581 support[1] = points[1]; 582 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 583 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); 584 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 585 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 586 support[0] = points[1]; 587 support[1] = points[2]; 588 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 589 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); 590 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 591 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 592 } 593 switch (Ner) { 594 case 2: Nt = 2;break; 595 case 3: Nt = 4;break; 596 case 4: Nt = 8;break; 597 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner); 598 } 599 for (t = 0; t < Nt; ++t) { 600 ierr = DMLabelSetValue(bodyLabel, cOff+t, b);CHKERRQ(ierr); 601 ierr = DMLabelSetValue(faceLabel, cOff+t, fid);CHKERRQ(ierr); 602 } 603 cOff += Nt; 604 } 605 EG_free(lobjs); 606 } 607 ierr = PetscHMapIDestroy(&edgeMap);CHKERRQ(ierr); 608 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 609 for (c = cStart; c < cEnd; ++c) { 610 PetscInt *closure = NULL; 611 PetscInt clSize, cl, bval, fval; 612 613 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 614 ierr = DMLabelGetValue(bodyLabel, c, &bval);CHKERRQ(ierr); 615 ierr = DMLabelGetValue(faceLabel, c, &fval);CHKERRQ(ierr); 616 for (cl = 0; cl < clSize*2; cl += 2) { 617 ierr = DMLabelSetValue(bodyLabel, closure[cl], bval);CHKERRQ(ierr); 618 ierr = DMLabelSetValue(faceLabel, closure[cl], fval);CHKERRQ(ierr); 619 } 620 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 621 } 622 *newdm = dm; 623 PetscFunctionReturn(0); 624 } 625 #endif 626 627 /*@C 628 DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADSLite file. 629 630 Collective 631 632 Input Parameters: 633 + comm - The MPI communicator 634 - filename - The name of the ExodusII file 635 636 Output Parameter: 637 . dm - The DM object representing the mesh 638 639 Level: beginner 640 641 .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS() 642 @*/ 643 PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm) 644 { 645 PetscMPIInt rank; 646 #if defined(PETSC_HAVE_EGADS) 647 ego context= NULL, model = NULL; 648 #endif 649 PetscBool printModel; 650 PetscErrorCode ierr; 651 652 PetscFunctionBegin; 653 PetscValidCharPointer(filename, 2); 654 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);CHKERRQ(ierr); 655 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 656 #if defined(PETSC_HAVE_EGADS) 657 if (!rank) { 658 659 ierr = EG_open(&context);CHKERRQ(ierr); 660 ierr = EG_loadModel(context, 0, filename, &model);CHKERRQ(ierr); 661 if (printModel) {ierr = DMPlexEGADSPrintModel(model);CHKERRQ(ierr);} 662 663 } 664 ierr = DMPlexCreateEGADS(comm, context, model, dm);CHKERRQ(ierr); 665 #else 666 SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads"); 667 #endif 668 PetscFunctionReturn(0); 669 } 670