Actual source code: plexegads.c
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/hashmapi.h>
4: #ifdef PETSC_HAVE_EGADS
5: #include <egads.h>
6: #endif
8: /* We need to understand how to natively parse STEP files. There seems to be only one open-source implementation of
9: the STEP parser contained in the OpenCASCADE package. It is enough to make a strong man weep:
11: https://github.com/tpaviot/oce/tree/master/src/STEPControl
13: The STEP, and inner EXPRESS, formats are ISO standards, so they are documented
15: https://stackoverflow.com/questions/26774037/documentation-or-specification-for-step-and-stp-files
16: http://stepmod.sourceforge.net/express_model_spec/
18: but again it seems that there has been a deliberate effort at obfuscation, probably to raise the bar for entrants.
19: */
21: #ifdef PETSC_HAVE_EGADS
22: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
23: PETSC_INTERN PetscErrorCode DMSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
25: PetscErrorCode DMSnapToGeomModel_EGADS_Internal(DM dm, PetscInt p, ego model, PetscInt bodyID, PetscInt faceID, PetscInt edgeID, const PetscScalar mcoords[], PetscScalar gcoords[])
26: {
27: DM cdm;
28: ego *bodies;
29: ego geom, body, obj;
30: /* result has to hold derivatives, along with the value */
31: double params[3], result[18], paramsV[16 * 3], resultV[16 * 3], range[4];
32: int Nb, oclass, mtype, *senses, peri;
33: Vec coordinatesLocal;
34: PetscScalar *coords = NULL;
35: PetscInt Nv, v, Np = 0, pm;
36: PetscInt dE, d;
38: PetscFunctionBeginHot;
39: PetscCall(DMGetCoordinateDM(dm, &cdm));
40: PetscCall(DMGetCoordinateDim(dm, &dE));
41: PetscCall(DMGetCoordinatesLocal(dm, &coordinatesLocal));
42: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
43: PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
44: body = bodies[bodyID];
46: if (edgeID >= 0) {
47: PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &obj));
48: Np = 1;
49: } else if (faceID >= 0) {
50: PetscCall(EG_objectBodyTopo(body, FACE, faceID, &obj));
51: Np = 2;
52: } else {
53: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: /* Calculate parameters (t or u,v) for vertices */
58: PetscCall(DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
59: Nv /= dE;
60: if (Nv == 1) {
61: PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
62: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
65: PetscCheck(Nv <= 16, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %" PetscInt_FMT " coordinates associated to point %" PetscInt_FMT, Nv, p);
67: /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */
68: PetscCall(EG_getRange(obj, range, &peri));
69: for (v = 0; v < Nv; ++v) {
70: PetscCall(EG_invEvaluate(obj, &coords[v * dE], ¶msV[v * 3], &resultV[v * 3]));
71: #if 1
72: if (peri > 0) {
73: if (paramsV[v * 3 + 0] + 1.e-4 < range[0]) {
74: paramsV[v * 3 + 0] += 2. * PETSC_PI;
75: } else if (paramsV[v * 3 + 0] - 1.e-4 > range[1]) {
76: paramsV[v * 3 + 0] -= 2. * PETSC_PI;
77: }
78: }
79: if (peri > 1) {
80: if (paramsV[v * 3 + 1] + 1.e-4 < range[2]) {
81: paramsV[v * 3 + 1] += 2. * PETSC_PI;
82: } else if (paramsV[v * 3 + 1] - 1.e-4 > range[3]) {
83: paramsV[v * 3 + 1] -= 2. * PETSC_PI;
84: }
85: }
86: #endif
87: }
88: PetscCall(DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords));
89: /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
90: for (pm = 0; pm < Np; ++pm) {
91: params[pm] = 0.;
92: for (v = 0; v < Nv; ++v) params[pm] += paramsV[v * 3 + pm];
93: params[pm] /= Nv;
94: }
95: PetscCheck(!(params[0] < range[0]) && !(params[0] > range[1]), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %" PetscInt_FMT " had bad interpolation", p);
96: 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);
97: /* Put coordinates for new vertex in result[] */
98: PetscCall(EG_evaluate(obj, params, result));
99: for (d = 0; d < dE; ++d) gcoords[d] = result[d];
100: PetscFunctionReturn(PETSC_SUCCESS);
101: }
102: #endif
104: PetscErrorCode DMSnapToGeomModel_EGADSLite(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
105: {
106: PetscFunctionBeginHot;
107: #ifdef PETSC_HAVE_EGADS
108: DMLabel bodyLabel, faceLabel, edgeLabel;
109: PetscInt bodyID, faceID, edgeID;
110: PetscContainer modelObj;
111: ego model;
113: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
114: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
115: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
116: PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite meshes must have body, face, and edge labels defined");
117: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj));
118: PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADSLite mesh missing model object");
120: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
121: PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID));
122: PetscCall(DMLabelGetValue(faceLabel, p, &faceID));
123: PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID));
124: /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
125: if (bodyID < 0) {
126: for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
129: PetscCall(DMSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
130: #endif
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: PetscErrorCode DMSnapToGeomModel_EGADS(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
135: {
136: PetscFunctionBeginHot;
137: #ifdef PETSC_HAVE_EGADS
138: DMLabel bodyLabel, faceLabel, edgeLabel;
139: PetscInt bodyID, faceID, edgeID;
140: PetscContainer modelObj;
141: ego model;
143: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
144: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
145: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
146: PetscCheck(bodyLabel && faceLabel && edgeLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS meshes must have body, face, and edge labels defined");
147: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
148: PetscCheck(modelObj, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "EGADS mesh missing model object");
150: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
151: PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID));
152: PetscCall(DMLabelGetValue(faceLabel, p, &faceID));
153: PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID));
154: /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
155: if (bodyID < 0) {
156: for (PetscInt d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
159: PetscCall(DMSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
160: #endif
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: #if defined(PETSC_HAVE_EGADS)
165: static PetscErrorCode DMPlexEGADSPrintModel_Internal(ego model)
166: {
167: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
168: int oclass, mtype, *senses;
169: int Nb, b;
171: PetscFunctionBeginUser;
172: /* test bodyTopo functions */
173: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
174: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb));
176: for (b = 0; b < Nb; ++b) {
177: ego body = bodies[b];
178: int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;
180: /* Output Basic Model Topology */
181: PetscCall(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs));
182: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh));
183: EG_free(objs);
185: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs));
186: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf));
187: EG_free(objs);
189: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
190: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl));
192: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs));
193: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne));
194: EG_free(objs);
196: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs));
197: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv));
198: EG_free(objs);
200: for (l = 0; l < Nl; ++l) {
201: ego loop = lobjs[l];
203: id = EG_indexBodyTopo(body, loop);
204: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id));
206: /* Get EDGE info which associated with the current LOOP */
207: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
209: for (e = 0; e < Ne; ++e) {
210: ego edge = objs[e];
211: double range[4] = {0., 0., 0., 0.};
212: double point[3] = {0., 0., 0.};
213: double params[3] = {0., 0., 0.};
214: double result[18];
215: int peri;
217: id = EG_indexBodyTopo(body, edge);
218: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e));
220: PetscCall(EG_getRange(edge, range, &peri));
221: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]));
223: /* Get NODE info which associated with the current EDGE */
224: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
225: if (mtype == DEGENERATE) {
226: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id));
227: } else {
228: params[0] = range[0];
229: PetscCall(EG_evaluate(edge, params, result));
230: PetscCall(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]));
231: params[0] = range[1];
232: PetscCall(EG_evaluate(edge, params, result));
233: PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]));
234: }
236: for (v = 0; v < Nv; ++v) {
237: ego vertex = nobjs[v];
238: double limits[4];
239: int dummy;
241: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
242: id = EG_indexBodyTopo(body, vertex);
243: PetscCall(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id));
244: PetscCall(PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]));
246: point[0] = point[0] + limits[0];
247: point[1] = point[1] + limits[1];
248: point[2] = point[2] + limits[2];
249: }
250: }
251: }
252: EG_free(lobjs);
253: }
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: static PetscErrorCode DMPlexEGADSDestroy_Private(void **context)
258: {
259: if (context) EG_close((ego)*context);
260: return 0;
261: }
263: static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
264: {
265: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
266: PetscInt cStart, cEnd, c;
267: /* EGADSLite variables */
268: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
269: int oclass, mtype, nbodies, *senses;
270: int b;
271: /* PETSc variables */
272: DM dm;
273: PetscHMapI edgeMap = NULL;
274: 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;
275: PetscInt *cells = NULL, *cone = NULL;
276: PetscReal *coords = NULL;
277: PetscMPIInt rank;
279: PetscFunctionBegin;
280: PetscCallMPI(MPI_Comm_rank(comm, &rank));
281: if (rank == 0) {
282: const PetscInt debug = 0;
284: /* ---------------------------------------------------------------------------------------------------
285: Generate Petsc Plex
286: Get all Nodes in model, record coordinates in a correctly formatted array
287: Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
288: We need to uniformly refine the initial geometry to guarantee a valid mesh
289: */
291: /* Calculate cell and vertex sizes */
292: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
293: PetscCall(PetscHMapICreate(&edgeMap));
294: numEdges = 0;
295: for (b = 0; b < nbodies; ++b) {
296: ego body = bodies[b];
297: int id, Nl, l, Nv, v;
299: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
300: for (l = 0; l < Nl; ++l) {
301: ego loop = lobjs[l];
302: int Ner = 0, Ne, e, Nc;
304: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
305: for (e = 0; e < Ne; ++e) {
306: ego edge = objs[e];
307: int Nv, id;
308: PetscHashIter iter;
309: PetscBool found;
311: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
312: if (mtype == DEGENERATE) continue;
313: id = EG_indexBodyTopo(body, edge);
314: PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found));
315: if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++));
316: ++Ner;
317: }
318: if (Ner == 2) {
319: Nc = 2;
320: } else if (Ner == 3) {
321: Nc = 4;
322: } else if (Ner == 4) {
323: Nc = 8;
324: ++numQuads;
325: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
326: numCells += Nc;
327: newCells += Nc - 1;
328: maxCorners = PetscMax(Ner * 2 + 1, maxCorners);
329: }
330: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
331: for (v = 0; v < Nv; ++v) {
332: ego vertex = nobjs[v];
334: id = EG_indexBodyTopo(body, vertex);
335: /* TODO: Instead of assuming contiguous ids, we could use a hash table */
336: numVertices = PetscMax(id, numVertices);
337: }
338: EG_free(lobjs);
339: EG_free(nobjs);
340: }
341: PetscCall(PetscHMapIGetSize(edgeMap, &numEdges));
342: newVertices = numEdges + numQuads;
343: numVertices += newVertices;
345: dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
346: cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
347: numCorners = 3; /* Split cells into triangles */
348: PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone));
350: /* Get vertex coordinates */
351: for (b = 0; b < nbodies; ++b) {
352: ego body = bodies[b];
353: int id, Nv, v;
355: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
356: for (v = 0; v < Nv; ++v) {
357: ego vertex = nobjs[v];
358: double limits[4];
359: int dummy;
361: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
362: id = EG_indexBodyTopo(body, vertex);
363: coords[(id - 1) * cdim + 0] = limits[0];
364: coords[(id - 1) * cdim + 1] = limits[1];
365: coords[(id - 1) * cdim + 2] = limits[2];
366: }
367: EG_free(nobjs);
368: }
369: PetscCall(PetscHMapIClear(edgeMap));
370: fOff = numVertices - newVertices + numEdges;
371: numEdges = 0;
372: numQuads = 0;
373: for (b = 0; b < nbodies; ++b) {
374: ego body = bodies[b];
375: int Nl, l;
377: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
378: for (l = 0; l < Nl; ++l) {
379: ego loop = lobjs[l];
380: int lid, Ner = 0, Ne, e;
382: lid = EG_indexBodyTopo(body, loop);
383: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
384: for (e = 0; e < Ne; ++e) {
385: ego edge = objs[e];
386: int eid, Nv;
387: PetscHashIter iter;
388: PetscBool found;
390: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
391: if (mtype == DEGENERATE) continue;
392: ++Ner;
393: eid = EG_indexBodyTopo(body, edge);
394: PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found));
395: if (!found) {
396: PetscInt v = numVertices - newVertices + numEdges;
397: double range[4], params[3] = {0., 0., 0.}, result[18];
398: int periodic[2];
400: PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++));
401: PetscCall(EG_getRange(edge, range, periodic));
402: params[0] = 0.5 * (range[0] + range[1]);
403: PetscCall(EG_evaluate(edge, params, result));
404: coords[v * cdim + 0] = result[0];
405: coords[v * cdim + 1] = result[1];
406: coords[v * cdim + 2] = result[2];
407: }
408: }
409: if (Ner == 4) {
410: PetscInt v = fOff + numQuads++;
411: ego *fobjs, face;
412: double range[4], params[3] = {0., 0., 0.}, result[18];
413: int Nf, fid, periodic[2];
415: PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
416: face = fobjs[0];
417: fid = EG_indexBodyTopo(body, face);
418: PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid);
419: PetscCall(EG_getRange(face, range, periodic));
420: params[0] = 0.5 * (range[0] + range[1]);
421: params[1] = 0.5 * (range[2] + range[3]);
422: PetscCall(EG_evaluate(face, params, result));
423: coords[v * cdim + 0] = result[0];
424: coords[v * cdim + 1] = result[1];
425: coords[v * cdim + 2] = result[2];
426: }
427: }
428: }
429: PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices);
431: /* Get cell vertices by traversing loops */
432: numQuads = 0;
433: cOff = 0;
434: for (b = 0; b < nbodies; ++b) {
435: ego body = bodies[b];
436: int id, Nl, l;
438: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
439: for (l = 0; l < Nl; ++l) {
440: ego loop = lobjs[l];
441: int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;
443: lid = EG_indexBodyTopo(body, loop);
444: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
446: for (e = 0; e < Ne; ++e) {
447: ego edge = objs[e];
448: int points[3];
449: int eid, Nv, v, tmp;
451: eid = EG_indexBodyTopo(body, edge);
452: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
453: if (mtype == DEGENERATE) continue;
454: else ++Ner;
455: PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);
457: for (v = 0; v < Nv; ++v) {
458: ego vertex = nobjs[v];
460: id = EG_indexBodyTopo(body, vertex);
461: points[v * 2] = id - 1;
462: }
463: {
464: PetscInt edgeNum;
466: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
467: points[1] = numVertices - newVertices + edgeNum;
468: }
469: /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
470: if (!nc) {
471: for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v];
472: } else {
473: if (cone[nc - 1] == points[0]) {
474: cone[nc++] = points[1];
475: if (cone[0] != points[2]) cone[nc++] = points[2];
476: } else if (cone[nc - 1] == points[2]) {
477: cone[nc++] = points[1];
478: if (cone[0] != points[0]) cone[nc++] = points[0];
479: } else if (cone[nc - 3] == points[0]) {
480: tmp = cone[nc - 3];
481: cone[nc - 3] = cone[nc - 1];
482: cone[nc - 1] = tmp;
483: cone[nc++] = points[1];
484: if (cone[0] != points[2]) cone[nc++] = points[2];
485: } else if (cone[nc - 3] == points[2]) {
486: tmp = cone[nc - 3];
487: cone[nc - 3] = cone[nc - 1];
488: cone[nc - 1] = tmp;
489: cone[nc++] = points[1];
490: if (cone[0] != points[0]) cone[nc++] = points[0];
491: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
492: }
493: }
494: PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner);
495: if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++;
496: PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners);
497: /* Triangulate the loop */
498: switch (Ner) {
499: case 2: /* Bi-Segment -> 2 triangles */
500: Nt = 2;
501: cells[cOff * numCorners + 0] = cone[0];
502: cells[cOff * numCorners + 1] = cone[1];
503: cells[cOff * numCorners + 2] = cone[2];
504: ++cOff;
505: cells[cOff * numCorners + 0] = cone[0];
506: cells[cOff * numCorners + 1] = cone[2];
507: cells[cOff * numCorners + 2] = cone[3];
508: ++cOff;
509: break;
510: case 3: /* Triangle -> 4 triangles */
511: Nt = 4;
512: cells[cOff * numCorners + 0] = cone[0];
513: cells[cOff * numCorners + 1] = cone[1];
514: cells[cOff * numCorners + 2] = cone[5];
515: ++cOff;
516: cells[cOff * numCorners + 0] = cone[1];
517: cells[cOff * numCorners + 1] = cone[2];
518: cells[cOff * numCorners + 2] = cone[3];
519: ++cOff;
520: cells[cOff * numCorners + 0] = cone[5];
521: cells[cOff * numCorners + 1] = cone[3];
522: cells[cOff * numCorners + 2] = cone[4];
523: ++cOff;
524: cells[cOff * numCorners + 0] = cone[1];
525: cells[cOff * numCorners + 1] = cone[3];
526: cells[cOff * numCorners + 2] = cone[5];
527: ++cOff;
528: break;
529: case 4: /* Quad -> 8 triangles */
530: Nt = 8;
531: cells[cOff * numCorners + 0] = cone[0];
532: cells[cOff * numCorners + 1] = cone[1];
533: cells[cOff * numCorners + 2] = cone[7];
534: ++cOff;
535: cells[cOff * numCorners + 0] = cone[1];
536: cells[cOff * numCorners + 1] = cone[2];
537: cells[cOff * numCorners + 2] = cone[3];
538: ++cOff;
539: cells[cOff * numCorners + 0] = cone[3];
540: cells[cOff * numCorners + 1] = cone[4];
541: cells[cOff * numCorners + 2] = cone[5];
542: ++cOff;
543: cells[cOff * numCorners + 0] = cone[5];
544: cells[cOff * numCorners + 1] = cone[6];
545: cells[cOff * numCorners + 2] = cone[7];
546: ++cOff;
547: cells[cOff * numCorners + 0] = cone[8];
548: cells[cOff * numCorners + 1] = cone[1];
549: cells[cOff * numCorners + 2] = cone[3];
550: ++cOff;
551: cells[cOff * numCorners + 0] = cone[8];
552: cells[cOff * numCorners + 1] = cone[3];
553: cells[cOff * numCorners + 2] = cone[5];
554: ++cOff;
555: cells[cOff * numCorners + 0] = cone[8];
556: cells[cOff * numCorners + 1] = cone[5];
557: cells[cOff * numCorners + 2] = cone[7];
558: ++cOff;
559: cells[cOff * numCorners + 0] = cone[8];
560: cells[cOff * numCorners + 1] = cone[7];
561: cells[cOff * numCorners + 2] = cone[1];
562: ++cOff;
563: break;
564: default:
565: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
566: }
567: if (debug) {
568: for (t = 0; t < Nt; ++t) {
569: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t));
570: for (c = 0; c < numCorners; ++c) {
571: if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
572: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c]));
573: }
574: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
575: }
576: }
577: }
578: EG_free(lobjs);
579: }
580: }
581: PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells);
582: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
583: PetscCall(PetscFree3(coords, cells, cone));
584: PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells));
585: PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices));
586: /* Embed EGADS model in DM */
587: {
588: PetscContainer modelObj, contextObj;
590: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
591: PetscCall(PetscContainerSetPointer(modelObj, model));
592: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
593: PetscCall(PetscContainerDestroy(&modelObj));
595: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
596: PetscCall(PetscContainerSetPointer(contextObj, context));
597: PetscCall(PetscContainerSetCtxDestroy(contextObj, DMPlexEGADSDestroy_Private));
598: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
599: PetscCall(PetscContainerDestroy(&contextObj));
600: }
601: /* Label points */
602: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
603: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
604: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
605: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
606: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
607: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
608: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
609: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
610: cOff = 0;
611: for (b = 0; b < nbodies; ++b) {
612: ego body = bodies[b];
613: int id, Nl, l;
615: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
616: for (l = 0; l < Nl; ++l) {
617: ego loop = lobjs[l];
618: ego *fobjs;
619: int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;
621: lid = EG_indexBodyTopo(body, loop);
622: PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
623: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
624: fid = EG_indexBodyTopo(body, fobjs[0]);
625: EG_free(fobjs);
626: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
627: for (e = 0; e < Ne; ++e) {
628: ego edge = objs[e];
629: int eid, Nv, v;
630: PetscInt points[3], support[2], numEdges, edgeNum;
631: const PetscInt *edges;
633: eid = EG_indexBodyTopo(body, edge);
634: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
635: if (mtype == DEGENERATE) continue;
636: else ++Ner;
637: for (v = 0; v < Nv; ++v) {
638: ego vertex = nobjs[v];
640: id = EG_indexBodyTopo(body, vertex);
641: PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid));
642: points[v * 2] = numCells + id - 1;
643: }
644: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
645: points[1] = numCells + numVertices - newVertices + edgeNum;
647: PetscCall(DMLabelSetValue(edgeLabel, points[1], eid));
648: support[0] = points[0];
649: support[1] = points[1];
650: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
651: 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);
652: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
653: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
654: support[0] = points[1];
655: support[1] = points[2];
656: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
657: 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);
658: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
659: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
660: }
661: switch (Ner) {
662: case 2:
663: Nt = 2;
664: break;
665: case 3:
666: Nt = 4;
667: break;
668: case 4:
669: Nt = 8;
670: break;
671: default:
672: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
673: }
674: for (t = 0; t < Nt; ++t) {
675: PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b));
676: PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid));
677: }
678: cOff += Nt;
679: }
680: EG_free(lobjs);
681: }
682: PetscCall(PetscHMapIDestroy(&edgeMap));
683: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
684: for (c = cStart; c < cEnd; ++c) {
685: PetscInt *closure = NULL;
686: PetscInt clSize, cl, bval, fval;
688: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
689: PetscCall(DMLabelGetValue(bodyLabel, c, &bval));
690: PetscCall(DMLabelGetValue(faceLabel, c, &fval));
691: for (cl = 0; cl < clSize * 2; cl += 2) {
692: PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval));
693: PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval));
694: }
695: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
696: }
697: *newdm = dm;
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
702: {
703: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
704: // EGADS/EGADSLite variables
705: ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs;
706: ego topRef, prev, next;
707: int oclass, mtype, nbodies, *senses, *lSenses, *eSenses;
708: int b;
709: // PETSc variables
710: DM dm;
711: PetscHMapI edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL;
712: PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0;
713: PetscInt cellCntr = 0, numPoints = 0;
714: PetscInt *cells = NULL;
715: const PetscInt *cone = NULL;
716: PetscReal *coords = NULL;
717: PetscMPIInt rank;
719: PetscFunctionBeginUser;
720: PetscCallMPI(MPI_Comm_rank(comm, &rank));
721: if (rank == 0) {
722: // ---------------------------------------------------------------------------------------------------
723: // Generate Petsc Plex
724: // Get all Nodes in model, record coordinates in a correctly formatted array
725: // Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
726: // We need to uniformly refine the initial geometry to guarantee a valid mesh
728: // Calculate cell and vertex sizes
729: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
731: PetscCall(PetscHMapICreate(&edgeMap));
732: PetscCall(PetscHMapICreate(&bodyIndexMap));
733: PetscCall(PetscHMapICreate(&bodyVertexMap));
734: PetscCall(PetscHMapICreate(&bodyEdgeMap));
735: PetscCall(PetscHMapICreate(&bodyEdgeGlobalMap));
736: PetscCall(PetscHMapICreate(&bodyFaceMap));
738: for (b = 0; b < nbodies; ++b) {
739: ego body = bodies[b];
740: int Nf, Ne, Nv;
741: PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter;
742: PetscBool BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound;
744: PetscCall(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound));
745: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
746: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
747: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
748: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
750: if (!BIfound) PetscCall(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices));
751: if (!BVfound) PetscCall(PetscHMapISet(bodyVertexMap, b, numVertices));
752: if (!BEfound) PetscCall(PetscHMapISet(bodyEdgeMap, b, numEdges));
753: if (!BEGfound) PetscCall(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr));
754: if (!BFfound) PetscCall(PetscHMapISet(bodyFaceMap, b, numFaces));
756: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
757: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
758: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
759: EG_free(fobjs);
760: EG_free(eobjs);
761: EG_free(nobjs);
763: // Remove DEGENERATE EDGES from Edge count
764: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
765: int Netemp = 0;
766: for (int e = 0; e < Ne; ++e) {
767: ego edge = eobjs[e];
768: int eid;
770: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
771: eid = EG_indexBodyTopo(body, edge);
773: PetscCall(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound));
774: if (mtype == DEGENERATE) {
775: if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1));
776: } else {
777: ++Netemp;
778: if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp));
779: }
780: }
781: EG_free(eobjs);
783: // Determine Number of Cells
784: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
785: for (int f = 0; f < Nf; ++f) {
786: ego face = fobjs[f];
787: int edgeTemp = 0;
789: PetscCall(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs));
790: for (int e = 0; e < Ne; ++e) {
791: ego edge = eobjs[e];
793: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
794: if (mtype != DEGENERATE) ++edgeTemp;
795: }
796: numCells += (2 * edgeTemp);
797: EG_free(eobjs);
798: }
799: EG_free(fobjs);
801: numFaces += Nf;
802: numEdges += Netemp;
803: numVertices += Nv;
804: edgeCntr += Ne;
805: }
807: // Set up basic DMPlex parameters
808: dim = 2; // Assumes 3D Models :: Need to handle 2D models in the future
809: cdim = 3; // Assumes 3D Models :: Need to update to handle 2D models in future
810: numCorners = 3; // Split Faces into triangles
811: numPoints = numVertices + numEdges + numFaces; // total number of coordinate points
813: PetscCall(PetscMalloc2(numPoints * cdim, &coords, numCells * numCorners, &cells));
815: // Get Vertex Coordinates and Set up Cells
816: for (b = 0; b < nbodies; ++b) {
817: ego body = bodies[b];
818: int Nf, Ne, Nv;
819: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
820: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
821: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
823: // Vertices on Current Body
824: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
826: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
827: PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
828: PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));
830: for (int v = 0; v < Nv; ++v) {
831: ego vertex = nobjs[v];
832: double limits[4];
833: int id, dummy;
835: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
836: id = EG_indexBodyTopo(body, vertex);
838: coords[(bodyVertexIndexStart + id - 1) * cdim + 0] = limits[0];
839: coords[(bodyVertexIndexStart + id - 1) * cdim + 1] = limits[1];
840: coords[(bodyVertexIndexStart + id - 1) * cdim + 2] = limits[2];
841: }
842: EG_free(nobjs);
844: // Edge Midpoint Vertices on Current Body
845: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
847: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
848: PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
849: PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));
851: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
852: PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
853: PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));
855: for (int e = 0; e < Ne; ++e) {
856: ego edge = eobjs[e];
857: double range[2], avgt[1], cntrPnt[9];
858: int eid, eOffset;
859: int periodic;
861: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
862: if (mtype == DEGENERATE) continue;
864: eid = EG_indexBodyTopo(body, edge);
866: // get relative offset from globalEdgeID Vector
867: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
868: PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1);
869: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
871: PetscCall(EG_getRange(edge, range, &periodic));
872: avgt[0] = (range[0] + range[1]) / 2.;
874: PetscCall(EG_evaluate(edge, avgt, cntrPnt));
875: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 0] = cntrPnt[0];
876: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 1] = cntrPnt[1];
877: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 2] = cntrPnt[2];
878: }
879: EG_free(eobjs);
881: // Face Midpoint Vertices on Current Body
882: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
884: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
885: PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
886: PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));
888: for (int f = 0; f < Nf; ++f) {
889: ego face = fobjs[f];
890: double range[4], avgUV[2], cntrPnt[18];
891: int peri, id;
893: id = EG_indexBodyTopo(body, face);
894: PetscCall(EG_getRange(face, range, &peri));
896: avgUV[0] = (range[0] + range[1]) / 2.;
897: avgUV[1] = (range[2] + range[3]) / 2.;
898: PetscCall(EG_evaluate(face, avgUV, cntrPnt));
900: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 0] = cntrPnt[0];
901: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 1] = cntrPnt[1];
902: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 2] = cntrPnt[2];
903: }
904: EG_free(fobjs);
906: // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity
907: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
908: for (int f = 0; f < Nf; ++f) {
909: ego face = fobjs[f];
910: int fID, midFaceID, midPntID, startID, endID, Nl;
912: fID = EG_indexBodyTopo(body, face);
913: midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1;
914: // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges.
915: // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are:
916: // 1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option.
917: // This will use a default EGADS tessellation as an initial surface mesh.
918: // 2) Create the initial surface mesh via a 2D mesher :: Currently not available (?future?)
919: // May I suggest the XXXX as a starting point?
921: PetscCall(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses));
923: PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Face has %d Loops. Can only handle Faces with 1 Loop. Please use --dm_plex_egads_with_tess = 1 Option", Nl);
924: for (int l = 0; l < Nl; ++l) {
925: ego loop = lobjs[l];
927: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
928: for (int e = 0; e < Ne; ++e) {
929: ego edge = eobjs[e];
930: int eid, eOffset;
932: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
933: eid = EG_indexBodyTopo(body, edge);
934: if (mtype == DEGENERATE) continue;
936: // get relative offset from globalEdgeID Vector
937: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
938: PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
939: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
941: midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1;
943: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
945: if (eSenses[e] > 0) {
946: startID = EG_indexBodyTopo(body, nobjs[0]);
947: endID = EG_indexBodyTopo(body, nobjs[1]);
948: } else {
949: startID = EG_indexBodyTopo(body, nobjs[1]);
950: endID = EG_indexBodyTopo(body, nobjs[0]);
951: }
953: // Define 2 Cells per Edge with correct orientation
954: cells[cellCntr * numCorners + 0] = midFaceID;
955: cells[cellCntr * numCorners + 1] = bodyVertexIndexStart + startID - 1;
956: cells[cellCntr * numCorners + 2] = midPntID;
958: cells[cellCntr * numCorners + 3] = midFaceID;
959: cells[cellCntr * numCorners + 4] = midPntID;
960: cells[cellCntr * numCorners + 5] = bodyVertexIndexStart + endID - 1;
962: cellCntr = cellCntr + 2;
963: }
964: }
965: }
966: EG_free(fobjs);
967: }
968: }
970: // Generate DMPlex
971: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
972: PetscCall(PetscFree2(coords, cells));
973: PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " \n", numCells));
974: PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " \n", numVertices));
976: // Embed EGADS model in DM
977: {
978: PetscContainer modelObj, contextObj;
980: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
981: PetscCall(PetscContainerSetPointer(modelObj, model));
982: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
983: PetscCall(PetscContainerDestroy(&modelObj));
985: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
986: PetscCall(PetscContainerSetPointer(contextObj, context));
987: PetscCall(PetscContainerSetCtxDestroy(contextObj, DMPlexEGADSDestroy_Private));
988: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
989: PetscCall(PetscContainerDestroy(&contextObj));
990: }
991: // Label points
992: PetscInt nStart, nEnd;
994: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
995: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
996: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
997: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
998: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
999: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1000: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
1001: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1003: PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));
1005: cellCntr = 0;
1006: for (b = 0; b < nbodies; ++b) {
1007: ego body = bodies[b];
1008: int Nv, Ne, Nf;
1009: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
1010: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
1011: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
1013: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
1014: PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
1015: PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));
1017: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
1018: PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
1019: PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));
1021: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
1022: PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
1023: PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));
1025: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
1026: PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
1027: PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));
1029: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1030: for (int f = 0; f < Nf; ++f) {
1031: ego face = fobjs[f];
1032: int fID, Nl;
1034: fID = EG_indexBodyTopo(body, face);
1036: PetscCall(EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs));
1037: for (int l = 0; l < Nl; ++l) {
1038: ego loop = lobjs[l];
1039: int lid;
1041: lid = EG_indexBodyTopo(body, loop);
1042: PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
1044: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
1045: for (int e = 0; e < Ne; ++e) {
1046: ego edge = eobjs[e];
1047: int eid, eOffset;
1049: // Skip DEGENERATE Edges
1050: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
1051: if (mtype == DEGENERATE) continue;
1052: eid = EG_indexBodyTopo(body, edge);
1054: // get relative offset from globalEdgeID Vector
1055: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
1056: PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
1057: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
1059: PetscCall(EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs));
1060: for (int v = 0; v < Nv; ++v) {
1061: ego vertex = nobjs[v];
1062: int vID;
1064: vID = EG_indexBodyTopo(body, vertex);
1065: PetscCall(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b));
1066: PetscCall(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID));
1067: }
1068: EG_free(nobjs);
1070: PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b));
1071: PetscCall(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid));
1073: // Define Cell faces
1074: for (int jj = 0; jj < 2; ++jj) {
1075: PetscCall(DMLabelSetValue(bodyLabel, cellCntr, b));
1076: PetscCall(DMLabelSetValue(faceLabel, cellCntr, fID));
1077: PetscCall(DMPlexGetCone(dm, cellCntr, &cone));
1079: PetscCall(DMLabelSetValue(bodyLabel, cone[0], b));
1080: PetscCall(DMLabelSetValue(faceLabel, cone[0], fID));
1082: PetscCall(DMLabelSetValue(bodyLabel, cone[1], b));
1083: PetscCall(DMLabelSetValue(edgeLabel, cone[1], eid));
1085: PetscCall(DMLabelSetValue(bodyLabel, cone[2], b));
1086: PetscCall(DMLabelSetValue(faceLabel, cone[2], fID));
1088: cellCntr = cellCntr + 1;
1089: }
1090: }
1091: }
1092: EG_free(lobjs);
1094: PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b));
1095: PetscCall(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID));
1096: }
1097: EG_free(fobjs);
1098: }
1100: PetscCall(PetscHMapIDestroy(&edgeMap));
1101: PetscCall(PetscHMapIDestroy(&bodyIndexMap));
1102: PetscCall(PetscHMapIDestroy(&bodyVertexMap));
1103: PetscCall(PetscHMapIDestroy(&bodyEdgeMap));
1104: PetscCall(PetscHMapIDestroy(&bodyEdgeGlobalMap));
1105: PetscCall(PetscHMapIDestroy(&bodyFaceMap));
1107: *newdm = dm;
1108: PetscFunctionReturn(PETSC_SUCCESS);
1109: }
1111: static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
1112: {
1113: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1114: /* EGADSLite variables */
1115: ego geom, *bodies, *fobjs;
1116: int b, oclass, mtype, nbodies, *senses;
1117: int totalNumTris = 0, totalNumPoints = 0;
1118: double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize;
1119: /* PETSc variables */
1120: DM dm;
1121: PetscHMapI pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL;
1122: PetscHMapI pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL;
1123: PetscInt dim = -1, cdim = -1, numCorners = 0, counter = 0;
1124: PetscInt *cells = NULL;
1125: const PetscInt *cone = NULL;
1126: PetscReal *coords = NULL;
1127: PetscMPIInt rank;
1129: PetscFunctionBeginUser;
1130: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1131: if (rank == 0) {
1132: // ---------------------------------------------------------------------------------------------------
1133: // Generate Petsc Plex from EGADSlite created Tessellation of geometry
1134: // ---------------------------------------------------------------------------------------------------
1136: // Calculate cell and vertex sizes
1137: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
1139: PetscCall(PetscHMapICreate(&pointIndexStartMap));
1140: PetscCall(PetscHMapICreate(&triIndexStartMap));
1141: PetscCall(PetscHMapICreate(&pTypeLabelMap));
1142: PetscCall(PetscHMapICreate(&pIndexLabelMap));
1143: PetscCall(PetscHMapICreate(&pBodyIndexLabelMap));
1144: PetscCall(PetscHMapICreate(&triFaceIDLabelMap));
1145: PetscCall(PetscHMapICreate(&triBodyIDLabelMap));
1147: /* Create Tessellation of Bodies */
1148: ego tessArray[nbodies];
1150: for (b = 0; b < nbodies; ++b) {
1151: ego body = bodies[b];
1152: double params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation
1153: int Nf, bodyNumPoints = 0, bodyNumTris = 0;
1154: PetscHashIter PISiter, TISiter;
1155: PetscBool PISfound, TISfound;
1157: /* Store Start Index for each Body's Point and Tris */
1158: PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
1159: PetscCall(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound));
1161: if (!PISfound) PetscCall(PetscHMapISet(pointIndexStartMap, b, totalNumPoints));
1162: if (!TISfound) PetscCall(PetscHMapISet(triIndexStartMap, b, totalNumTris));
1164: /* Calculate Tessellation parameters based on Bounding Box */
1165: /* Get Bounding Box Dimensions of the BODY */
1166: PetscCall(EG_getBoundingBox(body, boundBox));
1167: tessSize = boundBox[3] - boundBox[0];
1168: if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1];
1169: if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2];
1171: // TODO :: May want to give users tessellation parameter options //
1172: params[0] = 0.0250 * tessSize;
1173: params[1] = 0.0075 * tessSize;
1174: params[2] = 15.0;
1176: PetscCall(EG_makeTessBody(body, params, &tessArray[b]));
1178: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1180: for (int f = 0; f < Nf; ++f) {
1181: ego face = fobjs[f];
1182: int len, fID, ntris;
1183: const int *ptype, *pindex, *ptris, *ptric;
1184: const double *pxyz, *puv;
1186: // Get Face ID //
1187: fID = EG_indexBodyTopo(body, face);
1189: // Checkout the Surface Tessellation //
1190: PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));
1192: // Determine total number of triangle cells in the tessellation //
1193: bodyNumTris += (int)ntris;
1195: // Check out the point index and coordinate //
1196: for (int p = 0; p < len; ++p) {
1197: int global;
1199: PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));
1201: // Determine the total number of points in the tessellation //
1202: bodyNumPoints = PetscMax(bodyNumPoints, global);
1203: }
1204: }
1205: EG_free(fobjs);
1207: totalNumPoints += bodyNumPoints;
1208: totalNumTris += bodyNumTris;
1209: }
1210: //} - Original End of (rank == 0)
1212: dim = 2;
1213: cdim = 3;
1214: numCorners = 3;
1215: //PetscInt counter = 0;
1217: /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA */
1218: /* Fill in below and use to define DMLabels after DMPlex creation */
1219: PetscCall(PetscMalloc2(totalNumPoints * cdim, &coords, totalNumTris * numCorners, &cells));
1221: for (b = 0; b < nbodies; ++b) {
1222: ego body = bodies[b];
1223: int Nf;
1224: PetscInt pointIndexStart;
1225: PetscHashIter PISiter;
1226: PetscBool PISfound;
1228: PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
1229: PetscCheck(PISfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b);
1230: PetscCall(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart));
1232: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1234: for (int f = 0; f < Nf; ++f) {
1235: /* Get Face Object */
1236: ego face = fobjs[f];
1237: int len, fID, ntris;
1238: const int *ptype, *pindex, *ptris, *ptric;
1239: const double *pxyz, *puv;
1241: /* Get Face ID */
1242: fID = EG_indexBodyTopo(body, face);
1244: /* Checkout the Surface Tessellation */
1245: PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));
1247: /* Check out the point index and coordinate */
1248: for (int p = 0; p < len; ++p) {
1249: int global;
1250: PetscHashIter PTLiter, PILiter, PBLiter;
1251: PetscBool PTLfound, PILfound, PBLfound;
1253: PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));
1255: /* Set the coordinates array for DAG */
1256: coords[((global - 1 + pointIndexStart) * 3) + 0] = pxyz[(p * 3) + 0];
1257: coords[((global - 1 + pointIndexStart) * 3) + 1] = pxyz[(p * 3) + 1];
1258: coords[((global - 1 + pointIndexStart) * 3) + 2] = pxyz[(p * 3) + 2];
1260: /* Store Geometry Label Information for DMLabel assignment later */
1261: PetscCall(PetscHMapIFind(pTypeLabelMap, global - 1 + pointIndexStart, &PTLiter, &PTLfound));
1262: PetscCall(PetscHMapIFind(pIndexLabelMap, global - 1 + pointIndexStart, &PILiter, &PILfound));
1263: PetscCall(PetscHMapIFind(pBodyIndexLabelMap, global - 1 + pointIndexStart, &PBLiter, &PBLfound));
1265: if (!PTLfound) PetscCall(PetscHMapISet(pTypeLabelMap, global - 1 + pointIndexStart, ptype[p]));
1266: if (!PILfound) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, pindex[p]));
1267: if (!PBLfound) PetscCall(PetscHMapISet(pBodyIndexLabelMap, global - 1 + pointIndexStart, b));
1269: if (ptype[p] < 0) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, fID));
1270: }
1272: for (int t = 0; t < (int)ntris; ++t) {
1273: int global, globalA, globalB;
1274: PetscHashIter TFLiter, TBLiter;
1275: PetscBool TFLfound, TBLfound;
1277: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 0], &global));
1278: cells[(counter * 3) + 0] = global - 1 + pointIndexStart;
1280: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 1], &globalA));
1281: cells[(counter * 3) + 1] = globalA - 1 + pointIndexStart;
1283: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 2], &globalB));
1284: cells[(counter * 3) + 2] = globalB - 1 + pointIndexStart;
1286: PetscCall(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound));
1287: PetscCall(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound));
1289: if (!TFLfound) PetscCall(PetscHMapISet(triFaceIDLabelMap, counter, fID));
1290: if (!TBLfound) PetscCall(PetscHMapISet(triBodyIDLabelMap, counter, b));
1292: counter += 1;
1293: }
1294: }
1295: EG_free(fobjs);
1296: }
1297: }
1299: //Build DMPlex
1300: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
1301: PetscCall(PetscFree2(coords, cells));
1303: // Embed EGADS model in DM
1304: {
1305: PetscContainer modelObj, contextObj;
1307: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
1308: PetscCall(PetscContainerSetPointer(modelObj, model));
1309: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
1310: PetscCall(PetscContainerDestroy(&modelObj));
1312: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
1313: PetscCall(PetscContainerSetPointer(contextObj, context));
1314: PetscCall(PetscContainerSetCtxDestroy(contextObj, DMPlexEGADSDestroy_Private));
1315: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
1316: PetscCall(PetscContainerDestroy(&contextObj));
1317: }
1319: // Label Points
1320: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
1321: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
1322: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
1323: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
1324: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
1325: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1326: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
1327: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1329: /* Get Number of DAG Nodes at each level */
1330: int fStart, fEnd, eStart, eEnd, nStart, nEnd;
1332: PetscCall(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd));
1333: PetscCall(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd));
1334: PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));
1336: /* Set DMLabels for NODES */
1337: for (int n = nStart; n < nEnd; ++n) {
1338: int pTypeVal, pIndexVal, pBodyVal;
1339: PetscHashIter PTLiter, PILiter, PBLiter;
1340: PetscBool PTLfound, PILfound, PBLfound;
1342: //Converted to Hash Tables
1343: PetscCall(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound));
1344: PetscCheck(PTLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n);
1345: PetscCall(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal));
1347: PetscCall(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound));
1348: PetscCheck(PILfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n);
1349: PetscCall(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal));
1351: PetscCall(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound));
1352: PetscCheck(PBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n);
1353: PetscCall(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal));
1355: PetscCall(DMLabelSetValue(bodyLabel, n, pBodyVal));
1356: if (pTypeVal == 0) PetscCall(DMLabelSetValue(vertexLabel, n, pIndexVal));
1357: if (pTypeVal > 0) PetscCall(DMLabelSetValue(edgeLabel, n, pIndexVal));
1358: if (pTypeVal < 0) PetscCall(DMLabelSetValue(faceLabel, n, pIndexVal));
1359: }
1361: /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */
1362: for (int e = eStart; e < eEnd; ++e) {
1363: int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1;
1365: PetscCall(DMPlexGetCone(dm, e, &cone));
1366: PetscCall(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end?
1367: PetscCall(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0));
1368: PetscCall(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1));
1369: PetscCall(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0));
1370: PetscCall(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1));
1371: PetscCall(DMLabelGetValue(faceLabel, cone[0], &faceID_0));
1372: PetscCall(DMLabelGetValue(faceLabel, cone[1], &faceID_1));
1374: PetscCall(DMLabelSetValue(bodyLabel, e, bodyID_0));
1376: if (edgeID_0 == edgeID_1) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
1377: else if (vertexID_0 > 0 && edgeID_1 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_1));
1378: else if (vertexID_1 > 0 && edgeID_0 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
1379: else { /* Do Nothing */ }
1380: }
1382: /* Set DMLabels for Cells */
1383: for (int f = fStart; f < fEnd; ++f) {
1384: int edgeID_0;
1385: PetscInt triBodyVal, triFaceVal;
1386: PetscHashIter TFLiter, TBLiter;
1387: PetscBool TFLfound, TBLfound;
1389: // Convert to Hash Table
1390: PetscCall(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound));
1391: PetscCheck(TFLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f);
1392: PetscCall(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal));
1394: PetscCall(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound));
1395: PetscCheck(TBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f);
1396: PetscCall(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal));
1398: PetscCall(DMLabelSetValue(bodyLabel, f, triBodyVal));
1399: PetscCall(DMLabelSetValue(faceLabel, f, triFaceVal));
1401: /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */
1402: PetscCall(DMPlexGetCone(dm, f, &cone));
1404: for (int jj = 0; jj < 3; ++jj) {
1405: PetscCall(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0));
1407: if (edgeID_0 < 0) {
1408: PetscCall(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal));
1409: PetscCall(DMLabelSetValue(faceLabel, cone[jj], triFaceVal));
1410: }
1411: }
1412: }
1414: *newdm = dm;
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1417: #endif
1419: /*@
1420: DMPlexInflateToGeomModel - Snaps the vertex coordinates of a `DMPLEX` object representing the mesh to its geometry if some vertices depart from the model.
1421: This usually happens with non-conforming refinement.
1423: Collective
1425: Input Parameter:
1426: . dm - The uninflated `DM` object representing the mesh
1428: Level: intermediate
1430: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`
1431: @*/
1432: PetscErrorCode DMPlexInflateToGeomModel(DM dm)
1433: {
1434: #if defined(PETSC_HAVE_EGADS)
1435: /* EGADS Variables */
1436: ego model, geom, body, face, edge;
1437: ego *bodies;
1438: int Nb, oclass, mtype, *senses;
1439: double result[3];
1440: /* PETSc Variables */
1441: DM cdm;
1442: PetscContainer modelObj;
1443: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1444: Vec coordinates;
1445: PetscScalar *coords;
1446: PetscInt bodyID, faceID, edgeID, vertexID;
1447: PetscInt cdim, d, vStart, vEnd, v;
1448: #endif
1450: PetscFunctionBegin;
1451: #if defined(PETSC_HAVE_EGADS)
1452: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
1453: if (!modelObj) PetscFunctionReturn(PETSC_SUCCESS);
1454: PetscCall(DMGetCoordinateDim(dm, &cdim));
1455: PetscCall(DMGetCoordinateDM(dm, &cdm));
1456: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1457: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
1458: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
1459: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1460: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1462: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
1463: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
1465: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1466: PetscCall(VecGetArrayWrite(coordinates, &coords));
1467: for (v = vStart; v < vEnd; ++v) {
1468: PetscScalar *vcoords;
1470: PetscCall(DMLabelGetValue(bodyLabel, v, &bodyID));
1471: PetscCall(DMLabelGetValue(faceLabel, v, &faceID));
1472: PetscCall(DMLabelGetValue(edgeLabel, v, &edgeID));
1473: PetscCall(DMLabelGetValue(vertexLabel, v, &vertexID));
1475: PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
1476: body = bodies[bodyID];
1478: PetscCall(DMPlexPointLocalRef(cdm, v, coords, (void *)&vcoords));
1479: if (edgeID > 0) {
1480: /* Snap to EDGE at nearest location */
1481: double params[1];
1482: PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &edge));
1483: PetscCall(EG_invEvaluate(edge, vcoords, params, result)); // Get (x,y,z) of nearest point on EDGE
1484: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1485: } else if (faceID > 0) {
1486: /* Snap to FACE at nearest location */
1487: double params[2];
1488: PetscCall(EG_objectBodyTopo(body, FACE, faceID, &face));
1489: PetscCall(EG_invEvaluate(face, vcoords, params, result)); // Get (x,y,z) of nearest point on FACE
1490: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1491: }
1492: }
1493: PetscCall(VecRestoreArrayWrite(coordinates, &coords));
1494: /* Clear out global coordinates */
1495: PetscCall(VecDestroy(&dm->coordinates[0].x));
1496: #endif
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: /*@
1501: DMPlexCreateEGADSFromFile - Create a `DMPLEX` mesh from an EGADS, IGES, or STEP file.
1503: Collective
1505: Input Parameters:
1506: + comm - The MPI communicator
1507: - filename - The name of the EGADS, IGES, or STEP file
1509: Output Parameter:
1510: . dm - The `DM` object representing the mesh
1512: Level: beginner
1514: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSLiteFromFile()`
1515: @*/
1516: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
1517: {
1518: PetscMPIInt rank;
1519: #if defined(PETSC_HAVE_EGADS)
1520: ego context = NULL, model = NULL;
1521: #endif
1522: PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE;
1524: PetscFunctionBegin;
1525: PetscAssertPointer(filename, 2);
1526: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL));
1527: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL));
1528: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL));
1529: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1530: #if defined(PETSC_HAVE_EGADS)
1531: if (rank == 0) {
1532: PetscCall(EG_open(&context));
1533: PetscCall(EG_loadModel(context, 0, filename, &model));
1534: if (printModel) PetscCall(DMPlexEGADSPrintModel_Internal(model));
1535: }
1536: if (tessModel) PetscCall(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm));
1537: else if (newModel) PetscCall(DMPlexCreateEGADS(comm, context, model, dm));
1538: else PetscCall(DMPlexCreateEGADS_Internal(comm, context, model, dm));
1539: PetscFunctionReturn(PETSC_SUCCESS);
1540: #else
1541: SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads");
1542: #endif
1543: }