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; lex[2] = 1; lex[8] = 2; lex[6] = 3;
23: /* Edges */
24: lex[1] = 4; lex[5] = 5; lex[7] = 6; lex[3] = 7;
25: /* Cell */
26: lex[4] = -8;
27: }
28: return lex;
29: }
31: static int *GmshLexOrder_HEX_2_Serendipity(void)
32: {
33: static int Gmsh_LexOrder_HEX_2_Serendipity[27] = {-1};
34: int *lex = Gmsh_LexOrder_HEX_2_Serendipity;
35: if (lex[0] == -1) {
36: /* Vertices */
37: lex[ 0] = 0; lex[ 2] = 1; lex[ 8] = 2; lex[ 6] = 3;
38: lex[18] = 4; lex[20] = 5; lex[26] = 6; lex[24] = 7;
39: /* Edges */
40: lex[ 1] = 8; lex[ 3] = 9; lex[ 9] = 10; lex[ 5] = 11;
41: lex[11] = 12; lex[ 7] = 13; lex[17] = 14; lex[15] = 15;
42: lex[19] = 16; lex[21] = 17; lex[23] = 18; lex[25] = 19;
43: /* Faces */
44: lex[ 4] = -20; lex[10] = -21; lex[12] = -22; lex[14] = -23;
45: lex[16] = -24; lex[22] = -25;
46: /* Cell */
47: lex[13] = -26;
48: }
49: return lex;
50: }
52: #define GMSH_LEXORDER_LIST(T) \
53: GMSH_LEXORDER_ITEM(T, 1) \
54: GMSH_LEXORDER_ITEM(T, 2) \
55: GMSH_LEXORDER_ITEM(T, 3) \
56: GMSH_LEXORDER_ITEM(T, 4) \
57: GMSH_LEXORDER_ITEM(T, 5) \
58: GMSH_LEXORDER_ITEM(T, 6) \
59: GMSH_LEXORDER_ITEM(T, 7) \
60: GMSH_LEXORDER_ITEM(T, 8) \
61: GMSH_LEXORDER_ITEM(T, 9) \
62: GMSH_LEXORDER_ITEM(T, 10)
64: GMSH_LEXORDER_ITEM(VTX, 0)
65: GMSH_LEXORDER_LIST(SEG)
66: GMSH_LEXORDER_LIST(TRI)
67: GMSH_LEXORDER_LIST(QUA)
68: GMSH_LEXORDER_LIST(TET)
69: GMSH_LEXORDER_LIST(HEX)
70: GMSH_LEXORDER_LIST(PRI)
71: GMSH_LEXORDER_LIST(PYR)
73: typedef enum {
74: GMSH_VTX = 0,
75: GMSH_SEG = 1,
76: GMSH_TRI = 2,
77: GMSH_QUA = 3,
78: GMSH_TET = 4,
79: GMSH_HEX = 5,
80: GMSH_PRI = 6,
81: GMSH_PYR = 7,
82: GMSH_NUM_POLYTOPES = 8
83: } GmshPolytopeType;
85: typedef struct {
86: int cellType;
87: int polytope;
88: int dim;
89: int order;
90: int numVerts;
91: int numNodes;
92: int* (*lexorder)(void);
93: } GmshCellInfo;
95: #define GmshCellEntry(cellType, polytope, dim, order) \
96: {cellType, GMSH_##polytope, dim, order, \
97: GmshNumNodes_##polytope(1), \
98: GmshNumNodes_##polytope(order), \
99: GmshLexOrder_##polytope##_##order}
101: static const GmshCellInfo GmshCellTable[] = {
102: GmshCellEntry( 15, VTX, 0, 0),
104: GmshCellEntry( 1, SEG, 1, 1),
105: GmshCellEntry( 8, SEG, 1, 2),
106: GmshCellEntry( 26, SEG, 1, 3),
107: GmshCellEntry( 27, SEG, 1, 4),
108: GmshCellEntry( 28, SEG, 1, 5),
109: GmshCellEntry( 62, SEG, 1, 6),
110: GmshCellEntry( 63, SEG, 1, 7),
111: GmshCellEntry( 64, SEG, 1, 8),
112: GmshCellEntry( 65, SEG, 1, 9),
113: GmshCellEntry( 66, SEG, 1, 10),
115: GmshCellEntry( 2, TRI, 2, 1),
116: GmshCellEntry( 9, TRI, 2, 2),
117: GmshCellEntry( 21, TRI, 2, 3),
118: GmshCellEntry( 23, TRI, 2, 4),
119: GmshCellEntry( 25, TRI, 2, 5),
120: GmshCellEntry( 42, TRI, 2, 6),
121: GmshCellEntry( 43, TRI, 2, 7),
122: GmshCellEntry( 44, TRI, 2, 8),
123: GmshCellEntry( 45, TRI, 2, 9),
124: GmshCellEntry( 46, TRI, 2, 10),
126: GmshCellEntry( 3, QUA, 2, 1),
127: GmshCellEntry( 10, QUA, 2, 2),
128: {16, GMSH_QUA, 2, 2, 4, 8, GmshLexOrder_QUA_2_Serendipity},
129: GmshCellEntry( 36, QUA, 2, 3),
130: GmshCellEntry( 37, QUA, 2, 4),
131: GmshCellEntry( 38, QUA, 2, 5),
132: GmshCellEntry( 47, QUA, 2, 6),
133: GmshCellEntry( 48, QUA, 2, 7),
134: GmshCellEntry( 49, QUA, 2, 8),
135: GmshCellEntry( 50, QUA, 2, 9),
136: GmshCellEntry( 51, QUA, 2, 10),
138: GmshCellEntry( 4, TET, 3, 1),
139: GmshCellEntry( 11, TET, 3, 2),
140: GmshCellEntry( 29, TET, 3, 3),
141: GmshCellEntry( 30, TET, 3, 4),
142: GmshCellEntry( 31, TET, 3, 5),
143: GmshCellEntry( 71, TET, 3, 6),
144: GmshCellEntry( 72, TET, 3, 7),
145: GmshCellEntry( 73, TET, 3, 8),
146: GmshCellEntry( 74, TET, 3, 9),
147: GmshCellEntry( 75, TET, 3, 10),
149: GmshCellEntry( 5, HEX, 3, 1),
150: GmshCellEntry( 12, HEX, 3, 2),
151: {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity},
152: GmshCellEntry( 92, HEX, 3, 3),
153: GmshCellEntry( 93, HEX, 3, 4),
154: GmshCellEntry( 94, HEX, 3, 5),
155: GmshCellEntry( 95, HEX, 3, 6),
156: GmshCellEntry( 96, HEX, 3, 7),
157: GmshCellEntry( 97, HEX, 3, 8),
158: GmshCellEntry( 98, HEX, 3, 9),
159: GmshCellEntry( -1, HEX, 3, 10),
161: GmshCellEntry( 6, PRI, 3, 1),
162: GmshCellEntry( 13, PRI, 3, 2),
163: GmshCellEntry( 90, PRI, 3, 3),
164: GmshCellEntry( 91, PRI, 3, 4),
165: GmshCellEntry(106, PRI, 3, 5),
166: GmshCellEntry(107, PRI, 3, 6),
167: GmshCellEntry(108, PRI, 3, 7),
168: GmshCellEntry(109, PRI, 3, 8),
169: GmshCellEntry(110, PRI, 3, 9),
170: GmshCellEntry( -1, PRI, 3, 10),
172: GmshCellEntry( 7, PYR, 3, 1),
173: GmshCellEntry( 14, PYR, 3, 2),
174: GmshCellEntry(118, PYR, 3, 3),
175: GmshCellEntry(119, PYR, 3, 4),
176: GmshCellEntry(120, PYR, 3, 5),
177: GmshCellEntry(121, PYR, 3, 6),
178: GmshCellEntry(122, PYR, 3, 7),
179: GmshCellEntry(123, PYR, 3, 8),
180: GmshCellEntry(124, PYR, 3, 9),
181: GmshCellEntry( -1, PYR, 3, 10)
183: #if 0
184: {20, GMSH_TRI, 2, 3, 3, 9, NULL},
185: {18, GMSH_PRI, 3, 2, 6, 15, NULL},
186: {19, GMSH_PYR, 3, 2, 5, 13, NULL},
187: #endif
188: };
190: static GmshCellInfo GmshCellMap[150];
192: static PetscErrorCode GmshCellInfoSetUp(void)
193: {
194: size_t i, n;
195: static PetscBool called = PETSC_FALSE;
197: if (called) return 0;
198: called = PETSC_TRUE;
199: n = sizeof(GmshCellMap)/sizeof(GmshCellMap[0]);
200: for (i = 0; i < n; ++i) {
201: GmshCellMap[i].cellType = -1;
202: GmshCellMap[i].polytope = -1;
203: }
204: n = sizeof(GmshCellTable)/sizeof(GmshCellTable[0]);
205: for (i = 0; i < n; ++i) {
206: if (GmshCellTable[i].cellType <= 0) continue;
207: GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i];
208: }
209: return 0;
210: }
212: #define GmshCellTypeCheck(ct) PetscMacroReturnStandard( \
213: const int _ct_ = (int)ct; \
217: )
219: typedef struct {
220: PetscViewer viewer;
221: int fileFormat;
222: int dataSize;
223: PetscBool binary;
224: PetscBool byteSwap;
225: size_t wlen;
226: void *wbuf;
227: size_t slen;
228: void *sbuf;
229: PetscInt *nbuf;
230: PetscInt nodeStart;
231: PetscInt nodeEnd;
232: PetscInt *nodeMap;
233: } GmshFile;
235: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
236: {
237: size_t size = count * eltsize;
239: if (gmsh->wlen < size) {
240: PetscFree(gmsh->wbuf);
241: PetscMalloc(size, &gmsh->wbuf);
242: gmsh->wlen = size;
243: }
244: *(void**)buf = size ? gmsh->wbuf : NULL;
245: return 0;
246: }
248: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
249: {
250: size_t dataSize = (size_t)gmsh->dataSize;
251: size_t size = count * dataSize;
253: if (gmsh->slen < size) {
254: PetscFree(gmsh->sbuf);
255: PetscMalloc(size, &gmsh->sbuf);
256: gmsh->slen = size;
257: }
258: *(void**)buf = size ? gmsh->sbuf : NULL;
259: return 0;
260: }
262: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
263: {
264: PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
265: if (gmsh->byteSwap) PetscByteSwap(buf, dtype, count);
266: return 0;
267: }
269: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
270: {
271: PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
272: return 0;
273: }
275: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
276: {
277: PetscStrcmp(line, Section, match);
278: return 0;
279: }
281: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
282: {
283: PetscBool match;
285: GmshMatch(gmsh, Section, line, &match);
287: return 0;
288: }
290: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
291: {
292: PetscBool match;
294: while (PETSC_TRUE) {
295: GmshReadString(gmsh, line, 1);
296: GmshMatch(gmsh, "$Comments", line, &match);
297: if (!match) break;
298: while (PETSC_TRUE) {
299: GmshReadString(gmsh, line, 1);
300: GmshMatch(gmsh, "$EndComments", line, &match);
301: if (match) break;
302: }
303: }
304: return 0;
305: }
307: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
308: {
309: GmshReadString(gmsh, line, 1);
310: GmshExpect(gmsh, EndSection, line);
311: return 0;
312: }
314: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
315: {
316: PetscInt i;
317: size_t dataSize = (size_t)gmsh->dataSize;
319: if (dataSize == sizeof(PetscInt)) {
320: GmshRead(gmsh, buf, count, PETSC_INT);
321: } else if (dataSize == sizeof(int)) {
322: int *ibuf = NULL;
323: GmshBufferSizeGet(gmsh, count, &ibuf);
324: GmshRead(gmsh, ibuf, count, PETSC_ENUM);
325: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
326: } else if (dataSize == sizeof(long)) {
327: long *ibuf = NULL;
328: GmshBufferSizeGet(gmsh, count, &ibuf);
329: GmshRead(gmsh, ibuf, count, PETSC_LONG);
330: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
331: } else if (dataSize == sizeof(PetscInt64)) {
332: PetscInt64 *ibuf = NULL;
333: GmshBufferSizeGet(gmsh, count, &ibuf);
334: GmshRead(gmsh, ibuf, count, PETSC_INT64);
335: for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
336: }
337: return 0;
338: }
340: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
341: {
342: GmshRead(gmsh, buf, count, PETSC_ENUM);
343: return 0;
344: }
346: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
347: {
348: GmshRead(gmsh, buf, count, PETSC_DOUBLE);
349: return 0;
350: }
352: typedef struct {
353: PetscInt id; /* Entity ID */
354: PetscInt dim; /* Dimension */
355: double bbox[6]; /* Bounding box */
356: PetscInt numTags; /* Size of tag array */
357: int tags[4]; /* Tag array */
358: } GmshEntity;
360: typedef struct {
361: GmshEntity *entity[4];
362: PetscHMapI entityMap[4];
363: } GmshEntities;
365: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
366: {
367: PetscInt dim;
369: PetscNew(entities);
370: for (dim = 0; dim < 4; ++dim) {
371: PetscCalloc1(count[dim], &(*entities)->entity[dim]);
372: PetscHMapICreate(&(*entities)->entityMap[dim]);
373: }
374: return 0;
375: }
377: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
378: {
379: PetscInt dim;
381: if (!*entities) return 0;
382: for (dim = 0; dim < 4; ++dim) {
383: PetscFree((*entities)->entity[dim]);
384: PetscHMapIDestroy(&(*entities)->entityMap[dim]);
385: }
386: PetscFree((*entities));
387: return 0;
388: }
390: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity** entity)
391: {
392: PetscHMapISet(entities->entityMap[dim], eid, index);
393: entities->entity[dim][index].dim = dim;
394: entities->entity[dim][index].id = eid;
395: if (entity) *entity = &entities->entity[dim][index];
396: return 0;
397: }
399: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity** entity)
400: {
401: PetscInt index;
403: PetscHMapIGet(entities->entityMap[dim], eid, &index);
404: *entity = &entities->entity[dim][index];
405: return 0;
406: }
408: typedef struct {
409: PetscInt *id; /* Node IDs */
410: double *xyz; /* Coordinates */
411: PetscInt *tag; /* Physical tag */
412: } GmshNodes;
414: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
415: {
416: PetscNew(nodes);
417: PetscMalloc1(count*1, &(*nodes)->id);
418: PetscMalloc1(count*3, &(*nodes)->xyz);
419: PetscMalloc1(count*1, &(*nodes)->tag);
420: return 0;
421: }
423: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
424: {
425: if (!*nodes) return 0;
426: PetscFree((*nodes)->id);
427: PetscFree((*nodes)->xyz);
428: PetscFree((*nodes)->tag);
429: PetscFree((*nodes));
430: return 0;
431: }
433: typedef struct {
434: PetscInt id; /* Element ID */
435: PetscInt dim; /* Dimension */
436: PetscInt cellType; /* Cell type */
437: PetscInt numVerts; /* Size of vertex array */
438: PetscInt numNodes; /* Size of node array */
439: PetscInt *nodes; /* Vertex/Node array */
440: PetscInt numTags; /* Size of physical tag array */
441: int tags[4]; /* Physical tag array */
442: } GmshElement;
444: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
445: {
446: PetscCalloc1(count, elements);
447: return 0;
448: }
450: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
451: {
452: if (!*elements) return 0;
453: PetscFree(*elements);
454: return 0;
455: }
457: typedef struct {
458: PetscInt dim;
459: PetscInt order;
460: GmshEntities *entities;
461: PetscInt numNodes;
462: GmshNodes *nodelist;
463: PetscInt numElems;
464: GmshElement *elements;
465: PetscInt numVerts;
466: PetscInt numCells;
467: PetscInt *periodMap;
468: PetscInt *vertexMap;
469: PetscSegBuffer segbuf;
470: PetscInt numRegions;
471: PetscInt *regionTags;
472: char **regionNames;
473: } GmshMesh;
475: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
476: {
477: PetscNew(mesh);
478: PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);
479: return 0;
480: }
482: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
483: {
484: PetscInt r;
486: if (!*mesh) return 0;
487: GmshEntitiesDestroy(&(*mesh)->entities);
488: GmshNodesDestroy(&(*mesh)->nodelist);
489: GmshElementsDestroy(&(*mesh)->elements);
490: PetscFree((*mesh)->periodMap);
491: PetscFree((*mesh)->vertexMap);
492: PetscSegBufferDestroy(&(*mesh)->segbuf);
493: for (r = 0; r < (*mesh)->numRegions; ++r) PetscFree((*mesh)->regionNames[r]);
494: PetscFree2((*mesh)->regionTags, (*mesh)->regionNames);
495: PetscFree((*mesh));
496: return 0;
497: }
499: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
500: {
501: PetscViewer viewer = gmsh->viewer;
502: PetscBool byteSwap = gmsh->byteSwap;
503: char line[PETSC_MAX_PATH_LEN];
504: int n, num, nid, snum;
505: GmshNodes *nodes;
507: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
508: snum = sscanf(line, "%d", &num);
510: GmshNodesCreate(num, &nodes);
511: mesh->numNodes = num;
512: mesh->nodelist = nodes;
513: for (n = 0; n < num; ++n) {
514: double *xyz = nodes->xyz + n*3;
515: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
516: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
517: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
518: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
519: nodes->id[n] = nid;
520: nodes->tag[n] = -1;
521: }
522: return 0;
523: }
525: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
526: file contents multiple times to figure out the true number of cells and facets
527: in the given mesh. To make this more efficient we read the file contents only
528: once and store them in memory, while determining the true number of cells. */
529: static PetscErrorCode GmshReadElements_v22(GmshFile* gmsh, GmshMesh *mesh)
530: {
531: PetscViewer viewer = gmsh->viewer;
532: PetscBool binary = gmsh->binary;
533: PetscBool byteSwap = gmsh->byteSwap;
534: char line[PETSC_MAX_PATH_LEN];
535: int i, c, p, num, ibuf[1+4+1000], snum;
536: int cellType, numElem, numVerts, numNodes, numTags;
537: GmshElement *elements;
538: PetscInt *nodeMap = gmsh->nodeMap;
540: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
541: snum = sscanf(line, "%d", &num);
543: GmshElementsCreate(num, &elements);
544: mesh->numElems = num;
545: mesh->elements = elements;
546: for (c = 0; c < num;) {
547: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
548: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);
550: cellType = binary ? ibuf[0] : ibuf[1];
551: numElem = binary ? ibuf[1] : 1;
552: numTags = ibuf[2];
554: GmshCellTypeCheck(cellType);
555: numVerts = GmshCellMap[cellType].numVerts;
556: numNodes = GmshCellMap[cellType].numNodes;
558: for (i = 0; i < numElem; ++i, ++c) {
559: GmshElement *element = elements + c;
560: const int off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
561: const int *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
562: PetscViewerRead(viewer, ibuf+off, nint, NULL, PETSC_ENUM);
563: if (byteSwap) PetscByteSwap(ibuf+off, PETSC_ENUM, nint);
564: element->id = id[0];
565: element->dim = GmshCellMap[cellType].dim;
566: element->cellType = cellType;
567: element->numVerts = numVerts;
568: element->numNodes = numNodes;
569: element->numTags = PetscMin(numTags, 4);
570: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
571: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
572: for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
573: }
574: }
575: return 0;
576: }
578: /*
579: $Entities
580: numPoints(unsigned long) numCurves(unsigned long)
581: numSurfaces(unsigned long) numVolumes(unsigned long)
582: // points
583: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
584: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
585: numPhysicals(unsigned long) phyisicalTag[...](int)
586: ...
587: // curves
588: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
589: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
590: numPhysicals(unsigned long) physicalTag[...](int)
591: numBREPVert(unsigned long) tagBREPVert[...](int)
592: ...
593: // surfaces
594: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
595: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
596: numPhysicals(unsigned long) physicalTag[...](int)
597: numBREPCurve(unsigned long) tagBREPCurve[...](int)
598: ...
599: // volumes
600: tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
601: boxMaxX(double) boxMaxY(double) boxMaxZ(double)
602: numPhysicals(unsigned long) physicalTag[...](int)
603: numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
604: ...
605: $EndEntities
606: */
607: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
608: {
609: PetscViewer viewer = gmsh->viewer;
610: PetscBool byteSwap = gmsh->byteSwap;
611: long index, num, lbuf[4];
612: int dim, eid, numTags, *ibuf, t;
613: PetscInt count[4], i;
614: GmshEntity *entity = NULL;
616: PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
617: if (byteSwap) PetscByteSwap(lbuf, PETSC_LONG, 4);
618: for (i = 0; i < 4; ++i) count[i] = lbuf[i];
619: GmshEntitiesCreate(count, &mesh->entities);
620: for (dim = 0; dim < 4; ++dim) {
621: for (index = 0; index < count[dim]; ++index) {
622: PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
623: if (byteSwap) PetscByteSwap(&eid, PETSC_ENUM, 1);
624: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
625: PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
626: if (byteSwap) PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);
627: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
628: if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
629: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
630: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
631: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
632: entity->numTags = numTags = (int) PetscMin(num, 4);
633: for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
634: if (dim == 0) continue;
635: PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
636: if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
637: GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
638: PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
639: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
640: }
641: }
642: return 0;
643: }
645: /*
646: $Nodes
647: numEntityBlocks(unsigned long) numNodes(unsigned long)
648: tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
649: tag(int) x(double) y(double) z(double)
650: ...
651: ...
652: $EndNodes
653: */
654: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
655: {
656: PetscViewer viewer = gmsh->viewer;
657: PetscBool byteSwap = gmsh->byteSwap;
658: long block, node, n, numEntityBlocks, numTotalNodes, numNodes;
659: int info[3], nid;
660: GmshNodes *nodes;
662: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
663: if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
664: PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
665: if (byteSwap) PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);
666: GmshNodesCreate(numTotalNodes, &nodes);
667: mesh->numNodes = numTotalNodes;
668: mesh->nodelist = nodes;
669: for (n = 0, block = 0; block < numEntityBlocks; ++block) {
670: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
671: PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
672: if (byteSwap) PetscByteSwap(&numNodes, PETSC_LONG, 1);
673: if (gmsh->binary) {
674: size_t nbytes = sizeof(int) + 3*sizeof(double);
675: char *cbuf = NULL; /* dummy value to prevent warning from compiler about possible unitilized value */
676: GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
677: PetscViewerRead(viewer, cbuf, numNodes*nbytes, NULL, PETSC_CHAR);
678: for (node = 0; node < numNodes; ++node, ++n) {
679: char *cnid = cbuf + node*nbytes, *cxyz = cnid + sizeof(int);
680: double *xyz = nodes->xyz + n*3;
681: if (!PetscBinaryBigEndian()) PetscByteSwap(cnid, PETSC_ENUM, 1);
682: if (!PetscBinaryBigEndian()) PetscByteSwap(cxyz, PETSC_DOUBLE, 3);
683: PetscMemcpy(&nid, cnid, sizeof(int));
684: PetscMemcpy(xyz, cxyz, 3*sizeof(double));
685: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
686: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
687: nodes->id[n] = nid;
688: nodes->tag[n] = -1;
689: }
690: } else {
691: for (node = 0; node < numNodes; ++node, ++n) {
692: double *xyz = nodes->xyz + n*3;
693: PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
694: PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
695: if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
696: if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
697: nodes->id[n] = nid;
698: nodes->tag[n] = -1;
699: }
700: }
701: }
702: return 0;
703: }
705: /*
706: $Elements
707: numEntityBlocks(unsigned long) numElements(unsigned long)
708: tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
709: tag(int) numVert[...](int)
710: ...
711: ...
712: $EndElements
713: */
714: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
715: {
716: PetscViewer viewer = gmsh->viewer;
717: PetscBool byteSwap = gmsh->byteSwap;
718: long c, block, numEntityBlocks, numTotalElements, elem, numElements;
719: int p, info[3], *ibuf = NULL;
720: int eid, dim, cellType, numVerts, numNodes, numTags;
721: GmshEntity *entity = NULL;
722: GmshElement *elements;
723: PetscInt *nodeMap = gmsh->nodeMap;
725: PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
726: if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
727: PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
728: if (byteSwap) PetscByteSwap(&numTotalElements, PETSC_LONG, 1);
729: GmshElementsCreate(numTotalElements, &elements);
730: mesh->numElems = numTotalElements;
731: mesh->elements = elements;
732: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
733: PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
734: if (byteSwap) PetscByteSwap(info, PETSC_ENUM, 3);
735: eid = info[0]; dim = info[1]; cellType = info[2];
736: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
737: GmshCellTypeCheck(cellType);
738: numVerts = GmshCellMap[cellType].numVerts;
739: numNodes = GmshCellMap[cellType].numNodes;
740: numTags = entity->numTags;
741: PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
742: if (byteSwap) PetscByteSwap(&numElements, PETSC_LONG, 1);
743: GmshBufferGet(gmsh, (1+numNodes)*numElements, sizeof(int), &ibuf);
744: PetscViewerRead(viewer, ibuf, (1+numNodes)*numElements, NULL, PETSC_ENUM);
745: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, (1+numNodes)*numElements);
746: for (elem = 0; elem < numElements; ++elem, ++c) {
747: GmshElement *element = elements + c;
748: const int *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
749: element->id = id[0];
750: element->dim = dim;
751: element->cellType = cellType;
752: element->numVerts = numVerts;
753: element->numNodes = numNodes;
754: element->numTags = numTags;
755: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
756: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
757: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
758: }
759: }
760: return 0;
761: }
763: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
764: {
765: PetscViewer viewer = gmsh->viewer;
766: int fileFormat = gmsh->fileFormat;
767: PetscBool binary = gmsh->binary;
768: PetscBool byteSwap = gmsh->byteSwap;
769: int numPeriodic, snum, i;
770: char line[PETSC_MAX_PATH_LEN];
771: PetscInt *nodeMap = gmsh->nodeMap;
773: if (fileFormat == 22 || !binary) {
774: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
775: snum = sscanf(line, "%d", &numPeriodic);
777: } else {
778: PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
779: if (byteSwap) PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);
780: }
781: for (i = 0; i < numPeriodic; i++) {
782: int ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
783: long j, nNodes;
784: double affine[16];
786: if (fileFormat == 22 || !binary) {
787: PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
788: snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
790: } else {
791: PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
792: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);
793: correspondingDim = ibuf[0]; correspondingTag = ibuf[1]; primaryTag = ibuf[2];
794: }
795: (void)correspondingDim; (void)correspondingTag; (void)primaryTag; /* unused */
797: if (fileFormat == 22 || !binary) {
798: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
799: snum = sscanf(line, "%ld", &nNodes);
800: if (snum != 1) { /* discard transformation and try again */
801: PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
802: PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
803: snum = sscanf(line, "%ld", &nNodes);
805: }
806: } else {
807: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
808: if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
809: if (nNodes == -1) { /* discard transformation and try again */
810: PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
811: PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
812: if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
813: }
814: }
816: for (j = 0; j < nNodes; j++) {
817: if (fileFormat == 22 || !binary) {
818: PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
819: snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
821: } else {
822: PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
823: if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 2);
824: correspondingNode = ibuf[0]; primaryNode = ibuf[1];
825: }
826: correspondingNode = (int) nodeMap[correspondingNode];
827: primaryNode = (int) nodeMap[primaryNode];
828: periodicMap[correspondingNode] = primaryNode;
829: }
830: }
831: return 0;
832: }
834: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
835: $Entities
836: numPoints(size_t) numCurves(size_t)
837: numSurfaces(size_t) numVolumes(size_t)
838: pointTag(int) X(double) Y(double) Z(double)
839: numPhysicalTags(size_t) physicalTag(int) ...
840: ...
841: curveTag(int) minX(double) minY(double) minZ(double)
842: maxX(double) maxY(double) maxZ(double)
843: numPhysicalTags(size_t) physicalTag(int) ...
844: numBoundingPoints(size_t) pointTag(int) ...
845: ...
846: surfaceTag(int) minX(double) minY(double) minZ(double)
847: maxX(double) maxY(double) maxZ(double)
848: numPhysicalTags(size_t) physicalTag(int) ...
849: numBoundingCurves(size_t) curveTag(int) ...
850: ...
851: volumeTag(int) minX(double) minY(double) minZ(double)
852: maxX(double) maxY(double) maxZ(double)
853: numPhysicalTags(size_t) physicalTag(int) ...
854: numBoundngSurfaces(size_t) surfaceTag(int) ...
855: ...
856: $EndEntities
857: */
858: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
859: {
860: PetscInt count[4], index, numTags, i;
861: int dim, eid, *tags = NULL;
862: GmshEntity *entity = NULL;
864: GmshReadSize(gmsh, count, 4);
865: GmshEntitiesCreate(count, &mesh->entities);
866: for (dim = 0; dim < 4; ++dim) {
867: for (index = 0; index < count[dim]; ++index) {
868: GmshReadInt(gmsh, &eid, 1);
869: GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
870: GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
871: GmshReadSize(gmsh, &numTags, 1);
872: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
873: GmshReadInt(gmsh, tags, numTags);
874: entity->numTags = PetscMin(numTags, 4);
875: for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
876: if (dim == 0) continue;
877: GmshReadSize(gmsh, &numTags, 1);
878: GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
879: GmshReadInt(gmsh, tags, numTags);
880: /* Currently, we do not save the ids for the bounding entities */
881: }
882: }
883: return 0;
884: }
886: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
887: $Nodes
888: numEntityBlocks(size_t) numNodes(size_t)
889: minNodeTag(size_t) maxNodeTag(size_t)
890: entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
891: nodeTag(size_t)
892: ...
893: x(double) y(double) z(double)
894: < u(double; if parametric and entityDim = 1 or entityDim = 2) >
895: < v(double; if parametric and entityDim = 2) >
896: ...
897: ...
898: $EndNodes
899: */
900: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
901: {
902: int info[3], dim, eid, parametric;
903: PetscInt sizes[4], numEntityBlocks, numTags, numNodes, numNodesBlock = 0, block, node, n;
904: GmshEntity *entity = NULL;
905: GmshNodes *nodes;
907: GmshReadSize(gmsh, sizes, 4);
908: numEntityBlocks = sizes[0]; numNodes = sizes[1];
909: GmshNodesCreate(numNodes, &nodes);
910: mesh->numNodes = numNodes;
911: mesh->nodelist = nodes;
912: for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
913: GmshReadInt(gmsh, info, 3);
914: dim = info[0]; eid = info[1]; parametric = info[2];
915: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
916: numTags = PetscMin(1, entity->numTags);
917: if (entity->numTags > 1) PetscInfo(NULL, "Entity %d has more than %d physical tags, assigning only the first to nodes", eid, 1);
919: GmshReadSize(gmsh, &numNodesBlock, 1);
920: GmshReadSize(gmsh, nodes->id+node, numNodesBlock);
921: GmshReadDouble(gmsh, nodes->xyz+node*3, numNodesBlock*3);
922: for (n = 0; n < numNodesBlock; ++n) nodes->tag[node+n] = numTags ? entity->tags[0] : -1;
923: }
924: gmsh->nodeStart = sizes[2];
925: gmsh->nodeEnd = sizes[3]+1;
926: return 0;
927: }
929: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
930: $Elements
931: numEntityBlocks(size_t) numElements(size_t)
932: minElementTag(size_t) maxElementTag(size_t)
933: entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
934: elementTag(size_t) nodeTag(size_t) ...
935: ...
936: ...
937: $EndElements
938: */
939: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
940: {
941: int info[3], eid, dim, cellType;
942: PetscInt sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
943: GmshEntity *entity = NULL;
944: GmshElement *elements;
945: PetscInt *nodeMap = gmsh->nodeMap;
947: GmshReadSize(gmsh, sizes, 4);
948: numEntityBlocks = sizes[0]; numElements = sizes[1];
949: GmshElementsCreate(numElements, &elements);
950: mesh->numElems = numElements;
951: mesh->elements = elements;
952: for (c = 0, block = 0; block < numEntityBlocks; ++block) {
953: GmshReadInt(gmsh, info, 3);
954: dim = info[0]; eid = info[1]; cellType = info[2];
955: GmshEntitiesGet(mesh->entities, dim, eid, &entity);
956: GmshCellTypeCheck(cellType);
957: numVerts = GmshCellMap[cellType].numVerts;
958: numNodes = GmshCellMap[cellType].numNodes;
959: numTags = PetscMin(4, entity->numTags);
960: if (entity->numTags > 4) PetscInfo(NULL, "Entity %d has more then %d physical tags, assigning only the first to elements", eid, 4);
961: GmshReadSize(gmsh, &numBlockElements, 1);
962: GmshBufferGet(gmsh, (1+numNodes)*numBlockElements, sizeof(PetscInt), &ibuf);
963: GmshReadSize(gmsh, ibuf, (1+numNodes)*numBlockElements);
964: for (elem = 0; elem < numBlockElements; ++elem, ++c) {
965: GmshElement *element = elements + c;
966: const PetscInt *id = ibuf + elem*(1+numNodes), *nodes = id + 1;
967: element->id = id[0];
968: element->dim = dim;
969: element->cellType = cellType;
970: element->numVerts = numVerts;
971: element->numNodes = numNodes;
972: element->numTags = numTags;
973: PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
974: for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
975: for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
976: }
977: }
978: return 0;
979: }
981: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
982: $Periodic
983: numPeriodicLinks(size_t)
984: entityDim(int) entityTag(int) entityTagPrimary(int)
985: numAffine(size_t) value(double) ...
986: numCorrespondingNodes(size_t)
987: nodeTag(size_t) nodeTagPrimary(size_t)
988: ...
989: ...
990: $EndPeriodic
991: */
992: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
993: {
994: int info[3];
995: double dbuf[16];
996: PetscInt numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
997: PetscInt *nodeMap = gmsh->nodeMap;
999: GmshReadSize(gmsh, &numPeriodicLinks, 1);
1000: for (link = 0; link < numPeriodicLinks; ++link) {
1001: GmshReadInt(gmsh, info, 3);
1002: GmshReadSize(gmsh, &numAffine, 1);
1003: GmshReadDouble(gmsh, dbuf, numAffine);
1004: GmshReadSize(gmsh, &numCorrespondingNodes, 1);
1005: GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
1006: GmshReadSize(gmsh, nodeTags, numCorrespondingNodes*2);
1007: for (node = 0; node < numCorrespondingNodes; ++node) {
1008: PetscInt correspondingNode = nodeMap[nodeTags[node*2+0]];
1009: PetscInt primaryNode = nodeMap[nodeTags[node*2+1]];
1010: periodicMap[correspondingNode] = primaryNode;
1011: }
1012: }
1013: return 0;
1014: }
1016: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1017: $MeshFormat // same as MSH version 2
1018: version(ASCII double; currently 4.1)
1019: file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1020: data-size(ASCII int; sizeof(size_t))
1021: < int with value one; only in binary mode, to detect endianness >
1022: $EndMeshFormat
1023: */
1024: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1025: {
1026: char line[PETSC_MAX_PATH_LEN];
1027: int snum, fileType, fileFormat, dataSize, checkEndian;
1028: float version;
1030: GmshReadString(gmsh, line, 3);
1031: snum = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1038: fileFormat = (int)roundf(version*10);
1041: gmsh->fileFormat = fileFormat;
1042: gmsh->dataSize = dataSize;
1043: gmsh->byteSwap = PETSC_FALSE;
1044: if (gmsh->binary) {
1045: GmshReadInt(gmsh, &checkEndian, 1);
1046: if (checkEndian != 1) {
1047: PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
1049: gmsh->byteSwap = PETSC_TRUE;
1050: }
1051: }
1052: return 0;
1053: }
1055: /*
1056: PhysicalNames
1057: numPhysicalNames(ASCII int)
1058: dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1059: ...
1060: $EndPhysicalNames
1061: */
1062: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1063: {
1064: char line[PETSC_MAX_PATH_LEN], name[128+2], *p, *q;
1065: int snum, region, dim, tag;
1067: GmshReadString(gmsh, line, 1);
1068: snum = sscanf(line, "%d", ®ion);
1069: mesh->numRegions = region;
1071: PetscMalloc2(mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames);
1072: for (region = 0; region < mesh->numRegions; ++region) {
1073: GmshReadString(gmsh, line, 2);
1074: snum = sscanf(line, "%d %d", &dim, &tag);
1076: GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
1077: PetscStrchr(line, '"', &p);
1079: PetscStrrchr(line, '"', &q);
1081: PetscStrncpy(name, p+1, (size_t)(q-p-1));
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: GmshReadEntities_v41(gmsh, mesh); break;
1092: default: GmshReadEntities_v40(gmsh, mesh); break;
1093: }
1094: return 0;
1095: }
1097: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1098: {
1099: switch (gmsh->fileFormat) {
1100: case 41: GmshReadNodes_v41(gmsh, mesh); break;
1101: case 40: GmshReadNodes_v40(gmsh, mesh); break;
1102: default: GmshReadNodes_v22(gmsh, mesh); break;
1103: }
1105: { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1106: if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1107: PetscInt tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1108: GmshNodes *nodes = mesh->nodelist;
1109: for (n = 0; n < mesh->numNodes; ++n) {
1110: const PetscInt tag = nodes->id[n];
1111: tagMin = PetscMin(tag, tagMin);
1112: tagMax = PetscMax(tag, tagMax);
1113: }
1114: gmsh->nodeStart = tagMin;
1115: gmsh->nodeEnd = tagMax+1;
1116: }
1117: }
1119: { /* Support for sparse node tags */
1120: PetscInt n, t;
1121: GmshNodes *nodes = mesh->nodelist;
1122: PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);
1123: for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1124: gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1125: for (n = 0; n < mesh->numNodes; ++n) {
1126: const PetscInt tag = nodes->id[n];
1128: gmsh->nodeMap[tag] = n;
1129: }
1130: }
1131: return 0;
1132: }
1134: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1135: {
1136: switch (gmsh->fileFormat) {
1137: case 41: GmshReadElements_v41(gmsh, mesh); break;
1138: case 40: GmshReadElements_v40(gmsh, mesh); break;
1139: default: GmshReadElements_v22(gmsh, mesh); break;
1140: }
1142: { /* Reorder elements by codimension and polytope type */
1143: PetscInt ne = mesh->numElems;
1144: GmshElement *elements = mesh->elements;
1145: PetscInt keymap[GMSH_NUM_POLYTOPES], nk = 0;
1146: PetscInt offset[GMSH_NUM_POLYTOPES+1], e, k;
1148: for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1149: PetscMemzero(offset,sizeof(offset));
1151: keymap[GMSH_TET] = nk++;
1152: keymap[GMSH_HEX] = nk++;
1153: keymap[GMSH_PRI] = nk++;
1154: keymap[GMSH_PYR] = nk++;
1155: keymap[GMSH_TRI] = nk++;
1156: keymap[GMSH_QUA] = nk++;
1157: keymap[GMSH_SEG] = nk++;
1158: keymap[GMSH_VTX] = nk++;
1160: GmshElementsCreate(mesh->numElems, &mesh->elements);
1161: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1162: for (e = 0; e < ne; ++e) offset[1+key(e)]++;
1163: for (k = 1; k < nk; ++k) offset[k] += offset[k-1];
1164: for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1165: #undef key
1166: GmshElementsDestroy(&elements);
1167: }
1169: { /* Mesh dimension and order */
1170: GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1171: mesh->dim = elem ? GmshCellMap[elem->cellType].dim : 0;
1172: mesh->order = elem ? GmshCellMap[elem->cellType].order : 0;
1173: }
1175: {
1176: PetscBT vtx;
1177: PetscInt dim = mesh->dim, e, n, v;
1179: PetscBTCreate(mesh->numNodes, &vtx);
1181: /* Compute number of cells and set of vertices */
1182: mesh->numCells = 0;
1183: for (e = 0; e < mesh->numElems; ++e) {
1184: GmshElement *elem = mesh->elements + e;
1185: if (elem->dim == dim && dim > 0) mesh->numCells++;
1186: for (v = 0; v < elem->numVerts; v++) {
1187: PetscBTSet(vtx, elem->nodes[v]);
1188: }
1189: }
1191: /* Compute numbering for vertices */
1192: mesh->numVerts = 0;
1193: PetscMalloc1(mesh->numNodes, &mesh->vertexMap);
1194: for (n = 0; n < mesh->numNodes; ++n)
1195: mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;
1197: PetscBTDestroy(&vtx);
1198: }
1199: return 0;
1200: }
1202: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1203: {
1204: PetscInt n;
1206: PetscMalloc1(mesh->numNodes, &mesh->periodMap);
1207: for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1208: switch (gmsh->fileFormat) {
1209: case 41: GmshReadPeriodic_v41(gmsh, mesh->periodMap); break;
1210: default: GmshReadPeriodic_v40(gmsh, mesh->periodMap); break;
1211: }
1213: /* Find canonical primary nodes */
1214: for (n = 0; n < mesh->numNodes; ++n)
1215: while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]])
1216: mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];
1218: /* Renumber vertices (filter out correspondings) */
1219: mesh->numVerts = 0;
1220: for (n = 0; n < mesh->numNodes; ++n)
1221: if (mesh->vertexMap[n] >= 0) /* is vertex */
1222: if (mesh->periodMap[n] == n) /* is primary */
1223: mesh->vertexMap[n] = mesh->numVerts++;
1224: for (n = 0; n < mesh->numNodes; ++n)
1225: if (mesh->vertexMap[n] >= 0) /* is vertex */
1226: if (mesh->periodMap[n] != n) /* is corresponding */
1227: mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1228: return 0;
1229: }
1231: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1232: #define DM_POLYTOPE_PYRAMID DM_POLYTOPE_UNKNOWN
1233: static const DMPolytopeType DMPolytopeMap[] = {
1234: /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1235: /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1236: /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1237: /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1238: /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1239: /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1240: /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1241: /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,
1242: DM_POLYTOPE_UNKNOWN
1243: };
1245: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1246: {
1247: return DMPolytopeMap[GmshCellMap[cellType].polytope];
1248: }
1250: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1251: {
1252: DM K;
1253: PetscSpace P;
1254: PetscDualSpace Q;
1255: PetscQuadrature q, fq;
1256: PetscBool isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1257: PetscBool endpoint = PETSC_TRUE;
1258: char name[32];
1260: /* Create space */
1261: PetscSpaceCreate(comm, &P);
1262: PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1263: PetscSpacePolynomialSetTensor(P, isTensor);
1264: PetscSpaceSetNumComponents(P, Nc);
1265: PetscSpaceSetNumVariables(P, dim);
1266: PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1267: if (prefix) {
1268: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
1269: PetscSpaceSetFromOptions(P);
1270: PetscObjectSetOptionsPrefix((PetscObject) P, NULL);
1271: PetscSpaceGetDegree(P, &k, NULL);
1272: }
1273: PetscSpaceSetUp(P);
1274: /* Create dual space */
1275: PetscDualSpaceCreate(comm, &Q);
1276: PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1277: PetscDualSpaceLagrangeSetTensor(Q, isTensor);
1278: PetscDualSpaceLagrangeSetContinuity(Q, continuity);
1279: PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
1280: PetscDualSpaceSetNumComponents(Q, Nc);
1281: PetscDualSpaceSetOrder(Q, k);
1282: DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
1283: PetscDualSpaceSetDM(Q, K);
1284: DMDestroy(&K);
1285: if (prefix) {
1286: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
1287: PetscDualSpaceSetFromOptions(Q);
1288: PetscObjectSetOptionsPrefix((PetscObject) Q, NULL);
1289: }
1290: PetscDualSpaceSetUp(Q);
1291: /* Create quadrature */
1292: if (isSimplex) {
1293: PetscDTStroudConicalQuadrature(dim, 1, k+1, -1, +1, &q);
1294: PetscDTStroudConicalQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1295: } else {
1296: PetscDTGaussTensorQuadrature(dim, 1, k+1, -1, +1, &q);
1297: PetscDTGaussTensorQuadrature(dim-1, 1, k+1, -1, +1, &fq);
1298: }
1299: /* Create finite element */
1300: PetscFECreate(comm, fem);
1301: PetscFESetType(*fem, PETSCFEBASIC);
1302: PetscFESetNumComponents(*fem, Nc);
1303: PetscFESetBasisSpace(*fem, P);
1304: PetscFESetDualSpace(*fem, Q);
1305: PetscFESetQuadrature(*fem, q);
1306: PetscFESetFaceQuadrature(*fem, fq);
1307: if (prefix) {
1308: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
1309: PetscFESetFromOptions(*fem);
1310: PetscObjectSetOptionsPrefix((PetscObject) *fem, NULL);
1311: }
1312: PetscFESetUp(*fem);
1313: /* Cleanup */
1314: PetscSpaceDestroy(&P);
1315: PetscDualSpaceDestroy(&Q);
1316: PetscQuadratureDestroy(&q);
1317: PetscQuadratureDestroy(&fq);
1318: /* Set finite element name */
1319: PetscSNPrintf(name, sizeof(name), "%s%D", isSimplex? "P" : "Q", k);
1320: PetscFESetName(*fem, name);
1321: return 0;
1322: }
1324: /*@C
1325: DMPlexCreateGmshFromFile - Create a DMPlex mesh from a Gmsh file
1327: + comm - The MPI communicator
1328: . filename - Name of the Gmsh file
1329: - interpolate - Create faces and edges in the mesh
1331: Output Parameter:
1332: . dm - The DM object representing the mesh
1334: Level: beginner
1336: .seealso: DMPlexCreateFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
1337: @*/
1338: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1339: {
1340: PetscViewer viewer;
1341: PetscMPIInt rank;
1342: int fileType;
1343: PetscViewerType vtype;
1345: MPI_Comm_rank(comm, &rank);
1347: /* Determine Gmsh file type (ASCII or binary) from file header */
1348: if (rank == 0) {
1349: GmshFile gmsh[1];
1350: char line[PETSC_MAX_PATH_LEN];
1351: int snum;
1352: float version;
1354: PetscArrayzero(gmsh,1);
1355: PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1356: PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1357: PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1358: PetscViewerFileSetName(gmsh->viewer, filename);
1359: /* Read only the first two lines of the Gmsh file */
1360: GmshReadSection(gmsh, line);
1361: GmshExpect(gmsh, "$MeshFormat", line);
1362: GmshReadString(gmsh, line, 2);
1363: snum = sscanf(line, "%f %d", &version, &fileType);
1368: PetscViewerDestroy(&gmsh->viewer);
1369: }
1370: MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1371: vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;
1373: /* Create appropriate viewer and build plex */
1374: PetscViewerCreate(comm, &viewer);
1375: PetscViewerSetType(viewer, vtype);
1376: PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1377: PetscViewerFileSetName(viewer, filename);
1378: DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1379: PetscViewerDestroy(&viewer);
1380: return 0;
1381: }
1383: /*@
1384: DMPlexCreateGmsh - Create a DMPlex mesh from a Gmsh file viewer
1386: Collective
1388: Input Parameters:
1389: + comm - The MPI communicator
1390: . viewer - The Viewer associated with a Gmsh file
1391: - interpolate - Create faces and edges in the mesh
1393: Output Parameter:
1394: . dm - The DM object representing the mesh
1396: Note: http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
1398: Level: beginner
1400: .seealso: DMPLEX, DMCreate()
1401: @*/
1402: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1403: {
1404: GmshMesh *mesh = NULL;
1405: PetscViewer parentviewer = NULL;
1406: PetscBT periodicVerts = NULL;
1407: PetscBT periodicCells = NULL;
1408: DM cdm;
1409: PetscSection coordSection;
1410: Vec coordinates;
1411: DMLabel cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1412: PetscInt dim = 0, coordDim = -1, order = 0;
1413: PetscInt numNodes = 0, numElems = 0, numVerts = 0, numCells = 0;
1414: PetscInt cell, cone[8], e, n, v, d;
1415: PetscBool binary, usemarker = PETSC_FALSE, useregions = PETSC_FALSE, markvertices = PETSC_FALSE;
1416: PetscBool hybrid = interpolate, periodic = PETSC_TRUE;
1417: PetscBool highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1418: PetscBool isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1419: PetscMPIInt rank;
1422: MPI_Comm_rank(comm, &rank);
1423: PetscObjectOptionsBegin((PetscObject)viewer);
1424: PetscOptionsHead(PetscOptionsObject,"DMPlex Gmsh options");
1425: PetscOptionsBool("-dm_plex_gmsh_hybrid", "Generate hybrid cell bounds", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);
1426: PetscOptionsBool("-dm_plex_gmsh_periodic","Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1427: PetscOptionsBool("-dm_plex_gmsh_highorder","Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);
1428: PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);
1429: PetscOptionsBool("-dm_plex_gmsh_use_marker", "Generate marker label", "DMPlexCreateGmsh", usemarker, &usemarker, NULL);
1430: PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL);
1431: PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL);
1432: PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);
1433: PetscOptionsTail();
1434: PetscOptionsEnd();
1436: GmshCellInfoSetUp();
1438: DMCreate(comm, dm);
1439: DMSetType(*dm, DMPLEX);
1440: PetscLogEventBegin(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1442: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);
1444: /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1445: if (binary) {
1446: parentviewer = viewer;
1447: PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1448: }
1450: if (rank == 0) {
1451: GmshFile gmsh[1];
1452: char line[PETSC_MAX_PATH_LEN];
1453: PetscBool match;
1455: PetscArrayzero(gmsh,1);
1456: gmsh->viewer = viewer;
1457: gmsh->binary = binary;
1459: GmshMeshCreate(&mesh);
1461: /* Read mesh format */
1462: GmshReadSection(gmsh, line);
1463: GmshExpect(gmsh, "$MeshFormat", line);
1464: GmshReadMeshFormat(gmsh);
1465: GmshReadEndSection(gmsh, "$EndMeshFormat", line);
1467: /* OPTIONAL Read physical names */
1468: GmshReadSection(gmsh, line);
1469: GmshMatch(gmsh, "$PhysicalNames", line, &match);
1470: if (match) {
1471: GmshExpect(gmsh, "$PhysicalNames", line);
1472: GmshReadPhysicalNames(gmsh, mesh);
1473: GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1474: /* Initial read for entity section */
1475: GmshReadSection(gmsh, line);
1476: }
1478: /* Read entities */
1479: if (gmsh->fileFormat >= 40) {
1480: GmshExpect(gmsh, "$Entities", line);
1481: GmshReadEntities(gmsh, mesh);
1482: GmshReadEndSection(gmsh, "$EndEntities", line);
1483: /* Initial read for nodes section */
1484: GmshReadSection(gmsh, line);
1485: }
1487: /* Read nodes */
1488: GmshExpect(gmsh, "$Nodes", line);
1489: GmshReadNodes(gmsh, mesh);
1490: GmshReadEndSection(gmsh, "$EndNodes", line);
1492: /* Read elements */
1493: GmshReadSection(gmsh, line);
1494: GmshExpect(gmsh, "$Elements", line);
1495: GmshReadElements(gmsh, mesh);
1496: GmshReadEndSection(gmsh, "$EndElements", line);
1498: /* Read periodic section (OPTIONAL) */
1499: if (periodic) {
1500: GmshReadSection(gmsh, line);
1501: GmshMatch(gmsh, "$Periodic", line, &periodic);
1502: }
1503: if (periodic) {
1504: GmshExpect(gmsh, "$Periodic", line);
1505: GmshReadPeriodic(gmsh, mesh);
1506: GmshReadEndSection(gmsh, "$EndPeriodic", line);
1507: }
1509: PetscFree(gmsh->wbuf);
1510: PetscFree(gmsh->sbuf);
1511: PetscFree(gmsh->nbuf);
1513: dim = mesh->dim;
1514: order = mesh->order;
1515: numNodes = mesh->numNodes;
1516: numElems = mesh->numElems;
1517: numVerts = mesh->numVerts;
1518: numCells = mesh->numCells;
1520: {
1521: GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1522: GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1523: int ptA = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1524: int ptB = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1525: isSimplex = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1526: isHybrid = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1527: hasTetra = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1528: }
1529: }
1531: if (parentviewer) {
1532: PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1533: }
1535: {
1536: int buf[6];
1538: buf[0] = (int)dim;
1539: buf[1] = (int)order;
1540: buf[2] = periodic;
1541: buf[3] = isSimplex;
1542: buf[4] = isHybrid;
1543: buf[5] = hasTetra;
1545: MPI_Bcast(buf, 6, MPI_INT, 0, comm);
1547: dim = buf[0];
1548: order = buf[1];
1549: periodic = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1550: isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1551: isHybrid = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1552: hasTetra = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1553: }
1555: if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1558: /* We do not want this label automatically computed, instead we fill it here */
1559: DMCreateLabel(*dm, "celltype");
1561: /* Allocate the cell-vertex mesh */
1562: DMPlexSetChart(*dm, 0, numCells+numVerts);
1563: for (cell = 0; cell < numCells; ++cell) {
1564: GmshElement *elem = mesh->elements + cell;
1565: DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1566: if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1567: DMPlexSetConeSize(*dm, cell, elem->numVerts);
1568: DMPlexSetCellType(*dm, cell, ctype);
1569: }
1570: for (v = numCells; v < numCells+numVerts; ++v) {
1571: DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1572: }
1573: DMSetUp(*dm);
1575: /* Add cell-vertex connections */
1576: for (cell = 0; cell < numCells; ++cell) {
1577: GmshElement *elem = mesh->elements + cell;
1578: for (v = 0; v < elem->numVerts; ++v) {
1579: const PetscInt nn = elem->nodes[v];
1580: const PetscInt vv = mesh->vertexMap[nn];
1581: cone[v] = numCells + vv;
1582: }
1583: DMPlexReorderCell(*dm, cell, cone);
1584: DMPlexSetCone(*dm, cell, cone);
1585: }
1587: DMSetDimension(*dm, dim);
1588: DMPlexSymmetrize(*dm);
1589: DMPlexStratify(*dm);
1590: if (interpolate) {
1591: DM idm;
1593: DMPlexInterpolate(*dm, &idm);
1594: DMDestroy(dm);
1595: *dm = idm;
1596: }
1598: /* Create the label "marker" over the whole boundary */
1600: if (rank == 0 && usemarker) {
1601: PetscInt f, fStart, fEnd;
1603: DMCreateLabel(*dm, "marker");
1604: DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);
1605: for (f = fStart; f < fEnd; ++f) {
1606: PetscInt suppSize;
1608: DMPlexGetSupportSize(*dm, f, &suppSize);
1609: if (suppSize == 1) {
1610: PetscInt *cone = NULL, coneSize, p;
1612: DMPlexGetTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1613: for (p = 0; p < coneSize; p += 2) {
1614: DMSetLabelValue_Fast(*dm, &marker, "marker", cone[p], 1);
1615: }
1616: DMPlexRestoreTransitiveClosure(*dm, f, PETSC_TRUE, &coneSize, &cone);
1617: }
1618: }
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: const PetscInt tag = elem->tags[0];
1634: PetscInt r;
1636: if (!Nr) DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag);
1637: for (r = 0; r < Nr; ++r) {
1638: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], cell, tag);
1639: }
1640: }
1641: cell++;
1642: }
1644: /* Create face sets */
1645: if (interpolate && elem->dim == dim-1) {
1646: PetscInt joinSize;
1647: const PetscInt *join = NULL;
1648: const PetscInt tag = elem->tags[0];
1649: PetscInt r;
1651: /* Find the relevant facet with vertex joins */
1652: for (v = 0; v < elem->numVerts; ++v) {
1653: const PetscInt nn = elem->nodes[v];
1654: const PetscInt vv = mesh->vertexMap[nn];
1655: cone[v] = vStart + vv;
1656: }
1657: DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1659: if (!Nr) DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag);
1660: for (r = 0; r < Nr; ++r) {
1661: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], join[0], tag);
1662: }
1663: DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1664: }
1666: /* Create vertex sets */
1667: if (elem->dim == 0) {
1668: if (elem->numTags > 0) {
1669: const PetscInt nn = elem->nodes[0];
1670: const PetscInt vv = mesh->vertexMap[nn];
1671: const PetscInt tag = elem->tags[0];
1672: PetscInt r;
1674: if (!Nr) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1675: for (r = 0; r < Nr; ++r) {
1676: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag);
1677: }
1678: }
1679: }
1680: }
1681: if (markvertices) {
1682: for (v = 0; v < numNodes; ++v) {
1683: const PetscInt vv = mesh->vertexMap[v];
1684: const PetscInt tag = mesh->nodelist->tag[v];
1685: PetscInt r;
1687: if (tag != -1) {
1688: if (!Nr) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1689: for (r = 0; r < Nr; ++r) {
1690: if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, ®ionSets[r], mesh->regionNames[r], vStart + vv, tag);
1691: }
1692: }
1693: }
1694: }
1695: PetscFree(regionSets);
1696: }
1698: { /* Create Cell/Face/Vertex Sets labels at all processes */
1699: enum {n = 4};
1700: PetscBool flag[n];
1702: flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1703: flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1704: flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1705: flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1706: MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1707: if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1708: if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1709: if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1710: if (flag[3]) DMCreateLabel(*dm, "marker");
1711: }
1713: if (periodic) {
1714: PetscBTCreate(numVerts, &periodicVerts);
1715: for (n = 0; n < numNodes; ++n) {
1716: if (mesh->vertexMap[n] >= 0) {
1717: if (PetscUnlikely(mesh->periodMap[n] != n)) {
1718: PetscInt m = mesh->periodMap[n];
1719: PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1720: PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1721: }
1722: }
1723: }
1724: PetscBTCreate(numCells, &periodicCells);
1725: for (cell = 0; cell < numCells; ++cell) {
1726: GmshElement *elem = mesh->elements + cell;
1727: for (v = 0; v < elem->numVerts; ++v) {
1728: PetscInt nn = elem->nodes[v];
1729: PetscInt vv = mesh->vertexMap[nn];
1730: if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1731: PetscBTSet(periodicCells, cell); break;
1732: }
1733: }
1734: }
1735: }
1737: /* Setup coordinate DM */
1738: if (coordDim < 0) coordDim = dim;
1739: DMSetCoordinateDim(*dm, coordDim);
1740: DMGetCoordinateDM(*dm, &cdm);
1741: if (highOrder) {
1742: PetscFE fe;
1743: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1744: PetscDTNodeType nodeType = PETSCDTNODES_EQUISPACED;
1746: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1748: GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1749: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1750: DMSetField(cdm, 0, NULL, (PetscObject) fe);
1751: PetscFEDestroy(&fe);
1752: DMCreateDS(cdm);
1753: }
1755: /* Create coordinates */
1756: if (highOrder) {
1758: PetscInt maxDof = GmshNumNodes_HEX(order)*coordDim;
1759: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1760: PetscSection section;
1761: PetscScalar *cellCoords;
1763: DMSetLocalSection(cdm, NULL);
1764: DMGetLocalSection(cdm, &coordSection);
1765: PetscSectionClone(coordSection, §ion);
1766: DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */
1768: DMCreateLocalVector(cdm, &coordinates);
1769: PetscMalloc1(maxDof, &cellCoords);
1770: for (cell = 0; cell < numCells; ++cell) {
1771: GmshElement *elem = mesh->elements + cell;
1772: const int *lexorder = GmshCellMap[elem->cellType].lexorder();
1773: int s = 0;
1774: for (n = 0; n < elem->numNodes; ++n) {
1775: while (lexorder[n+s] < 0) ++s;
1776: const PetscInt node = elem->nodes[lexorder[n+s]];
1777: for (d = 0; d < coordDim; ++d) cellCoords[(n+s)*coordDim+d] = (PetscReal) coords[node*3+d];
1778: }
1779: if (s) {
1780: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1781: PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1782: /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1783: PetscReal hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25,
1784: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1785: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1786: PetscReal hexFrontWeights[27] = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1787: 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1788: -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1789: PetscReal hexLeftWeights[27] = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0,
1790: 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0,
1791: -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0};
1792: PetscReal hexRightWeights[27] = { 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25,
1793: 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5,
1794: 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25};
1795: PetscReal hexBackWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25,
1796: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5,
1797: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25};
1798: PetscReal hexTopWeights[27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1799: 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1800: -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1801: PetscReal hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25,
1802: 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25,
1803: -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25};
1804: PetscReal *sdWeights2[9] = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
1805: PetscReal *sdWeights3[27] = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL,
1806: NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
1807: NULL, NULL, NULL, NULL, hexTopWeights, NULL, NULL, NULL, NULL};
1808: PetscReal **sdWeights[4] = {NULL, NULL, sdWeights2, sdWeights3};
1810: /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
1811: for (n = 0; n < elem->numNodes+s; ++n) {
1812: if (lexorder[n] >= 0) continue;
1813: for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] = 0.0;
1814: for (int bn = 0; bn < elem->numNodes+s; ++bn) {
1815: if (lexorder[bn] < 0) continue;
1816: const PetscReal *weights = sdWeights[coordDim][n];
1817: const PetscInt bnode = elem->nodes[lexorder[bn]];
1818: for (d = 0; d < coordDim; ++d) cellCoords[n*coordDim+d] += weights[bn] * (PetscReal) coords[bnode*3+d];
1819: }
1820: }
1821: }
1822: DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1823: }
1824: PetscSectionDestroy(§ion);
1825: PetscFree(cellCoords);
1827: } else {
1829: PetscInt *nodeMap;
1830: double *coords = mesh ? mesh->nodelist->xyz : NULL;
1831: PetscScalar *pointCoords;
1833: DMGetLocalSection(cdm, &coordSection);
1834: PetscSectionSetNumFields(coordSection, 1);
1835: PetscSectionSetFieldComponents(coordSection, 0, coordDim);
1836: if (periodic) { /* we need to localize coordinates on cells */
1837: PetscSectionSetChart(coordSection, 0, numCells+numVerts);
1838: } else {
1839: PetscSectionSetChart(coordSection, numCells, numCells+numVerts);
1840: }
1841: if (periodic) {
1842: for (cell = 0; cell < numCells; ++cell) {
1843: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1844: GmshElement *elem = mesh->elements + cell;
1845: PetscInt dof = elem->numVerts * coordDim;
1846: PetscSectionSetDof(coordSection, cell, dof);
1847: PetscSectionSetFieldDof(coordSection, cell, 0, dof);
1848: }
1849: }
1850: }
1851: for (v = numCells; v < numCells+numVerts; ++v) {
1852: PetscSectionSetDof(coordSection, v, coordDim);
1853: PetscSectionSetFieldDof(coordSection, v, 0, coordDim);
1854: }
1855: PetscSectionSetUp(coordSection);
1857: DMCreateLocalVector(cdm, &coordinates);
1858: VecGetArray(coordinates, &pointCoords);
1859: if (periodic) {
1860: for (cell = 0; cell < numCells; ++cell) {
1861: if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1862: GmshElement *elem = mesh->elements + cell;
1863: PetscInt off, node;
1864: for (v = 0; v < elem->numVerts; ++v)
1865: cone[v] = elem->nodes[v];
1866: DMPlexReorderCell(cdm, cell, cone);
1867: PetscSectionGetOffset(coordSection, cell, &off);
1868: for (v = 0; v < elem->numVerts; ++v)
1869: for (node = cone[v], d = 0; d < coordDim; ++d)
1870: pointCoords[off++] = (PetscReal) coords[node*3+d];
1871: }
1872: }
1873: }
1874: PetscMalloc1(numVerts, &nodeMap);
1875: for (n = 0; n < numNodes; n++)
1876: if (mesh->vertexMap[n] >= 0)
1877: nodeMap[mesh->vertexMap[n]] = n;
1878: for (v = 0; v < numVerts; ++v) {
1879: PetscInt off, node = nodeMap[v];
1880: PetscSectionGetOffset(coordSection, numCells + v, &off);
1881: for (d = 0; d < coordDim; ++d)
1882: pointCoords[off+d] = (PetscReal) coords[node*3+d];
1883: }
1884: PetscFree(nodeMap);
1885: VecRestoreArray(coordinates, &pointCoords);
1887: }
1889: PetscObjectSetName((PetscObject) coordinates, "coordinates");
1890: VecSetBlockSize(coordinates, coordDim);
1891: DMSetCoordinatesLocal(*dm, coordinates);
1892: VecDestroy(&coordinates);
1893: DMSetPeriodicity(*dm, periodic, NULL, NULL, NULL);
1895: GmshMeshDestroy(&mesh);
1896: PetscBTDestroy(&periodicVerts);
1897: PetscBTDestroy(&periodicCells);
1899: if (highOrder && project) {
1900: PetscFE fe;
1901: const char prefix[] = "dm_plex_gmsh_project_";
1902: PetscBool continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1903: PetscDTNodeType nodeType = PETSCDTNODES_GAUSSJACOBI;
1905: if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
1907: GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1908: PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1909: DMProjectCoordinates(*dm, fe);
1910: PetscFEDestroy(&fe);
1911: }
1913: PetscLogEventEnd(DMPLEX_CreateGmsh,*dm,NULL,NULL,NULL);
1914: return 0;
1915: }