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