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

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

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

 84: static PetscErrorCode DMPlexEGADSLiteDestroy_Private(void *context)
 85: {
 86:   if (context) EGlite_close((ego)context);
 87:   return 0;
 88: }

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

106:   MPI_Comm_rank(comm, &rank);
107:   if (rank == 0) {
108:     const PetscInt debug = 0;

110:     /* ---------------------------------------------------------------------------------------------------
111:     Generate Petsc Plex
112:       Get all Nodes in model, record coordinates in a correctly formatted array
113:       Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
114:       We need to uniformly refine the initial geometry to guarantee a valid mesh
115:     */

117:     /* Calculate cell and vertex sizes */
118:     EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
119:     PetscHMapICreate(&edgeMap);
120:     numEdges = 0;
121:     for (b = 0; b < nbodies; ++b) {
122:       ego body = bodies[b];
123:       int id, Nl, l, Nv, v;

125:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
126:       for (l = 0; l < Nl; ++l) {
127:         ego loop = lobjs[l];
128:         int Ner  = 0, Ne, e, Nc;

130:         EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
131:         for (e = 0; e < Ne; ++e) {
132:           ego           edge = objs[e];
133:           int           Nv, id;
134:           PetscHashIter iter;
135:           PetscBool     found;

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

160:         id = EGlite_indexBodyTopo(body, vertex);
161:         /* TODO: Instead of assuming contiguous ids, we could use a hash table */
162:         numVertices = PetscMax(id, numVertices);
163:       }
164:       EGlite_free(lobjs);
165:       EGlite_free(nobjs);
166:     }
167:     PetscHMapIGetSize(edgeMap, &numEdges);
168:     newVertices = numEdges + numQuads;
169:     numVertices += newVertices;

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

176:     /* Get vertex coordinates */
177:     for (b = 0; b < nbodies; ++b) {
178:       ego body = bodies[b];
179:       int id, Nv, v;

181:       EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
182:       for (v = 0; v < Nv; ++v) {
183:         ego    vertex = nobjs[v];
184:         double limits[4];
185:         int    dummy;

187:         EGlite_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
188:         id                          = EGlite_indexBodyTopo(body, vertex);
189:         coords[(id - 1) * cdim + 0] = limits[0];
190:         coords[(id - 1) * cdim + 1] = limits[1];
191:         coords[(id - 1) * cdim + 2] = limits[2];
192:       }
193:       EGlite_free(nobjs);
194:     }
195:     PetscHMapIClear(edgeMap);
196:     fOff     = numVertices - newVertices + numEdges;
197:     numEdges = 0;
198:     numQuads = 0;
199:     for (b = 0; b < nbodies; ++b) {
200:       ego body = bodies[b];
201:       int Nl, l;

203:       EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
204:       for (l = 0; l < Nl; ++l) {
205:         ego loop = lobjs[l];
206:         int lid, Ner = 0, Ne, e;

208:         lid = EGlite_indexBodyTopo(body, loop);
209:         EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
210:         for (e = 0; e < Ne; ++e) {
211:           ego           edge = objs[e];
212:           int           eid, Nv;
213:           PetscHashIter iter;
214:           PetscBool     found;

216:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
217:           if (mtype == DEGENERATE) continue;
218:           ++Ner;
219:           eid = EGlite_indexBodyTopo(body, edge);
220:           PetscHMapIFind(edgeMap, eid - 1, &iter, &found);
221:           if (!found) {
222:             PetscInt v = numVertices - newVertices + numEdges;
223:             double   range[4], params[3] = {0., 0., 0.}, result[18];
224:             int      periodic[2];

226:             PetscHMapISet(edgeMap, eid - 1, numEdges++);
227:             EGlite_getRange(edge, range, periodic);
228:             params[0] = 0.5 * (range[0] + range[1]);
229:             EGlite_evaluate(edge, params, result);
230:             coords[v * cdim + 0] = result[0];
231:             coords[v * cdim + 1] = result[1];
232:             coords[v * cdim + 2] = result[2];
233:           }
234:         }
235:         if (Ner == 4) {
236:           PetscInt v = fOff + numQuads++;
237:           ego     *fobjs, face;
238:           double   range[4], params[3] = {0., 0., 0.}, result[18];
239:           int      Nf, fid, periodic[2];

241:           EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
242:           face = fobjs[0];
243:           fid  = EGlite_indexBodyTopo(body, face);
245:           EGlite_getRange(face, range, periodic);
246:           params[0] = 0.5 * (range[0] + range[1]);
247:           params[1] = 0.5 * (range[2] + range[3]);
248:           EGlite_evaluate(face, params, result);
249:           coords[v * cdim + 0] = result[0];
250:           coords[v * cdim + 1] = result[1];
251:           coords[v * cdim + 2] = result[2];
252:         }
253:       }
254:     }

257:     /* Get cell vertices by traversing loops */
258:     numQuads = 0;
259:     cOff     = 0;
260:     for (b = 0; b < nbodies; ++b) {
261:       ego body = bodies[b];
262:       int id, Nl, l;

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

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

272:         for (e = 0; e < Ne; ++e) {
273:           ego edge = objs[e];
274:           int points[3];
275:           int eid, Nv, v, tmp;

277:           eid = EGlite_indexBodyTopo(body, edge);
278:           EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
279:           if (mtype == DEGENERATE) continue;
280:           else ++Ner;

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

286:             id            = EGlite_indexBodyTopo(body, vertex);
287:             points[v * 2] = id - 1;
288:           }
289:           {
290:             PetscInt edgeNum;

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

416:     PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
417:     PetscContainerSetPointer(modelObj, model);
418:     PetscObjectCompose((PetscObject)dm, "EGADSLite Model", (PetscObject)modelObj);
419:     PetscContainerDestroy(&modelObj);

421:     PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
422:     PetscContainerSetPointer(contextObj, context);
423:     PetscContainerSetUserDestroy(contextObj, DMPlexEGADSLiteDestroy_Private);
424:     PetscObjectCompose((PetscObject)dm, "EGADSLite Context", (PetscObject)contextObj);
425:     PetscContainerDestroy(&contextObj);
426:   }
427:   /* Label points */
428:   DMCreateLabel(dm, "EGADS Body ID");
429:   DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
430:   DMCreateLabel(dm, "EGADS Face ID");
431:   DMGetLabel(dm, "EGADS Face ID", &faceLabel);
432:   DMCreateLabel(dm, "EGADS Edge ID");
433:   DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
434:   DMCreateLabel(dm, "EGADS Vertex ID");
435:   DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
436:   cOff = 0;
437:   for (b = 0; b < nbodies; ++b) {
438:     ego body = bodies[b];
439:     int id, Nl, l;

441:     EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
442:     for (l = 0; l < Nl; ++l) {
443:       ego  loop = lobjs[l];
444:       ego *fobjs;
445:       int  lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;

447:       lid = EGlite_indexBodyTopo(body, loop);
448:       EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
450:       fid = EGlite_indexBodyTopo(body, fobjs[0]);
451:       EGlite_free(fobjs);
452:       EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
453:       for (e = 0; e < Ne; ++e) {
454:         ego             edge = objs[e];
455:         int             eid, Nv, v;
456:         PetscInt        points[3], support[2], numEdges, edgeNum;
457:         const PetscInt *edges;

459:         eid = EGlite_indexBodyTopo(body, edge);
460:         EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
461:         if (mtype == DEGENERATE) continue;
462:         else ++Ner;
463:         for (v = 0; v < Nv; ++v) {
464:           ego vertex = nobjs[v];

466:           id = EGlite_indexBodyTopo(body, vertex);
467:           DMLabelSetValue(edgeLabel, numCells + id - 1, eid);
468:           points[v * 2] = numCells + id - 1;
469:         }
470:         PetscHMapIGet(edgeMap, eid - 1, &edgeNum);
471:         points[1] = numCells + numVertices - newVertices + edgeNum;

473:         DMLabelSetValue(edgeLabel, points[1], eid);
474:         support[0] = points[0];
475:         support[1] = points[1];
476:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
478:         DMLabelSetValue(edgeLabel, edges[0], eid);
479:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
480:         support[0] = points[1];
481:         support[1] = points[2];
482:         DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
484:         DMLabelSetValue(edgeLabel, edges[0], eid);
485:         DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
486:       }
487:       switch (Ner) {
488:       case 2:
489:         Nt = 2;
490:         break;
491:       case 3:
492:         Nt = 4;
493:         break;
494:       case 4:
495:         Nt = 8;
496:         break;
497:       default:
498:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
499:       }
500:       for (t = 0; t < Nt; ++t) {
501:         DMLabelSetValue(bodyLabel, cOff + t, b);
502:         DMLabelSetValue(faceLabel, cOff + t, fid);
503:       }
504:       cOff += Nt;
505:     }
506:     EGlite_free(lobjs);
507:   }
508:   PetscHMapIDestroy(&edgeMap);
509:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
510:   for (c = cStart; c < cEnd; ++c) {
511:     PetscInt *closure = NULL;
512:     PetscInt  clSize, cl, bval, fval;

514:     DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
515:     DMLabelGetValue(bodyLabel, c, &bval);
516:     DMLabelGetValue(faceLabel, c, &fval);
517:     for (cl = 0; cl < clSize * 2; cl += 2) {
518:       DMLabelSetValue(bodyLabel, closure[cl], bval);
519:       DMLabelSetValue(faceLabel, closure[cl], fval);
520:     }
521:     DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
522:   }
523:   *newdm = dm;
524:   return 0;
525: }

527: static PetscErrorCode DMPlexEGADSLitePrintModel_Internal(ego model)
528: {
529:   ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
530:   int oclass, mtype, *senses;
531:   int Nb, b;

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

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

542:     /* Output Basic Model Topology */
543:     EGlite_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
544:     PetscPrintf(PETSC_COMM_SELF, "   Number of SHELLS: %d \n", Nsh);
545:     EGlite_free(objs);

547:     EGlite_getBodyTopos(body, NULL, FACE, &Nf, &objs);
548:     PetscPrintf(PETSC_COMM_SELF, "   Number of FACES: %d \n", Nf);
549:     EGlite_free(objs);

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

554:     EGlite_getBodyTopos(body, NULL, EDGE, &Ne, &objs);
555:     PetscPrintf(PETSC_COMM_SELF, "   Number of EDGES: %d \n", Ne);
556:     EGlite_free(objs);

558:     EGlite_getBodyTopos(body, NULL, NODE, &Nv, &objs);
559:     PetscPrintf(PETSC_COMM_SELF, "   Number of NODES: %d \n", Nv);
560:     EGlite_free(objs);

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

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

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

571:       for (e = 0; e < Ne; ++e) {
572:         ego    edge      = objs[e];
573:         double range[4]  = {0., 0., 0., 0.};
574:         double point[3]  = {0., 0., 0.};
575:         double params[3] = {0., 0., 0.};
576:         double result[18];
577:         int    peri;

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

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

585:         /* Get NODE info which associated with the current EDGE */
586:         EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
587:         if (mtype == DEGENERATE) {
588:           PetscPrintf(PETSC_COMM_SELF, "  EDGE %d is DEGENERATE \n", id);
589:         } else {
590:           params[0] = range[0];
591:           EGlite_evaluate(edge, params, result);
592:           PetscPrintf(PETSC_COMM_SELF, "   between (%lf, %lf, %lf)", result[0], result[1], result[2]);
593:           params[0] = range[1];
594:           EGlite_evaluate(edge, params, result);
595:           PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]);
596:         }

598:         for (v = 0; v < Nv; ++v) {
599:           ego    vertex = nobjs[v];
600:           double limits[4];
601:           int    dummy;

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

608:           point[0] = point[0] + limits[0];
609:           point[1] = point[1] + limits[1];
610:           point[2] = point[2] + limits[2];
611:         }
612:       }
613:     }
614:     EGlite_free(lobjs);
615:   }
616:   return 0;
617: }
618: #endif

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

623:   Collective

625:   Input Parameters:
626: + comm     - The MPI communicator
627: - filename - The name of the EGADSLite file

629:   Output Parameter:
630: . dm       - The DM object representing the mesh

632:   Level: beginner

634: .seealso: `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSFromFile()`
635: @*/
636: PetscErrorCode DMPlexCreateEGADSLiteFromFile(MPI_Comm comm, const char filename[], DM *dm)
637: {
638:   PetscMPIInt rank;
639: #if defined(PETSC_HAVE_EGADS)
640:   ego context = NULL, model = NULL;
641: #endif
642:   PetscBool printModel = PETSC_FALSE;

645:   PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
646:   MPI_Comm_rank(comm, &rank);
647: #if defined(PETSC_HAVE_EGADS)
648:   if (rank == 0) {
649:     EGlite_open(&context);
650:     EGlite_loadModel(context, 0, filename, &model);
651:     if (printModel) DMPlexEGADSLitePrintModel_Internal(model);
652:   }
653:   DMPlexCreateEGADSLite_Internal(comm, context, model, dm);
654:   return 0;
655: #else
656:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
657: #endif
658: }