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;
40: DMGetCoordinateDM(dm, &cdm);
41: DMGetCoordinateDim(dm, &dE);
42: DMGetCoordinatesLocal(dm, &coordinatesLocal);
43: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
44: if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
45: body = bodies[bodyID];
47: if (edgeID >= 0) {EG_objectBodyTopo(body, EDGE, edgeID, &obj); Np = 1;}
48: else if (faceID >= 0) {EG_objectBodyTopo(body, FACE, faceID, &obj); Np = 2;}
49: else {
50: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
51: return(0);
52: }
54: /* Calculate parameters (t or u,v) for vertices */
55: DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
56: Nv /= dE;
57: if (Nv == 1) {
58: DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
59: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
60: return(0);
61: }
62: if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);
64: /* Correct EGADSlite 2pi bug when calculating nearest point on Periodic Surfaces */
65: EG_getRange(obj, range, &peri);
66: for (v = 0; v < Nv; ++v) {
67: EG_invEvaluate(obj, &coords[v*dE], ¶msV[v*3], &resultV[v*3]);
68: #if 1
69: if (peri > 0) {
70: if (paramsV[v*3+0] + 1.e-4 < range[0]) {paramsV[v*3+0] += 2. * PETSC_PI;}
71: else if (paramsV[v*3+0] - 1.e-4 > range[1]) {paramsV[v*3+0] -= 2. * PETSC_PI;}
72: }
73: if (peri > 1) {
74: if (paramsV[v*3+1] + 1.e-4 < range[2]) {paramsV[v*3+1] += 2. * PETSC_PI;}
75: else if (paramsV[v*3+1] - 1.e-4 > range[3]) {paramsV[v*3+1] -= 2. * PETSC_PI;}
76: }
77: #endif
78: }
79: DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
80: /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
81: for (pm = 0; pm < Np; ++pm) {
82: params[pm] = 0.;
83: for (v = 0; v < Nv; ++v) params[pm] += paramsV[v*3+pm];
84: params[pm] /= Nv;
85: }
86: if ((params[0] < range[0]) || (params[0] > range[1])) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D had bad interpolation", p);
87: if (Np > 1 && ((params[1] < range[2]) || (params[1] > range[3]))) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %D had bad interpolation", p);
88: /* Put coordinates for new vertex in result[] */
89: EG_evaluate(obj, params, result);
90: for (d = 0; d < dE; ++d) gcoords[d] = result[d];
91: return(0);
92: }
93: #endif
95: /*@
96: 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.
98: Not collective
100: Input Parameters:
101: + dm - The DMPlex object
102: . p - The mesh point
103: - mcoords - A coordinate point lying on the mesh point
105: Output Parameter:
106: . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point
108: Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS.
110: Level: intermediate
112: .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform()
113: @*/
114: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[])
115: {
116: PetscInt dE, d;
120: DMGetCoordinateDim(dm, &dE);
121: #ifdef PETSC_HAVE_EGADS
122: {
123: DM_Plex *plex = (DM_Plex *) dm->data;
124: DMLabel bodyLabel, faceLabel, edgeLabel;
125: PetscInt bodyID, faceID, edgeID;
126: PetscContainer modelObj;
127: ego model;
128: PetscBool islite = PETSC_FALSE;
130: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
131: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
132: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
133: if (!bodyLabel || !faceLabel || !edgeLabel || plex->ignoreModel) {
134: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
135: return(0);
136: }
137: PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
138: if (!modelObj) {
139: PetscObjectQuery((PetscObject) dm, "EGADSLite Model", (PetscObject *) &modelObj);
140: islite = PETSC_TRUE;
141: }
142: if (!modelObj) {
143: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
144: return(0);
145: }
146: PetscContainerGetPointer(modelObj, (void **) &model);
147: DMLabelGetValue(bodyLabel, p, &bodyID);
148: DMLabelGetValue(faceLabel, p, &faceID);
149: DMLabelGetValue(edgeLabel, p, &edgeID);
150: /* Allows for "Connective" Plex Edges present in models with multiple non-touching Entities */
151: if (bodyID < 0) {
152: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
153: return(0);
154: }
155: if (islite) {DMPlexSnapToGeomModel_EGADSLite_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords);}
156: else {DMPlexSnapToGeomModel_EGADS_Internal(dm, p, model, bodyID, faceID, edgeID, mcoords, gcoords);}
157: }
158: #else
159: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
160: #endif
161: return(0);
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;
173: /* test bodyTopo functions */
174: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
175: PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb);
177: for (b = 0; b < Nb; ++b) {
178: ego body = bodies[b];
179: int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;
181: /* Output Basic Model Topology */
182: EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
183: PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh);
184: EG_free(objs);
186: EG_getBodyTopos(body, NULL, FACE, &Nf, &objs);
187: PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf);
188: EG_free(objs);
190: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
191: PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl);
193: EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs);
194: PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne);
195: EG_free(objs);
197: EG_getBodyTopos(body, NULL, NODE, &Nv, &objs);
198: PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv);
199: EG_free(objs);
201: for (l = 0; l < Nl; ++l) {
202: ego loop = lobjs[l];
204: id = EG_indexBodyTopo(body, loop);
205: PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id);
207: /* Get EDGE info which associated with the current LOOP */
208: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
210: for (e = 0; e < Ne; ++e) {
211: ego edge = objs[e];
212: double range[4] = {0., 0., 0., 0.};
213: double point[3] = {0., 0., 0.};
214: double params[3] = {0., 0., 0.};
215: double result[18];
216: int peri;
218: id = EG_indexBodyTopo(body, edge);
219: PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e);
221: EG_getRange(edge, range, &peri);
222: PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]);
224: /* Get NODE info which associated with the current EDGE */
225: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
226: if (mtype == DEGENERATE) {
227: PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id);
228: } else {
229: params[0] = range[0];
230: EG_evaluate(edge, params, result);
231: PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]);
232: params[0] = range[1];
233: EG_evaluate(edge, params, result);
234: PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]);
235: }
237: for (v = 0; v < Nv; ++v) {
238: ego vertex = nobjs[v];
239: double limits[4];
240: int dummy;
242: EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
243: id = EG_indexBodyTopo(body, vertex);
244: PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);
245: PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]);
247: point[0] = point[0] + limits[0];
248: point[1] = point[1] + limits[1];
249: point[2] = point[2] + limits[2];
250: }
251: }
252: }
253: EG_free(lobjs);
254: }
255: return(0);
256: }
258: static PetscErrorCode DMPlexEGADSDestroy_Private(void *context)
259: {
260: if (context) EG_close((ego) context);
261: return 0;
262: }
264: static PetscErrorCode DMPlexCreateEGADS_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
265: {
266: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
267: PetscInt cStart, cEnd, c;
268: /* EGADSLite variables */
269: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
270: int oclass, mtype, nbodies, *senses;
271: int b;
272: /* PETSc variables */
273: DM dm;
274: PetscHMapI edgeMap = NULL;
275: 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;
276: PetscInt *cells = NULL, *cone = NULL;
277: PetscReal *coords = NULL;
278: PetscMPIInt rank;
282: MPI_Comm_rank(comm, &rank);
283: if (rank == 0) {
284: const PetscInt debug = 0;
286: /* ---------------------------------------------------------------------------------------------------
287: Generate Petsc Plex
288: Get all Nodes in model, record coordinates in a correctly formatted array
289: Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
290: We need to uniformly refine the initial geometry to guarantee a valid mesh
291: */
293: /* Calculate cell and vertex sizes */
294: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
295: PetscHMapICreate(&edgeMap);
296: numEdges = 0;
297: for (b = 0; b < nbodies; ++b) {
298: ego body = bodies[b];
299: int id, Nl, l, Nv, v;
301: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
302: for (l = 0; l < Nl; ++l) {
303: ego loop = lobjs[l];
304: int Ner = 0, Ne, e, Nc;
306: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
307: for (e = 0; e < Ne; ++e) {
308: ego edge = objs[e];
309: int Nv, id;
310: PetscHashIter iter;
311: PetscBool found;
313: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
314: if (mtype == DEGENERATE) continue;
315: id = EG_indexBodyTopo(body, edge);
316: PetscHMapIFind(edgeMap, id-1, &iter, &found);
317: if (!found) {PetscHMapISet(edgeMap, id-1, numEdges++);}
318: ++Ner;
319: }
320: if (Ner == 2) {Nc = 2;}
321: else if (Ner == 3) {Nc = 4;}
322: else if (Ner == 4) {Nc = 8; ++numQuads;}
323: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
324: numCells += Nc;
325: newCells += Nc-1;
326: maxCorners = PetscMax(Ner*2+1, maxCorners);
327: }
328: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
329: for (v = 0; v < Nv; ++v) {
330: ego vertex = nobjs[v];
332: id = EG_indexBodyTopo(body, vertex);
333: /* TODO: Instead of assuming contiguous ids, we could use a hash table */
334: numVertices = PetscMax(id, numVertices);
335: }
336: EG_free(lobjs);
337: EG_free(nobjs);
338: }
339: PetscHMapIGetSize(edgeMap, &numEdges);
340: newVertices = numEdges + numQuads;
341: numVertices += newVertices;
343: dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
344: cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
345: numCorners = 3; /* Split cells into triangles */
346: PetscMalloc3(numVertices*cdim, &coords, numCells*numCorners, &cells, maxCorners, &cone);
348: /* Get vertex coordinates */
349: for (b = 0; b < nbodies; ++b) {
350: ego body = bodies[b];
351: int id, Nv, v;
353: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
354: for (v = 0; v < Nv; ++v) {
355: ego vertex = nobjs[v];
356: double limits[4];
357: int dummy;
359: EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
360: id = EG_indexBodyTopo(body, vertex);
361: coords[(id-1)*cdim+0] = limits[0];
362: coords[(id-1)*cdim+1] = limits[1];
363: coords[(id-1)*cdim+2] = limits[2];
364: }
365: EG_free(nobjs);
366: }
367: PetscHMapIClear(edgeMap);
368: fOff = numVertices - newVertices + numEdges;
369: numEdges = 0;
370: numQuads = 0;
371: for (b = 0; b < nbodies; ++b) {
372: ego body = bodies[b];
373: int Nl, l;
375: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
376: for (l = 0; l < Nl; ++l) {
377: ego loop = lobjs[l];
378: int lid, Ner = 0, Ne, e;
380: lid = EG_indexBodyTopo(body, loop);
381: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
382: for (e = 0; e < Ne; ++e) {
383: ego edge = objs[e];
384: int eid, Nv;
385: PetscHashIter iter;
386: PetscBool found;
388: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
389: if (mtype == DEGENERATE) continue;
390: ++Ner;
391: eid = EG_indexBodyTopo(body, edge);
392: PetscHMapIFind(edgeMap, eid-1, &iter, &found);
393: if (!found) {
394: PetscInt v = numVertices - newVertices + numEdges;
395: double range[4], params[3] = {0., 0., 0.}, result[18];
396: int periodic[2];
398: PetscHMapISet(edgeMap, eid-1, numEdges++);
399: EG_getRange(edge, range, periodic);
400: params[0] = 0.5*(range[0] + range[1]);
401: EG_evaluate(edge, params, result);
402: coords[v*cdim+0] = result[0];
403: coords[v*cdim+1] = result[1];
404: coords[v*cdim+2] = result[2];
405: }
406: }
407: if (Ner == 4) {
408: PetscInt v = fOff + numQuads++;
409: ego *fobjs, face;
410: double range[4], params[3] = {0., 0., 0.}, result[18];
411: int Nf, fid, periodic[2];
413: EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
414: face = fobjs[0];
415: fid = EG_indexBodyTopo(body, face);
416: if (Nf != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid);
417: EG_getRange(face, range, periodic);
418: params[0] = 0.5*(range[0] + range[1]);
419: params[1] = 0.5*(range[2] + range[3]);
420: EG_evaluate(face, params, result);
421: coords[v*cdim+0] = result[0];
422: coords[v*cdim+1] = result[1];
423: coords[v*cdim+2] = result[2];
424: }
425: }
426: }
427: if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices);
429: /* Get cell vertices by traversing loops */
430: numQuads = 0;
431: cOff = 0;
432: for (b = 0; b < nbodies; ++b) {
433: ego body = bodies[b];
434: int id, Nl, l;
436: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
437: for (l = 0; l < Nl; ++l) {
438: ego loop = lobjs[l];
439: int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;
441: lid = EG_indexBodyTopo(body, loop);
442: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
444: for (e = 0; e < Ne; ++e) {
445: ego edge = objs[e];
446: int points[3];
447: int eid, Nv, v, tmp;
449: eid = EG_indexBodyTopo(body, edge);
450: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
451: if (mtype == DEGENERATE) continue;
452: else ++Ner;
453: if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);
455: for (v = 0; v < Nv; ++v) {
456: ego vertex = nobjs[v];
458: id = EG_indexBodyTopo(body, vertex);
459: points[v*2] = id-1;
460: }
461: {
462: PetscInt edgeNum;
464: PetscHMapIGet(edgeMap, eid-1, &edgeNum);
465: points[1] = numVertices - newVertices + edgeNum;
466: }
467: /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
468: if (!nc) {
469: for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v];
470: } else {
471: if (cone[nc-1] == points[0]) {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
472: else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
473: else if (cone[nc-3] == points[0]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
474: else if (cone[nc-3] == points[2]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
475: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
476: }
477: }
478: if (nc != 2*Ner) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner);
479: if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;}
480: if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners);
481: /* Triangulate the loop */
482: switch (Ner) {
483: case 2: /* Bi-Segment -> 2 triangles */
484: Nt = 2;
485: cells[cOff*numCorners+0] = cone[0];
486: cells[cOff*numCorners+1] = cone[1];
487: cells[cOff*numCorners+2] = cone[2];
488: ++cOff;
489: cells[cOff*numCorners+0] = cone[0];
490: cells[cOff*numCorners+1] = cone[2];
491: cells[cOff*numCorners+2] = cone[3];
492: ++cOff;
493: break;
494: case 3: /* Triangle -> 4 triangles */
495: Nt = 4;
496: cells[cOff*numCorners+0] = cone[0];
497: cells[cOff*numCorners+1] = cone[1];
498: cells[cOff*numCorners+2] = cone[5];
499: ++cOff;
500: cells[cOff*numCorners+0] = cone[1];
501: cells[cOff*numCorners+1] = cone[2];
502: cells[cOff*numCorners+2] = cone[3];
503: ++cOff;
504: cells[cOff*numCorners+0] = cone[5];
505: cells[cOff*numCorners+1] = cone[3];
506: cells[cOff*numCorners+2] = cone[4];
507: ++cOff;
508: cells[cOff*numCorners+0] = cone[1];
509: cells[cOff*numCorners+1] = cone[3];
510: cells[cOff*numCorners+2] = cone[5];
511: ++cOff;
512: break;
513: case 4: /* Quad -> 8 triangles */
514: Nt = 8;
515: cells[cOff*numCorners+0] = cone[0];
516: cells[cOff*numCorners+1] = cone[1];
517: cells[cOff*numCorners+2] = cone[7];
518: ++cOff;
519: cells[cOff*numCorners+0] = cone[1];
520: cells[cOff*numCorners+1] = cone[2];
521: cells[cOff*numCorners+2] = cone[3];
522: ++cOff;
523: cells[cOff*numCorners+0] = cone[3];
524: cells[cOff*numCorners+1] = cone[4];
525: cells[cOff*numCorners+2] = cone[5];
526: ++cOff;
527: cells[cOff*numCorners+0] = cone[5];
528: cells[cOff*numCorners+1] = cone[6];
529: cells[cOff*numCorners+2] = cone[7];
530: ++cOff;
531: cells[cOff*numCorners+0] = cone[8];
532: cells[cOff*numCorners+1] = cone[1];
533: cells[cOff*numCorners+2] = cone[3];
534: ++cOff;
535: cells[cOff*numCorners+0] = cone[8];
536: cells[cOff*numCorners+1] = cone[3];
537: cells[cOff*numCorners+2] = cone[5];
538: ++cOff;
539: cells[cOff*numCorners+0] = cone[8];
540: cells[cOff*numCorners+1] = cone[5];
541: cells[cOff*numCorners+2] = cone[7];
542: ++cOff;
543: cells[cOff*numCorners+0] = cone[8];
544: cells[cOff*numCorners+1] = cone[7];
545: cells[cOff*numCorners+2] = cone[1];
546: ++cOff;
547: break;
548: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
549: }
550: if (debug) {
551: for (t = 0; t < Nt; ++t) {
552: PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t);
553: for (c = 0; c < numCorners; ++c) {
554: if (c > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
555: PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);
556: }
557: PetscPrintf(PETSC_COMM_SELF, ")\n");
558: }
559: }
560: }
561: EG_free(lobjs);
562: }
563: }
564: if (cOff != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells);
565: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
566: PetscFree3(coords, cells, cone);
567: PetscInfo2(dm, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);
568: PetscInfo2(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices);
569: /* Embed EGADS model in DM */
570: {
571: PetscContainer modelObj, contextObj;
573: PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
574: PetscContainerSetPointer(modelObj, model);
575: PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
576: PetscContainerDestroy(&modelObj);
578: PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
579: PetscContainerSetPointer(contextObj, context);
580: PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
581: PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
582: PetscContainerDestroy(&contextObj);
583: }
584: /* Label points */
585: DMCreateLabel(dm, "EGADS Body ID");
586: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
587: DMCreateLabel(dm, "EGADS Face ID");
588: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
589: DMCreateLabel(dm, "EGADS Edge ID");
590: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
591: DMCreateLabel(dm, "EGADS Vertex ID");
592: DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
593: cOff = 0;
594: for (b = 0; b < nbodies; ++b) {
595: ego body = bodies[b];
596: int id, Nl, l;
598: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
599: for (l = 0; l < Nl; ++l) {
600: ego loop = lobjs[l];
601: ego *fobjs;
602: int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;
604: lid = EG_indexBodyTopo(body, loop);
605: EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
606: if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
607: fid = EG_indexBodyTopo(body, fobjs[0]);
608: EG_free(fobjs);
609: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
610: for (e = 0; e < Ne; ++e) {
611: ego edge = objs[e];
612: int eid, Nv, v;
613: PetscInt points[3], support[2], numEdges, edgeNum;
614: const PetscInt *edges;
616: eid = EG_indexBodyTopo(body, edge);
617: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
618: if (mtype == DEGENERATE) continue;
619: else ++Ner;
620: for (v = 0; v < Nv; ++v) {
621: ego vertex = nobjs[v];
623: id = EG_indexBodyTopo(body, vertex);
624: DMLabelSetValue(edgeLabel, numCells + id-1, eid);
625: points[v*2] = numCells + id-1;
626: }
627: PetscHMapIGet(edgeMap, eid-1, &edgeNum);
628: points[1] = numCells + numVertices - newVertices + edgeNum;
630: DMLabelSetValue(edgeLabel, points[1], eid);
631: support[0] = points[0];
632: support[1] = points[1];
633: DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
634: if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
635: DMLabelSetValue(edgeLabel, edges[0], eid);
636: DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
637: support[0] = points[1];
638: support[1] = points[2];
639: DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
640: if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges);
641: DMLabelSetValue(edgeLabel, edges[0], eid);
642: DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
643: }
644: switch (Ner) {
645: case 2: Nt = 2;break;
646: case 3: Nt = 4;break;
647: case 4: Nt = 8;break;
648: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
649: }
650: for (t = 0; t < Nt; ++t) {
651: DMLabelSetValue(bodyLabel, cOff+t, b);
652: DMLabelSetValue(faceLabel, cOff+t, fid);
653: }
654: cOff += Nt;
655: }
656: EG_free(lobjs);
657: }
658: PetscHMapIDestroy(&edgeMap);
659: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
660: for (c = cStart; c < cEnd; ++c) {
661: PetscInt *closure = NULL;
662: PetscInt clSize, cl, bval, fval;
664: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
665: DMLabelGetValue(bodyLabel, c, &bval);
666: DMLabelGetValue(faceLabel, c, &fval);
667: for (cl = 0; cl < clSize*2; cl += 2) {
668: DMLabelSetValue(bodyLabel, closure[cl], bval);
669: DMLabelSetValue(faceLabel, closure[cl], fval);
670: }
671: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
672: }
673: *newdm = dm;
674: return(0);
675: }
677: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
678: {
679: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
680: // EGADS/EGADSLite variables
681: ego geom, *bodies, *mobjs, *fobjs, *lobjs, *eobjs, *nobjs;
682: ego topRef, prev, next;
683: int oclass, mtype, nbodies, *senses, *lSenses, *eSenses;
684: int b;
685: // PETSc variables
686: DM dm;
687: PetscHMapI edgeMap = NULL, bodyIndexMap = NULL, bodyVertexMap = NULL, bodyEdgeMap = NULL, bodyFaceMap = NULL, bodyEdgeGlobalMap = NULL;
688: PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numEdges = 0, numFaces = 0, numCells = 0, edgeCntr = 0;
689: PetscInt cellCntr = 0, numPoints = 0;
690: PetscInt *cells = NULL;
691: const PetscInt *cone = NULL;
692: PetscReal *coords = NULL;
693: PetscMPIInt rank;
694: PetscErrorCode ierr;
697: MPI_Comm_rank(comm, &rank);
698: if (!rank) {
699: // ---------------------------------------------------------------------------------------------------
700: // Generate Petsc Plex
701: // Get all Nodes in model, record coordinates in a correctly formatted array
702: // Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
703: // We need to uniformly refine the initial geometry to guarantee a valid mesh
705: // Caluculate cell and vertex sizes
706: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
708: PetscHMapICreate(&edgeMap);
709: PetscHMapICreate(&bodyIndexMap);
710: PetscHMapICreate(&bodyVertexMap);
711: PetscHMapICreate(&bodyEdgeMap);
712: PetscHMapICreate(&bodyEdgeGlobalMap);
713: PetscHMapICreate(&bodyFaceMap);
715: for (b = 0; b < nbodies; ++b) {
716: ego body = bodies[b];
717: int Nf, Ne, Nv;
718: PetscHashIter BIiter, BViter, BEiter, BEGiter, BFiter, EMiter;
719: PetscBool BIfound, BVfound, BEfound, BEGfound, BFfound, EMfound;
721: PetscHMapIFind(bodyIndexMap, b, &BIiter, &BIfound);
722: PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
723: PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
724: PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
725: PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);
727: if (!BIfound) {PetscHMapISet(bodyIndexMap, b, numFaces + numEdges + numVertices);}
728: if (!BVfound) {PetscHMapISet(bodyVertexMap, b, numVertices);}
729: if (!BEfound) {PetscHMapISet(bodyEdgeMap, b, numEdges);}
730: if (!BEGfound) {PetscHMapISet(bodyEdgeGlobalMap, b, edgeCntr);}
731: if (!BFfound) {PetscHMapISet(bodyFaceMap, b, numFaces);}
733: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
734: EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);
735: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
736: EG_free(fobjs);
737: EG_free(eobjs);
738: EG_free(nobjs);
740: // Remove DEGENERATE EDGES from Edge count
741: EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);
742: int Netemp = 0;
743: for (int e = 0; e < Ne; ++e) {
744: ego edge = eobjs[e];
745: int eid;
747: EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
748: eid = EG_indexBodyTopo(body, edge);
750: PetscHMapIFind(edgeMap, edgeCntr + eid - 1, &EMiter, &EMfound);
751: if (mtype == DEGENERATE) {
752: if (!EMfound) {PetscHMapISet(edgeMap, edgeCntr + eid - 1, -1);}
753: }
754: else {
755: ++Netemp;
756: if (!EMfound) {PetscHMapISet(edgeMap, edgeCntr + eid - 1, Netemp);}
757: }
758: }
759: EG_free(eobjs);
761: // Determine Number of Cells
762: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
763: for (int f = 0; f < Nf; ++f) {
764: ego face = fobjs[f];
765: int edgeTemp = 0;
767: EG_getBodyTopos(body, face, EDGE, &Ne, &eobjs);
768: for (int e = 0; e < Ne; ++e) {
769: ego edge = eobjs[e];
771: EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
772: if (mtype != DEGENERATE) {++edgeTemp;}
773: }
774: numCells += (2 * edgeTemp);
775: EG_free(eobjs);
776: }
777: EG_free(fobjs);
779: numFaces += Nf;
780: numEdges += Netemp;
781: numVertices += Nv;
782: edgeCntr += Ne;
783: }
785: // Set up basic DMPlex parameters
786: dim = 2; // Assumes 3D Models :: Need to handle 2D modles in the future
787: cdim = 3; // Assumes 3D Models :: Need to update to handle 2D modles in future
788: numCorners = 3; // Split Faces into triangles
789: numPoints = numVertices + numEdges + numFaces; // total number of coordinate points
791: PetscMalloc2(numPoints*cdim, &coords, numCells*numCorners, &cells);
793: // Get Vertex Coordinates and Set up Cells
794: for (b = 0; b < nbodies; ++b) {
795: ego body = bodies[b];
796: int Nf, Ne, Nv;
797: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
798: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
799: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
801: // Vertices on Current Body
802: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
804: PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
805: if (!BVfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
806: PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart);
808: for (int v = 0; v < Nv; ++v) {
809: ego vertex = nobjs[v];
810: double limits[4];
811: int id, dummy;
813: EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
814: id = EG_indexBodyTopo(body, vertex);
816: coords[(bodyVertexIndexStart + id - 1)*cdim + 0] = limits[0];
817: coords[(bodyVertexIndexStart + id - 1)*cdim + 1] = limits[1];
818: coords[(bodyVertexIndexStart + id - 1)*cdim + 2] = limits[2];
819: }
820: EG_free(nobjs);
822: // Edge Midpoint Vertices on Current Body
823: EG_getBodyTopos(body, NULL, EDGE, &Ne, &eobjs);
825: PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
826: if (!BEfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
827: PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart);
829: PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
830: if (!BEGfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
831: PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart);
833: for (int e = 0; e < Ne; ++e) {
834: ego edge = eobjs[e];
835: double range[2], avgt[1], cntrPnt[9];
836: int eid, eOffset;
837: int periodic;
839: EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
840: if (mtype == DEGENERATE) {continue;}
842: eid = EG_indexBodyTopo(body, edge);
844: // get relative offset from globalEdgeID Vector
845: PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
846: if (!EMfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d not found in edgeMap", bodyEdgeGlobalIndexStart + eid - 1);
847: PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);
849: EG_getRange(edge, range, &periodic);
850: avgt[0] = (range[0] + range[1]) / 2.;
852: EG_evaluate(edge, avgt, cntrPnt);
853: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 0] = cntrPnt[0];
854: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 1] = cntrPnt[1];
855: coords[(numVertices + bodyEdgeIndexStart + eOffset - 1)*cdim + 2] = cntrPnt[2];
856: }
857: EG_free(eobjs);
859: // Face Midpoint Vertices on Current Body
860: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
862: PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);
863: if (!BFfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
864: PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart);
866: for (int f = 0; f < Nf; ++f) {
867: ego face = fobjs[f];
868: double range[4], avgUV[2], cntrPnt[18];
869: int peri, id;
871: id = EG_indexBodyTopo(body, face);
872: EG_getRange(face, range, &peri);
874: avgUV[0] = (range[0] + range[1]) / 2.;
875: avgUV[1] = (range[2] + range[3]) / 2.;
876: EG_evaluate(face, avgUV, cntrPnt);
878: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 0] = cntrPnt[0];
879: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 1] = cntrPnt[1];
880: coords[(numVertices + numEdges + bodyFaceIndexStart + id - 1)*cdim + 2] = cntrPnt[2];
881: }
882: EG_free(fobjs);
884: // Define Cells :: Note - This could be incorporated in the Face Midpoint Vertices Loop but was kept separate for clarity
885: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
886: for (int f = 0; f < Nf; ++f) {
887: ego face = fobjs[f];
888: int fID, midFaceID, midPntID, startID, endID, Nl;
890: fID = EG_indexBodyTopo(body, face);
891: midFaceID = numVertices + numEdges + bodyFaceIndexStart + fID - 1;
892: // Must Traverse Loop to ensure we have all necessary information like the sense (+/- 1) of the edges.
893: // TODO :: Only handles single loop faces (No holes). The choices for handling multiloop faces are:
894: // 1) Use the DMPlexCreateEGADSFromFile() with the -dm_plex_egads_with_tess = 1 option.
895: // This will use a default EGADS tessellation as an initial surface mesh.
896: // 2) Create the initial surface mesh via a 2D mesher :: Currently not availble (?future?)
897: // May I suggest the XXXX as a starting point?
899: EG_getTopology(face, &geom, &oclass, &mtype, NULL, &Nl, &lobjs, &lSenses);
901: if (Nl > 1) SETERRQ1(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);
902: for (int l = 0; l < Nl; ++l) {
903: ego loop = lobjs[l];
905: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses);
906: for (int e = 0; e < Ne; ++e) {
907: ego edge = eobjs[e];
908: int eid, eOffset;
910: EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
911: eid = EG_indexBodyTopo(body, edge);
912: if (mtype == DEGENERATE) { continue; }
914: // get relative offset from globalEdgeID Vector
915: PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
916: if (!EMfound) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
917: PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);
919: midPntID = numVertices + bodyEdgeIndexStart + eOffset - 1;
921: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
923: if (eSenses[e] > 0) { startID = EG_indexBodyTopo(body, nobjs[0]); endID = EG_indexBodyTopo(body, nobjs[1]); }
924: else { startID = EG_indexBodyTopo(body, nobjs[1]); endID = EG_indexBodyTopo(body, nobjs[0]); }
926: // Define 2 Cells per Edge with correct orientation
927: cells[cellCntr*numCorners + 0] = midFaceID;
928: cells[cellCntr*numCorners + 1] = bodyVertexIndexStart + startID - 1;
929: cells[cellCntr*numCorners + 2] = midPntID;
931: cells[cellCntr*numCorners + 3] = midFaceID;
932: cells[cellCntr*numCorners + 4] = midPntID;
933: cells[cellCntr*numCorners + 5] = bodyVertexIndexStart + endID - 1;
935: cellCntr = cellCntr + 2;
936: }
937: }
938: }
939: EG_free(fobjs);
940: }
941: }
943: // Generate DMPlex
944: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
945: PetscFree2(coords, cells);
946: PetscInfo1(dm, " Total Number of Unique Cells = %D \n", numCells);
947: PetscInfo1(dm, " Total Number of Unique Vertices = %D \n", numVertices);
949: // Embed EGADS model in DM
950: {
951: PetscContainer modelObj, contextObj;
953: PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
954: PetscContainerSetPointer(modelObj, model);
955: PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
956: PetscContainerDestroy(&modelObj);
958: PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
959: PetscContainerSetPointer(contextObj, context);
960: PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
961: PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
962: PetscContainerDestroy(&contextObj);
963: }
964: // Label points
965: PetscInt nStart, nEnd;
967: DMCreateLabel(dm, "EGADS Body ID");
968: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
969: DMCreateLabel(dm, "EGADS Face ID");
970: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
971: DMCreateLabel(dm, "EGADS Edge ID");
972: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
973: DMCreateLabel(dm, "EGADS Vertex ID");
974: DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
976: DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd);
978: cellCntr = 0;
979: for (b = 0; b < nbodies; ++b) {
980: ego body = bodies[b];
981: int Nv, Ne, Nf;
982: PetscInt bodyVertexIndexStart, bodyEdgeIndexStart, bodyEdgeGlobalIndexStart, bodyFaceIndexStart;
983: PetscHashIter BViter, BEiter, BEGiter, BFiter, EMiter;
984: PetscBool BVfound, BEfound, BEGfound, BFfound, EMfound;
986: PetscHMapIFind(bodyVertexMap, b, &BViter, &BVfound);
987: if (!BVfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyVertexMap", b);
988: PetscHMapIGet(bodyVertexMap, b, &bodyVertexIndexStart);
990: PetscHMapIFind(bodyEdgeMap, b, &BEiter, &BEfound);
991: if (!BEfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeMap", b);
992: PetscHMapIGet(bodyEdgeMap, b, &bodyEdgeIndexStart);
994: PetscHMapIFind(bodyFaceMap, b, &BFiter, &BFfound);
995: if (!BFfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyFaceMap", b);
996: PetscHMapIGet(bodyFaceMap, b, &bodyFaceIndexStart);
998: PetscHMapIFind(bodyEdgeGlobalMap, b, &BEGiter, &BEGfound);
999: if (!BEGfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in bodyEdgeGlobalMap", b);
1000: PetscHMapIGet(bodyEdgeGlobalMap, b, &bodyEdgeGlobalIndexStart);
1002: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
1003: for (int f = 0; f < Nf; ++f) {
1004: ego face = fobjs[f];
1005: int fID, Nl;
1007: fID = EG_indexBodyTopo(body, face);
1009: EG_getBodyTopos(body, face, LOOP, &Nl, &lobjs);
1010: for (int l = 0; l < Nl; ++l) {
1011: ego loop = lobjs[l];
1012: int lid;
1014: lid = EG_indexBodyTopo(body, loop);
1015: if (Nl > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
1017: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &eobjs, &eSenses);
1018: for (int e = 0; e < Ne; ++e) {
1019: ego edge = eobjs[e];
1020: int eid, eOffset;
1022: // Skip DEGENERATE Edges
1023: EG_getInfo(edge, &oclass, &mtype, &topRef, &prev, &next);
1024: if (mtype == DEGENERATE) {continue;}
1025: eid = EG_indexBodyTopo(body, edge);
1027: // get relative offset from globalEdgeID Vector
1028: PetscHMapIFind(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &EMiter, &EMfound);
1029: if (!EMfound) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_SUP, "Edge %d of Body %d not found in edgeMap. Global Edge ID :: %d", eid, b, bodyEdgeGlobalIndexStart + eid - 1);
1030: PetscHMapIGet(edgeMap, bodyEdgeGlobalIndexStart + eid - 1, &eOffset);
1032: EG_getBodyTopos(body, edge, NODE, &Nv, &nobjs);
1033: for (int v = 0; v < Nv; ++v){
1034: ego vertex = nobjs[v];
1035: int vID;
1037: vID = EG_indexBodyTopo(body, vertex);
1038: DMLabelSetValue(bodyLabel, nStart + bodyVertexIndexStart + vID - 1, b);
1039: DMLabelSetValue(vertexLabel, nStart + bodyVertexIndexStart + vID - 1, vID);
1040: }
1041: EG_free(nobjs);
1043: DMLabelSetValue(bodyLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, b);
1044: DMLabelSetValue(edgeLabel, nStart + numVertices + bodyEdgeIndexStart + eOffset - 1, eid);
1046: // Define Cell faces
1047: for (int jj = 0; jj < 2; ++jj){
1048: DMLabelSetValue(bodyLabel, cellCntr, b);
1049: DMLabelSetValue(faceLabel, cellCntr, fID);
1050: DMPlexGetCone(dm, cellCntr, &cone);
1052: DMLabelSetValue(bodyLabel, cone[0], b);
1053: DMLabelSetValue(faceLabel, cone[0], fID);
1055: DMLabelSetValue(bodyLabel, cone[1], b);
1056: DMLabelSetValue(edgeLabel, cone[1], eid);
1058: DMLabelSetValue(bodyLabel, cone[2], b);
1059: DMLabelSetValue(faceLabel, cone[2], fID);
1061: cellCntr = cellCntr + 1;
1062: }
1063: }
1064: }
1065: EG_free(lobjs);
1067: DMLabelSetValue(bodyLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, b);
1068: DMLabelSetValue(faceLabel, nStart + numVertices + numEdges + bodyFaceIndexStart + fID - 1, fID);
1069: }
1070: EG_free(fobjs);
1071: }
1073: PetscHMapIDestroy(&edgeMap);
1074: PetscHMapIDestroy(&bodyIndexMap);
1075: PetscHMapIDestroy(&bodyVertexMap);
1076: PetscHMapIDestroy(&bodyEdgeMap);
1077: PetscHMapIDestroy(&bodyEdgeGlobalMap);
1078: PetscHMapIDestroy(&bodyFaceMap);
1080: *newdm = dm;
1081: return(0);
1082: }
1084: static PetscErrorCode DMPlexCreateEGADS_Tess_Internal(MPI_Comm comm, ego context, ego model, DM *newdm)
1085: {
1086: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1087: /* EGADSLite variables */
1088: ego geom, *bodies, *fobjs;
1089: int b, oclass, mtype, nbodies, *senses;
1090: int totalNumTris = 0, totalNumPoints = 0;
1091: double boundBox[6] = {0., 0., 0., 0., 0., 0.}, tessSize;
1092: /* PETSc variables */
1093: DM dm;
1094: PetscHMapI pointIndexStartMap = NULL, triIndexStartMap = NULL, pTypeLabelMap = NULL, pIndexLabelMap = NULL;
1095: PetscHMapI pBodyIndexLabelMap = NULL, triFaceIDLabelMap = NULL, triBodyIDLabelMap = NULL;
1096: PetscInt dim = -1, cdim = -1, numCorners = 0, counter = 0;
1097: PetscInt *cells = NULL;
1098: const PetscInt *cone = NULL;
1099: PetscReal *coords = NULL;
1100: PetscMPIInt rank;
1101: PetscErrorCode ierr;
1104: MPI_Comm_rank(comm, &rank);
1105: if (!rank) {
1106: // ---------------------------------------------------------------------------------------------------
1107: // Generate Petsc Plex from EGADSlite created Tessellation of geometry
1108: // ---------------------------------------------------------------------------------------------------
1110: // Caluculate cell and vertex sizes
1111: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
1113: PetscHMapICreate(&pointIndexStartMap);
1114: PetscHMapICreate(&triIndexStartMap);
1115: PetscHMapICreate(&pTypeLabelMap);
1116: PetscHMapICreate(&pIndexLabelMap);
1117: PetscHMapICreate(&pBodyIndexLabelMap);
1118: PetscHMapICreate(&triFaceIDLabelMap);
1119: PetscHMapICreate(&triBodyIDLabelMap);
1121: /* Create Tessellation of Bodies */
1122: ego tessArray[nbodies];
1124: for (b = 0; b < nbodies; ++b) {
1125: ego body = bodies[b];
1126: double params[3] = {0.0, 0.0, 0.0}; // Parameters for Tessellation
1127: int Nf, bodyNumPoints = 0, bodyNumTris = 0;
1128: PetscHashIter PISiter, TISiter;
1129: PetscBool PISfound, TISfound;
1131: /* Store Start Index for each Body's Point and Tris */
1132: PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound);
1133: PetscHMapIFind(triIndexStartMap, b, &TISiter, &TISfound);
1135: if (!PISfound) {PetscHMapISet(pointIndexStartMap, b, totalNumPoints);}
1136: if (!TISfound) {PetscHMapISet(triIndexStartMap, b, totalNumTris);}
1138: /* Calculate Tessellation parameters based on Bounding Box */
1139: /* Get Bounding Box Dimensions of the BODY */
1140: EG_getBoundingBox(body, boundBox);
1141: tessSize = boundBox[3] - boundBox[0];
1142: if (tessSize < boundBox[4] - boundBox[1]) tessSize = boundBox[4] - boundBox[1];
1143: if (tessSize < boundBox[5] - boundBox[2]) tessSize = boundBox[5] - boundBox[2];
1145: // TODO :: May want to give users tessellation parameter options //
1146: params[0] = 0.0250 * tessSize;
1147: params[1] = 0.0075 * tessSize;
1148: params[2] = 15.0;
1150: EG_makeTessBody(body, params, &tessArray[b]);
1152: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
1154: for (int f = 0; f < Nf; ++f) {
1155: ego face = fobjs[f];
1156: int len, fID, ntris;
1157: const int *ptype, *pindex, *ptris, *ptric;
1158: const double *pxyz, *puv;
1160: // Get Face ID //
1161: fID = EG_indexBodyTopo(body, face);
1163: // Checkout the Surface Tessellation //
1164: EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric);
1166: // Determine total number of triangle cells in the tessellation //
1167: bodyNumTris += (int) ntris;
1169: // Check out the point index and coordinate //
1170: for (int p = 0; p < len; ++p) {
1171: int global;
1173: EG_localToGlobal(tessArray[b], fID, p+1, &global);
1175: // Determine the total number of points in the tessellation //
1176: bodyNumPoints = PetscMax(bodyNumPoints, global);
1177: }
1178: }
1179: EG_free(fobjs);
1181: totalNumPoints += bodyNumPoints;
1182: totalNumTris += bodyNumTris;
1183: }
1184: //} - Original End of (!rank)
1186: dim = 2;
1187: cdim = 3;
1188: numCorners = 3;
1189: //PetscInt counter = 0;
1191: /* NEED TO DEFINE MATRICES/VECTORS TO STORE GEOM REFERENCE DATA */
1192: /* Fill in below and use to define DMLabels after DMPlex creation */
1193: PetscMalloc2(totalNumPoints*cdim, &coords, totalNumTris*numCorners, &cells);
1195: for (b = 0; b < nbodies; ++b) {
1196: ego body = bodies[b];
1197: int Nf;
1198: PetscInt pointIndexStart;
1199: PetscHashIter PISiter;
1200: PetscBool PISfound;
1202: PetscHMapIFind(pointIndexStartMap, b, &PISiter, &PISfound);
1203: if (!PISfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Body %d not found in pointIndexStartMap", b);
1204: PetscHMapIGet(pointIndexStartMap, b, &pointIndexStart);
1206: EG_getBodyTopos(body, NULL, FACE, &Nf, &fobjs);
1208: for (int f = 0; f < Nf; ++f) {
1209: /* Get Face Object */
1210: ego face = fobjs[f];
1211: int len, fID, ntris;
1212: const int *ptype, *pindex, *ptris, *ptric;
1213: const double *pxyz, *puv;
1215: /* Get Face ID */
1216: fID = EG_indexBodyTopo(body, face);
1218: /* Checkout the Surface Tessellation */
1219: EG_getTessFace(tessArray[b], fID, &len, &pxyz, &puv, &ptype, &pindex, &ntris, &ptris, &ptric);
1221: /* Check out the point index and coordinate */
1222: for (int p = 0; p < len; ++p) {
1223: int global;
1224: PetscHashIter PTLiter, PILiter, PBLiter;
1225: PetscBool PTLfound, PILfound, PBLfound;
1227: EG_localToGlobal(tessArray[b], fID, p+1, &global);
1229: /* Set the coordinates array for DAG */
1230: coords[((global-1+pointIndexStart)*3) + 0] = pxyz[(p*3)+0];
1231: coords[((global-1+pointIndexStart)*3) + 1] = pxyz[(p*3)+1];
1232: coords[((global-1+pointIndexStart)*3) + 2] = pxyz[(p*3)+2];
1234: /* Store Geometry Label Information for DMLabel assignment later */
1235: PetscHMapIFind(pTypeLabelMap, global-1+pointIndexStart, &PTLiter, &PTLfound);
1236: PetscHMapIFind(pIndexLabelMap, global-1+pointIndexStart, &PILiter, &PILfound);
1237: PetscHMapIFind(pBodyIndexLabelMap, global-1+pointIndexStart, &PBLiter, &PBLfound);
1239: if (!PTLfound) {PetscHMapISet(pTypeLabelMap, global-1+pointIndexStart, ptype[p]);}
1240: if (!PILfound) {PetscHMapISet(pIndexLabelMap, global-1+pointIndexStart, pindex[p]);}
1241: if (!PBLfound) {PetscHMapISet(pBodyIndexLabelMap, global-1+pointIndexStart, b);}
1243: if (ptype[p] < 0) { PetscHMapISet(pIndexLabelMap, global-1+pointIndexStart, fID);}
1244: }
1246: for (int t = 0; t < (int) ntris; ++t){
1247: int global, globalA, globalB;
1248: PetscHashIter TFLiter, TBLiter;
1249: PetscBool TFLfound, TBLfound;
1251: EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 0], &global);
1252: cells[(counter*3) +0] = global-1+pointIndexStart;
1254: EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 1], &globalA);
1255: cells[(counter*3) +1] = globalA-1+pointIndexStart;
1257: EG_localToGlobal(tessArray[b], fID, ptris[(t*3) + 2], &globalB);
1258: cells[(counter*3) +2] = globalB-1+pointIndexStart;
1260: PetscHMapIFind(triFaceIDLabelMap, counter, &TFLiter, &TFLfound);
1261: PetscHMapIFind(triBodyIDLabelMap, counter, &TBLiter, &TBLfound);
1263: if (!TFLfound) {PetscHMapISet(triFaceIDLabelMap, counter, fID);}
1264: if (!TBLfound) {PetscHMapISet(triBodyIDLabelMap, counter, b);}
1266: counter += 1;
1267: }
1268: }
1269: EG_free(fobjs);
1270: }
1271: }
1273: //Build DMPlex
1274: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, totalNumTris, totalNumPoints, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
1275: PetscFree2(coords, cells);
1277: // Embed EGADS model in DM
1278: {
1279: PetscContainer modelObj, contextObj;
1281: PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
1282: PetscContainerSetPointer(modelObj, model);
1283: PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
1284: PetscContainerDestroy(&modelObj);
1286: PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
1287: PetscContainerSetPointer(contextObj, context);
1288: PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
1289: PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
1290: PetscContainerDestroy(&contextObj);
1291: }
1293: // Label Points
1294: DMCreateLabel(dm, "EGADS Body ID");
1295: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
1296: DMCreateLabel(dm, "EGADS Face ID");
1297: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
1298: DMCreateLabel(dm, "EGADS Edge ID");
1299: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
1300: DMCreateLabel(dm, "EGADS Vertex ID");
1301: DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
1303: /* Get Number of DAG Nodes at each level */
1304: int fStart, fEnd, eStart, eEnd, nStart, nEnd;
1306: DMPlexGetHeightStratum(dm, 0, &fStart, &fEnd);
1307: DMPlexGetHeightStratum(dm, 1, &eStart, &eEnd);
1308: DMPlexGetHeightStratum(dm, 2, &nStart, &nEnd);
1310: /* Set DMLabels for NODES */
1311: for (int n = nStart; n < nEnd; ++n) {
1312: int pTypeVal, pIndexVal, pBodyVal;
1313: PetscHashIter PTLiter, PILiter, PBLiter;
1314: PetscBool PTLfound, PILfound, PBLfound;
1316: //Converted to Hash Tables
1317: PetscHMapIFind(pTypeLabelMap, n - nStart, &PTLiter, &PTLfound);
1318: if (!PTLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pTypeLabelMap", n);
1319: PetscHMapIGet(pTypeLabelMap, n - nStart, &pTypeVal);
1321: PetscHMapIFind(pIndexLabelMap, n - nStart, &PILiter, &PILfound);
1322: if (!PILfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pIndexLabelMap", n);
1323: PetscHMapIGet(pIndexLabelMap, n - nStart, &pIndexVal);
1325: PetscHMapIFind(pBodyIndexLabelMap, n - nStart, &PBLiter, &PBLfound);
1326: if (!PBLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in pBodyLabelMap", n);
1327: PetscHMapIGet(pBodyIndexLabelMap, n - nStart, &pBodyVal);
1329: DMLabelSetValue(bodyLabel, n, pBodyVal);
1330: if (pTypeVal == 0) {DMLabelSetValue(vertexLabel, n, pIndexVal);}
1331: if (pTypeVal > 0) {DMLabelSetValue(edgeLabel, n, pIndexVal);}
1332: if (pTypeVal < 0) {DMLabelSetValue(faceLabel, n, pIndexVal);}
1333: }
1335: /* Set DMLabels for Edges - Based on the DMLabels of the EDGE's NODES */
1336: for (int e = eStart; e < eEnd; ++e) {
1337: int bodyID_0, vertexID_0, vertexID_1, edgeID_0, edgeID_1, faceID_0, faceID_1;
1339: DMPlexGetCone(dm, e, &cone);
1340: DMLabelGetValue(bodyLabel, cone[0], &bodyID_0); // Do I need to check the other end?
1341: DMLabelGetValue(vertexLabel, cone[0], &vertexID_0);
1342: DMLabelGetValue(vertexLabel, cone[1], &vertexID_1);
1343: DMLabelGetValue(edgeLabel, cone[0], &edgeID_0);
1344: DMLabelGetValue(edgeLabel, cone[1], &edgeID_1);
1345: DMLabelGetValue(faceLabel, cone[0], &faceID_0);
1346: DMLabelGetValue(faceLabel, cone[1], &faceID_1);
1348: DMLabelSetValue(bodyLabel, e, bodyID_0);
1350: if (edgeID_0 == edgeID_1) { DMLabelSetValue(edgeLabel, e, edgeID_0); }
1351: else if (vertexID_0 > 0 && edgeID_1 > 0) { DMLabelSetValue(edgeLabel, e, edgeID_1); }
1352: else if (vertexID_1 > 0 && edgeID_0 > 0) { DMLabelSetValue(edgeLabel, e, edgeID_0); }
1353: else { /* Do Nothing */ }
1354: }
1356: /* Set DMLabels for Cells */
1357: for (int f = fStart; f < fEnd; ++f){
1358: int edgeID_0;
1359: PetscInt triBodyVal, triFaceVal;
1360: PetscHashIter TFLiter, TBLiter;
1361: PetscBool TFLfound, TBLfound;
1363: // Convert to Hash Table
1364: PetscHMapIFind(triFaceIDLabelMap, f - fStart, &TFLiter, &TFLfound);
1365: if (!TFLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triFaceIDLabelMap", f);
1366: PetscHMapIGet(triFaceIDLabelMap, f - fStart, &triFaceVal);
1368: PetscHMapIFind(triBodyIDLabelMap, f - fStart, &TBLiter, &TBLfound);
1369: if (!TBLfound) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "DAG Point %d not found in triBodyIDLabelMap", f);
1370: PetscHMapIGet(triBodyIDLabelMap, f - fStart, &triBodyVal);
1372: DMLabelSetValue(bodyLabel, f, triBodyVal);
1373: DMLabelSetValue(faceLabel, f, triFaceVal);
1375: /* Finish Labeling previously unlabeled DMPlex Edges - Assumes Triangular Cell (3 Edges Max) */
1376: DMPlexGetCone(dm, f, &cone);
1378: for (int jj = 0; jj < 3; ++jj) {
1379: DMLabelGetValue(edgeLabel, cone[jj], &edgeID_0);
1381: if (edgeID_0 < 0) {
1382: DMLabelSetValue(bodyLabel, cone[jj], triBodyVal);
1383: DMLabelSetValue(faceLabel, cone[jj], triFaceVal);
1384: }
1385: }
1386: }
1388: *newdm = dm;
1389: return(0);
1390: }
1391: #endif
1393: /*@
1394: DMPlexInflateToGeomModel - Snaps the vertex coordinates of a DMPlex object representing the mesh to its geometry if some vertices depart from the model. This usually happens with non-conforming refinement.
1396: Collective on dm
1398: Input Parameter:
1399: . dm - The uninflated DM object representing the mesh
1401: Output Parameter:
1402: . dm - The inflated DM object representing the mesh
1404: Level: intermediate
1406: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS()
1407: @*/
1408: PetscErrorCode DMPlexInflateToGeomModel(DM dm)
1409: {
1410: #if defined(PETSC_HAVE_EGADS)
1411: /* EGADS Variables */
1412: ego model, geom, body, face, edge;
1413: ego *bodies;
1414: int Nb, oclass, mtype, *senses;
1415: double result[3];
1416: /* PETSc Variables */
1417: DM cdm;
1418: PetscContainer modelObj;
1419: DMLabel bodyLabel, faceLabel, edgeLabel, vertexLabel;
1420: Vec coordinates;
1421: PetscScalar *coords;
1422: PetscInt bodyID, faceID, edgeID, vertexID;
1423: PetscInt cdim, d, vStart, vEnd, v;
1425: #endif
1428: #if defined(PETSC_HAVE_EGADS)
1429: PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
1430: if (!modelObj) return(0);
1431: DMGetCoordinateDim(dm, &cdim);
1432: DMGetCoordinateDM(dm, &cdm);
1433: DMGetCoordinatesLocal(dm, &coordinates);
1434: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
1435: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
1436: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
1437: DMGetLabel(dm, "EGADS Vertex ID", &vertexLabel);
1439: PetscContainerGetPointer(modelObj, (void **) &model);
1440: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
1442: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1443: VecGetArrayWrite(coordinates, &coords);
1444: for (v = vStart; v < vEnd; ++v) {
1445: PetscScalar *vcoords;
1447: DMLabelGetValue(bodyLabel, v, &bodyID);
1448: DMLabelGetValue(faceLabel, v, &faceID);
1449: DMLabelGetValue(edgeLabel, v, &edgeID);
1450: DMLabelGetValue(vertexLabel, v, &vertexID);
1452: if (bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
1453: body = bodies[bodyID];
1455: DMPlexPointLocalRef(cdm, v, coords, (void *) &vcoords);
1456: if (edgeID > 0) {
1457: /* Snap to EDGE at nearest location */
1458: double params[1];
1459: EG_objectBodyTopo(body, EDGE, edgeID, &edge);
1460: EG_invEvaluate(edge, vcoords, params, result); // Get (x,y,z) of nearest point on EDGE
1461: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1462: } else if (faceID > 0) {
1463: /* Snap to FACE at nearest location */
1464: double params[2];
1465: EG_objectBodyTopo(body, FACE, faceID, &face);
1466: EG_invEvaluate(face, vcoords, params, result); // Get (x,y,z) of nearest point on FACE
1467: for (d = 0; d < cdim; ++d) vcoords[d] = result[d];
1468: }
1469: }
1470: VecRestoreArrayWrite(coordinates, &coords);
1471: /* Clear out global coordinates */
1472: VecDestroy(&dm->coordinates);
1473: #endif
1474: return(0);
1475: }
1477: /*@C
1478: DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADS, IGES, or STEP file.
1480: Collective
1482: Input Parameters:
1483: + comm - The MPI communicator
1484: - filename - The name of the EGADS, IGES, or STEP file
1486: Output Parameter:
1487: . dm - The DM object representing the mesh
1489: Level: beginner
1491: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS(), DMPlexCreateEGADSLiteFromFile()
1492: @*/
1493: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
1494: {
1495: PetscMPIInt rank;
1496: #if defined(PETSC_HAVE_EGADS)
1497: ego context= NULL, model = NULL;
1498: #endif
1499: PetscBool printModel = PETSC_FALSE, tessModel = PETSC_FALSE, newModel = PETSC_FALSE;
1504: PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
1505: PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_tess_model", &tessModel, NULL);
1506: PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_new_model", &newModel, NULL);
1507: MPI_Comm_rank(comm, &rank);
1508: #if defined(PETSC_HAVE_EGADS)
1509: if (rank == 0) {
1511: EG_open(&context);
1512: EG_loadModel(context, 0, filename, &model);
1513: if (printModel) {DMPlexEGADSPrintModel_Internal(model);}
1515: }
1516: if (tessModel) {DMPlexCreateEGADS_Tess_Internal(comm, context, model, dm);}
1517: else if (newModel) {DMPlexCreateEGADS(comm, context, model, dm);}
1518: else {DMPlexCreateEGADS_Internal(comm, context, model, dm);}
1519: return(0);
1520: #else
1521: SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADS support. Reconfigure using --download-egads");
1522: #endif
1523: }