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)
187: typedef struct {
188: PetscViewer viewer;
189: int fileFormat;
190: int dataSize;
191: PetscBool binary;
192: PetscBool byteSwap;
193: size_t wlen;
194: void *wbuf;
195: size_t slen;
196: void *sbuf;
197: PetscInt *nbuf;
198: PetscInt nodeStart;
199: PetscInt nodeEnd;
200: PetscInt *nodeMap;
201: } GmshFile;
203: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
204: {
205: size_t size = count * eltsize;
209: if (gmsh->wlen < size) {
210: PetscFree(gmsh->wbuf);
211: PetscMalloc(size, &gmsh->wbuf);
212: gmsh->wlen = size;
213: }
214: *(void**)buf = size ? gmsh->wbuf : NULL;
215: return(0);
216: }
218: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
219: {
220: size_t dataSize = (size_t)gmsh->dataSize;
221: size_t size = count * dataSize;
225: if (gmsh->slen < size) {
226: PetscFree(gmsh->sbuf);
227: PetscMalloc(size, &gmsh->sbuf);
228: gmsh->slen = size;
229: }
230: *(void**)buf = size ? gmsh->sbuf : NULL;
231: return(0);
232: }
234: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
235: {
238: PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
239: if (gmsh->byteSwap) {PetscByteSwap(buf, dtype, count);}
240: return(0);
241: }
243: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
244: {
247: PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
248: return(0);
249: }
251: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
252: {
255: PetscStrcmp(line, Section, match);
256: return(0);
257: }
259: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
260: {
261: PetscBool match;
265: GmshMatch(gmsh, Section, line, &match);
266: if (!match) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s",Section);
267: return(0);
268: }
270: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
271: {
272: PetscBool match;
276: while (PETSC_TRUE) {
277: GmshReadString(gmsh, line, 1);
278: GmshMatch(gmsh, "$Comments", line, &match);
279: if (!match) break;
280: while (PETSC_TRUE) {
281: GmshReadString(gmsh, line, 1);
282: GmshMatch(gmsh, "$EndComments", line, &match);
283: if (match) break;
284: }
285: }
286: return(0);
287: }
289: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
290: {
293: GmshReadString(gmsh, line, 1);
294: GmshExpect(gmsh, EndSection, line);
295: return(0);
296: }
298: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
299: {
300: PetscInt i;
301: size_t dataSize = (size_t)gmsh->dataSize;
305: if (dataSize == sizeof(PetscInt)) {
306: GmshRead(gmsh, buf, count, PETSC_INT);
307: } else if (dataSize == sizeof(int)) {
308: int *ibuf = NULL;
309: GmshBufferSizeGet(gmsh, count, &ibuf);
310: GmshRead(gmsh, ibuf, count, PETSC_ENUM);
311: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
312: } else if (dataSize == sizeof(long)) {
313: long *ibuf = NULL;
314: GmshBufferSizeGet(gmsh, count, &ibuf);
315: GmshRead(gmsh, ibuf, count, PETSC_LONG);
316: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
317: } else if (dataSize == sizeof(PetscInt64)) {
318: PetscInt64 *ibuf = NULL;
319: GmshBufferSizeGet(gmsh, count, &ibuf);
320: GmshRead(gmsh, ibuf, count, PETSC_INT64);
321: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
322: }
323: return(0);
324: }
326: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
327: {
330: GmshRead(gmsh, buf, count, PETSC_ENUM);
331: return(0);
332: }
334: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
335: {
338: GmshRead(gmsh, buf, count, PETSC_DOUBLE);
339: return(0);
340: }
342: typedef struct {
343: PetscInt id; /* Entity ID */
344: PetscInt dim; /* Dimension */
345: double bbox[6]; /* Bounding box */
346: PetscInt numTags; /* Size of tag array */
347: int tags[4]; /* Tag array */
348: } GmshEntity;
350: typedef struct {
351: GmshEntity *entity[4];
352: PetscHMapI entityMap[4];
353: } GmshEntities;
355: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
356: {
357: PetscInt dim;
361: PetscNew(entities);
362: for (dim = 0; dim < 4; ++dim) {
363: PetscCalloc1(count[dim], &(*entities)->entity[dim]);
364: PetscHMapICreate(&(*entities)->entityMap[dim]);
365: }
366: return(0);
367: }
369: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
370: {
371: PetscInt dim;
375: if (!*entities) return(0);
376: for (dim = 0; dim < 4; ++dim) {
377: PetscFree((*entities)->entity[dim]);
378: PetscHMapIDestroy(&(*entities)->entityMap[dim]);
379: }
380: PetscFree((*entities));
381: return(0);
382: }
384: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity)
385: {
389: PetscHMapISet(entities->entityMap[dim], eid, index);
390: entities->entity[dim][index].dim = dim;
391: entities->entity[dim][index].id = eid;
392: if (entity) *entity = &entities->entity[dim][index];
393: return(0);
394: }
396: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity)
397: {
398: PetscInt index;
402: PetscHMapIGet(entities->entityMap[dim], eid, &index);
403: *entity = &entities->entity[dim][index];
404: return(0);
405: }
407: typedef struct {
408: PetscInt *id; /* Node IDs */
409: double *xyz; /* Coordinates */
410: } GmshNodes;
412: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
413: {
417: PetscNew(nodes);
418: PetscMalloc1(count*1, &(*nodes)->id);
419: PetscMalloc1(count*3, &(*nodes)->xyz);
420: return(0);
421: }
423: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
424: {
427: if (!*nodes) return(0);
428: PetscFree((*nodes)->id);
429: PetscFree((*nodes)->xyz);
430: PetscFree((*nodes));
431: return(0);
432: }
434: typedef struct {
435: PetscInt id; /* Element ID */
436: PetscInt dim; /* Dimension */
437: PetscInt cellType; /* Cell type */
438: PetscInt numVerts; /* Size of vertex array */
439: PetscInt numNodes; /* Size of node array */
440: PetscInt *nodes; /* Vertex/Node array */
441: PetscInt numTags; /* Size of tag array */
442: int tags[4]; /* Tag array */
443: } GmshElement;
445: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
446: {
450: PetscCalloc1(count, elements);
451: return(0);
452: }
454: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
455: {
459: if (!*elements) return(0);
460: PetscFree(*elements);
461: return(0);
462: }
464: typedef struct {
465: PetscInt dim;
466: PetscInt order;
467: GmshEntities *entities;
468: PetscInt numNodes;
469: GmshNodes *nodelist;
470: PetscInt numElems;
471: GmshElement *elements;
472: PetscInt numVerts;
473: PetscInt numCells;
474: PetscInt *periodMap;
475: PetscInt *vertexMap;
476: PetscSegBuffer segbuf;
477: } GmshMesh;
479: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
480: {
484: PetscNew(mesh);
485: PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);
486: return(0);
487: }
489: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
490: {
494: if (!*mesh) return(0);
495: GmshEntitiesDestroy(&(*mesh)->entities);
496: GmshNodesDestroy(&(*mesh)->nodelist);
497: GmshElementsDestroy(&(*mesh)->elements);
498: PetscFree((*mesh)->periodMap);
499: PetscFree((*mesh)->vertexMap);
500: PetscSegBufferDestroy(&(*mesh)->segbuf);
501: PetscFree((*mesh));
502: return(0);
503: }
505: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
506: {
507: PetscViewer viewer = gmsh->viewer;
508: PetscBool byteSwap = gmsh->byteSwap;
509: char line[PETSC_MAX_PATH_LEN];
510: int n, num, nid, snum;
511: GmshNodes *nodes;
515: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
516: snum = sscanf(line, "%d", &num);
517: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
518: GmshNodesCreate(num, &nodes);
519: mesh->numNodes = num;
520: mesh->nodelist = nodes;
521: for (n = 0; n < num; ++n) {
522: double *xyz = nodes->xyz + n*3;
523: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
524: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
525: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
526: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
527: nodes->id[n] = nid;
528: }
529: return(0);
530: }
532: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
533: file contents multiple times to figure out the true number of cells and facets
534: in the given mesh. To make this more efficient we read the file contents only
535: once and store them in memory, while determining the true number of cells. */
536: static PetscErrorCode GmshReadElements_v22(GmshFile* gmsh, GmshMesh *mesh)
537: {
538: PetscViewer viewer = gmsh->viewer;
539: PetscBool binary = gmsh->binary;
540: PetscBool byteSwap = gmsh->byteSwap;
541: char line[PETSC_MAX_PATH_LEN];
542: int i, c, p, num, ibuf[1+4+1000], snum;
543: int cellType, numElem, numVerts, numNodes, numTags;
544: GmshElement *elements;
545: PetscInt *nodeMap = gmsh->nodeMap;
549: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
550: snum = sscanf(line, "%d", &num);
551: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
552: GmshElementsCreate(num, &elements);
553: mesh->numElems = num;
554: mesh->elements = elements;
555: for (c = 0; c < num;) {
556: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
557: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
559: cellType = binary ? ibuf[0] : ibuf[1];
560: numElem = binary ? ibuf[1] : 1;
561: numTags = ibuf[2];
563: GmshCellTypeCheck(cellType);
564: numVerts = GmshCellMap[cellType].numVerts;
565: numNodes = GmshCellMap[cellType].numNodes;
567: for (i = 0; i < numElem; ++i, ++c) {
568: GmshElement *element = elements + c;
569: const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
570: const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
571: PetscViewerRead(viewer, ibuf+off, nint, NULL, PETSC_ENUM);
572: if (byteSwap) {PetscByteSwap(ibuf+off, PETSC_ENUM, nint);}
573: element->id = id[0];
574: element->dim = GmshCellMap[cellType].dim;
575: element->cellType = cellType;
576: element->numVerts = numVerts;
577: element->numNodes = numNodes;
578: element->numTags = PetscMin(numTags, 4);
579: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
580: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
581: for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
582: }
583: }
584: return(0);
585: }
587: /*
588: $Entities
589: numPoints(unsigned long) numCurves(unsigned long)
590: numSurfaces(unsigned long) numVolumes(unsigned long)
591: // points
592: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
593: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
594: numPhysicals(unsigned long) phyisicalTag[...](int)
595: ...
596: // curves
597: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
598: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
599: numPhysicals(unsigned long) physicalTag[...](int)
600: numBREPVert(unsigned long) tagBREPVert[...](int)
601: ...
602: // surfaces
603: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
604: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
605: numPhysicals(unsigned long) physicalTag[...](int)
606: numBREPCurve(unsigned long) tagBREPCurve[...](int)
607: ...
608: // volumes
609: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
610: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
611: numPhysicals(unsigned long) physicalTag[...](int)
612: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
613: ...
614: $EndEntities
615: */
616: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
617: {
618: PetscViewer viewer = gmsh->viewer;
619: PetscBool byteSwap = gmsh->byteSwap;
620: long index, num, lbuf[4];
621: int dim, eid, numTags, *ibuf, t;
622: PetscInt count[4], i;
623: GmshEntity *entity = NULL;
627: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
628: if (byteSwap) {PetscByteSwap(lbuf, PETSC_LONG, 4);}
629: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
630: GmshEntitiesCreate(count, &mesh->entities);
631: for (dim = 0; dim < 4; ++dim) {
632: for (index = 0; index < count[dim]; ++index) {
633: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
634: if (byteSwap) {PetscByteSwap(&eid, PETSC_ENUM, 1);}
635: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
636: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
637: if (byteSwap) {PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);}
638: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
639: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
640: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
641: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
642: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
643: entity->numTags = numTags = (int) PetscMin(num, 4);
644: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
645: if (dim == 0) continue;
646: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
647: if (byteSwap) {PetscByteSwap(&num, PETSC_LONG, 1);}
648: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
649: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
650: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, num);}
651: }
652: }
653: return(0);
654: }
656: /*
657: $Nodes
658: numEntityBlocks(unsigned long) numNodes(unsigned long)
659: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
660: tag(int) x(double) y(double) z(double)
661: ...
662: ...
663: $EndNodes
664: */
665: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
666: {
667: PetscViewer viewer = gmsh->viewer;
668: PetscBool byteSwap = gmsh->byteSwap;
669: long block, node, n, numEntityBlocks, numTotalNodes, numNodes;
670: int info[3], nid;
671: GmshNodes *nodes;
675: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
676: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
677: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
678: if (byteSwap) {PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);}
679: GmshNodesCreate(numTotalNodes, &nodes);
680: mesh->numNodes = numTotalNodes;
681: mesh->nodelist = nodes;
682: for (n = 0, block = 0; block < numEntityBlocks; ++block) {
683: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
684: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
685: if (byteSwap) {PetscByteSwap(&numNodes, PETSC_LONG, 1);}
686: if (gmsh->binary) {
687: size_t nbytes = sizeof(int) + 3*sizeof(double);
688: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */
689: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
690: PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
691: for (node = 0; node < numNodes; ++node, ++n) {
692: char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
693: double *xyz = nodes->xyz + n*3;
694: if (!PetscBinaryBigEndian()) {PetscByteSwap(cnid, PETSC_ENUM, 1);}
695: if (!PetscBinaryBigEndian()) {PetscByteSwap(cxyz, PETSC_DOUBLE, 3);}
696: PetscMemcpy(&nid, cnid, sizeof(int));
697: PetscMemcpy(xyz, cxyz, 3*sizeof(double));
698: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
699: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
700: nodes->id[n] = nid;
701: }
702: } else {
703: for (node = 0; node < numNodes; ++node, ++n) {
704: double *xyz = nodes->xyz + n*3;
705: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
706: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
707: if (byteSwap) {PetscByteSwap(&nid, PETSC_ENUM, 1);}
708: if (byteSwap) {PetscByteSwap(xyz, PETSC_DOUBLE, 3);}
709: nodes->id[n] = nid;
710: }
711: }
712: }
713: return(0);
714: }
716: /*
717: $Elements
718: numEntityBlocks(unsigned long) numElements(unsigned long)
719: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
720: tag(int) numVert[...](int)
721: ...
722: ...
723: $EndElements
724: */
725: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
726: {
727: PetscViewer viewer = gmsh->viewer;
728: PetscBool byteSwap = gmsh->byteSwap;
729: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
730: int p, info[3], *ibuf = NULL;
731: int eid, dim, cellType, numVerts, numNodes, numTags;
732: GmshEntity *entity = NULL;
733: GmshElement *elements;
734: PetscInt *nodeMap = gmsh->nodeMap;
738: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
739: if (byteSwap) {PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);}
740: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
741: if (byteSwap) {PetscByteSwap(&numTotalElements, PETSC_LONG, 1);}
742: GmshElementsCreate(numTotalElements, &elements);
743: mesh->numElems = numTotalElements;
744: mesh->elements = elements;
745: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
746: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
747: if (byteSwap) {PetscByteSwap(info, PETSC_ENUM, 3);}
748: eid = info[0]; dim = info[1]; cellType = info[2];
749: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
750: GmshCellTypeCheck(cellType);
751: numVerts = GmshCellMap[cellType].numVerts;
752: numNodes = GmshCellMap[cellType].numNodes;
753: numTags = entity->numTags;
754: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
755: if (byteSwap) {PetscByteSwap(&numElements, PETSC_LONG, 1);}
756: GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
757: PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
758: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);}
759: for (elem = 0; elem < numElements; ++elem, ++c) {
760: GmshElement *element = elements + c;
761: const int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
762: element->id = id[0];
763: element->dim = dim;
764: element->cellType = cellType;
765: element->numVerts = numVerts;
766: element->numNodes = numNodes;
767: element->numTags = numTags;
768: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
769: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
770: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
771: }
772: }
773: return(0);
774: }
776: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
777: {
778: PetscViewer viewer = gmsh->viewer;
779: int fileFormat = gmsh->fileFormat;
780: PetscBool binary = gmsh->binary;
781: PetscBool byteSwap = gmsh->byteSwap;
782: int numPeriodic, snum, i;
783: char line[PETSC_MAX_PATH_LEN];
784: PetscInt *nodeMap = gmsh->nodeMap;
788: if (fileFormat == 22 || !binary) {
789: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
790: snum = sscanf(line, "%d", &numPeriodic);
791: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
792: } else {
793: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
794: if (byteSwap) {PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);}
795: }
796: for (i = 0; i < numPeriodic; i++) {
797: int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
798: long j, nNodes;
799: double affine[16];
801: if (fileFormat == 22 || !binary) {
802: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
803: snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
804: if (snum != 3) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
805: } else {
806: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
807: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 3);}
808: correspondingDim = ibuf[0]; correspondingTag = ibuf[1]; primaryTag = ibuf[2];
809: }
810: (void)correspondingDim; (void)correspondingTag; (void)primaryTag; /* unused */
812: if (fileFormat == 22 || !binary) {
813: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
814: snum = sscanf(line, "%ld", &nNodes);
815: if (snum != 1) { /* discard transformation and try again */
816: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
817: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
818: snum = sscanf(line, "%ld", &nNodes);
819: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
820: }
821: } else {
822: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
823: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
824: if (nNodes == -1) { /* discard transformation and try again */
825: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
826: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
827: if (byteSwap) {PetscByteSwap(&nNodes, PETSC_LONG, 1);}
828: }
829: }
831: for (j = 0; j < nNodes; j++) {
832: if (fileFormat == 22 || !binary) {
833: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
834: snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
835: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
836: } else {
837: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
838: if (byteSwap) {PetscByteSwap(ibuf, PETSC_ENUM, 2);}
839: correspondingNode = ibuf[0]; primaryNode = ibuf[1];
840: }
841: correspondingNode = (int) nodeMap[correspondingNode];
842: primaryNode = (int) nodeMap[primaryNode];
843: periodicMap[correspondingNode] = primaryNode;
844: }
845: }
846: return(0);
847: }
849: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
850: $Entities
851: numPoints(size_t) numCurves(size_t)
852: numSurfaces(size_t) numVolumes(size_t)
853: pointTag(int) X(double) Y(double) Z(double)
854: numPhysicalTags(size_t) physicalTag(int) ...
855: ...
856: curveTag(int) minX(double) minY(double) minZ(double)
857: maxX(double) maxY(double) maxZ(double)
858: numPhysicalTags(size_t) physicalTag(int) ...
859: numBoundingPoints(size_t) pointTag(int) ...
860: ...
861: surfaceTag(int) minX(double) minY(double) minZ(double)
862: maxX(double) maxY(double) maxZ(double)
863: numPhysicalTags(size_t) physicalTag(int) ...
864: numBoundingCurves(size_t) curveTag(int) ...
865: ...
866: volumeTag(int) minX(double) minY(double) minZ(double)
867: maxX(double) maxY(double) maxZ(double)
868: numPhysicalTags(size_t) physicalTag(int) ...
869: numBoundngSurfaces(size_t) surfaceTag(int) ...
870: ...
871: $EndEntities
872: */
873: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
874: {
875: PetscInt count[4], index, numTags, i;
876: int dim, eid, *tags = NULL;
877: GmshEntity *entity = NULL;
881: GmshReadSize(gmsh, count, 4);
882: GmshEntitiesCreate(count, &mesh->entities);
883: for (dim = 0; dim < 4; ++dim) {
884: for (index = 0; index < count[dim]; ++index) {
885: GmshReadInt(gmsh, &eid, 1);
886: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
887: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
888: GmshReadSize(gmsh, &numTags, 1);
889: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
890: GmshReadInt(gmsh, tags, numTags);
891: entity->numTags = PetscMin(numTags, 4);
892: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
893: if (dim == 0) continue;
894: GmshReadSize(gmsh, &numTags, 1);
895: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
896: GmshReadInt(gmsh, tags, numTags);
897: }
898: }
899: return(0);
900: }
902: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
903: $Nodes
904: numEntityBlocks(size_t) numNodes(size_t)
905: minNodeTag(size_t) maxNodeTag(size_t)
906: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
907: nodeTag(size_t)
908: ...
909: x(double) y(double) z(double)
910: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
911: < v(double; if parametric and entityDim = 2) >
912: ...
913: ...
914: $EndNodes
915: */
916: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
917: {
918: int info[3];
919: PetscInt sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0, block, node;
920: GmshNodes *nodes;
924: GmshReadSize(gmsh, sizes, 4);
925: numEntityBlocks = sizes[0]; numNodes = sizes[1];
926: GmshNodesCreate(numNodes, &nodes);
927: mesh->numNodes = numNodes;
928: mesh->nodelist = nodes;
929: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
930: GmshReadInt(gmsh, info, 3);
931: if (info[2] != 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
932: GmshReadSize(gmsh, &numNodesBlock, 1);
933: GmshReadSize(gmsh, nodes->id+node, numNodesBlock);
934: GmshReadDouble(gmsh, nodes->xyz+node*3, numNodesBlock*3);
935: }
936: gmsh->nodeStart = sizes[2];
937: gmsh->nodeEnd = sizes[3]+1;
938: return(0);
939: }
941: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
942: $Elements
943: numEntityBlocks(size_t) numElements(size_t)
944: minElementTag(size_t) maxElementTag(size_t)
945: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
946: elementTag(size_t) nodeTag(size_t) ...
947: ...
948: ...
949: $EndElements
950: */
951: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
952: {
953: int info[3], eid, dim, cellType;
954: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
955: GmshEntity *entity = NULL;
956: GmshElement *elements;
957: PetscInt *nodeMap = gmsh->nodeMap;
961: GmshReadSize(gmsh, sizes, 4);
962: numEntityBlocks = sizes[0]; numElements = sizes[1];
963: GmshElementsCreate(numElements, &elements);
964: mesh->numElems = numElements;
965: mesh->elements = elements;
966: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
967: GmshReadInt(gmsh, info, 3);
968: dim = info[0]; eid = info[1]; cellType = info[2];
969: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
970: GmshCellTypeCheck(cellType);
971: numVerts = GmshCellMap[cellType].numVerts;
972: numNodes = GmshCellMap[cellType].numNodes;
973: numTags = entity->numTags;
974: GmshReadSize(gmsh, &numBlockElements, 1);
975: GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
976: GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
977: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
978: GmshElement *element = elements + c;
979: const PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
980: element->id = id[0];
981: element->dim = dim;
982: element->cellType = cellType;
983: element->numVerts = numVerts;
984: element->numNodes = numNodes;
985: element->numTags = numTags;
986: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
987: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
988: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
989: }
990: }
991: return(0);
992: }
994: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
995: $Periodic
996: numPeriodicLinks(size_t)
997: entityDim(int) entityTag(int) entityTagPrimary(int)
998: numAffine(size_t) value(double) ...
999: numCorrespondingNodes(size_t)
1000: nodeTag(size_t) nodeTagPrimary(size_t)
1001: ...
1002: ...
1003: $EndPeriodic
1004: */
1005: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
1006: {
1007: int info[3];
1008: double dbuf[16];
1009: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
1010: PetscInt *nodeMap = gmsh->nodeMap;
1014: GmshReadSize(gmsh, &numPeriodicLinks, 1);
1015: for (link = 0; link < numPeriodicLinks; ++link) {
1016: GmshReadInt(gmsh, info, 3);
1017: GmshReadSize(gmsh, &numAffine, 1);
1018: GmshReadDouble(gmsh, dbuf, numAffine);
1019: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
1020: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
1021: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
1022: for (node = 0; node < numCorrespondingNodes; ++node) {
1023: PetscInt correspondingNode = nodeMap[nodeTags[node*2+0]];
1024: PetscInt primaryNode = nodeMap[nodeTags[node*2+1]];
1025: periodicMap[correspondingNode] = primaryNode;
1026: }
1027: }
1028: return(0);
1029: }
1031: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1032: $MeshFormat // same as MSH version 2
1033: version(ASCII double; currently 4.1)
1034: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1035: data-size(ASCII int; sizeof(size_t))
1036: < int with value one; only in binary mode, to detect endianness >
1037: $EndMeshFormat
1038: */
1039: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1040: {
1041: char line[PETSC_MAX_PATH_LEN];
1042: int snum, fileType, fileFormat, dataSize, checkEndian;
1043: float version;
1047: GmshReadString(gmsh, line, 3);
1048: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1049: if (snum != 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1050: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1051: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1052: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1053: if (gmsh->binary && !fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII");
1054: if (!gmsh->binary && fileType) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary");
1055: fileFormat = (int)roundf(version*10);
1056: 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);
1057: 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);
1058: gmsh->fileFormat = fileFormat;
1059: gmsh->dataSize = dataSize;
1060: gmsh->byteSwap = PETSC_FALSE;
1061: if (gmsh->binary) {
1062: GmshReadInt(gmsh, &checkEndian, 1);
1063: if (checkEndian != 1) {
1064: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
1065: if (checkEndian != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line);
1066: gmsh->byteSwap = PETSC_TRUE;
1067: }
1068: }
1069: return(0);
1070: }
1072: /*
1073: PhysicalNames
1074: numPhysicalNames(ASCII int)
1075: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1076: ...
1077: $EndPhysicalNames
1078: */
1079: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh)
1080: {
1081: char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
1082: int snum, numRegions, region, dim, tag;
1086: GmshReadString(gmsh, line, 1);
1087: snum = sscanf(line, "%d", &numRegions);
1088: if (snum != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1089: for (region = 0; region < numRegions; ++region) {
1090: GmshReadString(gmsh, line, 2);
1091: snum = sscanf(line, "%d %d", &dim, &tag);
1092: if (snum != 2) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1093: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
1094: PetscStrchr(line, '"', &p);
1095: if (!p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1096: PetscStrrchr(line, '"', &q);
1097: if (q == p) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1098: PetscStrncpy(name, p+1, (size_t)(q-p-1));
1099: }
1100: return(0);
1101: }
1103: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1104: {
1108: switch (gmsh->fileFormat) {
1109: case 41: GmshReadEntities_v41(gmsh, mesh); break;
1110: default: GmshReadEntities_v40(gmsh, mesh); break;
1111: }
1112: return(0);
1113: }
1115: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1116: {
1120: switch (gmsh->fileFormat) {
1121: case 41: GmshReadNodes_v41(gmsh, mesh); break;
1122: case 40: GmshReadNodes_v40(gmsh, mesh); break;
1123: default: GmshReadNodes_v22(gmsh, mesh); break;
1124: }
1126: { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1127: if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1128: PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1129: GmshNodes *nodes = mesh->nodelist;
1130: for (n = 0; n < mesh->numNodes; ++n) {
1131: const PetscInt tag = nodes->id[n];
1132: tagMin = PetscMin(tag, tagMin);
1133: tagMax = PetscMax(tag, tagMax);
1134: }
1135: gmsh->nodeStart = tagMin;
1136: gmsh->nodeEnd = tagMax+1;
1137: }
1138: }
1140: { /* Support for sparse node tags */
1141: PetscInt n, t;
1142: GmshNodes *nodes = mesh->nodelist;
1143: PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);
1144: for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1145: gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1146: for (n = 0; n < mesh->numNodes; ++n) {
1147: const PetscInt tag = nodes->id[n];
1148: if (gmsh->nodeMap[tag] >= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %D", tag);
1149: gmsh->nodeMap[tag] = n;
1150: }
1151: }
1152: return(0);
1153: }
1155: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1156: {
1160: switch (gmsh->fileFormat) {
1161: case 41: GmshReadElements_v41(gmsh, mesh); break;
1162: case 40: GmshReadElements_v40(gmsh, mesh); break;
1163: default: GmshReadElements_v22(gmsh, mesh); break;
1164: }
1166: { /* Reorder elements by codimension and polytope type */
1167: PetscInt ne = mesh->numElems;
1168: GmshElement *elements = mesh->elements;
1169: PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0;
1170: PetscInt offset[GMSH_NUM_POLYTOPES+1], e, k;
1172: for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1173: PetscMemzero(offset,sizeof(offset));
1175: keymap[GMSH_TET] = nk++;
1176: keymap[GMSH_HEX] = nk++;
1177: keymap[GMSH_PRI] = nk++;
1178: keymap[GMSH_PYR] = nk++;
1179: keymap[GMSH_TRI] = nk++;
1180: keymap[GMSH_QUA] = nk++;
1181: keymap[GMSH_SEG] = nk++;
1182: keymap[GMSH_VTX] = nk++;
1184: GmshElementsCreate(mesh->numElems, &mesh->elements);
1185: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1186: for (e = 0; e < ne; ++e) offset[1+key(e)]++;
1187: for (k = 1; k < nk; ++k) offset[k] += offset[k-1];
1188: for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1189: #undef key
1190: GmshElementsDestroy(&elements);
1191: }
1193: { /* Mesh dimension and order */
1194: GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1195: mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0;
1196: mesh->order = elem ? GmshCellMap[elem->cellType].order : 0;
1197: }
1199: {
1200: PetscBT vtx;
1201: PetscInt dim = mesh->dim, e, n, v;
1203: PetscBTCreate(mesh->numNodes, &vtx);
1205: /* Compute number of cells and set of vertices */
1206: mesh->numCells = 0;
1207: for (e = 0; e < mesh->numElems; ++e) {
1208: GmshElement *elem = mesh->elements + e;
1209: if (elem->dim == dim && dim > 0) mesh->numCells++;
1210: for (v = 0; v < elem->numVerts; v++) {
1211: PetscBTSet(vtx, elem->nodes[v]);
1212: }
1213: }
1215: /* Compute numbering for vertices */
1216: mesh->numVerts = 0;
1217: PetscMalloc1(mesh->numNodes, &mesh->vertexMap);
1218: for (n = 0; n < mesh->numNodes; ++n)
1219: mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;
1221: PetscBTDestroy(&vtx);
1222: }
1223: return(0);
1224: }
1226: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1227: {
1228: PetscInt n;
1232: PetscMalloc1(mesh->numNodes, &mesh->periodMap);
1233: for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1234: switch (gmsh->fileFormat) {
1235: case 41: GmshReadPeriodic_v41(gmsh, mesh->periodMap); break;
1236: default: GmshReadPeriodic_v40(gmsh, mesh->periodMap); break;
1237: }
1239: /* Find canonical primary nodes */
1240: for (n = 0; n < mesh->numNodes; ++n)
1241: while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]])
1242: mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];
1244: /* Renumber vertices (filter out correspondings) */
1245: mesh->numVerts = 0;
1246: for (n = 0; n < mesh->numNodes; ++n)
1247: if (mesh->vertexMap[n] >= 0) /* is vertex */
1248: if (mesh->periodMap[n] == n) /* is primary */
1249: mesh->vertexMap[n] = mesh->numVerts++;
1250: for (n = 0; n < mesh->numNodes; ++n)
1251: if (mesh->vertexMap[n] >= 0) /* is vertex */
1252: if (mesh->periodMap[n] != n) /* is corresponding */
1253: mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1254: return(0);
1255: }
1257: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1258: #define DM_POLYTOPE_PYRAMID DM_POLYTOPE_UNKNOWN
1259: static const DMPolytopeType DMPolytopeMap[] = {
1260: /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1261: /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1262: /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1263: /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1264: /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1265: /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1266: /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1267: /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,
1268: DM_POLYTOPE_UNKNOWN
1269: };
1271: PETSC_STATIC_INLINE DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1272: {
1273: return DMPolytopeMap[GmshCellMap[cellType].polytope];
1274: }
1276: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1277: {
1278: DM K;
1279: PetscSpace P;
1280: PetscDualSpace Q;
1281: PetscQuadrature q, fq;
1282: PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1283: PetscBool endpoint = PETSC_TRUE;
1284: char name[32];
1285: PetscErrorCode ierr;
1288: /* Create space */
1289: PetscSpaceCreate(comm, &P);
1290: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1291: PetscSpacePolynomialSetTensor(P, isTensor);
1292: PetscSpaceSetNumComponents(P, Nc);
1293: PetscSpaceSetNumVariables(P, dim);
1294: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1295: if (prefix) {
1296: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1297: PetscSpaceSetFromOptions(P);
1298: PetscObjectSetOptionsPrefix((PetscObject) P, NULL);
1299: PetscSpaceGetDegree(P, &k, NULL);
1300: }
1301: PetscSpaceSetUp(P);
1302: /* Create dual space */
1303: PetscDualSpaceCreate(comm, &Q);
1304: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1305: PetscDualSpaceLagrangeSetTensor(Q, isTensor);
1306: PetscDualSpaceLagrangeSetContinuity(Q, continuity);
1307: PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
1308: PetscDualSpaceSetNumComponents(Q, Nc);
1309: PetscDualSpaceSetOrder(Q, k);
1310: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
1311: PetscDualSpaceSetDM(Q, K);
1312: DMDestroy(&K);
1313: if (prefix) {
1314: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1315: PetscDualSpaceSetFromOptions(Q);
1316: PetscObjectSetOptionsPrefix((PetscObject) Q, NULL);
1317: }
1318: PetscDualSpaceSetUp(Q);
1319: /* Create quadrature */
1320: if (isSimplex) {
1321: PetscDTStroudConicalQuadrature(dim, 1, k+1, -1, +1, &q);
1322: PetscDTStroudConicalQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1323: } else {
1324: PetscDTGaussTensorQuadrature(dim, 1, k+1, -1, +1, &q);
1325: PetscDTGaussTensorQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1326: }
1327: /* Create finite element */
1328: PetscFECreate(comm, fem);
1329: PetscFESetType(*fem, PETSCFEBASIC);
1330: PetscFESetNumComponents(*fem, Nc);
1331: PetscFESetBasisSpace(*fem, P);
1332: PetscFESetDualSpace(*fem, Q);
1333: PetscFESetQuadrature(*fem, q);
1334: PetscFESetFaceQuadrature(*fem, fq);
1335: if (prefix) {
1336: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1337: PetscFESetFromOptions(*fem);
1338: PetscObjectSetOptionsPrefix((PetscObject) *fem, NULL);
1339: }
1340: PetscFESetUp(*fem);
1341: /* Cleanup */
1342: PetscSpaceDestroy(&P);
1343: PetscDualSpaceDestroy(&Q);
1344: PetscQuadratureDestroy(&q);
1345: PetscQuadratureDestroy(&fq);
1346: /* Set finite element name */
1347: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1348: PetscFESetName(*fem, name);
1349: return(0);
1350: }
1352: /*@C
1353: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
1355: + comm - The MPI communicator
1356: . filename - Name of the Gmsh file
1357: - interpolate - Create faces and edges in the mesh
1359: Output Parameter:
1360: . dm - The DM object representing the mesh
1362: Level: beginner
1364: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
1365: @*/
1366: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1367: {
1368: PetscViewer viewer;
1369: PetscMPIInt rank;
1370: int fileType;
1371: PetscViewerType vtype;
1372: PetscErrorCode ierr;
1375: MPI_Comm_rank(comm, &rank);
1377: /* Determine Gmsh file type (ASCII or binary) from file header */
1378: if (rank == 0) {
1379: GmshFile gmsh[1];
1380: char line[PETSC_MAX_PATH_LEN];
1381: int snum;
1382: float version;
1384: PetscArrayzero(gmsh,1);
1385: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1386: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1387: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1388: PetscViewerFileSetName(gmsh->viewer, filename);
1389: /* Read only the first two lines of the Gmsh file */
1390: GmshReadSection(gmsh, line);
1391: GmshExpect(gmsh, "$MeshFormat", line);
1392: GmshReadString(gmsh, line, 2);
1393: snum = sscanf(line, "%f %d", &version, &fileType);
1394: if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1395: if (version < 2.2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1396: if ((int)version == 3) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1397: if (version > 4.1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1398: PetscViewerDestroy(&gmsh->viewer);
1399: }
1400: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1401: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1403: /* Create appropriate viewer and build plex */
1404: PetscViewerCreate(comm, &viewer);
1405: PetscViewerSetType(viewer, vtype);
1406: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1407: PetscViewerFileSetName(viewer, filename);
1408: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1409: PetscViewerDestroy(&viewer);
1410: return(0);
1411: }
1413: /*@
1414: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
1416: Collective
1418: Input Parameters:
1419: + comm - The MPI communicator
1420: . viewer - The Viewer associated with a Gmsh file
1421: - interpolate - Create faces and edges in the mesh
1423: Output Parameter:
1424: . dm - The DM object representing the mesh
1426: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1428: Level: beginner
1430: .seealso: DMPLEX, DMCreate()
1431: @*/
1432: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1433: {
1434: GmshMesh *mesh = NULL;
1435: PetscViewer parentviewer = NULL;
1436: PetscBT periodicVerts = NULL;
1437: PetscBT periodicCells = NULL;
1438: DM cdm;
1439: PetscSection coordSection;
1440: Vec coordinates;
1441: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL;
1442: PetscInt dim = 0, coordDim = -1, order = 0;
1443: PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0;
1444: PetscInt cell, cone[8], e, n, v, d;
1445: PetscBool binary, usemarker = PETSC_FALSE;
1446: PetscBool hybrid = interpolate, periodic = PETSC_TRUE;
1447: PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1448: PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1449: PetscMPIInt rank;
1453: MPI_Comm_rank(comm, &rank);
1454: PetscObjectOptionsBegin((PetscObject)viewer);
1455: PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1456: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);
1457: PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1458: PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);
1459: PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);
1460: PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1461: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);
1462: PetscOptionsTail();
1463: PetscOptionsEnd();
1465: GmshCellInfoSetUp();
1467: DMCreate(comm, dm);
1468: DMSetType(*dm, DMPLEX);
1469: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1471: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1473: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1474: if (binary) {
1475: parentviewer = viewer;
1476: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1477: }
1479: if (rank == 0) {
1480: GmshFile gmsh[1];
1481: char line[PETSC_MAX_PATH_LEN];
1482: PetscBool match;
1484: PetscArrayzero(gmsh,1);
1485: gmsh->viewer = viewer;
1486: gmsh->binary = binary;
1488: GmshMeshCreate(&mesh);
1490: /* Read mesh format */
1491: GmshReadSection(gmsh, line);
1492: GmshExpect(gmsh, "$MeshFormat", line);
1493: GmshReadMeshFormat(gmsh);
1494: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1496: /* OPTIONAL Read physical names */
1497: GmshReadSection(gmsh, line);
1498: GmshMatch(gmsh, "$PhysicalNames", line, &match);
1499: if (match) {
1500: GmshExpect(gmsh, "$PhysicalNames", line);
1501: GmshReadPhysicalNames(gmsh);
1502: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1503: /* Initial read for entity section */
1504: GmshReadSection(gmsh, line);
1505: }
1507: /* Read entities */
1508: if (gmsh->fileFormat >= 40) {
1509: GmshExpect(gmsh, "$Entities", line);
1510: GmshReadEntities(gmsh, mesh);
1511: GmshReadEndSection(gmsh, "$EndEntities", line);
1512: /* Initial read for nodes section */
1513: GmshReadSection(gmsh, line);
1514: }
1516: /* Read nodes */
1517: GmshExpect(gmsh, "$Nodes", line);
1518: GmshReadNodes(gmsh, mesh);
1519: GmshReadEndSection(gmsh, "$EndNodes", line);
1521: /* Read elements */
1522: GmshReadSection(gmsh, line);
1523: GmshExpect(gmsh, "$Elements", line);
1524: GmshReadElements(gmsh, mesh);
1525: GmshReadEndSection(gmsh, "$EndElements", line);
1527: /* Read periodic section (OPTIONAL) */
1528: if (periodic) {
1529: GmshReadSection(gmsh, line);
1530: GmshMatch(gmsh, "$Periodic", line, &periodic);
1531: }
1532: if (periodic) {
1533: GmshExpect(gmsh, "$Periodic", line);
1534: GmshReadPeriodic(gmsh, mesh);
1535: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1536: }
1538: PetscFree(gmsh->wbuf);
1539: PetscFree(gmsh->sbuf);
1540: PetscFree(gmsh->nbuf);
1542: dim = mesh->dim;
1543: order = mesh->order;
1544: numNodes = mesh->numNodes;
1545: numElems = mesh->numElems;
1546: numVerts = mesh->numVerts;
1547: numCells = mesh->numCells;
1549: {
1550: GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1551: GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1552: int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1553: int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1554: isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1555: isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1556: hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1557: }
1558: }
1560: if (parentviewer) {
1561: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1562: }
1564: {
1565: int buf[6];
1567: buf[0] = (int)dim;
1568: buf[1] = (int)order;
1569: buf[2] = periodic;
1570: buf[3] = isSimplex;
1571: buf[4] = isHybrid;
1572: buf[5] = hasTetra;
1574: MPI_Bcast(buf, 6, MPI_INT, 0, comm);
1576: dim = buf[0];
1577: order = buf[1];
1578: periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1579: isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1580: isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1581: hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1582: }
1584: if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1585: if (highOrder && isHybrid) SETERRQ(comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet");
1587: /* We do not want this label automatically computed, instead we fill it here */
1588: DMCreateLabel(*dm, "celltype");
1590: /* Allocate the cell-vertex mesh */
1591: DMPlexSetChart(*dm, 0, numCells+numVerts);
1592: for (cell = 0; cell < numCells; ++cell) {
1593: GmshElement *elem = mesh->elements + cell;
1594: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1595: if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1596: DMPlexSetConeSize(*dm, cell, elem->numVerts);
1597: DMPlexSetCellType(*dm, cell, ctype);
1598: }
1599: for (v = numCells; v < numCells+numVerts; ++v) {
1600: DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1601: }
1602: DMSetUp(*dm);
1604: /* Add cell-vertex connections */
1605: for (cell = 0; cell < numCells; ++cell) {
1606: GmshElement *elem = mesh->elements + cell;
1607: for (v = 0; v < elem->numVerts; ++v) {
1608: const PetscInt nn = elem->nodes[v];
1609: const PetscInt vv = mesh->vertexMap[nn];
1610: cone[v] = numCells + vv;
1611: }
1612: DMPlexReorderCell(*dm, cell, cone);
1613: DMPlexSetCone(*dm, cell, cone);
1614: }
1616: DMSetDimension(*dm, dim);
1617: DMPlexSymmetrize(*dm);
1618: DMPlexStratify(*dm);
1619: if (interpolate) {
1620: DM idm;
1622: DMPlexInterpolate(*dm, &idm);
1623: DMDestroy(dm);
1624: *dm = idm;
1625: }
1627: if (usemarker && !interpolate && dim > 1) SETERRQ(comm,PETSC_ERR_SUP,"Cannot create marker label without interpolation");
1628: if (rank == 0 && usemarker) {
1629: PetscInt f, fStart, fEnd;
1631: DMCreateLabel(*dm, "marker");
1632: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1633: for (f = fStart; f < fEnd; ++f) {
1634: PetscInt suppSize;
1636: DMPlexGetSupportSize(*dm, f, &suppSize);
1637: if (suppSize == 1) {
1638: PetscInt *cone = NULL, coneSize, p;
1640: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1641: for (p = 0; p < coneSize; p += 2) {
1642: DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1);
1643: }
1644: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1645: }
1646: }
1647: }
1649: if (rank == 0) {
1650: PetscInt vStart, vEnd;
1652: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1653: for (cell = 0, e = 0; e < numElems; ++e) {
1654: GmshElement *elem = mesh->elements + e;
1656: /* Create cell sets */
1657: if (elem->dim == dim && dim > 0) {
1658: if (elem->numTags > 0) {
1659: DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, elem->tags[0]);
1660: }
1661: cell++;
1662: }
1664: /* Create face sets */
1665: if (interpolate && elem->dim == dim-1) {
1666: PetscInt joinSize;
1667: const PetscInt *join = NULL;
1668: /* Find the relevant facet with vertex joins */
1669: for (v = 0; v < elem->numVerts; ++v) {
1670: const PetscInt nn = elem->nodes[v];
1671: const PetscInt vv = mesh->vertexMap[nn];
1672: cone[v] = vStart + vv;
1673: }
1674: DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1675: 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);
1676: DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], elem->tags[0]);
1677: DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1678: }
1680: /* Create vertex sets */
1681: if (elem->dim == 0) {
1682: if (elem->numTags > 0) {
1683: const PetscInt nn = elem->nodes[0];
1684: const PetscInt vv = mesh->vertexMap[nn];
1685: DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, elem->tags[0]);
1686: }
1687: }
1688: }
1689: }
1691: { /* Create Cell/Face/Vertex Sets labels at all processes */
1692: enum {n = 4};
1693: PetscBool flag[n];
1695: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1696: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1697: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1698: flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1699: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1700: if (flag[0]) {DMCreateLabel(*dm, "Cell Sets");}
1701: if (flag[1]) {DMCreateLabel(*dm, "Face Sets");}
1702: if (flag[2]) {DMCreateLabel(*dm, "Vertex Sets");}
1703: if (flag[3]) {DMCreateLabel(*dm, "marker");}
1704: }
1706: if (periodic) {
1707: PetscBTCreate(numVerts, &periodicVerts);
1708: for (n = 0; n < numNodes; ++n) {
1709: if (mesh->vertexMap[n] >= 0) {
1710: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1711: PetscInt m = mesh->periodMap[n];
1712: PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1713: PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1714: }
1715: }
1716: }
1717: PetscBTCreate(numCells, &periodicCells);
1718: for (cell = 0; cell < numCells; ++cell) {
1719: GmshElement *elem = mesh->elements + cell;
1720: for (v = 0; v < elem->numVerts; ++v) {
1721: PetscInt nn = elem->nodes[v];
1722: PetscInt vv = mesh->vertexMap[nn];
1723: if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1724: PetscBTSet(periodicCells, cell); break;
1725: }
1726: }
1727: }
1728: }
1730: /* Setup coordinate DM */
1731: if (coordDim < 0) coordDim = dim;
1732: DMSetCoordinateDim(*dm, coordDim);
1733: DMGetCoordinateDM(*dm, &cdm);
1734: if (highOrder) {
1735: PetscFE fe;
1736: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1737: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1739: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1741: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1742: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1743: DMSetField(cdm, 0, NULL, (PetscObject) fe);
1744: PetscFEDestroy(&fe);
1745: DMCreateDS(cdm);
1746: }
1748: /* Create coordinates */
1749: if (highOrder) {
1751: PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim;
1752: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1753: PetscSection section;
1754: PetscScalar *cellCoords;
1756: DMSetLocalSection(cdm, NULL);
1757: DMGetLocalSection(cdm, &coordSection);
1758: PetscSectionClone(coordSection, §ion);
1759: DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1761: DMCreateLocalVector(cdm, &coordinates);
1762: PetscMalloc1(maxDof, &cellCoords);
1763: for (cell = 0; cell < numCells; ++cell) {
1764: GmshElement *elem = mesh->elements + cell;
1765: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1766: for (n = 0; n < elem->numNodes; ++n) {
1767: const PetscInt node = elem->nodes[lexorder[n]];
1768: for (d = 0; d < coordDim; ++d)
1769: cellCoords[n*coordDim+d] = (PetscReal) coords[node*3+d];
1770: }
1771: DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1772: }
1773: PetscSectionDestroy(§ion);
1774: PetscFree(cellCoords);
1776: } else {
1778: PetscInt *nodeMap;
1779: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1780: PetscScalar *pointCoords;
1782: DMGetLocalSection(cdm, &coordSection);
1783: PetscSectionSetNumFields(coordSection, 1);
1784: PetscSectionSetFieldComponents(coordSection, 0, coordDim);
1785: if (periodic) { /* we need to localize coordinates on cells */
1786: PetscSectionSetChart(coordSection, 0, numCells+numVerts);
1787: } else {
1788: PetscSectionSetChart(coordSection, numCells, numCells+numVerts);
1789: }
1790: if (periodic) {
1791: for (cell = 0; cell < numCells; ++cell) {
1792: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1793: GmshElement *elem = mesh->elements + cell;
1794: PetscInt dof = elem->numVerts * coordDim;
1795: PetscSectionSetDof(coordSection, cell, dof);
1796: PetscSectionSetFieldDof(coordSection, cell, 0, dof);
1797: }
1798: }
1799: }
1800: for (v = numCells; v < numCells+numVerts; ++v) {
1801: PetscSectionSetDof(coordSection, v, coordDim);
1802: PetscSectionSetFieldDof(coordSection, v, 0, coordDim);
1803: }
1804: PetscSectionSetUp(coordSection);
1806: DMCreateLocalVector(cdm, &coordinates);
1807: VecGetArray(coordinates, &pointCoords);
1808: if (periodic) {
1809: for (cell = 0; cell < numCells; ++cell) {
1810: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1811: GmshElement *elem = mesh->elements + cell;
1812: PetscInt off, node;
1813: for (v = 0; v < elem->numVerts; ++v)
1814: cone[v] = elem->nodes[v];
1815: DMPlexReorderCell(cdm, cell, cone);
1816: PetscSectionGetOffset(coordSection, cell, &off);
1817: for (v = 0; v < elem->numVerts; ++v)
1818: for (node = cone[v], d = 0; d < coordDim; ++d)
1819: pointCoords[off++] = (PetscReal) coords[node*3+d];
1820: }
1821: }
1822: }
1823: PetscMalloc1(numVerts, &nodeMap);
1824: for (n = 0; n < numNodes; n++)
1825: if (mesh->vertexMap[n] >= 0)
1826: nodeMap[mesh->vertexMap[n]] = n;
1827: for (v = 0; v < numVerts; ++v) {
1828: PetscInt off, node = nodeMap[v];
1829: PetscSectionGetOffset(coordSection, numCells + v, &off);
1830: for (d = 0; d < coordDim; ++d)
1831: pointCoords[off+d] = (PetscReal) coords[node*3+d];
1832: }
1833: PetscFree(nodeMap);
1834: VecRestoreArray(coordinates, &pointCoords);
1836: }
1838: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1839: VecSetBlockSize(coordinates, coordDim);
1840: DMSetCoordinatesLocal(*dm, coordinates);
1841: VecDestroy(&coordinates);
1842: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1844: GmshMeshDestroy(&mesh);
1845: PetscBTDestroy(&periodicVerts);
1846: PetscBTDestroy(&periodicCells);
1848: if (highOrder && project) {
1849: PetscFE fe;
1850: const char prefix[] = "dm_plex_gmsh_project_";
1851: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1852: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
1854: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1856: GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1857: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1858: DMProjectCoordinates(*dm, fe);
1859: PetscFEDestroy(&fe);
1860: }
1862: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1863: return(0);
1864: }