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 DMPlexSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
 23: PETSC_INTERN PetscErrorCode DMPlexSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);

 25: PetscErrorCode DMPlexSnapToGeomModel_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;

 40:   DMGetCoordinateDM(dm, &cdm);
 41:   DMGetCoordinateDim(dm, &dE);
 42:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
 43:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 44:   if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
 45:   body = bodies[bodyID];

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

 54:   /* Calculate parameters (t or u,v) for vertices */
 55:   DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 56:   Nv  /= dE;
 57:   if (Nv == 1) {
 58:     DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 59:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 60:     return(0);
 61:   }
 62:   if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);

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

 95: /*@
 96:   DMPlexSnapToGeomModel - Given a coordinate point 'mcoords' on the mesh point 'p', return the closest coordinate point 'gcoords' on the geometry model associated with that point.

 98:   Not collective

100:   Input Parameters:
101: + dm      - The DMPlex object
102: . p       - The mesh point
103: - mcoords - A coordinate point lying on the mesh point

105:   Output Parameter:
106: . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point

108:   Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS.

110:   Level: intermediate

112: .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform()
113: @*/
114: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[])
115: {
116:   PetscInt       dE, d;

120:   DMGetCoordinateDim(dm, &dE);
121: #ifdef PETSC_HAVE_EGADS
122:   {
123:     DM_Plex       *plex = (DM_Plex *) dm->data;
124:     DMLabel        bodyLabel, faceLabel, edgeLabel;
125:     PetscInt       bodyID, faceID, edgeID;
126:     PetscContainer modelObj;
127:     ego            model;
128:     PetscBool      islite = PETSC_FALSE;

130:     DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
131:     DMGetLabel(dm, "EGADS Face ID", &faceLabel);
132:     DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
133:     if (!bodyLabel || !faceLabel || !edgeLabel || plex->ignoreModel) {
134:       for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
135:       return(0);
136:     }
137:     PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
138:     if (!modelObj) {
139:       PetscObjectQuery((PetscObject) dm, "EGADSLite Model", (PetscObject *) &modelObj);
140:       islite = PETSC_TRUE;
141:     }
142:     if (!modelObj) {
143:       for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
144:       return(0);
145:     }
146:     PetscContainerGetPointer(modelObj, (void **) &model);
147:     DMLabelGetValue(bodyLabel, p, &bodyID);
148:     DMLabelGetValue(faceLabel, p, &faceID);
149:     DMLabelGetValue(edgeLabel, p, &edgeID);
150:     /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
151:     if (bodyID < 0) {
152:       for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
153:       return(0);
154:     }
155:     if (islite) {DMPlexSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords);}
156:     else        {DMPlexSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords);}
157:   }
158: #else
159:   for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
160: #endif
161:   return(0);
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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

282:   MPI_Comm_rank(comm, &rank);
283:   if (rank == 0) {
284:     const PetscInt debug = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

441:         lid  = EG_indexBodyTopo(body, loop);
442:         EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);

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

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

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

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

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

573:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
574:     PetscContainerSetPointer(modelObj, model);
575:     PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
576:     PetscContainerDestroy(&modelObj);

578:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
579:     PetscContainerSetPointer(contextObj, context);
580:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
581:     PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
582:     PetscContainerDestroy(&contextObj);
583:   }
584:   /* Label points */
585:   DMCreateLabel(dm, "EGADS Body ID");
586:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
587:   DMCreateLabel(dm, "EGADS Face ID");
588:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
589:   DMCreateLabel(dm, "EGADS Edge ID");
590:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
591:   DMCreateLabel(dm, "EGADS Vertex ID");
592:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
593:   cOff = 0;
594:   for (b = 0; b < nbodies; ++b) {
595:     ego body = bodies[b];
596:     int id, Nl, l;

598:     EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
599:     for (l = 0; l < Nl; ++l) {
600:       ego  loop = lobjs[l];
601:       ego *fobjs;
602:       int  lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;

604:       lid  = EG_indexBodyTopo(body, loop);
605:       EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
606:       if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
607:       fid  = EG_indexBodyTopo(body, fobjs[0]);
608:       EG_free(fobjs);
609:       EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
610:       for (e = 0; e < Ne; ++e) {
611:         ego             edge = objs[e];
612:         int             eid, Nv, v;
613:         PetscInt        points[3], support[2], numEdges, edgeNum;
614:         const PetscInt *edges;

616:         eid  = EG_indexBodyTopo(body, edge);
617:         EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
618:         if (mtype == DEGENERATE) continue;
619:         else                     ++Ner;
620:         for (v = 0; v < Nv; ++v) {
621:           ego vertex = nobjs[v];

623:           id   = EG_indexBodyTopo(body, vertex);
624:           DMLabelSetValue(edgeLabel, numCells + id-1, eid);
625:           points[v*2] = numCells + id-1;
626:         }
627:         PetscHMapIGet(edgeMap, eid-1, &edgeNum);
628:         points[1] = numCells + numVertices - newVertices + edgeNum;

630:         DMLabelSetValue(edgeLabel, points[1], eid);
631:         support[0] = points[0];
632:         support[1] = points[1];
633:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
634:         if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
635:         DMLabelSetValue(edgeLabel, edges[0], eid);
636:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
637:         support[0] = points[1];
638:         support[1] = points[2];
639:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
640:         if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
641:         DMLabelSetValue(edgeLabel, edges[0], eid);
642:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
643:       }
644:       switch (Ner) {
645:         case 2: Nt = 2;break;
646:         case 3: Nt = 4;break;
647:         case 4: Nt = 8;break;
648:         default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
649:       }
650:       for (t = 0; t < Nt; ++t) {
651:         DMLabelSetValue(bodyLabel, cOff+t, b);
652:         DMLabelSetValue(faceLabel, cOff+t, fid);
653:       }
654:       cOff += Nt;
655:     }
656:     EG_free(lobjs);
657:   }
658:   PetscHMapIDestroy(&edgeMap);
659:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
660:   for (c = cStart; c < cEnd; ++c) {
661:     PetscInt *closure = NULL;
662:     PetscInt  clSize, cl, bval, fval;

664:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
665:     DMLabelGetValue(bodyLabel, c, &bval);
666:     DMLabelGetValue(faceLabel, c, &fval);
667:     for (cl = 0; cl < clSize*2; cl += 2) {
668:       DMLabelSetValue(bodyLabel, closure[cl], bval);
669:       DMLabelSetValue(faceLabel, closure[cl], fval);
670:     }
671:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
672:   }
673:   *newdm = dm;
674:   return(0);
675: }

677: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
678: {
679:   DMLabel         bodyLabel, faceLabel, edgeLabel, vertexLabel;
680:   // EGADS/EGADSLite variables
681:   ego             geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs;
682:   ego             topRef, prev, next;
683:   int             oclass, mtype, nbodies, *senses, *lSenses, *eSenses;
684:   int             b;
685:   // PETSc variables
686:   DM              dm;
687:   PetscHMapI      edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL;
688:   PetscInt        dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0;
689:   PetscInt        cellCntr = 0, numPoints = 0;
690:   PetscInt        *cells  = NULL;
691:   const PetscInt  *cone = NULL;
692:   PetscReal       *coords = NULL;
693:   PetscMPIInt      rank;
694:   PetscErrorCode   ierr;

697:   MPI_Comm_rank(comm, &rank);
698:   if (!rank) {
699:     // ---------------------------------------------------------------------------------------------------
700:     // Generate Petsc Plex
701:     //  Get all Nodes in model, record coordinates in a correctly formatted array
702:     //  Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
703:     //  We need to uniformly refine the initial geometry to guarantee a valid mesh

705:   // Caluculate cell and vertex sizes
706:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);

708:     PetscHMapICreate(&edgeMap);
709:   PetscHMapICreate(&bodyIndexMap);
710:   PetscHMapICreate(&bodyVertexMap);
711:   PetscHMapICreate(&bodyEdgeMap);
712:   PetscHMapICreate(&bodyEdgeGlobalMap);
713:   PetscHMapICreate(&bodyFaceMap);

715:   for (b = 0; b < nbodies; ++b) {
716:       ego             body = bodies[b];
717:     int             Nf, Ne, Nv;
718:     PetscHashIter   BIiter, BViter, BEiter, BEGiter, BFiter, EMiter;
719:     PetscBool       BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound;

721:     PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound);
722:     PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
723:     PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
724:     PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
725:     PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);

727:     if (!BIfound)  {PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices);}
728:     if (!BVfound)  {PetscHMapISet(bodyVertexMap, b, numVertices);}
729:     if (!BEfound)  {PetscHMapISet(bodyEdgeMap, b, numEdges);}
730:     if (!BEGfound) {PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr);}
731:     if (!BFfound)  {PetscHMapISet(bodyFaceMap, b, numFaces);}

733:     EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
734:     EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);
735:     EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
736:     EG_free(fobjs);
737:     EG_free(eobjs);
738:     EG_free(nobjs);

740:     // Remove DEGENERATE EDGES from Edge count
741:     EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);
742:     int Netemp = 0;
743:     for (int e = 0; e < Ne; ++e) {
744:       ego     edge = eobjs[e];
745:       int     eid;

747:       EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
748:       eid = EG_indexBodyTopo(body, edge);

750:       PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound);
751:       if (mtype == DEGENERATE) {
752:         if (!EMfound) {PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1);}
753:       }
754:       else {
755:       ++Netemp;
756:         if (!EMfound) {PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp);}
757:       }
758:     }
759:     EG_free(eobjs);

761:     // Determine Number of Cells
762:     EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
763:     for (int f = 0; f < Nf; ++f) {
764:         ego     face = fobjs[f];
765:     int     edgeTemp = 0;

767:       EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs);
768:       for (int e = 0; e < Ne; ++e) {
769:         ego     edge = eobjs[e];

771:         EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
772:         if (mtype != DEGENERATE) {++edgeTemp;}
773:       }
774:       numCells += (2 * edgeTemp);
775:       EG_free(eobjs);
776:     }
777:     EG_free(fobjs);

779:     numFaces    += Nf;
780:     numEdges    += Netemp;
781:     numVertices += Nv;
782:     edgeCntr    += Ne;
783:   }

785:   // Set up basic DMPlex parameters
786:   dim        = 2;    // Assumes 3D Models :: Need to handle 2D modles in the future
787:   cdim       = 3;     // Assumes 3D Models :: Need to update to handle 2D modles in future
788:   numCorners = 3;     // Split Faces into triangles
789:     numPoints  = numVertices + numEdges + numFaces;   // total number of coordinate points

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

793:   // Get Vertex Coordinates and Set up Cells
794:   for (b = 0; b < nbodies; ++b) {
795:     ego             body = bodies[b];
796:     int             Nf, Ne, Nv;
797:     PetscInt        bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
798:     PetscHashIter   BViter, BEiter, BEGiter, BFiter, EMiter;
799:     PetscBool       BVfound, BEfound, BEGfound, BFfound, EMfound;

801:     // Vertices on Current Body
802:     EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);

804:     PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
805:     if (!BVfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
806:     PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart);

808:     for (int v = 0; v < Nv; ++v) {
809:       ego    vertex = nobjs[v];
810:     double limits[4];
811:     int    id, dummy;

813:     EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
814:     id = EG_indexBodyTopo(body, vertex);

816:     coords[(bodyVertexIndexStart + id - 1)*cdim + 0] = limits[0];
817:     coords[(bodyVertexIndexStart + id - 1)*cdim + 1] = limits[1];
818:     coords[(bodyVertexIndexStart + id - 1)*cdim + 2] = limits[2];
819:     }
820:     EG_free(nobjs);

822:     // Edge Midpoint Vertices on Current Body
823:     EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);

825:     PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
826:     if (!BEfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
827:     PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart);

829:     PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
830:     if (!BEGfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
831:     PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart);

833:     for (int e = 0; e < Ne; ++e) {
834:       ego          edge = eobjs[e];
835:     double       range[2], avgt[1], cntrPnt[9];
836:     int          eid, eOffset;
837:     int          periodic;

839:     EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
840:     if (mtype == DEGENERATE) {continue;}

842:     eid = EG_indexBodyTopo(body, edge);

844:     // get relative offset from globalEdgeID Vector
845:     PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
846:       if (!EMfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1);
847:       PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);

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

852:     EG_evaluate(edge, avgt, cntrPnt);
853:     coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 0] = cntrPnt[0];
854:         coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 1] = cntrPnt[1];
855:     coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 2] = cntrPnt[2];
856:     }
857:     EG_free(eobjs);

859:     // Face Midpoint Vertices on Current Body
860:     EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);

862:     PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);
863:     if (!BFfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
864:     PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart);

866:     for (int f = 0; f < Nf; ++f) {
867:     ego       face = fobjs[f];
868:     double    range[4], avgUV[2], cntrPnt[18];
869:     int       peri, id;

871:     id = EG_indexBodyTopo(body, face);
872:     EG_getRange(face, range, &peri);

874:     avgUV[0] = (range[0] + range[1]) / 2.;
875:     avgUV[1] = (range[2] + range[3]) / 2.;
876:     EG_evaluate(face, avgUV, cntrPnt);

878:     coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 0] = cntrPnt[0];
879:     coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 1] = cntrPnt[1];
880:     coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 2] = cntrPnt[2];
881:     }
882:     EG_free(fobjs);

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

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

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

901:       if (Nl > 1) SETERRQ1(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);
902:     for (int l = 0; l < Nl; ++l) {
903:           ego      loop = lobjs[l];

905:           EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses);
906:       for (int e = 0; e < Ne; ++e) {
907:         ego     edge = eobjs[e];
908:         int     eid, eOffset;

910:         EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
911:       eid = EG_indexBodyTopo(body, edge);
912:         if (mtype == DEGENERATE) { continue; }

914:         // get relative offset from globalEdgeID Vector
915:         PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
916:           if (!EMfound) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
917:           PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);

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

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

923:         if (eSenses[e] > 0) { startID = EG_indexBodyTopo(body, nobjs[0]); endID = EG_indexBodyTopo(body, nobjs[1]); }
924:         else { startID = EG_indexBodyTopo(body, nobjs[1]); endID = EG_indexBodyTopo(body, nobjs[0]); }

926:       // Define 2 Cells per Edge with correct orientation
927:       cells[cellCntr*numCorners + 0] = midFaceID;
928:       cells[cellCntr*numCorners + 1] = bodyVertexIndexStart + startID - 1;
929:       cells[cellCntr*numCorners + 2] = midPntID;

931:       cells[cellCntr*numCorners + 3] = midFaceID;
932:       cells[cellCntr*numCorners + 4] = midPntID;
933:       cells[cellCntr*numCorners + 5] = bodyVertexIndexStart + endID - 1;

935:       cellCntr = cellCntr + 2;
936:       }
937:     }
938:     }
939:     EG_free(fobjs);
940:   }
941:   }

943:   // Generate DMPlex
944:   DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
945:   PetscFree2(coords, cells);
946:   PetscInfo1(dm, " Total Number of Unique Cells    = %D \n", numCells);
947:   PetscInfo1(dm, " Total Number of Unique Vertices = %D \n", numVertices);

949:   // Embed EGADS model in DM
950:   {
951:     PetscContainer modelObj, contextObj;

953:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
954:     PetscContainerSetPointer(modelObj, model);
955:     PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
956:     PetscContainerDestroy(&modelObj);

958:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
959:     PetscContainerSetPointer(contextObj, context);
960:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
961:     PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
962:     PetscContainerDestroy(&contextObj);
963:   }
964:   // Label points
965:   PetscInt   nStart, nEnd;

967:   DMCreateLabel(dm, "EGADS Body ID");
968:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
969:   DMCreateLabel(dm, "EGADS Face ID");
970:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
971:   DMCreateLabel(dm, "EGADS Edge ID");
972:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
973:   DMCreateLabel(dm, "EGADS Vertex ID");
974:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);

976:   DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd);

978:   cellCntr = 0;
979:   for (b = 0; b < nbodies; ++b) {
980:     ego             body = bodies[b];
981:   int             Nv, Ne, Nf;
982:   PetscInt        bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
983:   PetscHashIter   BViter, BEiter, BEGiter, BFiter, EMiter;
984:   PetscBool       BVfound, BEfound, BEGfound, BFfound, EMfound;

986:   PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
987:   if (!BVfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
988:   PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart);

990:   PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
991:   if (!BEfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
992:   PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart);

994:     PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);
995:   if (!BFfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
996:   PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart);

998:     PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
999:     if (!BEGfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
1000:     PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart);

1002:   EG_getBodyTopos(body, NULL, FACE,  &Nf, &fobjs);
1003:   for (int f = 0; f < Nf; ++f) {
1004:     ego   face = fobjs[f];
1005:       int   fID, Nl;

1007:     fID  = EG_indexBodyTopo(body, face);

1009:     EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs);
1010:     for (int l = 0; l < Nl; ++l) {
1011:         ego  loop = lobjs[l];
1012:     int  lid;

1014:     lid  = EG_indexBodyTopo(body, loop);
1015:       if (Nl > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);

1017:     EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses);
1018:     for (int e = 0; e < Ne; ++e) {
1019:       ego     edge = eobjs[e];
1020:       int     eid, eOffset;

1022:       // Skip DEGENERATE Edges
1023:       EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
1024:       if (mtype == DEGENERATE) {continue;}
1025:       eid = EG_indexBodyTopo(body, edge);

1027:       // get relative offset from globalEdgeID Vector
1028:       PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
1029:       if (!EMfound) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
1030:       PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);

1032:       EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs);
1033:       for (int v = 0; v < Nv; ++v){
1034:         ego vertex = nobjs[v];
1035:         int vID;

1037:         vID = EG_indexBodyTopo(body, vertex);
1038:         DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b);
1039:         DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID);
1040:       }
1041:       EG_free(nobjs);

1043:       DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b);
1044:       DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid);

1046:       // Define Cell faces
1047:       for (int jj = 0; jj < 2; ++jj){
1048:         DMLabelSetValue(bodyLabel, cellCntr, b);
1049:         DMLabelSetValue(faceLabel, cellCntr, fID);
1050:         DMPlexGetCone(dm, cellCntr, &cone);

1052:         DMLabelSetValue(bodyLabel, cone[0], b);
1053:         DMLabelSetValue(faceLabel, cone[0], fID);

1055:         DMLabelSetValue(bodyLabel, cone[1], b);
1056:         DMLabelSetValue(edgeLabel, cone[1], eid);

1058:        DMLabelSetValue(bodyLabel, cone[2], b);
1059:        DMLabelSetValue(faceLabel, cone[2], fID);

1061:        cellCntr = cellCntr + 1;
1062:       }
1063:     }
1064:     }
1065:     EG_free(lobjs);

1067:     DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b);
1068:     DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID);
1069:   }
1070:   EG_free(fobjs);
1071:   }

1073:   PetscHMapIDestroy(&edgeMap);
1074:   PetscHMapIDestroy(&bodyIndexMap);
1075:   PetscHMapIDestroy(&bodyVertexMap);
1076:   PetscHMapIDestroy(&bodyEdgeMap);
1077:   PetscHMapIDestroy(&bodyEdgeGlobalMap);
1078:   PetscHMapIDestroy(&bodyFaceMap);

1080:   *newdm = dm;
1081:   return(0);
1082: }

1084: static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
1085: {
1086:   DMLabel              bodyLabel, faceLabel, edgeLabel, vertexLabel;
1087:   /* EGADSLite variables */
1088:   ego                  geom, *bodies, *fobjs;
1089:   int                  b, oclass, mtype, nbodies, *senses;
1090:   int                  totalNumTris = 0, totalNumPoints = 0;
1091:   double               boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize;
1092:   /* PETSc variables */
1093:   DM                   dm;
1094:   PetscHMapI           pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL;
1095:   PetscHMapI           pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL;
1096:   PetscInt             dim = -1, cdim = -1, numCorners = 0, counter = 0;
1097:   PetscInt            *cells  = NULL;
1098:   const PetscInt      *cone = NULL;
1099:   PetscReal           *coords = NULL;
1100:   PetscMPIInt          rank;
1101:   PetscErrorCode       ierr;

1104:   MPI_Comm_rank(comm, &rank);
1105:   if (!rank) {
1106:     // ---------------------------------------------------------------------------------------------------
1107:     // Generate Petsc Plex from EGADSlite created Tessellation of geometry
1108:     // ---------------------------------------------------------------------------------------------------

1110:   // Caluculate cell and vertex sizes
1111:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);

1113:   PetscHMapICreate(&pointIndexStartMap);
1114:   PetscHMapICreate(&triIndexStartMap);
1115:   PetscHMapICreate(&pTypeLabelMap);
1116:   PetscHMapICreate(&pIndexLabelMap);
1117:   PetscHMapICreate(&pBodyIndexLabelMap);
1118:   PetscHMapICreate(&triFaceIDLabelMap);
1119:   PetscHMapICreate(&triBodyIDLabelMap);

1121:   /* Create Tessellation of Bodies */
1122:   ego tessArray[nbodies];

1124:   for (b = 0; b < nbodies; ++b) {
1125:     ego             body = bodies[b];
1126:     double          params[3] = {0.0, 0.0, 0.0};    // Parameters for Tessellation
1127:     int             Nf, bodyNumPoints = 0, bodyNumTris = 0;
1128:     PetscHashIter   PISiter, TISiter;
1129:     PetscBool       PISfound, TISfound;

1131:     /* Store Start Index for each Body's Point and Tris */
1132:     PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound);
1133:     PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound);

1135:     if (!PISfound)  {PetscHMapISet(pointIndexStartMap, b, totalNumPoints);}
1136:     if (!TISfound)  {PetscHMapISet(triIndexStartMap, b, totalNumTris);}

1138:     /* Calculate Tessellation parameters based on Bounding Box */
1139:     /* Get Bounding Box Dimensions of the BODY */
1140:     EG_getBoundingBox(body, boundBox);
1141:     tessSize = boundBox[3] - boundBox[0];
1142:     if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1];
1143:     if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2];

1145:     // TODO :: May want to give users tessellation parameter options //
1146:     params[0] = 0.0250 * tessSize;
1147:     params[1] = 0.0075 * tessSize;
1148:     params[2] = 15.0;

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

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

1154:     for (int f = 0; f < Nf; ++f) {
1155:       ego             face = fobjs[f];
1156:     int             len, fID, ntris;
1157:     const int      *ptype, *pindex, *ptris, *ptric;
1158:     const double   *pxyz, *puv;

1160:     // Get Face ID //
1161:     fID = EG_indexBodyTopo(body, face);

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

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

1169:     // Check out the point index and coordinate //
1170:     for (int p = 0; p < len; ++p) {
1171:         int global;

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

1175:       // Determine the total number of points in the tessellation //
1176:         bodyNumPoints = PetscMax(bodyNumPoints, global);
1177:     }
1178:     }
1179:     EG_free(fobjs);

1181:     totalNumPoints += bodyNumPoints;
1182:     totalNumTris += bodyNumTris;
1183:     }
1184:   //}  - Original End of (!rank)

1186:   dim = 2;
1187:   cdim = 3;
1188:   numCorners = 3;
1189:   //PetscInt counter = 0;

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

1195:   for (b = 0; b < nbodies; ++b) {
1196:   ego             body = bodies[b];
1197:   int             Nf;
1198:   PetscInt        pointIndexStart;
1199:   PetscHashIter   PISiter;
1200:   PetscBool       PISfound;

1202:   PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound);
1203:   if (!PISfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b);
1204:   PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart);

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

1208:   for (int f = 0; f < Nf; ++f) {
1209:     /* Get Face Object */
1210:     ego              face = fobjs[f];
1211:     int              len, fID, ntris;
1212:     const int       *ptype, *pindex, *ptris, *ptric;
1213:     const double    *pxyz, *puv;

1215:     /* Get Face ID */
1216:     fID = EG_indexBodyTopo(body, face);

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

1221:     /* Check out the point index and coordinate */
1222:     for (int p = 0; p < len; ++p) {
1223:     int              global;
1224:     PetscHashIter    PTLiter, PILiter, PBLiter;
1225:     PetscBool        PTLfound, PILfound, PBLfound;

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

1229:     /* Set the coordinates array for DAG */
1230:     coords[((global-1+pointIndexStart)*3) + 0] = pxyz[(p*3)+0];
1231:     coords[((global-1+pointIndexStart)*3) + 1] = pxyz[(p*3)+1];
1232:     coords[((global-1+pointIndexStart)*3) + 2] = pxyz[(p*3)+2];

1234:     /* Store Geometry Label Information for DMLabel assignment later */
1235:     PetscHMapIFind(pTypeLabelMap, global-1+pointIndexStart, &PTLiter, &PTLfound);
1236:         PetscHMapIFind(pIndexLabelMap, global-1+pointIndexStart, &PILiter, &PILfound);
1237:         PetscHMapIFind(pBodyIndexLabelMap, global-1+pointIndexStart, &PBLiter, &PBLfound);

1239:         if (!PTLfound)  {PetscHMapISet(pTypeLabelMap, global-1+pointIndexStart, ptype[p]);}
1240:         if (!PILfound)  {PetscHMapISet(pIndexLabelMap, global-1+pointIndexStart, pindex[p]);}
1241:         if (!PBLfound)  {PetscHMapISet(pBodyIndexLabelMap, global-1+pointIndexStart, b);}

1243:     if (ptype[p] < 0) { PetscHMapISet(pIndexLabelMap, global-1+pointIndexStart, fID);}
1244:     }

1246:     for (int t = 0; t < (int) ntris; ++t){
1247:     int             global, globalA, globalB;
1248:     PetscHashIter   TFLiter, TBLiter;
1249:       PetscBool       TFLfound, TBLfound;

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

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

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

1260:     PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound);
1261:         PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound);

1263:     if (!TFLfound)  {PetscHMapISet(triFaceIDLabelMap, counter, fID);}
1264:         if (!TBLfound)  {PetscHMapISet(triBodyIDLabelMap, counter, b);}

1266:     counter += 1;
1267:     }
1268:   }
1269:   EG_free(fobjs);
1270:   }
1271: }

1273:   //Build DMPlex
1274:   DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
1275:   PetscFree2(coords, cells);

1277:   // Embed EGADS model in DM
1278:   {
1279:     PetscContainer modelObj, contextObj;

1281:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
1282:     PetscContainerSetPointer(modelObj, model);
1283:     PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
1284:     PetscContainerDestroy(&modelObj);

1286:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
1287:     PetscContainerSetPointer(contextObj, context);
1288:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
1289:     PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
1290:     PetscContainerDestroy(&contextObj);
1291:   }

1293:   // Label Points
1294:   DMCreateLabel(dm, "EGADS Body ID");
1295:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
1296:   DMCreateLabel(dm, "EGADS Face ID");
1297:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
1298:   DMCreateLabel(dm, "EGADS Edge ID");
1299:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
1300:   DMCreateLabel(dm, "EGADS Vertex ID");
1301:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);

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

1306:   DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd);
1307:   DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd);
1308:   DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd);

1310:   /* Set DMLabels for NODES */
1311:   for (int n = nStart; n < nEnd; ++n) {
1312:     int             pTypeVal, pIndexVal, pBodyVal;
1313:     PetscHashIter   PTLiter, PILiter, PBLiter;
1314:     PetscBool       PTLfound, PILfound, PBLfound;

1316:     //Converted to Hash Tables
1317:     PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound);
1318:     if (!PTLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n);
1319:     PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal);

1321:     PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound);
1322:     if (!PILfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n);
1323:     PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal);

1325:     PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound);
1326:     if (!PBLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n);
1327:     PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal);

1329:     DMLabelSetValue(bodyLabel, n, pBodyVal);
1330:     if (pTypeVal == 0) {DMLabelSetValue(vertexLabel, n, pIndexVal);}
1331:     if (pTypeVal >  0) {DMLabelSetValue(edgeLabel, n, pIndexVal);}
1332:     if (pTypeVal <  0) {DMLabelSetValue(faceLabel, n, pIndexVal);}
1333:   }

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

1339:   DMPlexGetCone(dm, e, &cone);
1340:   DMLabelGetValue(bodyLabel, cone[0], &bodyID_0);    // Do I need to check the other end?
1341:   DMLabelGetValue(vertexLabel, cone[0], &vertexID_0);
1342:   DMLabelGetValue(vertexLabel, cone[1], &vertexID_1);
1343:   DMLabelGetValue(edgeLabel, cone[0], &edgeID_0);
1344:   DMLabelGetValue(edgeLabel, cone[1], &edgeID_1);
1345:   DMLabelGetValue(faceLabel, cone[0], &faceID_0);
1346:   DMLabelGetValue(faceLabel, cone[1], &faceID_1);

1348:   DMLabelSetValue(bodyLabel, e, bodyID_0);

1350:   if (edgeID_0 == edgeID_1) { DMLabelSetValue(edgeLabel, e, edgeID_0); }
1351:   else if (vertexID_0 > 0 && edgeID_1 > 0) { DMLabelSetValue(edgeLabel, e, edgeID_1); }
1352:   else if (vertexID_1 > 0 && edgeID_0 > 0) { DMLabelSetValue(edgeLabel, e, edgeID_0); }
1353:   else { /* Do Nothing */ }
1354:   }

1356:   /* Set DMLabels for Cells */
1357:   for (int f = fStart; f < fEnd; ++f){
1358:   int             edgeID_0;
1359:   PetscInt        triBodyVal, triFaceVal;
1360:   PetscHashIter   TFLiter, TBLiter;
1361:   PetscBool       TFLfound, TBLfound;

1363:     // Convert to Hash Table
1364:   PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound);
1365:   if (!TFLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f);
1366:   PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal);

1368:   PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound);
1369:   if (!TBLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f);
1370:     PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal);

1372:   DMLabelSetValue(bodyLabel, f, triBodyVal);
1373:   DMLabelSetValue(faceLabel, f, triFaceVal);

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

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

1381:     if (edgeID_0 < 0) {
1382:     DMLabelSetValue(bodyLabel, cone[jj], triBodyVal);
1383:       DMLabelSetValue(faceLabel, cone[jj], triFaceVal);
1384:     }
1385:   }
1386:   }

1388:   *newdm = dm;
1389:   return(0);
1390: }
1391: #endif

1393: /*@
1394:   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.

1396:   Collective on dm

1398:   Input Parameter:
1399: . dm - The uninflated DM object representing the mesh

1401:   Output Parameter:
1402: . dm - The inflated DM object representing the mesh

1404:   Level: intermediate

1406: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS()
1407: @*/
1408: PetscErrorCode DMPlexInflateToGeomModel(DM dm)
1409: {
1410: #if defined(PETSC_HAVE_EGADS)
1411:   /* EGADS Variables */
1412:   ego            model, geom, body, face, edge;
1413:   ego           *bodies;
1414:   int            Nb, oclass, mtype, *senses;
1415:   double         result[3];
1416:   /* PETSc Variables */
1417:   DM             cdm;
1418:   PetscContainer modelObj;
1419:   DMLabel        bodyLabel, faceLabel, edgeLabel, vertexLabel;
1420:   Vec            coordinates;
1421:   PetscScalar   *coords;
1422:   PetscInt       bodyID, faceID, edgeID, vertexID;
1423:   PetscInt       cdim, d, vStart, vEnd, v;
1425: #endif

1428: #if defined(PETSC_HAVE_EGADS)
1429:   PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
1430:   if (!modelObj) return(0);
1431:   DMGetCoordinateDim(dm, &cdim);
1432:   DMGetCoordinateDM(dm, &cdm);
1433:   DMGetCoordinatesLocal(dm, &coordinates);
1434:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
1435:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
1436:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
1437:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);

1439:   PetscContainerGetPointer(modelObj, (void **) &model);
1440:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);

1442:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1443:   VecGetArrayWrite(coordinates, &coords);
1444:   for (v = vStart; v < vEnd; ++v) {
1445:     PetscScalar *vcoords;

1447:     DMLabelGetValue(bodyLabel, v, &bodyID);
1448:     DMLabelGetValue(faceLabel, v, &faceID);
1449:     DMLabelGetValue(edgeLabel, v, &edgeID);
1450:     DMLabelGetValue(vertexLabel, v, &vertexID);

1452:     if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
1453:     body = bodies[bodyID];

1455:     DMPlexPointLocalRef(cdm, v, coords, (void *) &vcoords);
1456:     if (edgeID > 0) {
1457:       /* Snap to EDGE at nearest location */
1458:       double params[1];
1459:       EG_objectBodyTopo(body, EDGE, edgeID, &edge);
1460:       EG_invEvaluate(edge, vcoords, params, result); // Get (x,y,z) of nearest point on EDGE
1461:       for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1462:     } else if (faceID > 0) {
1463:       /* Snap to FACE at nearest location */
1464:       double params[2];
1465:       EG_objectBodyTopo(body, FACE, faceID, &face);
1466:       EG_invEvaluate(face, vcoords, params, result); // Get (x,y,z) of nearest point on FACE
1467:       for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1468:     }
1469:   }
1470:   VecRestoreArrayWrite(coordinates, &coords);
1471:   /* Clear out global coordinates */
1472:   VecDestroy(&dm->coordinates);
1473: #endif
1474:   return(0);
1475: }

1477: /*@C
1478:   DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADS, IGES, or STEP file.

1480:   Collective

1482:   Input Parameters:
1483: + comm     - The MPI communicator
1484: - filename - The name of the EGADS, IGES, or STEP file

1486:   Output Parameter:
1487: . dm       - The DM object representing the mesh

1489:   Level: beginner

1491: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS(), DMPlexCreateEGADSLiteFromFile()
1492: @*/
1493: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
1494: {
1495:   PetscMPIInt    rank;
1496: #if defined(PETSC_HAVE_EGADS)
1497:   ego            context= NULL, model = NULL;
1498: #endif
1499:   PetscBool      printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE;

1504:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
1505:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL);
1506:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL);
1507:   MPI_Comm_rank(comm, &rank);
1508: #if defined(PETSC_HAVE_EGADS)
1509:   if (rank == 0) {

1511:     EG_open(&context);
1512:     EG_loadModel(context, 0, filename, &model);
1513:     if (printModel) {DMPlexEGADSPrintModel_Internal(model);}

1515:   }
1516:   if (tessModel)     {DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm);}
1517:   else if (newModel) {DMPlexCreateEGADS(comm, context, model, dm);}
1518:   else               {DMPlexCreateEGADS_Internal(comm, context, model, dm);}
1519:   return(0);
1520: #else
1521:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads");
1522: #endif
1523: }