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: */
22: /*@
23: 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.
25: Not collective
27: Input Parameters:
28: + dm - The DMPlex object
29: . p - The mesh point
30: - mcoords - A coordinate point lying on the mesh point
32: Output Parameter:
33: . gcoords - The closest coordinate point on the geometry model associated with 'p' to the given point
35: Note: Returns the original coordinates if no geometry model is found. Right now the only supported geometry model is EGADS.
37: Level: intermediate
39: .seealso: DMRefine(), DMPlexCreate(), DMPlexSetRefinementUniform()
40: @*/
41: PetscErrorCode DMPlexSnapToGeomModel(DM dm, PetscInt p, const PetscScalar mcoords[], PetscScalar gcoords[])
42: {
43: #ifdef PETSC_HAVE_EGADS
44: DM cdm;
45: DMLabel bodyLabel, faceLabel, edgeLabel;
46: PetscContainer modelObj;
47: PetscInt bodyID, faceID, edgeID;
48: ego *bodies;
49: ego model, geom, body, obj;
50: /* result has to hold derviatives, along with the value */
51: double params[3], result[18], paramsV[16*3], resultV[16*3];
52: int Nb, oclass, mtype, *senses;
53: Vec coordinatesLocal;
54: PetscScalar *coords = NULL;
55: PetscInt Nv, v, Np = 0, pm;
56: #endif
57: PetscInt dE, d;
61: DMGetCoordinateDim(dm, &dE);
62: #ifdef PETSC_HAVE_EGADS
63: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
64: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
65: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
66: if (!bodyLabel || !faceLabel || !edgeLabel) {
67: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
68: return(0);
69: }
70: DMGetCoordinateDM(dm, &cdm);
71: DMGetCoordinatesLocal(dm, &coordinatesLocal);
72: DMLabelGetValue(bodyLabel, p, &bodyID);
73: DMLabelGetValue(faceLabel, p, &faceID);
74: DMLabelGetValue(edgeLabel, p, &edgeID);
75: PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);
76: PetscContainerGetPointer(modelObj, (void **) &model);
77: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
78: if (bodyID < 0 || bodyID >= Nb) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Body %D is not in [0, %d)", bodyID, Nb);
79: body = bodies[bodyID];
81: if (edgeID >= 0) {EG_objectBodyTopo(body, EDGE, edgeID, &obj); Np = 1;}
82: else if (faceID >= 0) {EG_objectBodyTopo(body, FACE, faceID, &obj); Np = 2;}
83: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D is not in edge or face label for EGADS", p);
84: /* Calculate parameters (t or u,v) for vertices */
85: DMPlexVecGetClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
86: Nv /= dE;
87: if (Nv == 1) {
88: DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
89: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
90: return(0);
91: }
92: if (Nv > 16) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %D coordinates associated to point %D", Nv, p);
93: for (v = 0; v < Nv; ++v) {EG_invEvaluate(obj, &coords[v*dE], ¶msV[v*3], &resultV[v*3]);}
94: DMPlexVecRestoreClosure(cdm, NULL, coordinatesLocal, p, &Nv, &coords);
95: /* Calculate parameters (t or u,v) for new vertex at edge midpoint */
96: for (pm = 0; pm < Np; ++pm) {
97: params[pm] = 0.;
98: for (v = 0; v < Nv; ++v) {params[pm] += paramsV[v*3+pm];}
99: params[pm] /= Nv;
100: }
101: /* TODO Check
102: double range[4]; // [umin, umax, vmin, vmax]
103: int peri;
104: EG_getRange(face, range, &peri);
105: if ((paramsNew[0] < range[0]) || (paramsNew[0] > range[1]) || (paramsNew[1] < range[2]) || (paramsNew[1] > range[3])) SETERRQ();
106: */
107: /* Put coordinates for new vertex in result[] */
108: EG_evaluate(obj, params, result);
109: for (d = 0; d < dE; ++d) gcoords[d] = result[d];
110: #else
111: for (d = 0; d < dE; ++d) gcoords[d] = mcoords[d];
112: #endif
113: return(0);
114: }
116: #if defined(PETSC_HAVE_EGADS)
117: static PetscErrorCode DMPlexEGADSPrintModel(ego model)
118: {
119: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
120: int oclass, mtype, *senses;
121: int Nb, b;
125: /* test bodyTopo functions */
126: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);
127: PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb);
129: for (b = 0; b < Nb; ++b) {
130: ego body = bodies[b];
131: int id, Nsh, Nf, Nl, l, Ne, e, Nv, v;
133: /* Output Basic Model Topology */
134: EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);
135: PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh);
136: EG_free(objs);
138: EG_getBodyTopos(body, NULL, FACE, &Nf, &objs);
139: PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf);
140: EG_free(objs);
142: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
143: PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl);
145: EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs);
146: PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne);
147: EG_free(objs);
149: EG_getBodyTopos(body, NULL, NODE, &Nv, &objs);
150: PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv);
151: EG_free(objs);
153: for (l = 0; l < Nl; ++l) {
154: ego loop = lobjs[l];
156: id = EG_indexBodyTopo(body, loop);
157: PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id);
159: /* Get EDGE info which associated with the current LOOP */
160: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
162: for (e = 0; e < Ne; ++e) {
163: ego edge = objs[e];
164: double range[4] = {0., 0., 0., 0.};
165: double point[3] = {0., 0., 0.};
166: double params[3] = {0., 0., 0.};
167: double result[18];
168: int peri;
170: id = EG_indexBodyTopo(body, edge);
171: PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d (%d)\n", id, e);
173: EG_getRange(edge, range, &peri);
174: PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]);
176: /* Get NODE info which associated with the current EDGE */
177: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
178: if (mtype == DEGENERATE) {
179: PetscPrintf(PETSC_COMM_SELF, " EDGE %d is DEGENERATE \n", id);
180: } else {
181: params[0] = range[0];
182: EG_evaluate(edge, params, result);
183: PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]);
184: params[0] = range[1];
185: EG_evaluate(edge, params, result);
186: PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n", result[0], result[1], result[2]);
187: }
189: for (v = 0; v < Nv; ++v) {
190: ego vertex = nobjs[v];
191: double limits[4];
192: int dummy;
194: EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
195: id = EG_indexBodyTopo(body, vertex);
196: PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);
197: PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]);
199: point[0] = point[0] + limits[0];
200: point[1] = point[1] + limits[1];
201: point[2] = point[2] + limits[2];
202: }
203: }
204: }
205: EG_free(lobjs);
206: }
207: return(0);
208: }
210: static PetscErrorCode DMPlexEGADSDestroy_Private(void *context)
211: {
212: if (context) EG_close((ego) context);
213: return 0;
214: }
216: static PetscErrorCode DMPlexCreateEGADS(MPI_Comm comm, ego context, ego model, DM *newdm)
217: {
218: DMLabel bodyLabel, faceLabel, edgeLabel;
219: PetscInt cStart, cEnd, c;
220: /* EGADSLite variables */
221: ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs;
222: int oclass, mtype, nbodies, *senses;
223: int b;
224: /* PETSc variables */
225: DM dm;
226: PetscHMapI edgeMap = NULL;
227: 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;
228: PetscInt *cells = NULL, *cone = NULL;
229: PetscReal *coords = NULL;
230: PetscMPIInt rank;
234: MPI_Comm_rank(comm, &rank);
235: if (!rank) {
236: const PetscInt debug = 0;
238: /* ---------------------------------------------------------------------------------------------------
239: Generate Petsc Plex
240: Get all Nodes in model, record coordinates in a correctly formatted array
241: Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array
242: We need to uniformly refine the initial geometry to guarantee a valid mesh
243: */
245: /* Calculate cell and vertex sizes */
246: EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);
247: PetscHMapICreate(&edgeMap);
248: numEdges = 0;
249: for (b = 0; b < nbodies; ++b) {
250: ego body = bodies[b];
251: int id, Nl, l, Nv, v;
253: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
254: for (l = 0; l < Nl; ++l) {
255: ego loop = lobjs[l];
256: int Ner = 0, Ne, e, Nc;
258: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
259: for (e = 0; e < Ne; ++e) {
260: ego edge = objs[e];
261: int Nv, id;
262: PetscHashIter iter;
263: PetscBool found;
265: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
266: if (mtype == DEGENERATE) continue;
267: id = EG_indexBodyTopo(body, edge);
268: PetscHMapIFind(edgeMap, id-1, &iter, &found);
269: if (!found) {PetscHMapISet(edgeMap, id-1, numEdges++);}
270: ++Ner;
271: }
272: if (Ner == 2) {Nc = 2;}
273: else if (Ner == 3) {Nc = 4;}
274: else if (Ner == 4) {Nc = 8; ++numQuads;}
275: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner);
276: numCells += Nc;
277: newCells += Nc-1;
278: maxCorners = PetscMax(Ner*2+1, maxCorners);
279: }
280: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
281: for (v = 0; v < Nv; ++v) {
282: ego vertex = nobjs[v];
284: id = EG_indexBodyTopo(body, vertex);
285: /* TODO: Instead of assuming contiguous ids, we could use a hash table */
286: numVertices = PetscMax(id, numVertices);
287: }
288: EG_free(lobjs);
289: EG_free(nobjs);
290: }
291: PetscHMapIGetSize(edgeMap, &numEdges);
292: newVertices = numEdges + numQuads;
293: numVertices += newVertices;
295: dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
296: cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */
297: numCorners = 3; /* Split cells into triangles */
298: PetscMalloc3(numVertices*cdim, &coords, numCells*numCorners, &cells, maxCorners, &cone);
300: /* Get vertex coordinates */
301: for (b = 0; b < nbodies; ++b) {
302: ego body = bodies[b];
303: int id, Nv, v;
305: EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);
306: for (v = 0; v < Nv; ++v) {
307: ego vertex = nobjs[v];
308: double limits[4];
309: int dummy;
311: EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);
312: id = EG_indexBodyTopo(body, vertex);
313: coords[(id-1)*cdim+0] = limits[0];
314: coords[(id-1)*cdim+1] = limits[1];
315: coords[(id-1)*cdim+2] = limits[2];
316: }
317: EG_free(nobjs);
318: }
319: PetscHMapIClear(edgeMap);
320: fOff = numVertices - newVertices + numEdges;
321: numEdges = 0;
322: numQuads = 0;
323: for (b = 0; b < nbodies; ++b) {
324: ego body = bodies[b];
325: int Nl, l;
327: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
328: for (l = 0; l < Nl; ++l) {
329: ego loop = lobjs[l];
330: int lid, Ner = 0, Ne, e;
332: lid = EG_indexBodyTopo(body, loop);
333: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
334: for (e = 0; e < Ne; ++e) {
335: ego edge = objs[e];
336: int eid, Nv;
337: PetscHashIter iter;
338: PetscBool found;
340: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
341: if (mtype == DEGENERATE) continue;
342: ++Ner;
343: eid = EG_indexBodyTopo(body, edge);
344: PetscHMapIFind(edgeMap, eid-1, &iter, &found);
345: if (!found) {
346: PetscInt v = numVertices - newVertices + numEdges;
347: double range[4], params[3] = {0., 0., 0.}, result[18];
348: int periodic[2];
350: PetscHMapISet(edgeMap, eid-1, numEdges++);
351: EG_getRange(edge, range, periodic);
352: params[0] = 0.5*(range[0] + range[1]);
353: EG_evaluate(edge, params, result);
354: coords[v*cdim+0] = result[0];
355: coords[v*cdim+1] = result[1];
356: coords[v*cdim+2] = result[2];
357: }
358: }
359: if (Ner == 4) {
360: PetscInt v = fOff + numQuads++;
361: ego *fobjs, face;
362: double range[4], params[3] = {0., 0., 0.}, result[18];
363: int Nf, fid, periodic[2];
365: EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
366: face = fobjs[0];
367: fid = EG_indexBodyTopo(body, face);
368: if (Nf != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1 (%d)", lid-1, Nf, fid);
369: EG_getRange(face, range, periodic);
370: params[0] = 0.5*(range[0] + range[1]);
371: params[1] = 0.5*(range[2] + range[3]);
372: EG_evaluate(face, params, result);
373: coords[v*cdim+0] = result[0];
374: coords[v*cdim+1] = result[1];
375: coords[v*cdim+2] = result[2];
376: }
377: }
378: }
379: if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices);
380: CHKMEMQ;
382: /* Get cell vertices by traversing loops */
383: numQuads = 0;
384: cOff = 0;
385: for (b = 0; b < nbodies; ++b) {
386: ego body = bodies[b];
387: int id, Nl, l;
389: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
390: for (l = 0; l < Nl; ++l) {
391: ego loop = lobjs[l];
392: int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t;
394: lid = EG_indexBodyTopo(body, loop);
395: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
397: for (e = 0; e < Ne; ++e) {
398: ego edge = objs[e];
399: int points[3];
400: int eid, Nv, v, tmp;
402: eid = EG_indexBodyTopo(body, edge);
403: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
404: if (mtype == DEGENERATE) continue;
405: else ++Ner;
406: if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv);
408: for (v = 0; v < Nv; ++v) {
409: ego vertex = nobjs[v];
411: id = EG_indexBodyTopo(body, vertex);
412: points[v*2] = id-1;
413: }
414: {
415: PetscInt edgeNum;
417: PetscHMapIGet(edgeMap, eid-1, &edgeNum);
418: points[1] = numVertices - newVertices + edgeNum;
419: }
420: /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */
421: if (!nc) {
422: for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v];
423: } else {
424: if (cone[nc-1] == points[0]) {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];}
425: else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];}
426: 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];}
427: 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];}
428: else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid);
429: }
430: }
431: if (nc != 2*Ner) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner);
432: if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;}
433: if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners);
434: /* Triangulate the loop */
435: switch (Ner) {
436: case 2: /* Bi-Segment -> 2 triangles */
437: Nt = 2;
438: cells[cOff*numCorners+0] = cone[0];
439: cells[cOff*numCorners+1] = cone[1];
440: cells[cOff*numCorners+2] = cone[2];
441: ++cOff;
442: cells[cOff*numCorners+0] = cone[0];
443: cells[cOff*numCorners+1] = cone[2];
444: cells[cOff*numCorners+2] = cone[3];
445: ++cOff;
446: break;
447: case 3: /* Triangle -> 4 triangles */
448: Nt = 4;
449: cells[cOff*numCorners+0] = cone[0];
450: cells[cOff*numCorners+1] = cone[1];
451: cells[cOff*numCorners+2] = cone[5];
452: ++cOff;
453: cells[cOff*numCorners+0] = cone[1];
454: cells[cOff*numCorners+1] = cone[2];
455: cells[cOff*numCorners+2] = cone[3];
456: ++cOff;
457: cells[cOff*numCorners+0] = cone[5];
458: cells[cOff*numCorners+1] = cone[3];
459: cells[cOff*numCorners+2] = cone[4];
460: ++cOff;
461: cells[cOff*numCorners+0] = cone[1];
462: cells[cOff*numCorners+1] = cone[3];
463: cells[cOff*numCorners+2] = cone[5];
464: ++cOff;
465: break;
466: case 4: /* Quad -> 8 triangles */
467: Nt = 8;
468: cells[cOff*numCorners+0] = cone[0];
469: cells[cOff*numCorners+1] = cone[1];
470: cells[cOff*numCorners+2] = cone[7];
471: ++cOff;
472: cells[cOff*numCorners+0] = cone[1];
473: cells[cOff*numCorners+1] = cone[2];
474: cells[cOff*numCorners+2] = cone[3];
475: ++cOff;
476: cells[cOff*numCorners+0] = cone[3];
477: cells[cOff*numCorners+1] = cone[4];
478: cells[cOff*numCorners+2] = cone[5];
479: ++cOff;
480: cells[cOff*numCorners+0] = cone[5];
481: cells[cOff*numCorners+1] = cone[6];
482: cells[cOff*numCorners+2] = cone[7];
483: ++cOff;
484: cells[cOff*numCorners+0] = cone[8];
485: cells[cOff*numCorners+1] = cone[1];
486: cells[cOff*numCorners+2] = cone[3];
487: ++cOff;
488: cells[cOff*numCorners+0] = cone[8];
489: cells[cOff*numCorners+1] = cone[3];
490: cells[cOff*numCorners+2] = cone[5];
491: ++cOff;
492: cells[cOff*numCorners+0] = cone[8];
493: cells[cOff*numCorners+1] = cone[5];
494: cells[cOff*numCorners+2] = cone[7];
495: ++cOff;
496: cells[cOff*numCorners+0] = cone[8];
497: cells[cOff*numCorners+1] = cone[7];
498: cells[cOff*numCorners+2] = cone[1];
499: ++cOff;
500: break;
501: default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner);
502: }
503: if (debug) {
504: for (t = 0; t < Nt; ++t) {
505: PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t);
506: for (c = 0; c < numCorners; ++c) {
507: if (c > 0) {PetscPrintf(PETSC_COMM_SELF, ", ");}
508: PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]);
509: }
510: PetscPrintf(PETSC_COMM_SELF, ")\n");
511: }
512: }
513: }
514: EG_free(lobjs);
515: }
516: }
517: if (cOff != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells);
518: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);
519: PetscFree3(coords, cells, cone);
520: PetscInfo2(dm, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);
521: PetscInfo2(dm, " Total Number of Unique Vertices = %D (%D)\n", numVertices, newVertices);
522: /* Embed EGADS model in DM */
523: {
524: PetscContainer modelObj, contextObj;
526: PetscContainerCreate(PETSC_COMM_SELF, &modelObj);
527: PetscContainerSetPointer(modelObj, model);
528: PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);
529: PetscContainerDestroy(&modelObj);
531: PetscContainerCreate(PETSC_COMM_SELF, &contextObj);
532: PetscContainerSetPointer(contextObj, context);
533: PetscContainerSetUserDestroy(contextObj, DMPlexEGADSDestroy_Private);
534: PetscObjectCompose((PetscObject) dm, "EGADS Context", (PetscObject) contextObj);
535: PetscContainerDestroy(&contextObj);
536: }
537: /* Label points */
538: DMCreateLabel(dm, "EGADS Body ID");
539: DMGetLabel(dm, "EGADS Body ID", &bodyLabel);
540: DMCreateLabel(dm, "EGADS Face ID");
541: DMGetLabel(dm, "EGADS Face ID", &faceLabel);
542: DMCreateLabel(dm, "EGADS Edge ID");
543: DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);
544: cOff = 0;
545: for (b = 0; b < nbodies; ++b) {
546: ego body = bodies[b];
547: int id, Nl, l;
549: EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);
550: for (l = 0; l < Nl; ++l) {
551: ego loop = lobjs[l];
552: ego *fobjs;
553: int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t;
555: lid = EG_indexBodyTopo(body, loop);
556: EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);
557: if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);
558: fid = EG_indexBodyTopo(body, fobjs[0]);
559: EG_free(fobjs);
560: EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);
561: for (e = 0; e < Ne; ++e) {
562: ego edge = objs[e];
563: int eid, Nv, v;
564: PetscInt points[3], support[2], numEdges, edgeNum;
565: const PetscInt *edges;
567: eid = EG_indexBodyTopo(body, edge);
568: EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);
569: if (mtype == DEGENERATE) continue;
570: else ++Ner;
571: for (v = 0; v < Nv; ++v) {
572: ego vertex = nobjs[v];
574: id = EG_indexBodyTopo(body, vertex);
575: DMLabelSetValue(edgeLabel, numCells + id-1, eid);
576: points[v*2] = numCells + id-1;
577: }
578: PetscHMapIGet(edgeMap, eid-1, &edgeNum);
579: points[1] = numCells + numVertices - newVertices + edgeNum;
581: DMLabelSetValue(edgeLabel, points[1], eid);
582: support[0] = points[0];
583: support[1] = points[1];
584: DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
585: 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);
586: DMLabelSetValue(edgeLabel, edges[0], eid);
587: DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
588: support[0] = points[1];
589: support[1] = points[2];
590: DMPlexGetJoin(dm, 2, support, &numEdges, &edges);
591: 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);
592: DMLabelSetValue(edgeLabel, edges[0], eid);
593: DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);
594: }
595: switch (Ner) {
596: case 2: Nt = 2;break;
597: case 3: Nt = 4;break;
598: case 4: Nt = 8;break;
599: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner);
600: }
601: for (t = 0; t < Nt; ++t) {
602: DMLabelSetValue(bodyLabel, cOff+t, b);
603: DMLabelSetValue(faceLabel, cOff+t, fid);
604: }
605: cOff += Nt;
606: }
607: EG_free(lobjs);
608: }
609: PetscHMapIDestroy(&edgeMap);
610: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
611: for (c = cStart; c < cEnd; ++c) {
612: PetscInt *closure = NULL;
613: PetscInt clSize, cl, bval, fval;
615: DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
616: DMLabelGetValue(bodyLabel, c, &bval);
617: DMLabelGetValue(faceLabel, c, &fval);
618: for (cl = 0; cl < clSize*2; cl += 2) {
619: DMLabelSetValue(bodyLabel, closure[cl], bval);
620: DMLabelSetValue(faceLabel, closure[cl], fval);
621: }
622: DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);
623: }
624: *newdm = dm;
625: return(0);
626: }
627: #endif
629: /*@C
630: DMPlexCreateEGADSFromFile - Create a DMPlex mesh from an EGADSLite file.
632: Collective
634: Input Parameters:
635: + comm - The MPI communicator
636: - filename - The name of the ExodusII file
638: Output Parameter:
639: . dm - The DM object representing the mesh
641: Level: beginner
643: .seealso: DMPLEX, DMCreate(), DMPlexCreateEGADS()
644: @*/
645: PetscErrorCode DMPlexCreateEGADSFromFile(MPI_Comm comm, const char filename[], DM *dm)
646: {
647: PetscMPIInt rank;
648: #if defined(PETSC_HAVE_EGADS)
649: ego context= NULL, model = NULL;
650: #endif
651: PetscBool printModel = PETSC_FALSE;
656: PetscOptionsGetBool(NULL, NULL, "-dm_plex_egads_print_model", &printModel, NULL);
657: MPI_Comm_rank(comm, &rank);
658: #if defined(PETSC_HAVE_EGADS)
659: if (!rank) {
661: EG_open(&context);
662: EG_loadModel(context, 0, filename, &model);
663: if (printModel) {DMPlexEGADSPrintModel(model);}
665: }
666: DMPlexCreateEGADS(comm, context, model, dm);
667: return(0);
668: #else
669: SETERRQ(comm, PETSC_ERR_SUP, "This method requires EGADSLite support. Reconfigure using --download-egads");
670: #endif
671: }