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