Actual source code: plexegads.c

petsc-main 2021-04-20
Report Typos and Errors
  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: */


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

 25:   Not collective

 27:   Input Parameters:
 28: + dm      - The DMPlex object
 29: . p       - The mesh point
 30: - mcoords - A coordinate point lying on the mesh point

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

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

 37:   Level: intermediate

 39: .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform()
 40: @*/
 41: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[])
 42: {
 43: #ifdef PETSC_HAVE_EGADS
 44:   DM             cdm;
 45:   DMLabel        bodyLabel, faceLabel, edgeLabel;
 46:   PetscContainer modelObj;
 47:   PetscInt       bodyID, faceID, edgeID;
 48:   ego           *bodies;
 49:   ego            model, geom, body, obj;
 50:   /* result has to hold derviatives, along with the value */
 51:   double         params[3], result[18], paramsV[16*3], resultV[16*3];
 52:   int            Nb, oclass, mtype, *senses;
 53:   Vec            coordinatesLocal;
 54:   PetscScalar   *coords = NULL;
 55:   PetscInt       Nv, v, Np = 0, pm;
 56: #endif
 57:   PetscInt       dE, d;

 61:   DMGetCoordinateDim(dm, &dE);
 62: #ifdef PETSC_HAVE_EGADS
 63:   DMGetLabel(dm, "EGADS Body ID",   &bodyLabel);
 64:   DMGetLabel(dm, "EGADS Face ID",   &faceLabel);
 65:   DMGetLabel(dm, "EGADS Edge ID",   &edgeLabel);
 66:   if (!bodyLabel || !faceLabel || !edgeLabel) {
 67:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 68:     return(0);
 69:   }
 70:   DMGetCoordinateDM(dm, &cdm);
 71:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
 72:   DMLabelGetValue(bodyLabel,   p, &bodyID);
 73:   DMLabelGetValue(faceLabel,   p, &faceID);
 74:   DMLabelGetValue(edgeLabel,   p, &edgeID);
 75:   PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
 76:   PetscContainerGetPointer(modelObj, (void **) &model);
 77:   EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 78:   if (bodyID < 0 || bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
 79:   body = bodies[bodyID];

 81:   if (edgeID >= 0)      {EG_objectBodyTopo(body, EDGE, edgeID, &obj); Np = 1;}
 82:   else if (faceID >= 0) {EG_objectBodyTopo(body, FACE, faceID, &obj); Np = 2;}
 83:   else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D is not in edge or face label for EGADS", p);
 84:   /* Calculate parameters (t or u,v) for vertices */
 85:   DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 86:   Nv  /= dE;
 87:   if (Nv == 1) {
 88:     DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 89:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 90:     return(0);
 91:   }
 92:   if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);
 93:   for (v = 0; v < Nv; ++v) {EG_invEvaluate(obj, &coords[v*dE], &paramsV[v*3], &resultV[v*3]);}
 94:   DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 95:   /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
 96:   for (pm = 0; pm < Np; ++pm) {
 97:     params[pm] = 0.;
 98:     for (v = 0; v < Nv; ++v) {params[pm] += paramsV[v*3+pm];}
 99:     params[pm] /= Nv;
100:   }
101:   /* TODO Check
102:     double range[4]; // [umin, umax, vmin, vmax]
103:     int    peri;
104:     EG_getRange(face, range, &peri);
105:     if ((paramsNew[0] < range[0]) || (paramsNew[0] > range[1]) || (paramsNew[1] < range[2]) || (paramsNew[1] > range[3])) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE);
106:   */
107:   /* Put coordinates for new vertex in result[] */
108:   EG_evaluate(obj, params, result);
109:   for (d = 0; d < dE; ++d) gcoords[d] = result[d];
110: #else
111:   for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
112: #endif
113:   return(0);
114: }

116: #if defined(PETSC_HAVE_EGADS)
117: static PetscErrorCode DMPlexEGADSPrintModel(ego model)
118: {
119:   ego            geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
120:   int            oclass, mtype, *senses;
121:   int            Nb, b;

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

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

133:     /* Output Basic Model Topology */
134:     EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
135:     PetscPrintf(PETSC_COMM_SELF, "   Number of SHELLS: %d \n", Nsh);
136:     EG_free(objs);

138:     EG_getBodyTopos(body, NULL, FACE,  &Nf, &objs);
139:     PetscPrintf(PETSC_COMM_SELF, "   Number of FACES: %d \n", Nf);
140:     EG_free(objs);

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

145:     EG_getBodyTopos(body, NULL, EDGE,  &Ne, &objs);
146:     PetscPrintf(PETSC_COMM_SELF, "   Number of EDGES: %d \n", Ne);
147:     EG_free(objs);

149:     EG_getBodyTopos(body, NULL, NODE,  &Nv, &objs);
150:     PetscPrintf(PETSC_COMM_SELF, "   Number of NODES: %d \n", Nv);
151:     EG_free(objs);

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

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

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

162:       for (e = 0; e < Ne; ++e) {
163:         ego    edge      = objs[e];
164:         double range[4]  = {0., 0., 0., 0.};
165:         double point[3]  = {0., 0., 0.};
166:         double params[3] = {0., 0., 0.};
167:         double result[18];
168:         int    peri;

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

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

176:         /* Get NODE info which associated with the current EDGE */
177:         EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
178:         if (mtype == DEGENERATE) {
179:           PetscPrintf(PETSC_COMM_SELF, "  EDGE %d is DEGENERATE \n", id);
180:         } else {
181:           params[0] = range[0];
182:           EG_evaluate(edge, params, result);
183:           PetscPrintf(PETSC_COMM_SELF, "   between (%lf, %lf, %lf)", result[0], result[1], result[2]);
184:           params[0] = range[1];
185:           EG_evaluate(edge, params, result);
186:           PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]);
187:         }

189:         for (v = 0; v < Nv; ++v) {
190:           ego    vertex = nobjs[v];
191:           double limits[4];
192:           int    dummy;

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

199:           point[0] = point[0] + limits[0];
200:           point[1] = point[1] + limits[1];
201:           point[2] = point[2] + limits[2];
202:         }
203:       }
204:     }
205:     EG_free(lobjs);
206:   }
207:   return(0);
208: }

210: static PetscErrorCode DMPlexEGADSDestroy_Private(void *context)
211: {
212:   if (context) EG_close((ego) context);
213:   return 0;
214: }

216: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
217: {
218:   DMLabel        bodyLabel, faceLabel, edgeLabel;
219:   PetscInt       cStart, cEnd, c;
220:   /* EGADSLite variables */
221:   ego            geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
222:   int            oclass, mtype, nbodies, *senses;
223:   int            b;
224:   /* PETSc variables */
225:   DM             dm;
226:   PetscHMapI     edgeMap = NULL;
227:   PetscInt       dim = -1, cdim = -1, numCorners = 0, maxCorners = 0, numVertices = 0, newVertices = 0, numEdges = 0, numCells = 0, newCells = 0, numQuads = 0, cOff = 0, fOff = 0;
228:   PetscInt      *cells  = NULL, *cone = NULL;
229:   PetscReal     *coords = NULL;
230:   PetscMPIInt    rank;

234:   MPI_Comm_rank(comm, &rank);
235:   if (!rank) {
236:     const PetscInt debug = 0;

238:     /* ---------------------------------------------------------------------------------------------------
239:     Generate Petsc Plex
240:       Get all Nodes in model, record coordinates in a correctly formatted array
241:       Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
242:       We need to uniformly refine the initial geometry to guarantee a valid mesh
243:     */

245:     /* Calculate cell and vertex sizes */
246:     EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
247:     PetscHMapICreate(&edgeMap);
248:     numEdges = 0;
249:     for (b = 0; b < nbodies; ++b) {
250:       ego body = bodies[b];
251:       int id, Nl, l, Nv, v;

253:       EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
254:       for (l = 0; l < Nl; ++l) {
255:         ego loop = lobjs[l];
256:         int Ner  = 0, Ne, e, Nc;

258:         EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
259:         for (e = 0; e < Ne; ++e) {
260:           ego edge = objs[e];
261:           int Nv, id;
262:           PetscHashIter iter;
263:           PetscBool     found;

265:           EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
266:           if (mtype == DEGENERATE) continue;
267:           id   = EG_indexBodyTopo(body, edge);
268:           PetscHMapIFind(edgeMap, id-1, &iter, &found);
269:           if (!found) {PetscHMapISet(edgeMap, id-1, numEdges++);}
270:           ++Ner;
271:         }
272:         if (Ner == 2)      {Nc = 2;}
273:         else if (Ner == 3) {Nc = 4;}
274:         else if (Ner == 4) {Nc = 8; ++numQuads;}
275:         else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
276:         numCells += Nc;
277:         newCells += Nc-1;
278:         maxCorners = PetscMax(Ner*2+1, maxCorners);
279:       }
280:       EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
281:       for (v = 0; v < Nv; ++v) {
282:         ego vertex = nobjs[v];

284:         id = EG_indexBodyTopo(body, vertex);
285:         /* TODO: Instead of assuming contiguous ids, we could use a hash table */
286:         numVertices = PetscMax(id, numVertices);
287:       }
288:       EG_free(lobjs);
289:       EG_free(nobjs);
290:     }
291:     PetscHMapIGetSize(edgeMap, &numEdges);
292:     newVertices  = numEdges + numQuads;
293:     numVertices += newVertices;

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

300:     /* Get vertex coordinates */
301:     for (b = 0; b < nbodies; ++b) {
302:       ego body = bodies[b];
303:       int id, Nv, v;

305:       EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
306:       for (v = 0; v < Nv; ++v) {
307:         ego    vertex = nobjs[v];
308:         double limits[4];
309:         int    dummy;

311:         EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
312:         id   = EG_indexBodyTopo(body, vertex);
313:         coords[(id-1)*cdim+0] = limits[0];
314:         coords[(id-1)*cdim+1] = limits[1];
315:         coords[(id-1)*cdim+2] = limits[2];
316:       }
317:       EG_free(nobjs);
318:     }
319:     PetscHMapIClear(edgeMap);
320:     fOff     = numVertices - newVertices + numEdges;
321:     numEdges = 0;
322:     numQuads = 0;
323:     for (b = 0; b < nbodies; ++b) {
324:       ego body = bodies[b];
325:       int Nl, l;

327:       EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
328:       for (l = 0; l < Nl; ++l) {
329:         ego loop = lobjs[l];
330:         int lid, Ner = 0, Ne, e;

332:         lid  = EG_indexBodyTopo(body, loop);
333:         EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
334:         for (e = 0; e < Ne; ++e) {
335:           ego       edge = objs[e];
336:           int       eid, Nv;
337:           PetscHashIter iter;
338:           PetscBool     found;

340:           EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
341:           if (mtype == DEGENERATE) continue;
342:           ++Ner;
343:           eid  = EG_indexBodyTopo(body, edge);
344:           PetscHMapIFind(edgeMap, eid-1, &iter, &found);
345:           if (!found) {
346:             PetscInt v = numVertices - newVertices + numEdges;
347:             double range[4], params[3] = {0., 0., 0.}, result[18];
348:             int    periodic[2];

350:             PetscHMapISet(edgeMap, eid-1, numEdges++);
351:             EG_getRange(edge, range, periodic);
352:             params[0] = 0.5*(range[0] + range[1]);
353:             EG_evaluate(edge, params, result);
354:             coords[v*cdim+0] = result[0];
355:             coords[v*cdim+1] = result[1];
356:             coords[v*cdim+2] = result[2];
357:           }
358:         }
359:         if (Ner == 4) {
360:           PetscInt v = fOff + numQuads++;
361:           ego     *fobjs, face;
362:           double   range[4], params[3] = {0., 0., 0.}, result[18];
363:           int      Nf, fid, periodic[2];

365:           EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
366:           face = fobjs[0];
367:           fid  = EG_indexBodyTopo(body, face);
368:           if (Nf != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid);
369:           EG_getRange(face, range, periodic);
370:           params[0] = 0.5*(range[0] + range[1]);
371:           params[1] = 0.5*(range[2] + range[3]);
372:           EG_evaluate(face, params, result);
373:           coords[v*cdim+0] = result[0];
374:           coords[v*cdim+1] = result[1];
375:           coords[v*cdim+2] = result[2];
376:         }
377:       }
378:     }
379:     if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices);
380:     CHKMEMQ;

382:     /* Get cell vertices by traversing loops */
383:     numQuads = 0;
384:     cOff     = 0;
385:     for (b = 0; b < nbodies; ++b) {
386:       ego body = bodies[b];
387:       int id, Nl, l;

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

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

397:         for (e = 0; e < Ne; ++e) {
398:           ego edge = objs[e];
399:           int points[3];
400:           int eid, Nv, v, tmp;

402:           eid  = EG_indexBodyTopo(body, edge);
403:           EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
404:           if (mtype == DEGENERATE) continue;
405:           else                     ++Ner;
406:           if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);

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

411:             id = EG_indexBodyTopo(body, vertex);
412:             points[v*2] = id-1;
413:           }
414:           {
415:             PetscInt edgeNum;

417:             PetscHMapIGet(edgeMap, eid-1, &edgeNum);
418:             points[1] = numVertices - newVertices + edgeNum;
419:           }
420:           /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
421:           if (!nc) {
422:             for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v];
423:           } else {
424:             if (cone[nc-1] == points[0])      {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
425:             else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
426:             else if (cone[nc-3] == points[0]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
427:             else if (cone[nc-3] == points[2]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
428:             else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
429:           }
430:         }
431:         if (nc != 2*Ner)     SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner);
432:         if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;}
433:         if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners);
434:         /* Triangulate the loop */
435:         switch (Ner) {
436:           case 2: /* Bi-Segment -> 2 triangles */
437:             Nt = 2;
438:             cells[cOff*numCorners+0] = cone[0];
439:             cells[cOff*numCorners+1] = cone[1];
440:             cells[cOff*numCorners+2] = cone[2];
441:             ++cOff;
442:             cells[cOff*numCorners+0] = cone[0];
443:             cells[cOff*numCorners+1] = cone[2];
444:             cells[cOff*numCorners+2] = cone[3];
445:             ++cOff;
446:             break;
447:           case 3: /* Triangle   -> 4 triangles */
448:             Nt = 4;
449:             cells[cOff*numCorners+0] = cone[0];
450:             cells[cOff*numCorners+1] = cone[1];
451:             cells[cOff*numCorners+2] = cone[5];
452:             ++cOff;
453:             cells[cOff*numCorners+0] = cone[1];
454:             cells[cOff*numCorners+1] = cone[2];
455:             cells[cOff*numCorners+2] = cone[3];
456:             ++cOff;
457:             cells[cOff*numCorners+0] = cone[5];
458:             cells[cOff*numCorners+1] = cone[3];
459:             cells[cOff*numCorners+2] = cone[4];
460:             ++cOff;
461:             cells[cOff*numCorners+0] = cone[1];
462:             cells[cOff*numCorners+1] = cone[3];
463:             cells[cOff*numCorners+2] = cone[5];
464:             ++cOff;
465:             break;
466:           case 4: /* Quad       -> 8 triangles */
467:             Nt = 8;
468:             cells[cOff*numCorners+0] = cone[0];
469:             cells[cOff*numCorners+1] = cone[1];
470:             cells[cOff*numCorners+2] = cone[7];
471:             ++cOff;
472:             cells[cOff*numCorners+0] = cone[1];
473:             cells[cOff*numCorners+1] = cone[2];
474:             cells[cOff*numCorners+2] = cone[3];
475:             ++cOff;
476:             cells[cOff*numCorners+0] = cone[3];
477:             cells[cOff*numCorners+1] = cone[4];
478:             cells[cOff*numCorners+2] = cone[5];
479:             ++cOff;
480:             cells[cOff*numCorners+0] = cone[5];
481:             cells[cOff*numCorners+1] = cone[6];
482:             cells[cOff*numCorners+2] = cone[7];
483:             ++cOff;
484:             cells[cOff*numCorners+0] = cone[8];
485:             cells[cOff*numCorners+1] = cone[1];
486:             cells[cOff*numCorners+2] = cone[3];
487:             ++cOff;
488:             cells[cOff*numCorners+0] = cone[8];
489:             cells[cOff*numCorners+1] = cone[3];
490:             cells[cOff*numCorners+2] = cone[5];
491:             ++cOff;
492:             cells[cOff*numCorners+0] = cone[8];
493:             cells[cOff*numCorners+1] = cone[5];
494:             cells[cOff*numCorners+2] = cone[7];
495:             ++cOff;
496:             cells[cOff*numCorners+0] = cone[8];
497:             cells[cOff*numCorners+1] = cone[7];
498:             cells[cOff*numCorners+2] = cone[1];
499:             ++cOff;
500:             break;
501:           default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
502:         }
503:         if (debug) {
504:           for (t = 0; t < Nt; ++t) {
505:             PetscPrintf(PETSC_COMM_SELF, "  LOOP Corner NODEs Triangle %D (", t);
506:             for (c = 0; c < numCorners; ++c) {
507:               if (c > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
508:               PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);
509:             }
510:             PetscPrintf(PETSC_COMM_SELF, ")\n");
511:           }
512:         }
513:       }
514:       EG_free(lobjs);
515:     }
516:   }
517:   if (cOff != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells);
518:   DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
519:   PetscFree3(coords, cells, cone);
520:   PetscInfo2(dm, " Total Number of Unique Cells    = %D (%D)\n", numCells, newCells);
521:   PetscInfo2(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices);
522:   /* Embed EGADS model in DM */
523:   {
524:     PetscContainer modelObj, contextObj;

526:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
527:     PetscContainerSetPointer(modelObj, model);
528:     PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
529:     PetscContainerDestroy(&modelObj);

531:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
532:     PetscContainerSetPointer(contextObj, context);
533:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
534:     PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
535:     PetscContainerDestroy(&contextObj);
536:   }
537:   /* Label points */
538:   DMCreateLabel(dm, "EGADS Body ID");
539:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
540:   DMCreateLabel(dm, "EGADS Face ID");
541:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
542:   DMCreateLabel(dm, "EGADS Edge ID");
543:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
544:   cOff = 0;
545:   for (b = 0; b < nbodies; ++b) {
546:     ego body = bodies[b];
547:     int id, Nl, l;

549:     EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
550:     for (l = 0; l < Nl; ++l) {
551:       ego  loop = lobjs[l];
552:       ego *fobjs;
553:       int  lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;

555:       lid  = EG_indexBodyTopo(body, loop);
556:       EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
557:       if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
558:       fid  = EG_indexBodyTopo(body, fobjs[0]);
559:       EG_free(fobjs);
560:       EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
561:       for (e = 0; e < Ne; ++e) {
562:         ego             edge = objs[e];
563:         int             eid, Nv, v;
564:         PetscInt        points[3], support[2], numEdges, edgeNum;
565:         const PetscInt *edges;

567:         eid  = EG_indexBodyTopo(body, edge);
568:         EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
569:         if (mtype == DEGENERATE) continue;
570:         else                     ++Ner;
571:         for (v = 0; v < Nv; ++v) {
572:           ego vertex = nobjs[v];

574:           id   = EG_indexBodyTopo(body, vertex);
575:           DMLabelSetValue(edgeLabel, numCells + id-1, eid);
576:           points[v*2] = numCells + id-1;
577:         }
578:         PetscHMapIGet(edgeMap, eid-1, &edgeNum);
579:         points[1] = numCells + numVertices - newVertices + edgeNum;

581:         DMLabelSetValue(edgeLabel, points[1], eid);
582:         support[0] = points[0];
583:         support[1] = points[1];
584:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
585:         if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
586:         DMLabelSetValue(edgeLabel, edges[0], eid);
587:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
588:         support[0] = points[1];
589:         support[1] = points[2];
590:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
591:         if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
592:         DMLabelSetValue(edgeLabel, edges[0], eid);
593:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
594:       }
595:       switch (Ner) {
596:         case 2: Nt = 2;break;
597:         case 3: Nt = 4;break;
598:         case 4: Nt = 8;break;
599:         default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
600:       }
601:       for (t = 0; t < Nt; ++t) {
602:         DMLabelSetValue(bodyLabel, cOff+t, b);
603:         DMLabelSetValue(faceLabel, cOff+t, fid);
604:       }
605:       cOff += Nt;
606:     }
607:     EG_free(lobjs);
608:   }
609:   PetscHMapIDestroy(&edgeMap);
610:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
611:   for (c = cStart; c < cEnd; ++c) {
612:     PetscInt *closure = NULL;
613:     PetscInt  clSize, cl, bval, fval;

615:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
616:     DMLabelGetValue(bodyLabel, c, &bval);
617:     DMLabelGetValue(faceLabel, c, &fval);
618:     for (cl = 0; cl < clSize*2; cl += 2) {
619:       DMLabelSetValue(bodyLabel, closure[cl], bval);
620:       DMLabelSetValue(faceLabel, closure[cl], fval);
621:     }
622:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
623:   }
624:   *newdm = dm;
625:   return(0);
626: }
627: #endif

629: /*@C
630:   DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADSLite file.

632:   Collective

634:   Input Parameters:
635: + comm     - The MPI communicator
636: - filename - The name of the ExodusII file

638:   Output Parameter:
639: . dm       - The DM object representing the mesh

641:   Level: beginner

643: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS()
644: @*/
645: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
646: {
647:   PetscMPIInt    rank;
648: #if defined(PETSC_HAVE_EGADS)
649:   ego            context= NULL, model = NULL;
650: #endif
651:   PetscBool      printModel = PETSC_FALSE;

656:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
657:   MPI_Comm_rank(comm, &rank);
658: #if defined(PETSC_HAVE_EGADS)
659:   if (!rank) {

661:     EG_open(&context);
662:     EG_loadModel(context, 0, filename, &model);
663:     if (printModel) {DMPlexEGADSPrintModel(model);}

665:   }
666:   DMPlexCreateEGADS(comm, context, model, dm);
667:   return(0);
668: #else
669:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
670: #endif
671: }