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 double params[3] = {0., 0., 0.}; 167 double result[18]; 168 int peri; 169 170 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 171 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e);CHKERRQ(ierr); 172 173 ierr = EG_getRange(edge, range, &peri);CHKERRQ(ierr); 174 ierr = PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]); 175 176 /* Get NODE info which associated with the current EDGE */ 177 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);CHKERRQ(ierr); 178 if (mtype == DEGENERATE) { 179 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id);CHKERRQ(ierr); 180 } else { 181 params[0] = range[0]; 182 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 183 ierr = PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]); 184 params[0] = range[1]; 185 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 186 ierr = PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]); 187 } 188 189 for (v = 0; v < Nv; ++v) { 190 ego vertex = nobjs[v]; 191 double limits[4]; 192 int dummy; 193 194 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 195 id = EG_indexBodyTopo(body, vertex); 196 ierr = PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);CHKERRQ(ierr); 197 ierr = PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]); 198 199 point[0] = point[0] + limits[0]; 200 point[1] = point[1] + limits[1]; 201 point[2] = point[2] + limits[2]; 202 } 203 } 204 } 205 EG_free(lobjs); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 static PetscErrorCode DMPlexEGADSDestroy_Private(void *context) 211 { 212 if (context) EG_close((ego) context); 213 return 0; 214 } 215 216 static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm) 217 { 218 DMLabel bodyLabel, faceLabel, edgeLabel; 219 PetscInt cStart, cEnd, c; 220 /* EGADSLite variables */ 221 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 222 int oclass, mtype, nbodies, *senses; 223 int b; 224 /* PETSc variables */ 225 DM dm; 226 PetscHMapI edgeMap = NULL; 227 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; 228 PetscInt *cells = NULL, *cone = NULL; 229 PetscReal *coords = NULL; 230 PetscMPIInt rank; 231 PetscErrorCode ierr; 232 233 PetscFunctionBegin; 234 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 235 if (!rank) { 236 const PetscInt debug = 0; 237 238 /* --------------------------------------------------------------------------------------------------- 239 Generate Petsc Plex 240 Get all Nodes in model, record coordinates in a correctly formatted array 241 Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array 242 We need to uniformly refine the initial geometry to guarantee a valid mesh 243 */ 244 245 /* Calculate cell and vertex sizes */ 246 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);CHKERRQ(ierr); 247 ierr = PetscHMapICreate(&edgeMap);CHKERRQ(ierr); 248 numEdges = 0; 249 for (b = 0; b < nbodies; ++b) { 250 ego body = bodies[b]; 251 int id, Nl, l, Nv, v; 252 253 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 254 for (l = 0; l < Nl; ++l) { 255 ego loop = lobjs[l]; 256 int Ner = 0, Ne, e, Nc; 257 258 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 259 for (e = 0; e < Ne; ++e) { 260 ego edge = objs[e]; 261 int Nv, id; 262 PetscHashIter iter; 263 PetscBool found; 264 265 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 266 if (mtype == DEGENERATE) continue; 267 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 268 ierr = PetscHMapIFind(edgeMap, id-1, &iter, &found);CHKERRQ(ierr); 269 if (!found) {ierr = PetscHMapISet(edgeMap, id-1, numEdges++);CHKERRQ(ierr);} 270 ++Ner; 271 } 272 if (Ner == 2) {Nc = 2;} 273 else if (Ner == 3) {Nc = 4;} 274 else if (Ner == 4) {Nc = 8; ++numQuads;} 275 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner); 276 numCells += Nc; 277 newCells += Nc-1; 278 maxCorners = PetscMax(Ner*2+1, maxCorners); 279 } 280 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 281 for (v = 0; v < Nv; ++v) { 282 ego vertex = nobjs[v]; 283 284 id = EG_indexBodyTopo(body, vertex); 285 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 286 numVertices = PetscMax(id, numVertices); 287 } 288 EG_free(lobjs); 289 EG_free(nobjs); 290 } 291 ierr = PetscHMapIGetSize(edgeMap, &numEdges);CHKERRQ(ierr); 292 newVertices = numEdges + numQuads; 293 numVertices += newVertices; 294 295 dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 296 cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 297 numCorners = 3; /* Split cells into triangles */ 298 ierr = PetscMalloc3(numVertices*cdim, &coords, numCells*numCorners, &cells, maxCorners, &cone);CHKERRQ(ierr); 299 300 /* Get vertex coordinates */ 301 for (b = 0; b < nbodies; ++b) { 302 ego body = bodies[b]; 303 int id, Nv, v; 304 305 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 306 for (v = 0; v < Nv; ++v) { 307 ego vertex = nobjs[v]; 308 double limits[4]; 309 int dummy; 310 311 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 312 id = EG_indexBodyTopo(body, vertex);CHKERRQ(ierr); 313 coords[(id-1)*cdim+0] = limits[0]; 314 coords[(id-1)*cdim+1] = limits[1]; 315 coords[(id-1)*cdim+2] = limits[2]; 316 } 317 EG_free(nobjs); 318 } 319 ierr = PetscHMapIClear(edgeMap);CHKERRQ(ierr); 320 fOff = numVertices - newVertices + numEdges; 321 numEdges = 0; 322 numQuads = 0; 323 for (b = 0; b < nbodies; ++b) { 324 ego body = bodies[b]; 325 int Nl, l; 326 327 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 328 for (l = 0; l < Nl; ++l) { 329 ego loop = lobjs[l]; 330 int lid, Ner = 0, Ne, e; 331 332 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 333 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 334 for (e = 0; e < Ne; ++e) { 335 ego edge = objs[e]; 336 int eid, Nv; 337 PetscHashIter iter; 338 PetscBool found; 339 340 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 341 if (mtype == DEGENERATE) continue; 342 ++Ner; 343 eid = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 344 ierr = PetscHMapIFind(edgeMap, eid-1, &iter, &found);CHKERRQ(ierr); 345 if (!found) { 346 PetscInt v = numVertices - newVertices + numEdges; 347 double range[4], params[3] = {0., 0., 0.}, result[18]; 348 int periodic[2]; 349 350 ierr = PetscHMapISet(edgeMap, eid-1, numEdges++);CHKERRQ(ierr); 351 ierr = EG_getRange(edge, range, periodic);CHKERRQ(ierr); 352 params[0] = 0.5*(range[0] + range[1]); 353 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 354 coords[v*cdim+0] = result[0]; 355 coords[v*cdim+1] = result[1]; 356 coords[v*cdim+2] = result[2]; 357 } 358 } 359 if (Ner == 4) { 360 PetscInt v = fOff + numQuads++; 361 ego *fobjs, face; 362 double range[4], params[3] = {0., 0., 0.}, result[18]; 363 int Nf, fid, periodic[2]; 364 365 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 366 face = fobjs[0]; 367 fid = EG_indexBodyTopo(body, face);CHKERRQ(ierr); 368 if (Nf != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid); 369 ierr = EG_getRange(face, range, periodic);CHKERRQ(ierr); 370 params[0] = 0.5*(range[0] + range[1]); 371 params[1] = 0.5*(range[2] + range[3]); 372 ierr = EG_evaluate(face, params, result);CHKERRQ(ierr); 373 coords[v*cdim+0] = result[0]; 374 coords[v*cdim+1] = result[1]; 375 coords[v*cdim+2] = result[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 = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 396 397 for (e = 0; e < Ne; ++e) { 398 ego edge = objs[e]; 399 int points[3]; 400 int eid, Nv, v, tmp; 401 402 eid = EG_indexBodyTopo(body, edge); 403 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 404 if (mtype == DEGENERATE) continue; 405 else ++Ner; 406 if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv); 407 408 for (v = 0; v < Nv; ++v) { 409 ego vertex = nobjs[v]; 410 411 id = EG_indexBodyTopo(body, vertex); 412 points[v*2] = id-1; 413 } 414 { 415 PetscInt edgeNum; 416 417 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 418 points[1] = numVertices - newVertices + edgeNum; 419 } 420 /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */ 421 if (!nc) { 422 for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v]; 423 } else { 424 if (cone[nc-1] == points[0]) {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];} 425 else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];} 426 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];} 427 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];} 428 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); 429 } 430 } 431 if (nc != 2*Ner) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner); 432 if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;} 433 if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners); 434 /* Triangulate the loop */ 435 switch (Ner) { 436 case 2: /* Bi-Segment -> 2 triangles */ 437 Nt = 2; 438 cells[cOff*numCorners+0] = cone[0]; 439 cells[cOff*numCorners+1] = cone[1]; 440 cells[cOff*numCorners+2] = cone[2]; 441 ++cOff; 442 cells[cOff*numCorners+0] = cone[0]; 443 cells[cOff*numCorners+1] = cone[2]; 444 cells[cOff*numCorners+2] = cone[3]; 445 ++cOff; 446 break; 447 case 3: /* Triangle -> 4 triangles */ 448 Nt = 4; 449 cells[cOff*numCorners+0] = cone[0]; 450 cells[cOff*numCorners+1] = cone[1]; 451 cells[cOff*numCorners+2] = cone[5]; 452 ++cOff; 453 cells[cOff*numCorners+0] = cone[1]; 454 cells[cOff*numCorners+1] = cone[2]; 455 cells[cOff*numCorners+2] = cone[3]; 456 ++cOff; 457 cells[cOff*numCorners+0] = cone[5]; 458 cells[cOff*numCorners+1] = cone[3]; 459 cells[cOff*numCorners+2] = cone[4]; 460 ++cOff; 461 cells[cOff*numCorners+0] = cone[1]; 462 cells[cOff*numCorners+1] = cone[3]; 463 cells[cOff*numCorners+2] = cone[5]; 464 ++cOff; 465 break; 466 case 4: /* Quad -> 8 triangles */ 467 Nt = 8; 468 cells[cOff*numCorners+0] = cone[0]; 469 cells[cOff*numCorners+1] = cone[1]; 470 cells[cOff*numCorners+2] = cone[7]; 471 ++cOff; 472 cells[cOff*numCorners+0] = cone[1]; 473 cells[cOff*numCorners+1] = cone[2]; 474 cells[cOff*numCorners+2] = cone[3]; 475 ++cOff; 476 cells[cOff*numCorners+0] = cone[3]; 477 cells[cOff*numCorners+1] = cone[4]; 478 cells[cOff*numCorners+2] = cone[5]; 479 ++cOff; 480 cells[cOff*numCorners+0] = cone[5]; 481 cells[cOff*numCorners+1] = cone[6]; 482 cells[cOff*numCorners+2] = cone[7]; 483 ++cOff; 484 cells[cOff*numCorners+0] = cone[8]; 485 cells[cOff*numCorners+1] = cone[1]; 486 cells[cOff*numCorners+2] = cone[3]; 487 ++cOff; 488 cells[cOff*numCorners+0] = cone[8]; 489 cells[cOff*numCorners+1] = cone[3]; 490 cells[cOff*numCorners+2] = cone[5]; 491 ++cOff; 492 cells[cOff*numCorners+0] = cone[8]; 493 cells[cOff*numCorners+1] = cone[5]; 494 cells[cOff*numCorners+2] = cone[7]; 495 ++cOff; 496 cells[cOff*numCorners+0] = cone[8]; 497 cells[cOff*numCorners+1] = cone[7]; 498 cells[cOff*numCorners+2] = cone[1]; 499 ++cOff; 500 break; 501 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner); 502 } 503 if (debug) { 504 for (t = 0; t < Nt; ++t) { 505 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t);CHKERRQ(ierr); 506 for (c = 0; c < numCorners; ++c) { 507 if (c > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");CHKERRQ(ierr);} 508 ierr = PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);CHKERRQ(ierr); 509 } 510 ierr = PetscPrintf(PETSC_COMM_SELF, ")\n");CHKERRQ(ierr); 511 } 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 ierr = PetscInfo2(dm, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);CHKERRQ(ierr); 521 ierr = PetscInfo2(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices);CHKERRQ(ierr); 522 /* Embed EGADS model in DM */ 523 { 524 PetscContainer modelObj, contextObj; 525 526 ierr = PetscContainerCreate(PETSC_COMM_SELF, &modelObj);CHKERRQ(ierr); 527 ierr = PetscContainerSetPointer(modelObj, model);CHKERRQ(ierr); 528 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 529 ierr = PetscContainerDestroy(&modelObj);CHKERRQ(ierr); 530 531 ierr = PetscContainerCreate(PETSC_COMM_SELF, &contextObj);CHKERRQ(ierr); 532 ierr = PetscContainerSetPointer(contextObj, context);CHKERRQ(ierr); 533 ierr = PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);CHKERRQ(ierr); 534 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);CHKERRQ(ierr); 535 ierr = PetscContainerDestroy(&contextObj);CHKERRQ(ierr); 536 } 537 /* Label points */ 538 ierr = DMCreateLabel(dm, "EGADS Body ID");CHKERRQ(ierr); 539 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 540 ierr = DMCreateLabel(dm, "EGADS Face ID");CHKERRQ(ierr); 541 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 542 ierr = DMCreateLabel(dm, "EGADS Edge ID");CHKERRQ(ierr); 543 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 544 cOff = 0; 545 for (b = 0; b < nbodies; ++b) { 546 ego body = bodies[b]; 547 int id, Nl, l; 548 549 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 550 for (l = 0; l < Nl; ++l) { 551 ego loop = lobjs[l]; 552 ego *fobjs; 553 int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t; 554 555 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 556 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 557 if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);CHKERRQ(ierr); 558 fid = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 559 EG_free(fobjs); 560 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 561 for (e = 0; e < Ne; ++e) { 562 ego edge = objs[e]; 563 int eid, Nv, v; 564 PetscInt points[3], support[2], numEdges, edgeNum; 565 const PetscInt *edges; 566 567 eid = EG_indexBodyTopo(body, edge); 568 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 569 if (mtype == DEGENERATE) continue; 570 else ++Ner; 571 for (v = 0; v < Nv; ++v) { 572 ego vertex = nobjs[v]; 573 574 id = EG_indexBodyTopo(body, vertex); 575 ierr = DMLabelSetValue(edgeLabel, numCells + id-1, eid);CHKERRQ(ierr); 576 points[v*2] = numCells + id-1; 577 } 578 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 579 points[1] = numCells + numVertices - newVertices + edgeNum; 580 581 ierr = DMLabelSetValue(edgeLabel, points[1], eid);CHKERRQ(ierr); 582 support[0] = points[0]; 583 support[1] = points[1]; 584 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 585 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); 586 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 587 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 588 support[0] = points[1]; 589 support[1] = points[2]; 590 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 591 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); 592 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 593 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 594 } 595 switch (Ner) { 596 case 2: Nt = 2;break; 597 case 3: Nt = 4;break; 598 case 4: Nt = 8;break; 599 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner); 600 } 601 for (t = 0; t < Nt; ++t) { 602 ierr = DMLabelSetValue(bodyLabel, cOff+t, b);CHKERRQ(ierr); 603 ierr = DMLabelSetValue(faceLabel, cOff+t, fid);CHKERRQ(ierr); 604 } 605 cOff += Nt; 606 } 607 EG_free(lobjs); 608 } 609 ierr = PetscHMapIDestroy(&edgeMap);CHKERRQ(ierr); 610 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 611 for (c = cStart; c < cEnd; ++c) { 612 PetscInt *closure = NULL; 613 PetscInt clSize, cl, bval, fval; 614 615 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 616 ierr = DMLabelGetValue(bodyLabel, c, &bval);CHKERRQ(ierr); 617 ierr = DMLabelGetValue(faceLabel, c, &fval);CHKERRQ(ierr); 618 for (cl = 0; cl < clSize*2; cl += 2) { 619 ierr = DMLabelSetValue(bodyLabel, closure[cl], bval);CHKERRQ(ierr); 620 ierr = DMLabelSetValue(faceLabel, closure[cl], fval);CHKERRQ(ierr); 621 } 622 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 623 } 624 *newdm = dm; 625 PetscFunctionReturn(0); 626 } 627 #endif 628 629 /*@C 630 DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADSLite file. 631 632 Collective 633 634 Input Parameters: 635 + comm - The MPI communicator 636 - filename - The name of the ExodusII file 637 638 Output Parameter: 639 . dm - The DM object representing the mesh 640 641 Level: beginner 642 643 .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS() 644 @*/ 645 PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm) 646 { 647 PetscMPIInt rank; 648 #if defined(PETSC_HAVE_EGADS) 649 ego context= NULL, model = NULL; 650 #endif 651 PetscBool printModel = PETSC_FALSE; 652 PetscErrorCode ierr; 653 654 PetscFunctionBegin; 655 PetscValidCharPointer(filename, 2); 656 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);CHKERRQ(ierr); 657 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 658 #if defined(PETSC_HAVE_EGADS) 659 if (!rank) { 660 661 ierr = EG_open(&context);CHKERRQ(ierr); 662 ierr = EG_loadModel(context, 0, filename, &model);CHKERRQ(ierr); 663 if (printModel) {ierr = DMPlexEGADSPrintModel(model);CHKERRQ(ierr);} 664 665 } 666 ierr = DMPlexCreateEGADS(comm, context, model, dm);CHKERRQ(ierr); 667 PetscFunctionReturn(0); 668 #else 669 SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads"); 670 #endif 671 } 672