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], ¶msV[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: }