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 0;
178: called = PETSC_TRUE;
179: n = PETSC_STATIC_ARRAY_LENGTH(GmshCellMap);
180: for (i = 0; i < n; ++i) {
181: GmshCellMap[i].cellType = -1;
182: GmshCellMap[i].polytope = -1;
183: }
184: n = PETSC_STATIC_ARRAY_LENGTH(GmshCellTable);
185: for (i = 0; i < n; ++i) {
186: if (GmshCellTable[i].cellType <= 0) continue;
187: GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i];
188: }
189: return 0;
190: }
192: #define GmshCellTypeCheck(ct) \
196: typedef struct {
197: PetscViewer viewer;
198: int fileFormat;
199: int dataSize;
200: PetscBool binary;
201: PetscBool byteSwap;
202: size_t wlen;
203: void *wbuf;
204: size_t slen;
205: void *sbuf;
206: PetscInt *nbuf;
207: PetscInt nodeStart;
208: PetscInt nodeEnd;
209: PetscInt *nodeMap;
210: } GmshFile;
212: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
213: {
214: size_t size = count * eltsize;
216: if (gmsh->wlen < size) {
217: PetscFree(gmsh->wbuf);
218: PetscMalloc(size, &gmsh->wbuf);
219: gmsh->wlen = size;
220: }
221: *(void **)buf = size ? gmsh->wbuf : NULL;
222: return 0;
223: }
225: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
226: {
227: size_t dataSize = (size_t)gmsh->dataSize;
228: size_t size = count * dataSize;
230: if (gmsh->slen < size) {
231: PetscFree(gmsh->sbuf);
232: PetscMalloc(size, &gmsh->sbuf);
233: gmsh->slen = size;
234: }
235: *(void **)buf = size ? gmsh->sbuf : NULL;
236: return 0;
237: }
239: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
240: {
241: PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
242: if (gmsh->byteSwap) PetscByteSwap(buf, dtype, count);
243: return 0;
244: }
246: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
247: {
248: PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
249: return 0;
250: }
252: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
253: {
254: PetscStrcmp(line, Section, match);
255: return 0;
256: }
258: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
259: {
260: PetscBool match;
262: GmshMatch(gmsh, Section, line, &match);
264: return 0;
265: }
267: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
268: {
269: PetscBool match;
271: while (PETSC_TRUE) {
272: GmshReadString(gmsh, line, 1);
273: GmshMatch(gmsh, "$Comments", line, &match);
274: if (!match) break;
275: while (PETSC_TRUE) {
276: GmshReadString(gmsh, line, 1);
277: GmshMatch(gmsh, "$EndComments", line, &match);
278: if (match) break;
279: }
280: }
281: return 0;
282: }
284: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
285: {
286: GmshReadString(gmsh, line, 1);
287: GmshExpect(gmsh, EndSection, line);
288: return 0;
289: }
291: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
292: {
293: PetscInt i;
294: size_t dataSize = (size_t)gmsh->dataSize;
296: if (dataSize == sizeof(PetscInt)) {
297: GmshRead(gmsh, buf, count, PETSC_INT);
298: } else if (dataSize == sizeof(int)) {
299: int *ibuf = NULL;
300: GmshBufferSizeGet(gmsh, count, &ibuf);
301: GmshRead(gmsh, ibuf, count, PETSC_ENUM);
302: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
303: } else if (dataSize == sizeof(long)) {
304: long *ibuf = NULL;
305: GmshBufferSizeGet(gmsh, count, &ibuf);
306: GmshRead(gmsh, ibuf, count, PETSC_LONG);
307: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
308: } else if (dataSize == sizeof(PetscInt64)) {
309: PetscInt64 *ibuf = NULL;
310: GmshBufferSizeGet(gmsh, count, &ibuf);
311: GmshRead(gmsh, ibuf, count, PETSC_INT64);
312: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
313: }
314: return 0;
315: }
317: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
318: {
319: GmshRead(gmsh, buf, count, PETSC_ENUM);
320: return 0;
321: }
323: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
324: {
325: GmshRead(gmsh, buf, count, PETSC_DOUBLE);
326: return 0;
327: }
329: #define GMSH_MAX_TAGS 16
331: typedef struct {
332: PetscInt id; /* Entity ID */
333: PetscInt dim; /* Dimension */
334: double bbox[6]; /* Bounding box */
335: PetscInt numTags; /* Size of tag array */
336: int tags[GMSH_MAX_TAGS]; /* Tag array */
337: } GmshEntity;
339: typedef struct {
340: GmshEntity *entity[4];
341: PetscHMapI entityMap[4];
342: } GmshEntities;
344: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
345: {
346: PetscInt dim;
348: PetscNew(entities);
349: for (dim = 0; dim < 4; ++dim) {
350: PetscCalloc1(count[dim], &(*entities)->entity[dim]);
351: PetscHMapICreate(&(*entities)->entityMap[dim]);
352: }
353: return 0;
354: }
356: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
357: {
358: PetscInt dim;
360: if (!*entities) return 0;
361: for (dim = 0; dim < 4; ++dim) {
362: PetscFree((*entities)->entity[dim]);
363: PetscHMapIDestroy(&(*entities)->entityMap[dim]);
364: }
365: PetscFree((*entities));
366: return 0;
367: }
369: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity **entity)
370: {
371: PetscHMapISet(entities->entityMap[dim], eid, index);
372: entities->entity[dim][index].dim = dim;
373: entities->entity[dim][index].id = eid;
374: if (entity) *entity = &entities->entity[dim][index];
375: return 0;
376: }
378: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity **entity)
379: {
380: PetscInt index;
382: PetscHMapIGet(entities->entityMap[dim], eid, &index);
383: *entity = &entities->entity[dim][index];
384: return 0;
385: }
387: typedef struct {
388: PetscInt *id; /* Node IDs */
389: double *xyz; /* Coordinates */
390: PetscInt *tag; /* Physical tag */
391: } GmshNodes;
393: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
394: {
395: PetscNew(nodes);
396: PetscMalloc1(count * 1, &(*nodes)->id);
397: PetscMalloc1(count * 3, &(*nodes)->xyz);
398: PetscMalloc1(count * GMSH_MAX_TAGS, &(*nodes)->tag);
399: return 0;
400: }
402: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
403: {
404: if (!*nodes) return 0;
405: PetscFree((*nodes)->id);
406: PetscFree((*nodes)->xyz);
407: PetscFree((*nodes)->tag);
408: PetscFree((*nodes));
409: return 0;
410: }
412: typedef struct {
413: PetscInt id; /* Element ID */
414: PetscInt dim; /* Dimension */
415: PetscInt cellType; /* Cell type */
416: PetscInt numVerts; /* Size of vertex array */
417: PetscInt numNodes; /* Size of node array */
418: PetscInt *nodes; /* Vertex/Node array */
419: PetscInt numTags; /* Size of physical tag array */
420: int tags[GMSH_MAX_TAGS]; /* Physical tag array */
421: } GmshElement;
423: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
424: {
425: PetscCalloc1(count, elements);
426: return 0;
427: }
429: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
430: {
431: if (!*elements) return 0;
432: PetscFree(*elements);
433: return 0;
434: }
436: typedef struct {
437: PetscInt dim;
438: PetscInt order;
439: GmshEntities *entities;
440: PetscInt numNodes;
441: GmshNodes *nodelist;
442: PetscInt numElems;
443: GmshElement *elements;
444: PetscInt numVerts;
445: PetscInt numCells;
446: PetscInt *periodMap;
447: PetscInt *vertexMap;
448: PetscSegBuffer segbuf;
449: PetscInt numRegions;
450: PetscInt *regionDims;
451: PetscInt *regionTags;
452: char **regionNames;
453: } GmshMesh;
455: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
456: {
457: PetscNew(mesh);
458: PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);
459: return 0;
460: }
462: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
463: {
464: PetscInt r;
466: if (!*mesh) return 0;
467: GmshEntitiesDestroy(&(*mesh)->entities);
468: GmshNodesDestroy(&(*mesh)->nodelist);
469: GmshElementsDestroy(&(*mesh)->elements);
470: PetscFree((*mesh)->periodMap);
471: PetscFree((*mesh)->vertexMap);
472: PetscSegBufferDestroy(&(*mesh)->segbuf);
473: for (r = 0; r < (*mesh)->numRegions; ++r) PetscFree((*mesh)->regionNames[r]);
474: PetscFree3((*mesh)->regionDims, (*mesh)->regionTags, (*mesh)->regionNames);
475: PetscFree((*mesh));
476: return 0;
477: }
479: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
480: {
481: PetscViewer viewer = gmsh->viewer;
482: PetscBool byteSwap = gmsh->byteSwap;
483: char line[PETSC_MAX_PATH_LEN];
484: int n, t, num, nid, snum;
485: GmshNodes *nodes;
487: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
488: snum = sscanf(line, "%d", &num);
490: GmshNodesCreate(num, &nodes);
491: mesh->numNodes = num;
492: mesh->nodelist = nodes;
493: for (n = 0; n < num; ++n) {
494: double *xyz = nodes->xyz + n * 3;
495: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
496: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
497: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
498: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
499: nodes->id[n] = nid;
500: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
501: }
502: return 0;
503: }
505: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
506: file contents multiple times to figure out the true number of cells and facets
507: in the given mesh. To make this more efficient we read the file contents only
508: once and store them in memory, while determining the true number of cells. */
509: static PetscErrorCode GmshReadElements_v22(GmshFile *gmsh, GmshMesh *mesh)
510: {
511: PetscViewer viewer = gmsh->viewer;
512: PetscBool binary = gmsh->binary;
513: PetscBool byteSwap = gmsh->byteSwap;
514: char line[PETSC_MAX_PATH_LEN];
515: int i, c, p, num, ibuf[1 + 4 + 1000], snum;
516: int cellType, numElem, numVerts, numNodes, numTags;
517: GmshElement *elements;
518: PetscInt *nodeMap = gmsh->nodeMap;
520: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
521: snum = sscanf(line, "%d", &num);
523: GmshElementsCreate(num, &elements);
524: mesh->numElems = num;
525: mesh->elements = elements;
526: for (c = 0; c < num;) {
527: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
528: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);
530: cellType = binary ? ibuf[0] : ibuf[1];
531: numElem = binary ? ibuf[1] : 1;
532: numTags = ibuf[2];
534: GmshCellTypeCheck(cellType);
535: numVerts = GmshCellMap[cellType].numVerts;
536: numNodes = GmshCellMap[cellType].numNodes;
538: for (i = 0; i < numElem; ++i, ++c) {
539: GmshElement *element = elements + c;
540: const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
541: const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
542: PetscViewerRead(viewer, ibuf + off, nint, NULL, PETSC_ENUM);
543: if (byteSwap) PetscByteSwap(ibuf + off, PETSC_ENUM, nint);
544: element->id = id[0];
545: element->dim = GmshCellMap[cellType].dim;
546: element->cellType = cellType;
547: element->numVerts = numVerts;
548: element->numNodes = numNodes;
549: element->numTags = PetscMin(numTags, GMSH_MAX_TAGS);
550: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
551: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
552: for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
553: }
554: }
555: return 0;
556: }
558: /*
559: $Entities
560: numPoints(unsigned long) numCurves(unsigned long)
561: numSurfaces(unsigned long) numVolumes(unsigned long)
562: // points
563: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
564: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
565: numPhysicals(unsigned long) phyisicalTag[...](int)
566: ...
567: // curves
568: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
569: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
570: numPhysicals(unsigned long) physicalTag[...](int)
571: numBREPVert(unsigned long) tagBREPVert[...](int)
572: ...
573: // surfaces
574: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
575: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
576: numPhysicals(unsigned long) physicalTag[...](int)
577: numBREPCurve(unsigned long) tagBREPCurve[...](int)
578: ...
579: // volumes
580: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
581: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
582: numPhysicals(unsigned long) physicalTag[...](int)
583: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
584: ...
585: $EndEntities
586: */
587: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
588: {
589: PetscViewer viewer = gmsh->viewer;
590: PetscBool byteSwap = gmsh->byteSwap;
591: long index, num, lbuf[4];
592: int dim, eid, numTags, *ibuf, t;
593: PetscInt count[4], i;
594: GmshEntity *entity = NULL;
596: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
597: if (byteSwap) PetscByteSwap(lbuf, PETSC_LONG, 4);
598: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
599: GmshEntitiesCreate(count, &mesh->entities);
600: for (dim = 0; dim < 4; ++dim) {
601: for (index = 0; index < count[dim]; ++index) {
602: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
603: if (byteSwap) PetscByteSwap(&eid, PETSC_ENUM, 1);
604: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
605: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
606: if (byteSwap) PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);
607: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
608: if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
609: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
610: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
611: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
612: entity->numTags = numTags = (int)PetscMin(num, GMSH_MAX_TAGS);
613: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
614: if (dim == 0) continue;
615: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
616: if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
617: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
618: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
619: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
620: }
621: }
622: return 0;
623: }
625: /*
626: $Nodes
627: numEntityBlocks(unsigned long) numNodes(unsigned long)
628: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
629: tag(int) x(double) y(double) z(double)
630: ...
631: ...
632: $EndNodes
633: */
634: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
635: {
636: PetscViewer viewer = gmsh->viewer;
637: PetscBool byteSwap = gmsh->byteSwap;
638: long block, node, n, t, numEntityBlocks, numTotalNodes, numNodes;
639: int info[3], nid;
640: GmshNodes *nodes;
642: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
643: if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
644: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
645: if (byteSwap) PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);
646: GmshNodesCreate(numTotalNodes, &nodes);
647: mesh->numNodes = numTotalNodes;
648: mesh->nodelist = nodes;
649: for (n = 0, block = 0; block < numEntityBlocks; ++block) {
650: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
651: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
652: if (byteSwap) PetscByteSwap(&numNodes, PETSC_LONG, 1);
653: if (gmsh->binary) {
654: size_t nbytes = sizeof(int) + 3 * sizeof(double);
655: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible uninitialized value */
656: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
657: PetscViewerRead(viewer, cbuf, numNodes * nbytes, NULL, PETSC_CHAR);
658: for (node = 0; node < numNodes; ++node, ++n) {
659: char *cnid = cbuf + node * nbytes, *cxyz = cnid + sizeof(int);
660: double *xyz = nodes->xyz + n * 3;
661: if (!PetscBinaryBigEndian()) PetscByteSwap(cnid, PETSC_ENUM, 1);
662: if (!PetscBinaryBigEndian()) PetscByteSwap(cxyz, PETSC_DOUBLE, 3);
663: PetscMemcpy(&nid, cnid, sizeof(int));
664: PetscMemcpy(xyz, cxyz, 3 * sizeof(double));
665: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
666: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
667: nodes->id[n] = nid;
668: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
669: }
670: } else {
671: for (node = 0; node < numNodes; ++node, ++n) {
672: double *xyz = nodes->xyz + n * 3;
673: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
674: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
675: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
676: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
677: nodes->id[n] = nid;
678: for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
679: }
680: }
681: }
682: return 0;
683: }
685: /*
686: $Elements
687: numEntityBlocks(unsigned long) numElements(unsigned long)
688: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
689: tag(int) numVert[...](int)
690: ...
691: ...
692: $EndElements
693: */
694: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
695: {
696: PetscViewer viewer = gmsh->viewer;
697: PetscBool byteSwap = gmsh->byteSwap;
698: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
699: int p, info[3], *ibuf = NULL;
700: int eid, dim, cellType, numVerts, numNodes, numTags;
701: GmshEntity *entity = NULL;
702: GmshElement *elements;
703: PetscInt *nodeMap = gmsh->nodeMap;
705: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
706: if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
707: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
708: if (byteSwap) PetscByteSwap(&numTotalElements, PETSC_LONG, 1);
709: GmshElementsCreate(numTotalElements, &elements);
710: mesh->numElems = numTotalElements;
711: mesh->elements = elements;
712: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
713: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
714: if (byteSwap) PetscByteSwap(info, PETSC_ENUM, 3);
715: eid = info[0];
716: dim = info[1];
717: cellType = info[2];
718: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
719: GmshCellTypeCheck(cellType);
720: numVerts = GmshCellMap[cellType].numVerts;
721: numNodes = GmshCellMap[cellType].numNodes;
722: numTags = entity->numTags;
723: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
724: if (byteSwap) PetscByteSwap(&numElements, PETSC_LONG, 1);
725: GmshBufferGet(gmsh, (1 + numNodes) * numElements, sizeof(int), &ibuf);
726: PetscViewerRead(viewer, ibuf, (1 + numNodes) * numElements, NULL, PETSC_ENUM);
727: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, (1 + numNodes) * numElements);
728: for (elem = 0; elem < numElements; ++elem, ++c) {
729: GmshElement *element = elements + c;
730: const int *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
731: element->id = id[0];
732: element->dim = dim;
733: element->cellType = cellType;
734: element->numVerts = numVerts;
735: element->numNodes = numNodes;
736: element->numTags = numTags;
737: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
738: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
739: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
740: }
741: }
742: return 0;
743: }
745: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
746: {
747: PetscViewer viewer = gmsh->viewer;
748: int fileFormat = gmsh->fileFormat;
749: PetscBool binary = gmsh->binary;
750: PetscBool byteSwap = gmsh->byteSwap;
751: int numPeriodic, snum, i;
752: char line[PETSC_MAX_PATH_LEN];
753: PetscInt *nodeMap = gmsh->nodeMap;
755: if (fileFormat == 22 || !binary) {
756: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
757: snum = sscanf(line, "%d", &numPeriodic);
759: } else {
760: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
761: if (byteSwap) PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);
762: }
763: for (i = 0; i < numPeriodic; i++) {
764: int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
765: long j, nNodes;
766: double affine[16];
768: if (fileFormat == 22 || !binary) {
769: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
770: snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
772: } else {
773: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
774: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);
775: correspondingDim = ibuf[0];
776: correspondingTag = ibuf[1];
777: primaryTag = ibuf[2];
778: }
779: (void)correspondingDim;
780: (void)correspondingTag;
781: (void)primaryTag; /* unused */
783: if (fileFormat == 22 || !binary) {
784: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
785: snum = sscanf(line, "%ld", &nNodes);
786: if (snum != 1) { /* discard transformation and try again */
787: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
788: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
789: snum = sscanf(line, "%ld", &nNodes);
791: }
792: } else {
793: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
794: if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
795: if (nNodes == -1) { /* discard transformation and try again */
796: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
797: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
798: if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
799: }
800: }
802: for (j = 0; j < nNodes; j++) {
803: if (fileFormat == 22 || !binary) {
804: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
805: snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
807: } else {
808: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
809: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 2);
810: correspondingNode = ibuf[0];
811: primaryNode = ibuf[1];
812: }
813: correspondingNode = (int)nodeMap[correspondingNode];
814: primaryNode = (int)nodeMap[primaryNode];
815: periodicMap[correspondingNode] = primaryNode;
816: }
817: }
818: return 0;
819: }
821: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
822: $Entities
823: numPoints(size_t) numCurves(size_t)
824: numSurfaces(size_t) numVolumes(size_t)
825: pointTag(int) X(double) Y(double) Z(double)
826: numPhysicalTags(size_t) physicalTag(int) ...
827: ...
828: curveTag(int) minX(double) minY(double) minZ(double)
829: maxX(double) maxY(double) maxZ(double)
830: numPhysicalTags(size_t) physicalTag(int) ...
831: numBoundingPoints(size_t) pointTag(int) ...
832: ...
833: surfaceTag(int) minX(double) minY(double) minZ(double)
834: maxX(double) maxY(double) maxZ(double)
835: numPhysicalTags(size_t) physicalTag(int) ...
836: numBoundingCurves(size_t) curveTag(int) ...
837: ...
838: volumeTag(int) minX(double) minY(double) minZ(double)
839: maxX(double) maxY(double) maxZ(double)
840: numPhysicalTags(size_t) physicalTag(int) ...
841: numBoundngSurfaces(size_t) surfaceTag(int) ...
842: ...
843: $EndEntities
844: */
845: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
846: {
847: PetscInt count[4], index, numTags, i;
848: int dim, eid, *tags = NULL;
849: GmshEntity *entity = NULL;
851: GmshReadSize(gmsh, count, 4);
852: GmshEntitiesCreate(count, &mesh->entities);
853: for (dim = 0; dim < 4; ++dim) {
854: for (index = 0; index < count[dim]; ++index) {
855: GmshReadInt(gmsh, &eid, 1);
856: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
857: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
858: GmshReadSize(gmsh, &numTags, 1);
859: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
860: GmshReadInt(gmsh, tags, numTags);
862: entity->numTags = numTags;
863: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
864: if (dim == 0) continue;
865: GmshReadSize(gmsh, &numTags, 1);
866: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
867: GmshReadInt(gmsh, tags, numTags);
868: /* Currently, we do not save the ids for the bounding entities */
869: }
870: }
871: return 0;
872: }
874: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
875: $Nodes
876: numEntityBlocks(size_t) numNodes(size_t)
877: minNodeTag(size_t) maxNodeTag(size_t)
878: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
879: nodeTag(size_t)
880: ...
881: x(double) y(double) z(double)
882: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
883: < v(double; if parametric and entityDim = 2) >
884: ...
885: ...
886: $EndNodes
887: */
888: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
889: {
890: int info[3], dim, eid, parametric;
891: PetscInt sizes[4], numEntityBlocks, numTags, t, numNodes, numNodesBlock = 0, block, node, n;
892: GmshEntity *entity = NULL;
893: GmshNodes *nodes;
895: GmshReadSize(gmsh, sizes, 4);
896: numEntityBlocks = sizes[0];
897: numNodes = sizes[1];
898: GmshNodesCreate(numNodes, &nodes);
899: mesh->numNodes = numNodes;
900: mesh->nodelist = nodes;
901: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
902: GmshReadInt(gmsh, info, 3);
903: dim = info[0];
904: eid = info[1];
905: parametric = info[2];
906: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
907: numTags = entity->numTags;
909: GmshReadSize(gmsh, &numNodesBlock, 1);
910: GmshReadSize(gmsh, nodes->id + node, numNodesBlock);
911: GmshReadDouble(gmsh, nodes->xyz + node * 3, numNodesBlock * 3);
912: for (n = 0; n < numNodesBlock; ++n) {
913: PetscInt *tags = &nodes->tag[node * GMSH_MAX_TAGS];
915: for (t = 0; t < numTags; ++t) tags[n * GMSH_MAX_TAGS + t] = entity->tags[t];
916: for (t = numTags; t < GMSH_MAX_TAGS; ++t) tags[n * GMSH_MAX_TAGS + t] = -1;
917: }
918: }
919: gmsh->nodeStart = sizes[2];
920: gmsh->nodeEnd = sizes[3] + 1;
921: return 0;
922: }
924: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
925: $Elements
926: numEntityBlocks(size_t) numElements(size_t)
927: minElementTag(size_t) maxElementTag(size_t)
928: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
929: elementTag(size_t) nodeTag(size_t) ...
930: ...
931: ...
932: $EndElements
933: */
934: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
935: {
936: int info[3], eid, dim, cellType;
937: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
938: GmshEntity *entity = NULL;
939: GmshElement *elements;
940: PetscInt *nodeMap = gmsh->nodeMap;
942: GmshReadSize(gmsh, sizes, 4);
943: numEntityBlocks = sizes[0];
944: numElements = sizes[1];
945: GmshElementsCreate(numElements, &elements);
946: mesh->numElems = numElements;
947: mesh->elements = elements;
948: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
949: GmshReadInt(gmsh, info, 3);
950: dim = info[0];
951: eid = info[1];
952: cellType = info[2];
953: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
954: GmshCellTypeCheck(cellType);
955: numVerts = GmshCellMap[cellType].numVerts;
956: numNodes = GmshCellMap[cellType].numNodes;
957: numTags = entity->numTags;
958: GmshReadSize(gmsh, &numBlockElements, 1);
959: GmshBufferGet(gmsh, (1 + numNodes) * numBlockElements, sizeof(PetscInt), &ibuf);
960: GmshReadSize(gmsh, ibuf, (1 + numNodes) * numBlockElements);
961: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
962: GmshElement *element = elements + c;
963: const PetscInt *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
964: element->id = id[0];
965: element->dim = dim;
966: element->cellType = cellType;
967: element->numVerts = numVerts;
968: element->numNodes = numNodes;
969: element->numTags = numTags;
970: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
971: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
972: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
973: }
974: }
975: return 0;
976: }
978: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
979: $Periodic
980: numPeriodicLinks(size_t)
981: entityDim(int) entityTag(int) entityTagPrimary(int)
982: numAffine(size_t) value(double) ...
983: numCorrespondingNodes(size_t)
984: nodeTag(size_t) nodeTagPrimary(size_t)
985: ...
986: ...
987: $EndPeriodic
988: */
989: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
990: {
991: int info[3];
992: double dbuf[16];
993: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
994: PetscInt *nodeMap = gmsh->nodeMap;
996: GmshReadSize(gmsh, &numPeriodicLinks, 1);
997: for (link = 0; link < numPeriodicLinks; ++link) {
998: GmshReadInt(gmsh, info, 3);
999: GmshReadSize(gmsh, &numAffine, 1);
1000: GmshReadDouble(gmsh, dbuf, numAffine);
1001: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
1002: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
1003: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes * 2);
1004: for (node = 0; node < numCorrespondingNodes; ++node) {
1005: PetscInt correspondingNode = nodeMap[nodeTags[node * 2 + 0]];
1006: PetscInt primaryNode = nodeMap[nodeTags[node * 2 + 1]];
1007: periodicMap[correspondingNode] = primaryNode;
1008: }
1009: }
1010: return 0;
1011: }
1013: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1014: $MeshFormat // same as MSH version 2
1015: version(ASCII double; currently 4.1)
1016: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1017: data-size(ASCII int; sizeof(size_t))
1018: < int with value one; only in binary mode, to detect endianness >
1019: $EndMeshFormat
1020: */
1021: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1022: {
1023: char line[PETSC_MAX_PATH_LEN];
1024: int snum, fileType, fileFormat, dataSize, checkEndian;
1025: float version;
1027: GmshReadString(gmsh, line, 3);
1028: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1029: fileFormat = (int)roundf(version * 10);
1038: gmsh->fileFormat = fileFormat;
1039: gmsh->dataSize = dataSize;
1040: gmsh->byteSwap = PETSC_FALSE;
1041: if (gmsh->binary) {
1042: GmshReadInt(gmsh, &checkEndian, 1);
1043: if (checkEndian != 1) {
1044: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
1046: gmsh->byteSwap = PETSC_TRUE;
1047: }
1048: }
1049: return 0;
1050: }
1052: /*
1053: PhysicalNames
1054: numPhysicalNames(ASCII int)
1055: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1056: ...
1057: $EndPhysicalNames
1058: */
1059: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1060: {
1061: char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p, *q, *r;
1062: int snum, region, dim, tag;
1064: GmshReadString(gmsh, line, 1);
1065: snum = sscanf(line, "%d", ®ion);
1066: mesh->numRegions = region;
1068: PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames);
1069: for (region = 0; region < mesh->numRegions; ++region) {
1070: GmshReadString(gmsh, line, 2);
1071: snum = sscanf(line, "%d %d", &dim, &tag);
1073: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
1074: PetscStrchr(line, '"', &p);
1076: PetscStrrchr(line, '"', &q);
1078: PetscStrrchr(line, ':', &r);
1079: if (p != r) q = r;
1080: PetscStrncpy(name, p + 1, (size_t)(q - p - 1));
1081: mesh->regionDims[region] = dim;
1082: mesh->regionTags[region] = tag;
1083: PetscStrallocpy(name, &mesh->regionNames[region]);
1084: }
1085: return 0;
1086: }
1088: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1089: {
1090: switch (gmsh->fileFormat) {
1091: case 41:
1092: GmshReadEntities_v41(gmsh, mesh);
1093: break;
1094: default:
1095: GmshReadEntities_v40(gmsh, mesh);
1096: break;
1097: }
1098: return 0;
1099: }
1101: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1102: {
1103: switch (gmsh->fileFormat) {
1104: case 41:
1105: GmshReadNodes_v41(gmsh, mesh);
1106: break;
1107: case 40:
1108: GmshReadNodes_v40(gmsh, mesh);
1109: break;
1110: default:
1111: GmshReadNodes_v22(gmsh, mesh);
1112: break;
1113: }
1115: { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1116: if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1117: PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1118: GmshNodes *nodes = mesh->nodelist;
1119: for (n = 0; n < mesh->numNodes; ++n) {
1120: const PetscInt tag = nodes->id[n];
1121: tagMin = PetscMin(tag, tagMin);
1122: tagMax = PetscMax(tag, tagMax);
1123: }
1124: gmsh->nodeStart = tagMin;
1125: gmsh->nodeEnd = tagMax + 1;
1126: }
1127: }
1129: { /* Support for sparse node tags */
1130: PetscInt n, t;
1131: GmshNodes *nodes = mesh->nodelist;
1132: PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);
1133: for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1134: gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1135: for (n = 0; n < mesh->numNodes; ++n) {
1136: const PetscInt tag = nodes->id[n];
1138: gmsh->nodeMap[tag] = n;
1139: }
1140: }
1141: return 0;
1142: }
1144: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1145: {
1146: switch (gmsh->fileFormat) {
1147: case 41:
1148: GmshReadElements_v41(gmsh, mesh);
1149: break;
1150: case 40:
1151: GmshReadElements_v40(gmsh, mesh);
1152: break;
1153: default:
1154: GmshReadElements_v22(gmsh, mesh);
1155: break;
1156: }
1158: { /* Reorder elements by codimension and polytope type */
1159: PetscInt ne = mesh->numElems;
1160: GmshElement *elements = mesh->elements;
1161: PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0;
1162: PetscInt offset[GMSH_NUM_POLYTOPES + 1], e, k;
1164: for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1165: PetscMemzero(offset, sizeof(offset));
1167: keymap[GMSH_TET] = nk++;
1168: keymap[GMSH_HEX] = nk++;
1169: keymap[GMSH_PRI] = nk++;
1170: keymap[GMSH_PYR] = nk++;
1171: keymap[GMSH_TRI] = nk++;
1172: keymap[GMSH_QUA] = nk++;
1173: keymap[GMSH_SEG] = nk++;
1174: keymap[GMSH_VTX] = nk++;
1176: GmshElementsCreate(mesh->numElems, &mesh->elements);
1177: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1178: for (e = 0; e < ne; ++e) offset[1 + key(e)]++;
1179: for (k = 1; k < nk; ++k) offset[k] += offset[k - 1];
1180: for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1181: #undef key
1182: GmshElementsDestroy(&elements);
1183: }
1185: { /* Mesh dimension and order */
1186: GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1187: mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0;
1188: mesh->order = elem ? GmshCellMap[elem->cellType].order : 0;
1189: }
1191: {
1192: PetscBT vtx;
1193: PetscInt dim = mesh->dim, e, n, v;
1195: PetscBTCreate(mesh->numNodes, &vtx);
1197: /* Compute number of cells and set of vertices */
1198: mesh->numCells = 0;
1199: for (e = 0; e < mesh->numElems; ++e) {
1200: GmshElement *elem = mesh->elements + e;
1201: if (elem->dim == dim && dim > 0) mesh->numCells++;
1202: for (v = 0; v < elem->numVerts; v++) PetscBTSet(vtx, elem->nodes[v]);
1203: }
1205: /* Compute numbering for vertices */
1206: mesh->numVerts = 0;
1207: PetscMalloc1(mesh->numNodes, &mesh->vertexMap);
1208: for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;
1210: PetscBTDestroy(&vtx);
1211: }
1212: return 0;
1213: }
1215: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1216: {
1217: PetscInt n;
1219: PetscMalloc1(mesh->numNodes, &mesh->periodMap);
1220: for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1221: switch (gmsh->fileFormat) {
1222: case 41:
1223: GmshReadPeriodic_v41(gmsh, mesh->periodMap);
1224: break;
1225: default:
1226: GmshReadPeriodic_v40(gmsh, mesh->periodMap);
1227: break;
1228: }
1230: /* Find canonical primary nodes */
1231: for (n = 0; n < mesh->numNodes; ++n)
1232: while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];
1234: /* Renumber vertices (filter out correspondings) */
1235: mesh->numVerts = 0;
1236: for (n = 0; n < mesh->numNodes; ++n)
1237: if (mesh->vertexMap[n] >= 0) /* is vertex */
1238: if (mesh->periodMap[n] == n) /* is primary */
1239: mesh->vertexMap[n] = mesh->numVerts++;
1240: for (n = 0; n < mesh->numNodes; ++n)
1241: if (mesh->vertexMap[n] >= 0) /* is vertex */
1242: if (mesh->periodMap[n] != n) /* is corresponding */
1243: mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1244: return 0;
1245: }
1247: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1248: static const DMPolytopeType DMPolytopeMap[] = {
1249: /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1250: /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1251: /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1252: /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1253: /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1254: /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1255: /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1256: /* GMSH_PYR */ DM_POLYTOPE_PYRAMID, DM_POLYTOPE_UNKNOWN};
1258: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1259: {
1260: return DMPolytopeMap[GmshCellMap[cellType].polytope];
1261: }
1263: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1264: {
1265: DM K;
1266: PetscSpace P;
1267: PetscDualSpace Q;
1268: PetscQuadrature q, fq;
1269: PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1270: PetscBool endpoint = PETSC_TRUE;
1271: char name[32];
1273: /* Create space */
1274: PetscSpaceCreate(comm, &P);
1275: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1276: PetscSpacePolynomialSetTensor(P, isTensor);
1277: PetscSpaceSetNumComponents(P, Nc);
1278: PetscSpaceSetNumVariables(P, dim);
1279: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1280: if (prefix) {
1281: PetscObjectSetOptionsPrefix((PetscObject)P, prefix);
1282: PetscSpaceSetFromOptions(P);
1283: PetscObjectSetOptionsPrefix((PetscObject)P, NULL);
1284: PetscSpaceGetDegree(P, &k, NULL);
1285: }
1286: PetscSpaceSetUp(P);
1287: /* Create dual space */
1288: PetscDualSpaceCreate(comm, &Q);
1289: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1290: PetscDualSpaceLagrangeSetTensor(Q, isTensor);
1291: PetscDualSpaceLagrangeSetContinuity(Q, continuity);
1292: PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
1293: PetscDualSpaceSetNumComponents(Q, Nc);
1294: PetscDualSpaceSetOrder(Q, k);
1295: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
1296: PetscDualSpaceSetDM(Q, K);
1297: DMDestroy(&K);
1298: if (prefix) {
1299: PetscObjectSetOptionsPrefix((PetscObject)Q, prefix);
1300: PetscDualSpaceSetFromOptions(Q);
1301: PetscObjectSetOptionsPrefix((PetscObject)Q, NULL);
1302: }
1303: PetscDualSpaceSetUp(Q);
1304: /* Create quadrature */
1305: if (isSimplex) {
1306: PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q);
1307: PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq);
1308: } else {
1309: PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q);
1310: PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq);
1311: }
1312: /* Create finite element */
1313: PetscFECreate(comm, fem);
1314: PetscFESetType(*fem, PETSCFEBASIC);
1315: PetscFESetNumComponents(*fem, Nc);
1316: PetscFESetBasisSpace(*fem, P);
1317: PetscFESetDualSpace(*fem, Q);
1318: PetscFESetQuadrature(*fem, q);
1319: PetscFESetFaceQuadrature(*fem, fq);
1320: if (prefix) {
1321: PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix);
1322: PetscFESetFromOptions(*fem);
1323: PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL);
1324: }
1325: PetscFESetUp(*fem);
1326: /* Cleanup */
1327: PetscSpaceDestroy(&P);
1328: PetscDualSpaceDestroy(&Q);
1329: PetscQuadratureDestroy(&q);
1330: PetscQuadratureDestroy(&fq);
1331: /* Set finite element name */
1332: PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k);
1333: PetscFESetName(*fem, name);
1334: return 0;
1335: }
1337: /*@C
1338: DMPlexCreateGmshFromFile - Create a `DMPLEX` mesh from a Gmsh file
1340: + comm - The MPI communicator
1341: . filename - Name of the Gmsh file
1342: - interpolate - Create faces and edges in the mesh
1344: Output Parameter:
1345: . dm - The `DM` object representing the mesh
1347: Level: beginner
1349: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
1350: @*/
1351: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1352: {
1353: PetscViewer viewer;
1354: PetscMPIInt rank;
1355: int fileType;
1356: PetscViewerType vtype;
1358: MPI_Comm_rank(comm, &rank);
1360: /* Determine Gmsh file type (ASCII or binary) from file header */
1361: if (rank == 0) {
1362: GmshFile gmsh[1];
1363: char line[PETSC_MAX_PATH_LEN];
1364: int snum;
1365: float version;
1366: int fileFormat;
1368: PetscArrayzero(gmsh, 1);
1369: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1370: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1371: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1372: PetscViewerFileSetName(gmsh->viewer, filename);
1373: /* Read only the first two lines of the Gmsh file */
1374: GmshReadSection(gmsh, line);
1375: GmshExpect(gmsh, "$MeshFormat", line);
1376: GmshReadString(gmsh, line, 2);
1377: snum = sscanf(line, "%f %d", &version, &fileType);
1378: fileFormat = (int)roundf(version * 10);
1383: PetscViewerDestroy(&gmsh->viewer);
1384: }
1385: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1386: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1388: /* Create appropriate viewer and build plex */
1389: PetscViewerCreate(comm, &viewer);
1390: PetscViewerSetType(viewer, vtype);
1391: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1392: PetscViewerFileSetName(viewer, filename);
1393: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1394: PetscViewerDestroy(&viewer);
1395: return 0;
1396: }
1398: /*@
1399: DMPlexCreateGmsh - Create a `DMPLEX` mesh from a Gmsh file viewer
1401: Collective
1403: Input Parameters:
1404: + comm - The MPI communicator
1405: . viewer - The `PetscViewer` associated with a Gmsh file
1406: - interpolate - Create faces and edges in the mesh
1408: Output Parameter:
1409: . dm - The `DM` object representing the mesh
1411: Options Database Keys:
1412: + -dm_plex_gmsh_hybrid - Force triangular prisms to use tensor order
1413: . -dm_plex_gmsh_periodic - Read Gmsh periodic section and construct a periodic Plex
1414: . -dm_plex_gmsh_highorder - Generate high-order coordinates
1415: . -dm_plex_gmsh_project - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
1416: . -dm_plex_gmsh_use_regions - Generate labels with region names
1417: . -dm_plex_gmsh_mark_vertices - Add vertices to generated labels
1418: . -dm_plex_gmsh_multiple_tags - Allow multiple tags for default labels
1419: - -dm_plex_gmsh_spacedim <d> - Embedding space dimension, if different from topological dimension
1421: Note:
1422: The Gmsh file format is described in http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1424: 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.
1426: Level: beginner
1428: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1429: @*/
1430: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1431: {
1432: GmshMesh *mesh = NULL;
1433: PetscViewer parentviewer = NULL;
1434: PetscBT periodicVerts = NULL;
1435: PetscBT periodicCells = NULL;
1436: DM cdm, cdmCell = NULL;
1437: PetscSection cs, csCell = NULL;
1438: Vec coordinates, coordinatesCell;
1439: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1440: PetscInt dim = 0, coordDim = -1, order = 0;
1441: PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0;
1442: PetscInt cell, cone[8], e, n, v, d;
1443: PetscBool binary, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE;
1444: PetscBool hybrid = interpolate, periodic = PETSC_TRUE;
1445: PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1446: PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1447: PetscMPIInt rank;
1449: MPI_Comm_rank(comm, &rank);
1450: PetscObjectOptionsBegin((PetscObject)viewer);
1451: PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options");
1452: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);
1453: PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1454: PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);
1455: PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);
1456: PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL);
1457: PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL);
1458: PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL);
1459: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);
1460: PetscOptionsHeadEnd();
1461: PetscOptionsEnd();
1463: GmshCellInfoSetUp();
1465: DMCreate(comm, dm);
1466: DMSetType(*dm, DMPLEX);
1467: PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL);
1469: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1471: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1472: if (binary) {
1473: parentviewer = viewer;
1474: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1475: }
1477: if (rank == 0) {
1478: GmshFile gmsh[1];
1479: char line[PETSC_MAX_PATH_LEN];
1480: PetscBool match;
1482: PetscArrayzero(gmsh, 1);
1483: gmsh->viewer = viewer;
1484: gmsh->binary = binary;
1486: GmshMeshCreate(&mesh);
1488: /* Read mesh format */
1489: GmshReadSection(gmsh, line);
1490: GmshExpect(gmsh, "$MeshFormat", line);
1491: GmshReadMeshFormat(gmsh);
1492: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1494: /* OPTIONAL Read physical names */
1495: GmshReadSection(gmsh, line);
1496: GmshMatch(gmsh, "$PhysicalNames", line, &match);
1497: if (match) {
1498: GmshExpect(gmsh, "$PhysicalNames", line);
1499: GmshReadPhysicalNames(gmsh, mesh);
1500: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1501: /* Initial read for entity section */
1502: GmshReadSection(gmsh, line);
1503: }
1505: /* Read entities */
1506: if (gmsh->fileFormat >= 40) {
1507: GmshExpect(gmsh, "$Entities", line);
1508: GmshReadEntities(gmsh, mesh);
1509: GmshReadEndSection(gmsh, "$EndEntities", line);
1510: /* Initial read for nodes section */
1511: GmshReadSection(gmsh, line);
1512: }
1514: /* Read nodes */
1515: GmshExpect(gmsh, "$Nodes", line);
1516: GmshReadNodes(gmsh, mesh);
1517: GmshReadEndSection(gmsh, "$EndNodes", line);
1519: /* Read elements */
1520: GmshReadSection(gmsh, line);
1521: GmshExpect(gmsh, "$Elements", line);
1522: GmshReadElements(gmsh, mesh);
1523: GmshReadEndSection(gmsh, "$EndElements", line);
1525: /* Read periodic section (OPTIONAL) */
1526: if (periodic) {
1527: GmshReadSection(gmsh, line);
1528: GmshMatch(gmsh, "$Periodic", line, &periodic);
1529: }
1530: if (periodic) {
1531: GmshExpect(gmsh, "$Periodic", line);
1532: GmshReadPeriodic(gmsh, mesh);
1533: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1534: }
1536: PetscFree(gmsh->wbuf);
1537: PetscFree(gmsh->sbuf);
1538: PetscFree(gmsh->nbuf);
1540: dim = mesh->dim;
1541: order = mesh->order;
1542: numNodes = mesh->numNodes;
1543: numElems = mesh->numElems;
1544: numVerts = mesh->numVerts;
1545: numCells = mesh->numCells;
1547: {
1548: GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1549: GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1550: int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1551: int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1552: isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1553: isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1554: hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1555: }
1556: }
1558: if (parentviewer) PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1560: {
1561: int buf[6];
1563: buf[0] = (int)dim;
1564: buf[1] = (int)order;
1565: buf[2] = periodic;
1566: buf[3] = isSimplex;
1567: buf[4] = isHybrid;
1568: buf[5] = hasTetra;
1570: MPI_Bcast(buf, 6, MPI_INT, 0, comm);
1572: dim = buf[0];
1573: order = buf[1];
1574: periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1575: isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1576: isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1577: hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1578: }
1580: if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1583: /* We do not want this label automatically computed, instead we fill it here */
1584: DMCreateLabel(*dm, "celltype");
1586: /* Allocate the cell-vertex mesh */
1587: DMPlexSetChart(*dm, 0, numCells + numVerts);
1588: for (cell = 0; cell < numCells; ++cell) {
1589: GmshElement *elem = mesh->elements + cell;
1590: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1591: if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1592: DMPlexSetConeSize(*dm, cell, elem->numVerts);
1593: DMPlexSetCellType(*dm, cell, ctype);
1594: }
1595: for (v = numCells; v < numCells + numVerts; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1596: DMSetUp(*dm);
1598: /* Add cell-vertex connections */
1599: for (cell = 0; cell < numCells; ++cell) {
1600: GmshElement *elem = mesh->elements + cell;
1601: for (v = 0; v < elem->numVerts; ++v) {
1602: const PetscInt nn = elem->nodes[v];
1603: const PetscInt vv = mesh->vertexMap[nn];
1604: cone[v] = numCells + vv;
1605: }
1606: DMPlexReorderCell(*dm, cell, cone);
1607: DMPlexSetCone(*dm, cell, cone);
1608: }
1610: DMSetDimension(*dm, dim);
1611: DMPlexSymmetrize(*dm);
1612: DMPlexStratify(*dm);
1613: if (interpolate) {
1614: DM idm;
1616: DMPlexInterpolate(*dm, &idm);
1617: DMDestroy(dm);
1618: *dm = idm;
1619: }
1621: if (rank == 0) {
1622: const PetscInt Nr = useregions ? mesh->numRegions : 0;
1623: PetscInt vStart, vEnd;
1625: PetscCalloc1(Nr, ®ionSets);
1626: DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);
1627: for (cell = 0, e = 0; e < numElems; ++e) {
1628: GmshElement *elem = mesh->elements + e;
1630: /* Create cell sets */
1631: if (elem->dim == dim && dim > 0) {
1632: if (elem->numTags > 0) {
1633: PetscInt Nt = elem->numTags, t, r;
1635: for (t = 0; t < Nt; ++t) {
1636: const PetscInt tag = elem->tags[t];
1637: const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1639: if (generic) DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag);
1640: for (r = 0; r < Nr; ++r) {
1641: if (mesh->regionDims[r] != dim) continue;
1642: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag);
1643: }
1644: }
1645: }
1646: cell++;
1647: }
1649: /* Create face sets */
1650: if (interpolate && elem->dim == dim - 1) {
1651: PetscInt joinSize;
1652: const PetscInt *join = NULL;
1653: PetscInt Nt = elem->numTags, t, r;
1655: /* Find the relevant facet with vertex joins */
1656: for (v = 0; v < elem->numVerts; ++v) {
1657: const PetscInt nn = elem->nodes[v];
1658: const PetscInt vv = mesh->vertexMap[nn];
1659: cone[v] = vStart + vv;
1660: }
1661: DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1663: for (t = 0; t < Nt; ++t) {
1664: const PetscInt tag = elem->tags[t];
1665: const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1667: if (generic) DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag);
1668: for (r = 0; r < Nr; ++r) {
1669: if (mesh->regionDims[r] != dim - 1) continue;
1670: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag);
1671: }
1672: }
1673: DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1674: }
1676: /* Create vertex sets */
1677: if (elem->dim == 0 && markvertices) {
1678: if (elem->numTags > 0) {
1679: const PetscInt nn = elem->nodes[0];
1680: const PetscInt vv = mesh->vertexMap[nn];
1681: const PetscInt tag = elem->tags[0];
1682: PetscInt r;
1684: if (!Nr) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1685: for (r = 0; r < Nr; ++r) {
1686: if (mesh->regionDims[r] != 0) continue;
1687: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag);
1688: }
1689: }
1690: }
1691: }
1692: if (markvertices) {
1693: for (v = 0; v < numNodes; ++v) {
1694: const PetscInt vv = mesh->vertexMap[v];
1695: const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS];
1696: PetscInt r, t;
1698: for (t = 0; t < GMSH_MAX_TAGS; ++t) {
1699: const PetscInt tag = tags[t];
1700: const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;
1702: if (tag == -1) continue;
1703: if (generic) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1704: for (r = 0; r < Nr; ++r) {
1705: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag);
1706: }
1707: }
1708: }
1709: }
1710: PetscFree(regionSets);
1711: }
1713: { /* Create Cell/Face/Vertex Sets labels at all processes */
1714: enum {
1715: n = 4
1716: };
1717: PetscBool flag[n];
1719: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1720: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1721: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1722: flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1723: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1724: if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1725: if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1726: if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1727: if (flag[3]) DMCreateLabel(*dm, "marker");
1728: }
1730: if (periodic) {
1731: PetscBTCreate(numVerts, &periodicVerts);
1732: for (n = 0; n < numNodes; ++n) {
1733: if (mesh->vertexMap[n] >= 0) {
1734: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1735: PetscInt m = mesh->periodMap[n];
1736: PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1737: PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1738: }
1739: }
1740: }
1741: PetscBTCreate(numCells, &periodicCells);
1742: for (cell = 0; cell < numCells; ++cell) {
1743: GmshElement *elem = mesh->elements + cell;
1744: for (v = 0; v < elem->numVerts; ++v) {
1745: PetscInt nn = elem->nodes[v];
1746: PetscInt vv = mesh->vertexMap[nn];
1747: if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1748: PetscBTSet(periodicCells, cell);
1749: break;
1750: }
1751: }
1752: }
1753: }
1755: /* Setup coordinate DM */
1756: if (coordDim < 0) coordDim = dim;
1757: DMSetCoordinateDim(*dm, coordDim);
1758: DMGetCoordinateDM(*dm, &cdm);
1759: if (highOrder) {
1760: PetscFE fe;
1761: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1762: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1764: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1766: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1767: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1768: DMSetField(cdm, 0, NULL, (PetscObject)fe);
1769: PetscFEDestroy(&fe);
1770: DMCreateDS(cdm);
1771: }
1772: if (periodic) {
1773: DMClone(cdm, &cdmCell);
1774: DMSetCellCoordinateDM(*dm, cdmCell);
1775: }
1777: /* Create coordinates */
1778: if (highOrder) {
1779: PetscInt maxDof = GmshNumNodes_HEX(order) * coordDim;
1780: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1781: PetscSection section;
1782: PetscScalar *cellCoords;
1784: DMSetLocalSection(cdm, NULL);
1785: DMGetLocalSection(cdm, &cs);
1786: PetscSectionClone(cs, §ion);
1787: DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1789: DMCreateLocalVector(cdm, &coordinates);
1790: PetscMalloc1(maxDof, &cellCoords);
1791: for (cell = 0; cell < numCells; ++cell) {
1792: GmshElement *elem = mesh->elements + cell;
1793: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1794: int s = 0;
1795: for (n = 0; n < elem->numNodes; ++n) {
1796: while (lexorder[n + s] < 0) ++s;
1797: const PetscInt node = elem->nodes[lexorder[n + s]];
1798: for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d];
1799: }
1800: if (s) {
1801: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1802: PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1803: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1804: 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};
1805: 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};
1806: 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};
1807: 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};
1808: 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};
1809: 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};
1810: 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};
1811: PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
1812: PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
1813: NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL};
1814: PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3};
1816: /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
1817: for (n = 0; n < elem->numNodes + s; ++n) {
1818: if (lexorder[n] >= 0) continue;
1819: for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0;
1820: for (int bn = 0; bn < elem->numNodes + s; ++bn) {
1821: if (lexorder[bn] < 0) continue;
1822: const PetscReal *weights = sdWeights[coordDim][n];
1823: const PetscInt bnode = elem->nodes[lexorder[bn]];
1824: for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d];
1825: }
1826: }
1827: }
1828: DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1829: }
1830: PetscSectionDestroy(§ion);
1831: PetscFree(cellCoords);
1832: } else {
1833: PetscInt *nodeMap;
1834: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1835: PetscScalar *pointCoords;
1837: DMGetCoordinateSection(*dm, &cs);
1838: PetscSectionSetNumFields(cs, 1);
1839: PetscSectionSetFieldComponents(cs, 0, coordDim);
1840: PetscSectionSetChart(cs, numCells, numCells + numVerts);
1841: for (v = numCells; v < numCells + numVerts; ++v) {
1842: PetscSectionSetDof(cs, v, coordDim);
1843: PetscSectionSetFieldDof(cs, v, 0, coordDim);
1844: }
1845: PetscSectionSetUp(cs);
1847: /* We need to localize coordinates on cells */
1848: if (periodic) {
1849: PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell);
1850: PetscSectionSetNumFields(csCell, 1);
1851: PetscSectionSetFieldComponents(csCell, 0, coordDim);
1852: PetscSectionSetChart(csCell, 0, numCells);
1853: for (cell = 0; cell < numCells; ++cell) {
1854: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1855: GmshElement *elem = mesh->elements + cell;
1856: PetscInt dof = elem->numVerts * coordDim;
1858: PetscSectionSetDof(csCell, cell, dof);
1859: PetscSectionSetFieldDof(csCell, cell, 0, dof);
1860: }
1861: }
1862: PetscSectionSetUp(csCell);
1863: DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell);
1864: }
1866: DMCreateLocalVector(cdm, &coordinates);
1867: VecGetArray(coordinates, &pointCoords);
1868: PetscMalloc1(numVerts, &nodeMap);
1869: for (n = 0; n < numNodes; n++)
1870: if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n;
1871: for (v = 0; v < numVerts; ++v) {
1872: PetscInt off, node = nodeMap[v];
1874: PetscSectionGetOffset(cs, numCells + v, &off);
1875: for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
1876: }
1877: VecRestoreArray(coordinates, &pointCoords);
1878: PetscFree(nodeMap);
1880: if (periodic) {
1881: DMCreateLocalVector(cdmCell, &coordinatesCell);
1882: VecGetArray(coordinatesCell, &pointCoords);
1883: for (cell = 0; cell < numCells; ++cell) {
1884: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1885: GmshElement *elem = mesh->elements + cell;
1886: PetscInt off, node;
1887: for (v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v];
1888: DMPlexReorderCell(cdmCell, cell, cone);
1889: PetscSectionGetOffset(csCell, cell, &off);
1890: for (v = 0; v < elem->numVerts; ++v)
1891: for (node = cone[v], d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[node * 3 + d];
1892: }
1893: }
1894: VecSetBlockSize(coordinatesCell, coordDim);
1895: VecRestoreArray(coordinatesCell, &pointCoords);
1896: DMSetCellCoordinatesLocal(*dm, coordinatesCell);
1897: VecDestroy(&coordinatesCell);
1898: }
1899: PetscSectionDestroy(&csCell);
1900: DMDestroy(&cdmCell);
1901: }
1903: PetscObjectSetName((PetscObject)coordinates, "coordinates");
1904: VecSetBlockSize(coordinates, coordDim);
1905: DMSetCoordinatesLocal(*dm, coordinates);
1906: VecDestroy(&coordinates);
1908: GmshMeshDestroy(&mesh);
1909: PetscBTDestroy(&periodicVerts);
1910: PetscBTDestroy(&periodicCells);
1912: if (highOrder && project) {
1913: PetscFE fe;
1914: const char prefix[] = "dm_plex_gmsh_project_";
1915: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1916: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
1918: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1920: GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1921: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1922: DMProjectCoordinates(*dm, fe);
1923: PetscFEDestroy(&fe);
1924: }
1926: PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL);
1927: return 0;
1928: }