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: static int *GmshLexOrder_QUA_2_Serendipity(void)
17: {
18: static int Gmsh_LexOrder_QUA_2_Serendipity[9] = {-1};
19: int *lex = Gmsh_LexOrder_QUA_2_Serendipity;
20: if (lex[0] == -1) {
21: /* Vertices */
22: lex[0] = 0;
23: lex[2] = 1;
24: lex[8] = 2;
25: lex[6] = 3;
26: /* Edges */
27: lex[1] = 4;
28: lex[5] = 5;
29: lex[7] = 6;
30: lex[3] = 7;
31: /* Cell */
32: lex[4] = -8;
33: }
34: return lex;
35: }
37: static int *GmshLexOrder_HEX_2_Serendipity(void)
38: {
39: static int Gmsh_LexOrder_HEX_2_Serendipity[27] = {-1};
40: int *lex = Gmsh_LexOrder_HEX_2_Serendipity;
41: if (lex[0] == -1) {
42: /* Vertices */
43: lex[0] = 0;
44: lex[2] = 1;
45: lex[8] = 2;
46: lex[6] = 3;
47: lex[18] = 4;
48: lex[20] = 5;
49: lex[26] = 6;
50: lex[24] = 7;
51: /* Edges */
52: lex[1] = 8;
53: lex[3] = 9;
54: lex[9] = 10;
55: lex[5] = 11;
56: lex[11] = 12;
57: lex[7] = 13;
58: lex[17] = 14;
59: lex[15] = 15;
60: lex[19] = 16;
61: lex[21] = 17;
62: lex[23] = 18;
63: lex[25] = 19;
64: /* Faces */
65: lex[4] = -20;
66: lex[10] = -21;
67: lex[12] = -22;
68: lex[14] = -23;
69: lex[16] = -24;
70: lex[22] = -25;
71: /* Cell */
72: lex[13] = -26;
73: }
74: return lex;
75: }
77: #define GMSH_LEXORDER_LIST(T) \
78: GMSH_LEXORDER_ITEM(T, 1) \
79: GMSH_LEXORDER_ITEM(T, 2) \
80: GMSH_LEXORDER_ITEM(T, 3) \
81: GMSH_LEXORDER_ITEM(T, 4) \
82: GMSH_LEXORDER_ITEM(T, 5) \
83: GMSH_LEXORDER_ITEM(T, 6) \
84: GMSH_LEXORDER_ITEM(T, 7) \
85: GMSH_LEXORDER_ITEM(T, 8) \
86: GMSH_LEXORDER_ITEM(T, 9) \
87: GMSH_LEXORDER_ITEM(T, 10)
89: GMSH_LEXORDER_ITEM(VTX, 0)
90: GMSH_LEXORDER_LIST(SEG)
91: GMSH_LEXORDER_LIST(TRI)
92: GMSH_LEXORDER_LIST(QUA)
93: GMSH_LEXORDER_LIST(TET)
94: GMSH_LEXORDER_LIST(HEX)
95: GMSH_LEXORDER_LIST(PRI)
96: GMSH_LEXORDER_LIST(PYR)
98: typedef enum {
99: GMSH_VTX = 0,
100: GMSH_SEG = 1,
101: GMSH_TRI = 2,
102: GMSH_QUA = 3,
103: GMSH_TET = 4,
104: GMSH_HEX = 5,
105: GMSH_PRI = 6,
106: GMSH_PYR = 7,
107: GMSH_NUM_POLYTOPES = 8
108: } GmshPolytopeType;
110: typedef struct {
111: int cellType;
112: int polytope;
113: int dim;
114: int order;
115: int numVerts;
116: int numNodes;
117: int *(*lexorder)(void);
118: } GmshCellInfo;
120: // clang-format off
121: #define GmshCellEntry(cellType, polytope, dim, order) {cellType, GMSH_##polytope, dim, order, GmshNumNodes_##polytope(1), GmshNumNodes_##polytope(order), GmshLexOrder_##polytope##_##order}
122: // clang-format on
124: static const GmshCellInfo GmshCellTable[] = {
125: GmshCellEntry(15, VTX, 0, 0),
127: GmshCellEntry(1, SEG, 1, 1), GmshCellEntry(8, SEG, 1, 2), GmshCellEntry(26, SEG, 1, 3),
128: GmshCellEntry(27, SEG, 1, 4), GmshCellEntry(28, SEG, 1, 5), GmshCellEntry(62, SEG, 1, 6),
129: GmshCellEntry(63, SEG, 1, 7), GmshCellEntry(64, SEG, 1, 8), GmshCellEntry(65, SEG, 1, 9),
130: GmshCellEntry(66, SEG, 1, 10),
132: GmshCellEntry(2, TRI, 2, 1), GmshCellEntry(9, TRI, 2, 2), GmshCellEntry(21, TRI, 2, 3),
133: GmshCellEntry(23, TRI, 2, 4), GmshCellEntry(25, TRI, 2, 5), GmshCellEntry(42, TRI, 2, 6),
134: GmshCellEntry(43, TRI, 2, 7), GmshCellEntry(44, TRI, 2, 8), GmshCellEntry(45, TRI, 2, 9),
135: GmshCellEntry(46, TRI, 2, 10),
137: GmshCellEntry(3, QUA, 2, 1), GmshCellEntry(10, QUA, 2, 2), {16, GMSH_QUA, 2, 2, 4, 8, GmshLexOrder_QUA_2_Serendipity},
138: GmshCellEntry(36, QUA, 2, 3), GmshCellEntry(37, QUA, 2, 4), GmshCellEntry(38, QUA, 2, 5),
139: GmshCellEntry(47, QUA, 2, 6), GmshCellEntry(48, QUA, 2, 7), GmshCellEntry(49, QUA, 2, 8),
140: GmshCellEntry(50, QUA, 2, 9), GmshCellEntry(51, QUA, 2, 10),
142: GmshCellEntry(4, TET, 3, 1), GmshCellEntry(11, TET, 3, 2), GmshCellEntry(29, TET, 3, 3),
143: GmshCellEntry(30, TET, 3, 4), GmshCellEntry(31, TET, 3, 5), GmshCellEntry(71, TET, 3, 6),
144: GmshCellEntry(72, TET, 3, 7), GmshCellEntry(73, TET, 3, 8), GmshCellEntry(74, TET, 3, 9),
145: GmshCellEntry(75, TET, 3, 10),
147: GmshCellEntry(5, HEX, 3, 1), GmshCellEntry(12, HEX, 3, 2), {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity},
148: GmshCellEntry(92, HEX, 3, 3), GmshCellEntry(93, HEX, 3, 4), GmshCellEntry(94, HEX, 3, 5),
149: GmshCellEntry(95, HEX, 3, 6), GmshCellEntry(96, HEX, 3, 7), GmshCellEntry(97, HEX, 3, 8),
150: GmshCellEntry(98, HEX, 3, 9), GmshCellEntry(-1, HEX, 3, 10),
152: GmshCellEntry(6, PRI, 3, 1), GmshCellEntry(13, PRI, 3, 2), GmshCellEntry(90, PRI, 3, 3),
153: GmshCellEntry(91, PRI, 3, 4), GmshCellEntry(106, PRI, 3, 5), GmshCellEntry(107, PRI, 3, 6),
154: GmshCellEntry(108, PRI, 3, 7), GmshCellEntry(109, PRI, 3, 8), GmshCellEntry(110, PRI, 3, 9),
155: GmshCellEntry(-1, PRI, 3, 10),
157: GmshCellEntry(7, PYR, 3, 1), GmshCellEntry(14, PYR, 3, 2), GmshCellEntry(118, PYR, 3, 3),
158: GmshCellEntry(119, PYR, 3, 4), GmshCellEntry(120, PYR, 3, 5), GmshCellEntry(121, PYR, 3, 6),
159: GmshCellEntry(122, PYR, 3, 7), GmshCellEntry(123, PYR, 3, 8), GmshCellEntry(124, PYR, 3, 9),
160: GmshCellEntry(-1, PYR, 3, 10)
162: #if 0
163: {20, GMSH_TRI, 2, 3, 3, 9, NULL},
164: {18, GMSH_PRI, 3, 2, 6, 15, NULL},
165: {19, GMSH_PYR, 3, 2, 5, 13, NULL},
166: #endif
167: };
169: static GmshCellInfo GmshCellMap[150];
171: static PetscErrorCode GmshCellInfoSetUp(void)
172: {
173: size_t i, n;
174: static PetscBool called = PETSC_FALSE;
176: PetscFunctionBegin;
177: if (called) PetscFunctionReturn(PETSC_SUCCESS);
178: called = PETSC_TRUE;
179: n = PETSC_STATIC_ARRAY_LENGTH(GmshCellMap);
180: for (i = 0; i < n; ++i) {
181: GmshCellMap[i].cellType = -1;
182: GmshCellMap[i].polytope = -1;
183: }
184: n = PETSC_STATIC_ARRAY_LENGTH(GmshCellTable);
185: for (i = 0; i < n; ++i) {
186: if (GmshCellTable[i].cellType <= 0) continue;
187: GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i];
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: #define GmshCellTypeCheck(ct) \
193: PetscMacroReturnStandard(const int _ct_ = (int)ct; PetscCheck(_ct_ >= 0 && _ct_ < (int)PETSC_STATIC_ARRAY_LENGTH(GmshCellMap), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Gmsh element type %d", _ct_); PetscCheck(GmshCellMap[_ct_].cellType == _ct_, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \
194: PetscCheck(GmshCellMap[_ct_].polytope != -1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_);)
196: typedef struct {
197: PetscViewer viewer;
198: int fileFormat;
199: int dataSize;
200: PetscBool binary;
201: PetscBool byteSwap;
202: size_t wlen;
203: void *wbuf;
204: size_t slen;
205: void *sbuf;
206: PetscInt *nbuf;
207: PetscInt nodeStart;
208: PetscInt nodeEnd;
209: PetscInt *nodeMap;
210: } GmshFile;
212: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
213: {
214: size_t size = count * eltsize;
216: PetscFunctionBegin;
217: if (gmsh->wlen < size) {
218: PetscCall(PetscFree(gmsh->wbuf));
219: PetscCall(PetscMalloc(size, &gmsh->wbuf));
220: gmsh->wlen = size;
221: }
222: *(void **)buf = size ? gmsh->wbuf : NULL;
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
227: {
228: size_t dataSize = (size_t)gmsh->dataSize;
229: size_t size = count * dataSize;
231: PetscFunctionBegin;
232: if (gmsh->slen < size) {
233: PetscCall(PetscFree(gmsh->sbuf));
234: PetscCall(PetscMalloc(size, &gmsh->sbuf));
235: gmsh->slen = size;
236: }
237: *(void **)buf = size ? gmsh->sbuf : NULL;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
242: {
243: PetscFunctionBegin;
244: PetscCall(PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype));
245: if (gmsh->byteSwap) PetscCall(PetscByteSwap(buf, dtype, count));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
250: {
251: PetscFunctionBegin;
252: PetscCall(PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING));
253: PetscFunctionReturn(PETSC_SUCCESS);
254: }
256: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
257: {
258: PetscFunctionBegin;
259: PetscCall(PetscStrcmp(line, Section, match));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
264: {
265: PetscBool match;
267: PetscFunctionBegin;
268: PetscCall(GmshMatch(gmsh, Section, line, &match));
269: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s, not %s", Section, line);
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
274: {
275: PetscBool match;
277: PetscFunctionBegin;
278: while (PETSC_TRUE) {
279: PetscCall(GmshReadString(gmsh, line, 1));
280: PetscCall(GmshMatch(gmsh, "$Comments", line, &match));
281: if (!match) break;
282: while (PETSC_TRUE) {
283: PetscCall(GmshReadString(gmsh, line, 1));
284: PetscCall(GmshMatch(gmsh, "$EndComments", line, &match));
285: if (match) break;
286: }
287: }
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
292: {
293: PetscFunctionBegin;
294: PetscCall(GmshReadString(gmsh, line, 1));
295: PetscCall(GmshExpect(gmsh, EndSection, line));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
300: {
301: PetscInt i;
302: size_t dataSize = (size_t)gmsh->dataSize;
304: PetscFunctionBegin;
305: if (dataSize == sizeof(PetscInt)) {
306: PetscCall(GmshRead(gmsh, buf, count, PETSC_INT));
307: } else if (dataSize == sizeof(int)) {
308: int *ibuf = NULL;
309: PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
310: PetscCall(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: PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
315: PetscCall(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: PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
320: PetscCall(GmshRead(gmsh, ibuf, count, PETSC_INT64));
321: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
322: }
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
327: {
328: PetscFunctionBegin;
329: PetscCall(GmshRead(gmsh, buf, count, PETSC_ENUM));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
334: {
335: PetscFunctionBegin;
336: PetscCall(GmshRead(gmsh, buf, count, PETSC_DOUBLE));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: #define GMSH_MAX_TAGS 16
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[GMSH_MAX_TAGS]; /* 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;
359: PetscFunctionBegin;
360: PetscCall(PetscNew(entities));
361: for (dim = 0; dim < 4; ++dim) {
362: PetscCall(PetscCalloc1(count[dim], &(*entities)->entity[dim]));
363: PetscCall(PetscHMapICreate(&(*entities)->entityMap[dim]));
364: }
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
369: {
370: PetscInt dim;
372: PetscFunctionBegin;
373: if (!*entities) PetscFunctionReturn(PETSC_SUCCESS);
374: for (dim = 0; dim < 4; ++dim) {
375: PetscCall(PetscFree((*entities)->entity[dim]));
376: PetscCall(PetscHMapIDestroy(&(*entities)->entityMap[dim]));
377: }
378: PetscCall(PetscFree(*entities));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity **entity)
383: {
384: PetscFunctionBegin;
385: PetscCall(PetscHMapISet(entities->entityMap[dim], eid, index));
386: entities->entity[dim][index].dim = dim;
387: entities->entity[dim][index].id = eid;
388: if (entity) *entity = &entities->entity[dim][index];
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity **entity)
393: {
394: PetscInt index;
396: PetscFunctionBegin;
397: PetscCall(PetscHMapIGet(entities->entityMap[dim], eid, &index));
398: *entity = &entities->entity[dim][index];
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: typedef struct {
403: PetscInt *id; /* Node IDs */
404: double *xyz; /* Coordinates */
405: PetscInt *tag; /* Physical tag */
406: } GmshNodes;
408: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
409: {
410: PetscFunctionBegin;
411: PetscCall(PetscNew(nodes));
412: PetscCall(PetscMalloc1(count * 1, &(*nodes)->id));
413: PetscCall(PetscMalloc1(count * 3, &(*nodes)->xyz));
414: PetscCall(PetscMalloc1(count * GMSH_MAX_TAGS, &(*nodes)->tag));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
419: {
420: PetscFunctionBegin;
421: if (!*nodes) PetscFunctionReturn(PETSC_SUCCESS);
422: PetscCall(PetscFree((*nodes)->id));
423: PetscCall(PetscFree((*nodes)->xyz));
424: PetscCall(PetscFree((*nodes)->tag));
425: PetscCall(PetscFree(*nodes));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: typedef struct {
430: PetscInt id; /* Element ID */
431: PetscInt dim; /* Dimension */
432: PetscInt cellType; /* Cell type */
433: PetscInt numVerts; /* Size of vertex array */
434: PetscInt numNodes; /* Size of node array */
435: PetscInt *nodes; /* Vertex/Node array */
436: PetscInt numTags; /* Size of physical tag array */
437: int tags[GMSH_MAX_TAGS]; /* Physical tag array */
438: } GmshElement;
440: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
441: {
442: PetscFunctionBegin;
443: PetscCall(PetscCalloc1(count, elements));
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
448: {
449: PetscFunctionBegin;
450: if (!*elements) PetscFunctionReturn(PETSC_SUCCESS);
451: PetscCall(PetscFree(*elements));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: typedef struct {
456: PetscInt dim;
457: PetscInt order;
458: GmshEntities *entities;
459: PetscInt numNodes;
460: GmshNodes *nodelist;
461: PetscInt numElems;
462: GmshElement *elements;
463: PetscInt numVerts;
464: PetscInt numCells;
465: PetscInt *periodMap;
466: PetscInt *vertexMap;
467: PetscSegBuffer segbuf;
468: PetscInt numRegions;
469: PetscInt *regionDims;
470: PetscInt *regionTags;
471: char **regionNames;
472: } GmshMesh;
474: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
475: {
476: PetscFunctionBegin;
477: PetscCall(PetscNew(mesh));
478: PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf));
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
483: {
484: PetscInt r;
486: PetscFunctionBegin;
487: if (!*mesh) PetscFunctionReturn(PETSC_SUCCESS);
488: PetscCall(GmshEntitiesDestroy(&(*mesh)->entities));
489: PetscCall(GmshNodesDestroy(&(*mesh)->nodelist));
490: PetscCall(GmshElementsDestroy(&(*mesh)->elements));
491: PetscCall(PetscFree((*mesh)->periodMap));
492: PetscCall(PetscFree((*mesh)->vertexMap));
493: PetscCall(PetscSegBufferDestroy(&(*mesh)->segbuf));
494: for (r = 0; r < (*mesh)->numRegions; ++r) PetscCall(PetscFree((*mesh)->regionNames[r]));
495: PetscCall(PetscFree3((*mesh)->regionDims, (*mesh)->regionTags, (*mesh)->regionNames));
496: PetscCall(PetscFree(*mesh));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
501: {
502: PetscViewer viewer = gmsh->viewer;
503: PetscBool byteSwap = gmsh->byteSwap;
504: char line[PETSC_MAX_PATH_LEN];
505: int n, t, num, nid, snum;
506: GmshNodes *nodes;
508: PetscFunctionBegin;
509: PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
510: snum = sscanf(line, "%d", &num);
511: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
512: PetscCall(GmshNodesCreate(num, &nodes));
513: mesh->numNodes = num;
514: mesh->nodelist = nodes;
515: for (n = 0; n < num; ++n) {
516: double *xyz = nodes->xyz + n * 3;
517: PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
518: PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
519: if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
520: if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
521: nodes->id[n] = nid;
522: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
523: }
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
528: file contents multiple times to figure out the true number of cells and facets
529: in the given mesh. To make this more efficient we read the file contents only
530: once and store them in memory, while determining the true number of cells. */
531: static PetscErrorCode GmshReadElements_v22(GmshFile *gmsh, GmshMesh *mesh)
532: {
533: PetscViewer viewer = gmsh->viewer;
534: PetscBool binary = gmsh->binary;
535: PetscBool byteSwap = gmsh->byteSwap;
536: char line[PETSC_MAX_PATH_LEN];
537: int i, c, p, num, ibuf[1 + 4 + 1000], snum;
538: int cellType, numElem, numVerts, numNodes, numTags;
539: GmshElement *elements;
540: PetscInt *nodeMap = gmsh->nodeMap;
542: PetscFunctionBegin;
543: PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
544: snum = sscanf(line, "%d", &num);
545: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
546: PetscCall(GmshElementsCreate(num, &elements));
547: mesh->numElems = num;
548: mesh->elements = elements;
549: for (c = 0; c < num;) {
550: PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
551: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));
553: cellType = binary ? ibuf[0] : ibuf[1];
554: numElem = binary ? ibuf[1] : 1;
555: numTags = ibuf[2];
557: PetscCall(GmshCellTypeCheck(cellType));
558: numVerts = GmshCellMap[cellType].numVerts;
559: numNodes = GmshCellMap[cellType].numNodes;
561: for (i = 0; i < numElem; ++i, ++c) {
562: GmshElement *element = elements + c;
563: const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
564: const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
565: PetscCall(PetscViewerRead(viewer, ibuf + off, nint, NULL, PETSC_ENUM));
566: if (byteSwap) PetscCall(PetscByteSwap(ibuf + off, PETSC_ENUM, nint));
567: element->id = id[0];
568: element->dim = GmshCellMap[cellType].dim;
569: element->cellType = cellType;
570: element->numVerts = numVerts;
571: element->numNodes = numNodes;
572: element->numTags = PetscMin(numTags, GMSH_MAX_TAGS);
573: PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
574: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
575: for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
576: }
577: }
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*
582: $Entities
583: numPoints(unsigned long) numCurves(unsigned long)
584: numSurfaces(unsigned long) numVolumes(unsigned long)
585: // points
586: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
587: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
588: numPhysicals(unsigned long) phyisicalTag[...](int)
589: ...
590: // curves
591: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
592: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
593: numPhysicals(unsigned long) physicalTag[...](int)
594: numBREPVert(unsigned long) tagBREPVert[...](int)
595: ...
596: // surfaces
597: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
598: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
599: numPhysicals(unsigned long) physicalTag[...](int)
600: numBREPCurve(unsigned long) tagBREPCurve[...](int)
601: ...
602: // volumes
603: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
604: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
605: numPhysicals(unsigned long) physicalTag[...](int)
606: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
607: ...
608: $EndEntities
609: */
610: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
611: {
612: PetscViewer viewer = gmsh->viewer;
613: PetscBool byteSwap = gmsh->byteSwap;
614: long index, num, lbuf[4];
615: int dim, eid, numTags, *ibuf, t;
616: PetscInt count[4], i;
617: GmshEntity *entity = NULL;
619: PetscFunctionBegin;
620: PetscCall(PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG));
621: if (byteSwap) PetscCall(PetscByteSwap(lbuf, PETSC_LONG, 4));
622: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
623: PetscCall(GmshEntitiesCreate(count, &mesh->entities));
624: for (dim = 0; dim < 4; ++dim) {
625: for (index = 0; index < count[dim]; ++index) {
626: PetscCall(PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM));
627: if (byteSwap) PetscCall(PetscByteSwap(&eid, PETSC_ENUM, 1));
628: PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
629: PetscCall(PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE));
630: if (byteSwap) PetscCall(PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6));
631: PetscCall(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG));
632: if (byteSwap) PetscCall(PetscByteSwap(&num, PETSC_LONG, 1));
633: PetscCall(GmshBufferGet(gmsh, num, sizeof(int), &ibuf));
634: PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
635: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
636: entity->numTags = numTags = (int)PetscMin(num, GMSH_MAX_TAGS);
637: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
638: if (dim == 0) continue;
639: PetscCall(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG));
640: if (byteSwap) PetscCall(PetscByteSwap(&num, PETSC_LONG, 1));
641: PetscCall(GmshBufferGet(gmsh, num, sizeof(int), &ibuf));
642: PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
643: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
644: }
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*
650: $Nodes
651: numEntityBlocks(unsigned long) numNodes(unsigned long)
652: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
653: tag(int) x(double) y(double) z(double)
654: ...
655: ...
656: $EndNodes
657: */
658: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
659: {
660: PetscViewer viewer = gmsh->viewer;
661: PetscBool byteSwap = gmsh->byteSwap;
662: long block, node, n, t, numEntityBlocks, numTotalNodes, numNodes;
663: int info[3], nid;
664: GmshNodes *nodes;
666: PetscFunctionBegin;
667: PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
668: if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
669: PetscCall(PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG));
670: if (byteSwap) PetscCall(PetscByteSwap(&numTotalNodes, PETSC_LONG, 1));
671: PetscCall(GmshNodesCreate(numTotalNodes, &nodes));
672: mesh->numNodes = numTotalNodes;
673: mesh->nodelist = nodes;
674: for (n = 0, block = 0; block < numEntityBlocks; ++block) {
675: PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
676: PetscCall(PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG));
677: if (byteSwap) PetscCall(PetscByteSwap(&numNodes, PETSC_LONG, 1));
678: if (gmsh->binary) {
679: size_t nbytes = sizeof(int) + 3 * sizeof(double);
680: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible uninitialized value */
681: PetscCall(GmshBufferGet(gmsh, numNodes, nbytes, &cbuf));
682: PetscCall(PetscViewerRead(viewer, cbuf, numNodes * nbytes, NULL, PETSC_CHAR));
683: for (node = 0; node < numNodes; ++node, ++n) {
684: char *cnid = cbuf + node * nbytes, *cxyz = cnid + sizeof(int);
685: double *xyz = nodes->xyz + n * 3;
686: if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cnid, PETSC_ENUM, 1));
687: if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cxyz, PETSC_DOUBLE, 3));
688: PetscCall(PetscMemcpy(&nid, cnid, sizeof(int)));
689: PetscCall(PetscMemcpy(xyz, cxyz, 3 * sizeof(double)));
690: if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
691: if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
692: nodes->id[n] = nid;
693: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
694: }
695: } else {
696: for (node = 0; node < numNodes; ++node, ++n) {
697: double *xyz = nodes->xyz + n * 3;
698: PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
699: PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
700: if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
701: if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
702: nodes->id[n] = nid;
703: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
704: }
705: }
706: }
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: /*
711: $Elements
712: numEntityBlocks(unsigned long) numElements(unsigned long)
713: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
714: tag(int) numVert[...](int)
715: ...
716: ...
717: $EndElements
718: */
719: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
720: {
721: PetscViewer viewer = gmsh->viewer;
722: PetscBool byteSwap = gmsh->byteSwap;
723: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
724: int p, info[3], *ibuf = NULL;
725: int eid, dim, cellType, numVerts, numNodes, numTags;
726: GmshEntity *entity = NULL;
727: GmshElement *elements;
728: PetscInt *nodeMap = gmsh->nodeMap;
730: PetscFunctionBegin;
731: PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
732: if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
733: PetscCall(PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG));
734: if (byteSwap) PetscCall(PetscByteSwap(&numTotalElements, PETSC_LONG, 1));
735: PetscCall(GmshElementsCreate(numTotalElements, &elements));
736: mesh->numElems = numTotalElements;
737: mesh->elements = elements;
738: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
739: PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
740: if (byteSwap) PetscCall(PetscByteSwap(info, PETSC_ENUM, 3));
741: eid = info[0];
742: dim = info[1];
743: cellType = info[2];
744: PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
745: PetscCall(GmshCellTypeCheck(cellType));
746: numVerts = GmshCellMap[cellType].numVerts;
747: numNodes = GmshCellMap[cellType].numNodes;
748: numTags = entity->numTags;
749: PetscCall(PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG));
750: if (byteSwap) PetscCall(PetscByteSwap(&numElements, PETSC_LONG, 1));
751: PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numElements, sizeof(int), &ibuf));
752: PetscCall(PetscViewerRead(viewer, ibuf, (1 + numNodes) * numElements, NULL, PETSC_ENUM));
753: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, (1 + numNodes) * numElements));
754: for (elem = 0; elem < numElements; ++elem, ++c) {
755: GmshElement *element = elements + c;
756: const int *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
757: element->id = id[0];
758: element->dim = dim;
759: element->cellType = cellType;
760: element->numVerts = numVerts;
761: element->numNodes = numNodes;
762: element->numTags = numTags;
763: PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
764: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
765: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
766: }
767: }
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
772: {
773: PetscViewer viewer = gmsh->viewer;
774: int fileFormat = gmsh->fileFormat;
775: PetscBool binary = gmsh->binary;
776: PetscBool byteSwap = gmsh->byteSwap;
777: int numPeriodic, snum, i;
778: char line[PETSC_MAX_PATH_LEN];
779: PetscInt *nodeMap = gmsh->nodeMap;
781: PetscFunctionBegin;
782: if (fileFormat == 22 || !binary) {
783: PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
784: snum = sscanf(line, "%d", &numPeriodic);
785: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
786: } else {
787: PetscCall(PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM));
788: if (byteSwap) PetscCall(PetscByteSwap(&numPeriodic, PETSC_ENUM, 1));
789: }
790: for (i = 0; i < numPeriodic; i++) {
791: int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
792: long j, nNodes;
793: double affine[16];
795: if (fileFormat == 22 || !binary) {
796: PetscCall(PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING));
797: snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
798: PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
799: } else {
800: PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
801: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));
802: correspondingDim = ibuf[0];
803: correspondingTag = ibuf[1];
804: primaryTag = ibuf[2];
805: }
806: (void)correspondingDim;
807: (void)correspondingTag;
808: (void)primaryTag; /* unused */
810: if (fileFormat == 22 || !binary) {
811: PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
812: snum = sscanf(line, "%ld", &nNodes);
813: if (snum != 1) { /* discard transformation and try again */
814: PetscCall(PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING));
815: PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
816: snum = sscanf(line, "%ld", &nNodes);
817: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
818: }
819: } else {
820: PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
821: if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
822: if (nNodes == -1) { /* discard transformation and try again */
823: PetscCall(PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE));
824: PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
825: if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
826: }
827: }
829: for (j = 0; j < nNodes; j++) {
830: if (fileFormat == 22 || !binary) {
831: PetscCall(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING));
832: snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
833: PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
834: } else {
835: PetscCall(PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM));
836: if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 2));
837: correspondingNode = ibuf[0];
838: primaryNode = ibuf[1];
839: }
840: correspondingNode = (int)nodeMap[correspondingNode];
841: primaryNode = (int)nodeMap[primaryNode];
842: periodicMap[correspondingNode] = primaryNode;
843: }
844: }
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
849: $Entities
850: numPoints(size_t) numCurves(size_t)
851: numSurfaces(size_t) numVolumes(size_t)
852: pointTag(int) X(double) Y(double) Z(double)
853: numPhysicalTags(size_t) physicalTag(int) ...
854: ...
855: curveTag(int) minX(double) minY(double) minZ(double)
856: maxX(double) maxY(double) maxZ(double)
857: numPhysicalTags(size_t) physicalTag(int) ...
858: numBoundingPoints(size_t) pointTag(int) ...
859: ...
860: surfaceTag(int) minX(double) minY(double) minZ(double)
861: maxX(double) maxY(double) maxZ(double)
862: numPhysicalTags(size_t) physicalTag(int) ...
863: numBoundingCurves(size_t) curveTag(int) ...
864: ...
865: volumeTag(int) minX(double) minY(double) minZ(double)
866: maxX(double) maxY(double) maxZ(double)
867: numPhysicalTags(size_t) physicalTag(int) ...
868: numBoundngSurfaces(size_t) surfaceTag(int) ...
869: ...
870: $EndEntities
871: */
872: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
873: {
874: PetscInt count[4], index, numTags, i;
875: int dim, eid, *tags = NULL;
876: GmshEntity *entity = NULL;
878: PetscFunctionBegin;
879: PetscCall(GmshReadSize(gmsh, count, 4));
880: PetscCall(GmshEntitiesCreate(count, &mesh->entities));
881: for (dim = 0; dim < 4; ++dim) {
882: for (index = 0; index < count[dim]; ++index) {
883: PetscCall(GmshReadInt(gmsh, &eid, 1));
884: PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
885: PetscCall(GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6));
886: PetscCall(GmshReadSize(gmsh, &numTags, 1));
887: PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
888: PetscCall(GmshReadInt(gmsh, tags, numTags));
889: PetscCheck(numTags <= GMSH_MAX_TAGS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PETSc currently supports up to %" PetscInt_FMT " tags per entity, not %" PetscInt_FMT, (PetscInt)GMSH_MAX_TAGS, numTags);
890: entity->numTags = numTags;
891: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
892: if (dim == 0) continue;
893: PetscCall(GmshReadSize(gmsh, &numTags, 1));
894: PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
895: PetscCall(GmshReadInt(gmsh, tags, numTags));
896: /* Currently, we do not save the ids for the bounding entities */
897: }
898: }
899: PetscFunctionReturn(PETSC_SUCCESS);
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], dim, eid, parametric;
919: PetscInt sizes[4], numEntityBlocks, numTags, t, numNodes, numNodesBlock = 0, block, node, n;
920: GmshEntity *entity = NULL;
921: GmshNodes *nodes;
923: PetscFunctionBegin;
924: PetscCall(GmshReadSize(gmsh, sizes, 4));
925: numEntityBlocks = sizes[0];
926: numNodes = sizes[1];
927: PetscCall(GmshNodesCreate(numNodes, &nodes));
928: mesh->numNodes = numNodes;
929: mesh->nodelist = nodes;
930: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
931: PetscCall(GmshReadInt(gmsh, info, 3));
932: dim = info[0];
933: eid = info[1];
934: parametric = info[2];
935: PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
936: numTags = entity->numTags;
937: PetscCheck(!parametric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
938: PetscCall(GmshReadSize(gmsh, &numNodesBlock, 1));
939: PetscCall(GmshReadSize(gmsh, nodes->id + node, numNodesBlock));
940: PetscCall(GmshReadDouble(gmsh, nodes->xyz + node * 3, numNodesBlock * 3));
941: for (n = 0; n < numNodesBlock; ++n) {
942: PetscInt *tags = &nodes->tag[node * GMSH_MAX_TAGS];
944: for (t = 0; t < numTags; ++t) tags[n * GMSH_MAX_TAGS + t] = entity->tags[t];
945: for (t = numTags; t < GMSH_MAX_TAGS; ++t) tags[n * GMSH_MAX_TAGS + t] = -1;
946: }
947: }
948: gmsh->nodeStart = sizes[2];
949: gmsh->nodeEnd = sizes[3] + 1;
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
954: $Elements
955: numEntityBlocks(size_t) numElements(size_t)
956: minElementTag(size_t) maxElementTag(size_t)
957: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
958: elementTag(size_t) nodeTag(size_t) ...
959: ...
960: ...
961: $EndElements
962: */
963: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
964: {
965: int info[3], eid, dim, cellType;
966: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
967: GmshEntity *entity = NULL;
968: GmshElement *elements;
969: PetscInt *nodeMap = gmsh->nodeMap;
971: PetscFunctionBegin;
972: PetscCall(GmshReadSize(gmsh, sizes, 4));
973: numEntityBlocks = sizes[0];
974: numElements = sizes[1];
975: PetscCall(GmshElementsCreate(numElements, &elements));
976: mesh->numElems = numElements;
977: mesh->elements = elements;
978: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
979: PetscCall(GmshReadInt(gmsh, info, 3));
980: dim = info[0];
981: eid = info[1];
982: cellType = info[2];
983: PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
984: PetscCall(GmshCellTypeCheck(cellType));
985: numVerts = GmshCellMap[cellType].numVerts;
986: numNodes = GmshCellMap[cellType].numNodes;
987: numTags = entity->numTags;
988: PetscCall(GmshReadSize(gmsh, &numBlockElements, 1));
989: PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numBlockElements, sizeof(PetscInt), &ibuf));
990: PetscCall(GmshReadSize(gmsh, ibuf, (1 + numNodes) * numBlockElements));
991: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
992: GmshElement *element = elements + c;
993: const PetscInt *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
994: element->id = id[0];
995: element->dim = dim;
996: element->cellType = cellType;
997: element->numVerts = numVerts;
998: element->numNodes = numNodes;
999: element->numTags = numTags;
1000: PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
1001: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
1002: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
1003: }
1004: }
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1009: $Periodic
1010: numPeriodicLinks(size_t)
1011: entityDim(int) entityTag(int) entityTagPrimary(int)
1012: numAffine(size_t) value(double) ...
1013: numCorrespondingNodes(size_t)
1014: nodeTag(size_t) nodeTagPrimary(size_t)
1015: ...
1016: ...
1017: $EndPeriodic
1018: */
1019: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
1020: {
1021: int info[3];
1022: double dbuf[16];
1023: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
1024: PetscInt *nodeMap = gmsh->nodeMap;
1026: PetscFunctionBegin;
1027: PetscCall(GmshReadSize(gmsh, &numPeriodicLinks, 1));
1028: for (link = 0; link < numPeriodicLinks; ++link) {
1029: PetscCall(GmshReadInt(gmsh, info, 3));
1030: PetscCall(GmshReadSize(gmsh, &numAffine, 1));
1031: PetscCall(GmshReadDouble(gmsh, dbuf, numAffine));
1032: PetscCall(GmshReadSize(gmsh, &numCorrespondingNodes, 1));
1033: PetscCall(GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags));
1034: PetscCall(GmshReadSize(gmsh, nodeTags, numCorrespondingNodes * 2));
1035: for (node = 0; node < numCorrespondingNodes; ++node) {
1036: PetscInt correspondingNode = nodeMap[nodeTags[node * 2 + 0]];
1037: PetscInt primaryNode = nodeMap[nodeTags[node * 2 + 1]];
1038: periodicMap[correspondingNode] = primaryNode;
1039: }
1040: }
1041: PetscFunctionReturn(PETSC_SUCCESS);
1042: }
1044: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1045: $MeshFormat // same as MSH version 2
1046: version(ASCII double; currently 4.1)
1047: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1048: data-size(ASCII int; sizeof(size_t))
1049: < int with value one; only in binary mode, to detect endianness >
1050: $EndMeshFormat
1051: */
1052: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1053: {
1054: char line[PETSC_MAX_PATH_LEN];
1055: int snum, fileType, fileFormat, dataSize, checkEndian;
1056: float version;
1058: PetscFunctionBegin;
1059: PetscCall(GmshReadString(gmsh, line, 3));
1060: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1061: fileFormat = (int)roundf(version * 10);
1062: PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1063: PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1064: PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1065: PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1066: PetscCheck(!gmsh->binary || fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII");
1067: PetscCheck(gmsh->binary || !fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary");
1068: PetscCheck(fileFormat > 40 || dataSize == sizeof(double), PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1069: PetscCheck(fileFormat < 41 || dataSize == sizeof(int) || dataSize == sizeof(PetscInt64), PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1070: gmsh->fileFormat = fileFormat;
1071: gmsh->dataSize = dataSize;
1072: gmsh->byteSwap = PETSC_FALSE;
1073: if (gmsh->binary) {
1074: PetscCall(GmshReadInt(gmsh, &checkEndian, 1));
1075: if (checkEndian != 1) {
1076: PetscCall(PetscByteSwap(&checkEndian, PETSC_ENUM, 1));
1077: PetscCheck(checkEndian == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line);
1078: gmsh->byteSwap = PETSC_TRUE;
1079: }
1080: }
1081: PetscFunctionReturn(PETSC_SUCCESS);
1082: }
1084: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1085: Neper: https://neper.info/ adds this section
1087: $MeshVersion
1088: <major>.<minor>,<patch>
1089: $EndMeshVersion
1090: */
1091: static PetscErrorCode GmshReadMeshVersion(GmshFile *gmsh)
1092: {
1093: char line[PETSC_MAX_PATH_LEN];
1094: int snum, major, minor, patch;
1096: PetscFunctionBegin;
1097: PetscCall(GmshReadString(gmsh, line, 1));
1098: snum = sscanf(line, "%d.%d.%d", &major, &minor, &patch);
1099: PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1100: PetscFunctionReturn(PETSC_SUCCESS);
1101: }
1103: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1104: Neper: https://neper.info/ adds this section
1106: $Domain
1107: <shape>
1108: $EndDomain
1109: */
1110: static PetscErrorCode GmshReadMeshDomain(GmshFile *gmsh)
1111: {
1112: char line[PETSC_MAX_PATH_LEN];
1114: PetscFunctionBegin;
1115: PetscCall(GmshReadString(gmsh, line, 1));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: /*
1120: PhysicalNames
1121: numPhysicalNames(ASCII int)
1122: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1123: ...
1124: $EndPhysicalNames
1125: */
1126: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1127: {
1128: char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p = NULL, *q = NULL, *r = NULL;
1129: int snum, region, dim, tag;
1131: PetscFunctionBegin;
1132: PetscCall(GmshReadString(gmsh, line, 1));
1133: snum = sscanf(line, "%d", ®ion);
1134: mesh->numRegions = region;
1135: PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1136: PetscCall(PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames));
1137: for (region = 0; region < mesh->numRegions; ++region) {
1138: PetscCall(GmshReadString(gmsh, line, 2));
1139: snum = sscanf(line, "%d %d", &dim, &tag);
1140: PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1141: PetscCall(GmshReadString(gmsh, line, -(PetscInt)sizeof(line)));
1142: PetscCall(PetscStrchr(line, '"', &p));
1143: PetscCheck(p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1144: PetscCall(PetscStrrchr(line, '"', &q));
1145: PetscCheck(q != p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1146: PetscCall(PetscStrrchr(line, ':', &r));
1147: if (p != r) q = r;
1148: PetscCall(PetscStrncpy(name, p + 1, (size_t)(q - p - 1)));
1149: mesh->regionDims[region] = dim;
1150: mesh->regionTags[region] = tag;
1151: PetscCall(PetscStrallocpy(name, &mesh->regionNames[region]));
1152: }
1153: PetscFunctionReturn(PETSC_SUCCESS);
1154: }
1156: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1157: {
1158: PetscFunctionBegin;
1159: switch (gmsh->fileFormat) {
1160: case 41:
1161: PetscCall(GmshReadEntities_v41(gmsh, mesh));
1162: break;
1163: default:
1164: PetscCall(GmshReadEntities_v40(gmsh, mesh));
1165: break;
1166: }
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1171: {
1172: PetscFunctionBegin;
1173: switch (gmsh->fileFormat) {
1174: case 41:
1175: PetscCall(GmshReadNodes_v41(gmsh, mesh));
1176: break;
1177: case 40:
1178: PetscCall(GmshReadNodes_v40(gmsh, mesh));
1179: break;
1180: default:
1181: PetscCall(GmshReadNodes_v22(gmsh, mesh));
1182: break;
1183: }
1185: { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1186: if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1187: PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1188: GmshNodes *nodes = mesh->nodelist;
1189: for (n = 0; n < mesh->numNodes; ++n) {
1190: const PetscInt tag = nodes->id[n];
1191: tagMin = PetscMin(tag, tagMin);
1192: tagMax = PetscMax(tag, tagMax);
1193: }
1194: gmsh->nodeStart = tagMin;
1195: gmsh->nodeEnd = tagMax + 1;
1196: }
1197: }
1199: { /* Support for sparse node tags */
1200: PetscInt n, t;
1201: GmshNodes *nodes = mesh->nodelist;
1202: PetscCall(PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf));
1203: for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1204: gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1205: for (n = 0; n < mesh->numNodes; ++n) {
1206: const PetscInt tag = nodes->id[n];
1207: PetscCheck(gmsh->nodeMap[tag] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %" PetscInt_FMT, tag);
1208: gmsh->nodeMap[tag] = n;
1209: }
1210: }
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1215: {
1216: PetscFunctionBegin;
1217: switch (gmsh->fileFormat) {
1218: case 41:
1219: PetscCall(GmshReadElements_v41(gmsh, mesh));
1220: break;
1221: case 40:
1222: PetscCall(GmshReadElements_v40(gmsh, mesh));
1223: break;
1224: default:
1225: PetscCall(GmshReadElements_v22(gmsh, mesh));
1226: break;
1227: }
1229: { /* Reorder elements by codimension and polytope type */
1230: PetscInt ne = mesh->numElems;
1231: GmshElement *elements = mesh->elements;
1232: PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0;
1233: PetscInt offset[GMSH_NUM_POLYTOPES + 1], e, k;
1235: for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1236: PetscCall(PetscMemzero(offset, sizeof(offset)));
1238: keymap[GMSH_TET] = nk++;
1239: keymap[GMSH_HEX] = nk++;
1240: keymap[GMSH_PRI] = nk++;
1241: keymap[GMSH_PYR] = nk++;
1242: keymap[GMSH_TRI] = nk++;
1243: keymap[GMSH_QUA] = nk++;
1244: keymap[GMSH_SEG] = nk++;
1245: keymap[GMSH_VTX] = nk++;
1247: PetscCall(GmshElementsCreate(mesh->numElems, &mesh->elements));
1248: #define key(eid) keymap[GmshCellMap[elements[eid].cellType].polytope]
1249: for (e = 0; e < ne; ++e) offset[1 + key(e)]++;
1250: for (k = 1; k < nk; ++k) offset[k] += offset[k - 1];
1251: for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1252: #undef key
1253: PetscCall(GmshElementsDestroy(&elements));
1254: }
1256: { /* Mesh dimension and order */
1257: GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1258: mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0;
1259: mesh->order = elem ? GmshCellMap[elem->cellType].order : 0;
1260: }
1262: {
1263: PetscBT vtx;
1264: PetscInt dim = mesh->dim, e, n, v;
1266: PetscCall(PetscBTCreate(mesh->numNodes, &vtx));
1268: /* Compute number of cells and set of vertices */
1269: mesh->numCells = 0;
1270: for (e = 0; e < mesh->numElems; ++e) {
1271: GmshElement *elem = mesh->elements + e;
1272: if (elem->dim == dim && dim > 0) mesh->numCells++;
1273: for (v = 0; v < elem->numVerts; v++) PetscCall(PetscBTSet(vtx, elem->nodes[v]));
1274: }
1276: /* Compute numbering for vertices */
1277: mesh->numVerts = 0;
1278: PetscCall(PetscMalloc1(mesh->numNodes, &mesh->vertexMap));
1279: for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;
1281: PetscCall(PetscBTDestroy(&vtx));
1282: }
1283: PetscFunctionReturn(PETSC_SUCCESS);
1284: }
1286: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1287: {
1288: PetscInt n;
1290: PetscFunctionBegin;
1291: PetscCall(PetscMalloc1(mesh->numNodes, &mesh->periodMap));
1292: for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1293: switch (gmsh->fileFormat) {
1294: case 41:
1295: PetscCall(GmshReadPeriodic_v41(gmsh, mesh->periodMap));
1296: break;
1297: default:
1298: PetscCall(GmshReadPeriodic_v40(gmsh, mesh->periodMap));
1299: break;
1300: }
1302: /* Find canonical primary nodes */
1303: for (n = 0; n < mesh->numNodes; ++n)
1304: while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];
1306: /* Renumber vertices (filter out correspondings) */
1307: mesh->numVerts = 0;
1308: for (n = 0; n < mesh->numNodes; ++n)
1309: if (mesh->vertexMap[n] >= 0) /* is vertex */
1310: if (mesh->periodMap[n] == n) /* is primary */
1311: mesh->vertexMap[n] = mesh->numVerts++;
1312: for (n = 0; n < mesh->numNodes; ++n)
1313: if (mesh->vertexMap[n] >= 0) /* is vertex */
1314: if (mesh->periodMap[n] != n) /* is corresponding */
1315: mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1316: PetscFunctionReturn(PETSC_SUCCESS);
1317: }
1319: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1320: static const DMPolytopeType DMPolytopeMap[] = {
1321: /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1322: /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1323: /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1324: /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1325: /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1326: /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1327: /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1328: /* GMSH_PYR */ DM_POLYTOPE_PYRAMID, DM_POLYTOPE_UNKNOWN};
1330: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1331: {
1332: return DMPolytopeMap[GmshCellMap[cellType].polytope];
1333: }
1335: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1336: {
1337: DM K;
1338: PetscSpace P;
1339: PetscDualSpace Q;
1340: PetscQuadrature q, fq;
1341: PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1342: PetscBool endpoint = PETSC_TRUE;
1343: char name[32];
1345: PetscFunctionBegin;
1346: /* Create space */
1347: PetscCall(PetscSpaceCreate(comm, &P));
1348: PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
1349: PetscCall(PetscSpacePolynomialSetTensor(P, isTensor));
1350: PetscCall(PetscSpaceSetNumComponents(P, Nc));
1351: PetscCall(PetscSpaceSetNumVariables(P, dim));
1352: PetscCall(PetscSpaceSetDegree(P, k, PETSC_DETERMINE));
1353: if (prefix) {
1354: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
1355: PetscCall(PetscSpaceSetFromOptions(P));
1356: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, NULL));
1357: PetscCall(PetscSpaceGetDegree(P, &k, NULL));
1358: }
1359: PetscCall(PetscSpaceSetUp(P));
1360: /* Create dual space */
1361: PetscCall(PetscDualSpaceCreate(comm, &Q));
1362: PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
1363: PetscCall(PetscDualSpaceLagrangeSetTensor(Q, isTensor));
1364: PetscCall(PetscDualSpaceLagrangeSetContinuity(Q, continuity));
1365: PetscCall(PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0));
1366: PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
1367: PetscCall(PetscDualSpaceSetOrder(Q, k));
1368: PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K));
1369: PetscCall(PetscDualSpaceSetDM(Q, K));
1370: PetscCall(DMDestroy(&K));
1371: if (prefix) {
1372: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
1373: PetscCall(PetscDualSpaceSetFromOptions(Q));
1374: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, NULL));
1375: }
1376: PetscCall(PetscDualSpaceSetUp(Q));
1377: /* Create quadrature */
1378: if (isSimplex) {
1379: PetscCall(PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q));
1380: PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1381: } else {
1382: PetscCall(PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q));
1383: PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1384: }
1385: /* Create finite element */
1386: PetscCall(PetscFECreate(comm, fem));
1387: PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
1388: PetscCall(PetscFESetNumComponents(*fem, Nc));
1389: PetscCall(PetscFESetBasisSpace(*fem, P));
1390: PetscCall(PetscFESetDualSpace(*fem, Q));
1391: PetscCall(PetscFESetQuadrature(*fem, q));
1392: PetscCall(PetscFESetFaceQuadrature(*fem, fq));
1393: if (prefix) {
1394: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
1395: PetscCall(PetscFESetFromOptions(*fem));
1396: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL));
1397: }
1398: PetscCall(PetscFESetUp(*fem));
1399: /* Cleanup */
1400: PetscCall(PetscSpaceDestroy(&P));
1401: PetscCall(PetscDualSpaceDestroy(&Q));
1402: PetscCall(PetscQuadratureDestroy(&q));
1403: PetscCall(PetscQuadratureDestroy(&fq));
1404: /* Set finite element name */
1405: PetscCall(PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k));
1406: PetscCall(PetscFESetName(*fem, name));
1407: PetscFunctionReturn(PETSC_SUCCESS);
1408: }
1410: /*@C
1411: DMPlexCreateGmshFromFile - Create a `DMPLEX` mesh from a Gmsh file
1413: Input Parameters:
1414: + comm - The MPI communicator
1415: . filename - Name of the Gmsh file
1416: - interpolate - Create faces and edges in the mesh
1418: Output Parameter:
1419: . dm - The `DM` object representing the mesh
1421: Level: beginner
1423: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
1424: @*/
1425: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1426: {
1427: PetscViewer viewer;
1428: PetscMPIInt rank;
1429: int fileType;
1430: PetscViewerType vtype;
1432: PetscFunctionBegin;
1433: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1435: /* Determine Gmsh file type (ASCII or binary) from file header */
1436: if (rank == 0) {
1437: GmshFile gmsh[1];
1438: char line[PETSC_MAX_PATH_LEN];
1439: int snum;
1440: float version;
1441: int fileFormat;
1443: PetscCall(PetscArrayzero(gmsh, 1));
1444: PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer));
1445: PetscCall(PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII));
1446: PetscCall(PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ));
1447: PetscCall(PetscViewerFileSetName(gmsh->viewer, filename));
1448: /* Read only the first two lines of the Gmsh file */
1449: PetscCall(GmshReadSection(gmsh, line));
1450: PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1451: PetscCall(GmshReadString(gmsh, line, 2));
1452: snum = sscanf(line, "%f %d", &version, &fileType);
1453: fileFormat = (int)roundf(version * 10);
1454: PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1455: PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1456: PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1457: PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1458: PetscCall(PetscViewerDestroy(&gmsh->viewer));
1459: }
1460: PetscCallMPI(MPI_Bcast(&fileType, 1, MPI_INT, 0, comm));
1461: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1463: /* Create appropriate viewer and build plex */
1464: PetscCall(PetscViewerCreate(comm, &viewer));
1465: PetscCall(PetscViewerSetType(viewer, vtype));
1466: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
1467: PetscCall(PetscViewerFileSetName(viewer, filename));
1468: PetscCall(DMPlexCreateGmsh(comm, viewer, interpolate, dm));
1469: PetscCall(PetscViewerDestroy(&viewer));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /*@
1474: DMPlexCreateGmsh - Create a `DMPLEX` mesh from a Gmsh file viewer
1476: Collective
1478: Input Parameters:
1479: + comm - The MPI communicator
1480: . viewer - The `PetscViewer` associated with a Gmsh file
1481: - interpolate - Create faces and edges in the mesh
1483: Output Parameter:
1484: . dm - The `DM` object representing the mesh
1486: Options Database Keys:
1487: + -dm_plex_gmsh_hybrid - Force triangular prisms to use tensor order
1488: . -dm_plex_gmsh_periodic - Read Gmsh periodic section and construct a periodic Plex
1489: . -dm_plex_gmsh_highorder - Generate high-order coordinates
1490: . -dm_plex_gmsh_project - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
1491: . -dm_plex_gmsh_use_generic - Generate generic labels, i.e. Cell Sets, Face Sets, etc.
1492: . -dm_plex_gmsh_use_regions - Generate labels with region names
1493: . -dm_plex_gmsh_mark_vertices - Add vertices to generated labels
1494: . -dm_plex_gmsh_multiple_tags - Allow multiple tags for default labels
1495: - -dm_plex_gmsh_spacedim <d> - Embedding space dimension, if different from topological dimension
1497: Level: beginner
1499: Notes:
1500: The Gmsh file format is described in <http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format>
1502: By default, the "Cell Sets", "Face Sets", and "Vertex Sets" labels are created, and only insert the first tag on a point. By using -dm_plex_gmsh_multiple_tags, all tags can be inserted. Alternatively, -dm_plex_gmsh_use_regions creates labels based on the region names from the PhysicalNames section, and all tags are used, but you can retain the generic labels using -dm_plex_gmsh_use_generic.
1504: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1505: @*/
1506: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1507: {
1508: GmshMesh *mesh = NULL;
1509: PetscViewer parentviewer = NULL;
1510: PetscBT periodicVerts = NULL;
1511: PetscBT *periodicCells = NULL;
1512: DM cdm, cdmCell = NULL;
1513: PetscSection cs, csCell = NULL;
1514: Vec coordinates, coordinatesCell;
1515: DMLabel cellSets = NULL, faceSets = NULL, edgeSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1516: PetscInt dim = 0, coordDim = -1, order = 0, maxHeight = 0;
1517: PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0, vStart, vEnd;
1518: PetscInt cell, cone[8], e, n, v, d;
1519: PetscBool usegeneric = PETSC_TRUE, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE;
1520: PetscBool flg, binary, hybrid = interpolate, periodic = PETSC_TRUE;
1521: PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1522: PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1523: PetscMPIInt rank;
1525: PetscFunctionBegin;
1526: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1527: PetscObjectOptionsBegin((PetscObject)viewer);
1528: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options");
1529: PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL));
1530: PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL));
1531: PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet));
1532: PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL));
1533: PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL));
1534: PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_generic", "Generate generic labels, i.e. Cell Sets, Face Sets, etc", "DMPlexCreateGmsh", usegeneric, &usegeneric, &flg));
1535: if (!flg && useregions) usegeneric = PETSC_FALSE;
1536: PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL));
1537: PetscCall(PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL));
1538: PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE));
1539: PetscCall(PetscOptionsBoundedInt("-dm_localize_height", "Localize edges and faces in addition to cells", "", maxHeight, &maxHeight, NULL, 0));
1540: PetscOptionsHeadEnd();
1541: PetscOptionsEnd();
1543: PetscCall(GmshCellInfoSetUp());
1545: PetscCall(DMCreate(comm, dm));
1546: PetscCall(DMSetType(*dm, DMPLEX));
1547: PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));
1549: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary));
1551: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1552: if (binary) {
1553: parentviewer = viewer;
1554: PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer));
1555: }
1557: if (rank == 0) {
1558: GmshFile gmsh[1];
1559: char line[PETSC_MAX_PATH_LEN];
1560: PetscBool match;
1562: PetscCall(PetscArrayzero(gmsh, 1));
1563: gmsh->viewer = viewer;
1564: gmsh->binary = binary;
1566: PetscCall(GmshMeshCreate(&mesh));
1568: /* Read mesh format */
1569: PetscCall(GmshReadSection(gmsh, line));
1570: PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1571: PetscCall(GmshReadMeshFormat(gmsh));
1572: PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line));
1574: /* OPTIONAL Read mesh version (Neper only) */
1575: PetscCall(GmshReadSection(gmsh, line));
1576: PetscCall(GmshMatch(gmsh, "$MeshVersion", line, &match));
1577: if (match) {
1578: PetscCall(GmshExpect(gmsh, "$MeshVersion", line));
1579: PetscCall(GmshReadMeshVersion(gmsh));
1580: PetscCall(GmshReadEndSection(gmsh, "$EndMeshVersion", line));
1581: /* Initial read for entity section */
1582: PetscCall(GmshReadSection(gmsh, line));
1583: }
1585: /* OPTIONAL Read mesh domain (Neper only) */
1586: PetscCall(GmshMatch(gmsh, "$Domain", line, &match));
1587: if (match) {
1588: PetscCall(GmshExpect(gmsh, "$Domain", line));
1589: PetscCall(GmshReadMeshDomain(gmsh));
1590: PetscCall(GmshReadEndSection(gmsh, "$EndDomain", line));
1591: /* Initial read for entity section */
1592: PetscCall(GmshReadSection(gmsh, line));
1593: }
1595: /* OPTIONAL Read physical names */
1596: PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match));
1597: if (match) {
1598: PetscCall(GmshExpect(gmsh, "$PhysicalNames", line));
1599: PetscCall(GmshReadPhysicalNames(gmsh, mesh));
1600: PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line));
1601: /* Initial read for entity section */
1602: PetscCall(GmshReadSection(gmsh, line));
1603: }
1605: /* Read entities */
1606: if (gmsh->fileFormat >= 40) {
1607: PetscCall(GmshExpect(gmsh, "$Entities", line));
1608: PetscCall(GmshReadEntities(gmsh, mesh));
1609: PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line));
1610: /* Initial read for nodes section */
1611: PetscCall(GmshReadSection(gmsh, line));
1612: }
1614: /* Read nodes */
1615: PetscCall(GmshExpect(gmsh, "$Nodes", line));
1616: PetscCall(GmshReadNodes(gmsh, mesh));
1617: PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line));
1619: /* Read elements */
1620: PetscCall(GmshReadSection(gmsh, line));
1621: PetscCall(GmshExpect(gmsh, "$Elements", line));
1622: PetscCall(GmshReadElements(gmsh, mesh));
1623: PetscCall(GmshReadEndSection(gmsh, "$EndElements", line));
1625: /* Read periodic section (OPTIONAL) */
1626: if (periodic) {
1627: PetscCall(GmshReadSection(gmsh, line));
1628: PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic));
1629: }
1630: if (periodic) {
1631: PetscCall(GmshExpect(gmsh, "$Periodic", line));
1632: PetscCall(GmshReadPeriodic(gmsh, mesh));
1633: PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line));
1634: }
1636: PetscCall(PetscFree(gmsh->wbuf));
1637: PetscCall(PetscFree(gmsh->sbuf));
1638: PetscCall(PetscFree(gmsh->nbuf));
1640: dim = mesh->dim;
1641: order = mesh->order;
1642: numNodes = mesh->numNodes;
1643: numElems = mesh->numElems;
1644: numVerts = mesh->numVerts;
1645: numCells = mesh->numCells;
1647: {
1648: GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1649: GmshElement *elemB = PetscSafePointerPlusOffset(elemA, mesh->numCells - 1);
1650: int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1651: int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1652: isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1653: isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1654: hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1655: }
1656: }
1658: if (parentviewer) PetscCall(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer));
1660: {
1661: int buf[6];
1663: buf[0] = (int)dim;
1664: buf[1] = (int)order;
1665: buf[2] = periodic;
1666: buf[3] = isSimplex;
1667: buf[4] = isHybrid;
1668: buf[5] = hasTetra;
1670: PetscCallMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm));
1672: dim = buf[0];
1673: order = buf[1];
1674: periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1675: isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1676: isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1677: hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1678: }
1680: if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1681: PetscCheck(!highOrder || !isHybrid, comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet");
1683: /* We do not want this label automatically computed, instead we fill it here */
1684: PetscCall(DMCreateLabel(*dm, "celltype"));
1686: /* Allocate the cell-vertex mesh */
1687: PetscCall(DMPlexSetChart(*dm, 0, numCells + numVerts));
1688: for (cell = 0; cell < numCells; ++cell) {
1689: GmshElement *elem = mesh->elements + cell;
1690: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1691: if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1692: PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts));
1693: PetscCall(DMPlexSetCellType(*dm, cell, ctype));
1694: }
1695: for (v = numCells; v < numCells + numVerts; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1696: PetscCall(DMSetUp(*dm));
1698: /* Add cell-vertex connections */
1699: for (cell = 0; cell < numCells; ++cell) {
1700: GmshElement *elem = mesh->elements + cell;
1701: for (v = 0; v < elem->numVerts; ++v) {
1702: const PetscInt nn = elem->nodes[v];
1703: const PetscInt vv = mesh->vertexMap[nn];
1704: cone[v] = numCells + vv;
1705: }
1706: PetscCall(DMPlexReorderCell(*dm, cell, cone));
1707: PetscCall(DMPlexSetCone(*dm, cell, cone));
1708: }
1710: PetscCall(DMSetDimension(*dm, dim));
1711: PetscCall(DMPlexSymmetrize(*dm));
1712: PetscCall(DMPlexStratify(*dm));
1713: if (interpolate) {
1714: DM idm;
1716: PetscCall(DMPlexInterpolate(*dm, &idm));
1717: PetscCall(DMDestroy(dm));
1718: *dm = idm;
1719: }
1720: PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd));
1722: if (rank == 0) {
1723: const PetscInt Nr = useregions ? mesh->numRegions : 0;
1725: PetscCall(PetscCalloc1(Nr, ®ionSets));
1726: for (cell = 0, e = 0; e < numElems; ++e) {
1727: GmshElement *elem = mesh->elements + e;
1729: /* Create cell sets */
1730: if (elem->dim == dim && dim > 0) {
1731: if (elem->numTags > 0) {
1732: PetscInt Nt = elem->numTags, t, r;
1734: for (t = 0; t < Nt; ++t) {
1735: const PetscInt tag = elem->tags[t];
1736: const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1738: if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag));
1739: for (r = 0; r < Nr; ++r) {
1740: if (mesh->regionDims[r] != dim) continue;
1741: if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag));
1742: }
1743: }
1744: }
1745: cell++;
1746: }
1748: /* Create face sets */
1749: if (elem->numTags && interpolate && elem->dim == dim - 1) {
1750: PetscInt joinSize;
1751: const PetscInt *join = NULL;
1752: PetscInt Nt = elem->numTags, pdepth;
1754: /* Find the relevant facet with vertex joins */
1755: for (v = 0; v < elem->numVerts; ++v) {
1756: const PetscInt nn = elem->nodes[v];
1757: const PetscInt vv = mesh->vertexMap[nn];
1758: cone[v] = vStart + vv;
1759: }
1760: PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1761: PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex facet for Gmsh element %" PetscInt_FMT " (Plex cell %" PetscInt_FMT ")", elem->id, e);
1762: PetscCall(DMPlexGetPointDepth(*dm, join[0], &pdepth));
1763: PetscCheck(pdepth == dim - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Plex facet %" PetscInt_FMT " for Gmsh element %" PetscInt_FMT " had depth %" PetscInt_FMT " != %" PetscInt_FMT, join[0], elem->id, pdepth, dim - 1);
1764: for (PetscInt t = 0; t < Nt; ++t) {
1765: const PetscInt tag = elem->tags[t];
1766: const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1768: if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag));
1769: for (PetscInt r = 0; r < Nr; ++r) {
1770: if (mesh->regionDims[r] != dim - 1) continue;
1771: if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag));
1772: }
1773: }
1774: PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1775: }
1777: /* Create edge sets */
1778: if (elem->numTags && interpolate && dim > 2 && elem->dim == 1) {
1779: PetscInt joinSize;
1780: const PetscInt *join = NULL;
1781: PetscInt Nt = elem->numTags, t, r;
1783: /* Find the relevant edge with vertex joins */
1784: for (v = 0; v < elem->numVerts; ++v) {
1785: const PetscInt nn = elem->nodes[v];
1786: const PetscInt vv = mesh->vertexMap[nn];
1787: cone[v] = vStart + vv;
1788: }
1789: PetscCall(DMPlexGetJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1790: PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex edge for Gmsh element %" PetscInt_FMT " (Plex cell %" PetscInt_FMT ")", elem->id, e);
1791: for (t = 0; t < Nt; ++t) {
1792: const PetscInt tag = elem->tags[t];
1793: const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1795: if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &edgeSets, "Edge Sets", join[0], tag));
1796: for (r = 0; r < Nr; ++r) {
1797: if (mesh->regionDims[r] != 1) continue;
1798: if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag));
1799: }
1800: }
1801: PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1802: }
1804: /* Create vertex sets */
1805: if (elem->numTags && elem->dim == 0 && markvertices) {
1806: const PetscInt nn = elem->nodes[0];
1807: const PetscInt vv = mesh->vertexMap[nn];
1808: const PetscInt tag = elem->tags[0];
1809: PetscInt r;
1811: if (vv < 0) continue;
1812: if (usegeneric) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1813: for (r = 0; r < Nr; ++r) {
1814: if (mesh->regionDims[r] != 0) continue;
1815: if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag));
1816: }
1817: }
1818: }
1819: if (markvertices) {
1820: for (v = 0; v < numNodes; ++v) {
1821: const PetscInt vv = mesh->vertexMap[v];
1822: const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS];
1823: PetscInt r, t;
1825: if (vv < 0) continue;
1826: for (t = 0; t < GMSH_MAX_TAGS; ++t) {
1827: const PetscInt tag = tags[t];
1828: const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1830: if (tag == -1) continue;
1831: if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1832: for (r = 0; r < Nr; ++r) {
1833: if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag));
1834: }
1835: }
1836: }
1837: }
1838: PetscCall(PetscFree(regionSets));
1839: }
1841: { /* Create Cell/Face/Edge/Vertex Sets labels at all processes */
1842: enum {
1843: n = 5
1844: };
1845: PetscBool flag[n];
1847: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1848: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1849: flag[2] = edgeSets ? PETSC_TRUE : PETSC_FALSE;
1850: flag[3] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1851: flag[4] = marker ? PETSC_TRUE : PETSC_FALSE;
1852: PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
1853: if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1854: if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1855: if (flag[2]) PetscCall(DMCreateLabel(*dm, "Edge Sets"));
1856: if (flag[3]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1857: if (flag[4]) PetscCall(DMCreateLabel(*dm, "marker"));
1858: }
1860: if (periodic) {
1861: PetscCall(PetscBTCreate(numVerts, &periodicVerts));
1862: for (n = 0; n < numNodes; ++n) {
1863: if (mesh->vertexMap[n] >= 0) {
1864: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1865: PetscInt m = mesh->periodMap[n];
1866: PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n]));
1867: PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m]));
1868: }
1869: }
1870: }
1871: PetscCall(DMGetCoordinateDM(*dm, &cdm));
1872: PetscCall(PetscMalloc1(maxHeight + 1, &periodicCells));
1873: for (PetscInt h = 0; h <= maxHeight; ++h) {
1874: PetscInt pStart, pEnd;
1876: PetscCall(DMPlexGetHeightStratum(*dm, h, &pStart, &pEnd));
1877: PetscCall(PetscBTCreate(pEnd - pStart, &periodicCells[h]));
1878: for (PetscInt p = pStart; p < pEnd; ++p) {
1879: PetscInt *closure = NULL;
1880: PetscInt Ncl;
1882: PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1883: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1884: if (closure[cl] >= vStart && closure[cl] < vEnd) {
1885: if (PetscUnlikely(PetscBTLookup(periodicVerts, closure[cl] - vStart))) {
1886: PetscCall(PetscBTSet(periodicCells[h], p - pStart));
1887: break;
1888: }
1889: }
1890: }
1891: PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1892: }
1893: }
1894: }
1896: /* Setup coordinate DM */
1897: if (coordDim < 0) coordDim = dim;
1898: PetscCall(DMSetCoordinateDim(*dm, coordDim));
1899: PetscCall(DMGetCoordinateDM(*dm, &cdm));
1900: if (highOrder) {
1901: PetscFE fe;
1902: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1903: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1905: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1907: PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
1908: PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view"));
1909: PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe));
1910: PetscCall(PetscFEDestroy(&fe));
1911: PetscCall(DMCreateDS(cdm));
1912: }
1913: if (periodic) {
1914: PetscCall(DMClone(cdm, &cdmCell));
1915: PetscCall(DMSetCellCoordinateDM(*dm, cdmCell));
1916: }
1918: /* Create coordinates */
1919: if (highOrder) {
1920: PetscInt maxDof = GmshNumNodes_HEX(order) * coordDim;
1921: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1922: PetscSection section;
1923: PetscScalar *cellCoords;
1925: PetscCall(DMSetLocalSection(cdm, NULL));
1926: PetscCall(DMGetLocalSection(cdm, &cs));
1927: PetscCall(PetscSectionClone(cs, §ion));
1928: PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1930: PetscCall(DMCreateLocalVector(cdm, &coordinates));
1931: PetscCall(PetscMalloc1(maxDof, &cellCoords));
1932: for (cell = 0; cell < numCells; ++cell) {
1933: GmshElement *elem = mesh->elements + cell;
1934: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1935: int s = 0;
1936: for (n = 0; n < elem->numNodes; ++n) {
1937: while (lexorder[n + s] < 0) ++s;
1938: const PetscInt node = elem->nodes[lexorder[n + s]];
1939: for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d];
1940: }
1941: if (s) {
1942: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1943: PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1944: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1945: PetscReal hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1946: PetscReal hexFrontWeights[27] = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1947: PetscReal hexLeftWeights[27] = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0};
1948: PetscReal hexRightWeights[27] = {0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25};
1949: PetscReal hexBackWeights[27] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25};
1950: PetscReal hexTopWeights[27] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1951: PetscReal hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25};
1952: PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
1953: PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
1954: NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL};
1955: PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3};
1957: /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
1958: for (n = 0; n < elem->numNodes + s; ++n) {
1959: if (lexorder[n] >= 0) continue;
1960: for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0;
1961: for (int bn = 0; bn < elem->numNodes + s; ++bn) {
1962: if (lexorder[bn] < 0) continue;
1963: const PetscReal *weights = sdWeights[coordDim][n];
1964: const PetscInt bnode = elem->nodes[lexorder[bn]];
1965: for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d];
1966: }
1967: }
1968: }
1969: PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES));
1970: }
1971: PetscCall(PetscSectionDestroy(§ion));
1972: PetscCall(PetscFree(cellCoords));
1973: } else {
1974: PetscInt *nodeMap;
1975: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1976: PetscScalar *pointCoords;
1978: PetscCall(DMGetCoordinateSection(*dm, &cs));
1979: PetscCall(PetscSectionSetNumFields(cs, 1));
1980: PetscCall(PetscSectionSetFieldComponents(cs, 0, coordDim));
1981: PetscCall(PetscSectionSetChart(cs, numCells, numCells + numVerts));
1982: for (v = numCells; v < numCells + numVerts; ++v) {
1983: PetscCall(PetscSectionSetDof(cs, v, coordDim));
1984: PetscCall(PetscSectionSetFieldDof(cs, v, 0, coordDim));
1985: }
1986: PetscCall(PetscSectionSetUp(cs));
1988: /* We need to localize coordinates on cells */
1989: if (periodic) {
1990: PetscInt newStart = PETSC_MAX_INT, newEnd = -1, pStart, pEnd;
1992: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell));
1993: PetscCall(PetscSectionSetNumFields(csCell, 1));
1994: PetscCall(PetscSectionSetFieldComponents(csCell, 0, coordDim));
1995: for (PetscInt h = 0; h <= maxHeight; h++) {
1996: PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
1997: newStart = PetscMin(newStart, pStart);
1998: newEnd = PetscMax(newEnd, pEnd);
1999: }
2000: PetscCall(PetscSectionSetChart(csCell, newStart, newEnd));
2001: for (PetscInt h = 0; h <= maxHeight; h++) {
2002: PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
2003: for (PetscInt p = pStart; p < pEnd; ++p) {
2004: PetscInt *closure = NULL;
2005: PetscInt Ncl, Nv = 0;
2007: if (PetscUnlikely(PetscBTLookup(periodicCells[h], p - pStart))) {
2008: PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
2009: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2010: if (closure[cl] >= vStart && closure[cl] < vEnd) ++Nv;
2011: }
2012: PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
2013: PetscCall(PetscSectionSetDof(csCell, p, Nv * coordDim));
2014: PetscCall(PetscSectionSetFieldDof(csCell, p, 0, Nv * coordDim));
2015: }
2016: }
2017: }
2018: PetscCall(PetscSectionSetUp(csCell));
2019: PetscCall(DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell));
2020: }
2022: PetscCall(DMCreateLocalVector(cdm, &coordinates));
2023: PetscCall(VecGetArray(coordinates, &pointCoords));
2024: PetscCall(PetscMalloc1(numVerts, &nodeMap));
2025: for (n = 0; n < numNodes; n++)
2026: if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n;
2027: for (v = 0; v < numVerts; ++v) {
2028: PetscInt off, node = nodeMap[v];
2030: PetscCall(PetscSectionGetOffset(cs, numCells + v, &off));
2031: for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
2032: }
2033: PetscCall(VecRestoreArray(coordinates, &pointCoords));
2035: if (periodic) {
2036: PetscInt cStart, cEnd;
2038: PetscCall(DMPlexGetHeightStratum(cdmCell, 0, &cStart, &cEnd));
2039: PetscCall(DMCreateLocalVector(cdmCell, &coordinatesCell));
2040: PetscCall(VecGetArray(coordinatesCell, &pointCoords));
2041: for (PetscInt c = cStart; c < cEnd; ++c) {
2042: GmshElement *elem = mesh->elements + c - cStart;
2043: PetscInt *closure = NULL;
2044: PetscInt verts[8];
2045: PetscInt Ncl, Nv = 0;
2047: for (PetscInt v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v];
2048: PetscCall(DMPlexReorderCell(cdmCell, c, cone));
2049: PetscCall(DMPlexGetTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
2050: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2051: if (closure[cl] >= vStart && closure[cl] < vEnd) verts[Nv++] = closure[cl];
2052: }
2053: PetscCheck(Nv == elem->numVerts, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of vertices %" PetscInt_FMT " in closure does not match number of vertices %" PetscInt_FMT " in Gmsh cell", Nv, elem->numVerts);
2054: for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2055: const PetscInt point = closure[cl];
2056: PetscInt hStart, h;
2058: PetscCall(DMPlexGetPointHeight(cdmCell, point, &h));
2059: if (h > maxHeight) continue;
2060: PetscCall(DMPlexGetHeightStratum(cdmCell, h, &hStart, NULL));
2061: if (PetscUnlikely(PetscBTLookup(periodicCells[h], point - hStart))) {
2062: PetscInt *pclosure = NULL;
2063: PetscInt Npcl, off, v;
2065: PetscCall(PetscSectionGetOffset(csCell, point, &off));
2066: PetscCall(DMPlexGetTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2067: for (PetscInt pcl = 0; pcl < Npcl * 2; pcl += 2) {
2068: if (pclosure[pcl] >= vStart && pclosure[pcl] < vEnd) {
2069: for (v = 0; v < Nv; ++v)
2070: if (verts[v] == pclosure[pcl]) break;
2071: PetscCheck(v < Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find vertex %" PetscInt_FMT " in closure of cell %" PetscInt_FMT, pclosure[pcl], c);
2072: for (PetscInt d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[cone[v] * 3 + d];
2073: ++v;
2074: }
2075: }
2076: PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2077: }
2078: }
2079: PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
2080: }
2081: PetscCall(VecSetBlockSize(coordinatesCell, coordDim));
2082: PetscCall(VecRestoreArray(coordinatesCell, &pointCoords));
2083: PetscCall(DMSetCellCoordinatesLocal(*dm, coordinatesCell));
2084: PetscCall(VecDestroy(&coordinatesCell));
2085: }
2086: PetscCall(PetscFree(nodeMap));
2087: PetscCall(PetscSectionDestroy(&csCell));
2088: PetscCall(DMDestroy(&cdmCell));
2089: }
2091: PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
2092: PetscCall(VecSetBlockSize(coordinates, coordDim));
2093: PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
2094: PetscCall(VecDestroy(&coordinates));
2096: PetscCall(GmshMeshDestroy(&mesh));
2097: PetscCall(PetscBTDestroy(&periodicVerts));
2098: if (periodic) {
2099: for (PetscInt h = 0; h <= maxHeight; ++h) PetscCall(PetscBTDestroy(periodicCells + h));
2100: PetscCall(PetscFree(periodicCells));
2101: }
2103: if (highOrder && project) {
2104: PetscFE fe;
2105: const char prefix[] = "dm_plex_gmsh_project_";
2106: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
2107: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
2109: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
2110: PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
2111: PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view"));
2112: PetscCall(DMSetCoordinateDisc(*dm, fe, PETSC_TRUE));
2113: PetscCall(PetscFEDestroy(&fe));
2114: }
2116: PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));
2117: PetscFunctionReturn(PETSC_SUCCESS);
2118: }