Actual source code: plexgmsh.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>
  3: #include <petsc/private/hashmapi.h>

  5: #include <../src/dm/impls/plex/gmshlex.h>

  7: #define GMSH_LEXORDER_ITEM(T, p) \
  8:   static int *GmshLexOrder_##T##_##p(void) \
  9:   { \
 10:     static int Gmsh_LexOrder_##T##_##p[GmshNumNodes_##T(p)] = {-1}; \
 11:     int       *lex                                          = Gmsh_LexOrder_##T##_##p; \
 12:     if (lex[0] == -1) (void)GmshLexOrder_##T(p, lex, 0); \
 13:     return lex; \
 14:   }

 16: static int *GmshLexOrder_QUA_2_Serendipity(void)
 17: {
 18:   static int Gmsh_LexOrder_QUA_2_Serendipity[9] = {-1};
 19:   int       *lex                                = Gmsh_LexOrder_QUA_2_Serendipity;
 20:   if (lex[0] == -1) {
 21:     /* Vertices */
 22:     lex[0] = 0;
 23:     lex[2] = 1;
 24:     lex[8] = 2;
 25:     lex[6] = 3;
 26:     /* Edges */
 27:     lex[1] = 4;
 28:     lex[5] = 5;
 29:     lex[7] = 6;
 30:     lex[3] = 7;
 31:     /* Cell */
 32:     lex[4] = -8;
 33:   }
 34:   return lex;
 35: }

 37: static int *GmshLexOrder_HEX_2_Serendipity(void)
 38: {
 39:   static int Gmsh_LexOrder_HEX_2_Serendipity[27] = {-1};
 40:   int       *lex                                 = Gmsh_LexOrder_HEX_2_Serendipity;
 41:   if (lex[0] == -1) {
 42:     /* Vertices */
 43:     lex[0]  = 0;
 44:     lex[2]  = 1;
 45:     lex[8]  = 2;
 46:     lex[6]  = 3;
 47:     lex[18] = 4;
 48:     lex[20] = 5;
 49:     lex[26] = 6;
 50:     lex[24] = 7;
 51:     /* Edges */
 52:     lex[1]  = 8;
 53:     lex[3]  = 9;
 54:     lex[9]  = 10;
 55:     lex[5]  = 11;
 56:     lex[11] = 12;
 57:     lex[7]  = 13;
 58:     lex[17] = 14;
 59:     lex[15] = 15;
 60:     lex[19] = 16;
 61:     lex[21] = 17;
 62:     lex[23] = 18;
 63:     lex[25] = 19;
 64:     /* Faces */
 65:     lex[4]  = -20;
 66:     lex[10] = -21;
 67:     lex[12] = -22;
 68:     lex[14] = -23;
 69:     lex[16] = -24;
 70:     lex[22] = -25;
 71:     /* Cell */
 72:     lex[13] = -26;
 73:   }
 74:   return lex;
 75: }

 77: #define GMSH_LEXORDER_LIST(T) \
 78:   GMSH_LEXORDER_ITEM(T, 1) \
 79:   GMSH_LEXORDER_ITEM(T, 2) \
 80:   GMSH_LEXORDER_ITEM(T, 3) \
 81:   GMSH_LEXORDER_ITEM(T, 4) \
 82:   GMSH_LEXORDER_ITEM(T, 5) \
 83:   GMSH_LEXORDER_ITEM(T, 6) \
 84:   GMSH_LEXORDER_ITEM(T, 7) \
 85:   GMSH_LEXORDER_ITEM(T, 8) \
 86:   GMSH_LEXORDER_ITEM(T, 9) \
 87:   GMSH_LEXORDER_ITEM(T, 10)

 89: GMSH_LEXORDER_ITEM(VTX, 0)
 90: GMSH_LEXORDER_LIST(SEG)
 91: GMSH_LEXORDER_LIST(TRI)
 92: GMSH_LEXORDER_LIST(QUA)
 93: GMSH_LEXORDER_LIST(TET)
 94: GMSH_LEXORDER_LIST(HEX)
 95: GMSH_LEXORDER_LIST(PRI)
 96: GMSH_LEXORDER_LIST(PYR)

 98: typedef enum {
 99:   GMSH_VTX           = 0,
100:   GMSH_SEG           = 1,
101:   GMSH_TRI           = 2,
102:   GMSH_QUA           = 3,
103:   GMSH_TET           = 4,
104:   GMSH_HEX           = 5,
105:   GMSH_PRI           = 6,
106:   GMSH_PYR           = 7,
107:   GMSH_NUM_POLYTOPES = 8
108: } GmshPolytopeType;

110: typedef struct {
111:   int cellType;
112:   int polytope;
113:   int dim;
114:   int order;
115:   int numVerts;
116:   int numNodes;
117:   int *(*lexorder)(void);
118: } GmshCellInfo;

120: #define GmshCellEntry(cellType, polytope, dim, order) \
121:   { \
122:     cellType, GMSH_##polytope, dim, order, GmshNumNodes_##polytope(1), GmshNumNodes_##polytope(order), GmshLexOrder_##polytope##_##order \
123:   }

125: static const GmshCellInfo GmshCellTable[] = {
126:   GmshCellEntry(15, VTX, 0, 0),

128:   GmshCellEntry(1, SEG, 1, 1),   GmshCellEntry(8, SEG, 1, 2),   GmshCellEntry(26, SEG, 1, 3),
129:   GmshCellEntry(27, SEG, 1, 4),  GmshCellEntry(28, SEG, 1, 5),  GmshCellEntry(62, SEG, 1, 6),
130:   GmshCellEntry(63, SEG, 1, 7),  GmshCellEntry(64, SEG, 1, 8),  GmshCellEntry(65, SEG, 1, 9),
131:   GmshCellEntry(66, SEG, 1, 10),

133:   GmshCellEntry(2, TRI, 2, 1),   GmshCellEntry(9, TRI, 2, 2),   GmshCellEntry(21, TRI, 2, 3),
134:   GmshCellEntry(23, TRI, 2, 4),  GmshCellEntry(25, TRI, 2, 5),  GmshCellEntry(42, TRI, 2, 6),
135:   GmshCellEntry(43, TRI, 2, 7),  GmshCellEntry(44, TRI, 2, 8),  GmshCellEntry(45, TRI, 2, 9),
136:   GmshCellEntry(46, TRI, 2, 10),

138:   GmshCellEntry(3, QUA, 2, 1),   GmshCellEntry(10, QUA, 2, 2),  {16, GMSH_QUA, 2, 2, 4, 8,  GmshLexOrder_QUA_2_Serendipity},
139:   GmshCellEntry(36, QUA, 2, 3),  GmshCellEntry(37, QUA, 2, 4),  GmshCellEntry(38, QUA, 2, 5),
140:   GmshCellEntry(47, QUA, 2, 6),  GmshCellEntry(48, QUA, 2, 7),  GmshCellEntry(49, QUA, 2, 8),
141:   GmshCellEntry(50, QUA, 2, 9),  GmshCellEntry(51, QUA, 2, 10),

143:   GmshCellEntry(4, TET, 3, 1),   GmshCellEntry(11, TET, 3, 2),  GmshCellEntry(29, TET, 3, 3),
144:   GmshCellEntry(30, TET, 3, 4),  GmshCellEntry(31, TET, 3, 5),  GmshCellEntry(71, TET, 3, 6),
145:   GmshCellEntry(72, TET, 3, 7),  GmshCellEntry(73, TET, 3, 8),  GmshCellEntry(74, TET, 3, 9),
146:   GmshCellEntry(75, TET, 3, 10),

148:   GmshCellEntry(5, HEX, 3, 1),   GmshCellEntry(12, HEX, 3, 2),  {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity},
149:   GmshCellEntry(92, HEX, 3, 3),  GmshCellEntry(93, HEX, 3, 4),  GmshCellEntry(94, HEX, 3, 5),
150:   GmshCellEntry(95, HEX, 3, 6),  GmshCellEntry(96, HEX, 3, 7),  GmshCellEntry(97, HEX, 3, 8),
151:   GmshCellEntry(98, HEX, 3, 9),  GmshCellEntry(-1, HEX, 3, 10),

153:   GmshCellEntry(6, PRI, 3, 1),   GmshCellEntry(13, PRI, 3, 2),  GmshCellEntry(90, PRI, 3, 3),
154:   GmshCellEntry(91, PRI, 3, 4),  GmshCellEntry(106, PRI, 3, 5), GmshCellEntry(107, PRI, 3, 6),
155:   GmshCellEntry(108, PRI, 3, 7), GmshCellEntry(109, PRI, 3, 8), GmshCellEntry(110, PRI, 3, 9),
156:   GmshCellEntry(-1, PRI, 3, 10),

158:   GmshCellEntry(7, PYR, 3, 1),   GmshCellEntry(14, PYR, 3, 2),  GmshCellEntry(118, PYR, 3, 3),
159:   GmshCellEntry(119, PYR, 3, 4), GmshCellEntry(120, PYR, 3, 5), GmshCellEntry(121, PYR, 3, 6),
160:   GmshCellEntry(122, PYR, 3, 7), GmshCellEntry(123, PYR, 3, 8), GmshCellEntry(124, PYR, 3, 9),
161:   GmshCellEntry(-1, PYR, 3, 10)

163: #if 0
164:   {20, GMSH_TRI, 2, 3, 3,  9, NULL},
165:   {18, GMSH_PRI, 3, 2, 6, 15, NULL},
166:   {19, GMSH_PYR, 3, 2, 5, 13, NULL},
167: #endif
168: };

170: static GmshCellInfo GmshCellMap[150];

172: static PetscErrorCode GmshCellInfoSetUp(void)
173: {
174:   size_t           i, n;
175:   static PetscBool called = PETSC_FALSE;

177:   if (called) return 0;
178:   called = PETSC_TRUE;
179:   n      = PETSC_STATIC_ARRAY_LENGTH(GmshCellMap);
180:   for (i = 0; i < n; ++i) {
181:     GmshCellMap[i].cellType = -1;
182:     GmshCellMap[i].polytope = -1;
183:   }
184:   n = PETSC_STATIC_ARRAY_LENGTH(GmshCellTable);
185:   for (i = 0; i < n; ++i) {
186:     if (GmshCellTable[i].cellType <= 0) continue;
187:     GmshCellMap[GmshCellTable[i].cellType] = GmshCellTable[i];
188:   }
189:   return 0;
190: }

192: #define GmshCellTypeCheck(ct) \

196: typedef struct {
197:   PetscViewer viewer;
198:   int         fileFormat;
199:   int         dataSize;
200:   PetscBool   binary;
201:   PetscBool   byteSwap;
202:   size_t      wlen;
203:   void       *wbuf;
204:   size_t      slen;
205:   void       *sbuf;
206:   PetscInt   *nbuf;
207:   PetscInt    nodeStart;
208:   PetscInt    nodeEnd;
209:   PetscInt   *nodeMap;
210: } GmshFile;

212: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, size_t count, size_t eltsize, void *buf)
213: {
214:   size_t size = count * eltsize;

216:   if (gmsh->wlen < size) {
217:     PetscFree(gmsh->wbuf);
218:     PetscMalloc(size, &gmsh->wbuf);
219:     gmsh->wlen = size;
220:   }
221:   *(void **)buf = size ? gmsh->wbuf : NULL;
222:   return 0;
223: }

225: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, size_t count, void *buf)
226: {
227:   size_t dataSize = (size_t)gmsh->dataSize;
228:   size_t size     = count * dataSize;

230:   if (gmsh->slen < size) {
231:     PetscFree(gmsh->sbuf);
232:     PetscMalloc(size, &gmsh->sbuf);
233:     gmsh->slen = size;
234:   }
235:   *(void **)buf = size ? gmsh->sbuf : NULL;
236:   return 0;
237: }

239: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
240: {
241:   PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype);
242:   if (gmsh->byteSwap) PetscByteSwap(buf, dtype, count);
243:   return 0;
244: }

246: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
247: {
248:   PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING);
249:   return 0;
250: }

252: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
253: {
254:   PetscStrcmp(line, Section, match);
255:   return 0;
256: }

258: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
259: {
260:   PetscBool match;

262:   GmshMatch(gmsh, Section, line, &match);
264:   return 0;
265: }

267: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
268: {
269:   PetscBool match;

271:   while (PETSC_TRUE) {
272:     GmshReadString(gmsh, line, 1);
273:     GmshMatch(gmsh, "$Comments", line, &match);
274:     if (!match) break;
275:     while (PETSC_TRUE) {
276:       GmshReadString(gmsh, line, 1);
277:       GmshMatch(gmsh, "$EndComments", line, &match);
278:       if (match) break;
279:     }
280:   }
281:   return 0;
282: }

284: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
285: {
286:   GmshReadString(gmsh, line, 1);
287:   GmshExpect(gmsh, EndSection, line);
288:   return 0;
289: }

291: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
292: {
293:   PetscInt i;
294:   size_t   dataSize = (size_t)gmsh->dataSize;

296:   if (dataSize == sizeof(PetscInt)) {
297:     GmshRead(gmsh, buf, count, PETSC_INT);
298:   } else if (dataSize == sizeof(int)) {
299:     int *ibuf = NULL;
300:     GmshBufferSizeGet(gmsh, count, &ibuf);
301:     GmshRead(gmsh, ibuf, count, PETSC_ENUM);
302:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
303:   } else if (dataSize == sizeof(long)) {
304:     long *ibuf = NULL;
305:     GmshBufferSizeGet(gmsh, count, &ibuf);
306:     GmshRead(gmsh, ibuf, count, PETSC_LONG);
307:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
308:   } else if (dataSize == sizeof(PetscInt64)) {
309:     PetscInt64 *ibuf = NULL;
310:     GmshBufferSizeGet(gmsh, count, &ibuf);
311:     GmshRead(gmsh, ibuf, count, PETSC_INT64);
312:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
313:   }
314:   return 0;
315: }

317: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
318: {
319:   GmshRead(gmsh, buf, count, PETSC_ENUM);
320:   return 0;
321: }

323: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
324: {
325:   GmshRead(gmsh, buf, count, PETSC_DOUBLE);
326:   return 0;
327: }

329: #define GMSH_MAX_TAGS 16

331: typedef struct {
332:   PetscInt id;                  /* Entity ID */
333:   PetscInt dim;                 /* Dimension */
334:   double   bbox[6];             /* Bounding box */
335:   PetscInt numTags;             /* Size of tag array */
336:   int      tags[GMSH_MAX_TAGS]; /* Tag array */
337: } GmshEntity;

339: typedef struct {
340:   GmshEntity *entity[4];
341:   PetscHMapI  entityMap[4];
342: } GmshEntities;

344: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
345: {
346:   PetscInt dim;

348:   PetscNew(entities);
349:   for (dim = 0; dim < 4; ++dim) {
350:     PetscCalloc1(count[dim], &(*entities)->entity[dim]);
351:     PetscHMapICreate(&(*entities)->entityMap[dim]);
352:   }
353:   return 0;
354: }

356: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
357: {
358:   PetscInt dim;

360:   if (!*entities) return 0;
361:   for (dim = 0; dim < 4; ++dim) {
362:     PetscFree((*entities)->entity[dim]);
363:     PetscHMapIDestroy(&(*entities)->entityMap[dim]);
364:   }
365:   PetscFree((*entities));
366:   return 0;
367: }

369: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity **entity)
370: {
371:   PetscHMapISet(entities->entityMap[dim], eid, index);
372:   entities->entity[dim][index].dim = dim;
373:   entities->entity[dim][index].id  = eid;
374:   if (entity) *entity = &entities->entity[dim][index];
375:   return 0;
376: }

378: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity **entity)
379: {
380:   PetscInt index;

382:   PetscHMapIGet(entities->entityMap[dim], eid, &index);
383:   *entity = &entities->entity[dim][index];
384:   return 0;
385: }

387: typedef struct {
388:   PetscInt *id;  /* Node IDs */
389:   double   *xyz; /* Coordinates */
390:   PetscInt *tag; /* Physical tag */
391: } GmshNodes;

393: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
394: {
395:   PetscNew(nodes);
396:   PetscMalloc1(count * 1, &(*nodes)->id);
397:   PetscMalloc1(count * 3, &(*nodes)->xyz);
398:   PetscMalloc1(count * GMSH_MAX_TAGS, &(*nodes)->tag);
399:   return 0;
400: }

402: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
403: {
404:   if (!*nodes) return 0;
405:   PetscFree((*nodes)->id);
406:   PetscFree((*nodes)->xyz);
407:   PetscFree((*nodes)->tag);
408:   PetscFree((*nodes));
409:   return 0;
410: }

412: typedef struct {
413:   PetscInt  id;                  /* Element ID */
414:   PetscInt  dim;                 /* Dimension */
415:   PetscInt  cellType;            /* Cell type */
416:   PetscInt  numVerts;            /* Size of vertex array */
417:   PetscInt  numNodes;            /* Size of node array */
418:   PetscInt *nodes;               /* Vertex/Node array */
419:   PetscInt  numTags;             /* Size of physical tag array */
420:   int       tags[GMSH_MAX_TAGS]; /* Physical tag array */
421: } GmshElement;

423: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
424: {
425:   PetscCalloc1(count, elements);
426:   return 0;
427: }

429: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
430: {
431:   if (!*elements) return 0;
432:   PetscFree(*elements);
433:   return 0;
434: }

436: typedef struct {
437:   PetscInt       dim;
438:   PetscInt       order;
439:   GmshEntities  *entities;
440:   PetscInt       numNodes;
441:   GmshNodes     *nodelist;
442:   PetscInt       numElems;
443:   GmshElement   *elements;
444:   PetscInt       numVerts;
445:   PetscInt       numCells;
446:   PetscInt      *periodMap;
447:   PetscInt      *vertexMap;
448:   PetscSegBuffer segbuf;
449:   PetscInt       numRegions;
450:   PetscInt      *regionDims;
451:   PetscInt      *regionTags;
452:   char         **regionNames;
453: } GmshMesh;

455: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
456: {
457:   PetscNew(mesh);
458:   PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf);
459:   return 0;
460: }

462: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
463: {
464:   PetscInt r;

466:   if (!*mesh) return 0;
467:   GmshEntitiesDestroy(&(*mesh)->entities);
468:   GmshNodesDestroy(&(*mesh)->nodelist);
469:   GmshElementsDestroy(&(*mesh)->elements);
470:   PetscFree((*mesh)->periodMap);
471:   PetscFree((*mesh)->vertexMap);
472:   PetscSegBufferDestroy(&(*mesh)->segbuf);
473:   for (r = 0; r < (*mesh)->numRegions; ++r) PetscFree((*mesh)->regionNames[r]);
474:   PetscFree3((*mesh)->regionDims, (*mesh)->regionTags, (*mesh)->regionNames);
475:   PetscFree((*mesh));
476:   return 0;
477: }

479: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
480: {
481:   PetscViewer viewer   = gmsh->viewer;
482:   PetscBool   byteSwap = gmsh->byteSwap;
483:   char        line[PETSC_MAX_PATH_LEN];
484:   int         n, t, num, nid, snum;
485:   GmshNodes  *nodes;

487:   PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
488:   snum = sscanf(line, "%d", &num);
490:   GmshNodesCreate(num, &nodes);
491:   mesh->numNodes = num;
492:   mesh->nodelist = nodes;
493:   for (n = 0; n < num; ++n) {
494:     double *xyz = nodes->xyz + n * 3;
495:     PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
496:     PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
497:     if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
498:     if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
499:     nodes->id[n] = nid;
500:     for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
501:   }
502:   return 0;
503: }

505: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
506:    file contents multiple times to figure out the true number of cells and facets
507:    in the given mesh. To make this more efficient we read the file contents only
508:    once and store them in memory, while determining the true number of cells. */
509: static PetscErrorCode GmshReadElements_v22(GmshFile *gmsh, GmshMesh *mesh)
510: {
511:   PetscViewer  viewer   = gmsh->viewer;
512:   PetscBool    binary   = gmsh->binary;
513:   PetscBool    byteSwap = gmsh->byteSwap;
514:   char         line[PETSC_MAX_PATH_LEN];
515:   int          i, c, p, num, ibuf[1 + 4 + 1000], snum;
516:   int          cellType, numElem, numVerts, numNodes, numTags;
517:   GmshElement *elements;
518:   PetscInt    *nodeMap = gmsh->nodeMap;

520:   PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
521:   snum = sscanf(line, "%d", &num);
523:   GmshElementsCreate(num, &elements);
524:   mesh->numElems = num;
525:   mesh->elements = elements;
526:   for (c = 0; c < num;) {
527:     PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
528:     if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);

530:     cellType = binary ? ibuf[0] : ibuf[1];
531:     numElem  = binary ? ibuf[1] : 1;
532:     numTags  = ibuf[2];

534:     GmshCellTypeCheck(cellType);
535:     numVerts = GmshCellMap[cellType].numVerts;
536:     numNodes = GmshCellMap[cellType].numNodes;

538:     for (i = 0; i < numElem; ++i, ++c) {
539:       GmshElement *element = elements + c;
540:       const int    off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
541:       const int   *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
542:       PetscViewerRead(viewer, ibuf + off, nint, NULL, PETSC_ENUM);
543:       if (byteSwap) PetscByteSwap(ibuf + off, PETSC_ENUM, nint);
544:       element->id       = id[0];
545:       element->dim      = GmshCellMap[cellType].dim;
546:       element->cellType = cellType;
547:       element->numVerts = numVerts;
548:       element->numNodes = numNodes;
549:       element->numTags  = PetscMin(numTags, GMSH_MAX_TAGS);
550:       PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
551:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
552:       for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
553:     }
554:   }
555:   return 0;
556: }

558: /*
559: $Entities
560:   numPoints(unsigned long) numCurves(unsigned long)
561:     numSurfaces(unsigned long) numVolumes(unsigned long)
562:   // points
563:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
564:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
565:     numPhysicals(unsigned long) phyisicalTag[...](int)
566:   ...
567:   // curves
568:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
569:      boxMaxX(double) boxMaxY(double) boxMaxZ(double)
570:      numPhysicals(unsigned long) physicalTag[...](int)
571:      numBREPVert(unsigned long) tagBREPVert[...](int)
572:   ...
573:   // surfaces
574:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
575:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
576:     numPhysicals(unsigned long) physicalTag[...](int)
577:     numBREPCurve(unsigned long) tagBREPCurve[...](int)
578:   ...
579:   // volumes
580:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
581:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
582:     numPhysicals(unsigned long) physicalTag[...](int)
583:     numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
584:   ...
585: $EndEntities
586: */
587: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
588: {
589:   PetscViewer viewer   = gmsh->viewer;
590:   PetscBool   byteSwap = gmsh->byteSwap;
591:   long        index, num, lbuf[4];
592:   int         dim, eid, numTags, *ibuf, t;
593:   PetscInt    count[4], i;
594:   GmshEntity *entity = NULL;

596:   PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG);
597:   if (byteSwap) PetscByteSwap(lbuf, PETSC_LONG, 4);
598:   for (i = 0; i < 4; ++i) count[i] = lbuf[i];
599:   GmshEntitiesCreate(count, &mesh->entities);
600:   for (dim = 0; dim < 4; ++dim) {
601:     for (index = 0; index < count[dim]; ++index) {
602:       PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM);
603:       if (byteSwap) PetscByteSwap(&eid, PETSC_ENUM, 1);
604:       GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
605:       PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE);
606:       if (byteSwap) PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6);
607:       PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
608:       if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
609:       GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
610:       PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
611:       if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
612:       entity->numTags = numTags = (int)PetscMin(num, GMSH_MAX_TAGS);
613:       for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
614:       if (dim == 0) continue;
615:       PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG);
616:       if (byteSwap) PetscByteSwap(&num, PETSC_LONG, 1);
617:       GmshBufferGet(gmsh, num, sizeof(int), &ibuf);
618:       PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM);
619:       if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, num);
620:     }
621:   }
622:   return 0;
623: }

625: /*
626: $Nodes
627:   numEntityBlocks(unsigned long) numNodes(unsigned long)
628:   tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
629:     tag(int) x(double) y(double) z(double)
630:     ...
631:   ...
632: $EndNodes
633: */
634: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
635: {
636:   PetscViewer viewer   = gmsh->viewer;
637:   PetscBool   byteSwap = gmsh->byteSwap;
638:   long        block, node, n, t, numEntityBlocks, numTotalNodes, numNodes;
639:   int         info[3], nid;
640:   GmshNodes  *nodes;

642:   PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
643:   if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
644:   PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG);
645:   if (byteSwap) PetscByteSwap(&numTotalNodes, PETSC_LONG, 1);
646:   GmshNodesCreate(numTotalNodes, &nodes);
647:   mesh->numNodes = numTotalNodes;
648:   mesh->nodelist = nodes;
649:   for (n = 0, block = 0; block < numEntityBlocks; ++block) {
650:     PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
651:     PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG);
652:     if (byteSwap) PetscByteSwap(&numNodes, PETSC_LONG, 1);
653:     if (gmsh->binary) {
654:       size_t nbytes = sizeof(int) + 3 * sizeof(double);
655:       char  *cbuf   = NULL; /* dummy value to prevent warning from compiler about possible uninitialized value */
656:       GmshBufferGet(gmsh, numNodes, nbytes, &cbuf);
657:       PetscViewerRead(viewer, cbuf, numNodes * nbytes, NULL, PETSC_CHAR);
658:       for (node = 0; node < numNodes; ++node, ++n) {
659:         char   *cnid = cbuf + node * nbytes, *cxyz = cnid + sizeof(int);
660:         double *xyz = nodes->xyz + n * 3;
661:         if (!PetscBinaryBigEndian()) PetscByteSwap(cnid, PETSC_ENUM, 1);
662:         if (!PetscBinaryBigEndian()) PetscByteSwap(cxyz, PETSC_DOUBLE, 3);
663:         PetscMemcpy(&nid, cnid, sizeof(int));
664:         PetscMemcpy(xyz, cxyz, 3 * sizeof(double));
665:         if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
666:         if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
667:         nodes->id[n] = nid;
668:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
669:       }
670:     } else {
671:       for (node = 0; node < numNodes; ++node, ++n) {
672:         double *xyz = nodes->xyz + n * 3;
673:         PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM);
674:         PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE);
675:         if (byteSwap) PetscByteSwap(&nid, PETSC_ENUM, 1);
676:         if (byteSwap) PetscByteSwap(xyz, PETSC_DOUBLE, 3);
677:         nodes->id[n] = nid;
678:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
679:       }
680:     }
681:   }
682:   return 0;
683: }

685: /*
686: $Elements
687:   numEntityBlocks(unsigned long) numElements(unsigned long)
688:   tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
689:     tag(int) numVert[...](int)
690:     ...
691:   ...
692: $EndElements
693: */
694: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
695: {
696:   PetscViewer  viewer   = gmsh->viewer;
697:   PetscBool    byteSwap = gmsh->byteSwap;
698:   long         c, block, numEntityBlocks, numTotalElements, elem, numElements;
699:   int          p, info[3], *ibuf = NULL;
700:   int          eid, dim, cellType, numVerts, numNodes, numTags;
701:   GmshEntity  *entity = NULL;
702:   GmshElement *elements;
703:   PetscInt    *nodeMap = gmsh->nodeMap;

705:   PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG);
706:   if (byteSwap) PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1);
707:   PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG);
708:   if (byteSwap) PetscByteSwap(&numTotalElements, PETSC_LONG, 1);
709:   GmshElementsCreate(numTotalElements, &elements);
710:   mesh->numElems = numTotalElements;
711:   mesh->elements = elements;
712:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
713:     PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM);
714:     if (byteSwap) PetscByteSwap(info, PETSC_ENUM, 3);
715:     eid      = info[0];
716:     dim      = info[1];
717:     cellType = info[2];
718:     GmshEntitiesGet(mesh->entities, dim, eid, &entity);
719:     GmshCellTypeCheck(cellType);
720:     numVerts = GmshCellMap[cellType].numVerts;
721:     numNodes = GmshCellMap[cellType].numNodes;
722:     numTags  = entity->numTags;
723:     PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG);
724:     if (byteSwap) PetscByteSwap(&numElements, PETSC_LONG, 1);
725:     GmshBufferGet(gmsh, (1 + numNodes) * numElements, sizeof(int), &ibuf);
726:     PetscViewerRead(viewer, ibuf, (1 + numNodes) * numElements, NULL, PETSC_ENUM);
727:     if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, (1 + numNodes) * numElements);
728:     for (elem = 0; elem < numElements; ++elem, ++c) {
729:       GmshElement *element = elements + c;
730:       const int   *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
731:       element->id       = id[0];
732:       element->dim      = dim;
733:       element->cellType = cellType;
734:       element->numVerts = numVerts;
735:       element->numNodes = numNodes;
736:       element->numTags  = numTags;
737:       PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
738:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
739:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
740:     }
741:   }
742:   return 0;
743: }

745: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
746: {
747:   PetscViewer viewer     = gmsh->viewer;
748:   int         fileFormat = gmsh->fileFormat;
749:   PetscBool   binary     = gmsh->binary;
750:   PetscBool   byteSwap   = gmsh->byteSwap;
751:   int         numPeriodic, snum, i;
752:   char        line[PETSC_MAX_PATH_LEN];
753:   PetscInt   *nodeMap = gmsh->nodeMap;

755:   if (fileFormat == 22 || !binary) {
756:     PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
757:     snum = sscanf(line, "%d", &numPeriodic);
759:   } else {
760:     PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM);
761:     if (byteSwap) PetscByteSwap(&numPeriodic, PETSC_ENUM, 1);
762:   }
763:   for (i = 0; i < numPeriodic; i++) {
764:     int    ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
765:     long   j, nNodes;
766:     double affine[16];

768:     if (fileFormat == 22 || !binary) {
769:       PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING);
770:       snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
772:     } else {
773:       PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM);
774:       if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 3);
775:       correspondingDim = ibuf[0];
776:       correspondingTag = ibuf[1];
777:       primaryTag       = ibuf[2];
778:     }
779:     (void)correspondingDim;
780:     (void)correspondingTag;
781:     (void)primaryTag; /* unused */

783:     if (fileFormat == 22 || !binary) {
784:       PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
785:       snum = sscanf(line, "%ld", &nNodes);
786:       if (snum != 1) { /* discard transformation and try again */
787:         PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING);
788:         PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING);
789:         snum = sscanf(line, "%ld", &nNodes);
791:       }
792:     } else {
793:       PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
794:       if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
795:       if (nNodes == -1) { /* discard transformation and try again */
796:         PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE);
797:         PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG);
798:         if (byteSwap) PetscByteSwap(&nNodes, PETSC_LONG, 1);
799:       }
800:     }

802:     for (j = 0; j < nNodes; j++) {
803:       if (fileFormat == 22 || !binary) {
804:         PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);
805:         snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
807:       } else {
808:         PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM);
809:         if (byteSwap) PetscByteSwap(ibuf, PETSC_ENUM, 2);
810:         correspondingNode = ibuf[0];
811:         primaryNode       = ibuf[1];
812:       }
813:       correspondingNode              = (int)nodeMap[correspondingNode];
814:       primaryNode                    = (int)nodeMap[primaryNode];
815:       periodicMap[correspondingNode] = primaryNode;
816:     }
817:   }
818:   return 0;
819: }

821: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
822: $Entities
823:   numPoints(size_t) numCurves(size_t)
824:     numSurfaces(size_t) numVolumes(size_t)
825:   pointTag(int) X(double) Y(double) Z(double)
826:     numPhysicalTags(size_t) physicalTag(int) ...
827:   ...
828:   curveTag(int) minX(double) minY(double) minZ(double)
829:     maxX(double) maxY(double) maxZ(double)
830:     numPhysicalTags(size_t) physicalTag(int) ...
831:     numBoundingPoints(size_t) pointTag(int) ...
832:   ...
833:   surfaceTag(int) minX(double) minY(double) minZ(double)
834:     maxX(double) maxY(double) maxZ(double)
835:     numPhysicalTags(size_t) physicalTag(int) ...
836:     numBoundingCurves(size_t) curveTag(int) ...
837:   ...
838:   volumeTag(int) minX(double) minY(double) minZ(double)
839:     maxX(double) maxY(double) maxZ(double)
840:     numPhysicalTags(size_t) physicalTag(int) ...
841:     numBoundngSurfaces(size_t) surfaceTag(int) ...
842:   ...
843: $EndEntities
844: */
845: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
846: {
847:   PetscInt    count[4], index, numTags, i;
848:   int         dim, eid, *tags = NULL;
849:   GmshEntity *entity = NULL;

851:   GmshReadSize(gmsh, count, 4);
852:   GmshEntitiesCreate(count, &mesh->entities);
853:   for (dim = 0; dim < 4; ++dim) {
854:     for (index = 0; index < count[dim]; ++index) {
855:       GmshReadInt(gmsh, &eid, 1);
856:       GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity);
857:       GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6);
858:       GmshReadSize(gmsh, &numTags, 1);
859:       GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
860:       GmshReadInt(gmsh, tags, numTags);
862:       entity->numTags = numTags;
863:       for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
864:       if (dim == 0) continue;
865:       GmshReadSize(gmsh, &numTags, 1);
866:       GmshBufferGet(gmsh, numTags, sizeof(int), &tags);
867:       GmshReadInt(gmsh, tags, numTags);
868:       /* Currently, we do not save the ids for the bounding entities */
869:     }
870:   }
871:   return 0;
872: }

874: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
875: $Nodes
876:   numEntityBlocks(size_t) numNodes(size_t)
877:     minNodeTag(size_t) maxNodeTag(size_t)
878:   entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
879:     nodeTag(size_t)
880:     ...
881:     x(double) y(double) z(double)
882:        < u(double; if parametric and entityDim = 1 or entityDim = 2) >
883:        < v(double; if parametric and entityDim = 2) >
884:     ...
885:   ...
886: $EndNodes
887: */
888: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
889: {
890:   int         info[3], dim, eid, parametric;
891:   PetscInt    sizes[4], numEntityBlocks, numTags, t, numNodes, numNodesBlock = 0, block, node, n;
892:   GmshEntity *entity = NULL;
893:   GmshNodes  *nodes;

895:   GmshReadSize(gmsh, sizes, 4);
896:   numEntityBlocks = sizes[0];
897:   numNodes        = sizes[1];
898:   GmshNodesCreate(numNodes, &nodes);
899:   mesh->numNodes = numNodes;
900:   mesh->nodelist = nodes;
901:   for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
902:     GmshReadInt(gmsh, info, 3);
903:     dim        = info[0];
904:     eid        = info[1];
905:     parametric = info[2];
906:     GmshEntitiesGet(mesh->entities, dim, eid, &entity);
907:     numTags = entity->numTags;
909:     GmshReadSize(gmsh, &numNodesBlock, 1);
910:     GmshReadSize(gmsh, nodes->id + node, numNodesBlock);
911:     GmshReadDouble(gmsh, nodes->xyz + node * 3, numNodesBlock * 3);
912:     for (n = 0; n < numNodesBlock; ++n) {
913:       PetscInt *tags = &nodes->tag[node * GMSH_MAX_TAGS];

915:       for (t = 0; t < numTags; ++t) tags[n * GMSH_MAX_TAGS + t] = entity->tags[t];
916:       for (t = numTags; t < GMSH_MAX_TAGS; ++t) tags[n * GMSH_MAX_TAGS + t] = -1;
917:     }
918:   }
919:   gmsh->nodeStart = sizes[2];
920:   gmsh->nodeEnd   = sizes[3] + 1;
921:   return 0;
922: }

924: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
925: $Elements
926:   numEntityBlocks(size_t) numElements(size_t)
927:     minElementTag(size_t) maxElementTag(size_t)
928:   entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
929:     elementTag(size_t) nodeTag(size_t) ...
930:     ...
931:   ...
932: $EndElements
933: */
934: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
935: {
936:   int          info[3], eid, dim, cellType;
937:   PetscInt     sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
938:   GmshEntity  *entity = NULL;
939:   GmshElement *elements;
940:   PetscInt    *nodeMap = gmsh->nodeMap;

942:   GmshReadSize(gmsh, sizes, 4);
943:   numEntityBlocks = sizes[0];
944:   numElements     = sizes[1];
945:   GmshElementsCreate(numElements, &elements);
946:   mesh->numElems = numElements;
947:   mesh->elements = elements;
948:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
949:     GmshReadInt(gmsh, info, 3);
950:     dim      = info[0];
951:     eid      = info[1];
952:     cellType = info[2];
953:     GmshEntitiesGet(mesh->entities, dim, eid, &entity);
954:     GmshCellTypeCheck(cellType);
955:     numVerts = GmshCellMap[cellType].numVerts;
956:     numNodes = GmshCellMap[cellType].numNodes;
957:     numTags  = entity->numTags;
958:     GmshReadSize(gmsh, &numBlockElements, 1);
959:     GmshBufferGet(gmsh, (1 + numNodes) * numBlockElements, sizeof(PetscInt), &ibuf);
960:     GmshReadSize(gmsh, ibuf, (1 + numNodes) * numBlockElements);
961:     for (elem = 0; elem < numBlockElements; ++elem, ++c) {
962:       GmshElement    *element = elements + c;
963:       const PetscInt *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
964:       element->id       = id[0];
965:       element->dim      = dim;
966:       element->cellType = cellType;
967:       element->numVerts = numVerts;
968:       element->numNodes = numNodes;
969:       element->numTags  = numTags;
970:       PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes);
971:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
972:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
973:     }
974:   }
975:   return 0;
976: }

978: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
979: $Periodic
980:   numPeriodicLinks(size_t)
981:   entityDim(int) entityTag(int) entityTagPrimary(int)
982:   numAffine(size_t) value(double) ...
983:   numCorrespondingNodes(size_t)
984:     nodeTag(size_t) nodeTagPrimary(size_t)
985:     ...
986:   ...
987: $EndPeriodic
988: */
989: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
990: {
991:   int       info[3];
992:   double    dbuf[16];
993:   PetscInt  numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
994:   PetscInt *nodeMap = gmsh->nodeMap;

996:   GmshReadSize(gmsh, &numPeriodicLinks, 1);
997:   for (link = 0; link < numPeriodicLinks; ++link) {
998:     GmshReadInt(gmsh, info, 3);
999:     GmshReadSize(gmsh, &numAffine, 1);
1000:     GmshReadDouble(gmsh, dbuf, numAffine);
1001:     GmshReadSize(gmsh, &numCorrespondingNodes, 1);
1002:     GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags);
1003:     GmshReadSize(gmsh, nodeTags, numCorrespondingNodes * 2);
1004:     for (node = 0; node < numCorrespondingNodes; ++node) {
1005:       PetscInt correspondingNode     = nodeMap[nodeTags[node * 2 + 0]];
1006:       PetscInt primaryNode           = nodeMap[nodeTags[node * 2 + 1]];
1007:       periodicMap[correspondingNode] = primaryNode;
1008:     }
1009:   }
1010:   return 0;
1011: }

1013: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1014: $MeshFormat // same as MSH version 2
1015:   version(ASCII double; currently 4.1)
1016:   file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1017:   data-size(ASCII int; sizeof(size_t))
1018:   < int with value one; only in binary mode, to detect endianness >
1019: $EndMeshFormat
1020: */
1021: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1022: {
1023:   char  line[PETSC_MAX_PATH_LEN];
1024:   int   snum, fileType, fileFormat, dataSize, checkEndian;
1025:   float version;

1027:   GmshReadString(gmsh, line, 3);
1028:   snum       = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1029:   fileFormat = (int)roundf(version * 10);
1038:   gmsh->fileFormat = fileFormat;
1039:   gmsh->dataSize   = dataSize;
1040:   gmsh->byteSwap   = PETSC_FALSE;
1041:   if (gmsh->binary) {
1042:     GmshReadInt(gmsh, &checkEndian, 1);
1043:     if (checkEndian != 1) {
1044:       PetscByteSwap(&checkEndian, PETSC_ENUM, 1);
1046:       gmsh->byteSwap = PETSC_TRUE;
1047:     }
1048:   }
1049:   return 0;
1050: }

1052: /*
1053: PhysicalNames
1054:   numPhysicalNames(ASCII int)
1055:   dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1056:   ...
1057: $EndPhysicalNames
1058: */
1059: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1060: {
1061:   char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p, *q, *r;
1062:   int  snum, region, dim, tag;

1064:   GmshReadString(gmsh, line, 1);
1065:   snum             = sscanf(line, "%d", &region);
1066:   mesh->numRegions = region;
1068:   PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames);
1069:   for (region = 0; region < mesh->numRegions; ++region) {
1070:     GmshReadString(gmsh, line, 2);
1071:     snum = sscanf(line, "%d %d", &dim, &tag);
1073:     GmshReadString(gmsh, line, -(PetscInt)sizeof(line));
1074:     PetscStrchr(line, '"', &p);
1076:     PetscStrrchr(line, '"', &q);
1078:     PetscStrrchr(line, ':', &r);
1079:     if (p != r) q = r;
1080:     PetscStrncpy(name, p + 1, (size_t)(q - p - 1));
1081:     mesh->regionDims[region] = dim;
1082:     mesh->regionTags[region] = tag;
1083:     PetscStrallocpy(name, &mesh->regionNames[region]);
1084:   }
1085:   return 0;
1086: }

1088: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1089: {
1090:   switch (gmsh->fileFormat) {
1091:   case 41:
1092:     GmshReadEntities_v41(gmsh, mesh);
1093:     break;
1094:   default:
1095:     GmshReadEntities_v40(gmsh, mesh);
1096:     break;
1097:   }
1098:   return 0;
1099: }

1101: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1102: {
1103:   switch (gmsh->fileFormat) {
1104:   case 41:
1105:     GmshReadNodes_v41(gmsh, mesh);
1106:     break;
1107:   case 40:
1108:     GmshReadNodes_v40(gmsh, mesh);
1109:     break;
1110:   default:
1111:     GmshReadNodes_v22(gmsh, mesh);
1112:     break;
1113:   }

1115:   { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1116:     if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1117:       PetscInt   tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1118:       GmshNodes *nodes = mesh->nodelist;
1119:       for (n = 0; n < mesh->numNodes; ++n) {
1120:         const PetscInt tag = nodes->id[n];
1121:         tagMin             = PetscMin(tag, tagMin);
1122:         tagMax             = PetscMax(tag, tagMax);
1123:       }
1124:       gmsh->nodeStart = tagMin;
1125:       gmsh->nodeEnd   = tagMax + 1;
1126:     }
1127:   }

1129:   { /* Support for sparse node tags */
1130:     PetscInt   n, t;
1131:     GmshNodes *nodes = mesh->nodelist;
1132:     PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf);
1133:     for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1134:     gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1135:     for (n = 0; n < mesh->numNodes; ++n) {
1136:       const PetscInt tag = nodes->id[n];
1138:       gmsh->nodeMap[tag] = n;
1139:     }
1140:   }
1141:   return 0;
1142: }

1144: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1145: {
1146:   switch (gmsh->fileFormat) {
1147:   case 41:
1148:     GmshReadElements_v41(gmsh, mesh);
1149:     break;
1150:   case 40:
1151:     GmshReadElements_v40(gmsh, mesh);
1152:     break;
1153:   default:
1154:     GmshReadElements_v22(gmsh, mesh);
1155:     break;
1156:   }

1158:   { /* Reorder elements by codimension and polytope type */
1159:     PetscInt     ne       = mesh->numElems;
1160:     GmshElement *elements = mesh->elements;
1161:     PetscInt     keymap[GMSH_NUM_POLYTOPES], nk = 0;
1162:     PetscInt     offset[GMSH_NUM_POLYTOPES + 1], e, k;

1164:     for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1165:     PetscMemzero(offset, sizeof(offset));

1167:     keymap[GMSH_TET] = nk++;
1168:     keymap[GMSH_HEX] = nk++;
1169:     keymap[GMSH_PRI] = nk++;
1170:     keymap[GMSH_PYR] = nk++;
1171:     keymap[GMSH_TRI] = nk++;
1172:     keymap[GMSH_QUA] = nk++;
1173:     keymap[GMSH_SEG] = nk++;
1174:     keymap[GMSH_VTX] = nk++;

1176:     GmshElementsCreate(mesh->numElems, &mesh->elements);
1177: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1178:     for (e = 0; e < ne; ++e) offset[1 + key(e)]++;
1179:     for (k = 1; k < nk; ++k) offset[k] += offset[k - 1];
1180:     for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1181: #undef key
1182:     GmshElementsDestroy(&elements);
1183:   }

1185:   { /* Mesh dimension and order */
1186:     GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1187:     mesh->dim         = elem ? GmshCellMap[elem->cellType].dim : 0;
1188:     mesh->order       = elem ? GmshCellMap[elem->cellType].order : 0;
1189:   }

1191:   {
1192:     PetscBT  vtx;
1193:     PetscInt dim = mesh->dim, e, n, v;

1195:     PetscBTCreate(mesh->numNodes, &vtx);

1197:     /* Compute number of cells and set of vertices */
1198:     mesh->numCells = 0;
1199:     for (e = 0; e < mesh->numElems; ++e) {
1200:       GmshElement *elem = mesh->elements + e;
1201:       if (elem->dim == dim && dim > 0) mesh->numCells++;
1202:       for (v = 0; v < elem->numVerts; v++) PetscBTSet(vtx, elem->nodes[v]);
1203:     }

1205:     /* Compute numbering for vertices */
1206:     mesh->numVerts = 0;
1207:     PetscMalloc1(mesh->numNodes, &mesh->vertexMap);
1208:     for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;

1210:     PetscBTDestroy(&vtx);
1211:   }
1212:   return 0;
1213: }

1215: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1216: {
1217:   PetscInt n;

1219:   PetscMalloc1(mesh->numNodes, &mesh->periodMap);
1220:   for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1221:   switch (gmsh->fileFormat) {
1222:   case 41:
1223:     GmshReadPeriodic_v41(gmsh, mesh->periodMap);
1224:     break;
1225:   default:
1226:     GmshReadPeriodic_v40(gmsh, mesh->periodMap);
1227:     break;
1228:   }

1230:   /* Find canonical primary nodes */
1231:   for (n = 0; n < mesh->numNodes; ++n)
1232:     while (mesh->periodMap[n] != mesh->periodMap[mesh->periodMap[n]]) mesh->periodMap[n] = mesh->periodMap[mesh->periodMap[n]];

1234:   /* Renumber vertices (filter out correspondings) */
1235:   mesh->numVerts = 0;
1236:   for (n = 0; n < mesh->numNodes; ++n)
1237:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1238:       if (mesh->periodMap[n] == n) /* is primary */
1239:         mesh->vertexMap[n] = mesh->numVerts++;
1240:   for (n = 0; n < mesh->numNodes; ++n)
1241:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1242:       if (mesh->periodMap[n] != n) /* is corresponding  */
1243:         mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1244:   return 0;
1245: }

1247: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1248: static const DMPolytopeType DMPolytopeMap[] = {
1249:   /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1250:   /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1251:   /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1252:   /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1253:   /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1254:   /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1255:   /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1256:   /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,       DM_POLYTOPE_UNKNOWN};

1258: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1259: {
1260:   return DMPolytopeMap[GmshCellMap[cellType].polytope];
1261: }

1263: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1264: {
1265:   DM              K;
1266:   PetscSpace      P;
1267:   PetscDualSpace  Q;
1268:   PetscQuadrature q, fq;
1269:   PetscBool       isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1270:   PetscBool       endpoint = PETSC_TRUE;
1271:   char            name[32];

1273:   /* Create space */
1274:   PetscSpaceCreate(comm, &P);
1275:   PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL);
1276:   PetscSpacePolynomialSetTensor(P, isTensor);
1277:   PetscSpaceSetNumComponents(P, Nc);
1278:   PetscSpaceSetNumVariables(P, dim);
1279:   PetscSpaceSetDegree(P, k, PETSC_DETERMINE);
1280:   if (prefix) {
1281:     PetscObjectSetOptionsPrefix((PetscObject)P, prefix);
1282:     PetscSpaceSetFromOptions(P);
1283:     PetscObjectSetOptionsPrefix((PetscObject)P, NULL);
1284:     PetscSpaceGetDegree(P, &k, NULL);
1285:   }
1286:   PetscSpaceSetUp(P);
1287:   /* Create dual space */
1288:   PetscDualSpaceCreate(comm, &Q);
1289:   PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE);
1290:   PetscDualSpaceLagrangeSetTensor(Q, isTensor);
1291:   PetscDualSpaceLagrangeSetContinuity(Q, continuity);
1292:   PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0);
1293:   PetscDualSpaceSetNumComponents(Q, Nc);
1294:   PetscDualSpaceSetOrder(Q, k);
1295:   DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K);
1296:   PetscDualSpaceSetDM(Q, K);
1297:   DMDestroy(&K);
1298:   if (prefix) {
1299:     PetscObjectSetOptionsPrefix((PetscObject)Q, prefix);
1300:     PetscDualSpaceSetFromOptions(Q);
1301:     PetscObjectSetOptionsPrefix((PetscObject)Q, NULL);
1302:   }
1303:   PetscDualSpaceSetUp(Q);
1304:   /* Create quadrature */
1305:   if (isSimplex) {
1306:     PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q);
1307:     PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq);
1308:   } else {
1309:     PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q);
1310:     PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq);
1311:   }
1312:   /* Create finite element */
1313:   PetscFECreate(comm, fem);
1314:   PetscFESetType(*fem, PETSCFEBASIC);
1315:   PetscFESetNumComponents(*fem, Nc);
1316:   PetscFESetBasisSpace(*fem, P);
1317:   PetscFESetDualSpace(*fem, Q);
1318:   PetscFESetQuadrature(*fem, q);
1319:   PetscFESetFaceQuadrature(*fem, fq);
1320:   if (prefix) {
1321:     PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix);
1322:     PetscFESetFromOptions(*fem);
1323:     PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL);
1324:   }
1325:   PetscFESetUp(*fem);
1326:   /* Cleanup */
1327:   PetscSpaceDestroy(&P);
1328:   PetscDualSpaceDestroy(&Q);
1329:   PetscQuadratureDestroy(&q);
1330:   PetscQuadratureDestroy(&fq);
1331:   /* Set finite element name */
1332:   PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k);
1333:   PetscFESetName(*fem, name);
1334:   return 0;
1335: }

1337: /*@C
1338:   DMPlexCreateGmshFromFile - Create a `DMPLEX` mesh from a Gmsh file

1340: + comm        - The MPI communicator
1341: . filename    - Name of the Gmsh file
1342: - interpolate - Create faces and edges in the mesh

1344:   Output Parameter:
1345: . dm  - The `DM` object representing the mesh

1347:   Level: beginner

1349: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
1350: @*/
1351: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1352: {
1353:   PetscViewer     viewer;
1354:   PetscMPIInt     rank;
1355:   int             fileType;
1356:   PetscViewerType vtype;

1358:   MPI_Comm_rank(comm, &rank);

1360:   /* Determine Gmsh file type (ASCII or binary) from file header */
1361:   if (rank == 0) {
1362:     GmshFile gmsh[1];
1363:     char     line[PETSC_MAX_PATH_LEN];
1364:     int      snum;
1365:     float    version;
1366:     int      fileFormat;

1368:     PetscArrayzero(gmsh, 1);
1369:     PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer);
1370:     PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII);
1371:     PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ);
1372:     PetscViewerFileSetName(gmsh->viewer, filename);
1373:     /* Read only the first two lines of the Gmsh file */
1374:     GmshReadSection(gmsh, line);
1375:     GmshExpect(gmsh, "$MeshFormat", line);
1376:     GmshReadString(gmsh, line, 2);
1377:     snum       = sscanf(line, "%f %d", &version, &fileType);
1378:     fileFormat = (int)roundf(version * 10);
1383:     PetscViewerDestroy(&gmsh->viewer);
1384:   }
1385:   MPI_Bcast(&fileType, 1, MPI_INT, 0, comm);
1386:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

1388:   /* Create appropriate viewer and build plex */
1389:   PetscViewerCreate(comm, &viewer);
1390:   PetscViewerSetType(viewer, vtype);
1391:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
1392:   PetscViewerFileSetName(viewer, filename);
1393:   DMPlexCreateGmsh(comm, viewer, interpolate, dm);
1394:   PetscViewerDestroy(&viewer);
1395:   return 0;
1396: }

1398: /*@
1399:   DMPlexCreateGmsh - Create a `DMPLEX` mesh from a Gmsh file viewer

1401:   Collective

1403:   Input Parameters:
1404: + comm  - The MPI communicator
1405: . viewer - The `PetscViewer` associated with a Gmsh file
1406: - interpolate - Create faces and edges in the mesh

1408:   Output Parameter:
1409: . dm  - The `DM` object representing the mesh

1411:   Options Database Keys:
1412: + -dm_plex_gmsh_hybrid        - Force triangular prisms to use tensor order
1413: . -dm_plex_gmsh_periodic      - Read Gmsh periodic section and construct a periodic Plex
1414: . -dm_plex_gmsh_highorder     - Generate high-order coordinates
1415: . -dm_plex_gmsh_project       - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
1416: . -dm_plex_gmsh_use_regions   - Generate labels with region names
1417: . -dm_plex_gmsh_mark_vertices - Add vertices to generated labels
1418: . -dm_plex_gmsh_multiple_tags - Allow multiple tags for default labels
1419: - -dm_plex_gmsh_spacedim <d>  - Embedding space dimension, if different from topological dimension

1421:   Note:
1422:   The Gmsh file format is described in http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format

1424:   By default, the "Cell Sets", "Face Sets", and "Vertex Sets" labels are created, and only insert the first tag on a point. By using -dm_plex_gmsh_multiple_tags, all tags can be inserted. Instead, -dm_plex_gmsh_use_regions creates labels based on the region names from the PhysicalNames section, and all tags are used.

1426:   Level: beginner

1428: .seealso: [](chapter_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1429: @*/
1430: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1431: {
1432:   GmshMesh    *mesh          = NULL;
1433:   PetscViewer  parentviewer  = NULL;
1434:   PetscBT      periodicVerts = NULL;
1435:   PetscBT      periodicCells = NULL;
1436:   DM           cdm, cdmCell = NULL;
1437:   PetscSection cs, csCell   = NULL;
1438:   Vec          coordinates, coordinatesCell;
1439:   DMLabel      cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1440:   PetscInt     dim = 0, coordDim = -1, order = 0;
1441:   PetscInt     numNodes = 0, numElems = 0, numVerts = 0, numCells = 0;
1442:   PetscInt     cell, cone[8], e, n, v, d;
1443:   PetscBool    binary, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE;
1444:   PetscBool    hybrid = interpolate, periodic = PETSC_TRUE;
1445:   PetscBool    highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1446:   PetscBool    isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1447:   PetscMPIInt  rank;

1449:   MPI_Comm_rank(comm, &rank);
1450:   PetscObjectOptionsBegin((PetscObject)viewer);
1451:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options");
1452:   PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL);
1453:   PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL);
1454:   PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet);
1455:   PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL);
1456:   PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL);
1457:   PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL);
1458:   PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL);
1459:   PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE);
1460:   PetscOptionsHeadEnd();
1461:   PetscOptionsEnd();

1463:   GmshCellInfoSetUp();

1465:   DMCreate(comm, dm);
1466:   DMSetType(*dm, DMPLEX);
1467:   PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL);

1469:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary);

1471:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1472:   if (binary) {
1473:     parentviewer = viewer;
1474:     PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);
1475:   }

1477:   if (rank == 0) {
1478:     GmshFile  gmsh[1];
1479:     char      line[PETSC_MAX_PATH_LEN];
1480:     PetscBool match;

1482:     PetscArrayzero(gmsh, 1);
1483:     gmsh->viewer = viewer;
1484:     gmsh->binary = binary;

1486:     GmshMeshCreate(&mesh);

1488:     /* Read mesh format */
1489:     GmshReadSection(gmsh, line);
1490:     GmshExpect(gmsh, "$MeshFormat", line);
1491:     GmshReadMeshFormat(gmsh);
1492:     GmshReadEndSection(gmsh, "$EndMeshFormat", line);

1494:     /* OPTIONAL Read physical names */
1495:     GmshReadSection(gmsh, line);
1496:     GmshMatch(gmsh, "$PhysicalNames", line, &match);
1497:     if (match) {
1498:       GmshExpect(gmsh, "$PhysicalNames", line);
1499:       GmshReadPhysicalNames(gmsh, mesh);
1500:       GmshReadEndSection(gmsh, "$EndPhysicalNames", line);
1501:       /* Initial read for entity section */
1502:       GmshReadSection(gmsh, line);
1503:     }

1505:     /* Read entities */
1506:     if (gmsh->fileFormat >= 40) {
1507:       GmshExpect(gmsh, "$Entities", line);
1508:       GmshReadEntities(gmsh, mesh);
1509:       GmshReadEndSection(gmsh, "$EndEntities", line);
1510:       /* Initial read for nodes section */
1511:       GmshReadSection(gmsh, line);
1512:     }

1514:     /* Read nodes */
1515:     GmshExpect(gmsh, "$Nodes", line);
1516:     GmshReadNodes(gmsh, mesh);
1517:     GmshReadEndSection(gmsh, "$EndNodes", line);

1519:     /* Read elements */
1520:     GmshReadSection(gmsh, line);
1521:     GmshExpect(gmsh, "$Elements", line);
1522:     GmshReadElements(gmsh, mesh);
1523:     GmshReadEndSection(gmsh, "$EndElements", line);

1525:     /* Read periodic section (OPTIONAL) */
1526:     if (periodic) {
1527:       GmshReadSection(gmsh, line);
1528:       GmshMatch(gmsh, "$Periodic", line, &periodic);
1529:     }
1530:     if (periodic) {
1531:       GmshExpect(gmsh, "$Periodic", line);
1532:       GmshReadPeriodic(gmsh, mesh);
1533:       GmshReadEndSection(gmsh, "$EndPeriodic", line);
1534:     }

1536:     PetscFree(gmsh->wbuf);
1537:     PetscFree(gmsh->sbuf);
1538:     PetscFree(gmsh->nbuf);

1540:     dim      = mesh->dim;
1541:     order    = mesh->order;
1542:     numNodes = mesh->numNodes;
1543:     numElems = mesh->numElems;
1544:     numVerts = mesh->numVerts;
1545:     numCells = mesh->numCells;

1547:     {
1548:       GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1549:       GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1550:       int          ptA   = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1551:       int          ptB   = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1552:       isSimplex          = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1553:       isHybrid           = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1554:       hasTetra           = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1555:     }
1556:   }

1558:   if (parentviewer) PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer);

1560:   {
1561:     int buf[6];

1563:     buf[0] = (int)dim;
1564:     buf[1] = (int)order;
1565:     buf[2] = periodic;
1566:     buf[3] = isSimplex;
1567:     buf[4] = isHybrid;
1568:     buf[5] = hasTetra;

1570:     MPI_Bcast(buf, 6, MPI_INT, 0, comm);

1572:     dim       = buf[0];
1573:     order     = buf[1];
1574:     periodic  = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1575:     isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1576:     isHybrid  = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1577:     hasTetra  = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1578:   }

1580:   if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;

1583:   /* We do not want this label automatically computed, instead we fill it here */
1584:   DMCreateLabel(*dm, "celltype");

1586:   /* Allocate the cell-vertex mesh */
1587:   DMPlexSetChart(*dm, 0, numCells + numVerts);
1588:   for (cell = 0; cell < numCells; ++cell) {
1589:     GmshElement   *elem  = mesh->elements + cell;
1590:     DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1591:     if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1592:     DMPlexSetConeSize(*dm, cell, elem->numVerts);
1593:     DMPlexSetCellType(*dm, cell, ctype);
1594:   }
1595:   for (v = numCells; v < numCells + numVerts; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
1596:   DMSetUp(*dm);

1598:   /* Add cell-vertex connections */
1599:   for (cell = 0; cell < numCells; ++cell) {
1600:     GmshElement *elem = mesh->elements + cell;
1601:     for (v = 0; v < elem->numVerts; ++v) {
1602:       const PetscInt nn = elem->nodes[v];
1603:       const PetscInt vv = mesh->vertexMap[nn];
1604:       cone[v]           = numCells + vv;
1605:     }
1606:     DMPlexReorderCell(*dm, cell, cone);
1607:     DMPlexSetCone(*dm, cell, cone);
1608:   }

1610:   DMSetDimension(*dm, dim);
1611:   DMPlexSymmetrize(*dm);
1612:   DMPlexStratify(*dm);
1613:   if (interpolate) {
1614:     DM idm;

1616:     DMPlexInterpolate(*dm, &idm);
1617:     DMDestroy(dm);
1618:     *dm = idm;
1619:   }

1621:   if (rank == 0) {
1622:     const PetscInt Nr = useregions ? mesh->numRegions : 0;
1623:     PetscInt       vStart, vEnd;

1625:     PetscCalloc1(Nr, &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:           PetscInt Nt = elem->numTags, t, r;

1635:           for (t = 0; t < Nt; ++t) {
1636:             const PetscInt  tag     = elem->tags[t];
1637:             const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1639:             if (generic) DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag);
1640:             for (r = 0; r < Nr; ++r) {
1641:               if (mesh->regionDims[r] != dim) continue;
1642:               if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], cell, tag);
1643:             }
1644:           }
1645:         }
1646:         cell++;
1647:       }

1649:       /* Create face sets */
1650:       if (interpolate && elem->dim == dim - 1) {
1651:         PetscInt        joinSize;
1652:         const PetscInt *join = NULL;
1653:         PetscInt        Nt   = elem->numTags, t, r;

1655:         /* Find the relevant facet with vertex joins */
1656:         for (v = 0; v < elem->numVerts; ++v) {
1657:           const PetscInt nn = elem->nodes[v];
1658:           const PetscInt vv = mesh->vertexMap[nn];
1659:           cone[v]           = vStart + vv;
1660:         }
1661:         DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1663:         for (t = 0; t < Nt; ++t) {
1664:           const PetscInt  tag     = elem->tags[t];
1665:           const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1667:           if (generic) DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag);
1668:           for (r = 0; r < Nr; ++r) {
1669:             if (mesh->regionDims[r] != dim - 1) continue;
1670:             if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], join[0], tag);
1671:           }
1672:         }
1673:         DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join);
1674:       }

1676:       /* Create vertex sets */
1677:       if (elem->dim == 0 && markvertices) {
1678:         if (elem->numTags > 0) {
1679:           const PetscInt nn  = elem->nodes[0];
1680:           const PetscInt vv  = mesh->vertexMap[nn];
1681:           const PetscInt tag = elem->tags[0];
1682:           PetscInt       r;

1684:           if (!Nr) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1685:           for (r = 0; r < Nr; ++r) {
1686:             if (mesh->regionDims[r] != 0) continue;
1687:             if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag);
1688:           }
1689:         }
1690:       }
1691:     }
1692:     if (markvertices) {
1693:       for (v = 0; v < numNodes; ++v) {
1694:         const PetscInt  vv   = mesh->vertexMap[v];
1695:         const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS];
1696:         PetscInt        r, t;

1698:         for (t = 0; t < GMSH_MAX_TAGS; ++t) {
1699:           const PetscInt  tag     = tags[t];
1700:           const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1702:           if (tag == -1) continue;
1703:           if (generic) DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag);
1704:           for (r = 0; r < Nr; ++r) {
1705:             if (mesh->regionTags[r] == tag) DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag);
1706:           }
1707:         }
1708:       }
1709:     }
1710:     PetscFree(regionSets);
1711:   }

1713:   { /* Create Cell/Face/Vertex Sets labels at all processes */
1714:     enum {
1715:       n = 4
1716:     };
1717:     PetscBool flag[n];

1719:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1720:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1721:     flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1722:     flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1723:     MPI_Bcast(flag, n, MPIU_BOOL, 0, comm);
1724:     if (flag[0]) DMCreateLabel(*dm, "Cell Sets");
1725:     if (flag[1]) DMCreateLabel(*dm, "Face Sets");
1726:     if (flag[2]) DMCreateLabel(*dm, "Vertex Sets");
1727:     if (flag[3]) DMCreateLabel(*dm, "marker");
1728:   }

1730:   if (periodic) {
1731:     PetscBTCreate(numVerts, &periodicVerts);
1732:     for (n = 0; n < numNodes; ++n) {
1733:       if (mesh->vertexMap[n] >= 0) {
1734:         if (PetscUnlikely(mesh->periodMap[n] != n)) {
1735:           PetscInt m = mesh->periodMap[n];
1736:           PetscBTSet(periodicVerts, mesh->vertexMap[n]);
1737:           PetscBTSet(periodicVerts, mesh->vertexMap[m]);
1738:         }
1739:       }
1740:     }
1741:     PetscBTCreate(numCells, &periodicCells);
1742:     for (cell = 0; cell < numCells; ++cell) {
1743:       GmshElement *elem = mesh->elements + cell;
1744:       for (v = 0; v < elem->numVerts; ++v) {
1745:         PetscInt nn = elem->nodes[v];
1746:         PetscInt vv = mesh->vertexMap[nn];
1747:         if (PetscUnlikely(PetscBTLookup(periodicVerts, vv))) {
1748:           PetscBTSet(periodicCells, cell);
1749:           break;
1750:         }
1751:       }
1752:     }
1753:   }

1755:   /* Setup coordinate DM */
1756:   if (coordDim < 0) coordDim = dim;
1757:   DMSetCoordinateDim(*dm, coordDim);
1758:   DMGetCoordinateDM(*dm, &cdm);
1759:   if (highOrder) {
1760:     PetscFE         fe;
1761:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1762:     PetscDTNodeType nodeType   = PETSCDTNODES_EQUISPACED;

1764:     if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */

1766:     GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1767:     PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view");
1768:     DMSetField(cdm, 0, NULL, (PetscObject)fe);
1769:     PetscFEDestroy(&fe);
1770:     DMCreateDS(cdm);
1771:   }
1772:   if (periodic) {
1773:     DMClone(cdm, &cdmCell);
1774:     DMSetCellCoordinateDM(*dm, cdmCell);
1775:   }

1777:   /* Create coordinates */
1778:   if (highOrder) {
1779:     PetscInt     maxDof = GmshNumNodes_HEX(order) * coordDim;
1780:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
1781:     PetscSection section;
1782:     PetscScalar *cellCoords;

1784:     DMSetLocalSection(cdm, NULL);
1785:     DMGetLocalSection(cdm, &cs);
1786:     PetscSectionClone(cs, &section);
1787:     DMPlexSetClosurePermutationTensor(cdm, 0, section); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */

1789:     DMCreateLocalVector(cdm, &coordinates);
1790:     PetscMalloc1(maxDof, &cellCoords);
1791:     for (cell = 0; cell < numCells; ++cell) {
1792:       GmshElement *elem     = mesh->elements + cell;
1793:       const int   *lexorder = GmshCellMap[elem->cellType].lexorder();
1794:       int          s        = 0;
1795:       for (n = 0; n < elem->numNodes; ++n) {
1796:         while (lexorder[n + s] < 0) ++s;
1797:         const PetscInt node = elem->nodes[lexorder[n + s]];
1798:         for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d];
1799:       }
1800:       if (s) {
1801:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1802:         PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1803:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1804:         PetscReal   hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1805:         PetscReal   hexFrontWeights[27]  = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1806:         PetscReal   hexLeftWeights[27]   = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0};
1807:         PetscReal   hexRightWeights[27]  = {0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25};
1808:         PetscReal   hexBackWeights[27]   = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25};
1809:         PetscReal   hexTopWeights[27]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1810:         PetscReal   hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25};
1811:         PetscReal  *sdWeights2[9]        = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
1812:         PetscReal  *sdWeights3[27]       = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
1813:                                             NULL, NULL, NULL, NULL, hexTopWeights,    NULL, NULL, NULL, NULL};
1814:         PetscReal **sdWeights[4]         = {NULL, NULL, sdWeights2, sdWeights3};

1816:         /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
1817:         for (n = 0; n < elem->numNodes + s; ++n) {
1818:           if (lexorder[n] >= 0) continue;
1819:           for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0;
1820:           for (int bn = 0; bn < elem->numNodes + s; ++bn) {
1821:             if (lexorder[bn] < 0) continue;
1822:             const PetscReal *weights = sdWeights[coordDim][n];
1823:             const PetscInt   bnode   = elem->nodes[lexorder[bn]];
1824:             for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d];
1825:           }
1826:         }
1827:       }
1828:       DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES);
1829:     }
1830:     PetscSectionDestroy(&section);
1831:     PetscFree(cellCoords);
1832:   } else {
1833:     PetscInt    *nodeMap;
1834:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
1835:     PetscScalar *pointCoords;

1837:     DMGetCoordinateSection(*dm, &cs);
1838:     PetscSectionSetNumFields(cs, 1);
1839:     PetscSectionSetFieldComponents(cs, 0, coordDim);
1840:     PetscSectionSetChart(cs, numCells, numCells + numVerts);
1841:     for (v = numCells; v < numCells + numVerts; ++v) {
1842:       PetscSectionSetDof(cs, v, coordDim);
1843:       PetscSectionSetFieldDof(cs, v, 0, coordDim);
1844:     }
1845:     PetscSectionSetUp(cs);

1847:     /* We need to localize coordinates on cells */
1848:     if (periodic) {
1849:       PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell);
1850:       PetscSectionSetNumFields(csCell, 1);
1851:       PetscSectionSetFieldComponents(csCell, 0, coordDim);
1852:       PetscSectionSetChart(csCell, 0, numCells);
1853:       for (cell = 0; cell < numCells; ++cell) {
1854:         if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1855:           GmshElement *elem = mesh->elements + cell;
1856:           PetscInt     dof  = elem->numVerts * coordDim;

1858:           PetscSectionSetDof(csCell, cell, dof);
1859:           PetscSectionSetFieldDof(csCell, cell, 0, dof);
1860:         }
1861:       }
1862:       PetscSectionSetUp(csCell);
1863:       DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell);
1864:     }

1866:     DMCreateLocalVector(cdm, &coordinates);
1867:     VecGetArray(coordinates, &pointCoords);
1868:     PetscMalloc1(numVerts, &nodeMap);
1869:     for (n = 0; n < numNodes; n++)
1870:       if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n;
1871:     for (v = 0; v < numVerts; ++v) {
1872:       PetscInt off, node = nodeMap[v];

1874:       PetscSectionGetOffset(cs, numCells + v, &off);
1875:       for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
1876:     }
1877:     VecRestoreArray(coordinates, &pointCoords);
1878:     PetscFree(nodeMap);

1880:     if (periodic) {
1881:       DMCreateLocalVector(cdmCell, &coordinatesCell);
1882:       VecGetArray(coordinatesCell, &pointCoords);
1883:       for (cell = 0; cell < numCells; ++cell) {
1884:         if (PetscUnlikely(PetscBTLookup(periodicCells, cell))) {
1885:           GmshElement *elem = mesh->elements + cell;
1886:           PetscInt     off, node;
1887:           for (v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v];
1888:           DMPlexReorderCell(cdmCell, cell, cone);
1889:           PetscSectionGetOffset(csCell, cell, &off);
1890:           for (v = 0; v < elem->numVerts; ++v)
1891:             for (node = cone[v], d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[node * 3 + d];
1892:         }
1893:       }
1894:       VecSetBlockSize(coordinatesCell, coordDim);
1895:       VecRestoreArray(coordinatesCell, &pointCoords);
1896:       DMSetCellCoordinatesLocal(*dm, coordinatesCell);
1897:       VecDestroy(&coordinatesCell);
1898:     }
1899:     PetscSectionDestroy(&csCell);
1900:     DMDestroy(&cdmCell);
1901:   }

1903:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
1904:   VecSetBlockSize(coordinates, coordDim);
1905:   DMSetCoordinatesLocal(*dm, coordinates);
1906:   VecDestroy(&coordinates);

1908:   GmshMeshDestroy(&mesh);
1909:   PetscBTDestroy(&periodicVerts);
1910:   PetscBTDestroy(&periodicCells);

1912:   if (highOrder && project) {
1913:     PetscFE         fe;
1914:     const char      prefix[]   = "dm_plex_gmsh_project_";
1915:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1916:     PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;

1918:     if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */

1920:     GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe);
1921:     PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view");
1922:     DMProjectCoordinates(*dm, fe);
1923:     PetscFEDestroy(&fe);
1924:   }

1926:   PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL);
1927:   return 0;
1928: }