Actual source code: plexgmsh.c
petsc-3.14.6 2021-03-30
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], slaveDim = -1, slaveTag = -1, masterTag = -1, slaveNode, masterNode;
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", &slaveDim, &slaveTag, &masterTag);
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: slaveDim = ibuf[0]; slaveTag = ibuf[1]; masterTag = ibuf[2];
810: }
811: (void)slaveDim; (void)slaveTag; (void)masterTag; /* 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", &slaveNode, &masterNode);
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: slaveNode = ibuf[0]; masterNode = ibuf[1];
841: }
842: slaveNode = (int) nodeMap[slaveNode];
843: masterNode = (int) nodeMap[masterNode];
844: periodicMap[slaveNode] = masterNode;
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) entityTagMaster(int)
999: numAffine(size_t) value(double) ...
1000: numCorrespondingNodes(size_t)
1001: nodeTag(size_t) nodeTagMaster(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 slaveNode = nodeMap[nodeTags[node*2+0]];
1025: PetscInt masterNode = nodeMap[nodeTags[node*2+1]];
1026: periodicMap[slaveNode] = masterNode;
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)(version*10); /* XXX Should use (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 master 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 slaves) */
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 master */
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 slave */
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: 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) {
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 && 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(*dm, "marker", cone[p], 1);
1643: }
1644: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1645: }
1646: }
1647: }
1649: if (!rank) {
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(*dm, "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(*dm, "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(*dm, "Vertex Sets", vStart + vv, elem->tags[0]);
1686: }
1687: }
1688: }
1689: }
1691: if (periodic) {
1692: PetscBTCreate(numVerts, &periodicVerts);
1693: for (n = 0; n < numNodes; ++n) {
1694: if (mesh->vertexMap[n] >= 0) {
1695: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1696: PetscInt m = mesh->periodMap[n];
1697: ierr= PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1698: ierr= PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1699: }
1700: }
1701: }
1702: PetscBTCreate(numCells, &periodicCells);
1703: for (cell = 0; cell < numCells; ++cell) {
1704: GmshElement *elem = mesh->elements + cell;
1705: for (v = 0; v < elem->numVerts; ++v) {
1706: PetscInt nn = elem->nodes[v];
1707: PetscInt vv = mesh->vertexMap[nn];
1708: if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1709: PetscBTSet(periodicCells, cell); break;
1710: }
1711: }
1712: }
1713: }
1715: /* Setup coordinate DM */
1716: if (coordDim < 0) coordDim = dim;
1717: DMSetCoordinateDim(*dm, coordDim);
1718: DMGetCoordinateDM(*dm, &cdm);
1719: if (highOrder) {
1720: PetscFE fe;
1721: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1722: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1724: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1726: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1727: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1728: DMSetField(cdm, 0, NULL, (PetscObject) fe);
1729: PetscFEDestroy(&fe);
1730: DMCreateDS(cdm);
1731: }
1733: /* Create coordinates */
1734: if (highOrder) {
1736: PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim;
1737: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1738: PetscSection section;
1739: PetscScalar *cellCoords;
1741: DMSetLocalSection(cdm, NULL);
1742: DMGetLocalSection(cdm, &coordSection);
1743: PetscSectionClone(coordSection, §ion);
1744: DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1746: DMCreateLocalVector(cdm, &coordinates);
1747: PetscMalloc1(maxDof, &cellCoords);
1748: for (cell = 0; cell < numCells; ++cell) {
1749: GmshElement *elem = mesh->elements + cell;
1750: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1751: for (n = 0; n < elem->numNodes; ++n) {
1752: const PetscInt node = elem->nodes[lexorder[n]];
1753: for (d = 0; d < coordDim; ++d)
1754: cellCoords[n*coordDim+d] = (PetscReal) coords[node*3+d];
1755: }
1756: DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1757: }
1758: PetscSectionDestroy(§ion);
1759: PetscFree(cellCoords);
1761: } else {
1763: PetscInt *nodeMap;
1764: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1765: PetscScalar *pointCoords;
1767: DMGetLocalSection(cdm, &coordSection);
1768: PetscSectionSetNumFields(coordSection, 1);
1769: PetscSectionSetFieldComponents(coordSection, 0, coordDim);
1770: if (periodic) { /* we need to localize coordinates on cells */
1771: PetscSectionSetChart(coordSection, 0, numCells+numVerts);
1772: } else {
1773: PetscSectionSetChart(coordSection, numCells, numCells+numVerts);
1774: }
1775: if (periodic) {
1776: for (cell = 0; cell < numCells; ++cell) {
1777: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1778: GmshElement *elem = mesh->elements + cell;
1779: PetscInt dof = elem->numVerts * coordDim;
1780: PetscSectionSetDof(coordSection, cell, dof);
1781: PetscSectionSetFieldDof(coordSection, cell, 0, dof);
1782: }
1783: }
1784: }
1785: for (v = numCells; v < numCells+numVerts; ++v) {
1786: PetscSectionSetDof(coordSection, v, coordDim);
1787: PetscSectionSetFieldDof(coordSection, v, 0, coordDim);
1788: }
1789: PetscSectionSetUp(coordSection);
1791: DMCreateLocalVector(cdm, &coordinates);
1792: VecGetArray(coordinates, &pointCoords);
1793: if (periodic) {
1794: for (cell = 0; cell < numCells; ++cell) {
1795: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1796: GmshElement *elem = mesh->elements + cell;
1797: PetscInt off, node;
1798: for (v = 0; v < elem->numVerts; ++v)
1799: cone[v] = elem->nodes[v];
1800: DMPlexReorderCell(cdm, cell, cone);
1801: PetscSectionGetOffset(coordSection, cell, &off);
1802: for (v = 0; v < elem->numVerts; ++v)
1803: for (node = cone[v], d = 0; d < coordDim; ++d)
1804: pointCoords[off++] = (PetscReal) coords[node*3+d];
1805: }
1806: }
1807: }
1808: PetscMalloc1(numVerts, &nodeMap);
1809: for (n = 0; n < numNodes; n++)
1810: if (mesh->vertexMap[n] >= 0)
1811: nodeMap[mesh->vertexMap[n]] = n;
1812: for (v = 0; v < numVerts; ++v) {
1813: PetscInt off, node = nodeMap[v];
1814: PetscSectionGetOffset(coordSection, numCells + v, &off);
1815: for (d = 0; d < coordDim; ++d)
1816: pointCoords[off+d] = (PetscReal) coords[node*3+d];
1817: }
1818: PetscFree(nodeMap);
1819: VecRestoreArray(coordinates, &pointCoords);
1821: }
1823: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1824: VecSetBlockSize(coordinates, coordDim);
1825: DMSetCoordinatesLocal(*dm, coordinates);
1826: VecDestroy(&coordinates);
1827: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1829: GmshMeshDestroy(&mesh);
1830: PetscBTDestroy(&periodicVerts);
1831: PetscBTDestroy(&periodicCells);
1833: if (highOrder && project) {
1834: PetscFE fe;
1835: const char prefix[] = "dm_plex_gmsh_project_";
1836: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1837: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
1839: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1841: GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1842: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1843: DMProjectCoordinates(*dm, fe);
1844: PetscFEDestroy(&fe);
1845: }
1847: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1848: return(0);
1849: }