Actual source code: plexegadslite.c

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

  4: #ifdef PETSC_HAVE_EGADS
  5: #include <egads_lite.h>

  7: PetscErrorCode DMPlexSnapToGeomModel_EGADSLite_Internal(DM dm, PetscInt p, PetscInt dE, ego model, PetscInt bodyID, PetscInt faceID, PetscInt edgeID, const PetscScalar mcoords[], PetscScalar gcoords[])
  8: {
  9:   DM             cdm;
 10:   ego           *bodies;
 11:   ego            geom, body, obj;
 12:   /* result has to hold derivatives, along with the value */
 13:   double         params[3], result[18], paramsV[16*3], resultV[16*3], range[4];
 14:   int            Nb, oclass, mtype, *senses, peri;
 15:   Vec            coordinatesLocal;
 16:   PetscScalar   *coords = NULL;
 17:   PetscInt       Nv, v, Np = 0, pm;
 18:   PetscInt       d;

 21:   DMGetCoordinateDM(dm, &cdm);
 22:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
 23:   EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 25:   body = bodies[bodyID];

 27:   if (edgeID >= 0)      {EGlite_objectBodyTopo(body, EDGE, edgeID, &obj); Np = 1;}
 28:   else if (faceID >= 0) {EGlite_objectBodyTopo(body, FACE, faceID, &obj); Np = 2;}
 29:   else {
 30:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 31:     return 0;
 32:   }

 34:   /* Calculate parameters (t or u,v) for vertices */
 35:   DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 36:   Nv  /= dE;
 37:   if (Nv == 1) {
 38:     DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 39:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 40:     return 0;
 41:   }

 44:   /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */
 45:   EGlite_getRange(obj, range, &peri);
 46:   for (v = 0; v < Nv; ++v) {
 47:     EGlite_invEvaluate(obj, &coords[v*dE], &paramsV[v*3], &resultV[v*3]);
 48: #if 1
 49:     if (peri > 0) {
 50:       if      (paramsV[v*3+0] + 1.e-4 < range[0]) {paramsV[v*3+0] += 2. * PETSC_PI;}
 51:       else if (paramsV[v*3+0] - 1.e-4 > range[1]) {paramsV[v*3+0] -= 2. * PETSC_PI;}
 52:     }
 53:     if (peri > 1) {
 54:       if      (paramsV[v*3+1] + 1.e-4 < range[2]) {paramsV[v*3+1] += 2. * PETSC_PI;}
 55:       else if (paramsV[v*3+1] - 1.e-4 > range[3]) {paramsV[v*3+1] -= 2. * PETSC_PI;}
 56:     }
 57: #endif
 58:   }
 59:   DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 60:   /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
 61:   for (pm = 0; pm < Np; ++pm) {
 62:     params[pm] = 0.;
 63:     for (v = 0; v < Nv; ++v) params[pm] += paramsV[v*3+pm];
 64:     params[pm] /= Nv;
 65:   }
 68:   /* Put coordinates for new vertex in result[] */
 69:   EGlite_evaluate(obj, params, result);
 70:   for (d = 0; d < dE; ++d) gcoords[d] = result[d];
 71:   return 0;
 72: }

 74: static PetscErrorCode DMPlexEGADSLiteDestroy_Private(void *context)
 75: {
 76:   if (context) EGlite_close((ego) context);
 77:   return 0;
 78: }

 80: static PetscErrorCode DMPlexCreateEGADSLite_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
 81: {
 82:   DMLabel        bodyLabel, faceLabel, edgeLabel, vertexLabel;
 83:   PetscInt       cStart, cEnd, c;
 84:   /* EGADSLite variables */
 85:   ego            geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
 86:   int            oclass, mtype, nbodies, *senses;
 87:   int            b;
 88:   /* PETSc variables */
 89:   DM             dm;
 90:   PetscHMapI     edgeMap = NULL;
 91:   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;
 92:   PetscInt      *cells  = NULL, *cone = NULL;
 93:   PetscReal     *coords = NULL;
 94:   PetscMPIInt    rank;

 96:   MPI_Comm_rank(comm, &rank);
 97:   if (!rank) {
 98:     const PetscInt debug = 0;

100:     /* ---------------------------------------------------------------------------------------------------
101:     Generate Petsc Plex
102:       Get all Nodes in model, record coordinates in a correctly formatted array
103:       Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
104:       We need to uniformly refine the initial geometry to guarantee a valid mesh
105:     */

107:     /* Calculate cell and vertex sizes */
108:     EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
109:     PetscHMapICreate(&edgeMap);
110:     numEdges = 0;
111:     for (b = 0; b < nbodies; ++b) {
112:       ego body = bodies[b];
113:       int id, Nl, l, Nv, v;

115:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
116:       for (l = 0; l < Nl; ++l) {
117:         ego loop = lobjs[l];
118:         int Ner  = 0, Ne, e, Nc;

120:         EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
121:         for (e = 0; e < Ne; ++e) {
122:           ego edge = objs[e];
123:           int Nv, id;
124:           PetscHashIter iter;
125:           PetscBool     found;

127:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
128:           if (mtype == DEGENERATE) continue;
129:           id   = EGlite_indexBodyTopo(body, edge);
130:           PetscHMapIFind(edgeMap, id-1, &iter, &found);
131:           if (!found) PetscHMapISet(edgeMap, id-1, numEdges++);
132:           ++Ner;
133:         }
134:         if (Ner == 2)      {Nc = 2;}
135:         else if (Ner == 3) {Nc = 4;}
136:         else if (Ner == 4) {Nc = 8; ++numQuads;}
137:         else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
138:         numCells += Nc;
139:         newCells += Nc-1;
140:         maxCorners = PetscMax(Ner*2+1, maxCorners);
141:       }
142:       EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
143:       for (v = 0; v < Nv; ++v) {
144:         ego vertex = nobjs[v];

146:         id = EGlite_indexBodyTopo(body, vertex);
147:         /* TODO: Instead of assuming contiguous ids, we could use a hash table */
148:         numVertices = PetscMax(id, numVertices);
149:       }
150:       EGlite_free(lobjs);
151:       EGlite_free(nobjs);
152:     }
153:     PetscHMapIGetSize(edgeMap, &numEdges);
154:     newVertices  = numEdges + numQuads;
155:     numVertices += newVertices;

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

162:     /* Get vertex coordinates */
163:     for (b = 0; b < nbodies; ++b) {
164:       ego body = bodies[b];
165:       int id, Nv, v;

167:       EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
168:       for (v = 0; v < Nv; ++v) {
169:         ego    vertex = nobjs[v];
170:         double limits[4];
171:         int    dummy;

173:         EGlite_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
174:         id   = EGlite_indexBodyTopo(body, vertex);
175:         coords[(id-1)*cdim+0] = limits[0];
176:         coords[(id-1)*cdim+1] = limits[1];
177:         coords[(id-1)*cdim+2] = limits[2];
178:       }
179:       EGlite_free(nobjs);
180:     }
181:     PetscHMapIClear(edgeMap);
182:     fOff     = numVertices - newVertices + numEdges;
183:     numEdges = 0;
184:     numQuads = 0;
185:     for (b = 0; b < nbodies; ++b) {
186:       ego body = bodies[b];
187:       int Nl, l;

189:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
190:       for (l = 0; l < Nl; ++l) {
191:         ego loop = lobjs[l];
192:         int lid, Ner = 0, Ne, e;

194:         lid  = EGlite_indexBodyTopo(body, loop);
195:         EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
196:         for (e = 0; e < Ne; ++e) {
197:           ego       edge = objs[e];
198:           int       eid, Nv;
199:           PetscHashIter iter;
200:           PetscBool     found;

202:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
203:           if (mtype == DEGENERATE) continue;
204:           ++Ner;
205:           eid  = EGlite_indexBodyTopo(body, edge);
206:           PetscHMapIFind(edgeMap, eid-1, &iter, &found);
207:           if (!found) {
208:             PetscInt v = numVertices - newVertices + numEdges;
209:             double range[4], params[3] = {0., 0., 0.}, result[18];
210:             int    periodic[2];

212:             PetscHMapISet(edgeMap, eid-1, numEdges++);
213:             EGlite_getRange(edge, range, periodic);
214:             params[0] = 0.5*(range[0] + range[1]);
215:             EGlite_evaluate(edge, params, result);
216:             coords[v*cdim+0] = result[0];
217:             coords[v*cdim+1] = result[1];
218:             coords[v*cdim+2] = result[2];
219:           }
220:         }
221:         if (Ner == 4) {
222:           PetscInt v = fOff + numQuads++;
223:           ego     *fobjs, face;
224:           double   range[4], params[3] = {0., 0., 0.}, result[18];
225:           int      Nf, fid, periodic[2];

227:           EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
228:           face = fobjs[0];
229:           fid  = EGlite_indexBodyTopo(body, face);
231:           EGlite_getRange(face, range, periodic);
232:           params[0] = 0.5*(range[0] + range[1]);
233:           params[1] = 0.5*(range[2] + range[3]);
234:           EGlite_evaluate(face, params, result);
235:           coords[v*cdim+0] = result[0];
236:           coords[v*cdim+1] = result[1];
237:           coords[v*cdim+2] = result[2];
238:         }
239:       }
240:     }

243:     /* Get cell vertices by traversing loops */
244:     numQuads = 0;
245:     cOff     = 0;
246:     for (b = 0; b < nbodies; ++b) {
247:       ego body = bodies[b];
248:       int id, Nl, l;

250:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
251:       for (l = 0; l < Nl; ++l) {
252:         ego loop = lobjs[l];
253:         int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;

255:         lid  = EGlite_indexBodyTopo(body, loop);
256:         EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);

258:         for (e = 0; e < Ne; ++e) {
259:           ego edge = objs[e];
260:           int points[3];
261:           int eid, Nv, v, tmp;

263:           eid  = EGlite_indexBodyTopo(body, edge);
264:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
265:           if (mtype == DEGENERATE) continue;
266:           else                     ++Ner;

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

272:             id = EGlite_indexBodyTopo(body, vertex);
273:             points[v*2] = id-1;
274:           }
275:           {
276:             PetscInt edgeNum;

278:             PetscHMapIGet(edgeMap, eid-1, &edgeNum);
279:             points[1] = numVertices - newVertices + edgeNum;
280:           }
281:           /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
282:           if (!nc) {
283:             for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v];
284:           } else {
285:             if (cone[nc-1] == points[0])      {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
286:             else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
287:             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];}
288:             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];}
289:             else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
290:           }
291:         }
293:         if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;}
295:         /* Triangulate the loop */
296:         switch (Ner) {
297:           case 2: /* Bi-Segment -> 2 triangles */
298:             Nt = 2;
299:             cells[cOff*numCorners+0] = cone[0];
300:             cells[cOff*numCorners+1] = cone[1];
301:             cells[cOff*numCorners+2] = cone[2];
302:             ++cOff;
303:             cells[cOff*numCorners+0] = cone[0];
304:             cells[cOff*numCorners+1] = cone[2];
305:             cells[cOff*numCorners+2] = cone[3];
306:             ++cOff;
307:             break;
308:           case 3: /* Triangle   -> 4 triangles */
309:             Nt = 4;
310:             cells[cOff*numCorners+0] = cone[0];
311:             cells[cOff*numCorners+1] = cone[1];
312:             cells[cOff*numCorners+2] = cone[5];
313:             ++cOff;
314:             cells[cOff*numCorners+0] = cone[1];
315:             cells[cOff*numCorners+1] = cone[2];
316:             cells[cOff*numCorners+2] = cone[3];
317:             ++cOff;
318:             cells[cOff*numCorners+0] = cone[5];
319:             cells[cOff*numCorners+1] = cone[3];
320:             cells[cOff*numCorners+2] = cone[4];
321:             ++cOff;
322:             cells[cOff*numCorners+0] = cone[1];
323:             cells[cOff*numCorners+1] = cone[3];
324:             cells[cOff*numCorners+2] = cone[5];
325:             ++cOff;
326:             break;
327:           case 4: /* Quad       -> 8 triangles */
328:             Nt = 8;
329:             cells[cOff*numCorners+0] = cone[0];
330:             cells[cOff*numCorners+1] = cone[1];
331:             cells[cOff*numCorners+2] = cone[7];
332:             ++cOff;
333:             cells[cOff*numCorners+0] = cone[1];
334:             cells[cOff*numCorners+1] = cone[2];
335:             cells[cOff*numCorners+2] = cone[3];
336:             ++cOff;
337:             cells[cOff*numCorners+0] = cone[3];
338:             cells[cOff*numCorners+1] = cone[4];
339:             cells[cOff*numCorners+2] = cone[5];
340:             ++cOff;
341:             cells[cOff*numCorners+0] = cone[5];
342:             cells[cOff*numCorners+1] = cone[6];
343:             cells[cOff*numCorners+2] = cone[7];
344:             ++cOff;
345:             cells[cOff*numCorners+0] = cone[8];
346:             cells[cOff*numCorners+1] = cone[1];
347:             cells[cOff*numCorners+2] = cone[3];
348:             ++cOff;
349:             cells[cOff*numCorners+0] = cone[8];
350:             cells[cOff*numCorners+1] = cone[3];
351:             cells[cOff*numCorners+2] = cone[5];
352:             ++cOff;
353:             cells[cOff*numCorners+0] = cone[8];
354:             cells[cOff*numCorners+1] = cone[5];
355:             cells[cOff*numCorners+2] = cone[7];
356:             ++cOff;
357:             cells[cOff*numCorners+0] = cone[8];
358:             cells[cOff*numCorners+1] = cone[7];
359:             cells[cOff*numCorners+2] = cone[1];
360:             ++cOff;
361:             break;
362:           default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
363:         }
364:         if (debug) {
365:           for (t = 0; t < Nt; ++t) {
366:             PetscPrintf(PETSC_COMM_SELF, "  LOOP Corner NODEs Triangle %D (", t);
367:             for (c = 0; c < numCorners; ++c) {
368:               if (c > 0) PetscPrintf(PETSC_COMM_SELF, ", ");
369:               PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);
370:             }
371:             PetscPrintf(PETSC_COMM_SELF, ")\n");
372:           }
373:         }
374:       }
375:       EGlite_free(lobjs);
376:     }
377:   }
379:   DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
380:   PetscFree3(coords, cells, cone);
381:   PetscInfo(dm, " Total Number of Unique Cells    = %D (%D)\n", numCells, newCells);
382:   PetscInfo(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices);
383:   /* Embed EGADS model in DM */
384:   {
385:     PetscContainer modelObj, contextObj;

387:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
388:     PetscContainerSetPointer(modelObj, model);
389:     PetscObjectCompose((PetscObject) dm, "EGADSLite Model", (PetscObject) modelObj);
390:     PetscContainerDestroy(&modelObj);

392:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
393:     PetscContainerSetPointer(contextObj, context);
394:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSLiteDestroy_Private);
395:     PetscObjectCompose((PetscObject) dm, "EGADSLite Context", (PetscObject) contextObj);
396:     PetscContainerDestroy(&contextObj);
397:   }
398:   /* Label points */
399:   DMCreateLabel(dm, "EGADS Body ID");
400:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
401:   DMCreateLabel(dm, "EGADS Face ID");
402:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
403:   DMCreateLabel(dm, "EGADS Edge ID");
404:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
405:   DMCreateLabel(dm, "EGADS Vertex ID");
406:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
407:   cOff = 0;
408:   for (b = 0; b < nbodies; ++b) {
409:     ego body = bodies[b];
410:     int id, Nl, l;

412:     EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
413:     for (l = 0; l < Nl; ++l) {
414:       ego  loop = lobjs[l];
415:       ego *fobjs;
416:       int  lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;

418:       lid  = EGlite_indexBodyTopo(body, loop);
419:       EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
421:       fid  = EGlite_indexBodyTopo(body, fobjs[0]);
422:       EGlite_free(fobjs);
423:       EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
424:       for (e = 0; e < Ne; ++e) {
425:         ego             edge = objs[e];
426:         int             eid, Nv, v;
427:         PetscInt        points[3], support[2], numEdges, edgeNum;
428:         const PetscInt *edges;

430:         eid  = EGlite_indexBodyTopo(body, edge);
431:         EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
432:         if (mtype == DEGENERATE) continue;
433:         else                     ++Ner;
434:         for (v = 0; v < Nv; ++v) {
435:           ego vertex = nobjs[v];

437:           id   = EGlite_indexBodyTopo(body, vertex);
438:           DMLabelSetValue(edgeLabel, numCells + id-1, eid);
439:           points[v*2] = numCells + id-1;
440:         }
441:         PetscHMapIGet(edgeMap, eid-1, &edgeNum);
442:         points[1] = numCells + numVertices - newVertices + edgeNum;

444:         DMLabelSetValue(edgeLabel, points[1], eid);
445:         support[0] = points[0];
446:         support[1] = points[1];
447:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
449:         DMLabelSetValue(edgeLabel, edges[0], eid);
450:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
451:         support[0] = points[1];
452:         support[1] = points[2];
453:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
455:         DMLabelSetValue(edgeLabel, edges[0], eid);
456:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
457:       }
458:       switch (Ner) {
459:         case 2: Nt = 2;break;
460:         case 3: Nt = 4;break;
461:         case 4: Nt = 8;break;
462:         default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
463:       }
464:       for (t = 0; t < Nt; ++t) {
465:         DMLabelSetValue(bodyLabel, cOff+t, b);
466:         DMLabelSetValue(faceLabel, cOff+t, fid);
467:       }
468:       cOff += Nt;
469:     }
470:     EGlite_free(lobjs);
471:   }
472:   PetscHMapIDestroy(&edgeMap);
473:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
474:   for (c = cStart; c < cEnd; ++c) {
475:     PetscInt *closure = NULL;
476:     PetscInt  clSize, cl, bval, fval;

478:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
479:     DMLabelGetValue(bodyLabel, c, &bval);
480:     DMLabelGetValue(faceLabel, c, &fval);
481:     for (cl = 0; cl < clSize*2; cl += 2) {
482:       DMLabelSetValue(bodyLabel, closure[cl], bval);
483:       DMLabelSetValue(faceLabel, closure[cl], fval);
484:     }
485:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
486:   }
487:   *newdm = dm;
488:   return 0;
489: }

491: static PetscErrorCode DMPlexEGADSLitePrintModel_Internal(ego model)
492: {
493:   ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
494:   int oclass, mtype, *senses;
495:   int Nb, b;

498:   /* test bodyTopo functions */
499:   EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
500:   PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb);

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

506:     /* Output Basic Model Topology */
507:     EGlite_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
508:     PetscPrintf(PETSC_COMM_SELF, "   Number of SHELLS: %d \n", Nsh);
509:     EGlite_free(objs);

511:     EGlite_getBodyTopos(body, NULL, FACE,  &Nf, &objs);
512:     PetscPrintf(PETSC_COMM_SELF, "   Number of FACES: %d \n", Nf);
513:     EGlite_free(objs);

515:     EGlite_getBodyTopos(body, NULL, LOOP,  &Nl, &lobjs);
516:     PetscPrintf(PETSC_COMM_SELF, "   Number of LOOPS: %d \n", Nl);

518:     EGlite_getBodyTopos(body, NULL, EDGE,  &Ne, &objs);
519:     PetscPrintf(PETSC_COMM_SELF, "   Number of EDGES: %d \n", Ne);
520:     EGlite_free(objs);

522:     EGlite_getBodyTopos(body, NULL, NODE,  &Nv, &objs);
523:     PetscPrintf(PETSC_COMM_SELF, "   Number of NODES: %d \n", Nv);
524:     EGlite_free(objs);

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

529:       id   = EGlite_indexBodyTopo(body, loop);
530:       PetscPrintf(PETSC_COMM_SELF, "          LOOP ID: %d\n", id);

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

535:       for (e = 0; e < Ne; ++e) {
536:         ego    edge      = objs[e];
537:         double range[4]  = {0., 0., 0., 0.};
538:         double point[3]  = {0., 0., 0.};
539:         double params[3] = {0., 0., 0.};
540:         double result[18];
541:         int    peri;

543:         EGlite_indexBodyTopo(body, edge);
544:         PetscPrintf(PETSC_COMM_SELF, "            EDGE ID: %d (%d)\n", id, e);

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

549:         /* Get NODE info which associated with the current EDGE */
550:         EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
551:         if (mtype == DEGENERATE) {
552:           PetscPrintf(PETSC_COMM_SELF, "  EDGE %d is DEGENERATE \n", id);
553:         } else {
554:           params[0] = range[0];
555:           EGlite_evaluate(edge, params, result);
556:           PetscPrintf(PETSC_COMM_SELF, "   between (%lf, %lf, %lf)", result[0], result[1], result[2]);
557:           params[0] = range[1];
558:           EGlite_evaluate(edge, params, result);
559:           PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]);
560:         }

562:         for (v = 0; v < Nv; ++v) {
563:           ego    vertex = nobjs[v];
564:           double limits[4];
565:           int    dummy;

567:           EGlite_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
568:           EGlite_indexBodyTopo(body, vertex);
569:           PetscPrintf(PETSC_COMM_SELF, "              NODE ID: %d \n", id);
570:           PetscPrintf(PETSC_COMM_SELF, "                 (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]);

572:           point[0] = point[0] + limits[0];
573:           point[1] = point[1] + limits[1];
574:           point[2] = point[2] + limits[2];
575:         }
576:       }
577:     }
578:     EGlite_free(lobjs);
579:   }
580:   return 0;
581: }
582: #endif

584: /*@C
585:   DMPlexCreateEGADSLiteFromFile - Create a DMPlex mesh from an EGADSLite file.

587:   Collective

589:   Input Parameters:
590: + comm     - The MPI communicator
591: - filename - The name of the EGADSLite file

593:   Output Parameter:
594: . dm       - The DM object representing the mesh

596:   Level: beginner

598: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS(), DMPlexCreateEGADSFromFile()
599: @*/
600: PetscErrorCode DMPlexCreateEGADSLiteFromFile(MPI_Comm comm, const char filename[], DM *dm)
601: {
602:   PetscMPIInt    rank;
603: #if defined(PETSC_HAVE_EGADS)
604:   ego            context= NULL, model = NULL;
605: #endif
606:   PetscBool      printModel = PETSC_FALSE;

609:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
610:   MPI_Comm_rank(comm, &rank);
611: #if defined(PETSC_HAVE_EGADS)
612:   if (!rank) {

614:     EGlite_open(&context);
615:     EGlite_loadModel(context, 0, filename, &model);
616:     if (printModel) DMPlexEGADSLitePrintModel_Internal(model);

618:   }
619:   DMPlexCreateEGADSLite_Internal(comm, context, model, dm);
620:   return 0;
621: #else
622:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
623: #endif
624: }