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