#include /*I "petscdmplex.h" I*/ #include #ifdef PETSC_HAVE_EGADS #include #endif /* We need to understand how to natively parse STEP files. There seems to be only one open-source implementation of the STEP parser contained in the OpenCASCADE package. It is enough to make a strong man weep: https://github.com/tpaviot/oce/tree/master/src/STEPControl The STEP, and inner EXPRESS, formats are ISO standards, so they are documented https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files http://stepmod.sourceforge.net/express_model_spec/ but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants. */ #ifdef PETSC_HAVE_EGADS PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]); PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]); PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM dm, PetscInt p, ego model, PetscInt bodyID, PetscInt faceID, PetscInt edgeID, const PetscScalar mcoords[], PetscScalar gcoords[]) { DM cdm; ego *bodies; ego geom, body, obj; /* result has to hold derivatives, along with the value */ double params[3], result[18], paramsV[16 * 3], resultV[16 * 3], range[4]; int Nb, oclass, mtype, *senses, peri; Vec coordinatesLocal; PetscScalar *coords = NULL; PetscInt Nv, v, Np = 0, pm; PetscInt dE, d; PetscFunctionBeginHot; PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMGetCoordinateDim(dm, &dE)); PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal)); PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb); body = bodies[bodyID]; if (edgeID >= 0) { PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &obj)); Np = 1; } else if (faceID >= 0) { PetscCall(EG_objectBodyTopo(body, FACE, faceID, &obj)); Np = 2; } else { for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; PetscFunctionReturn(PETSC_SUCCESS); } /* Calculate parameters (t or u,v) for vertices */ PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); Nv /= dE; if (Nv == 1) { PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; PetscFunctionReturn(PETSC_SUCCESS); } PetscCheck(Nv <= 16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %" PetscInt_FMT " coordinates associated to point %" PetscInt_FMT, Nv, p); /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */ PetscCall(EG_getRange(obj, range, &peri)); for (v = 0; v < Nv; ++v) { PetscCall(EG_invEvaluate(obj, &coords[v * dE], ¶msV[v * 3], &resultV[v * 3])); #if 1 if (peri > 0) { if (paramsV[v * 3 + 0] + 1.e-4 < range[0]) { paramsV[v * 3 + 0] += 2. * PETSC_PI; } else if (paramsV[v * 3 + 0] - 1.e-4 > range[1]) { paramsV[v * 3 + 0] -= 2. * PETSC_PI; } } if (peri > 1) { if (paramsV[v * 3 + 1] + 1.e-4 < range[2]) { paramsV[v * 3 + 1] += 2. * PETSC_PI; } else if (paramsV[v * 3 + 1] - 1.e-4 > range[3]) { paramsV[v * 3 + 1] -= 2. * PETSC_PI; } } #endif } PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords)); /* Calculate parameters (t or u,v) for new vertex at edge midpoint */ for (pm = 0; pm < Np; ++pm) { params[pm] = 0.; for (v = 0; v < Nv; ++v) params[pm] += paramsV[v * 3 + pm]; params[pm] /= Nv; } PetscCheck(!(params[0] < range[0]) && !(params[0] > range[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p); PetscCheck(Np <= 1 || (params[1] >= range[2] && params[1] <= range[3]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p); /* Put coordinates for new vertex in result[] */ PetscCall(EG_evaluate(obj, params, result)); for (d = 0; d < dE; ++d) gcoords[d] = result[d]; PetscFunctionReturn(PETSC_SUCCESS); } #endif PetscErrorCode DMSnapToGeomModel_EGADSLite(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[]) { PetscFunctionBeginHot; #ifdef PETSC_HAVE_EGADS DMLabel bodyLabel, faceLabel, edgeLabel; PetscInt bodyID, faceID, edgeID; PetscContainer modelObj; ego model; PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite meshes must have body, face, and edge labels defined"); PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj)); PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite mesh missing model object"); PetscCall(PetscContainerGetPointer(modelObj, (void **)&model)); PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID)); PetscCall(DMLabelGetValue(faceLabel, p, &faceID)); PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID)); /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */ if (bodyID < 0) { for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(DMSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords)); #endif PetscFunctionReturn(PETSC_SUCCESS); } PetscErrorCode DMSnapToGeomModel_EGADS(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[]) { PetscFunctionBeginHot; #ifdef PETSC_HAVE_EGADS DMLabel bodyLabel, faceLabel, edgeLabel; PetscInt bodyID, faceID, edgeID; PetscContainer modelObj; ego model; PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS meshes must have body, face, and edge labels defined"); PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj)); PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS mesh missing model object"); PetscCall(PetscContainerGetPointer(modelObj, (void **)&model)); PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID)); PetscCall(DMLabelGetValue(faceLabel, p, &faceID)); PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID)); /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */ if (bodyID < 0) { for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d]; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(DMSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords)); #endif PetscFunctionReturn(PETSC_SUCCESS); } #if defined(PETSC_HAVE_EGADS) static PetscErrorCode DMPlexEGADSPrintModel_Internal(ego model) { ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; int oclass, mtype, *senses; int Nb, b; PetscFunctionBeginUser; /* test bodyTopo functions */ PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb)); for (b = 0; b < Nb; ++b) { ego body = bodies[b]; int id, Nsh, Nf, Nl, l, Ne, e, Nv, v; /* Output Basic Model Topology */ PetscCall(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh)); EG_free(objs); PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf)); EG_free(objs); PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl)); PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne)); EG_free(objs); PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv)); EG_free(objs); for (l = 0; l < Nl; ++l) { ego loop = lobjs[l]; id = EG_indexBodyTopo(body, loop); PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id)); /* Get EDGE info which associated with the current LOOP */ PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); for (e = 0; e < Ne; ++e) { ego edge = objs[e]; double range[4] = {0., 0., 0., 0.}; double point[3] = {0., 0., 0.}; double params[3] = {0., 0., 0.}; double result[18]; int peri; id = EG_indexBodyTopo(body, edge); PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e)); PetscCall(EG_getRange(edge, range, &peri)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3])); /* Get NODE info which associated with the current EDGE */ PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (mtype == DEGENERATE) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id)); } else { params[0] = range[0]; PetscCall(EG_evaluate(edge, params, result)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2])); params[0] = range[1]; PetscCall(EG_evaluate(edge, params, result)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2])); } for (v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; double limits[4]; int dummy; PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); id = EG_indexBodyTopo(body, vertex); PetscCall(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id)); PetscCall(PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2])); point[0] = point[0] + limits[0]; point[1] = point[1] + limits[1]; point[2] = point[2] + limits[2]; } } } EG_free(lobjs); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexEGADSDestroy_Private(void *context) { if (context) EG_close((ego)context); return 0; } static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm) { DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; PetscInt cStart, cEnd, c; /* EGADSLite variables */ ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; int oclass, mtype, nbodies, *senses; int b; /* PETSc variables */ DM dm; PetscHMapI edgeMap = NULL; 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; PetscInt *cells = NULL, *cone = NULL; PetscReal *coords = NULL; PetscMPIInt rank; PetscFunctionBegin; PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (rank == 0) { const PetscInt debug = 0; /* --------------------------------------------------------------------------------------------------- Generate Petsc Plex Get all Nodes in model, record coordinates in a correctly formatted array Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array We need to uniformly refine the initial geometry to guarantee a valid mesh */ /* Calculate cell and vertex sizes */ PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); PetscCall(PetscHMapICreate(&edgeMap)); numEdges = 0; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int id, Nl, l, Nv, v; PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); for (l = 0; l < Nl; ++l) { ego loop = lobjs[l]; int Ner = 0, Ne, e, Nc; PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); for (e = 0; e < Ne; ++e) { ego edge = objs[e]; int Nv, id; PetscHashIter iter; PetscBool found; PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (mtype == DEGENERATE) continue; id = EG_indexBodyTopo(body, edge); PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found)); if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++)); ++Ner; } if (Ner == 2) { Nc = 2; } else if (Ner == 3) { Nc = 4; } else if (Ner == 4) { Nc = 8; ++numQuads; } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner); numCells += Nc; newCells += Nc - 1; maxCorners = PetscMax(Ner * 2 + 1, maxCorners); } PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); for (v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; id = EG_indexBodyTopo(body, vertex); /* TODO: Instead of assuming contiguous ids, we could use a hash table */ numVertices = PetscMax(id, numVertices); } EG_free(lobjs); EG_free(nobjs); } PetscCall(PetscHMapIGetSize(edgeMap, &numEdges)); newVertices = numEdges + numQuads; numVertices += newVertices; dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ numCorners = 3; /* Split cells into triangles */ PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone)); /* Get vertex coordinates */ for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int id, Nv, v; PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); for (v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; double limits[4]; int dummy; PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); id = EG_indexBodyTopo(body, vertex); coords[(id - 1) * cdim + 0] = limits[0]; coords[(id - 1) * cdim + 1] = limits[1]; coords[(id - 1) * cdim + 2] = limits[2]; } EG_free(nobjs); } PetscCall(PetscHMapIClear(edgeMap)); fOff = numVertices - newVertices + numEdges; numEdges = 0; numQuads = 0; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int Nl, l; PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); for (l = 0; l < Nl; ++l) { ego loop = lobjs[l]; int lid, Ner = 0, Ne, e; lid = EG_indexBodyTopo(body, loop); PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); for (e = 0; e < Ne; ++e) { ego edge = objs[e]; int eid, Nv; PetscHashIter iter; PetscBool found; PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (mtype == DEGENERATE) continue; ++Ner; eid = EG_indexBodyTopo(body, edge); PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found)); if (!found) { PetscInt v = numVertices - newVertices + numEdges; double range[4], params[3] = {0., 0., 0.}, result[18]; int periodic[2]; PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++)); PetscCall(EG_getRange(edge, range, periodic)); params[0] = 0.5 * (range[0] + range[1]); PetscCall(EG_evaluate(edge, params, result)); coords[v * cdim + 0] = result[0]; coords[v * cdim + 1] = result[1]; coords[v * cdim + 2] = result[2]; } } if (Ner == 4) { PetscInt v = fOff + numQuads++; ego *fobjs, face; double range[4], params[3] = {0., 0., 0.}, result[18]; int Nf, fid, periodic[2]; PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); face = fobjs[0]; fid = EG_indexBodyTopo(body, face); PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid); PetscCall(EG_getRange(face, range, periodic)); params[0] = 0.5 * (range[0] + range[1]); params[1] = 0.5 * (range[2] + range[3]); PetscCall(EG_evaluate(face, params, result)); coords[v * cdim + 0] = result[0]; coords[v * cdim + 1] = result[1]; coords[v * cdim + 2] = result[2]; } } } PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices); /* Get cell vertices by traversing loops */ numQuads = 0; cOff = 0; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int id, Nl, l; PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); for (l = 0; l < Nl; ++l) { ego loop = lobjs[l]; int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t; lid = EG_indexBodyTopo(body, loop); PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); for (e = 0; e < Ne; ++e) { ego edge = objs[e]; int points[3]; int eid, Nv, v, tmp; eid = EG_indexBodyTopo(body, edge); PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (mtype == DEGENERATE) continue; else ++Ner; PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv); for (v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; id = EG_indexBodyTopo(body, vertex); points[v * 2] = id - 1; } { PetscInt edgeNum; PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum)); points[1] = numVertices - newVertices + edgeNum; } /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */ if (!nc) { for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v]; } else { if (cone[nc - 1] == points[0]) { cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2]; } else if (cone[nc - 1] == points[2]) { cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0]; } 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]; } 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]; } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); } } PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner); if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++; PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners); /* Triangulate the loop */ switch (Ner) { case 2: /* Bi-Segment -> 2 triangles */ Nt = 2; cells[cOff * numCorners + 0] = cone[0]; cells[cOff * numCorners + 1] = cone[1]; cells[cOff * numCorners + 2] = cone[2]; ++cOff; cells[cOff * numCorners + 0] = cone[0]; cells[cOff * numCorners + 1] = cone[2]; cells[cOff * numCorners + 2] = cone[3]; ++cOff; break; case 3: /* Triangle -> 4 triangles */ Nt = 4; cells[cOff * numCorners + 0] = cone[0]; cells[cOff * numCorners + 1] = cone[1]; cells[cOff * numCorners + 2] = cone[5]; ++cOff; cells[cOff * numCorners + 0] = cone[1]; cells[cOff * numCorners + 1] = cone[2]; cells[cOff * numCorners + 2] = cone[3]; ++cOff; cells[cOff * numCorners + 0] = cone[5]; cells[cOff * numCorners + 1] = cone[3]; cells[cOff * numCorners + 2] = cone[4]; ++cOff; cells[cOff * numCorners + 0] = cone[1]; cells[cOff * numCorners + 1] = cone[3]; cells[cOff * numCorners + 2] = cone[5]; ++cOff; break; case 4: /* Quad -> 8 triangles */ Nt = 8; cells[cOff * numCorners + 0] = cone[0]; cells[cOff * numCorners + 1] = cone[1]; cells[cOff * numCorners + 2] = cone[7]; ++cOff; cells[cOff * numCorners + 0] = cone[1]; cells[cOff * numCorners + 1] = cone[2]; cells[cOff * numCorners + 2] = cone[3]; ++cOff; cells[cOff * numCorners + 0] = cone[3]; cells[cOff * numCorners + 1] = cone[4]; cells[cOff * numCorners + 2] = cone[5]; ++cOff; cells[cOff * numCorners + 0] = cone[5]; cells[cOff * numCorners + 1] = cone[6]; cells[cOff * numCorners + 2] = cone[7]; ++cOff; cells[cOff * numCorners + 0] = cone[8]; cells[cOff * numCorners + 1] = cone[1]; cells[cOff * numCorners + 2] = cone[3]; ++cOff; cells[cOff * numCorners + 0] = cone[8]; cells[cOff * numCorners + 1] = cone[3]; cells[cOff * numCorners + 2] = cone[5]; ++cOff; cells[cOff * numCorners + 0] = cone[8]; cells[cOff * numCorners + 1] = cone[5]; cells[cOff * numCorners + 2] = cone[7]; ++cOff; cells[cOff * numCorners + 0] = cone[8]; cells[cOff * numCorners + 1] = cone[7]; cells[cOff * numCorners + 2] = cone[1]; ++cOff; break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner); } if (debug) { for (t = 0; t < Nt; ++t) { PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t)); for (c = 0; c < numCorners; ++c) { if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", ")); PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c])); } PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n")); } } } EG_free(lobjs); } } PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells); PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); PetscCall(PetscFree3(coords, cells, cone)); PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells)); PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices)); /* Embed EGADS model in DM */ { PetscContainer modelObj, contextObj; PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); PetscCall(PetscContainerSetPointer(modelObj, model)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); PetscCall(PetscContainerDestroy(&modelObj)); PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); PetscCall(PetscContainerSetPointer(contextObj, context)); PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); PetscCall(PetscContainerDestroy(&contextObj)); } /* Label points */ PetscCall(DMCreateLabel(dm, "EGADS Body ID")); PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMCreateLabel(dm, "EGADS Face ID")); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); cOff = 0; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int id, Nl, l; PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs)); for (l = 0; l < Nl; ++l) { ego loop = lobjs[l]; ego *fobjs; int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t; lid = EG_indexBodyTopo(body, loop); PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs)); PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf); fid = EG_indexBodyTopo(body, fobjs[0]); EG_free(fobjs); PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses)); for (e = 0; e < Ne; ++e) { ego edge = objs[e]; int eid, Nv, v; PetscInt points[3], support[2], numEdges, edgeNum; const PetscInt *edges; eid = EG_indexBodyTopo(body, edge); PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (mtype == DEGENERATE) continue; else ++Ner; for (v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; id = EG_indexBodyTopo(body, vertex); PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid)); points[v * 2] = numCells + id - 1; } PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum)); points[1] = numCells + numVertices - newVertices + edgeNum; PetscCall(DMLabelSetValue(edgeLabel, points[1], eid)); support[0] = points[0]; support[1] = points[1]; PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges); PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid)); PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges)); support[0] = points[1]; support[1] = points[2]; PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges)); PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges); PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid)); PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges)); } switch (Ner) { case 2: Nt = 2; break; case 3: Nt = 4; break; case 4: Nt = 8; break; default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner); } for (t = 0; t < Nt; ++t) { PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b)); PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid)); } cOff += Nt; } EG_free(lobjs); } PetscCall(PetscHMapIDestroy(&edgeMap)); PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); for (c = cStart; c < cEnd; ++c) { PetscInt *closure = NULL; PetscInt clSize, cl, bval, fval; PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); PetscCall(DMLabelGetValue(bodyLabel, c, &bval)); PetscCall(DMLabelGetValue(faceLabel, c, &fval)); for (cl = 0; cl < clSize * 2; cl += 2) { PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval)); PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval)); } PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure)); } *newdm = dm; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm) { DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; // EGADS/EGADSLite variables ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs; ego topRef, prev, next; int oclass, mtype, nbodies, *senses, *lSenses, *eSenses; int b; // PETSc variables DM dm; PetscHMapI edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL; PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0; PetscInt cellCntr = 0, numPoints = 0; PetscInt *cells = NULL; const PetscInt *cone = NULL; PetscReal *coords = NULL; PetscMPIInt rank; PetscFunctionBeginUser; PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (rank == 0) { // --------------------------------------------------------------------------------------------------- // Generate Petsc Plex // Get all Nodes in model, record coordinates in a correctly formatted array // Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array // We need to uniformly refine the initial geometry to guarantee a valid mesh // Calculate cell and vertex sizes PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); PetscCall(PetscHMapICreate(&edgeMap)); PetscCall(PetscHMapICreate(&bodyIndexMap)); PetscCall(PetscHMapICreate(&bodyVertexMap)); PetscCall(PetscHMapICreate(&bodyEdgeMap)); PetscCall(PetscHMapICreate(&bodyEdgeGlobalMap)); PetscCall(PetscHMapICreate(&bodyFaceMap)); for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int Nf, Ne, Nv; PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter; PetscBool BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound; PetscCall(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound)); PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); if (!BIfound) PetscCall(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices)); if (!BVfound) PetscCall(PetscHMapISet(bodyVertexMap, b, numVertices)); if (!BEfound) PetscCall(PetscHMapISet(bodyEdgeMap, b, numEdges)); if (!BEGfound) PetscCall(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr)); if (!BFfound) PetscCall(PetscHMapISet(bodyFaceMap, b, numFaces)); PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); EG_free(fobjs); EG_free(eobjs); EG_free(nobjs); // Remove DEGENERATE EDGES from Edge count PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); int Netemp = 0; for (int e = 0; e < Ne; ++e) { ego edge = eobjs[e]; int eid; PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); eid = EG_indexBodyTopo(body, edge); PetscCall(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound)); if (mtype == DEGENERATE) { if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1)); } else { ++Netemp; if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp)); } } EG_free(eobjs); // Determine Number of Cells PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); for (int f = 0; f < Nf; ++f) { ego face = fobjs[f]; int edgeTemp = 0; PetscCall(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs)); for (int e = 0; e < Ne; ++e) { ego edge = eobjs[e]; PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); if (mtype != DEGENERATE) ++edgeTemp; } numCells += (2 * edgeTemp); EG_free(eobjs); } EG_free(fobjs); numFaces += Nf; numEdges += Netemp; numVertices += Nv; edgeCntr += Ne; } // Set up basic DMPlex parameters dim = 2; // Assumes 3D Models :: Need to handle 2D models in the future cdim = 3; // Assumes 3D Models :: Need to update to handle 2D models in future numCorners = 3; // Split Faces into triangles numPoints = numVertices + numEdges + numFaces; // total number of coordinate points PetscCall(PetscMalloc2(numPoints * cdim, &coords, numCells * numCorners, &cells)); // Get Vertex Coordinates and Set up Cells for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int Nf, Ne, Nv; PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart; PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter; PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound; // Vertices on Current Body PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs)); PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart)); for (int v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; double limits[4]; int id, dummy; PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses)); id = EG_indexBodyTopo(body, vertex); coords[(bodyVertexIndexStart + id - 1) * cdim + 0] = limits[0]; coords[(bodyVertexIndexStart + id - 1) * cdim + 1] = limits[1]; coords[(bodyVertexIndexStart + id - 1) * cdim + 2] = limits[2]; } EG_free(nobjs); // Edge Midpoint Vertices on Current Body PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs)); PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart)); for (int e = 0; e < Ne; ++e) { ego edge = eobjs[e]; double range[2], avgt[1], cntrPnt[9]; int eid, eOffset; int periodic; PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); if (mtype == DEGENERATE) continue; eid = EG_indexBodyTopo(body, edge); // get relative offset from globalEdgeID Vector PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1); PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); PetscCall(EG_getRange(edge, range, &periodic)); avgt[0] = (range[0] + range[1]) / 2.; PetscCall(EG_evaluate(edge, avgt, cntrPnt)); coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 0] = cntrPnt[0]; coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 1] = cntrPnt[1]; coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 2] = cntrPnt[2]; } EG_free(eobjs); // Face Midpoint Vertices on Current Body PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart)); for (int f = 0; f < Nf; ++f) { ego face = fobjs[f]; double range[4], avgUV[2], cntrPnt[18]; int peri, id; id = EG_indexBodyTopo(body, face); PetscCall(EG_getRange(face, range, &peri)); avgUV[0] = (range[0] + range[1]) / 2.; avgUV[1] = (range[2] + range[3]) / 2.; PetscCall(EG_evaluate(face, avgUV, cntrPnt)); coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 0] = cntrPnt[0]; coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 1] = cntrPnt[1]; coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 2] = cntrPnt[2]; } EG_free(fobjs); // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); for (int f = 0; f < Nf; ++f) { ego face = fobjs[f]; int fID, midFaceID, midPntID, startID, endID, Nl; fID = EG_indexBodyTopo(body, face); midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1; // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges. // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are: // 1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option. // This will use a default EGADS tessellation as an initial surface mesh. // 2) Create the initial surface mesh via a 2D mesher :: Currently not available (?future?) // May I suggest the XXXX as a starting point? PetscCall(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses)); PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face has %d Loops. Can only handle Faces with 1 Loop. Please use --dm_plex_egads_with_tess = 1 Option", Nl); for (int l = 0; l < Nl; ++l) { ego loop = lobjs[l]; PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses)); for (int e = 0; e < Ne; ++e) { ego edge = eobjs[e]; int eid, eOffset; PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); eid = EG_indexBodyTopo(body, edge); if (mtype == DEGENERATE) continue; // get relative offset from globalEdgeID Vector PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1); PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1; PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses)); if (eSenses[e] > 0) { startID = EG_indexBodyTopo(body, nobjs[0]); endID = EG_indexBodyTopo(body, nobjs[1]); } else { startID = EG_indexBodyTopo(body, nobjs[1]); endID = EG_indexBodyTopo(body, nobjs[0]); } // Define 2 Cells per Edge with correct orientation cells[cellCntr * numCorners + 0] = midFaceID; cells[cellCntr * numCorners + 1] = bodyVertexIndexStart + startID - 1; cells[cellCntr * numCorners + 2] = midPntID; cells[cellCntr * numCorners + 3] = midFaceID; cells[cellCntr * numCorners + 4] = midPntID; cells[cellCntr * numCorners + 5] = bodyVertexIndexStart + endID - 1; cellCntr = cellCntr + 2; } } } EG_free(fobjs); } } // Generate DMPlex PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); PetscCall(PetscFree2(coords, cells)); PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " \n", numCells)); PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " \n", numVertices)); // Embed EGADS model in DM { PetscContainer modelObj, contextObj; PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); PetscCall(PetscContainerSetPointer(modelObj, model)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); PetscCall(PetscContainerDestroy(&modelObj)); PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); PetscCall(PetscContainerSetPointer(contextObj, context)); PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); PetscCall(PetscContainerDestroy(&contextObj)); } // Label points PetscInt nStart, nEnd; PetscCall(DMCreateLabel(dm, "EGADS Body ID")); PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMCreateLabel(dm, "EGADS Face ID")); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd)); cellCntr = 0; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int Nv, Ne, Nf; PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart; PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter; PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound; PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound)); PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b); PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart)); PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound)); PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b); PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart)); PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound)); PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b); PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart)); PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound)); PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b); PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart)); PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); for (int f = 0; f < Nf; ++f) { ego face = fobjs[f]; int fID, Nl; fID = EG_indexBodyTopo(body, face); PetscCall(EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs)); for (int l = 0; l < Nl; ++l) { ego loop = lobjs[l]; int lid; lid = EG_indexBodyTopo(body, loop); PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf); PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses)); for (int e = 0; e < Ne; ++e) { ego edge = eobjs[e]; int eid, eOffset; // Skip DEGENERATE Edges PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next)); if (mtype == DEGENERATE) continue; eid = EG_indexBodyTopo(body, edge); // get relative offset from globalEdgeID Vector PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound)); PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1); PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset)); PetscCall(EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs)); for (int v = 0; v < Nv; ++v) { ego vertex = nobjs[v]; int vID; vID = EG_indexBodyTopo(body, vertex); PetscCall(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b)); PetscCall(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID)); } EG_free(nobjs); PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b)); PetscCall(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid)); // Define Cell faces for (int jj = 0; jj < 2; ++jj) { PetscCall(DMLabelSetValue(bodyLabel, cellCntr, b)); PetscCall(DMLabelSetValue(faceLabel, cellCntr, fID)); PetscCall(DMPlexGetCone(dm, cellCntr, &cone)); PetscCall(DMLabelSetValue(bodyLabel, cone[0], b)); PetscCall(DMLabelSetValue(faceLabel, cone[0], fID)); PetscCall(DMLabelSetValue(bodyLabel, cone[1], b)); PetscCall(DMLabelSetValue(edgeLabel, cone[1], eid)); PetscCall(DMLabelSetValue(bodyLabel, cone[2], b)); PetscCall(DMLabelSetValue(faceLabel, cone[2], fID)); cellCntr = cellCntr + 1; } } } EG_free(lobjs); PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b)); PetscCall(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID)); } EG_free(fobjs); } PetscCall(PetscHMapIDestroy(&edgeMap)); PetscCall(PetscHMapIDestroy(&bodyIndexMap)); PetscCall(PetscHMapIDestroy(&bodyVertexMap)); PetscCall(PetscHMapIDestroy(&bodyEdgeMap)); PetscCall(PetscHMapIDestroy(&bodyEdgeGlobalMap)); PetscCall(PetscHMapIDestroy(&bodyFaceMap)); *newdm = dm; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm) { DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; /* EGADSLite variables */ ego geom, *bodies, *fobjs; int b, oclass, mtype, nbodies, *senses; int totalNumTris = 0, totalNumPoints = 0; double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize; /* PETSc variables */ DM dm; PetscHMapI pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL; PetscHMapI pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL; PetscInt dim = -1, cdim = -1, numCorners = 0, counter = 0; PetscInt *cells = NULL; const PetscInt *cone = NULL; PetscReal *coords = NULL; PetscMPIInt rank; PetscFunctionBeginUser; PetscCallMPI(MPI_Comm_rank(comm, &rank)); if (rank == 0) { // --------------------------------------------------------------------------------------------------- // Generate Petsc Plex from EGADSlite created Tessellation of geometry // --------------------------------------------------------------------------------------------------- // Calculate cell and vertex sizes PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses)); PetscCall(PetscHMapICreate(&pointIndexStartMap)); PetscCall(PetscHMapICreate(&triIndexStartMap)); PetscCall(PetscHMapICreate(&pTypeLabelMap)); PetscCall(PetscHMapICreate(&pIndexLabelMap)); PetscCall(PetscHMapICreate(&pBodyIndexLabelMap)); PetscCall(PetscHMapICreate(&triFaceIDLabelMap)); PetscCall(PetscHMapICreate(&triBodyIDLabelMap)); /* Create Tessellation of Bodies */ ego tessArray[nbodies]; for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; double params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation int Nf, bodyNumPoints = 0, bodyNumTris = 0; PetscHashIter PISiter, TISiter; PetscBool PISfound, TISfound; /* Store Start Index for each Body's Point and Tris */ PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); PetscCall(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound)); if (!PISfound) PetscCall(PetscHMapISet(pointIndexStartMap, b, totalNumPoints)); if (!TISfound) PetscCall(PetscHMapISet(triIndexStartMap, b, totalNumTris)); /* Calculate Tessellation parameters based on Bounding Box */ /* Get Bounding Box Dimensions of the BODY */ PetscCall(EG_getBoundingBox(body, boundBox)); tessSize = boundBox[3] - boundBox[0]; if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1]; if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2]; // TODO :: May want to give users tessellation parameter options // params[0] = 0.0250 * tessSize; params[1] = 0.0075 * tessSize; params[2] = 15.0; PetscCall(EG_makeTessBody(body, params, &tessArray[b])); PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); for (int f = 0; f < Nf; ++f) { ego face = fobjs[f]; int len, fID, ntris; const int *ptype, *pindex, *ptris, *ptric; const double *pxyz, *puv; // Get Face ID // fID = EG_indexBodyTopo(body, face); // Checkout the Surface Tessellation // PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric)); // Determine total number of triangle cells in the tessellation // bodyNumTris += (int)ntris; // Check out the point index and coordinate // for (int p = 0; p < len; ++p) { int global; PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global)); // Determine the total number of points in the tessellation // bodyNumPoints = PetscMax(bodyNumPoints, global); } } EG_free(fobjs); totalNumPoints += bodyNumPoints; totalNumTris += bodyNumTris; } //} - Original End of (rank == 0) dim = 2; cdim = 3; numCorners = 3; //PetscInt counter = 0; /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA */ /* Fill in below and use to define DMLabels after DMPlex creation */ PetscCall(PetscMalloc2(totalNumPoints * cdim, &coords, totalNumTris * numCorners, &cells)); for (b = 0; b < nbodies; ++b) { ego body = bodies[b]; int Nf; PetscInt pointIndexStart; PetscHashIter PISiter; PetscBool PISfound; PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound)); PetscCheck(PISfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b); PetscCall(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart)); PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs)); for (int f = 0; f < Nf; ++f) { /* Get Face Object */ ego face = fobjs[f]; int len, fID, ntris; const int *ptype, *pindex, *ptris, *ptric; const double *pxyz, *puv; /* Get Face ID */ fID = EG_indexBodyTopo(body, face); /* Checkout the Surface Tessellation */ PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric)); /* Check out the point index and coordinate */ for (int p = 0; p < len; ++p) { int global; PetscHashIter PTLiter, PILiter, PBLiter; PetscBool PTLfound, PILfound, PBLfound; PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global)); /* Set the coordinates array for DAG */ coords[((global - 1 + pointIndexStart) * 3) + 0] = pxyz[(p * 3) + 0]; coords[((global - 1 + pointIndexStart) * 3) + 1] = pxyz[(p * 3) + 1]; coords[((global - 1 + pointIndexStart) * 3) + 2] = pxyz[(p * 3) + 2]; /* Store Geometry Label Information for DMLabel assignment later */ PetscCall(PetscHMapIFind(pTypeLabelMap, global - 1 + pointIndexStart, &PTLiter, &PTLfound)); PetscCall(PetscHMapIFind(pIndexLabelMap, global - 1 + pointIndexStart, &PILiter, &PILfound)); PetscCall(PetscHMapIFind(pBodyIndexLabelMap, global - 1 + pointIndexStart, &PBLiter, &PBLfound)); if (!PTLfound) PetscCall(PetscHMapISet(pTypeLabelMap, global - 1 + pointIndexStart, ptype[p])); if (!PILfound) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, pindex[p])); if (!PBLfound) PetscCall(PetscHMapISet(pBodyIndexLabelMap, global - 1 + pointIndexStart, b)); if (ptype[p] < 0) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, fID)); } for (int t = 0; t < (int)ntris; ++t) { int global, globalA, globalB; PetscHashIter TFLiter, TBLiter; PetscBool TFLfound, TBLfound; PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 0], &global)); cells[(counter * 3) + 0] = global - 1 + pointIndexStart; PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 1], &globalA)); cells[(counter * 3) + 1] = globalA - 1 + pointIndexStart; PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 2], &globalB)); cells[(counter * 3) + 2] = globalB - 1 + pointIndexStart; PetscCall(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound)); PetscCall(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound)); if (!TFLfound) PetscCall(PetscHMapISet(triFaceIDLabelMap, counter, fID)); if (!TBLfound) PetscCall(PetscHMapISet(triBodyIDLabelMap, counter, b)); counter += 1; } } EG_free(fobjs); } } //Build DMPlex PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm)); PetscCall(PetscFree2(coords, cells)); // Embed EGADS model in DM { PetscContainer modelObj, contextObj; PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj)); PetscCall(PetscContainerSetPointer(modelObj, model)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj)); PetscCall(PetscContainerDestroy(&modelObj)); PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj)); PetscCall(PetscContainerSetPointer(contextObj, context)); PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private)); PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj)); PetscCall(PetscContainerDestroy(&contextObj)); } // Label Points PetscCall(DMCreateLabel(dm, "EGADS Body ID")); PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMCreateLabel(dm, "EGADS Face ID")); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMCreateLabel(dm, "EGADS Edge ID")); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCall(DMCreateLabel(dm, "EGADS Vertex ID")); PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); /* Get Number of DAG Nodes at each level */ int fStart, fEnd, eStart, eEnd, nStart, nEnd; PetscCall(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd)); PetscCall(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd)); PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd)); /* Set DMLabels for NODES */ for (int n = nStart; n < nEnd; ++n) { int pTypeVal, pIndexVal, pBodyVal; PetscHashIter PTLiter, PILiter, PBLiter; PetscBool PTLfound, PILfound, PBLfound; //Converted to Hash Tables PetscCall(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound)); PetscCheck(PTLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n); PetscCall(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal)); PetscCall(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound)); PetscCheck(PILfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n); PetscCall(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal)); PetscCall(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound)); PetscCheck(PBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n); PetscCall(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal)); PetscCall(DMLabelSetValue(bodyLabel, n, pBodyVal)); if (pTypeVal == 0) PetscCall(DMLabelSetValue(vertexLabel, n, pIndexVal)); if (pTypeVal > 0) PetscCall(DMLabelSetValue(edgeLabel, n, pIndexVal)); if (pTypeVal < 0) PetscCall(DMLabelSetValue(faceLabel, n, pIndexVal)); } /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */ for (int e = eStart; e < eEnd; ++e) { int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1; PetscCall(DMPlexGetCone(dm, e, &cone)); PetscCall(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end? PetscCall(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0)); PetscCall(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1)); PetscCall(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0)); PetscCall(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1)); PetscCall(DMLabelGetValue(faceLabel, cone[0], &faceID_0)); PetscCall(DMLabelGetValue(faceLabel, cone[1], &faceID_1)); PetscCall(DMLabelSetValue(bodyLabel, e, bodyID_0)); if (edgeID_0 == edgeID_1) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0)); else if (vertexID_0 > 0 && edgeID_1 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_1)); else if (vertexID_1 > 0 && edgeID_0 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0)); else { /* Do Nothing */ } } /* Set DMLabels for Cells */ for (int f = fStart; f < fEnd; ++f) { int edgeID_0; PetscInt triBodyVal, triFaceVal; PetscHashIter TFLiter, TBLiter; PetscBool TFLfound, TBLfound; // Convert to Hash Table PetscCall(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound)); PetscCheck(TFLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f); PetscCall(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal)); PetscCall(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound)); PetscCheck(TBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f); PetscCall(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal)); PetscCall(DMLabelSetValue(bodyLabel, f, triBodyVal)); PetscCall(DMLabelSetValue(faceLabel, f, triFaceVal)); /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */ PetscCall(DMPlexGetCone(dm, f, &cone)); for (int jj = 0; jj < 3; ++jj) { PetscCall(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0)); if (edgeID_0 < 0) { PetscCall(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal)); PetscCall(DMLabelSetValue(faceLabel, cone[jj], triFaceVal)); } } } *newdm = dm; PetscFunctionReturn(PETSC_SUCCESS); } #endif /*@ 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. Collective Input Parameter: . dm - The uninflated `DM` object representing the mesh Level: intermediate .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()` @*/ PetscErrorCode DMPlexInflateToGeomModel(DM dm) { #if defined(PETSC_HAVE_EGADS) /* EGADS Variables */ ego model, geom, body, face, edge; ego *bodies; int Nb, oclass, mtype, *senses; double result[3]; /* PETSc Variables */ DM cdm; PetscContainer modelObj; DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel; Vec coordinates; PetscScalar *coords; PetscInt bodyID, faceID, edgeID, vertexID; PetscInt cdim, d, vStart, vEnd, v; #endif PetscFunctionBegin; #if defined(PETSC_HAVE_EGADS) PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj)); if (!modelObj) PetscFunctionReturn(PETSC_SUCCESS); PetscCall(DMGetCoordinateDim(dm, &cdim)); PetscCall(DMGetCoordinateDM(dm, &cdm)); PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel)); PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel)); PetscCall(PetscContainerGetPointer(modelObj, (void **)&model)); PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses)); PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); PetscCall(VecGetArrayWrite(coordinates, &coords)); for (v = vStart; v < vEnd; ++v) { PetscScalar *vcoords; PetscCall(DMLabelGetValue(bodyLabel, v, &bodyID)); PetscCall(DMLabelGetValue(faceLabel, v, &faceID)); PetscCall(DMLabelGetValue(edgeLabel, v, &edgeID)); PetscCall(DMLabelGetValue(vertexLabel, v, &vertexID)); PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb); body = bodies[bodyID]; PetscCall(DMPlexPointLocalRef(cdm, v, coords, (void *)&vcoords)); if (edgeID > 0) { /* Snap to EDGE at nearest location */ double params[1]; PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &edge)); PetscCall(EG_invEvaluate(edge, vcoords, params, result)); // Get (x,y,z) of nearest point on EDGE for (d = 0; d < cdim; ++d) vcoords[d] = result[d]; } else if (faceID > 0) { /* Snap to FACE at nearest location */ double params[2]; PetscCall(EG_objectBodyTopo(body, FACE, faceID, &face)); PetscCall(EG_invEvaluate(face, vcoords, params, result)); // Get (x,y,z) of nearest point on FACE for (d = 0; d < cdim; ++d) vcoords[d] = result[d]; } } PetscCall(VecRestoreArrayWrite(coordinates, &coords)); /* Clear out global coordinates */ PetscCall(VecDestroy(&dm->coordinates[0].x)); #endif PetscFunctionReturn(PETSC_SUCCESS); } /*@ DMPlexCreateEGADSFromFile - Create a `DMPLEX` mesh from an EGADS, IGES, or STEP file. Collective Input Parameters: + comm - The MPI communicator - filename - The name of the EGADS, IGES, or STEP file Output Parameter: . dm - The `DM` object representing the mesh Level: beginner .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSLiteFromFile()` @*/ PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm) { PetscMPIInt rank; #if defined(PETSC_HAVE_EGADS) ego context = NULL, model = NULL; #endif PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE; PetscFunctionBegin; PetscAssertPointer(filename, 2); PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL)); PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL)); PetscCallMPI(MPI_Comm_rank(comm, &rank)); #if defined(PETSC_HAVE_EGADS) if (rank == 0) { PetscCall(EG_open(&context)); PetscCall(EG_loadModel(context, 0, filename, &model)); if (printModel) PetscCall(DMPlexEGADSPrintModel_Internal(model)); } if (tessModel) PetscCall(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm)); else if (newModel) PetscCall(DMPlexCreateEGADS(comm, context, model, dm)); else PetscCall(DMPlexCreateEGADS_Internal(comm, context, model, dm)); PetscFunctionReturn(PETSC_SUCCESS); #else SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads"); #endif }