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;
20: PetscFunctionBeginHot;
21: PetscCall(DMGetCoordinateDM(dm, &cdm));
22: PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
23: PetscCall(EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
24: PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
25: body = bodies[bodyID];
27: if (edgeID >= 0) {
28: PetscCall(EGlite_objectBodyTopo(body, EDGE, edgeID, &obj));
29: Np = 1;
30: } else if (faceID >= 0) {
31: PetscCall(EGlite_objectBodyTopo(body, FACE, faceID, &obj));
32: Np = 2;
33: } else {
34: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: /* Calculate parameters (t or u,v) for vertices */
39: PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
40: Nv /= dE;
41: if (Nv == 1) {
42: PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
43: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
46: PetscCheck(Nv <= 16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %" PetscInt_FMT " coordinates associated to point %" PetscInt_FMT, Nv, p);
48: /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */
49: PetscCall(EGlite_getRange(obj, range, &peri));
50: for (v = 0; v < Nv; ++v) {
51: PetscCall(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: PetscCall(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: }
76: PetscCheck(!(params[0] < range[0]) && !(params[0] > range[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p);
77: PetscCheck(Np <= 1 || (params[1] >= range[2] && params[1] <= range[3]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p);
78: /* Put coordinates for new vertex in result[] */
79: PetscCall(EGlite_evaluate(obj, params, result));
80: for (d = 0; d < dE; ++d) gcoords[d] = result[d];
81: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
107: PetscCallMPI(MPI_Comm_rank(comm, &rank));
108: if (rank == 0) {
109: const PetscInt debug = 0;
111: /* ---------------------------------------------------------------------------------------------------
112: Generate Petsc Plex
113: Get all Nodes in model, record coordinates in a correctly formatted array
114: Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
115: We need to uniformly refine the initial geometry to guarantee a valid mesh
116: */
118: /* Calculate cell and vertex sizes */
119: PetscCall(EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
120: PetscCall(PetscHMapICreate(&edgeMap));
121: numEdges = 0;
122: for (b = 0; b < nbodies; ++b) {
123: ego body = bodies[b];
124: int id, Nl, l, Nv, v;
126: PetscCall(EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
127: for (l = 0; l < Nl; ++l) {
128: ego loop = lobjs[l];
129: int Ner = 0, Ne, e, Nc;
131: PetscCall(EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
132: for (e = 0; e < Ne; ++e) {
133: ego edge = objs[e];
134: int Nv, id;
135: PetscHashIter iter;
136: PetscBool found;
138: PetscCall(EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
139: if (mtype == DEGENERATE) continue;
140: id = EGlite_indexBodyTopo(body, edge);
141: PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found));
142: if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++));
143: ++Ner;
144: }
145: if (Ner == 2) {
146: Nc = 2;
147: } else if (Ner == 3) {
148: Nc = 4;
149: } else if (Ner == 4) {
150: Nc = 8;
151: ++numQuads;
152: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
153: numCells += Nc;
154: newCells += Nc - 1;
155: maxCorners = PetscMax(Ner * 2 + 1, maxCorners);
156: }
157: PetscCall(EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
158: for (v = 0; v < Nv; ++v) {
159: ego vertex = nobjs[v];
161: id = EGlite_indexBodyTopo(body, vertex);
162: /* TODO: Instead of assuming contiguous ids, we could use a hash table */
163: numVertices = PetscMax(id, numVertices);
164: }
165: EGlite_free(lobjs);
166: EGlite_free(nobjs);
167: }
168: PetscCall(PetscHMapIGetSize(edgeMap, &numEdges));
169: newVertices = numEdges + numQuads;
170: numVertices += newVertices;
172: dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
173: cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
174: numCorners = 3; /* Split cells into triangles */
175: PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone));
177: /* Get vertex coordinates */
178: for (b = 0; b < nbodies; ++b) {
179: ego body = bodies[b];
180: int id, Nv, v;
182: PetscCall(EGlite_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
183: for (v = 0; v < Nv; ++v) {
184: ego vertex = nobjs[v];
185: double limits[4];
186: int dummy;
188: PetscCall(EGlite_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
189: id = EGlite_indexBodyTopo(body, vertex);
190: coords[(id - 1) * cdim + 0] = limits[0];
191: coords[(id - 1) * cdim + 1] = limits[1];
192: coords[(id - 1) * cdim + 2] = limits[2];
193: }
194: EGlite_free(nobjs);
195: }
196: PetscCall(PetscHMapIClear(edgeMap));
197: fOff = numVertices - newVertices + numEdges;
198: numEdges = 0;
199: numQuads = 0;
200: for (b = 0; b < nbodies; ++b) {
201: ego body = bodies[b];
202: int Nl, l;
204: PetscCall(EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
205: for (l = 0; l < Nl; ++l) {
206: ego loop = lobjs[l];
207: int lid, Ner = 0, Ne, e;
209: lid = EGlite_indexBodyTopo(body, loop);
210: PetscCall(EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
211: for (e = 0; e < Ne; ++e) {
212: ego edge = objs[e];
213: int eid, Nv;
214: PetscHashIter iter;
215: PetscBool found;
217: PetscCall(EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
218: if (mtype == DEGENERATE) continue;
219: ++Ner;
220: eid = EGlite_indexBodyTopo(body, edge);
221: PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found));
222: if (!found) {
223: PetscInt v = numVertices - newVertices + numEdges;
224: double range[4], params[3] = {0., 0., 0.}, result[18];
225: int periodic[2];
227: PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++));
228: PetscCall(EGlite_getRange(edge, range, periodic));
229: params[0] = 0.5 * (range[0] + range[1]);
230: PetscCall(EGlite_evaluate(edge, params, result));
231: coords[v * cdim + 0] = result[0];
232: coords[v * cdim + 1] = result[1];
233: coords[v * cdim + 2] = result[2];
234: }
235: }
236: if (Ner == 4) {
237: PetscInt v = fOff + numQuads++;
238: ego *fobjs, face;
239: double range[4], params[3] = {0., 0., 0.}, result[18];
240: int Nf, fid, periodic[2];
242: PetscCall(EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
243: face = fobjs[0];
244: fid = EGlite_indexBodyTopo(body, face);
245: PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid);
246: PetscCall(EGlite_getRange(face, range, periodic));
247: params[0] = 0.5 * (range[0] + range[1]);
248: params[1] = 0.5 * (range[2] + range[3]);
249: PetscCall(EGlite_evaluate(face, params, result));
250: coords[v * cdim + 0] = result[0];
251: coords[v * cdim + 1] = result[1];
252: coords[v * cdim + 2] = result[2];
253: }
254: }
255: }
256: PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices);
258: /* Get cell vertices by traversing loops */
259: numQuads = 0;
260: cOff = 0;
261: for (b = 0; b < nbodies; ++b) {
262: ego body = bodies[b];
263: int id, Nl, l;
265: PetscCall(EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
266: for (l = 0; l < Nl; ++l) {
267: ego loop = lobjs[l];
268: int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;
270: lid = EGlite_indexBodyTopo(body, loop);
271: PetscCall(EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
273: for (e = 0; e < Ne; ++e) {
274: ego edge = objs[e];
275: int points[3];
276: int eid, Nv, v, tmp;
278: eid = EGlite_indexBodyTopo(body, edge);
279: PetscCall(EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
280: if (mtype == DEGENERATE) continue;
281: else ++Ner;
282: PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);
284: for (v = 0; v < Nv; ++v) {
285: ego vertex = nobjs[v];
287: id = EGlite_indexBodyTopo(body, vertex);
288: points[v * 2] = id - 1;
289: }
290: {
291: PetscInt edgeNum;
293: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
294: points[1] = numVertices - newVertices + edgeNum;
295: }
296: /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
297: if (!nc) {
298: for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v];
299: } else {
300: if (cone[nc - 1] == points[0]) {
301: cone[nc++] = points[1];
302: if (cone[0] != points[2]) cone[nc++] = points[2];
303: } else if (cone[nc - 1] == points[2]) {
304: cone[nc++] = points[1];
305: if (cone[0] != points[0]) cone[nc++] = points[0];
306: } else if (cone[nc - 3] == points[0]) {
307: tmp = cone[nc - 3];
308: cone[nc - 3] = cone[nc - 1];
309: cone[nc - 1] = tmp;
310: cone[nc++] = points[1];
311: if (cone[0] != points[2]) cone[nc++] = points[2];
312: } else if (cone[nc - 3] == points[2]) {
313: tmp = cone[nc - 3];
314: cone[nc - 3] = cone[nc - 1];
315: cone[nc - 1] = tmp;
316: cone[nc++] = points[1];
317: if (cone[0] != points[0]) cone[nc++] = points[0];
318: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
319: }
320: }
321: PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner);
322: if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++;
323: PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners);
324: /* Triangulate the loop */
325: switch (Ner) {
326: case 2: /* Bi-Segment -> 2 triangles */
327: Nt = 2;
328: cells[cOff * numCorners + 0] = cone[0];
329: cells[cOff * numCorners + 1] = cone[1];
330: cells[cOff * numCorners + 2] = cone[2];
331: ++cOff;
332: cells[cOff * numCorners + 0] = cone[0];
333: cells[cOff * numCorners + 1] = cone[2];
334: cells[cOff * numCorners + 2] = cone[3];
335: ++cOff;
336: break;
337: case 3: /* Triangle -> 4 triangles */
338: Nt = 4;
339: cells[cOff * numCorners + 0] = cone[0];
340: cells[cOff * numCorners + 1] = cone[1];
341: cells[cOff * numCorners + 2] = cone[5];
342: ++cOff;
343: cells[cOff * numCorners + 0] = cone[1];
344: cells[cOff * numCorners + 1] = cone[2];
345: cells[cOff * numCorners + 2] = cone[3];
346: ++cOff;
347: cells[cOff * numCorners + 0] = cone[5];
348: cells[cOff * numCorners + 1] = cone[3];
349: cells[cOff * numCorners + 2] = cone[4];
350: ++cOff;
351: cells[cOff * numCorners + 0] = cone[1];
352: cells[cOff * numCorners + 1] = cone[3];
353: cells[cOff * numCorners + 2] = cone[5];
354: ++cOff;
355: break;
356: case 4: /* Quad -> 8 triangles */
357: Nt = 8;
358: cells[cOff * numCorners + 0] = cone[0];
359: cells[cOff * numCorners + 1] = cone[1];
360: cells[cOff * numCorners + 2] = cone[7];
361: ++cOff;
362: cells[cOff * numCorners + 0] = cone[1];
363: cells[cOff * numCorners + 1] = cone[2];
364: cells[cOff * numCorners + 2] = cone[3];
365: ++cOff;
366: cells[cOff * numCorners + 0] = cone[3];
367: cells[cOff * numCorners + 1] = cone[4];
368: cells[cOff * numCorners + 2] = cone[5];
369: ++cOff;
370: cells[cOff * numCorners + 0] = cone[5];
371: cells[cOff * numCorners + 1] = cone[6];
372: cells[cOff * numCorners + 2] = cone[7];
373: ++cOff;
374: cells[cOff * numCorners + 0] = cone[8];
375: cells[cOff * numCorners + 1] = cone[1];
376: cells[cOff * numCorners + 2] = cone[3];
377: ++cOff;
378: cells[cOff * numCorners + 0] = cone[8];
379: cells[cOff * numCorners + 1] = cone[3];
380: cells[cOff * numCorners + 2] = cone[5];
381: ++cOff;
382: cells[cOff * numCorners + 0] = cone[8];
383: cells[cOff * numCorners + 1] = cone[5];
384: cells[cOff * numCorners + 2] = cone[7];
385: ++cOff;
386: cells[cOff * numCorners + 0] = cone[8];
387: cells[cOff * numCorners + 1] = cone[7];
388: cells[cOff * numCorners + 2] = cone[1];
389: ++cOff;
390: break;
391: default:
392: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
393: }
394: if (debug) {
395: for (t = 0; t < Nt; ++t) {
396: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t));
397: for (c = 0; c < numCorners; ++c) {
398: if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
399: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c]));
400: }
401: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
402: }
403: }
404: }
405: EGlite_free(lobjs);
406: }
407: }
408: PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells);
409: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
410: PetscCall(PetscFree3(coords, cells, cone));
411: PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells));
412: PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices));
413: /* Embed EGADS model in DM */
414: {
415: PetscContainer modelObj, contextObj;
417: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
418: PetscCall(PetscContainerSetPointer(modelObj, model));
419: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADSLite Model", (PetscObject)modelObj));
420: PetscCall(PetscContainerDestroy(&modelObj));
422: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
423: PetscCall(PetscContainerSetPointer(contextObj, context));
424: PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSLiteDestroy_Private));
425: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADSLite Context", (PetscObject)contextObj));
426: PetscCall(PetscContainerDestroy(&contextObj));
427: }
428: /* Label points */
429: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
430: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
431: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
432: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
433: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
434: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
435: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
436: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
437: cOff = 0;
438: for (b = 0; b < nbodies; ++b) {
439: ego body = bodies[b];
440: int id, Nl, l;
442: PetscCall(EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
443: for (l = 0; l < Nl; ++l) {
444: ego loop = lobjs[l];
445: ego *fobjs;
446: int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;
448: lid = EGlite_indexBodyTopo(body, loop);
449: PetscCall(EGlite_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
450: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
451: fid = EGlite_indexBodyTopo(body, fobjs[0]);
452: EGlite_free(fobjs);
453: PetscCall(EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
454: for (e = 0; e < Ne; ++e) {
455: ego edge = objs[e];
456: int eid, Nv, v;
457: PetscInt points[3], support[2], numEdges, edgeNum;
458: const PetscInt *edges;
460: eid = EGlite_indexBodyTopo(body, edge);
461: PetscCall(EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
462: if (mtype == DEGENERATE) continue;
463: else ++Ner;
464: for (v = 0; v < Nv; ++v) {
465: ego vertex = nobjs[v];
467: id = EGlite_indexBodyTopo(body, vertex);
468: PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid));
469: points[v * 2] = numCells + id - 1;
470: }
471: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
472: points[1] = numCells + numVertices - newVertices + edgeNum;
474: PetscCall(DMLabelSetValue(edgeLabel, points[1], eid));
475: support[0] = points[0];
476: support[1] = points[1];
477: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
478: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges);
479: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
480: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
481: support[0] = points[1];
482: support[1] = points[2];
483: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
484: PetscCheck(numEdges == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%" PetscInt_FMT ", %" PetscInt_FMT ") should only bound 1 edge, not %" PetscInt_FMT, support[0], support[1], numEdges);
485: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
486: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
487: }
488: switch (Ner) {
489: case 2:
490: Nt = 2;
491: break;
492: case 3:
493: Nt = 4;
494: break;
495: case 4:
496: Nt = 8;
497: break;
498: default:
499: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
500: }
501: for (t = 0; t < Nt; ++t) {
502: PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b));
503: PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid));
504: }
505: cOff += Nt;
506: }
507: EGlite_free(lobjs);
508: }
509: PetscCall(PetscHMapIDestroy(&edgeMap));
510: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
511: for (c = cStart; c < cEnd; ++c) {
512: PetscInt *closure = NULL;
513: PetscInt clSize, cl, bval, fval;
515: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
516: PetscCall(DMLabelGetValue(bodyLabel, c, &bval));
517: PetscCall(DMLabelGetValue(faceLabel, c, &fval));
518: for (cl = 0; cl < clSize * 2; cl += 2) {
519: PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval));
520: PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval));
521: }
522: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
523: }
524: *newdm = dm;
525: PetscFunctionReturn(PETSC_SUCCESS);
526: }
528: static PetscErrorCode DMPlexEGADSLitePrintModel_Internal(ego model)
529: {
530: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
531: int oclass, mtype, *senses;
532: int Nb, b;
534: PetscFunctionBeginUser;
535: /* test bodyTopo functions */
536: PetscCall(EGlite_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
537: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb));
539: for (b = 0; b < Nb; ++b) {
540: ego body = bodies[b];
541: int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;
543: /* Output Basic Model Topology */
544: PetscCall(EGlite_getBodyTopos(body, NULL, SHELL, &Nsh, &objs));
545: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh));
546: EGlite_free(objs);
548: PetscCall(EGlite_getBodyTopos(body, NULL, FACE, &Nf, &objs));
549: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf));
550: EGlite_free(objs);
552: PetscCall(EGlite_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
553: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl));
555: PetscCall(EGlite_getBodyTopos(body, NULL, EDGE, &Ne, &objs));
556: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne));
557: EGlite_free(objs);
559: PetscCall(EGlite_getBodyTopos(body, NULL, NODE, &Nv, &objs));
560: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv));
561: EGlite_free(objs);
563: for (l = 0; l < Nl; ++l) {
564: ego loop = lobjs[l];
566: id = EGlite_indexBodyTopo(body, loop);
567: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id));
569: /* Get EDGE info which associated with the current LOOP */
570: PetscCall(EGlite_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
572: for (e = 0; e < Ne; ++e) {
573: ego edge = objs[e];
574: double range[4] = {0., 0., 0., 0.};
575: double point[3] = {0., 0., 0.};
576: double params[3] = {0., 0., 0.};
577: double result[18];
578: int peri;
580: PetscCall(EGlite_indexBodyTopo(body, edge));
581: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e));
583: PetscCall(EGlite_getRange(edge, range, &peri));
584: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]));
586: /* Get NODE info which associated with the current EDGE */
587: PetscCall(EGlite_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
588: if (mtype == DEGENERATE) {
589: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id));
590: } else {
591: params[0] = range[0];
592: PetscCall(EGlite_evaluate(edge, params, result));
593: PetscCall(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]));
594: params[0] = range[1];
595: PetscCall(EGlite_evaluate(edge, params, result));
596: PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]));
597: }
599: for (v = 0; v < Nv; ++v) {
600: ego vertex = nobjs[v];
601: double limits[4];
602: int dummy;
604: PetscCall(EGlite_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
605: PetscCall(EGlite_indexBodyTopo(body, vertex));
606: PetscCall(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id));
607: PetscCall(PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]));
609: point[0] = point[0] + limits[0];
610: point[1] = point[1] + limits[1];
611: point[2] = point[2] + limits[2];
612: }
613: }
614: }
615: EGlite_free(lobjs);
616: }
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
619: #endif
621: /*@C
622: DMPlexCreateEGADSLiteFromFile - Create a DMPlex mesh from an EGADSLite file.
624: Collective
626: Input Parameters:
627: + comm - The MPI communicator
628: - filename - The name of the EGADSLite file
630: Output Parameter:
631: . dm - The DM object representing the mesh
633: Level: beginner
635: .seealso: `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSFromFile()`
636: @*/
637: PetscErrorCode DMPlexCreateEGADSLiteFromFile(MPI_Comm comm, const char filename[], DM *dm)
638: {
639: PetscMPIInt rank;
640: #if defined(PETSC_HAVE_EGADS)
641: ego context = NULL, model = NULL;
642: #endif
643: PetscBool printModel = PETSC_FALSE;
645: PetscFunctionBegin;
646: PetscAssertPointer(filename, 2);
647: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL));
648: PetscCallMPI(MPI_Comm_rank(comm, &rank));
649: #if defined(PETSC_HAVE_EGADS)
650: if (rank == 0) {
651: PetscCall(EGlite_open(&context));
652: PetscCall(EGlite_loadModel(context, 0, filename, &model));
653: if (printModel) PetscCall(DMPlexEGADSLitePrintModel_Internal(model));
654: }
655: PetscCall(DMPlexCreateEGADSLite_Internal(comm, context, model, dm));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: #else
658: SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
659: #endif
660: }