Actual source code: plexgmsh.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmpleximpl.h>
3: #include <petsc/private/hashmapi.h>
5: #include <../src/dm/impls/plex/gmshlex.h>
7: #define GMSH_LEXORDER_ITEM(T, p) \
8: static int *GmshLexOrder_##T##_##p(void) \
9: { \
10: static int Gmsh_LexOrder_##T##_##p[GmshNumNodes_##T(p)] = {-1}; \
11: int *lex = Gmsh_LexOrder_##T##_##p; \
12: if (lex[0] == -1) (void)GmshLexOrder_##T(p, lex, 0); \
13: return lex; \
14: }
16: #define GMSH_LEXORDER_LIST(T) \
17: GMSH_LEXORDER_ITEM(T, 1) \
18: GMSH_LEXORDER_ITEM(T, 2) \
19: GMSH_LEXORDER_ITEM(T, 3) \
20: GMSH_LEXORDER_ITEM(T, 4) \
21: GMSH_LEXORDER_ITEM(T, 5) \
22: GMSH_LEXORDER_ITEM(T, 6) \
23: GMSH_LEXORDER_ITEM(T, 7) \
24: GMSH_LEXORDER_ITEM(T, 8) \
25: GMSH_LEXORDER_ITEM(T, 9) \
26: GMSH_LEXORDER_ITEM(T, 10)
28: GMSH_LEXORDER_ITEM(VTX, 0)
29: GMSH_LEXORDER_LIST(SEG)
30: GMSH_LEXORDER_LIST(TRI)
31: GMSH_LEXORDER_LIST(QUA)
32: GMSH_LEXORDER_LIST(TET)
33: GMSH_LEXORDER_LIST(HEX)
34: GMSH_LEXORDER_LIST(PRI)
35: GMSH_LEXORDER_LIST(PYR)
37: typedef enum {
38: GMSH_VTX = 0,
39: GMSH_SEG = 1,
40: GMSH_TRI = 2,
41: GMSH_QUA = 3,
42: GMSH_TET = 4,
43: GMSH_HEX = 5,
44: GMSH_PRI = 6,
45: GMSH_PYR = 7,
46: GMSH_NUM_POLYTOPES = 8
47: } GmshPolytopeType;
49: typedef struct {
50: int cellType;
51: int polytope;
52: int dim;
53: int order;
54: int numVerts;
55: int numNodes;
56: int* (*lexorder)(void);
57: } GmshCellInfo;
59: #define GmshCellEntry(cellType, polytope, dim, order) \
60: {cellType, GMSH_##polytope, dim, order, \
61: GmshNumNodes_##polytope(1), \
62: GmshNumNodes_##polytope(order), \
63: GmshLexOrder_##polytope##_##order}
65: static const GmshCellInfo GmshCellTable[] = {
66: GmshCellEntry( 15, VTX, 0, 0),
68: GmshCellEntry( 1, SEG, 1, 1),
69: GmshCellEntry( 8, SEG, 1, 2),
70: GmshCellEntry( 26, SEG, 1, 3),
71: GmshCellEntry( 27, SEG, 1, 4),
72: GmshCellEntry( 28, SEG, 1, 5),
73: GmshCellEntry( 62, SEG, 1, 6),
74: GmshCellEntry( 63, SEG, 1, 7),
75: GmshCellEntry( 64, SEG, 1, 8),
76: GmshCellEntry( 65, SEG, 1, 9),
77: GmshCellEntry( 66, SEG, 1, 10),
79: GmshCellEntry( 2, TRI, 2, 1),
80: GmshCellEntry( 9, TRI, 2, 2),
81: GmshCellEntry( 21, TRI, 2, 3),
82: GmshCellEntry( 23, TRI, 2, 4),
83: GmshCellEntry( 25, TRI, 2, 5),
84: GmshCellEntry( 42, TRI, 2, 6),
85: GmshCellEntry( 43, TRI, 2, 7),
86: GmshCellEntry( 44, TRI, 2, 8),
87: GmshCellEntry( 45, TRI, 2, 9),
88: GmshCellEntry( 46, TRI, 2, 10),
90: GmshCellEntry( 3, QUA, 2, 1),
91: GmshCellEntry( 10, QUA, 2, 2),
92: GmshCellEntry( 36, QUA, 2, 3),
93: GmshCellEntry( 37, QUA, 2, 4),
94: GmshCellEntry( 38, QUA, 2, 5),
95: GmshCellEntry( 47, QUA, 2, 6),
96: GmshCellEntry( 48, QUA, 2, 7),
97: GmshCellEntry( 49, QUA, 2, 8),
98: GmshCellEntry( 50, QUA, 2, 9),
99: GmshCellEntry( 51, QUA, 2, 10),
101: GmshCellEntry( 4, TET, 3, 1),
102: GmshCellEntry( 11, TET, 3, 2),
103: GmshCellEntry( 29, TET, 3, 3),
104: GmshCellEntry( 30, TET, 3, 4),
105: GmshCellEntry( 31, TET, 3, 5),
106: GmshCellEntry( 71, TET, 3, 6),
107: GmshCellEntry( 72, TET, 3, 7),
108: GmshCellEntry( 73, TET, 3, 8),
109: GmshCellEntry( 74, TET, 3, 9),
110: GmshCellEntry( 75, TET, 3, 10),
112: GmshCellEntry( 5, HEX, 3, 1),
113: GmshCellEntry( 12, HEX, 3, 2),
114: GmshCellEntry( 92, HEX, 3, 3),
115: GmshCellEntry( 93, HEX, 3, 4),
116: GmshCellEntry( 94, HEX, 3, 5),
117: GmshCellEntry( 95, HEX, 3, 6),
118: GmshCellEntry( 96, HEX, 3, 7),
119: GmshCellEntry( 97, HEX, 3, 8),
120: GmshCellEntry( 98, HEX, 3, 9),
121: GmshCellEntry( -1, HEX, 3, 10),
123: GmshCellEntry( 6, PRI, 3, 1),
124: GmshCellEntry( 13, PRI, 3, 2),
125: GmshCellEntry( 90, PRI, 3, 3),
126: GmshCellEntry( 91, PRI, 3, 4),
127: GmshCellEntry(106, PRI, 3, 5),
128: GmshCellEntry(107, PRI, 3, 6),
129: GmshCellEntry(108, PRI, 3, 7),
130: GmshCellEntry(109, PRI, 3, 8),
131: GmshCellEntry(110, PRI, 3, 9),
132: GmshCellEntry( -1, PRI, 3, 10),
134: GmshCellEntry( 7, PYR, 3, 1),
135: GmshCellEntry( 14, PYR, 3, 2),
136: GmshCellEntry(118, PYR, 3, 3),
137: GmshCellEntry(119, PYR, 3, 4),
138: GmshCellEntry(120, PYR, 3, 5),
139: GmshCellEntry(121, PYR, 3, 6),
140: GmshCellEntry(122, PYR, 3, 7),
141: GmshCellEntry(123, PYR, 3, 8),
142: GmshCellEntry(124, PYR, 3, 9),
143: GmshCellEntry( -1, PYR, 3, 10)
145: #if 0
146: {20, GMSH_TRI, 2, 3, 3, 9, NULL},
147: {16, GMSH_QUA, 2, 2, 4, 8, NULL},
148: {17, GMSH_HEX, 3, 2, 8, 20, NULL},
149: {18, GMSH_PRI, 3, 2, 6, 15, NULL},
150: {19, GMSH_PYR, 3, 2, 5, 13, NULL},
151: #endif
152: };
154: static GmshCellInfo GmshCellMap[150];
156: static PetscErrorCode GmshCellInfoSetUp(void)
157: {
158: size_t i, n;
159: static PetscBool called = PETSC_FALSE;
161: if (called) return 0;
163: called = PETSC_TRUE;
164: n = sizeof(GmshCellMap)/sizeof(GmshCellMap[0]);
165: for (i = 0; i < n; ++i) {
166: GmshCellMap[i].cellType = -1;
167: GmshCellMap[i].polytope = -1;
168: }
169: n = sizeof(GmshCellTable)/sizeof(GmshCellTable[0]);
170: for (i = 0; i < n; ++i) {
171: if (GmshCellTable[i].cellType <= 0) continue;
172: GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i];
173: }
174: return(0);
175: }
177: #define GmshCellTypeCheck(ct) 0; do { \
178: const int _ct_ = (int)ct; \
179: if (_ct_ < 0 || _ct_ >= (int)(sizeof(GmshCellMap)/sizeof(GmshCellMap[0]))) \
180: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Gmsh element type %d", _ct_); \
181: if (GmshCellMap[_ct_].cellType != _ct_) \
182: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \
183: if (GmshCellMap[_ct_].polytope == -1) \
184: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \
185: } while (0)
188: typedef struct {
189: PetscViewer viewer;
190: int fileFormat;
191: int dataSize;
192: PetscBool binary;
193: PetscBool byteSwap;
194: size_t wlen;
195: void *wbuf;
196: size_t slen;
197: void *sbuf;
198: PetscInt *nbuf;
199: PetscInt nodeStart;
200: PetscInt nodeEnd;
201: PetscInt *nodeMap;
202: } GmshFile;
204: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
205: {
206: size_t size = count * eltsize;
210: if (gmsh->wlen < size) {
211: PetscFree(gmsh->wbuf);
212: PetscMalloc(size, &gmsh->wbuf);
213: gmsh->wlen = size;
214: }
215: *(void**)buf = size ? gmsh->wbuf : NULL;
216: return(0);
217: }
219: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
220: {
221: size_t dataSize = (size_t)gmsh->dataSize;
222: size_t size = count * dataSize;
226: if (gmsh->slen < size) {
227: PetscFree(gmsh->sbuf);
228: PetscMalloc(size, &gmsh->sbuf);
229: gmsh->slen = size;
230: }
231: *(void**)buf = size ? gmsh->sbuf : NULL;
232: return(0);
233: }
235: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
236: {
239: PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
240: if (gmsh->byteSwap) {PetscByteSwap(buf, dtype, count);}
241: return(0);
242: }
244: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
245: {
248: PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
249: return(0);
250: }
252: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
253: {
256: PetscStrcmp(line, Section, match);
257: return(0);
258: }
260: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
261: {
262: PetscBool match;
266: GmshMatch(gmsh, Section, line, &match);
267: if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s",Section);
268: return(0);
269: }
271: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
272: {
273: PetscBool match;
277: while (PETSC_TRUE) {
278: GmshReadString(gmsh, line, 1);
279: GmshMatch(gmsh, "$Comments", line, &match);
280: if (!match) break;
281: while (PETSC_TRUE) {
282: GmshReadString(gmsh, line, 1);
283: GmshMatch(gmsh, "$EndComments", line, &match);
284: if (match) break;
285: }
286: }
287: return(0);
288: }
290: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
291: {
294: GmshReadString(gmsh, line, 1);
295: GmshExpect(gmsh, EndSection, line);
296: return(0);
297: }
299: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
300: {
301: PetscInt i;
302: size_t dataSize = (size_t)gmsh->dataSize;
306: if (dataSize == sizeof(PetscInt)) {
307: GmshRead(gmsh, buf, count, PETSC_INT);
308: } else if (dataSize == sizeof(int)) {
309: int *ibuf = NULL;
310: GmshBufferSizeGet(gmsh, count, &ibuf);
311: GmshRead(gmsh, ibuf, count, PETSC_ENUM);
312: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
313: } else if (dataSize == sizeof(long)) {
314: long *ibuf = NULL;
315: GmshBufferSizeGet(gmsh, count, &ibuf);
316: GmshRead(gmsh, ibuf, count, PETSC_LONG);
317: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
318: } else if (dataSize == sizeof(PetscInt64)) {
319: PetscInt64 *ibuf = NULL;
320: GmshBufferSizeGet(gmsh, count, &ibuf);
321: GmshRead(gmsh, ibuf, count, PETSC_INT64);
322: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
323: }
324: return(0);
325: }
327: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
328: {
331: GmshRead(gmsh, buf, count, PETSC_ENUM);
332: return(0);
333: }
335: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
336: {
339: GmshRead(gmsh, buf, count, PETSC_DOUBLE);
340: return(0);
341: }
343: typedef struct {
344: PetscInt id; /* Entity ID */
345: PetscInt dim; /* Dimension */
346: double bbox[6]; /* Bounding box */
347: PetscInt numTags; /* Size of tag array */
348: int tags[4]; /* Tag array */
349: } GmshEntity;
351: typedef struct {
352: GmshEntity *entity[4];
353: PetscHMapI entityMap[4];
354: } GmshEntities;
356: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
357: {
358: PetscInt dim;
362: PetscNew(entities);
363: for (dim = 0; dim < 4; ++dim) {
364: PetscCalloc1(count[dim], &(*entities)->entity[dim]);
365: PetscHMapICreate(&(*entities)->entityMap[dim]);
366: }
367: return(0);
368: }
370: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
371: {
372: PetscInt dim;
376: if (!*entities) return(0);
377: for (dim = 0; dim < 4; ++dim) {
378: PetscFree((*entities)->entity[dim]);
379: PetscHMapIDestroy(&(*entities)->entityMap[dim]);
380: }
381: PetscFree((*entities));
382: return(0);
383: }
385: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity)
386: {
390: PetscHMapISet(entities->entityMap[dim], eid, index);
391: entities->entity[dim][index].dim = dim;
392: entities->entity[dim][index].id = eid;
393: if (entity) *entity = &entities->entity[dim][index];
394: return(0);
395: }
397: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity)
398: {
399: PetscInt index;
403: PetscHMapIGet(entities->entityMap[dim], eid, &index);
404: *entity = &entities->entity[dim][index];
405: return(0);
406: }
408: typedef struct {
409: PetscInt *id; /* Node IDs */
410: double *xyz; /* Coordinates */
411: } GmshNodes;
413: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
414: {
418: PetscNew(nodes);
419: PetscMalloc1(count*1, &(*nodes)->id);
420: PetscMalloc1(count*3, &(*nodes)->xyz);
421: return(0);
422: }
424: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
425: {
428: if (!*nodes) return(0);
429: PetscFree((*nodes)->id);
430: PetscFree((*nodes)->xyz);
431: PetscFree((*nodes));
432: return(0);
433: }
435: typedef struct {
436: PetscInt id; /* Element ID */
437: PetscInt dim; /* Dimension */
438: PetscInt cellType; /* Cell type */
439: PetscInt numVerts; /* Size of vertex array */
440: PetscInt numNodes; /* Size of node array */
441: PetscInt *nodes; /* Vertex/Node array */
442: PetscInt numTags; /* Size of tag array */
443: int tags[4]; /* Tag array */
444: } GmshElement;
446: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
447: {
451: PetscCalloc1(count, elements);
452: return(0);
453: }
455: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
456: {
460: if (!*elements) return(0);
461: PetscFree(*elements);
462: return(0);
463: }
465: typedef struct {
466: PetscInt dim;
467: PetscInt order;
468: GmshEntities *entities;
469: PetscInt numNodes;
470: GmshNodes *nodelist;
471: PetscInt numElems;
472: GmshElement *elements;
473: PetscInt numVerts;
474: PetscInt numCells;
475: PetscInt *periodMap;
476: PetscInt *vertexMap;
477: PetscSegBuffer segbuf;
478: } GmshMesh;
480: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
481: {
485: PetscNew(mesh);
486: PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);
487: return(0);
488: }
490: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
491: {
495: if (!*mesh) return(0);
496: GmshEntitiesDestroy(&(*mesh)->entities);
497: GmshNodesDestroy(&(*mesh)->nodelist);
498: GmshElementsDestroy(&(*mesh)->elements);
499: PetscFree((*mesh)->periodMap);
500: PetscFree((*mesh)->vertexMap);
501: PetscSegBufferDestroy(&(*mesh)->segbuf);
502: PetscFree((*mesh));
503: return(0);
504: }
506: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
507: {
508: PetscViewer viewer = gmsh->viewer;
509: PetscBool byteSwap = gmsh->byteSwap;
510: char line[PETSC_MAX_PATH_LEN];
511: int n, num, nid, snum;
512: GmshNodes *nodes;
516: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
517: snum = sscanf(line, "%d", &num);
518: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
519: GmshNodesCreate(num, &nodes);
520: mesh->numNodes = num;
521: mesh->nodelist = nodes;
522: for (n = 0; n < num; ++n) {
523: double *xyz = nodes->xyz + n*3;
524: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
525: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
526: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
527: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
528: nodes->id[n] = nid;
529: }
530: return(0);
531: }
533: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
534: file contents multiple times to figure out the true number of cells and facets
535: in the given mesh. To make this more efficient we read the file contents only
536: once and store them in memory, while determining the true number of cells. */
537: static PetscErrorCode GmshReadElements_v22(GmshFile* gmsh, GmshMesh *mesh)
538: {
539: PetscViewer viewer = gmsh->viewer;
540: PetscBool binary = gmsh->binary;
541: PetscBool byteSwap = gmsh->byteSwap;
542: char line[PETSC_MAX_PATH_LEN];
543: int i, c, p, num, ibuf[1+4+1000], snum;
544: int cellType, numElem, numVerts, numNodes, numTags;
545: GmshElement *elements;
546: PetscInt *nodeMap = gmsh->nodeMap;
550: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
551: snum = sscanf(line, "%d", &num);
552: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
553: GmshElementsCreate(num, &elements);
554: mesh->numElems = num;
555: mesh->elements = elements;
556: for (c = 0; c < num;) {
557: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
558: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
560: cellType = binary ? ibuf[0] : ibuf[1];
561: numElem = binary ? ibuf[1] : 1;
562: numTags = ibuf[2];
564: GmshCellTypeCheck(cellType);
565: numVerts = GmshCellMap[cellType].numVerts;
566: numNodes = GmshCellMap[cellType].numNodes;
568: for (i = 0; i < numElem; ++i, ++c) {
569: GmshElement *element = elements + c;
570: const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
571: const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
572: PetscViewerRead(viewer, ibuf+off, nint, NULL, PETSC_ENUM);
573: if (byteSwap) {PetscByteSwap(ibuf+off, PETSC_ENUM, nint);}
574: element->id = id[0];
575: element->dim = GmshCellMap[cellType].dim;
576: element->cellType = cellType;
577: element->numVerts = numVerts;
578: element->numNodes = numNodes;
579: element->numTags = PetscMin(numTags, 4);
580: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
581: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
582: for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
583: }
584: }
585: return(0);
586: }
588: /*
589: $Entities
590: numPoints(unsigned long) numCurves(unsigned long)
591: numSurfaces(unsigned long) numVolumes(unsigned long)
592: // points
593: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
594: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
595: numPhysicals(unsigned long) phyisicalTag[...](int)
596: ...
597: // curves
598: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
599: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
600: numPhysicals(unsigned long) physicalTag[...](int)
601: numBREPVert(unsigned long) tagBREPVert[...](int)
602: ...
603: // surfaces
604: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
605: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
606: numPhysicals(unsigned long) physicalTag[...](int)
607: numBREPCurve(unsigned long) tagBREPCurve[...](int)
608: ...
609: // volumes
610: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
611: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
612: numPhysicals(unsigned long) physicalTag[...](int)
613: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
614: ...
615: $EndEntities
616: */
617: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
618: {
619: PetscViewer viewer = gmsh->viewer;
620: PetscBool byteSwap = gmsh->byteSwap;
621: long index, num, lbuf[4];
622: int dim, eid, numTags, *ibuf, t;
623: PetscInt count[4], i;
624: GmshEntity *entity = NULL;
628: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
629: if (byteSwap) {PetscByteSwap(lbuf, PETSC_LONG, 4);}
630: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
631: GmshEntitiesCreate(count, &mesh->entities);
632: for (dim = 0; dim < 4; ++dim) {
633: for (index = 0; index < count[dim]; ++index) {
634: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
635: if (byteSwap) {PetscByteSwap(&eid, PETSC_ENUM, 1);}
636: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
637: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
638: if (byteSwap) {PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);}
639: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
640: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
641: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
642: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
643: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
644: entity->numTags = numTags = (int) PetscMin(num, 4);
645: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
646: if (dim == 0) continue;
647: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
648: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
649: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
650: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
651: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
652: }
653: }
654: return(0);
655: }
657: /*
658: $Nodes
659: numEntityBlocks(unsigned long) numNodes(unsigned long)
660: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
661: tag(int) x(double) y(double) z(double)
662: ...
663: ...
664: $EndNodes
665: */
666: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
667: {
668: PetscViewer viewer = gmsh->viewer;
669: PetscBool byteSwap = gmsh->byteSwap;
670: long block, node, n, numEntityBlocks, numTotalNodes, numNodes;
671: int info[3], nid;
672: GmshNodes *nodes;
676: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
677: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
678: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
679: if (byteSwap) {PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);}
680: GmshNodesCreate(numTotalNodes, &nodes);
681: mesh->numNodes = numTotalNodes;
682: mesh->nodelist = nodes;
683: for (n = 0, block = 0; block < numEntityBlocks; ++block) {
684: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
685: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
686: if (byteSwap) {PetscByteSwap(&numNodes, PETSC_LONG, 1);}
687: if (gmsh->binary) {
688: size_t nbytes = sizeof(int) + 3*sizeof(double);
689: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */
690: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
691: PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
692: for (node = 0; node < numNodes; ++node, ++n) {
693: char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
694: double *xyz = nodes->xyz + n*3;
695: if (!PetscBinaryBigEndian()) {PetscByteSwap(cnid, PETSC_ENUM, 1);}
696: if (!PetscBinaryBigEndian()) {PetscByteSwap(cxyz, PETSC_DOUBLE, 3);}
697: PetscMemcpy(&nid, cnid, sizeof(int));
698: PetscMemcpy(xyz, cxyz, 3*sizeof(double));
699: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
700: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
701: nodes->id[n] = nid;
702: }
703: } else {
704: for (node = 0; node < numNodes; ++node, ++n) {
705: double *xyz = nodes->xyz + n*3;
706: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
707: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
708: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
709: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
710: nodes->id[n] = nid;
711: }
712: }
713: }
714: return(0);
715: }
717: /*
718: $Elements
719: numEntityBlocks(unsigned long) numElements(unsigned long)
720: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
721: tag(int) numVert[...](int)
722: ...
723: ...
724: $EndElements
725: */
726: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
727: {
728: PetscViewer viewer = gmsh->viewer;
729: PetscBool byteSwap = gmsh->byteSwap;
730: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
731: int p, info[3], *ibuf = NULL;
732: int eid, dim, cellType, numVerts, numNodes, numTags;
733: GmshEntity *entity = NULL;
734: GmshElement *elements;
735: PetscInt *nodeMap = gmsh->nodeMap;
739: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
740: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
741: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
742: if (byteSwap) {PetscByteSwap(&numTotalElements, PETSC_LONG, 1);}
743: GmshElementsCreate(numTotalElements, &elements);
744: mesh->numElems = numTotalElements;
745: mesh->elements = elements;
746: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
747: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
748: if (byteSwap) {PetscByteSwap(info, PETSC_ENUM, 3);}
749: eid = info[0]; dim = info[1]; cellType = info[2];
750: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
751: GmshCellTypeCheck(cellType);
752: numVerts = GmshCellMap[cellType].numVerts;
753: numNodes = GmshCellMap[cellType].numNodes;
754: numTags = entity->numTags;
755: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
756: if (byteSwap) {PetscByteSwap(&numElements, PETSC_LONG, 1);}
757: GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
758: PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
759: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);}
760: for (elem = 0; elem < numElements; ++elem, ++c) {
761: GmshElement *element = elements + c;
762: const int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
763: element->id = id[0];
764: element->dim = dim;
765: element->cellType = cellType;
766: element->numVerts = numVerts;
767: element->numNodes = numNodes;
768: element->numTags = numTags;
769: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
770: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
771: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
772: }
773: }
774: return(0);
775: }
777: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
778: {
779: PetscViewer viewer = gmsh->viewer;
780: int fileFormat = gmsh->fileFormat;
781: PetscBool binary = gmsh->binary;
782: PetscBool byteSwap = gmsh->byteSwap;
783: int numPeriodic, snum, i;
784: char line[PETSC_MAX_PATH_LEN];
785: PetscInt *nodeMap = gmsh->nodeMap;
789: if (fileFormat == 22 || !binary) {
790: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
791: snum = sscanf(line, "%d", &numPeriodic);
792: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
793: } else {
794: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
795: if (byteSwap) {PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);}
796: }
797: for (i = 0; i < numPeriodic; i++) {
798: int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
799: long j, nNodes;
800: double affine[16];
802: if (fileFormat == 22 || !binary) {
803: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
804: snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
805: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
806: } else {
807: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
808: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
809: correspondingDim = ibuf[0]; correspondingTag = ibuf[1]; primaryTag = ibuf[2];
810: }
811: (void)correspondingDim; (void)correspondingTag; (void)primaryTag; /* unused */
813: if (fileFormat == 22 || !binary) {
814: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
815: snum = sscanf(line, "%ld", &nNodes);
816: if (snum != 1) { /* discard transformation and try again */
817: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
818: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
819: snum = sscanf(line, "%ld", &nNodes);
820: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
821: }
822: } else {
823: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
824: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
825: if (nNodes == -1) { /* discard transformation and try again */
826: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
827: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
828: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
829: }
830: }
832: for (j = 0; j < nNodes; j++) {
833: if (fileFormat == 22 || !binary) {
834: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
835: snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
836: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
837: } else {
838: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
839: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 2);}
840: correspondingNode = ibuf[0]; primaryNode = ibuf[1];
841: }
842: correspondingNode = (int) nodeMap[correspondingNode];
843: primaryNode = (int) nodeMap[primaryNode];
844: periodicMap[correspondingNode] = primaryNode;
845: }
846: }
847: return(0);
848: }
850: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
851: $Entities
852: numPoints(size_t) numCurves(size_t)
853: numSurfaces(size_t) numVolumes(size_t)
854: pointTag(int) X(double) Y(double) Z(double)
855: numPhysicalTags(size_t) physicalTag(int) ...
856: ...
857: curveTag(int) minX(double) minY(double) minZ(double)
858: maxX(double) maxY(double) maxZ(double)
859: numPhysicalTags(size_t) physicalTag(int) ...
860: numBoundingPoints(size_t) pointTag(int) ...
861: ...
862: surfaceTag(int) minX(double) minY(double) minZ(double)
863: maxX(double) maxY(double) maxZ(double)
864: numPhysicalTags(size_t) physicalTag(int) ...
865: numBoundingCurves(size_t) curveTag(int) ...
866: ...
867: volumeTag(int) minX(double) minY(double) minZ(double)
868: maxX(double) maxY(double) maxZ(double)
869: numPhysicalTags(size_t) physicalTag(int) ...
870: numBoundngSurfaces(size_t) surfaceTag(int) ...
871: ...
872: $EndEntities
873: */
874: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
875: {
876: PetscInt count[4], index, numTags, i;
877: int dim, eid, *tags = NULL;
878: GmshEntity *entity = NULL;
882: GmshReadSize(gmsh, count, 4);
883: GmshEntitiesCreate(count, &mesh->entities);
884: for (dim = 0; dim < 4; ++dim) {
885: for (index = 0; index < count[dim]; ++index) {
886: GmshReadInt(gmsh, &eid, 1);
887: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
888: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
889: GmshReadSize(gmsh, &numTags, 1);
890: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
891: GmshReadInt(gmsh, tags, numTags);
892: entity->numTags = PetscMin(numTags, 4);
893: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
894: if (dim == 0) continue;
895: GmshReadSize(gmsh, &numTags, 1);
896: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
897: GmshReadInt(gmsh, tags, numTags);
898: }
899: }
900: return(0);
901: }
903: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
904: $Nodes
905: numEntityBlocks(size_t) numNodes(size_t)
906: minNodeTag(size_t) maxNodeTag(size_t)
907: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
908: nodeTag(size_t)
909: ...
910: x(double) y(double) z(double)
911: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
912: < v(double; if parametric and entityDim = 2) >
913: ...
914: ...
915: $EndNodes
916: */
917: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
918: {
919: int info[3];
920: PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, block, node;
921: GmshNodes *nodes;
925: GmshReadSize(gmsh, sizes, 4);
926: numEntityBlocks = sizes[0]; numNodes = sizes[1];
927: GmshNodesCreate(numNodes, &nodes);
928: mesh->numNodes = numNodes;
929: mesh->nodelist = nodes;
930: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
931: GmshReadInt(gmsh, info, 3);
932: if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
933: GmshReadSize(gmsh, &numNodesBlock, 1);
934: GmshReadSize(gmsh, nodes->id+node, numNodesBlock);
935: GmshReadDouble(gmsh, nodes->xyz+node*3, numNodesBlock*3);
936: }
937: gmsh->nodeStart = sizes[2];
938: gmsh->nodeEnd = sizes[3]+1;
939: return(0);
940: }
942: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
943: $Elements
944: numEntityBlocks(size_t) numElements(size_t)
945: minElementTag(size_t) maxElementTag(size_t)
946: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
947: elementTag(size_t) nodeTag(size_t) ...
948: ...
949: ...
950: $EndElements
951: */
952: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
953: {
954: int info[3], eid, dim, cellType;
955: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
956: GmshEntity *entity = NULL;
957: GmshElement *elements;
958: PetscInt *nodeMap = gmsh->nodeMap;
962: GmshReadSize(gmsh, sizes, 4);
963: numEntityBlocks = sizes[0]; numElements = sizes[1];
964: GmshElementsCreate(numElements, &elements);
965: mesh->numElems = numElements;
966: mesh->elements = elements;
967: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
968: GmshReadInt(gmsh, info, 3);
969: dim = info[0]; eid = info[1]; cellType = info[2];
970: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
971: GmshCellTypeCheck(cellType);
972: numVerts = GmshCellMap[cellType].numVerts;
973: numNodes = GmshCellMap[cellType].numNodes;
974: numTags = entity->numTags;
975: GmshReadSize(gmsh, &numBlockElements, 1);
976: GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
977: GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
978: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
979: GmshElement *element = elements + c;
980: const PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
981: element->id = id[0];
982: element->dim = dim;
983: element->cellType = cellType;
984: element->numVerts = numVerts;
985: element->numNodes = numNodes;
986: element->numTags = numTags;
987: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
988: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
989: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
990: }
991: }
992: return(0);
993: }
995: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
996: $Periodic
997: numPeriodicLinks(size_t)
998: entityDim(int) entityTag(int) entityTagPrimary(int)
999: numAffine(size_t) value(double) ...
1000: numCorrespondingNodes(size_t)
1001: nodeTag(size_t) nodeTagPrimary(size_t)
1002: ...
1003: ...
1004: $EndPeriodic
1005: */
1006: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
1007: {
1008: int info[3];
1009: double dbuf[16];
1010: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
1011: PetscInt *nodeMap = gmsh->nodeMap;
1015: GmshReadSize(gmsh, &numPeriodicLinks, 1);
1016: for (link = 0; link < numPeriodicLinks; ++link) {
1017: GmshReadInt(gmsh, info, 3);
1018: GmshReadSize(gmsh, &numAffine, 1);
1019: GmshReadDouble(gmsh, dbuf, numAffine);
1020: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
1021: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
1022: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
1023: for (node = 0; node < numCorrespondingNodes; ++node) {
1024: PetscInt correspondingNode = nodeMap[nodeTags[node*2+0]];
1025: PetscInt primaryNode = nodeMap[nodeTags[node*2+1]];
1026: periodicMap[correspondingNode] = primaryNode;
1027: }
1028: }
1029: return(0);
1030: }
1032: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1033: $MeshFormat // same as MSH version 2
1034: version(ASCII double; currently 4.1)
1035: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1036: data-size(ASCII int; sizeof(size_t))
1037: < int with value one; only in binary mode, to detect endianness >
1038: $EndMeshFormat
1039: */
1040: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1041: {
1042: char line[PETSC_MAX_PATH_LEN];
1043: int snum, fileType, fileFormat, dataSize, checkEndian;
1044: float version;
1048: GmshReadString(gmsh, line, 3);
1049: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1050: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1051: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1052: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1053: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1054: if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII");
1055: if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary");
1056: fileFormat = (int)roundf(version*10);
1057: if (fileFormat <= 40 && dataSize != sizeof(double)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1058: if (fileFormat >= 41 && dataSize != sizeof(int) && dataSize != sizeof(PetscInt64)) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1059: gmsh->fileFormat = fileFormat;
1060: gmsh->dataSize = dataSize;
1061: gmsh->byteSwap = PETSC_FALSE;
1062: if (gmsh->binary) {
1063: GmshReadInt(gmsh, &checkEndian, 1);
1064: if (checkEndian != 1) {
1065: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
1066: if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line);
1067: gmsh->byteSwap = PETSC_TRUE;
1068: }
1069: }
1070: return(0);
1071: }
1073: /*
1074: PhysicalNames
1075: numPhysicalNames(ASCII int)
1076: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1077: ...
1078: $EndPhysicalNames
1079: */
1080: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh)
1081: {
1082: char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
1083: int snum, numRegions, region, dim, tag;
1087: GmshReadString(gmsh, line, 1);
1088: snum = sscanf(line, "%d", &numRegions);
1089: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1090: for (region = 0; region < numRegions; ++region) {
1091: GmshReadString(gmsh, line, 2);
1092: snum = sscanf(line, "%d %d", &dim, &tag);
1093: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1094: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
1095: PetscStrchr(line, '"', &p);
1096: if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1097: PetscStrrchr(line, '"', &q);
1098: if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1099: PetscStrncpy(name, p+1, (size_t)(q-p-1));
1100: }
1101: return(0);
1102: }
1104: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1105: {
1109: switch (gmsh->fileFormat) {
1110: case 41: GmshReadEntities_v41(gmsh, mesh); break;
1111: default: GmshReadEntities_v40(gmsh, mesh); break;
1112: }
1113: return(0);
1114: }
1116: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1117: {
1121: switch (gmsh->fileFormat) {
1122: case 41: GmshReadNodes_v41(gmsh, mesh); break;
1123: case 40: GmshReadNodes_v40(gmsh, mesh); break;
1124: default: GmshReadNodes_v22(gmsh, mesh); break;
1125: }
1127: { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1128: if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1129: PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1130: GmshNodes *nodes = mesh->nodelist;
1131: for (n = 0; n < mesh->numNodes; ++n) {
1132: const PetscInt tag = nodes->id[n];
1133: tagMin = PetscMin(tag, tagMin);
1134: tagMax = PetscMax(tag, tagMax);
1135: }
1136: gmsh->nodeStart = tagMin;
1137: gmsh->nodeEnd = tagMax+1;
1138: }
1139: }
1141: { /* Support for sparse node tags */
1142: PetscInt n, t;
1143: GmshNodes *nodes = mesh->nodelist;
1144: PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);
1145: for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1146: gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1147: for (n = 0; n < mesh->numNodes; ++n) {
1148: const PetscInt tag = nodes->id[n];
1149: if (gmsh->nodeMap[tag] >= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %D", tag);
1150: gmsh->nodeMap[tag] = n;
1151: }
1152: }
1153: return(0);
1154: }
1156: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1157: {
1161: switch (gmsh->fileFormat) {
1162: case 41: GmshReadElements_v41(gmsh, mesh); break;
1163: case 40: GmshReadElements_v40(gmsh, mesh); break;
1164: default: GmshReadElements_v22(gmsh, mesh); break;
1165: }
1167: { /* Reorder elements by codimension and polytope type */
1168: PetscInt ne = mesh->numElems;
1169: GmshElement *elements = mesh->elements;
1170: PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0;
1171: PetscInt offset[GMSH_NUM_POLYTOPES+1], e, k;
1173: for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1174: PetscMemzero(offset,sizeof(offset));
1176: keymap[GMSH_TET] = nk++;
1177: keymap[GMSH_HEX] = nk++;
1178: keymap[GMSH_PRI] = nk++;
1179: keymap[GMSH_PYR] = nk++;
1180: keymap[GMSH_TRI] = nk++;
1181: keymap[GMSH_QUA] = nk++;
1182: keymap[GMSH_SEG] = nk++;
1183: keymap[GMSH_VTX] = nk++;
1185: GmshElementsCreate(mesh->numElems, &mesh->elements);
1186: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1187: for (e = 0; e < ne; ++e) offset[1+key(e)]++;
1188: for (k = 1; k < nk; ++k) offset[k] += offset[k-1];
1189: for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1190: #undef key
1191: GmshElementsDestroy(&elements);
1192: }
1194: { /* Mesh dimension and order */
1195: GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1196: mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0;
1197: mesh->order = elem ? GmshCellMap[elem->cellType].order : 0;
1198: }
1200: {
1201: PetscBT vtx;
1202: PetscInt dim = mesh->dim, e, n, v;
1204: PetscBTCreate(mesh->numNodes, &vtx);
1206: /* Compute number of cells and set of vertices */
1207: mesh->numCells = 0;
1208: for (e = 0; e < mesh->numElems; ++e) {
1209: GmshElement *elem = mesh->elements + e;
1210: if (elem->dim == dim && dim > 0) mesh->numCells++;
1211: for (v = 0; v < elem->numVerts; v++) {
1212: PetscBTSet(vtx, elem->nodes[v]);
1213: }
1214: }
1216: /* Compute numbering for vertices */
1217: mesh->numVerts = 0;
1218: PetscMalloc1(mesh->numNodes, &mesh->vertexMap);
1219: for (n = 0; n < mesh->numNodes; ++n)
1220: mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;
1222: PetscBTDestroy(&vtx);
1223: }
1224: return(0);
1225: }
1227: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1228: {
1229: PetscInt n;
1233: PetscMalloc1(mesh->numNodes, &mesh->periodMap);
1234: for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1235: switch (gmsh->fileFormat) {
1236: case 41: GmshReadPeriodic_v41(gmsh, mesh->periodMap); break;
1237: default: GmshReadPeriodic_v40(gmsh, mesh->periodMap); break;
1238: }
1240: /* Find canonical primary nodes */
1241: for (n = 0; n < mesh->numNodes; ++n)
1242: while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]])
1243: mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];
1245: /* Renumber vertices (filter out correspondings) */
1246: mesh->numVerts = 0;
1247: for (n = 0; n < mesh->numNodes; ++n)
1248: if (mesh->vertexMap[n] >= 0) /* is vertex */
1249: if (mesh->periodMap[n] == n) /* is primary */
1250: mesh->vertexMap[n] = mesh->numVerts++;
1251: for (n = 0; n < mesh->numNodes; ++n)
1252: if (mesh->vertexMap[n] >= 0) /* is vertex */
1253: if (mesh->periodMap[n] != n) /* is corresponding */
1254: mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1255: return(0);
1256: }
1258: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1259: #define DM_POLYTOPE_PYRAMID DM_POLYTOPE_UNKNOWN
1260: static const DMPolytopeType DMPolytopeMap[] = {
1261: /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1262: /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1263: /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1264: /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1265: /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1266: /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1267: /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1268: /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,
1269: DM_POLYTOPE_UNKNOWN
1270: };
1272: PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1273: {
1274: return DMPolytopeMap[GmshCellMap[cellType].polytope];
1275: }
1277: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1278: {
1279: DM K;
1280: PetscSpace P;
1281: PetscDualSpace Q;
1282: PetscQuadrature q, fq;
1283: PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1284: PetscBool endpoint = PETSC_TRUE;
1285: char name[32];
1286: PetscErrorCode ierr;
1289: /* Create space */
1290: PetscSpaceCreate(comm, &P);
1291: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1292: PetscSpacePolynomialSetTensor(P, isTensor);
1293: PetscSpaceSetNumComponents(P, Nc);
1294: PetscSpaceSetNumVariables(P, dim);
1295: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1296: if (prefix) {
1297: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1298: PetscSpaceSetFromOptions(P);
1299: PetscObjectSetOptionsPrefix((PetscObject) P, NULL);
1300: PetscSpaceGetDegree(P, &k, NULL);
1301: }
1302: PetscSpaceSetUp(P);
1303: /* Create dual space */
1304: PetscDualSpaceCreate(comm, &Q);
1305: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1306: PetscDualSpaceLagrangeSetTensor(Q, isTensor);
1307: PetscDualSpaceLagrangeSetContinuity(Q, continuity);
1308: PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
1309: PetscDualSpaceSetNumComponents(Q, Nc);
1310: PetscDualSpaceSetOrder(Q, k);
1311: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1312: PetscDualSpaceSetDM(Q, K);
1313: DMDestroy(&K);
1314: if (prefix) {
1315: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1316: PetscDualSpaceSetFromOptions(Q);
1317: PetscObjectSetOptionsPrefix((PetscObject) Q, NULL);
1318: }
1319: PetscDualSpaceSetUp(Q);
1320: /* Create quadrature */
1321: if (isSimplex) {
1322: PetscDTStroudConicalQuadrature(dim, 1, k+1, -1, +1, &q);
1323: PetscDTStroudConicalQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1324: } else {
1325: PetscDTGaussTensorQuadrature(dim, 1, k+1, -1, +1, &q);
1326: PetscDTGaussTensorQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1327: }
1328: /* Create finite element */
1329: PetscFECreate(comm, fem);
1330: PetscFESetType(*fem, PETSCFEBASIC);
1331: PetscFESetNumComponents(*fem, Nc);
1332: PetscFESetBasisSpace(*fem, P);
1333: PetscFESetDualSpace(*fem, Q);
1334: PetscFESetQuadrature(*fem, q);
1335: PetscFESetFaceQuadrature(*fem, fq);
1336: if (prefix) {
1337: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1338: PetscFESetFromOptions(*fem);
1339: PetscObjectSetOptionsPrefix((PetscObject) *fem, NULL);
1340: }
1341: PetscFESetUp(*fem);
1342: /* Cleanup */
1343: PetscSpaceDestroy(&P);
1344: PetscDualSpaceDestroy(&Q);
1345: PetscQuadratureDestroy(&q);
1346: PetscQuadratureDestroy(&fq);
1347: /* Set finite element name */
1348: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1349: PetscFESetName(*fem, name);
1350: return(0);
1351: }
1353: /*@C
1354: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
1356: + comm - The MPI communicator
1357: . filename - Name of the Gmsh file
1358: - interpolate - Create faces and edges in the mesh
1360: Output Parameter:
1361: . dm - The DM object representing the mesh
1363: Level: beginner
1365: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
1366: @*/
1367: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1368: {
1369: PetscViewer viewer;
1370: PetscMPIInt rank;
1371: int fileType;
1372: PetscViewerType vtype;
1373: PetscErrorCode ierr;
1376: MPI_Comm_rank(comm, &rank);
1378: /* Determine Gmsh file type (ASCII or binary) from file header */
1379: if (!rank) {
1380: GmshFile gmsh[1];
1381: char line[PETSC_MAX_PATH_LEN];
1382: int snum;
1383: float version;
1385: PetscArrayzero(gmsh,1);
1386: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1387: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1388: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1389: PetscViewerFileSetName(gmsh->viewer, filename);
1390: /* Read only the first two lines of the Gmsh file */
1391: GmshReadSection(gmsh, line);
1392: GmshExpect(gmsh, "$MeshFormat", line);
1393: GmshReadString(gmsh, line, 2);
1394: snum = sscanf(line, "%f %d", &version, &fileType);
1395: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1396: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1397: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1398: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1399: PetscViewerDestroy(&gmsh->viewer);
1400: }
1401: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1402: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1404: /* Create appropriate viewer and build plex */
1405: PetscViewerCreate(comm, &viewer);
1406: PetscViewerSetType(viewer, vtype);
1407: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1408: PetscViewerFileSetName(viewer, filename);
1409: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1410: PetscViewerDestroy(&viewer);
1411: return(0);
1412: }
1414: /*@
1415: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
1417: Collective
1419: Input Parameters:
1420: + comm - The MPI communicator
1421: . viewer - The Viewer associated with a Gmsh file
1422: - interpolate - Create faces and edges in the mesh
1424: Output Parameter:
1425: . dm - The DM object representing the mesh
1427: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1429: Level: beginner
1431: .seealso: DMPLEX, DMCreate()
1432: @*/
1433: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1434: {
1435: GmshMesh *mesh = NULL;
1436: PetscViewer parentviewer = NULL;
1437: PetscBT periodicVerts = NULL;
1438: PetscBT periodicCells = NULL;
1439: DM cdm;
1440: PetscSection coordSection;
1441: Vec coordinates;
1442: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL;
1443: PetscInt dim = 0, coordDim = -1, order = 0;
1444: PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0;
1445: PetscInt cell, cone[8], e, n, v, d;
1446: PetscBool binary, usemarker = PETSC_FALSE;
1447: PetscBool hybrid = interpolate, periodic = PETSC_TRUE;
1448: PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1449: PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1450: PetscMPIInt rank;
1454: MPI_Comm_rank(comm, &rank);
1455: PetscObjectOptionsBegin((PetscObject)viewer);
1456: PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1457: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);
1458: PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1459: PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);
1460: PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);
1461: PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1462: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);
1463: PetscOptionsTail();
1464: PetscOptionsEnd();
1466: GmshCellInfoSetUp();
1468: DMCreate(comm, dm);
1469: DMSetType(*dm, DMPLEX);
1470: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1472: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1474: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1475: if (binary) {
1476: parentviewer = viewer;
1477: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1478: }
1480: if (!rank) {
1481: GmshFile gmsh[1];
1482: char line[PETSC_MAX_PATH_LEN];
1483: PetscBool match;
1485: PetscArrayzero(gmsh,1);
1486: gmsh->viewer = viewer;
1487: gmsh->binary = binary;
1489: GmshMeshCreate(&mesh);
1491: /* Read mesh format */
1492: GmshReadSection(gmsh, line);
1493: GmshExpect(gmsh, "$MeshFormat", line);
1494: GmshReadMeshFormat(gmsh);
1495: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1497: /* OPTIONAL Read physical names */
1498: GmshReadSection(gmsh, line);
1499: GmshMatch(gmsh, "$PhysicalNames", line, &match);
1500: if (match) {
1501: GmshExpect(gmsh, "$PhysicalNames", line);
1502: GmshReadPhysicalNames(gmsh);
1503: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1504: /* Initial read for entity section */
1505: GmshReadSection(gmsh, line);
1506: }
1508: /* Read entities */
1509: if (gmsh->fileFormat >= 40) {
1510: GmshExpect(gmsh, "$Entities", line);
1511: GmshReadEntities(gmsh, mesh);
1512: GmshReadEndSection(gmsh, "$EndEntities", line);
1513: /* Initial read for nodes section */
1514: GmshReadSection(gmsh, line);
1515: }
1517: /* Read nodes */
1518: GmshExpect(gmsh, "$Nodes", line);
1519: GmshReadNodes(gmsh, mesh);
1520: GmshReadEndSection(gmsh, "$EndNodes", line);
1522: /* Read elements */
1523: GmshReadSection(gmsh, line);
1524: GmshExpect(gmsh, "$Elements", line);
1525: GmshReadElements(gmsh, mesh);
1526: GmshReadEndSection(gmsh, "$EndElements", line);
1528: /* Read periodic section (OPTIONAL) */
1529: if (periodic) {
1530: GmshReadSection(gmsh, line);
1531: GmshMatch(gmsh, "$Periodic", line, &periodic);
1532: }
1533: if (periodic) {
1534: GmshExpect(gmsh, "$Periodic", line);
1535: GmshReadPeriodic(gmsh, mesh);
1536: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1537: }
1539: PetscFree(gmsh->wbuf);
1540: PetscFree(gmsh->sbuf);
1541: PetscFree(gmsh->nbuf);
1543: dim = mesh->dim;
1544: order = mesh->order;
1545: numNodes = mesh->numNodes;
1546: numElems = mesh->numElems;
1547: numVerts = mesh->numVerts;
1548: numCells = mesh->numCells;
1550: {
1551: GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1552: GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1553: int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1554: int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1555: isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1556: isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1557: hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1558: }
1559: }
1561: if (parentviewer) {
1562: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1563: }
1565: {
1566: int buf[6];
1568: buf[0] = (int)dim;
1569: buf[1] = (int)order;
1570: buf[2] = periodic;
1571: buf[3] = isSimplex;
1572: buf[4] = isHybrid;
1573: buf[5] = hasTetra;
1575: MPI_Bcast(buf, 6, MPI_INT, 0, comm);
1577: dim = buf[0];
1578: order = buf[1];
1579: periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1580: isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1581: isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1582: hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1583: }
1585: if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1586: if (highOrder && isHybrid) SETERRQ(comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet");
1588: /* We do not want this label automatically computed, instead we fill it here */
1589: DMCreateLabel(*dm, "celltype");
1591: /* Allocate the cell-vertex mesh */
1592: DMPlexSetChart(*dm, 0, numCells+numVerts);
1593: for (cell = 0; cell < numCells; ++cell) {
1594: GmshElement *elem = mesh->elements + cell;
1595: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1596: if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1597: DMPlexSetConeSize(*dm, cell, elem->numVerts);
1598: DMPlexSetCellType(*dm, cell, ctype);
1599: }
1600: for (v = numCells; v < numCells+numVerts; ++v) {
1601: DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1602: }
1603: DMSetUp(*dm);
1605: /* Add cell-vertex connections */
1606: for (cell = 0; cell < numCells; ++cell) {
1607: GmshElement *elem = mesh->elements + cell;
1608: for (v = 0; v < elem->numVerts; ++v) {
1609: const PetscInt nn = elem->nodes[v];
1610: const PetscInt vv = mesh->vertexMap[nn];
1611: cone[v] = numCells + vv;
1612: }
1613: DMPlexReorderCell(*dm, cell, cone);
1614: DMPlexSetCone(*dm, cell, cone);
1615: }
1617: DMSetDimension(*dm, dim);
1618: DMPlexSymmetrize(*dm);
1619: DMPlexStratify(*dm);
1620: if (interpolate) {
1621: DM idm;
1623: DMPlexInterpolate(*dm, &idm);
1624: DMDestroy(dm);
1625: *dm = idm;
1626: }
1628: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1629: if (!rank && usemarker) {
1630: PetscInt f, fStart, fEnd;
1632: DMCreateLabel(*dm, "marker");
1633: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1634: for (f = fStart; f < fEnd; ++f) {
1635: PetscInt suppSize;
1637: DMPlexGetSupportSize(*dm, f, &suppSize);
1638: if (suppSize == 1) {
1639: PetscInt *cone = NULL, coneSize, p;
1641: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1642: for (p = 0; p < coneSize; p += 2) {
1643: DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1);
1644: }
1645: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1646: }
1647: }
1648: }
1650: if (!rank) {
1651: PetscInt vStart, vEnd;
1653: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1654: for (cell = 0, e = 0; e < numElems; ++e) {
1655: GmshElement *elem = mesh->elements + e;
1657: /* Create cell sets */
1658: if (elem->dim == dim && dim > 0) {
1659: if (elem->numTags > 0) {
1660: DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, elem->tags[0]);
1661: }
1662: cell++;
1663: }
1665: /* Create face sets */
1666: if (interpolate && elem->dim == dim-1) {
1667: PetscInt joinSize;
1668: const PetscInt *join = NULL;
1669: /* Find the relevant facet with vertex joins */
1670: for (v = 0; v < elem->numVerts; ++v) {
1671: const PetscInt nn = elem->nodes[v];
1672: const PetscInt vv = mesh->vertexMap[nn];
1673: cone[v] = vStart + vv;
1674: }
1675: DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1676: if (joinSize != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex facet for Gmsh element %D (Plex cell %D)", elem->id, e);
1677: DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], elem->tags[0]);
1678: DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1679: }
1681: /* Create vertex sets */
1682: if (elem->dim == 0) {
1683: if (elem->numTags > 0) {
1684: const PetscInt nn = elem->nodes[0];
1685: const PetscInt vv = mesh->vertexMap[nn];
1686: DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, elem->tags[0]);
1687: }
1688: }
1689: }
1690: }
1692: { /* Create Cell/Face/Vertex Sets labels at all processes */
1693: enum {n = 4};
1694: PetscBool flag[n];
1696: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1697: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1698: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1699: flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1700: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1701: if (flag[0]) {DMCreateLabel(*dm, "Cell Sets");}
1702: if (flag[1]) {DMCreateLabel(*dm, "Face Sets");}
1703: if (flag[2]) {DMCreateLabel(*dm, "Vertex Sets");}
1704: if (flag[3]) {DMCreateLabel(*dm, "marker");}
1705: }
1707: if (periodic) {
1708: PetscBTCreate(numVerts, &periodicVerts);
1709: for (n = 0; n < numNodes; ++n) {
1710: if (mesh->vertexMap[n] >= 0) {
1711: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1712: PetscInt m = mesh->periodMap[n];
1713: ierr= PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1714: ierr= PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1715: }
1716: }
1717: }
1718: PetscBTCreate(numCells, &periodicCells);
1719: for (cell = 0; cell < numCells; ++cell) {
1720: GmshElement *elem = mesh->elements + cell;
1721: for (v = 0; v < elem->numVerts; ++v) {
1722: PetscInt nn = elem->nodes[v];
1723: PetscInt vv = mesh->vertexMap[nn];
1724: if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1725: PetscBTSet(periodicCells, cell); break;
1726: }
1727: }
1728: }
1729: }
1731: /* Setup coordinate DM */
1732: if (coordDim < 0) coordDim = dim;
1733: DMSetCoordinateDim(*dm, coordDim);
1734: DMGetCoordinateDM(*dm, &cdm);
1735: if (highOrder) {
1736: PetscFE fe;
1737: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1738: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1740: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1742: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1743: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1744: DMSetField(cdm, 0, NULL, (PetscObject) fe);
1745: PetscFEDestroy(&fe);
1746: DMCreateDS(cdm);
1747: }
1749: /* Create coordinates */
1750: if (highOrder) {
1752: PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim;
1753: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1754: PetscSection section;
1755: PetscScalar *cellCoords;
1757: DMSetLocalSection(cdm, NULL);
1758: DMGetLocalSection(cdm, &coordSection);
1759: PetscSectionClone(coordSection, §ion);
1760: DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1762: DMCreateLocalVector(cdm, &coordinates);
1763: PetscMalloc1(maxDof, &cellCoords);
1764: for (cell = 0; cell < numCells; ++cell) {
1765: GmshElement *elem = mesh->elements + cell;
1766: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1767: for (n = 0; n < elem->numNodes; ++n) {
1768: const PetscInt node = elem->nodes[lexorder[n]];
1769: for (d = 0; d < coordDim; ++d)
1770: cellCoords[n*coordDim+d] = (PetscReal) coords[node*3+d];
1771: }
1772: DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1773: }
1774: PetscSectionDestroy(§ion);
1775: PetscFree(cellCoords);
1777: } else {
1779: PetscInt *nodeMap;
1780: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1781: PetscScalar *pointCoords;
1783: DMGetLocalSection(cdm, &coordSection);
1784: PetscSectionSetNumFields(coordSection, 1);
1785: PetscSectionSetFieldComponents(coordSection, 0, coordDim);
1786: if (periodic) { /* we need to localize coordinates on cells */
1787: PetscSectionSetChart(coordSection, 0, numCells+numVerts);
1788: } else {
1789: PetscSectionSetChart(coordSection, numCells, numCells+numVerts);
1790: }
1791: if (periodic) {
1792: for (cell = 0; cell < numCells; ++cell) {
1793: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1794: GmshElement *elem = mesh->elements + cell;
1795: PetscInt dof = elem->numVerts * coordDim;
1796: PetscSectionSetDof(coordSection, cell, dof);
1797: PetscSectionSetFieldDof(coordSection, cell, 0, dof);
1798: }
1799: }
1800: }
1801: for (v = numCells; v < numCells+numVerts; ++v) {
1802: PetscSectionSetDof(coordSection, v, coordDim);
1803: PetscSectionSetFieldDof(coordSection, v, 0, coordDim);
1804: }
1805: PetscSectionSetUp(coordSection);
1807: DMCreateLocalVector(cdm, &coordinates);
1808: VecGetArray(coordinates, &pointCoords);
1809: if (periodic) {
1810: for (cell = 0; cell < numCells; ++cell) {
1811: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1812: GmshElement *elem = mesh->elements + cell;
1813: PetscInt off, node;
1814: for (v = 0; v < elem->numVerts; ++v)
1815: cone[v] = elem->nodes[v];
1816: DMPlexReorderCell(cdm, cell, cone);
1817: PetscSectionGetOffset(coordSection, cell, &off);
1818: for (v = 0; v < elem->numVerts; ++v)
1819: for (node = cone[v], d = 0; d < coordDim; ++d)
1820: pointCoords[off++] = (PetscReal) coords[node*3+d];
1821: }
1822: }
1823: }
1824: PetscMalloc1(numVerts, &nodeMap);
1825: for (n = 0; n < numNodes; n++)
1826: if (mesh->vertexMap[n] >= 0)
1827: nodeMap[mesh->vertexMap[n]] = n;
1828: for (v = 0; v < numVerts; ++v) {
1829: PetscInt off, node = nodeMap[v];
1830: PetscSectionGetOffset(coordSection, numCells + v, &off);
1831: for (d = 0; d < coordDim; ++d)
1832: pointCoords[off+d] = (PetscReal) coords[node*3+d];
1833: }
1834: PetscFree(nodeMap);
1835: VecRestoreArray(coordinates, &pointCoords);
1837: }
1839: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1840: VecSetBlockSize(coordinates, coordDim);
1841: DMSetCoordinatesLocal(*dm, coordinates);
1842: VecDestroy(&coordinates);
1843: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1845: GmshMeshDestroy(&mesh);
1846: PetscBTDestroy(&periodicVerts);
1847: PetscBTDestroy(&periodicCells);
1849: if (highOrder && project) {
1850: PetscFE fe;
1851: const char prefix[] = "dm_plex_gmsh_project_";
1852: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1853: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
1855: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1857: GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1858: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1859: DMProjectCoordinates(*dm, fe);
1860: PetscFEDestroy(&fe);
1861: }
1863: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1864: return(0);
1865: }