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