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