Actual source code: plexegads.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/hashmapi.h>

      4: #ifdef PETSC_HAVE_EGADS
      5:   #include <egads.h>
      6: #endif

      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:

     11:      https://github.com/tpaviot/oce/tree/master/src/STEPControl

     13:    The STEP, and inner EXPRESS, formats are ISO standards, so they are documented

     15:      https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files
     16:      http://stepmod.sourceforge.net/express_model_spec/

     18:    but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants.
     19: */

     21: #ifdef PETSC_HAVE_EGADS
     22: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
     23: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);

     25: PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM dm, PetscInt p, ego model, PetscInt bodyID, PetscInt faceID, PetscInt edgeID, const PetscScalar mcoords[], PetscScalar gcoords[])
     26: {
     27:   DM   cdm;
     28:   ego *bodies;
     29:   ego  geom, body, obj;
     30:   /* result has to hold derivatives, along with the value */
     31:   double       params[3], result[18], paramsV[16 * 3], resultV[16 * 3], range[4];
     32:   int          Nb, oclass, mtype, *senses, peri;
     33:   Vec          coordinatesLocal;
     34:   PetscScalar *coords = NULL;
     35:   PetscInt     Nv, v, Np = 0, pm;
     36:   PetscInt     dE, d;

     38:   PetscFunctionBeginHot;
     39:   PetscCall(DMGetCoordinateDM(dm, &cdm));
     40:   PetscCall(DMGetCoordinateDim(dm, &dE));
     41:   PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
     42:   PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
     43:   PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
     44:   body = bodies[bodyID];

     46:   if (edgeID >= 0) {
     47:     PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &obj));
     48:     Np = 1;
     49:   } else if (faceID >= 0) {
     50:     PetscCall(EG_objectBodyTopo(body, FACE, faceID, &obj));
     51:     Np = 2;
     52:   } else {
     53:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
     54:     PetscFunctionReturn(PETSC_SUCCESS);
     55:   }

     57:   /* Calculate parameters (t or u,v) for vertices */
     58:   PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
     59:   Nv /= dE;
     60:   if (Nv == 1) {
     61:     PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
     62:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
     63:     PetscFunctionReturn(PETSC_SUCCESS);
     64:   }
     65:   PetscCheck(Nv <= 16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %" PetscInt_FMT " coordinates associated to point %" PetscInt_FMT, Nv, p);

     67:   /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */
     68:   PetscCall(EG_getRange(obj, range, &peri));
     69:   for (v = 0; v < Nv; ++v) {
     70:     PetscCall(EG_invEvaluate(obj, &coords[v * dE], &paramsV[v * 3], &resultV[v * 3]));
     71:   #if 1
     72:     if (peri > 0) {
     73:       if (paramsV[v * 3 + 0] + 1.e-4 < range[0]) {
     74:         paramsV[v * 3 + 0] += 2. * PETSC_PI;
     75:       } else if (paramsV[v * 3 + 0] - 1.e-4 > range[1]) {
     76:         paramsV[v * 3 + 0] -= 2. * PETSC_PI;
     77:       }
     78:     }
     79:     if (peri > 1) {
     80:       if (paramsV[v * 3 + 1] + 1.e-4 < range[2]) {
     81:         paramsV[v * 3 + 1] += 2. * PETSC_PI;
     82:       } else if (paramsV[v * 3 + 1] - 1.e-4 > range[3]) {
     83:         paramsV[v * 3 + 1] -= 2. * PETSC_PI;
     84:       }
     85:     }
     86:   #endif
     87:   }
     88:   PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
     89:   /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
     90:   for (pm = 0; pm < Np; ++pm) {
     91:     params[pm] = 0.;
     92:     for (v = 0; v < Nv; ++v) params[pm] += paramsV[v * 3 + pm];
     93:     params[pm] /= Nv;
     94:   }
     95:   PetscCheck(!(params[0] < range[0]) && !(params[0] > range[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p);
     96:   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);
     97:   /* Put coordinates for new vertex in result[] */
     98:   PetscCall(EG_evaluate(obj, params, result));
     99:   for (d = 0; d < dE; ++d) gcoords[d] = result[d];
    100:   PetscFunctionReturn(PETSC_SUCCESS);
    101: }
    102: #endif

    104: PetscErrorCode DMSnapToGeomModel_EGADSLite(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
    105: {
    106:   PetscFunctionBeginHot;
    107: #ifdef PETSC_HAVE_EGADS
    108:   DMLabel        bodyLabel, faceLabel, edgeLabel;
    109:   PetscInt       bodyID, faceID, edgeID;
    110:   PetscContainer modelObj;
    111:   ego            model;

    113:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    114:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    115:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    116:   PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite meshes must have body, face, and edge labels defined");
    117:   PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj));
    118:   PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite mesh missing model object");

    120:   PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
    121:   PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID));
    122:   PetscCall(DMLabelGetValue(faceLabel, p, &faceID));
    123:   PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID));
    124:   /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
    125:   if (bodyID < 0) {
    126:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
    127:     PetscFunctionReturn(PETSC_SUCCESS);
    128:   }
    129:   PetscCall(DMSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
    130: #endif
    131:   PetscFunctionReturn(PETSC_SUCCESS);
    132: }

    134: PetscErrorCode DMSnapToGeomModel_EGADS(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
    135: {
    136:   PetscFunctionBeginHot;
    137: #ifdef PETSC_HAVE_EGADS
    138:   DMLabel        bodyLabel, faceLabel, edgeLabel;
    139:   PetscInt       bodyID, faceID, edgeID;
    140:   PetscContainer modelObj;
    141:   ego            model;

    143:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    144:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    145:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    146:   PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS meshes must have body, face, and edge labels defined");
    147:   PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
    148:   PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS mesh missing model object");

    150:   PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
    151:   PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID));
    152:   PetscCall(DMLabelGetValue(faceLabel, p, &faceID));
    153:   PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID));
    154:   /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
    155:   if (bodyID < 0) {
    156:     for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
    157:     PetscFunctionReturn(PETSC_SUCCESS);
    158:   }
    159:   PetscCall(DMSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
    160: #endif
    161:   PetscFunctionReturn(PETSC_SUCCESS);
    162: }

    164: #if defined(PETSC_HAVE_EGADS)
    165: static PetscErrorCode DMPlexEGADSPrintModel_Internal(ego model)
    166: {
    167:   ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
    168:   int oclass, mtype, *senses;
    169:   int Nb, b;

    171:   PetscFunctionBeginUser;
    172:   /* test bodyTopo functions */
    173:   PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
    174:   PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb));

    176:   for (b = 0; b < Nb; ++b) {
    177:     ego body = bodies[b];
    178:     int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;

    180:     /* Output Basic Model Topology */
    181:     PetscCall(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs));
    182:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "   Number of SHELLS: %d \n", Nsh));
    183:     EG_free(objs);

    185:     PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs));
    186:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "   Number of FACES: %d \n", Nf));
    187:     EG_free(objs);

    189:     PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
    190:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "   Number of LOOPS: %d \n", Nl));

    192:     PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs));
    193:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "   Number of EDGES: %d \n", Ne));
    194:     EG_free(objs);

    196:     PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs));
    197:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "   Number of NODES: %d \n", Nv));
    198:     EG_free(objs);

    200:     for (l = 0; l < Nl; ++l) {
    201:       ego loop = lobjs[l];

    203:       id = EG_indexBodyTopo(body, loop);
    204:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "          LOOP ID: %d\n", id));

    206:       /* Get EDGE info which associated with the current LOOP */
    207:       PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));

    209:       for (e = 0; e < Ne; ++e) {
    210:         ego    edge      = objs[e];
    211:         double range[4]  = {0., 0., 0., 0.};
    212:         double point[3]  = {0., 0., 0.};
    213:         double params[3] = {0., 0., 0.};
    214:         double result[18];
    215:         int    peri;

    217:         id = EG_indexBodyTopo(body, edge);
    218:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "            EDGE ID: %d (%d)\n", id, e));

    220:         PetscCall(EG_getRange(edge, range, &peri));
    221:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]));

    223:         /* Get NODE info which associated with the current EDGE */
    224:         PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
    225:         if (mtype == DEGENERATE) {
    226:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  EDGE %d is DEGENERATE \n", id));
    227:         } else {
    228:           params[0] = range[0];
    229:           PetscCall(EG_evaluate(edge, params, result));
    230:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "   between (%lf, %lf, %lf)", result[0], result[1], result[2]));
    231:           params[0] = range[1];
    232:           PetscCall(EG_evaluate(edge, params, result));
    233:           PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]));
    234:         }

    236:         for (v = 0; v < Nv; ++v) {
    237:           ego    vertex = nobjs[v];
    238:           double limits[4];
    239:           int    dummy;

    241:           PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
    242:           id = EG_indexBodyTopo(body, vertex);
    243:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "              NODE ID: %d \n", id));
    244:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "                 (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]));

    246:           point[0] = point[0] + limits[0];
    247:           point[1] = point[1] + limits[1];
    248:           point[2] = point[2] + limits[2];
    249:         }
    250:       }
    251:     }
    252:     EG_free(lobjs);
    253:   }
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

    257: static PetscErrorCode DMPlexEGADSDestroy_Private(void *context)
    258: {
    259:   if (context) EG_close((ego)context);
    260:   return 0;
    261: }

    263: static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
    264: {
    265:   DMLabel  bodyLabel, faceLabel, edgeLabel, vertexLabel;
    266:   PetscInt cStart, cEnd, c;
    267:   /* EGADSLite variables */
    268:   ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
    269:   int oclass, mtype, nbodies, *senses;
    270:   int b;
    271:   /* PETSc variables */
    272:   DM          dm;
    273:   PetscHMapI  edgeMap = NULL;
    274:   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;
    275:   PetscInt   *cells = NULL, *cone = NULL;
    276:   PetscReal  *coords = NULL;
    277:   PetscMPIInt rank;

    279:   PetscFunctionBegin;
    280:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    281:   if (rank == 0) {
    282:     const PetscInt debug = 0;

    284:     /* ---------------------------------------------------------------------------------------------------
    285:     Generate Petsc Plex
    286:       Get all Nodes in model, record coordinates in a correctly formatted array
    287:       Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
    288:       We need to uniformly refine the initial geometry to guarantee a valid mesh
    289:     */

    291:     /* Calculate cell and vertex sizes */
    292:     PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
    293:     PetscCall(PetscHMapICreate(&edgeMap));
    294:     numEdges = 0;
    295:     for (b = 0; b < nbodies; ++b) {
    296:       ego body = bodies[b];
    297:       int id, Nl, l, Nv, v;

    299:       PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
    300:       for (l = 0; l < Nl; ++l) {
    301:         ego loop = lobjs[l];
    302:         int Ner  = 0, Ne, e, Nc;

    304:         PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
    305:         for (e = 0; e < Ne; ++e) {
    306:           ego           edge = objs[e];
    307:           int           Nv, id;
    308:           PetscHashIter iter;
    309:           PetscBool     found;

    311:           PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
    312:           if (mtype == DEGENERATE) continue;
    313:           id = EG_indexBodyTopo(body, edge);
    314:           PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found));
    315:           if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++));
    316:           ++Ner;
    317:         }
    318:         if (Ner == 2) {
    319:           Nc = 2;
    320:         } else if (Ner == 3) {
    321:           Nc = 4;
    322:         } else if (Ner == 4) {
    323:           Nc = 8;
    324:           ++numQuads;
    325:         } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
    326:         numCells += Nc;
    327:         newCells += Nc - 1;
    328:         maxCorners = PetscMax(Ner * 2 + 1, maxCorners);
    329:       }
    330:       PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
    331:       for (v = 0; v < Nv; ++v) {
    332:         ego vertex = nobjs[v];

    334:         id = EG_indexBodyTopo(body, vertex);
    335:         /* TODO: Instead of assuming contiguous ids, we could use a hash table */
    336:         numVertices = PetscMax(id, numVertices);
    337:       }
    338:       EG_free(lobjs);
    339:       EG_free(nobjs);
    340:     }
    341:     PetscCall(PetscHMapIGetSize(edgeMap, &numEdges));
    342:     newVertices = numEdges + numQuads;
    343:     numVertices += newVertices;

    345:     dim        = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
    346:     cdim       = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
    347:     numCorners = 3; /* Split cells into triangles */
    348:     PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone));

    350:     /* Get vertex coordinates */
    351:     for (b = 0; b < nbodies; ++b) {
    352:       ego body = bodies[b];
    353:       int id, Nv, v;

    355:       PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
    356:       for (v = 0; v < Nv; ++v) {
    357:         ego    vertex = nobjs[v];
    358:         double limits[4];
    359:         int    dummy;

    361:         PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
    362:         id                          = EG_indexBodyTopo(body, vertex);
    363:         coords[(id - 1) * cdim + 0] = limits[0];
    364:         coords[(id - 1) * cdim + 1] = limits[1];
    365:         coords[(id - 1) * cdim + 2] = limits[2];
    366:       }
    367:       EG_free(nobjs);
    368:     }
    369:     PetscCall(PetscHMapIClear(edgeMap));
    370:     fOff     = numVertices - newVertices + numEdges;
    371:     numEdges = 0;
    372:     numQuads = 0;
    373:     for (b = 0; b < nbodies; ++b) {
    374:       ego body = bodies[b];
    375:       int Nl, l;

    377:       PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
    378:       for (l = 0; l < Nl; ++l) {
    379:         ego loop = lobjs[l];
    380:         int lid, Ner = 0, Ne, e;

    382:         lid = EG_indexBodyTopo(body, loop);
    383:         PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
    384:         for (e = 0; e < Ne; ++e) {
    385:           ego           edge = objs[e];
    386:           int           eid, Nv;
    387:           PetscHashIter iter;
    388:           PetscBool     found;

    390:           PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
    391:           if (mtype == DEGENERATE) continue;
    392:           ++Ner;
    393:           eid = EG_indexBodyTopo(body, edge);
    394:           PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found));
    395:           if (!found) {
    396:             PetscInt v = numVertices - newVertices + numEdges;
    397:             double   range[4], params[3] = {0., 0., 0.}, result[18];
    398:             int      periodic[2];

    400:             PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++));
    401:             PetscCall(EG_getRange(edge, range, periodic));
    402:             params[0] = 0.5 * (range[0] + range[1]);
    403:             PetscCall(EG_evaluate(edge, params, result));
    404:             coords[v * cdim + 0] = result[0];
    405:             coords[v * cdim + 1] = result[1];
    406:             coords[v * cdim + 2] = result[2];
    407:           }
    408:         }
    409:         if (Ner == 4) {
    410:           PetscInt v = fOff + numQuads++;
    411:           ego     *fobjs, face;
    412:           double   range[4], params[3] = {0., 0., 0.}, result[18];
    413:           int      Nf, fid, periodic[2];

    415:           PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
    416:           face = fobjs[0];
    417:           fid  = EG_indexBodyTopo(body, face);
    418:           PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid);
    419:           PetscCall(EG_getRange(face, range, periodic));
    420:           params[0] = 0.5 * (range[0] + range[1]);
    421:           params[1] = 0.5 * (range[2] + range[3]);
    422:           PetscCall(EG_evaluate(face, params, result));
    423:           coords[v * cdim + 0] = result[0];
    424:           coords[v * cdim + 1] = result[1];
    425:           coords[v * cdim + 2] = result[2];
    426:         }
    427:       }
    428:     }
    429:     PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices);

    431:     /* Get cell vertices by traversing loops */
    432:     numQuads = 0;
    433:     cOff     = 0;
    434:     for (b = 0; b < nbodies; ++b) {
    435:       ego body = bodies[b];
    436:       int id, Nl, l;

    438:       PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
    439:       for (l = 0; l < Nl; ++l) {
    440:         ego loop = lobjs[l];
    441:         int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;

    443:         lid = EG_indexBodyTopo(body, loop);
    444:         PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));

    446:         for (e = 0; e < Ne; ++e) {
    447:           ego edge = objs[e];
    448:           int points[3];
    449:           int eid, Nv, v, tmp;

    451:           eid = EG_indexBodyTopo(body, edge);
    452:           PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
    453:           if (mtype == DEGENERATE) continue;
    454:           else ++Ner;
    455:           PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);

    457:           for (v = 0; v < Nv; ++v) {
    458:             ego vertex = nobjs[v];

    460:             id            = EG_indexBodyTopo(body, vertex);
    461:             points[v * 2] = id - 1;
    462:           }
    463:           {
    464:             PetscInt edgeNum;

    466:             PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
    467:             points[1] = numVertices - newVertices + edgeNum;
    468:           }
    469:           /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
    470:           if (!nc) {
    471:             for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v];
    472:           } else {
    473:             if (cone[nc - 1] == points[0]) {
    474:               cone[nc++] = points[1];
    475:               if (cone[0] != points[2]) cone[nc++] = points[2];
    476:             } else if (cone[nc - 1] == points[2]) {
    477:               cone[nc++] = points[1];
    478:               if (cone[0] != points[0]) cone[nc++] = points[0];
    479:             } else if (cone[nc - 3] == points[0]) {
    480:               tmp          = cone[nc - 3];
    481:               cone[nc - 3] = cone[nc - 1];
    482:               cone[nc - 1] = tmp;
    483:               cone[nc++]   = points[1];
    484:               if (cone[0] != points[2]) cone[nc++] = points[2];
    485:             } else if (cone[nc - 3] == points[2]) {
    486:               tmp          = cone[nc - 3];
    487:               cone[nc - 3] = cone[nc - 1];
    488:               cone[nc - 1] = tmp;
    489:               cone[nc++]   = points[1];
    490:               if (cone[0] != points[0]) cone[nc++] = points[0];
    491:             } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
    492:           }
    493:         }
    494:         PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner);
    495:         if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++;
    496:         PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners);
    497:         /* Triangulate the loop */
    498:         switch (Ner) {
    499:         case 2: /* Bi-Segment -> 2 triangles */
    500:           Nt                           = 2;
    501:           cells[cOff * numCorners + 0] = cone[0];
    502:           cells[cOff * numCorners + 1] = cone[1];
    503:           cells[cOff * numCorners + 2] = cone[2];
    504:           ++cOff;
    505:           cells[cOff * numCorners + 0] = cone[0];
    506:           cells[cOff * numCorners + 1] = cone[2];
    507:           cells[cOff * numCorners + 2] = cone[3];
    508:           ++cOff;
    509:           break;
    510:         case 3: /* Triangle   -> 4 triangles */
    511:           Nt                           = 4;
    512:           cells[cOff * numCorners + 0] = cone[0];
    513:           cells[cOff * numCorners + 1] = cone[1];
    514:           cells[cOff * numCorners + 2] = cone[5];
    515:           ++cOff;
    516:           cells[cOff * numCorners + 0] = cone[1];
    517:           cells[cOff * numCorners + 1] = cone[2];
    518:           cells[cOff * numCorners + 2] = cone[3];
    519:           ++cOff;
    520:           cells[cOff * numCorners + 0] = cone[5];
    521:           cells[cOff * numCorners + 1] = cone[3];
    522:           cells[cOff * numCorners + 2] = cone[4];
    523:           ++cOff;
    524:           cells[cOff * numCorners + 0] = cone[1];
    525:           cells[cOff * numCorners + 1] = cone[3];
    526:           cells[cOff * numCorners + 2] = cone[5];
    527:           ++cOff;
    528:           break;
    529:         case 4: /* Quad       -> 8 triangles */
    530:           Nt                           = 8;
    531:           cells[cOff * numCorners + 0] = cone[0];
    532:           cells[cOff * numCorners + 1] = cone[1];
    533:           cells[cOff * numCorners + 2] = cone[7];
    534:           ++cOff;
    535:           cells[cOff * numCorners + 0] = cone[1];
    536:           cells[cOff * numCorners + 1] = cone[2];
    537:           cells[cOff * numCorners + 2] = cone[3];
    538:           ++cOff;
    539:           cells[cOff * numCorners + 0] = cone[3];
    540:           cells[cOff * numCorners + 1] = cone[4];
    541:           cells[cOff * numCorners + 2] = cone[5];
    542:           ++cOff;
    543:           cells[cOff * numCorners + 0] = cone[5];
    544:           cells[cOff * numCorners + 1] = cone[6];
    545:           cells[cOff * numCorners + 2] = cone[7];
    546:           ++cOff;
    547:           cells[cOff * numCorners + 0] = cone[8];
    548:           cells[cOff * numCorners + 1] = cone[1];
    549:           cells[cOff * numCorners + 2] = cone[3];
    550:           ++cOff;
    551:           cells[cOff * numCorners + 0] = cone[8];
    552:           cells[cOff * numCorners + 1] = cone[3];
    553:           cells[cOff * numCorners + 2] = cone[5];
    554:           ++cOff;
    555:           cells[cOff * numCorners + 0] = cone[8];
    556:           cells[cOff * numCorners + 1] = cone[5];
    557:           cells[cOff * numCorners + 2] = cone[7];
    558:           ++cOff;
    559:           cells[cOff * numCorners + 0] = cone[8];
    560:           cells[cOff * numCorners + 1] = cone[7];
    561:           cells[cOff * numCorners + 2] = cone[1];
    562:           ++cOff;
    563:           break;
    564:         default:
    565:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
    566:         }
    567:         if (debug) {
    568:           for (t = 0; t < Nt; ++t) {
    569:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "  LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t));
    570:             for (c = 0; c < numCorners; ++c) {
    571:               if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    572:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c]));
    573:             }
    574:             PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
    575:           }
    576:         }
    577:       }
    578:       EG_free(lobjs);
    579:     }
    580:   }
    581:   PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells);
    582:   PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
    583:   PetscCall(PetscFree3(coords, cells, cone));
    584:   PetscCall(PetscInfo(dm, " Total Number of Unique Cells    = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells));
    585:   PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices));
    586:   /* Embed EGADS model in DM */
    587:   {
    588:     PetscContainer modelObj, contextObj;

    590:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
    591:     PetscCall(PetscContainerSetPointer(modelObj, model));
    592:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
    593:     PetscCall(PetscContainerDestroy(&modelObj));

    595:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
    596:     PetscCall(PetscContainerSetPointer(contextObj, context));
    597:     PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
    598:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
    599:     PetscCall(PetscContainerDestroy(&contextObj));
    600:   }
    601:   /* Label points */
    602:   PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
    603:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    604:   PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
    605:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    606:   PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
    607:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    608:   PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
    609:   PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
    610:   cOff = 0;
    611:   for (b = 0; b < nbodies; ++b) {
    612:     ego body = bodies[b];
    613:     int id, Nl, l;

    615:     PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
    616:     for (l = 0; l < Nl; ++l) {
    617:       ego  loop = lobjs[l];
    618:       ego *fobjs;
    619:       int  lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;

    621:       lid = EG_indexBodyTopo(body, loop);
    622:       PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
    623:       PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
    624:       fid = EG_indexBodyTopo(body, fobjs[0]);
    625:       EG_free(fobjs);
    626:       PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
    627:       for (e = 0; e < Ne; ++e) {
    628:         ego             edge = objs[e];
    629:         int             eid, Nv, v;
    630:         PetscInt        points[3], support[2], numEdges, edgeNum;
    631:         const PetscInt *edges;

    633:         eid = EG_indexBodyTopo(body, edge);
    634:         PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
    635:         if (mtype == DEGENERATE) continue;
    636:         else ++Ner;
    637:         for (v = 0; v < Nv; ++v) {
    638:           ego vertex = nobjs[v];

    640:           id = EG_indexBodyTopo(body, vertex);
    641:           PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid));
    642:           points[v * 2] = numCells + id - 1;
    643:         }
    644:         PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
    645:         points[1] = numCells + numVertices - newVertices + edgeNum;

    647:         PetscCall(DMLabelSetValue(edgeLabel, points[1], eid));
    648:         support[0] = points[0];
    649:         support[1] = points[1];
    650:         PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
    651:         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);
    652:         PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
    653:         PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
    654:         support[0] = points[1];
    655:         support[1] = points[2];
    656:         PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
    657:         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);
    658:         PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
    659:         PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
    660:       }
    661:       switch (Ner) {
    662:       case 2:
    663:         Nt = 2;
    664:         break;
    665:       case 3:
    666:         Nt = 4;
    667:         break;
    668:       case 4:
    669:         Nt = 8;
    670:         break;
    671:       default:
    672:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
    673:       }
    674:       for (t = 0; t < Nt; ++t) {
    675:         PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b));
    676:         PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid));
    677:       }
    678:       cOff += Nt;
    679:     }
    680:     EG_free(lobjs);
    681:   }
    682:   PetscCall(PetscHMapIDestroy(&edgeMap));
    683:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
    684:   for (c = cStart; c < cEnd; ++c) {
    685:     PetscInt *closure = NULL;
    686:     PetscInt  clSize, cl, bval, fval;

    688:     PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
    689:     PetscCall(DMLabelGetValue(bodyLabel, c, &bval));
    690:     PetscCall(DMLabelGetValue(faceLabel, c, &fval));
    691:     for (cl = 0; cl < clSize * 2; cl += 2) {
    692:       PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval));
    693:       PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval));
    694:     }
    695:     PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
    696:   }
    697:   *newdm = dm;
    698:   PetscFunctionReturn(PETSC_SUCCESS);
    699: }

    701: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
    702: {
    703:   DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
    704:   // EGADS/EGADSLite variables
    705:   ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs;
    706:   ego topRef, prev, next;
    707:   int oclass, mtype, nbodies, *senses, *lSenses, *eSenses;
    708:   int b;
    709:   // PETSc variables
    710:   DM              dm;
    711:   PetscHMapI      edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL;
    712:   PetscInt        dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0;
    713:   PetscInt        cellCntr = 0, numPoints = 0;
    714:   PetscInt       *cells  = NULL;
    715:   const PetscInt *cone   = NULL;
    716:   PetscReal      *coords = NULL;
    717:   PetscMPIInt     rank;

    719:   PetscFunctionBeginUser;
    720:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    721:   if (rank == 0) {
    722:     // ---------------------------------------------------------------------------------------------------
    723:     // Generate Petsc Plex
    724:     //  Get all Nodes in model, record coordinates in a correctly formatted array
    725:     //  Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
    726:     //  We need to uniformly refine the initial geometry to guarantee a valid mesh

    728:     // Calculate cell and vertex sizes
    729:     PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));

    731:     PetscCall(PetscHMapICreate(&edgeMap));
    732:     PetscCall(PetscHMapICreate(&bodyIndexMap));
    733:     PetscCall(PetscHMapICreate(&bodyVertexMap));
    734:     PetscCall(PetscHMapICreate(&bodyEdgeMap));
    735:     PetscCall(PetscHMapICreate(&bodyEdgeGlobalMap));
    736:     PetscCall(PetscHMapICreate(&bodyFaceMap));

    738:     for (b = 0; b < nbodies; ++b) {
    739:       ego           body = bodies[b];
    740:       int           Nf, Ne, Nv;
    741:       PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter;
    742:       PetscBool     BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound;

    744:       PetscCall(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound));
    745:       PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
    746:       PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
    747:       PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
    748:       PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));

    750:       if (!BIfound) PetscCall(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices));
    751:       if (!BVfound) PetscCall(PetscHMapISet(bodyVertexMap, b, numVertices));
    752:       if (!BEfound) PetscCall(PetscHMapISet(bodyEdgeMap, b, numEdges));
    753:       if (!BEGfound) PetscCall(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr));
    754:       if (!BFfound) PetscCall(PetscHMapISet(bodyFaceMap, b, numFaces));

    756:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
    757:       PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
    758:       PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
    759:       EG_free(fobjs);
    760:       EG_free(eobjs);
    761:       EG_free(nobjs);

    763:       // Remove DEGENERATE EDGES from Edge count
    764:       PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
    765:       int Netemp = 0;
    766:       for (int e = 0; e < Ne; ++e) {
    767:         ego edge = eobjs[e];
    768:         int eid;

    770:         PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
    771:         eid = EG_indexBodyTopo(body, edge);

    773:         PetscCall(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound));
    774:         if (mtype == DEGENERATE) {
    775:           if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1));
    776:         } else {
    777:           ++Netemp;
    778:           if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp));
    779:         }
    780:       }
    781:       EG_free(eobjs);

    783:       // Determine Number of Cells
    784:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
    785:       for (int f = 0; f < Nf; ++f) {
    786:         ego face     = fobjs[f];
    787:         int edgeTemp = 0;

    789:         PetscCall(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs));
    790:         for (int e = 0; e < Ne; ++e) {
    791:           ego edge = eobjs[e];

    793:           PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
    794:           if (mtype != DEGENERATE) ++edgeTemp;
    795:         }
    796:         numCells += (2 * edgeTemp);
    797:         EG_free(eobjs);
    798:       }
    799:       EG_free(fobjs);

    801:       numFaces += Nf;
    802:       numEdges += Netemp;
    803:       numVertices += Nv;
    804:       edgeCntr += Ne;
    805:     }

    807:     // Set up basic DMPlex parameters
    808:     dim        = 2;                                 // Assumes 3D Models :: Need to handle 2D models in the future
    809:     cdim       = 3;                                 // Assumes 3D Models :: Need to update to handle 2D models in future
    810:     numCorners = 3;                                 // Split Faces into triangles
    811:     numPoints  = numVertices + numEdges + numFaces; // total number of coordinate points

    813:     PetscCall(PetscMalloc2(numPoints * cdim, &coords, numCells * numCorners, &cells));

    815:     // Get Vertex Coordinates and Set up Cells
    816:     for (b = 0; b < nbodies; ++b) {
    817:       ego           body = bodies[b];
    818:       int           Nf, Ne, Nv;
    819:       PetscInt      bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
    820:       PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
    821:       PetscBool     BVfound, BEfound, BEGfound, BFfound, EMfound;

    823:       // Vertices on Current Body
    824:       PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));

    826:       PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
    827:       PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
    828:       PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));

    830:       for (int v = 0; v < Nv; ++v) {
    831:         ego    vertex = nobjs[v];
    832:         double limits[4];
    833:         int    id, dummy;

    835:         PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
    836:         id = EG_indexBodyTopo(body, vertex);

    838:         coords[(bodyVertexIndexStart + id - 1) * cdim + 0] = limits[0];
    839:         coords[(bodyVertexIndexStart + id - 1) * cdim + 1] = limits[1];
    840:         coords[(bodyVertexIndexStart + id - 1) * cdim + 2] = limits[2];
    841:       }
    842:       EG_free(nobjs);

    844:       // Edge Midpoint Vertices on Current Body
    845:       PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));

    847:       PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
    848:       PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
    849:       PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));

    851:       PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
    852:       PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
    853:       PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));

    855:       for (int e = 0; e < Ne; ++e) {
    856:         ego    edge = eobjs[e];
    857:         double range[2], avgt[1], cntrPnt[9];
    858:         int    eid, eOffset;
    859:         int    periodic;

    861:         PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
    862:         if (mtype == DEGENERATE) continue;

    864:         eid = EG_indexBodyTopo(body, edge);

    866:         // get relative offset from globalEdgeID Vector
    867:         PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
    868:         PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1);
    869:         PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));

    871:         PetscCall(EG_getRange(edge, range, &periodic));
    872:         avgt[0] = (range[0] + range[1]) / 2.;

    874:         PetscCall(EG_evaluate(edge, avgt, cntrPnt));
    875:         coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 0] = cntrPnt[0];
    876:         coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 1] = cntrPnt[1];
    877:         coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 2] = cntrPnt[2];
    878:       }
    879:       EG_free(eobjs);

    881:       // Face Midpoint Vertices on Current Body
    882:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));

    884:       PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
    885:       PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
    886:       PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));

    888:       for (int f = 0; f < Nf; ++f) {
    889:         ego    face = fobjs[f];
    890:         double range[4], avgUV[2], cntrPnt[18];
    891:         int    peri, id;

    893:         id = EG_indexBodyTopo(body, face);
    894:         PetscCall(EG_getRange(face, range, &peri));

    896:         avgUV[0] = (range[0] + range[1]) / 2.;
    897:         avgUV[1] = (range[2] + range[3]) / 2.;
    898:         PetscCall(EG_evaluate(face, avgUV, cntrPnt));

    900:         coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 0] = cntrPnt[0];
    901:         coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 1] = cntrPnt[1];
    902:         coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 2] = cntrPnt[2];
    903:       }
    904:       EG_free(fobjs);

    906:       // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity
    907:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
    908:       for (int f = 0; f < Nf; ++f) {
    909:         ego face = fobjs[f];
    910:         int fID, midFaceID, midPntID, startID, endID, Nl;

    912:         fID       = EG_indexBodyTopo(body, face);
    913:         midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1;
    914:         // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges.
    915:         // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are:
    916:         //            1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option.
    917:         //               This will use a default EGADS tessellation as an initial surface mesh.
    918:         //            2) Create the initial surface mesh via a 2D mesher :: Currently not available (?future?)
    919:         //               May I suggest the XXXX as a starting point?

    921:         PetscCall(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses));

    923:         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);
    924:         for (int l = 0; l < Nl; ++l) {
    925:           ego loop = lobjs[l];

    927:           PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
    928:           for (int e = 0; e < Ne; ++e) {
    929:             ego edge = eobjs[e];
    930:             int eid, eOffset;

    932:             PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
    933:             eid = EG_indexBodyTopo(body, edge);
    934:             if (mtype == DEGENERATE) continue;

    936:             // get relative offset from globalEdgeID Vector
    937:             PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
    938:             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);
    939:             PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));

    941:             midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1;

    943:             PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));

    945:             if (eSenses[e] > 0) {
    946:               startID = EG_indexBodyTopo(body, nobjs[0]);
    947:               endID   = EG_indexBodyTopo(body, nobjs[1]);
    948:             } else {
    949:               startID = EG_indexBodyTopo(body, nobjs[1]);
    950:               endID   = EG_indexBodyTopo(body, nobjs[0]);
    951:             }

    953:             // Define 2 Cells per Edge with correct orientation
    954:             cells[cellCntr * numCorners + 0] = midFaceID;
    955:             cells[cellCntr * numCorners + 1] = bodyVertexIndexStart + startID - 1;
    956:             cells[cellCntr * numCorners + 2] = midPntID;

    958:             cells[cellCntr * numCorners + 3] = midFaceID;
    959:             cells[cellCntr * numCorners + 4] = midPntID;
    960:             cells[cellCntr * numCorners + 5] = bodyVertexIndexStart + endID - 1;

    962:             cellCntr = cellCntr + 2;
    963:           }
    964:         }
    965:       }
    966:       EG_free(fobjs);
    967:     }
    968:   }

    970:   // Generate DMPlex
    971:   PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
    972:   PetscCall(PetscFree2(coords, cells));
    973:   PetscCall(PetscInfo(dm, " Total Number of Unique Cells    = %" PetscInt_FMT " \n", numCells));
    974:   PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " \n", numVertices));

    976:   // Embed EGADS model in DM
    977:   {
    978:     PetscContainer modelObj, contextObj;

    980:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
    981:     PetscCall(PetscContainerSetPointer(modelObj, model));
    982:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
    983:     PetscCall(PetscContainerDestroy(&modelObj));

    985:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
    986:     PetscCall(PetscContainerSetPointer(contextObj, context));
    987:     PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
    988:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
    989:     PetscCall(PetscContainerDestroy(&contextObj));
    990:   }
    991:   // Label points
    992:   PetscInt nStart, nEnd;

    994:   PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
    995:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    996:   PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
    997:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    998:   PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
    999:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    1000:   PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
    1001:   PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));

    1003:   PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));

    1005:   cellCntr = 0;
    1006:   for (b = 0; b < nbodies; ++b) {
    1007:     ego           body = bodies[b];
    1008:     int           Nv, Ne, Nf;
    1009:     PetscInt      bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
    1010:     PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
    1011:     PetscBool     BVfound, BEfound, BEGfound, BFfound, EMfound;

    1013:     PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
    1014:     PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
    1015:     PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));

    1017:     PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
    1018:     PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
    1019:     PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));

    1021:     PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
    1022:     PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
    1023:     PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));

    1025:     PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
    1026:     PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
    1027:     PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));

    1029:     PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
    1030:     for (int f = 0; f < Nf; ++f) {
    1031:       ego face = fobjs[f];
    1032:       int fID, Nl;

    1034:       fID = EG_indexBodyTopo(body, face);

    1036:       PetscCall(EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs));
    1037:       for (int l = 0; l < Nl; ++l) {
    1038:         ego loop = lobjs[l];
    1039:         int lid;

    1041:         lid = EG_indexBodyTopo(body, loop);
    1042:         PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);

    1044:         PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
    1045:         for (int e = 0; e < Ne; ++e) {
    1046:           ego edge = eobjs[e];
    1047:           int eid, eOffset;

    1049:           // Skip DEGENERATE Edges
    1050:           PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
    1051:           if (mtype == DEGENERATE) continue;
    1052:           eid = EG_indexBodyTopo(body, edge);

    1054:           // get relative offset from globalEdgeID Vector
    1055:           PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
    1056:           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);
    1057:           PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));

    1059:           PetscCall(EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs));
    1060:           for (int v = 0; v < Nv; ++v) {
    1061:             ego vertex = nobjs[v];
    1062:             int vID;

    1064:             vID = EG_indexBodyTopo(body, vertex);
    1065:             PetscCall(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b));
    1066:             PetscCall(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID));
    1067:           }
    1068:           EG_free(nobjs);

    1070:           PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b));
    1071:           PetscCall(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid));

    1073:           // Define Cell faces
    1074:           for (int jj = 0; jj < 2; ++jj) {
    1075:             PetscCall(DMLabelSetValue(bodyLabel, cellCntr, b));
    1076:             PetscCall(DMLabelSetValue(faceLabel, cellCntr, fID));
    1077:             PetscCall(DMPlexGetCone(dm, cellCntr, &cone));

    1079:             PetscCall(DMLabelSetValue(bodyLabel, cone[0], b));
    1080:             PetscCall(DMLabelSetValue(faceLabel, cone[0], fID));

    1082:             PetscCall(DMLabelSetValue(bodyLabel, cone[1], b));
    1083:             PetscCall(DMLabelSetValue(edgeLabel, cone[1], eid));

    1085:             PetscCall(DMLabelSetValue(bodyLabel, cone[2], b));
    1086:             PetscCall(DMLabelSetValue(faceLabel, cone[2], fID));

    1088:             cellCntr = cellCntr + 1;
    1089:           }
    1090:         }
    1091:       }
    1092:       EG_free(lobjs);

    1094:       PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b));
    1095:       PetscCall(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID));
    1096:     }
    1097:     EG_free(fobjs);
    1098:   }

    1100:   PetscCall(PetscHMapIDestroy(&edgeMap));
    1101:   PetscCall(PetscHMapIDestroy(&bodyIndexMap));
    1102:   PetscCall(PetscHMapIDestroy(&bodyVertexMap));
    1103:   PetscCall(PetscHMapIDestroy(&bodyEdgeMap));
    1104:   PetscCall(PetscHMapIDestroy(&bodyEdgeGlobalMap));
    1105:   PetscCall(PetscHMapIDestroy(&bodyFaceMap));

    1107:   *newdm = dm;
    1108:   PetscFunctionReturn(PETSC_SUCCESS);
    1109: }

    1111: static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
    1112: {
    1113:   DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
    1114:   /* EGADSLite variables */
    1115:   ego    geom, *bodies, *fobjs;
    1116:   int    b, oclass, mtype, nbodies, *senses;
    1117:   int    totalNumTris = 0, totalNumPoints = 0;
    1118:   double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize;
    1119:   /* PETSc variables */
    1120:   DM              dm;
    1121:   PetscHMapI      pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL;
    1122:   PetscHMapI      pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL;
    1123:   PetscInt        dim = -1, cdim = -1, numCorners = 0, counter = 0;
    1124:   PetscInt       *cells  = NULL;
    1125:   const PetscInt *cone   = NULL;
    1126:   PetscReal      *coords = NULL;
    1127:   PetscMPIInt     rank;

    1129:   PetscFunctionBeginUser;
    1130:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1131:   if (rank == 0) {
    1132:     // ---------------------------------------------------------------------------------------------------
    1133:     // Generate Petsc Plex from EGADSlite created Tessellation of geometry
    1134:     // ---------------------------------------------------------------------------------------------------

    1136:     // Calculate cell and vertex sizes
    1137:     PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));

    1139:     PetscCall(PetscHMapICreate(&pointIndexStartMap));
    1140:     PetscCall(PetscHMapICreate(&triIndexStartMap));
    1141:     PetscCall(PetscHMapICreate(&pTypeLabelMap));
    1142:     PetscCall(PetscHMapICreate(&pIndexLabelMap));
    1143:     PetscCall(PetscHMapICreate(&pBodyIndexLabelMap));
    1144:     PetscCall(PetscHMapICreate(&triFaceIDLabelMap));
    1145:     PetscCall(PetscHMapICreate(&triBodyIDLabelMap));

    1147:     /* Create Tessellation of Bodies */
    1148:     ego tessArray[nbodies];

    1150:     for (b = 0; b < nbodies; ++b) {
    1151:       ego           body      = bodies[b];
    1152:       double        params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation
    1153:       int           Nf, bodyNumPoints = 0, bodyNumTris = 0;
    1154:       PetscHashIter PISiter, TISiter;
    1155:       PetscBool     PISfound, TISfound;

    1157:       /* Store Start Index for each Body's Point and Tris */
    1158:       PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
    1159:       PetscCall(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound));

    1161:       if (!PISfound) PetscCall(PetscHMapISet(pointIndexStartMap, b, totalNumPoints));
    1162:       if (!TISfound) PetscCall(PetscHMapISet(triIndexStartMap, b, totalNumTris));

    1164:       /* Calculate Tessellation parameters based on Bounding Box */
    1165:       /* Get Bounding Box Dimensions of the BODY */
    1166:       PetscCall(EG_getBoundingBox(body, boundBox));
    1167:       tessSize = boundBox[3] - boundBox[0];
    1168:       if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1];
    1169:       if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2];

    1171:       // TODO :: May want to give users tessellation parameter options //
    1172:       params[0] = 0.0250 * tessSize;
    1173:       params[1] = 0.0075 * tessSize;
    1174:       params[2] = 15.0;

    1176:       PetscCall(EG_makeTessBody(body, params, &tessArray[b]));

    1178:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));

    1180:       for (int f = 0; f < Nf; ++f) {
    1181:         ego           face = fobjs[f];
    1182:         int           len, fID, ntris;
    1183:         const int    *ptype, *pindex, *ptris, *ptric;
    1184:         const double *pxyz, *puv;

    1186:         // Get Face ID //
    1187:         fID = EG_indexBodyTopo(body, face);

    1189:         // Checkout the Surface Tessellation //
    1190:         PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));

    1192:         // Determine total number of triangle cells in the tessellation //
    1193:         bodyNumTris += (int)ntris;

    1195:         // Check out the point index and coordinate //
    1196:         for (int p = 0; p < len; ++p) {
    1197:           int global;

    1199:           PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));

    1201:           // Determine the total number of points in the tessellation //
    1202:           bodyNumPoints = PetscMax(bodyNumPoints, global);
    1203:         }
    1204:       }
    1205:       EG_free(fobjs);

    1207:       totalNumPoints += bodyNumPoints;
    1208:       totalNumTris += bodyNumTris;
    1209:     }
    1210:     //}  - Original End of (rank == 0)

    1212:     dim        = 2;
    1213:     cdim       = 3;
    1214:     numCorners = 3;
    1215:     //PetscInt counter = 0;

    1217:     /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA   */
    1218:     /* Fill in below and use to define DMLabels after DMPlex creation */
    1219:     PetscCall(PetscMalloc2(totalNumPoints * cdim, &coords, totalNumTris * numCorners, &cells));

    1221:     for (b = 0; b < nbodies; ++b) {
    1222:       ego           body = bodies[b];
    1223:       int           Nf;
    1224:       PetscInt      pointIndexStart;
    1225:       PetscHashIter PISiter;
    1226:       PetscBool     PISfound;

    1228:       PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
    1229:       PetscCheck(PISfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b);
    1230:       PetscCall(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart));

    1232:       PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));

    1234:       for (int f = 0; f < Nf; ++f) {
    1235:         /* Get Face Object */
    1236:         ego           face = fobjs[f];
    1237:         int           len, fID, ntris;
    1238:         const int    *ptype, *pindex, *ptris, *ptric;
    1239:         const double *pxyz, *puv;

    1241:         /* Get Face ID */
    1242:         fID = EG_indexBodyTopo(body, face);

    1244:         /* Checkout the Surface Tessellation */
    1245:         PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));

    1247:         /* Check out the point index and coordinate */
    1248:         for (int p = 0; p < len; ++p) {
    1249:           int           global;
    1250:           PetscHashIter PTLiter, PILiter, PBLiter;
    1251:           PetscBool     PTLfound, PILfound, PBLfound;

    1253:           PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));

    1255:           /* Set the coordinates array for DAG */
    1256:           coords[((global - 1 + pointIndexStart) * 3) + 0] = pxyz[(p * 3) + 0];
    1257:           coords[((global - 1 + pointIndexStart) * 3) + 1] = pxyz[(p * 3) + 1];
    1258:           coords[((global - 1 + pointIndexStart) * 3) + 2] = pxyz[(p * 3) + 2];

    1260:           /* Store Geometry Label Information for DMLabel assignment later */
    1261:           PetscCall(PetscHMapIFind(pTypeLabelMap, global - 1 + pointIndexStart, &PTLiter, &PTLfound));
    1262:           PetscCall(PetscHMapIFind(pIndexLabelMap, global - 1 + pointIndexStart, &PILiter, &PILfound));
    1263:           PetscCall(PetscHMapIFind(pBodyIndexLabelMap, global - 1 + pointIndexStart, &PBLiter, &PBLfound));

    1265:           if (!PTLfound) PetscCall(PetscHMapISet(pTypeLabelMap, global - 1 + pointIndexStart, ptype[p]));
    1266:           if (!PILfound) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, pindex[p]));
    1267:           if (!PBLfound) PetscCall(PetscHMapISet(pBodyIndexLabelMap, global - 1 + pointIndexStart, b));

    1269:           if (ptype[p] < 0) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, fID));
    1270:         }

    1272:         for (int t = 0; t < (int)ntris; ++t) {
    1273:           int           global, globalA, globalB;
    1274:           PetscHashIter TFLiter, TBLiter;
    1275:           PetscBool     TFLfound, TBLfound;

    1277:           PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 0], &global));
    1278:           cells[(counter * 3) + 0] = global - 1 + pointIndexStart;

    1280:           PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 1], &globalA));
    1281:           cells[(counter * 3) + 1] = globalA - 1 + pointIndexStart;

    1283:           PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 2], &globalB));
    1284:           cells[(counter * 3) + 2] = globalB - 1 + pointIndexStart;

    1286:           PetscCall(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound));
    1287:           PetscCall(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound));

    1289:           if (!TFLfound) PetscCall(PetscHMapISet(triFaceIDLabelMap, counter, fID));
    1290:           if (!TBLfound) PetscCall(PetscHMapISet(triBodyIDLabelMap, counter, b));

    1292:           counter += 1;
    1293:         }
    1294:       }
    1295:       EG_free(fobjs);
    1296:     }
    1297:   }

    1299:   //Build DMPlex
    1300:   PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
    1301:   PetscCall(PetscFree2(coords, cells));

    1303:   // Embed EGADS model in DM
    1304:   {
    1305:     PetscContainer modelObj, contextObj;

    1307:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
    1308:     PetscCall(PetscContainerSetPointer(modelObj, model));
    1309:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
    1310:     PetscCall(PetscContainerDestroy(&modelObj));

    1312:     PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
    1313:     PetscCall(PetscContainerSetPointer(contextObj, context));
    1314:     PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
    1315:     PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
    1316:     PetscCall(PetscContainerDestroy(&contextObj));
    1317:   }

    1319:   // Label Points
    1320:   PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
    1321:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    1322:   PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
    1323:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    1324:   PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
    1325:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    1326:   PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
    1327:   PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));

    1329:   /* Get Number of DAG Nodes at each level */
    1330:   int fStart, fEnd, eStart, eEnd, nStart, nEnd;

    1332:   PetscCall(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd));
    1333:   PetscCall(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd));
    1334:   PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));

    1336:   /* Set DMLabels for NODES */
    1337:   for (int n = nStart; n < nEnd; ++n) {
    1338:     int           pTypeVal, pIndexVal, pBodyVal;
    1339:     PetscHashIter PTLiter, PILiter, PBLiter;
    1340:     PetscBool     PTLfound, PILfound, PBLfound;

    1342:     //Converted to Hash Tables
    1343:     PetscCall(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound));
    1344:     PetscCheck(PTLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n);
    1345:     PetscCall(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal));

    1347:     PetscCall(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound));
    1348:     PetscCheck(PILfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n);
    1349:     PetscCall(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal));

    1351:     PetscCall(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound));
    1352:     PetscCheck(PBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n);
    1353:     PetscCall(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal));

    1355:     PetscCall(DMLabelSetValue(bodyLabel, n, pBodyVal));
    1356:     if (pTypeVal == 0) PetscCall(DMLabelSetValue(vertexLabel, n, pIndexVal));
    1357:     if (pTypeVal > 0) PetscCall(DMLabelSetValue(edgeLabel, n, pIndexVal));
    1358:     if (pTypeVal < 0) PetscCall(DMLabelSetValue(faceLabel, n, pIndexVal));
    1359:   }

    1361:   /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */
    1362:   for (int e = eStart; e < eEnd; ++e) {
    1363:     int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1;

    1365:     PetscCall(DMPlexGetCone(dm, e, &cone));
    1366:     PetscCall(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end?
    1367:     PetscCall(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0));
    1368:     PetscCall(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1));
    1369:     PetscCall(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0));
    1370:     PetscCall(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1));
    1371:     PetscCall(DMLabelGetValue(faceLabel, cone[0], &faceID_0));
    1372:     PetscCall(DMLabelGetValue(faceLabel, cone[1], &faceID_1));

    1374:     PetscCall(DMLabelSetValue(bodyLabel, e, bodyID_0));

    1376:     if (edgeID_0 == edgeID_1) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
    1377:     else if (vertexID_0 > 0 && edgeID_1 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_1));
    1378:     else if (vertexID_1 > 0 && edgeID_0 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
    1379:     else { /* Do Nothing */ }
    1380:   }

    1382:   /* Set DMLabels for Cells */
    1383:   for (int f = fStart; f < fEnd; ++f) {
    1384:     int           edgeID_0;
    1385:     PetscInt      triBodyVal, triFaceVal;
    1386:     PetscHashIter TFLiter, TBLiter;
    1387:     PetscBool     TFLfound, TBLfound;

    1389:     // Convert to Hash Table
    1390:     PetscCall(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound));
    1391:     PetscCheck(TFLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f);
    1392:     PetscCall(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal));

    1394:     PetscCall(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound));
    1395:     PetscCheck(TBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f);
    1396:     PetscCall(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal));

    1398:     PetscCall(DMLabelSetValue(bodyLabel, f, triBodyVal));
    1399:     PetscCall(DMLabelSetValue(faceLabel, f, triFaceVal));

    1401:     /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */
    1402:     PetscCall(DMPlexGetCone(dm, f, &cone));

    1404:     for (int jj = 0; jj < 3; ++jj) {
    1405:       PetscCall(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0));

    1407:       if (edgeID_0 < 0) {
    1408:         PetscCall(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal));
    1409:         PetscCall(DMLabelSetValue(faceLabel, cone[jj], triFaceVal));
    1410:       }
    1411:     }
    1412:   }

    1414:   *newdm = dm;
    1415:   PetscFunctionReturn(PETSC_SUCCESS);
    1416: }
    1417: #endif

    1419: /*@
    1420:   DMPlexInflateToGeomModel - Snaps the vertex coordinates of a `DMPLEX` object representing the mesh to its geometry if some vertices depart from the model.
    1421:   This usually happens with non-conforming refinement.

    1423:   Collective

    1425:   Input Parameter:
    1426: . dm - The uninflated `DM` object representing the mesh

    1428:   Level: intermediate

    1430: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`
    1431: @*/
    1432: PetscErrorCode DMPlexInflateToGeomModel(DM dm)
    1433: {
    1434: #if defined(PETSC_HAVE_EGADS)
    1435:   /* EGADS Variables */
    1436:   ego    model, geom, body, face, edge;
    1437:   ego   *bodies;
    1438:   int    Nb, oclass, mtype, *senses;
    1439:   double result[3];
    1440:   /* PETSc Variables */
    1441:   DM             cdm;
    1442:   PetscContainer modelObj;
    1443:   DMLabel        bodyLabel, faceLabel, edgeLabel, vertexLabel;
    1444:   Vec            coordinates;
    1445:   PetscScalar   *coords;
    1446:   PetscInt       bodyID, faceID, edgeID, vertexID;
    1447:   PetscInt       cdim, d, vStart, vEnd, v;
    1448: #endif

    1450:   PetscFunctionBegin;
    1451: #if defined(PETSC_HAVE_EGADS)
    1452:   PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
    1453:   if (!modelObj) PetscFunctionReturn(PETSC_SUCCESS);
    1454:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    1455:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1456:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
    1457:   PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
    1458:   PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
    1459:   PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
    1460:   PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));

    1462:   PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
    1463:   PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));

    1465:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
    1466:   PetscCall(VecGetArrayWrite(coordinates, &coords));
    1467:   for (v = vStart; v < vEnd; ++v) {
    1468:     PetscScalar *vcoords;

    1470:     PetscCall(DMLabelGetValue(bodyLabel, v, &bodyID));
    1471:     PetscCall(DMLabelGetValue(faceLabel, v, &faceID));
    1472:     PetscCall(DMLabelGetValue(edgeLabel, v, &edgeID));
    1473:     PetscCall(DMLabelGetValue(vertexLabel, v, &vertexID));

    1475:     PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
    1476:     body = bodies[bodyID];

    1478:     PetscCall(DMPlexPointLocalRef(cdm, v, coords, (void *)&vcoords));
    1479:     if (edgeID > 0) {
    1480:       /* Snap to EDGE at nearest location */
    1481:       double params[1];
    1482:       PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &edge));
    1483:       PetscCall(EG_invEvaluate(edge, vcoords, params, result)); // Get (x,y,z) of nearest point on EDGE
    1484:       for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
    1485:     } else if (faceID > 0) {
    1486:       /* Snap to FACE at nearest location */
    1487:       double params[2];
    1488:       PetscCall(EG_objectBodyTopo(body, FACE, faceID, &face));
    1489:       PetscCall(EG_invEvaluate(face, vcoords, params, result)); // Get (x,y,z) of nearest point on FACE
    1490:       for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
    1491:     }
    1492:   }
    1493:   PetscCall(VecRestoreArrayWrite(coordinates, &coords));
    1494:   /* Clear out global coordinates */
    1495:   PetscCall(VecDestroy(&dm->coordinates[0].x));
    1496: #endif
    1497:   PetscFunctionReturn(PETSC_SUCCESS);
    1498: }

    1500: /*@
    1501:   DMPlexCreateEGADSFromFile - Create a `DMPLEX` mesh from an EGADS, IGES, or STEP file.

    1503:   Collective

    1505:   Input Parameters:
    1506: + comm     - The MPI communicator
    1507: - filename - The name of the EGADS, IGES, or STEP file

    1509:   Output Parameter:
    1510: . dm - The `DM` object representing the mesh

    1512:   Level: beginner

    1514: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSLiteFromFile()`
    1515: @*/
    1516: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
    1517: {
    1518:   PetscMPIInt rank;
    1519: #if defined(PETSC_HAVE_EGADS)
    1520:   ego context = NULL, model = NULL;
    1521: #endif
    1522:   PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE;

    1524:   PetscFunctionBegin;
    1525:   PetscAssertPointer(filename, 2);
    1526:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL));
    1527:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL));
    1528:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL));
    1529:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
    1530: #if defined(PETSC_HAVE_EGADS)
    1531:   if (rank == 0) {
    1532:     PetscCall(EG_open(&context));
    1533:     PetscCall(EG_loadModel(context, 0, filename, &model));
    1534:     if (printModel) PetscCall(DMPlexEGADSPrintModel_Internal(model));
    1535:   }
    1536:   if (tessModel) PetscCall(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm));
    1537:   else if (newModel) PetscCall(DMPlexCreateEGADS(comm, context, model, dm));
    1538:   else PetscCall(DMPlexCreateEGADS_Internal(comm, context, model, dm));
    1539:   PetscFunctionReturn(PETSC_SUCCESS);
    1540: #else
    1541:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads");
    1542: #endif
    1543: }