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, 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       dE, d;

 22:   DMGetCoordinateDM(dm, &cdm);
 23:   DMGetCoordinateDim(dm, &dE);
 24:   DMGetCoordinatesLocal(dm, &coordinatesLocal);
 25:   EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
 26:   if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
 27:   body = bodies[bodyID];

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

 36:   /* Calculate parameters (t or u,v) for vertices */
 37:   DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 38:   Nv  /= dE;
 39:   if (Nv == 1) {
 40:     DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
 41:     for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
 42:     return(0);
 43:   }
 44:   if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);

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

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

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

100:   MPI_Comm_rank(comm, &rank);
101:   if (!rank) {
102:     const PetscInt debug = 0;

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

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

119:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
120:       for (l = 0; l < Nl; ++l) {
121:         ego loop = lobjs[l];
122:         int Ner  = 0, Ne, e, Nc;

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

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

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

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

166:     /* Get vertex coordinates */
167:     for (b = 0; b < nbodies; ++b) {
168:       ego body = bodies[b];
169:       int id, Nv, v;

171:       EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
172:       for (v = 0; v < Nv; ++v) {
173:         ego    vertex = nobjs[v];
174:         double limits[4];
175:         int    dummy;

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

193:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
194:       for (l = 0; l < Nl; ++l) {
195:         ego loop = lobjs[l];
196:         int lid, Ner = 0, Ne, e;

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

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

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

231:           EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
232:           face = fobjs[0];
233:           fid  = EGlite_indexBodyTopo(body, face);
234:           if (Nf != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid);
235:           EGlite_getRange(face, range, periodic);
236:           params[0] = 0.5*(range[0] + range[1]);
237:           params[1] = 0.5*(range[2] + range[3]);
238:           EGlite_evaluate(face, params, result);
239:           coords[v*cdim+0] = result[0];
240:           coords[v*cdim+1] = result[1];
241:           coords[v*cdim+2] = result[2];
242:         }
243:       }
244:     }
245:     if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices);

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

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

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

262:         for (e = 0; e < Ne; ++e) {
263:           ego edge = objs[e];
264:           int points[3];
265:           int eid, Nv, v, tmp;

267:           eid  = EGlite_indexBodyTopo(body, edge);
268:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
269:           if (mtype == DEGENERATE) continue;
270:           else                     ++Ner;
271:           if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);

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

276:             id = EGlite_indexBodyTopo(body, vertex);
277:             points[v*2] = id-1;
278:           }
279:           {
280:             PetscInt edgeNum;

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

391:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
392:     PetscContainerSetPointer(modelObj, model);
393:     PetscObjectCompose((PetscObject) dm, "EGADSLite Model", (PetscObject) modelObj);
394:     PetscContainerDestroy(&modelObj);

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

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

422:       lid  = EGlite_indexBodyTopo(body, loop);
423:       EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
424:       if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
425:       fid  = EGlite_indexBodyTopo(body, fobjs[0]);
426:       EGlite_free(fobjs);
427:       EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
428:       for (e = 0; e < Ne; ++e) {
429:         ego             edge = objs[e];
430:         int             eid, Nv, v;
431:         PetscInt        points[3], support[2], numEdges, edgeNum;
432:         const PetscInt *edges;

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

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

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

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

495: static PetscErrorCode DMPlexEGADSLitePrintModel_Internal(ego model)
496: {
497:   ego            geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
498:   int            oclass, mtype, *senses;
499:   int            Nb, b;

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

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

511:     /* Output Basic Model Topology */
512:     EGlite_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
513:     PetscPrintf(PETSC_COMM_SELF, "   Number of SHELLS: %d \n", Nsh);
514:     EGlite_free(objs);

516:     EGlite_getBodyTopos(body, NULL, FACE,  &Nf, &objs);
517:     PetscPrintf(PETSC_COMM_SELF, "   Number of FACES: %d \n", Nf);
518:     EGlite_free(objs);

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

523:     EGlite_getBodyTopos(body, NULL, EDGE,  &Ne, &objs);
524:     PetscPrintf(PETSC_COMM_SELF, "   Number of EDGES: %d \n", Ne);
525:     EGlite_free(objs);

527:     EGlite_getBodyTopos(body, NULL, NODE,  &Nv, &objs);
528:     PetscPrintf(PETSC_COMM_SELF, "   Number of NODES: %d \n", Nv);
529:     EGlite_free(objs);

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

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

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

540:       for (e = 0; e < Ne; ++e) {
541:         ego    edge      = objs[e];
542:         double range[4]  = {0., 0., 0., 0.};
543:         double point[3]  = {0., 0., 0.};
544:         double params[3] = {0., 0., 0.};
545:         double result[18];
546:         int    peri;

548:         id   = EGlite_indexBodyTopo(body, edge);
549:         PetscPrintf(PETSC_COMM_SELF, "            EDGE ID: %d (%d)\n", id, e);

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

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

567:         for (v = 0; v < Nv; ++v) {
568:           ego    vertex = nobjs[v];
569:           double limits[4];
570:           int    dummy;

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

577:           point[0] = point[0] + limits[0];
578:           point[1] = point[1] + limits[1];
579:           point[2] = point[2] + limits[2];
580:         }
581:       }
582:     }
583:     EGlite_free(lobjs);
584:   }
585:   return(0);
586: }
587: #endif

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

592:   Collective

594:   Input Parameters:
595: + comm     - The MPI communicator
596: - filename - The name of the EGADSLite file

598:   Output Parameter:
599: . dm       - The DM object representing the mesh

601:   Level: beginner

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

616:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
617:   MPI_Comm_rank(comm, &rank);
618: #if defined(PETSC_HAVE_EGADS)
619:   if (!rank) {

621:     EGlite_open(&context);
622:     EGlite_loadModel(context, 0, filename, &model);
623:     if (printModel) {DMPlexEGADSLitePrintModel_Internal(model);}

625:   }
626:   DMPlexCreateEGADSLite_Internal(comm, context, model, dm);
627:   return(0);
628: #else
629:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
630: #endif
631: }