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