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", &region);
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, &regionSets);
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, &regionSets[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, &regionSets[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, &regionSets[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, &regionSets[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, &section);
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(&section);
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: }