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 DMPlexSnapToGeomModel_EGADS_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
23: PETSC_INTERN PetscErrorCode DMPlexSnapToGeomModel_EGADSLite_Internal(DM, PetscInt, ego, PetscInt, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
25: PetscErrorCode DMPlexSnapToGeomModel_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: /*@
105: DMPlexSnapToGeomModel - Given a coordinate point 'mcoords' on the mesh point 'p', return the closest coordinate point 'gcoords' on the geometry model associated with that point.
107: Not Collective
109: Input Parameters:
110: + dm - The `DMPLEX` object
111: . p - The mesh point
112: . dE - The coordinate dimension
113: - mcoords - A coordinate point lying on the mesh point
115: Output Parameter:
116: . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point
118: Level: intermediate
120: Note:
121: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS.
123: The coordinate dimension may be different from the coordinate dimension of the `dm`, for example if the transformation is extrusion.
125: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexCreate()`, `DMPlexSetRefinementUniform()`
126: @*/
127: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, PetscInt dE, const PetscScalar mcoords[], PetscScalar gcoords[])
128: {
129: PetscInt d;
131: PetscFunctionBeginHot;
132: #ifdef PETSC_HAVE_EGADS
133: {
134: DM_Plex *plex = (DM_Plex *)dm->data;
135: DMLabel bodyLabel, faceLabel, edgeLabel;
136: PetscInt bodyID, faceID, edgeID;
137: PetscContainer modelObj;
138: ego model;
139: PetscBool islite = PETSC_FALSE;
141: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
142: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
143: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
144: if (!bodyLabel || !faceLabel || !edgeLabel || plex->ignoreModel) {
145: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
148: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
149: if (!modelObj) {
150: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSLite Model", (PetscObject *)&modelObj));
151: islite = PETSC_TRUE;
152: }
153: if (!modelObj) {
154: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
157: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
158: PetscCall(DMLabelGetValue(bodyLabel, p, &bodyID));
159: PetscCall(DMLabelGetValue(faceLabel, p, &faceID));
160: PetscCall(DMLabelGetValue(edgeLabel, p, &edgeID));
161: /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
162: if (bodyID < 0) {
163: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
164: PetscFunctionReturn(PETSC_SUCCESS);
165: }
166: if (islite) PetscCall(DMPlexSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
167: else PetscCall(DMPlexSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords));
168: }
169: #else
170: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
171: #endif
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: #if defined(PETSC_HAVE_EGADS)
176: static PetscErrorCode DMPlexEGADSPrintModel_Internal(ego model)
177: {
178: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
179: int oclass, mtype, *senses;
180: int Nb, b;
182: PetscFunctionBeginUser;
183: /* test bodyTopo functions */
184: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
185: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb));
187: for (b = 0; b < Nb; ++b) {
188: ego body = bodies[b];
189: int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;
191: /* Output Basic Model Topology */
192: PetscCall(EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs));
193: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh));
194: EG_free(objs);
196: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &objs));
197: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf));
198: EG_free(objs);
200: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
201: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl));
203: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs));
204: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne));
205: EG_free(objs);
207: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &objs));
208: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv));
209: EG_free(objs);
211: for (l = 0; l < Nl; ++l) {
212: ego loop = lobjs[l];
214: id = EG_indexBodyTopo(body, loop);
215: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id));
217: /* Get EDGE info which associated with the current LOOP */
218: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
220: for (e = 0; e < Ne; ++e) {
221: ego edge = objs[e];
222: double range[4] = {0., 0., 0., 0.};
223: double point[3] = {0., 0., 0.};
224: double params[3] = {0., 0., 0.};
225: double result[18];
226: int peri;
228: id = EG_indexBodyTopo(body, edge);
229: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e));
231: PetscCall(EG_getRange(edge, range, &peri));
232: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]));
234: /* Get NODE info which associated with the current EDGE */
235: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
236: if (mtype == DEGENERATE) {
237: PetscCall(PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id));
238: } else {
239: params[0] = range[0];
240: PetscCall(EG_evaluate(edge, params, result));
241: PetscCall(PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]));
242: params[0] = range[1];
243: PetscCall(EG_evaluate(edge, params, result));
244: PetscCall(PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]));
245: }
247: for (v = 0; v < Nv; ++v) {
248: ego vertex = nobjs[v];
249: double limits[4];
250: int dummy;
252: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
253: id = EG_indexBodyTopo(body, vertex);
254: PetscCall(PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id));
255: PetscCall(PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]));
257: point[0] = point[0] + limits[0];
258: point[1] = point[1] + limits[1];
259: point[2] = point[2] + limits[2];
260: }
261: }
262: }
263: EG_free(lobjs);
264: }
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: static PetscErrorCode DMPlexEGADSDestroy_Private(void *context)
269: {
270: if (context) EG_close((ego)context);
271: return 0;
272: }
274: static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
275: {
276: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
277: PetscInt cStart, cEnd, c;
278: /* EGADSLite variables */
279: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
280: int oclass, mtype, nbodies, *senses;
281: int b;
282: /* PETSc variables */
283: DM dm;
284: PetscHMapI edgeMap = NULL;
285: 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;
286: PetscInt *cells = NULL, *cone = NULL;
287: PetscReal *coords = NULL;
288: PetscMPIInt rank;
290: PetscFunctionBegin;
291: PetscCallMPI(MPI_Comm_rank(comm, &rank));
292: if (rank == 0) {
293: const PetscInt debug = 0;
295: /* ---------------------------------------------------------------------------------------------------
296: Generate Petsc Plex
297: Get all Nodes in model, record coordinates in a correctly formatted array
298: Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
299: We need to uniformly refine the initial geometry to guarantee a valid mesh
300: */
302: /* Calculate cell and vertex sizes */
303: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
304: PetscCall(PetscHMapICreate(&edgeMap));
305: numEdges = 0;
306: for (b = 0; b < nbodies; ++b) {
307: ego body = bodies[b];
308: int id, Nl, l, Nv, v;
310: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
311: for (l = 0; l < Nl; ++l) {
312: ego loop = lobjs[l];
313: int Ner = 0, Ne, e, Nc;
315: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
316: for (e = 0; e < Ne; ++e) {
317: ego edge = objs[e];
318: int Nv, id;
319: PetscHashIter iter;
320: PetscBool found;
322: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
323: if (mtype == DEGENERATE) continue;
324: id = EG_indexBodyTopo(body, edge);
325: PetscCall(PetscHMapIFind(edgeMap, id - 1, &iter, &found));
326: if (!found) PetscCall(PetscHMapISet(edgeMap, id - 1, numEdges++));
327: ++Ner;
328: }
329: if (Ner == 2) {
330: Nc = 2;
331: } else if (Ner == 3) {
332: Nc = 4;
333: } else if (Ner == 4) {
334: Nc = 8;
335: ++numQuads;
336: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
337: numCells += Nc;
338: newCells += Nc - 1;
339: maxCorners = PetscMax(Ner * 2 + 1, maxCorners);
340: }
341: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
342: for (v = 0; v < Nv; ++v) {
343: ego vertex = nobjs[v];
345: id = EG_indexBodyTopo(body, vertex);
346: /* TODO: Instead of assuming contiguous ids, we could use a hash table */
347: numVertices = PetscMax(id, numVertices);
348: }
349: EG_free(lobjs);
350: EG_free(nobjs);
351: }
352: PetscCall(PetscHMapIGetSize(edgeMap, &numEdges));
353: newVertices = numEdges + numQuads;
354: numVertices += newVertices;
356: dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
357: cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
358: numCorners = 3; /* Split cells into triangles */
359: PetscCall(PetscMalloc3(numVertices * cdim, &coords, numCells * numCorners, &cells, maxCorners, &cone));
361: /* Get vertex coordinates */
362: for (b = 0; b < nbodies; ++b) {
363: ego body = bodies[b];
364: int id, Nv, v;
366: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
367: for (v = 0; v < Nv; ++v) {
368: ego vertex = nobjs[v];
369: double limits[4];
370: int dummy;
372: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
373: id = EG_indexBodyTopo(body, vertex);
374: coords[(id - 1) * cdim + 0] = limits[0];
375: coords[(id - 1) * cdim + 1] = limits[1];
376: coords[(id - 1) * cdim + 2] = limits[2];
377: }
378: EG_free(nobjs);
379: }
380: PetscCall(PetscHMapIClear(edgeMap));
381: fOff = numVertices - newVertices + numEdges;
382: numEdges = 0;
383: numQuads = 0;
384: for (b = 0; b < nbodies; ++b) {
385: ego body = bodies[b];
386: int Nl, l;
388: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
389: for (l = 0; l < Nl; ++l) {
390: ego loop = lobjs[l];
391: int lid, Ner = 0, Ne, e;
393: lid = EG_indexBodyTopo(body, loop);
394: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
395: for (e = 0; e < Ne; ++e) {
396: ego edge = objs[e];
397: int eid, Nv;
398: PetscHashIter iter;
399: PetscBool found;
401: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
402: if (mtype == DEGENERATE) continue;
403: ++Ner;
404: eid = EG_indexBodyTopo(body, edge);
405: PetscCall(PetscHMapIFind(edgeMap, eid - 1, &iter, &found));
406: if (!found) {
407: PetscInt v = numVertices - newVertices + numEdges;
408: double range[4], params[3] = {0., 0., 0.}, result[18];
409: int periodic[2];
411: PetscCall(PetscHMapISet(edgeMap, eid - 1, numEdges++));
412: PetscCall(EG_getRange(edge, range, periodic));
413: params[0] = 0.5 * (range[0] + range[1]);
414: PetscCall(EG_evaluate(edge, params, result));
415: coords[v * cdim + 0] = result[0];
416: coords[v * cdim + 1] = result[1];
417: coords[v * cdim + 2] = result[2];
418: }
419: }
420: if (Ner == 4) {
421: PetscInt v = fOff + numQuads++;
422: ego *fobjs, face;
423: double range[4], params[3] = {0., 0., 0.}, result[18];
424: int Nf, fid, periodic[2];
426: PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
427: face = fobjs[0];
428: fid = EG_indexBodyTopo(body, face);
429: PetscCheck(Nf == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid - 1, Nf, fid);
430: PetscCall(EG_getRange(face, range, periodic));
431: params[0] = 0.5 * (range[0] + range[1]);
432: params[1] = 0.5 * (range[2] + range[3]);
433: PetscCall(EG_evaluate(face, params, result));
434: coords[v * cdim + 0] = result[0];
435: coords[v * cdim + 1] = result[1];
436: coords[v * cdim + 2] = result[2];
437: }
438: }
439: }
440: PetscCheck(numEdges + numQuads == newVertices, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %" PetscInt_FMT " != %" PetscInt_FMT " previous count", numEdges + numQuads, newVertices);
442: /* Get cell vertices by traversing loops */
443: numQuads = 0;
444: cOff = 0;
445: for (b = 0; b < nbodies; ++b) {
446: ego body = bodies[b];
447: int id, Nl, l;
449: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
450: for (l = 0; l < Nl; ++l) {
451: ego loop = lobjs[l];
452: int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;
454: lid = EG_indexBodyTopo(body, loop);
455: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
457: for (e = 0; e < Ne; ++e) {
458: ego edge = objs[e];
459: int points[3];
460: int eid, Nv, v, tmp;
462: eid = EG_indexBodyTopo(body, edge);
463: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
464: if (mtype == DEGENERATE) continue;
465: else ++Ner;
466: PetscCheck(Nv == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);
468: for (v = 0; v < Nv; ++v) {
469: ego vertex = nobjs[v];
471: id = EG_indexBodyTopo(body, vertex);
472: points[v * 2] = id - 1;
473: }
474: {
475: PetscInt edgeNum;
477: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
478: points[1] = numVertices - newVertices + edgeNum;
479: }
480: /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
481: if (!nc) {
482: for (v = 0; v < Nv + 1; ++v) cone[nc++] = points[v];
483: } else {
484: if (cone[nc - 1] == points[0]) {
485: cone[nc++] = points[1];
486: if (cone[0] != points[2]) cone[nc++] = points[2];
487: } else if (cone[nc - 1] == points[2]) {
488: cone[nc++] = points[1];
489: if (cone[0] != points[0]) cone[nc++] = points[0];
490: } else if (cone[nc - 3] == points[0]) {
491: tmp = cone[nc - 3];
492: cone[nc - 3] = cone[nc - 1];
493: cone[nc - 1] = tmp;
494: cone[nc++] = points[1];
495: if (cone[0] != points[2]) cone[nc++] = points[2];
496: } else if (cone[nc - 3] == points[2]) {
497: tmp = cone[nc - 3];
498: cone[nc - 3] = cone[nc - 1];
499: cone[nc - 1] = tmp;
500: cone[nc++] = points[1];
501: if (cone[0] != points[0]) cone[nc++] = points[0];
502: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
503: }
504: }
505: PetscCheck(nc == 2 * Ner, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " != %" PetscInt_FMT, nc, 2 * Ner);
506: if (Ner == 4) cone[nc++] = numVertices - newVertices + numEdges + numQuads++;
507: PetscCheck(nc <= maxCorners, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %" PetscInt_FMT " > %" PetscInt_FMT " max", nc, maxCorners);
508: /* Triangulate the loop */
509: switch (Ner) {
510: case 2: /* Bi-Segment -> 2 triangles */
511: Nt = 2;
512: cells[cOff * numCorners + 0] = cone[0];
513: cells[cOff * numCorners + 1] = cone[1];
514: cells[cOff * numCorners + 2] = cone[2];
515: ++cOff;
516: cells[cOff * numCorners + 0] = cone[0];
517: cells[cOff * numCorners + 1] = cone[2];
518: cells[cOff * numCorners + 2] = cone[3];
519: ++cOff;
520: break;
521: case 3: /* Triangle -> 4 triangles */
522: Nt = 4;
523: cells[cOff * numCorners + 0] = cone[0];
524: cells[cOff * numCorners + 1] = cone[1];
525: cells[cOff * numCorners + 2] = cone[5];
526: ++cOff;
527: cells[cOff * numCorners + 0] = cone[1];
528: cells[cOff * numCorners + 1] = cone[2];
529: cells[cOff * numCorners + 2] = cone[3];
530: ++cOff;
531: cells[cOff * numCorners + 0] = cone[5];
532: cells[cOff * numCorners + 1] = cone[3];
533: cells[cOff * numCorners + 2] = cone[4];
534: ++cOff;
535: cells[cOff * numCorners + 0] = cone[1];
536: cells[cOff * numCorners + 1] = cone[3];
537: cells[cOff * numCorners + 2] = cone[5];
538: ++cOff;
539: break;
540: case 4: /* Quad -> 8 triangles */
541: Nt = 8;
542: cells[cOff * numCorners + 0] = cone[0];
543: cells[cOff * numCorners + 1] = cone[1];
544: cells[cOff * numCorners + 2] = cone[7];
545: ++cOff;
546: cells[cOff * numCorners + 0] = cone[1];
547: cells[cOff * numCorners + 1] = cone[2];
548: cells[cOff * numCorners + 2] = cone[3];
549: ++cOff;
550: cells[cOff * numCorners + 0] = cone[3];
551: cells[cOff * numCorners + 1] = cone[4];
552: cells[cOff * numCorners + 2] = cone[5];
553: ++cOff;
554: cells[cOff * numCorners + 0] = cone[5];
555: cells[cOff * numCorners + 1] = cone[6];
556: cells[cOff * numCorners + 2] = cone[7];
557: ++cOff;
558: cells[cOff * numCorners + 0] = cone[8];
559: cells[cOff * numCorners + 1] = cone[1];
560: cells[cOff * numCorners + 2] = cone[3];
561: ++cOff;
562: cells[cOff * numCorners + 0] = cone[8];
563: cells[cOff * numCorners + 1] = cone[3];
564: cells[cOff * numCorners + 2] = cone[5];
565: ++cOff;
566: cells[cOff * numCorners + 0] = cone[8];
567: cells[cOff * numCorners + 1] = cone[5];
568: cells[cOff * numCorners + 2] = cone[7];
569: ++cOff;
570: cells[cOff * numCorners + 0] = cone[8];
571: cells[cOff * numCorners + 1] = cone[7];
572: cells[cOff * numCorners + 2] = cone[1];
573: ++cOff;
574: break;
575: default:
576: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
577: }
578: if (debug) {
579: for (t = 0; t < Nt; ++t) {
580: PetscCall(PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %" PetscInt_FMT " (", t));
581: for (c = 0; c < numCorners; ++c) {
582: if (c > 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
583: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%" PetscInt_FMT, cells[(cOff - Nt + t) * numCorners + c]));
584: }
585: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
586: }
587: }
588: }
589: EG_free(lobjs);
590: }
591: }
592: PetscCheck(cOff == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %" PetscInt_FMT " != %" PetscInt_FMT " previous count", cOff, numCells);
593: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
594: PetscCall(PetscFree3(coords, cells, cone));
595: PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numCells, newCells));
596: PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " (%" PetscInt_FMT ")\n", numVertices, newVertices));
597: /* Embed EGADS model in DM */
598: {
599: PetscContainer modelObj, contextObj;
601: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
602: PetscCall(PetscContainerSetPointer(modelObj, model));
603: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
604: PetscCall(PetscContainerDestroy(&modelObj));
606: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
607: PetscCall(PetscContainerSetPointer(contextObj, context));
608: PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
609: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
610: PetscCall(PetscContainerDestroy(&contextObj));
611: }
612: /* Label points */
613: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
614: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
615: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
616: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
617: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
618: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
619: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
620: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
621: cOff = 0;
622: for (b = 0; b < nbodies; ++b) {
623: ego body = bodies[b];
624: int id, Nl, l;
626: PetscCall(EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs));
627: for (l = 0; l < Nl; ++l) {
628: ego loop = lobjs[l];
629: ego *fobjs;
630: int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;
632: lid = EG_indexBodyTopo(body, loop);
633: PetscCall(EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs));
634: PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
635: fid = EG_indexBodyTopo(body, fobjs[0]);
636: EG_free(fobjs);
637: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses));
638: for (e = 0; e < Ne; ++e) {
639: ego edge = objs[e];
640: int eid, Nv, v;
641: PetscInt points[3], support[2], numEdges, edgeNum;
642: const PetscInt *edges;
644: eid = EG_indexBodyTopo(body, edge);
645: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
646: if (mtype == DEGENERATE) continue;
647: else ++Ner;
648: for (v = 0; v < Nv; ++v) {
649: ego vertex = nobjs[v];
651: id = EG_indexBodyTopo(body, vertex);
652: PetscCall(DMLabelSetValue(edgeLabel, numCells + id - 1, eid));
653: points[v * 2] = numCells + id - 1;
654: }
655: PetscCall(PetscHMapIGet(edgeMap, eid - 1, &edgeNum));
656: points[1] = numCells + numVertices - newVertices + edgeNum;
658: PetscCall(DMLabelSetValue(edgeLabel, points[1], eid));
659: support[0] = points[0];
660: support[1] = points[1];
661: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
662: 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);
663: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
664: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
665: support[0] = points[1];
666: support[1] = points[2];
667: PetscCall(DMPlexGetJoin(dm, 2, support, &numEdges, &edges));
668: 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);
669: PetscCall(DMLabelSetValue(edgeLabel, edges[0], eid));
670: PetscCall(DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges));
671: }
672: switch (Ner) {
673: case 2:
674: Nt = 2;
675: break;
676: case 3:
677: Nt = 4;
678: break;
679: case 4:
680: Nt = 8;
681: break;
682: default:
683: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
684: }
685: for (t = 0; t < Nt; ++t) {
686: PetscCall(DMLabelSetValue(bodyLabel, cOff + t, b));
687: PetscCall(DMLabelSetValue(faceLabel, cOff + t, fid));
688: }
689: cOff += Nt;
690: }
691: EG_free(lobjs);
692: }
693: PetscCall(PetscHMapIDestroy(&edgeMap));
694: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
695: for (c = cStart; c < cEnd; ++c) {
696: PetscInt *closure = NULL;
697: PetscInt clSize, cl, bval, fval;
699: PetscCall(DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
700: PetscCall(DMLabelGetValue(bodyLabel, c, &bval));
701: PetscCall(DMLabelGetValue(faceLabel, c, &fval));
702: for (cl = 0; cl < clSize * 2; cl += 2) {
703: PetscCall(DMLabelSetValue(bodyLabel, closure[cl], bval));
704: PetscCall(DMLabelSetValue(faceLabel, closure[cl], fval));
705: }
706: PetscCall(DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure));
707: }
708: *newdm = dm;
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
713: {
714: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
715: // EGADS/EGADSLite variables
716: ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs;
717: ego topRef, prev, next;
718: int oclass, mtype, nbodies, *senses, *lSenses, *eSenses;
719: int b;
720: // PETSc variables
721: DM dm;
722: PetscHMapI edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL;
723: PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0;
724: PetscInt cellCntr = 0, numPoints = 0;
725: PetscInt *cells = NULL;
726: const PetscInt *cone = NULL;
727: PetscReal *coords = NULL;
728: PetscMPIInt rank;
730: PetscFunctionBeginUser;
731: PetscCallMPI(MPI_Comm_rank(comm, &rank));
732: if (rank == 0) {
733: // ---------------------------------------------------------------------------------------------------
734: // Generate Petsc Plex
735: // Get all Nodes in model, record coordinates in a correctly formatted array
736: // Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
737: // We need to uniformly refine the initial geometry to guarantee a valid mesh
739: // Calculate cell and vertex sizes
740: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
742: PetscCall(PetscHMapICreate(&edgeMap));
743: PetscCall(PetscHMapICreate(&bodyIndexMap));
744: PetscCall(PetscHMapICreate(&bodyVertexMap));
745: PetscCall(PetscHMapICreate(&bodyEdgeMap));
746: PetscCall(PetscHMapICreate(&bodyEdgeGlobalMap));
747: PetscCall(PetscHMapICreate(&bodyFaceMap));
749: for (b = 0; b < nbodies; ++b) {
750: ego body = bodies[b];
751: int Nf, Ne, Nv;
752: PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter;
753: PetscBool BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound;
755: PetscCall(PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound));
756: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
757: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
758: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
759: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
761: if (!BIfound) PetscCall(PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices));
762: if (!BVfound) PetscCall(PetscHMapISet(bodyVertexMap, b, numVertices));
763: if (!BEfound) PetscCall(PetscHMapISet(bodyEdgeMap, b, numEdges));
764: if (!BEGfound) PetscCall(PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr));
765: if (!BFfound) PetscCall(PetscHMapISet(bodyFaceMap, b, numFaces));
767: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
768: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
769: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
770: EG_free(fobjs);
771: EG_free(eobjs);
772: EG_free(nobjs);
774: // Remove DEGENERATE EDGES from Edge count
775: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
776: int Netemp = 0;
777: for (int e = 0; e < Ne; ++e) {
778: ego edge = eobjs[e];
779: int eid;
781: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
782: eid = EG_indexBodyTopo(body, edge);
784: PetscCall(PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound));
785: if (mtype == DEGENERATE) {
786: if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1));
787: } else {
788: ++Netemp;
789: if (!EMfound) PetscCall(PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp));
790: }
791: }
792: EG_free(eobjs);
794: // Determine Number of Cells
795: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
796: for (int f = 0; f < Nf; ++f) {
797: ego face = fobjs[f];
798: int edgeTemp = 0;
800: PetscCall(EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs));
801: for (int e = 0; e < Ne; ++e) {
802: ego edge = eobjs[e];
804: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
805: if (mtype != DEGENERATE) ++edgeTemp;
806: }
807: numCells += (2 * edgeTemp);
808: EG_free(eobjs);
809: }
810: EG_free(fobjs);
812: numFaces += Nf;
813: numEdges += Netemp;
814: numVertices += Nv;
815: edgeCntr += Ne;
816: }
818: // Set up basic DMPlex parameters
819: dim = 2; // Assumes 3D Models :: Need to handle 2D models in the future
820: cdim = 3; // Assumes 3D Models :: Need to update to handle 2D models in future
821: numCorners = 3; // Split Faces into triangles
822: numPoints = numVertices + numEdges + numFaces; // total number of coordinate points
824: PetscCall(PetscMalloc2(numPoints * cdim, &coords, numCells * numCorners, &cells));
826: // Get Vertex Coordinates and Set up Cells
827: for (b = 0; b < nbodies; ++b) {
828: ego body = bodies[b];
829: int Nf, Ne, Nv;
830: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
831: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
832: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
834: // Vertices on Current Body
835: PetscCall(EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs));
837: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
838: PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
839: PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));
841: for (int v = 0; v < Nv; ++v) {
842: ego vertex = nobjs[v];
843: double limits[4];
844: int id, dummy;
846: PetscCall(EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses));
847: id = EG_indexBodyTopo(body, vertex);
849: coords[(bodyVertexIndexStart + id - 1) * cdim + 0] = limits[0];
850: coords[(bodyVertexIndexStart + id - 1) * cdim + 1] = limits[1];
851: coords[(bodyVertexIndexStart + id - 1) * cdim + 2] = limits[2];
852: }
853: EG_free(nobjs);
855: // Edge Midpoint Vertices on Current Body
856: PetscCall(EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs));
858: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
859: PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
860: PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));
862: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
863: PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
864: PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));
866: for (int e = 0; e < Ne; ++e) {
867: ego edge = eobjs[e];
868: double range[2], avgt[1], cntrPnt[9];
869: int eid, eOffset;
870: int periodic;
872: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
873: if (mtype == DEGENERATE) continue;
875: eid = EG_indexBodyTopo(body, edge);
877: // get relative offset from globalEdgeID Vector
878: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
879: PetscCheck(EMfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1);
880: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
882: PetscCall(EG_getRange(edge, range, &periodic));
883: avgt[0] = (range[0] + range[1]) / 2.;
885: PetscCall(EG_evaluate(edge, avgt, cntrPnt));
886: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 0] = cntrPnt[0];
887: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 1] = cntrPnt[1];
888: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1) * cdim + 2] = cntrPnt[2];
889: }
890: EG_free(eobjs);
892: // Face Midpoint Vertices on Current Body
893: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
895: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
896: PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
897: PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));
899: for (int f = 0; f < Nf; ++f) {
900: ego face = fobjs[f];
901: double range[4], avgUV[2], cntrPnt[18];
902: int peri, id;
904: id = EG_indexBodyTopo(body, face);
905: PetscCall(EG_getRange(face, range, &peri));
907: avgUV[0] = (range[0] + range[1]) / 2.;
908: avgUV[1] = (range[2] + range[3]) / 2.;
909: PetscCall(EG_evaluate(face, avgUV, cntrPnt));
911: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 0] = cntrPnt[0];
912: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 1] = cntrPnt[1];
913: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1) * cdim + 2] = cntrPnt[2];
914: }
915: EG_free(fobjs);
917: // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity
918: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
919: for (int f = 0; f < Nf; ++f) {
920: ego face = fobjs[f];
921: int fID, midFaceID, midPntID, startID, endID, Nl;
923: fID = EG_indexBodyTopo(body, face);
924: midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1;
925: // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges.
926: // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are:
927: // 1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option.
928: // This will use a default EGADS tessellation as an initial surface mesh.
929: // 2) Create the initial surface mesh via a 2D mesher :: Currently not available (?future?)
930: // May I suggest the XXXX as a starting point?
932: PetscCall(EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses));
934: 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);
935: for (int l = 0; l < Nl; ++l) {
936: ego loop = lobjs[l];
938: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
939: for (int e = 0; e < Ne; ++e) {
940: ego edge = eobjs[e];
941: int eid, eOffset;
943: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
944: eid = EG_indexBodyTopo(body, edge);
945: if (mtype == DEGENERATE) continue;
947: // get relative offset from globalEdgeID Vector
948: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
949: 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);
950: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
952: midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1;
954: PetscCall(EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses));
956: if (eSenses[e] > 0) {
957: startID = EG_indexBodyTopo(body, nobjs[0]);
958: endID = EG_indexBodyTopo(body, nobjs[1]);
959: } else {
960: startID = EG_indexBodyTopo(body, nobjs[1]);
961: endID = EG_indexBodyTopo(body, nobjs[0]);
962: }
964: // Define 2 Cells per Edge with correct orientation
965: cells[cellCntr * numCorners + 0] = midFaceID;
966: cells[cellCntr * numCorners + 1] = bodyVertexIndexStart + startID - 1;
967: cells[cellCntr * numCorners + 2] = midPntID;
969: cells[cellCntr * numCorners + 3] = midFaceID;
970: cells[cellCntr * numCorners + 4] = midPntID;
971: cells[cellCntr * numCorners + 5] = bodyVertexIndexStart + endID - 1;
973: cellCntr = cellCntr + 2;
974: }
975: }
976: }
977: EG_free(fobjs);
978: }
979: }
981: // Generate DMPlex
982: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
983: PetscCall(PetscFree2(coords, cells));
984: PetscCall(PetscInfo(dm, " Total Number of Unique Cells = %" PetscInt_FMT " \n", numCells));
985: PetscCall(PetscInfo(dm, " Total Number of Unique Vertices = %" PetscInt_FMT " \n", numVertices));
987: // Embed EGADS model in DM
988: {
989: PetscContainer modelObj, contextObj;
991: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
992: PetscCall(PetscContainerSetPointer(modelObj, model));
993: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
994: PetscCall(PetscContainerDestroy(&modelObj));
996: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
997: PetscCall(PetscContainerSetPointer(contextObj, context));
998: PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
999: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
1000: PetscCall(PetscContainerDestroy(&contextObj));
1001: }
1002: // Label points
1003: PetscInt nStart, nEnd;
1005: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
1006: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
1007: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
1008: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
1009: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
1010: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1011: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
1012: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1014: PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));
1016: cellCntr = 0;
1017: for (b = 0; b < nbodies; ++b) {
1018: ego body = bodies[b];
1019: int Nv, Ne, Nf;
1020: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
1021: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
1022: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
1024: PetscCall(PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound));
1025: PetscCheck(BVfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
1026: PetscCall(PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart));
1028: PetscCall(PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound));
1029: PetscCheck(BEfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
1030: PetscCall(PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart));
1032: PetscCall(PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound));
1033: PetscCheck(BFfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
1034: PetscCall(PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart));
1036: PetscCall(PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound));
1037: PetscCheck(BEGfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
1038: PetscCall(PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart));
1040: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1041: for (int f = 0; f < Nf; ++f) {
1042: ego face = fobjs[f];
1043: int fID, Nl;
1045: fID = EG_indexBodyTopo(body, face);
1047: PetscCall(EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs));
1048: for (int l = 0; l < Nl; ++l) {
1049: ego loop = lobjs[l];
1050: int lid;
1052: lid = EG_indexBodyTopo(body, loop);
1053: PetscCheck(Nl <= 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
1055: PetscCall(EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses));
1056: for (int e = 0; e < Ne; ++e) {
1057: ego edge = eobjs[e];
1058: int eid, eOffset;
1060: // Skip DEGENERATE Edges
1061: PetscCall(EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next));
1062: if (mtype == DEGENERATE) continue;
1063: eid = EG_indexBodyTopo(body, edge);
1065: // get relative offset from globalEdgeID Vector
1066: PetscCall(PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound));
1067: 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);
1068: PetscCall(PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset));
1070: PetscCall(EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs));
1071: for (int v = 0; v < Nv; ++v) {
1072: ego vertex = nobjs[v];
1073: int vID;
1075: vID = EG_indexBodyTopo(body, vertex);
1076: PetscCall(DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b));
1077: PetscCall(DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID));
1078: }
1079: EG_free(nobjs);
1081: PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b));
1082: PetscCall(DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid));
1084: // Define Cell faces
1085: for (int jj = 0; jj < 2; ++jj) {
1086: PetscCall(DMLabelSetValue(bodyLabel, cellCntr, b));
1087: PetscCall(DMLabelSetValue(faceLabel, cellCntr, fID));
1088: PetscCall(DMPlexGetCone(dm, cellCntr, &cone));
1090: PetscCall(DMLabelSetValue(bodyLabel, cone[0], b));
1091: PetscCall(DMLabelSetValue(faceLabel, cone[0], fID));
1093: PetscCall(DMLabelSetValue(bodyLabel, cone[1], b));
1094: PetscCall(DMLabelSetValue(edgeLabel, cone[1], eid));
1096: PetscCall(DMLabelSetValue(bodyLabel, cone[2], b));
1097: PetscCall(DMLabelSetValue(faceLabel, cone[2], fID));
1099: cellCntr = cellCntr + 1;
1100: }
1101: }
1102: }
1103: EG_free(lobjs);
1105: PetscCall(DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b));
1106: PetscCall(DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID));
1107: }
1108: EG_free(fobjs);
1109: }
1111: PetscCall(PetscHMapIDestroy(&edgeMap));
1112: PetscCall(PetscHMapIDestroy(&bodyIndexMap));
1113: PetscCall(PetscHMapIDestroy(&bodyVertexMap));
1114: PetscCall(PetscHMapIDestroy(&bodyEdgeMap));
1115: PetscCall(PetscHMapIDestroy(&bodyEdgeGlobalMap));
1116: PetscCall(PetscHMapIDestroy(&bodyFaceMap));
1118: *newdm = dm;
1119: PetscFunctionReturn(PETSC_SUCCESS);
1120: }
1122: static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
1123: {
1124: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1125: /* EGADSLite variables */
1126: ego geom, *bodies, *fobjs;
1127: int b, oclass, mtype, nbodies, *senses;
1128: int totalNumTris = 0, totalNumPoints = 0;
1129: double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize;
1130: /* PETSc variables */
1131: DM dm;
1132: PetscHMapI pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL;
1133: PetscHMapI pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL;
1134: PetscInt dim = -1, cdim = -1, numCorners = 0, counter = 0;
1135: PetscInt *cells = NULL;
1136: const PetscInt *cone = NULL;
1137: PetscReal *coords = NULL;
1138: PetscMPIInt rank;
1140: PetscFunctionBeginUser;
1141: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1142: if (rank == 0) {
1143: // ---------------------------------------------------------------------------------------------------
1144: // Generate Petsc Plex from EGADSlite created Tessellation of geometry
1145: // ---------------------------------------------------------------------------------------------------
1147: // Calculate cell and vertex sizes
1148: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses));
1150: PetscCall(PetscHMapICreate(&pointIndexStartMap));
1151: PetscCall(PetscHMapICreate(&triIndexStartMap));
1152: PetscCall(PetscHMapICreate(&pTypeLabelMap));
1153: PetscCall(PetscHMapICreate(&pIndexLabelMap));
1154: PetscCall(PetscHMapICreate(&pBodyIndexLabelMap));
1155: PetscCall(PetscHMapICreate(&triFaceIDLabelMap));
1156: PetscCall(PetscHMapICreate(&triBodyIDLabelMap));
1158: /* Create Tessellation of Bodies */
1159: ego tessArray[nbodies];
1161: for (b = 0; b < nbodies; ++b) {
1162: ego body = bodies[b];
1163: double params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation
1164: int Nf, bodyNumPoints = 0, bodyNumTris = 0;
1165: PetscHashIter PISiter, TISiter;
1166: PetscBool PISfound, TISfound;
1168: /* Store Start Index for each Body's Point and Tris */
1169: PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
1170: PetscCall(PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound));
1172: if (!PISfound) PetscCall(PetscHMapISet(pointIndexStartMap, b, totalNumPoints));
1173: if (!TISfound) PetscCall(PetscHMapISet(triIndexStartMap, b, totalNumTris));
1175: /* Calculate Tessellation parameters based on Bounding Box */
1176: /* Get Bounding Box Dimensions of the BODY */
1177: PetscCall(EG_getBoundingBox(body, boundBox));
1178: tessSize = boundBox[3] - boundBox[0];
1179: if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1];
1180: if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2];
1182: // TODO :: May want to give users tessellation parameter options //
1183: params[0] = 0.0250 * tessSize;
1184: params[1] = 0.0075 * tessSize;
1185: params[2] = 15.0;
1187: PetscCall(EG_makeTessBody(body, params, &tessArray[b]));
1189: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1191: for (int f = 0; f < Nf; ++f) {
1192: ego face = fobjs[f];
1193: int len, fID, ntris;
1194: const int *ptype, *pindex, *ptris, *ptric;
1195: const double *pxyz, *puv;
1197: // Get Face ID //
1198: fID = EG_indexBodyTopo(body, face);
1200: // Checkout the Surface Tessellation //
1201: PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));
1203: // Determine total number of triangle cells in the tessellation //
1204: bodyNumTris += (int)ntris;
1206: // Check out the point index and coordinate //
1207: for (int p = 0; p < len; ++p) {
1208: int global;
1210: PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));
1212: // Determine the total number of points in the tessellation //
1213: bodyNumPoints = PetscMax(bodyNumPoints, global);
1214: }
1215: }
1216: EG_free(fobjs);
1218: totalNumPoints += bodyNumPoints;
1219: totalNumTris += bodyNumTris;
1220: }
1221: //} - Original End of (rank == 0)
1223: dim = 2;
1224: cdim = 3;
1225: numCorners = 3;
1226: //PetscInt counter = 0;
1228: /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA */
1229: /* Fill in below and use to define DMLabels after DMPlex creation */
1230: PetscCall(PetscMalloc2(totalNumPoints * cdim, &coords, totalNumTris * numCorners, &cells));
1232: for (b = 0; b < nbodies; ++b) {
1233: ego body = bodies[b];
1234: int Nf;
1235: PetscInt pointIndexStart;
1236: PetscHashIter PISiter;
1237: PetscBool PISfound;
1239: PetscCall(PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound));
1240: PetscCheck(PISfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b);
1241: PetscCall(PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart));
1243: PetscCall(EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs));
1245: for (int f = 0; f < Nf; ++f) {
1246: /* Get Face Object */
1247: ego face = fobjs[f];
1248: int len, fID, ntris;
1249: const int *ptype, *pindex, *ptris, *ptric;
1250: const double *pxyz, *puv;
1252: /* Get Face ID */
1253: fID = EG_indexBodyTopo(body, face);
1255: /* Checkout the Surface Tessellation */
1256: PetscCall(EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric));
1258: /* Check out the point index and coordinate */
1259: for (int p = 0; p < len; ++p) {
1260: int global;
1261: PetscHashIter PTLiter, PILiter, PBLiter;
1262: PetscBool PTLfound, PILfound, PBLfound;
1264: PetscCall(EG_localToGlobal(tessArray[b], fID, p + 1, &global));
1266: /* Set the coordinates array for DAG */
1267: coords[((global - 1 + pointIndexStart) * 3) + 0] = pxyz[(p * 3) + 0];
1268: coords[((global - 1 + pointIndexStart) * 3) + 1] = pxyz[(p * 3) + 1];
1269: coords[((global - 1 + pointIndexStart) * 3) + 2] = pxyz[(p * 3) + 2];
1271: /* Store Geometry Label Information for DMLabel assignment later */
1272: PetscCall(PetscHMapIFind(pTypeLabelMap, global - 1 + pointIndexStart, &PTLiter, &PTLfound));
1273: PetscCall(PetscHMapIFind(pIndexLabelMap, global - 1 + pointIndexStart, &PILiter, &PILfound));
1274: PetscCall(PetscHMapIFind(pBodyIndexLabelMap, global - 1 + pointIndexStart, &PBLiter, &PBLfound));
1276: if (!PTLfound) PetscCall(PetscHMapISet(pTypeLabelMap, global - 1 + pointIndexStart, ptype[p]));
1277: if (!PILfound) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, pindex[p]));
1278: if (!PBLfound) PetscCall(PetscHMapISet(pBodyIndexLabelMap, global - 1 + pointIndexStart, b));
1280: if (ptype[p] < 0) PetscCall(PetscHMapISet(pIndexLabelMap, global - 1 + pointIndexStart, fID));
1281: }
1283: for (int t = 0; t < (int)ntris; ++t) {
1284: int global, globalA, globalB;
1285: PetscHashIter TFLiter, TBLiter;
1286: PetscBool TFLfound, TBLfound;
1288: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 0], &global));
1289: cells[(counter * 3) + 0] = global - 1 + pointIndexStart;
1291: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 1], &globalA));
1292: cells[(counter * 3) + 1] = globalA - 1 + pointIndexStart;
1294: PetscCall(EG_localToGlobal(tessArray[b], fID, ptris[(t * 3) + 2], &globalB));
1295: cells[(counter * 3) + 2] = globalB - 1 + pointIndexStart;
1297: PetscCall(PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound));
1298: PetscCall(PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound));
1300: if (!TFLfound) PetscCall(PetscHMapISet(triFaceIDLabelMap, counter, fID));
1301: if (!TBLfound) PetscCall(PetscHMapISet(triBodyIDLabelMap, counter, b));
1303: counter += 1;
1304: }
1305: }
1306: EG_free(fobjs);
1307: }
1308: }
1310: //Build DMPlex
1311: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm));
1312: PetscCall(PetscFree2(coords, cells));
1314: // Embed EGADS model in DM
1315: {
1316: PetscContainer modelObj, contextObj;
1318: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &modelObj));
1319: PetscCall(PetscContainerSetPointer(modelObj, model));
1320: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Model", (PetscObject)modelObj));
1321: PetscCall(PetscContainerDestroy(&modelObj));
1323: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &contextObj));
1324: PetscCall(PetscContainerSetPointer(contextObj, context));
1325: PetscCall(PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private));
1326: PetscCall(PetscObjectCompose((PetscObject)dm, "EGADS Context", (PetscObject)contextObj));
1327: PetscCall(PetscContainerDestroy(&contextObj));
1328: }
1330: // Label Points
1331: PetscCall(DMCreateLabel(dm, "EGADS Body ID"));
1332: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
1333: PetscCall(DMCreateLabel(dm, "EGADS Face ID"));
1334: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
1335: PetscCall(DMCreateLabel(dm, "EGADS Edge ID"));
1336: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1337: PetscCall(DMCreateLabel(dm, "EGADS Vertex ID"));
1338: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1340: /* Get Number of DAG Nodes at each level */
1341: int fStart, fEnd, eStart, eEnd, nStart, nEnd;
1343: PetscCall(DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd));
1344: PetscCall(DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd));
1345: PetscCall(DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd));
1347: /* Set DMLabels for NODES */
1348: for (int n = nStart; n < nEnd; ++n) {
1349: int pTypeVal, pIndexVal, pBodyVal;
1350: PetscHashIter PTLiter, PILiter, PBLiter;
1351: PetscBool PTLfound, PILfound, PBLfound;
1353: //Converted to Hash Tables
1354: PetscCall(PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound));
1355: PetscCheck(PTLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n);
1356: PetscCall(PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal));
1358: PetscCall(PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound));
1359: PetscCheck(PILfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n);
1360: PetscCall(PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal));
1362: PetscCall(PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound));
1363: PetscCheck(PBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n);
1364: PetscCall(PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal));
1366: PetscCall(DMLabelSetValue(bodyLabel, n, pBodyVal));
1367: if (pTypeVal == 0) PetscCall(DMLabelSetValue(vertexLabel, n, pIndexVal));
1368: if (pTypeVal > 0) PetscCall(DMLabelSetValue(edgeLabel, n, pIndexVal));
1369: if (pTypeVal < 0) PetscCall(DMLabelSetValue(faceLabel, n, pIndexVal));
1370: }
1372: /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */
1373: for (int e = eStart; e < eEnd; ++e) {
1374: int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1;
1376: PetscCall(DMPlexGetCone(dm, e, &cone));
1377: PetscCall(DMLabelGetValue(bodyLabel, cone[0], &bodyID_0)); // Do I need to check the other end?
1378: PetscCall(DMLabelGetValue(vertexLabel, cone[0], &vertexID_0));
1379: PetscCall(DMLabelGetValue(vertexLabel, cone[1], &vertexID_1));
1380: PetscCall(DMLabelGetValue(edgeLabel, cone[0], &edgeID_0));
1381: PetscCall(DMLabelGetValue(edgeLabel, cone[1], &edgeID_1));
1382: PetscCall(DMLabelGetValue(faceLabel, cone[0], &faceID_0));
1383: PetscCall(DMLabelGetValue(faceLabel, cone[1], &faceID_1));
1385: PetscCall(DMLabelSetValue(bodyLabel, e, bodyID_0));
1387: if (edgeID_0 == edgeID_1) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
1388: else if (vertexID_0 > 0 && edgeID_1 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_1));
1389: else if (vertexID_1 > 0 && edgeID_0 > 0) PetscCall(DMLabelSetValue(edgeLabel, e, edgeID_0));
1390: else { /* Do Nothing */ }
1391: }
1393: /* Set DMLabels for Cells */
1394: for (int f = fStart; f < fEnd; ++f) {
1395: int edgeID_0;
1396: PetscInt triBodyVal, triFaceVal;
1397: PetscHashIter TFLiter, TBLiter;
1398: PetscBool TFLfound, TBLfound;
1400: // Convert to Hash Table
1401: PetscCall(PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound));
1402: PetscCheck(TFLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f);
1403: PetscCall(PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal));
1405: PetscCall(PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound));
1406: PetscCheck(TBLfound, PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f);
1407: PetscCall(PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal));
1409: PetscCall(DMLabelSetValue(bodyLabel, f, triBodyVal));
1410: PetscCall(DMLabelSetValue(faceLabel, f, triFaceVal));
1412: /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */
1413: PetscCall(DMPlexGetCone(dm, f, &cone));
1415: for (int jj = 0; jj < 3; ++jj) {
1416: PetscCall(DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0));
1418: if (edgeID_0 < 0) {
1419: PetscCall(DMLabelSetValue(bodyLabel, cone[jj], triBodyVal));
1420: PetscCall(DMLabelSetValue(faceLabel, cone[jj], triFaceVal));
1421: }
1422: }
1423: }
1425: *newdm = dm;
1426: PetscFunctionReturn(PETSC_SUCCESS);
1427: }
1428: #endif
1430: /*@
1431: DMPlexInflateToGeomModel - Snaps the vertex coordinates of a `DMPLEX` object representing the mesh to its geometry if some vertices depart from the model.
1432: This usually happens with non-conforming refinement.
1434: Collective
1436: Input Parameter:
1437: . dm - The uninflated `DM` object representing the mesh
1439: Level: intermediate
1441: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`
1442: @*/
1443: PetscErrorCode DMPlexInflateToGeomModel(DM dm)
1444: {
1445: #if defined(PETSC_HAVE_EGADS)
1446: /* EGADS Variables */
1447: ego model, geom, body, face, edge;
1448: ego *bodies;
1449: int Nb, oclass, mtype, *senses;
1450: double result[3];
1451: /* PETSc Variables */
1452: DM cdm;
1453: PetscContainer modelObj;
1454: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1455: Vec coordinates;
1456: PetscScalar *coords;
1457: PetscInt bodyID, faceID, edgeID, vertexID;
1458: PetscInt cdim, d, vStart, vEnd, v;
1459: #endif
1461: PetscFunctionBegin;
1462: #if defined(PETSC_HAVE_EGADS)
1463: PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
1464: if (!modelObj) PetscFunctionReturn(PETSC_SUCCESS);
1465: PetscCall(DMGetCoordinateDim(dm, &cdim));
1466: PetscCall(DMGetCoordinateDM(dm, &cdm));
1467: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1468: PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
1469: PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
1470: PetscCall(DMGetLabel(dm, "EGADS Edge ID", &edgeLabel));
1471: PetscCall(DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel));
1473: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
1474: PetscCall(EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses));
1476: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1477: PetscCall(VecGetArrayWrite(coordinates, &coords));
1478: for (v = vStart; v < vEnd; ++v) {
1479: PetscScalar *vcoords;
1481: PetscCall(DMLabelGetValue(bodyLabel, v, &bodyID));
1482: PetscCall(DMLabelGetValue(faceLabel, v, &faceID));
1483: PetscCall(DMLabelGetValue(edgeLabel, v, &edgeID));
1484: PetscCall(DMLabelGetValue(vertexLabel, v, &vertexID));
1486: PetscCheck(bodyID < Nb, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %" PetscInt_FMT " is not in [0, %d)", bodyID, Nb);
1487: body = bodies[bodyID];
1489: PetscCall(DMPlexPointLocalRef(cdm, v, coords, (void *)&vcoords));
1490: if (edgeID > 0) {
1491: /* Snap to EDGE at nearest location */
1492: double params[1];
1493: PetscCall(EG_objectBodyTopo(body, EDGE, edgeID, &edge));
1494: PetscCall(EG_invEvaluate(edge, vcoords, params, result)); // Get (x,y,z) of nearest point on EDGE
1495: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1496: } else if (faceID > 0) {
1497: /* Snap to FACE at nearest location */
1498: double params[2];
1499: PetscCall(EG_objectBodyTopo(body, FACE, faceID, &face));
1500: PetscCall(EG_invEvaluate(face, vcoords, params, result)); // Get (x,y,z) of nearest point on FACE
1501: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1502: }
1503: }
1504: PetscCall(VecRestoreArrayWrite(coordinates, &coords));
1505: /* Clear out global coordinates */
1506: PetscCall(VecDestroy(&dm->coordinates[0].x));
1507: #endif
1508: PetscFunctionReturn(PETSC_SUCCESS);
1509: }
1511: /*@C
1512: DMPlexCreateEGADSFromFile - Create a `DMPLEX` mesh from an EGADS, IGES, or STEP file.
1514: Collective
1516: Input Parameters:
1517: + comm - The MPI communicator
1518: - filename - The name of the EGADS, IGES, or STEP file
1520: Output Parameter:
1521: . dm - The `DM` object representing the mesh
1523: Level: beginner
1525: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`, `DMPlexCreateEGADS()`, `DMPlexCreateEGADSLiteFromFile()`
1526: @*/
1527: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
1528: {
1529: PetscMPIInt rank;
1530: #if defined(PETSC_HAVE_EGADS)
1531: ego context = NULL, model = NULL;
1532: #endif
1533: PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE;
1535: PetscFunctionBegin;
1536: PetscAssertPointer(filename, 2);
1537: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL));
1538: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL));
1539: PetscCall(PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL));
1540: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1541: #if defined(PETSC_HAVE_EGADS)
1542: if (rank == 0) {
1543: PetscCall(EG_open(&context));
1544: PetscCall(EG_loadModel(context, 0, filename, &model));
1545: if (printModel) PetscCall(DMPlexEGADSPrintModel_Internal(model));
1546: }
1547: if (tessModel) PetscCall(DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm));
1548: else if (newModel) PetscCall(DMPlexCreateEGADS(comm, context, model, dm));
1549: else PetscCall(DMPlexCreateEGADS_Internal(comm, context, model, dm));
1550: PetscFunctionReturn(PETSC_SUCCESS);
1551: #else
1552: SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads");
1553: #endif
1554: }