Actual source code: plexgmsh.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

170: static GmshCellInfo GmshCellMap[150];

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

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

193: #define GmshCellTypeCheck(ct) \
194:   PetscMacroReturnStandard(const int _ct_ = (int)ct; PetscCheck(_ct_ >= 0 && _ct_ < (int)PETSC_STATIC_ARRAY_LENGTH(GmshCellMap), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid Gmsh element type %d", _ct_); PetscCheck(GmshCellMap[_ct_].cellType == _ct_, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_); \
195:                            PetscCheck(GmshCellMap[_ct_].polytope != -1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_);)

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

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

217:   PetscFunctionBegin;
218:   if (gmsh->wlen < size) {
219:     PetscCall(PetscFree(gmsh->wbuf));
220:     PetscCall(PetscMalloc(size, &gmsh->wbuf));
221:     gmsh->wlen = size;
222:   }
223:   *(void **)buf = size ? gmsh->wbuf : NULL;
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

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

232:   PetscFunctionBegin;
233:   if (gmsh->slen < size) {
234:     PetscCall(PetscFree(gmsh->sbuf));
235:     PetscCall(PetscMalloc(size, &gmsh->sbuf));
236:     gmsh->slen = size;
237:   }
238:   *(void **)buf = size ? gmsh->sbuf : NULL;
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscInt count, PetscDataType dtype)
243: {
244:   PetscFunctionBegin;
245:   PetscCall(PetscViewerRead(gmsh->viewer, buf, count, NULL, dtype));
246:   if (gmsh->byteSwap) PetscCall(PetscByteSwap(buf, dtype, count));
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscInt count)
251: {
252:   PetscFunctionBegin;
253:   PetscCall(PetscViewerRead(gmsh->viewer, buf, count, NULL, PETSC_STRING));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
258: {
259:   PetscFunctionBegin;
260:   PetscCall(PetscStrcmp(line, Section, match));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
265: {
266:   PetscBool match;

268:   PetscFunctionBegin;
269:   PetscCall(GmshMatch(gmsh, Section, line, &match));
270:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s", Section);
271:   PetscFunctionReturn(PETSC_SUCCESS);
272: }

274: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
275: {
276:   PetscBool match;

278:   PetscFunctionBegin;
279:   while (PETSC_TRUE) {
280:     PetscCall(GmshReadString(gmsh, line, 1));
281:     PetscCall(GmshMatch(gmsh, "$Comments", line, &match));
282:     if (!match) break;
283:     while (PETSC_TRUE) {
284:       PetscCall(GmshReadString(gmsh, line, 1));
285:       PetscCall(GmshMatch(gmsh, "$EndComments", line, &match));
286:       if (match) break;
287:     }
288:   }
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
293: {
294:   PetscFunctionBegin;
295:   PetscCall(GmshReadString(gmsh, line, 1));
296:   PetscCall(GmshExpect(gmsh, EndSection, line));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: static PetscErrorCode GmshReadSize(GmshFile *gmsh, PetscInt *buf, PetscInt count)
301: {
302:   PetscInt i;
303:   size_t   dataSize = (size_t)gmsh->dataSize;

305:   PetscFunctionBegin;
306:   if (dataSize == sizeof(PetscInt)) {
307:     PetscCall(GmshRead(gmsh, buf, count, PETSC_INT));
308:   } else if (dataSize == sizeof(int)) {
309:     int *ibuf = NULL;
310:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
311:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_ENUM));
312:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
313:   } else if (dataSize == sizeof(long)) {
314:     long *ibuf = NULL;
315:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
316:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_LONG));
317:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
318:   } else if (dataSize == sizeof(PetscInt64)) {
319:     PetscInt64 *ibuf = NULL;
320:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
321:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_INT64));
322:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
323:   }
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscInt count)
328: {
329:   PetscFunctionBegin;
330:   PetscCall(GmshRead(gmsh, buf, count, PETSC_ENUM));
331:   PetscFunctionReturn(PETSC_SUCCESS);
332: }

334: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscInt count)
335: {
336:   PetscFunctionBegin;
337:   PetscCall(GmshRead(gmsh, buf, count, PETSC_DOUBLE));
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: #define GMSH_MAX_TAGS 16

343: typedef struct {
344:   PetscInt id;                  /* Entity ID */
345:   PetscInt dim;                 /* Dimension */
346:   double   bbox[6];             /* Bounding box */
347:   PetscInt numTags;             /* Size of tag array */
348:   int      tags[GMSH_MAX_TAGS]; /* Tag array */
349: } GmshEntity;

351: typedef struct {
352:   GmshEntity *entity[4];
353:   PetscHMapI  entityMap[4];
354: } GmshEntities;

356: static PetscErrorCode GmshEntitiesCreate(PetscInt count[4], GmshEntities **entities)
357: {
358:   PetscInt dim;

360:   PetscFunctionBegin;
361:   PetscCall(PetscNew(entities));
362:   for (dim = 0; dim < 4; ++dim) {
363:     PetscCall(PetscCalloc1(count[dim], &(*entities)->entity[dim]));
364:     PetscCall(PetscHMapICreate(&(*entities)->entityMap[dim]));
365:   }
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
370: {
371:   PetscInt dim;

373:   PetscFunctionBegin;
374:   if (!*entities) PetscFunctionReturn(PETSC_SUCCESS);
375:   for (dim = 0; dim < 4; ++dim) {
376:     PetscCall(PetscFree((*entities)->entity[dim]));
377:     PetscCall(PetscHMapIDestroy(&(*entities)->entityMap[dim]));
378:   }
379:   PetscCall(PetscFree((*entities)));
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

383: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity **entity)
384: {
385:   PetscFunctionBegin;
386:   PetscCall(PetscHMapISet(entities->entityMap[dim], eid, index));
387:   entities->entity[dim][index].dim = dim;
388:   entities->entity[dim][index].id  = eid;
389:   if (entity) *entity = &entities->entity[dim][index];
390:   PetscFunctionReturn(PETSC_SUCCESS);
391: }

393: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity **entity)
394: {
395:   PetscInt index;

397:   PetscFunctionBegin;
398:   PetscCall(PetscHMapIGet(entities->entityMap[dim], eid, &index));
399:   *entity = &entities->entity[dim][index];
400:   PetscFunctionReturn(PETSC_SUCCESS);
401: }

403: typedef struct {
404:   PetscInt *id;  /* Node IDs */
405:   double   *xyz; /* Coordinates */
406:   PetscInt *tag; /* Physical tag */
407: } GmshNodes;

409: static PetscErrorCode GmshNodesCreate(PetscInt count, GmshNodes **nodes)
410: {
411:   PetscFunctionBegin;
412:   PetscCall(PetscNew(nodes));
413:   PetscCall(PetscMalloc1(count * 1, &(*nodes)->id));
414:   PetscCall(PetscMalloc1(count * 3, &(*nodes)->xyz));
415:   PetscCall(PetscMalloc1(count * GMSH_MAX_TAGS, &(*nodes)->tag));
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
420: {
421:   PetscFunctionBegin;
422:   if (!*nodes) PetscFunctionReturn(PETSC_SUCCESS);
423:   PetscCall(PetscFree((*nodes)->id));
424:   PetscCall(PetscFree((*nodes)->xyz));
425:   PetscCall(PetscFree((*nodes)->tag));
426:   PetscCall(PetscFree((*nodes)));
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: typedef struct {
431:   PetscInt  id;                  /* Element ID */
432:   PetscInt  dim;                 /* Dimension */
433:   PetscInt  cellType;            /* Cell type */
434:   PetscInt  numVerts;            /* Size of vertex array */
435:   PetscInt  numNodes;            /* Size of node array */
436:   PetscInt *nodes;               /* Vertex/Node array */
437:   PetscInt  numTags;             /* Size of physical tag array */
438:   int       tags[GMSH_MAX_TAGS]; /* Physical tag array */
439: } GmshElement;

441: static PetscErrorCode GmshElementsCreate(PetscInt count, GmshElement **elements)
442: {
443:   PetscFunctionBegin;
444:   PetscCall(PetscCalloc1(count, elements));
445:   PetscFunctionReturn(PETSC_SUCCESS);
446: }

448: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
449: {
450:   PetscFunctionBegin;
451:   if (!*elements) PetscFunctionReturn(PETSC_SUCCESS);
452:   PetscCall(PetscFree(*elements));
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: typedef struct {
457:   PetscInt       dim;
458:   PetscInt       order;
459:   GmshEntities  *entities;
460:   PetscInt       numNodes;
461:   GmshNodes     *nodelist;
462:   PetscInt       numElems;
463:   GmshElement   *elements;
464:   PetscInt       numVerts;
465:   PetscInt       numCells;
466:   PetscInt      *periodMap;
467:   PetscInt      *vertexMap;
468:   PetscSegBuffer segbuf;
469:   PetscInt       numRegions;
470:   PetscInt      *regionDims;
471:   PetscInt      *regionTags;
472:   char         **regionNames;
473: } GmshMesh;

475: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
476: {
477:   PetscFunctionBegin;
478:   PetscCall(PetscNew(mesh));
479:   PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf));
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

483: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
484: {
485:   PetscInt r;

487:   PetscFunctionBegin;
488:   if (!*mesh) PetscFunctionReturn(PETSC_SUCCESS);
489:   PetscCall(GmshEntitiesDestroy(&(*mesh)->entities));
490:   PetscCall(GmshNodesDestroy(&(*mesh)->nodelist));
491:   PetscCall(GmshElementsDestroy(&(*mesh)->elements));
492:   PetscCall(PetscFree((*mesh)->periodMap));
493:   PetscCall(PetscFree((*mesh)->vertexMap));
494:   PetscCall(PetscSegBufferDestroy(&(*mesh)->segbuf));
495:   for (r = 0; r < (*mesh)->numRegions; ++r) PetscCall(PetscFree((*mesh)->regionNames[r]));
496:   PetscCall(PetscFree3((*mesh)->regionDims, (*mesh)->regionTags, (*mesh)->regionNames));
497:   PetscCall(PetscFree((*mesh)));
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
502: {
503:   PetscViewer viewer   = gmsh->viewer;
504:   PetscBool   byteSwap = gmsh->byteSwap;
505:   char        line[PETSC_MAX_PATH_LEN];
506:   int         n, t, num, nid, snum;
507:   GmshNodes  *nodes;

509:   PetscFunctionBegin;
510:   PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
511:   snum = sscanf(line, "%d", &num);
512:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
513:   PetscCall(GmshNodesCreate(num, &nodes));
514:   mesh->numNodes = num;
515:   mesh->nodelist = nodes;
516:   for (n = 0; n < num; ++n) {
517:     double *xyz = nodes->xyz + n * 3;
518:     PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
519:     PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
520:     if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
521:     if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
522:     nodes->id[n] = nid;
523:     for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
524:   }
525:   PetscFunctionReturn(PETSC_SUCCESS);
526: }

528: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
529:    file contents multiple times to figure out the true number of cells and facets
530:    in the given mesh. To make this more efficient we read the file contents only
531:    once and store them in memory, while determining the true number of cells. */
532: static PetscErrorCode GmshReadElements_v22(GmshFile *gmsh, GmshMesh *mesh)
533: {
534:   PetscViewer  viewer   = gmsh->viewer;
535:   PetscBool    binary   = gmsh->binary;
536:   PetscBool    byteSwap = gmsh->byteSwap;
537:   char         line[PETSC_MAX_PATH_LEN];
538:   int          i, c, p, num, ibuf[1 + 4 + 1000], snum;
539:   int          cellType, numElem, numVerts, numNodes, numTags;
540:   GmshElement *elements;
541:   PetscInt    *nodeMap = gmsh->nodeMap;

543:   PetscFunctionBegin;
544:   PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
545:   snum = sscanf(line, "%d", &num);
546:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
547:   PetscCall(GmshElementsCreate(num, &elements));
548:   mesh->numElems = num;
549:   mesh->elements = elements;
550:   for (c = 0; c < num;) {
551:     PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
552:     if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));

554:     cellType = binary ? ibuf[0] : ibuf[1];
555:     numElem  = binary ? ibuf[1] : 1;
556:     numTags  = ibuf[2];

558:     PetscCall(GmshCellTypeCheck(cellType));
559:     numVerts = GmshCellMap[cellType].numVerts;
560:     numNodes = GmshCellMap[cellType].numNodes;

562:     for (i = 0; i < numElem; ++i, ++c) {
563:       GmshElement *element = elements + c;
564:       const int    off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
565:       const int   *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
566:       PetscCall(PetscViewerRead(viewer, ibuf + off, nint, NULL, PETSC_ENUM));
567:       if (byteSwap) PetscCall(PetscByteSwap(ibuf + off, PETSC_ENUM, nint));
568:       element->id       = id[0];
569:       element->dim      = GmshCellMap[cellType].dim;
570:       element->cellType = cellType;
571:       element->numVerts = numVerts;
572:       element->numNodes = numNodes;
573:       element->numTags  = PetscMin(numTags, GMSH_MAX_TAGS);
574:       PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
575:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
576:       for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
577:     }
578:   }
579:   PetscFunctionReturn(PETSC_SUCCESS);
580: }

582: /*
583: $Entities
584:   numPoints(unsigned long) numCurves(unsigned long)
585:     numSurfaces(unsigned long) numVolumes(unsigned long)
586:   // points
587:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
588:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
589:     numPhysicals(unsigned long) phyisicalTag[...](int)
590:   ...
591:   // curves
592:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
593:      boxMaxX(double) boxMaxY(double) boxMaxZ(double)
594:      numPhysicals(unsigned long) physicalTag[...](int)
595:      numBREPVert(unsigned long) tagBREPVert[...](int)
596:   ...
597:   // surfaces
598:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
599:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
600:     numPhysicals(unsigned long) physicalTag[...](int)
601:     numBREPCurve(unsigned long) tagBREPCurve[...](int)
602:   ...
603:   // volumes
604:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
605:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
606:     numPhysicals(unsigned long) physicalTag[...](int)
607:     numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
608:   ...
609: $EndEntities
610: */
611: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
612: {
613:   PetscViewer viewer   = gmsh->viewer;
614:   PetscBool   byteSwap = gmsh->byteSwap;
615:   long        index, num, lbuf[4];
616:   int         dim, eid, numTags, *ibuf, t;
617:   PetscInt    count[4], i;
618:   GmshEntity *entity = NULL;

620:   PetscFunctionBegin;
621:   PetscCall(PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG));
622:   if (byteSwap) PetscCall(PetscByteSwap(lbuf, PETSC_LONG, 4));
623:   for (i = 0; i < 4; ++i) count[i] = lbuf[i];
624:   PetscCall(GmshEntitiesCreate(count, &mesh->entities));
625:   for (dim = 0; dim < 4; ++dim) {
626:     for (index = 0; index < count[dim]; ++index) {
627:       PetscCall(PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM));
628:       if (byteSwap) PetscCall(PetscByteSwap(&eid, PETSC_ENUM, 1));
629:       PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
630:       PetscCall(PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE));
631:       if (byteSwap) PetscCall(PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6));
632:       PetscCall(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG));
633:       if (byteSwap) PetscCall(PetscByteSwap(&num, PETSC_LONG, 1));
634:       PetscCall(GmshBufferGet(gmsh, num, sizeof(int), &ibuf));
635:       PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
636:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
637:       entity->numTags = numTags = (int)PetscMin(num, GMSH_MAX_TAGS);
638:       for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
639:       if (dim == 0) continue;
640:       PetscCall(PetscViewerRead(viewer, &num, 1, NULL, PETSC_LONG));
641:       if (byteSwap) PetscCall(PetscByteSwap(&num, PETSC_LONG, 1));
642:       PetscCall(GmshBufferGet(gmsh, num, sizeof(int), &ibuf));
643:       PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
644:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
645:     }
646:   }
647:   PetscFunctionReturn(PETSC_SUCCESS);
648: }

650: /*
651: $Nodes
652:   numEntityBlocks(unsigned long) numNodes(unsigned long)
653:   tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
654:     tag(int) x(double) y(double) z(double)
655:     ...
656:   ...
657: $EndNodes
658: */
659: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
660: {
661:   PetscViewer viewer   = gmsh->viewer;
662:   PetscBool   byteSwap = gmsh->byteSwap;
663:   long        block, node, n, t, numEntityBlocks, numTotalNodes, numNodes;
664:   int         info[3], nid;
665:   GmshNodes  *nodes;

667:   PetscFunctionBegin;
668:   PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
669:   if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
670:   PetscCall(PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG));
671:   if (byteSwap) PetscCall(PetscByteSwap(&numTotalNodes, PETSC_LONG, 1));
672:   PetscCall(GmshNodesCreate(numTotalNodes, &nodes));
673:   mesh->numNodes = numTotalNodes;
674:   mesh->nodelist = nodes;
675:   for (n = 0, block = 0; block < numEntityBlocks; ++block) {
676:     PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
677:     PetscCall(PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG));
678:     if (byteSwap) PetscCall(PetscByteSwap(&numNodes, PETSC_LONG, 1));
679:     if (gmsh->binary) {
680:       size_t nbytes = sizeof(int) + 3 * sizeof(double);
681:       char  *cbuf   = NULL; /* dummy value to prevent warning from compiler about possible uninitialized value */
682:       PetscCall(GmshBufferGet(gmsh, numNodes, nbytes, &cbuf));
683:       PetscCall(PetscViewerRead(viewer, cbuf, numNodes * nbytes, NULL, PETSC_CHAR));
684:       for (node = 0; node < numNodes; ++node, ++n) {
685:         char   *cnid = cbuf + node * nbytes, *cxyz = cnid + sizeof(int);
686:         double *xyz = nodes->xyz + n * 3;
687:         if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cnid, PETSC_ENUM, 1));
688:         if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cxyz, PETSC_DOUBLE, 3));
689:         PetscCall(PetscMemcpy(&nid, cnid, sizeof(int)));
690:         PetscCall(PetscMemcpy(xyz, cxyz, 3 * sizeof(double)));
691:         if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
692:         if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
693:         nodes->id[n] = nid;
694:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
695:       }
696:     } else {
697:       for (node = 0; node < numNodes; ++node, ++n) {
698:         double *xyz = nodes->xyz + n * 3;
699:         PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
700:         PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
701:         if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
702:         if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
703:         nodes->id[n] = nid;
704:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
705:       }
706:     }
707:   }
708:   PetscFunctionReturn(PETSC_SUCCESS);
709: }

711: /*
712: $Elements
713:   numEntityBlocks(unsigned long) numElements(unsigned long)
714:   tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
715:     tag(int) numVert[...](int)
716:     ...
717:   ...
718: $EndElements
719: */
720: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
721: {
722:   PetscViewer  viewer   = gmsh->viewer;
723:   PetscBool    byteSwap = gmsh->byteSwap;
724:   long         c, block, numEntityBlocks, numTotalElements, elem, numElements;
725:   int          p, info[3], *ibuf = NULL;
726:   int          eid, dim, cellType, numVerts, numNodes, numTags;
727:   GmshEntity  *entity = NULL;
728:   GmshElement *elements;
729:   PetscInt    *nodeMap = gmsh->nodeMap;

731:   PetscFunctionBegin;
732:   PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
733:   if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
734:   PetscCall(PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG));
735:   if (byteSwap) PetscCall(PetscByteSwap(&numTotalElements, PETSC_LONG, 1));
736:   PetscCall(GmshElementsCreate(numTotalElements, &elements));
737:   mesh->numElems = numTotalElements;
738:   mesh->elements = elements;
739:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
740:     PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
741:     if (byteSwap) PetscCall(PetscByteSwap(info, PETSC_ENUM, 3));
742:     eid      = info[0];
743:     dim      = info[1];
744:     cellType = info[2];
745:     PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
746:     PetscCall(GmshCellTypeCheck(cellType));
747:     numVerts = GmshCellMap[cellType].numVerts;
748:     numNodes = GmshCellMap[cellType].numNodes;
749:     numTags  = entity->numTags;
750:     PetscCall(PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG));
751:     if (byteSwap) PetscCall(PetscByteSwap(&numElements, PETSC_LONG, 1));
752:     PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numElements, sizeof(int), &ibuf));
753:     PetscCall(PetscViewerRead(viewer, ibuf, (1 + numNodes) * numElements, NULL, PETSC_ENUM));
754:     if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, (1 + numNodes) * numElements));
755:     for (elem = 0; elem < numElements; ++elem, ++c) {
756:       GmshElement *element = elements + c;
757:       const int   *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
758:       element->id       = id[0];
759:       element->dim      = dim;
760:       element->cellType = cellType;
761:       element->numVerts = numVerts;
762:       element->numNodes = numNodes;
763:       element->numTags  = numTags;
764:       PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
765:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
766:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
767:     }
768:   }
769:   PetscFunctionReturn(PETSC_SUCCESS);
770: }

772: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
773: {
774:   PetscViewer viewer     = gmsh->viewer;
775:   int         fileFormat = gmsh->fileFormat;
776:   PetscBool   binary     = gmsh->binary;
777:   PetscBool   byteSwap   = gmsh->byteSwap;
778:   int         numPeriodic, snum, i;
779:   char        line[PETSC_MAX_PATH_LEN];
780:   PetscInt   *nodeMap = gmsh->nodeMap;

782:   PetscFunctionBegin;
783:   if (fileFormat == 22 || !binary) {
784:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
785:     snum = sscanf(line, "%d", &numPeriodic);
786:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
787:   } else {
788:     PetscCall(PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM));
789:     if (byteSwap) PetscCall(PetscByteSwap(&numPeriodic, PETSC_ENUM, 1));
790:   }
791:   for (i = 0; i < numPeriodic; i++) {
792:     int    ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
793:     long   j, nNodes;
794:     double affine[16];

796:     if (fileFormat == 22 || !binary) {
797:       PetscCall(PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING));
798:       snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
799:       PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
800:     } else {
801:       PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
802:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));
803:       correspondingDim = ibuf[0];
804:       correspondingTag = ibuf[1];
805:       primaryTag       = ibuf[2];
806:     }
807:     (void)correspondingDim;
808:     (void)correspondingTag;
809:     (void)primaryTag; /* unused */

811:     if (fileFormat == 22 || !binary) {
812:       PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
813:       snum = sscanf(line, "%ld", &nNodes);
814:       if (snum != 1) { /* discard transformation and try again */
815:         PetscCall(PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING));
816:         PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
817:         snum = sscanf(line, "%ld", &nNodes);
818:         PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
819:       }
820:     } else {
821:       PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
822:       if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
823:       if (nNodes == -1) { /* discard transformation and try again */
824:         PetscCall(PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE));
825:         PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
826:         if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
827:       }
828:     }

830:     for (j = 0; j < nNodes; j++) {
831:       if (fileFormat == 22 || !binary) {
832:         PetscCall(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING));
833:         snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
834:         PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
835:       } else {
836:         PetscCall(PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM));
837:         if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 2));
838:         correspondingNode = ibuf[0];
839:         primaryNode       = ibuf[1];
840:       }
841:       correspondingNode              = (int)nodeMap[correspondingNode];
842:       primaryNode                    = (int)nodeMap[primaryNode];
843:       periodicMap[correspondingNode] = primaryNode;
844:     }
845:   }
846:   PetscFunctionReturn(PETSC_SUCCESS);
847: }

849: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
850: $Entities
851:   numPoints(size_t) numCurves(size_t)
852:     numSurfaces(size_t) numVolumes(size_t)
853:   pointTag(int) X(double) Y(double) Z(double)
854:     numPhysicalTags(size_t) physicalTag(int) ...
855:   ...
856:   curveTag(int) minX(double) minY(double) minZ(double)
857:     maxX(double) maxY(double) maxZ(double)
858:     numPhysicalTags(size_t) physicalTag(int) ...
859:     numBoundingPoints(size_t) pointTag(int) ...
860:   ...
861:   surfaceTag(int) minX(double) minY(double) minZ(double)
862:     maxX(double) maxY(double) maxZ(double)
863:     numPhysicalTags(size_t) physicalTag(int) ...
864:     numBoundingCurves(size_t) curveTag(int) ...
865:   ...
866:   volumeTag(int) minX(double) minY(double) minZ(double)
867:     maxX(double) maxY(double) maxZ(double)
868:     numPhysicalTags(size_t) physicalTag(int) ...
869:     numBoundngSurfaces(size_t) surfaceTag(int) ...
870:   ...
871: $EndEntities
872: */
873: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
874: {
875:   PetscInt    count[4], index, numTags, i;
876:   int         dim, eid, *tags = NULL;
877:   GmshEntity *entity = NULL;

879:   PetscFunctionBegin;
880:   PetscCall(GmshReadSize(gmsh, count, 4));
881:   PetscCall(GmshEntitiesCreate(count, &mesh->entities));
882:   for (dim = 0; dim < 4; ++dim) {
883:     for (index = 0; index < count[dim]; ++index) {
884:       PetscCall(GmshReadInt(gmsh, &eid, 1));
885:       PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
886:       PetscCall(GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6));
887:       PetscCall(GmshReadSize(gmsh, &numTags, 1));
888:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
889:       PetscCall(GmshReadInt(gmsh, tags, numTags));
890:       PetscCheck(numTags <= GMSH_MAX_TAGS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PETSc currently supports up to %" PetscInt_FMT " tags per entity, not %" PetscInt_FMT, (PetscInt)GMSH_MAX_TAGS, numTags);
891:       entity->numTags = numTags;
892:       for (i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
893:       if (dim == 0) continue;
894:       PetscCall(GmshReadSize(gmsh, &numTags, 1));
895:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
896:       PetscCall(GmshReadInt(gmsh, tags, numTags));
897:       /* Currently, we do not save the ids for the bounding entities */
898:     }
899:   }
900:   PetscFunctionReturn(PETSC_SUCCESS);
901: }

903: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
904: $Nodes
905:   numEntityBlocks(size_t) numNodes(size_t)
906:     minNodeTag(size_t) maxNodeTag(size_t)
907:   entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
908:     nodeTag(size_t)
909:     ...
910:     x(double) y(double) z(double)
911:        < u(double; if parametric and entityDim = 1 or entityDim = 2) >
912:        < v(double; if parametric and entityDim = 2) >
913:     ...
914:   ...
915: $EndNodes
916: */
917: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
918: {
919:   int         info[3], dim, eid, parametric;
920:   PetscInt    sizes[4], numEntityBlocks, numTags, t, numNodes, numNodesBlock = 0, block, node, n;
921:   GmshEntity *entity = NULL;
922:   GmshNodes  *nodes;

924:   PetscFunctionBegin;
925:   PetscCall(GmshReadSize(gmsh, sizes, 4));
926:   numEntityBlocks = sizes[0];
927:   numNodes        = sizes[1];
928:   PetscCall(GmshNodesCreate(numNodes, &nodes));
929:   mesh->numNodes = numNodes;
930:   mesh->nodelist = nodes;
931:   for (block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
932:     PetscCall(GmshReadInt(gmsh, info, 3));
933:     dim        = info[0];
934:     eid        = info[1];
935:     parametric = info[2];
936:     PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
937:     numTags = entity->numTags;
938:     PetscCheck(!parametric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
939:     PetscCall(GmshReadSize(gmsh, &numNodesBlock, 1));
940:     PetscCall(GmshReadSize(gmsh, nodes->id + node, numNodesBlock));
941:     PetscCall(GmshReadDouble(gmsh, nodes->xyz + node * 3, numNodesBlock * 3));
942:     for (n = 0; n < numNodesBlock; ++n) {
943:       PetscInt *tags = &nodes->tag[node * GMSH_MAX_TAGS];

945:       for (t = 0; t < numTags; ++t) tags[n * GMSH_MAX_TAGS + t] = entity->tags[t];
946:       for (t = numTags; t < GMSH_MAX_TAGS; ++t) tags[n * GMSH_MAX_TAGS + t] = -1;
947:     }
948:   }
949:   gmsh->nodeStart = sizes[2];
950:   gmsh->nodeEnd   = sizes[3] + 1;
951:   PetscFunctionReturn(PETSC_SUCCESS);
952: }

954: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
955: $Elements
956:   numEntityBlocks(size_t) numElements(size_t)
957:     minElementTag(size_t) maxElementTag(size_t)
958:   entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
959:     elementTag(size_t) nodeTag(size_t) ...
960:     ...
961:   ...
962: $EndElements
963: */
964: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
965: {
966:   int          info[3], eid, dim, cellType;
967:   PetscInt     sizes[4], *ibuf = NULL, numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
968:   GmshEntity  *entity = NULL;
969:   GmshElement *elements;
970:   PetscInt    *nodeMap = gmsh->nodeMap;

972:   PetscFunctionBegin;
973:   PetscCall(GmshReadSize(gmsh, sizes, 4));
974:   numEntityBlocks = sizes[0];
975:   numElements     = sizes[1];
976:   PetscCall(GmshElementsCreate(numElements, &elements));
977:   mesh->numElems = numElements;
978:   mesh->elements = elements;
979:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
980:     PetscCall(GmshReadInt(gmsh, info, 3));
981:     dim      = info[0];
982:     eid      = info[1];
983:     cellType = info[2];
984:     PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
985:     PetscCall(GmshCellTypeCheck(cellType));
986:     numVerts = GmshCellMap[cellType].numVerts;
987:     numNodes = GmshCellMap[cellType].numNodes;
988:     numTags  = entity->numTags;
989:     PetscCall(GmshReadSize(gmsh, &numBlockElements, 1));
990:     PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numBlockElements, sizeof(PetscInt), &ibuf));
991:     PetscCall(GmshReadSize(gmsh, ibuf, (1 + numNodes) * numBlockElements));
992:     for (elem = 0; elem < numBlockElements; ++elem, ++c) {
993:       GmshElement    *element = elements + c;
994:       const PetscInt *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
995:       element->id       = id[0];
996:       element->dim      = dim;
997:       element->cellType = cellType;
998:       element->numVerts = numVerts;
999:       element->numNodes = numNodes;
1000:       element->numTags  = numTags;
1001:       PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
1002:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
1003:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
1004:     }
1005:   }
1006:   PetscFunctionReturn(PETSC_SUCCESS);
1007: }

1009: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1010: $Periodic
1011:   numPeriodicLinks(size_t)
1012:   entityDim(int) entityTag(int) entityTagPrimary(int)
1013:   numAffine(size_t) value(double) ...
1014:   numCorrespondingNodes(size_t)
1015:     nodeTag(size_t) nodeTagPrimary(size_t)
1016:     ...
1017:   ...
1018: $EndPeriodic
1019: */
1020: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
1021: {
1022:   int       info[3];
1023:   double    dbuf[16];
1024:   PetscInt  numPeriodicLinks, numAffine, numCorrespondingNodes, *nodeTags = NULL, link, node;
1025:   PetscInt *nodeMap = gmsh->nodeMap;

1027:   PetscFunctionBegin;
1028:   PetscCall(GmshReadSize(gmsh, &numPeriodicLinks, 1));
1029:   for (link = 0; link < numPeriodicLinks; ++link) {
1030:     PetscCall(GmshReadInt(gmsh, info, 3));
1031:     PetscCall(GmshReadSize(gmsh, &numAffine, 1));
1032:     PetscCall(GmshReadDouble(gmsh, dbuf, numAffine));
1033:     PetscCall(GmshReadSize(gmsh, &numCorrespondingNodes, 1));
1034:     PetscCall(GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags));
1035:     PetscCall(GmshReadSize(gmsh, nodeTags, numCorrespondingNodes * 2));
1036:     for (node = 0; node < numCorrespondingNodes; ++node) {
1037:       PetscInt correspondingNode     = nodeMap[nodeTags[node * 2 + 0]];
1038:       PetscInt primaryNode           = nodeMap[nodeTags[node * 2 + 1]];
1039:       periodicMap[correspondingNode] = primaryNode;
1040:     }
1041:   }
1042:   PetscFunctionReturn(PETSC_SUCCESS);
1043: }

1045: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1046: $MeshFormat // same as MSH version 2
1047:   version(ASCII double; currently 4.1)
1048:   file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1049:   data-size(ASCII int; sizeof(size_t))
1050:   < int with value one; only in binary mode, to detect endianness >
1051: $EndMeshFormat
1052: */
1053: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1054: {
1055:   char  line[PETSC_MAX_PATH_LEN];
1056:   int   snum, fileType, fileFormat, dataSize, checkEndian;
1057:   float version;

1059:   PetscFunctionBegin;
1060:   PetscCall(GmshReadString(gmsh, line, 3));
1061:   snum       = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1062:   fileFormat = (int)roundf(version * 10);
1063:   PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1064:   PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1065:   PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1066:   PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1067:   PetscCheck(!gmsh->binary || fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII");
1068:   PetscCheck(gmsh->binary || !fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary");
1069:   PetscCheck(fileFormat > 40 || dataSize == sizeof(double), PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1070:   PetscCheck(fileFormat < 41 || dataSize == sizeof(int) || dataSize == sizeof(PetscInt64), PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1071:   gmsh->fileFormat = fileFormat;
1072:   gmsh->dataSize   = dataSize;
1073:   gmsh->byteSwap   = PETSC_FALSE;
1074:   if (gmsh->binary) {
1075:     PetscCall(GmshReadInt(gmsh, &checkEndian, 1));
1076:     if (checkEndian != 1) {
1077:       PetscCall(PetscByteSwap(&checkEndian, PETSC_ENUM, 1));
1078:       PetscCheck(checkEndian == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line);
1079:       gmsh->byteSwap = PETSC_TRUE;
1080:     }
1081:   }
1082:   PetscFunctionReturn(PETSC_SUCCESS);
1083: }

1085: /*
1086: PhysicalNames
1087:   numPhysicalNames(ASCII int)
1088:   dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1089:   ...
1090: $EndPhysicalNames
1091: */
1092: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1093: {
1094:   char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p = NULL, *q = NULL, *r = NULL;
1095:   int  snum, region, dim, tag;

1097:   PetscFunctionBegin;
1098:   PetscCall(GmshReadString(gmsh, line, 1));
1099:   snum             = sscanf(line, "%d", &region);
1100:   mesh->numRegions = region;
1101:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1102:   PetscCall(PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames));
1103:   for (region = 0; region < mesh->numRegions; ++region) {
1104:     PetscCall(GmshReadString(gmsh, line, 2));
1105:     snum = sscanf(line, "%d %d", &dim, &tag);
1106:     PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1107:     PetscCall(GmshReadString(gmsh, line, -(PetscInt)sizeof(line)));
1108:     PetscCall(PetscStrchr(line, '"', &p));
1109:     PetscCheck(p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1110:     PetscCall(PetscStrrchr(line, '"', &q));
1111:     PetscCheck(q != p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1112:     PetscCall(PetscStrrchr(line, ':', &r));
1113:     if (p != r) q = r;
1114:     PetscCall(PetscStrncpy(name, p + 1, (size_t)(q - p - 1)));
1115:     mesh->regionDims[region] = dim;
1116:     mesh->regionTags[region] = tag;
1117:     PetscCall(PetscStrallocpy(name, &mesh->regionNames[region]));
1118:   }
1119:   PetscFunctionReturn(PETSC_SUCCESS);
1120: }

1122: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1123: {
1124:   PetscFunctionBegin;
1125:   switch (gmsh->fileFormat) {
1126:   case 41:
1127:     PetscCall(GmshReadEntities_v41(gmsh, mesh));
1128:     break;
1129:   default:
1130:     PetscCall(GmshReadEntities_v40(gmsh, mesh));
1131:     break;
1132:   }
1133:   PetscFunctionReturn(PETSC_SUCCESS);
1134: }

1136: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1137: {
1138:   PetscFunctionBegin;
1139:   switch (gmsh->fileFormat) {
1140:   case 41:
1141:     PetscCall(GmshReadNodes_v41(gmsh, mesh));
1142:     break;
1143:   case 40:
1144:     PetscCall(GmshReadNodes_v40(gmsh, mesh));
1145:     break;
1146:   default:
1147:     PetscCall(GmshReadNodes_v22(gmsh, mesh));
1148:     break;
1149:   }

1151:   { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1152:     if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1153:       PetscInt   tagMin = PETSC_MAX_INT, tagMax = PETSC_MIN_INT, n;
1154:       GmshNodes *nodes = mesh->nodelist;
1155:       for (n = 0; n < mesh->numNodes; ++n) {
1156:         const PetscInt tag = nodes->id[n];
1157:         tagMin             = PetscMin(tag, tagMin);
1158:         tagMax             = PetscMax(tag, tagMax);
1159:       }
1160:       gmsh->nodeStart = tagMin;
1161:       gmsh->nodeEnd   = tagMax + 1;
1162:     }
1163:   }

1165:   { /* Support for sparse node tags */
1166:     PetscInt   n, t;
1167:     GmshNodes *nodes = mesh->nodelist;
1168:     PetscCall(PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf));
1169:     for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_MIN_INT;
1170:     gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1171:     for (n = 0; n < mesh->numNodes; ++n) {
1172:       const PetscInt tag = nodes->id[n];
1173:       PetscCheck(gmsh->nodeMap[tag] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %" PetscInt_FMT, tag);
1174:       gmsh->nodeMap[tag] = n;
1175:     }
1176:   }
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1181: {
1182:   PetscFunctionBegin;
1183:   switch (gmsh->fileFormat) {
1184:   case 41:
1185:     PetscCall(GmshReadElements_v41(gmsh, mesh));
1186:     break;
1187:   case 40:
1188:     PetscCall(GmshReadElements_v40(gmsh, mesh));
1189:     break;
1190:   default:
1191:     PetscCall(GmshReadElements_v22(gmsh, mesh));
1192:     break;
1193:   }

1195:   { /* Reorder elements by codimension and polytope type */
1196:     PetscInt     ne       = mesh->numElems;
1197:     GmshElement *elements = mesh->elements;
1198:     PetscInt     keymap[GMSH_NUM_POLYTOPES], nk = 0;
1199:     PetscInt     offset[GMSH_NUM_POLYTOPES + 1], e, k;

1201:     for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_MIN_INT;
1202:     PetscCall(PetscMemzero(offset, sizeof(offset)));

1204:     keymap[GMSH_TET] = nk++;
1205:     keymap[GMSH_HEX] = nk++;
1206:     keymap[GMSH_PRI] = nk++;
1207:     keymap[GMSH_PYR] = nk++;
1208:     keymap[GMSH_TRI] = nk++;
1209:     keymap[GMSH_QUA] = nk++;
1210:     keymap[GMSH_SEG] = nk++;
1211:     keymap[GMSH_VTX] = nk++;

1213:     PetscCall(GmshElementsCreate(mesh->numElems, &mesh->elements));
1214: #define key(eid) keymap[GmshCellMap[elements[(eid)].cellType].polytope]
1215:     for (e = 0; e < ne; ++e) offset[1 + key(e)]++;
1216:     for (k = 1; k < nk; ++k) offset[k] += offset[k - 1];
1217:     for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1218: #undef key
1219:     PetscCall(GmshElementsDestroy(&elements));
1220:   }

1222:   { /* Mesh dimension and order */
1223:     GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1224:     mesh->dim         = elem ? GmshCellMap[elem->cellType].dim : 0;
1225:     mesh->order       = elem ? GmshCellMap[elem->cellType].order : 0;
1226:   }

1228:   {
1229:     PetscBT  vtx;
1230:     PetscInt dim = mesh->dim, e, n, v;

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

1234:     /* Compute number of cells and set of vertices */
1235:     mesh->numCells = 0;
1236:     for (e = 0; e < mesh->numElems; ++e) {
1237:       GmshElement *elem = mesh->elements + e;
1238:       if (elem->dim == dim && dim > 0) mesh->numCells++;
1239:       for (v = 0; v < elem->numVerts; v++) PetscCall(PetscBTSet(vtx, elem->nodes[v]));
1240:     }

1242:     /* Compute numbering for vertices */
1243:     mesh->numVerts = 0;
1244:     PetscCall(PetscMalloc1(mesh->numNodes, &mesh->vertexMap));
1245:     for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_MIN_INT;

1247:     PetscCall(PetscBTDestroy(&vtx));
1248:   }
1249:   PetscFunctionReturn(PETSC_SUCCESS);
1250: }

1252: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1253: {
1254:   PetscInt n;

1256:   PetscFunctionBegin;
1257:   PetscCall(PetscMalloc1(mesh->numNodes, &mesh->periodMap));
1258:   for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1259:   switch (gmsh->fileFormat) {
1260:   case 41:
1261:     PetscCall(GmshReadPeriodic_v41(gmsh, mesh->periodMap));
1262:     break;
1263:   default:
1264:     PetscCall(GmshReadPeriodic_v40(gmsh, mesh->periodMap));
1265:     break;
1266:   }

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

1272:   /* Renumber vertices (filter out correspondings) */
1273:   mesh->numVerts = 0;
1274:   for (n = 0; n < mesh->numNodes; ++n)
1275:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1276:       if (mesh->periodMap[n] == n) /* is primary */
1277:         mesh->vertexMap[n] = mesh->numVerts++;
1278:   for (n = 0; n < mesh->numNodes; ++n)
1279:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1280:       if (mesh->periodMap[n] != n) /* is corresponding  */
1281:         mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1282:   PetscFunctionReturn(PETSC_SUCCESS);
1283: }

1285: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1286: static const DMPolytopeType DMPolytopeMap[] = {
1287:   /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1288:   /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1289:   /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1290:   /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1291:   /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1292:   /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1293:   /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1294:   /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,       DM_POLYTOPE_UNKNOWN};

1296: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1297: {
1298:   return DMPolytopeMap[GmshCellMap[cellType].polytope];
1299: }

1301: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1302: {
1303:   DM              K;
1304:   PetscSpace      P;
1305:   PetscDualSpace  Q;
1306:   PetscQuadrature q, fq;
1307:   PetscBool       isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1308:   PetscBool       endpoint = PETSC_TRUE;
1309:   char            name[32];

1311:   PetscFunctionBegin;
1312:   /* Create space */
1313:   PetscCall(PetscSpaceCreate(comm, &P));
1314:   PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
1315:   PetscCall(PetscSpacePolynomialSetTensor(P, isTensor));
1316:   PetscCall(PetscSpaceSetNumComponents(P, Nc));
1317:   PetscCall(PetscSpaceSetNumVariables(P, dim));
1318:   PetscCall(PetscSpaceSetDegree(P, k, PETSC_DETERMINE));
1319:   if (prefix) {
1320:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
1321:     PetscCall(PetscSpaceSetFromOptions(P));
1322:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, NULL));
1323:     PetscCall(PetscSpaceGetDegree(P, &k, NULL));
1324:   }
1325:   PetscCall(PetscSpaceSetUp(P));
1326:   /* Create dual space */
1327:   PetscCall(PetscDualSpaceCreate(comm, &Q));
1328:   PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
1329:   PetscCall(PetscDualSpaceLagrangeSetTensor(Q, isTensor));
1330:   PetscCall(PetscDualSpaceLagrangeSetContinuity(Q, continuity));
1331:   PetscCall(PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0));
1332:   PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
1333:   PetscCall(PetscDualSpaceSetOrder(Q, k));
1334:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K));
1335:   PetscCall(PetscDualSpaceSetDM(Q, K));
1336:   PetscCall(DMDestroy(&K));
1337:   if (prefix) {
1338:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
1339:     PetscCall(PetscDualSpaceSetFromOptions(Q));
1340:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, NULL));
1341:   }
1342:   PetscCall(PetscDualSpaceSetUp(Q));
1343:   /* Create quadrature */
1344:   if (isSimplex) {
1345:     PetscCall(PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q));
1346:     PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1347:   } else {
1348:     PetscCall(PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q));
1349:     PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1350:   }
1351:   /* Create finite element */
1352:   PetscCall(PetscFECreate(comm, fem));
1353:   PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
1354:   PetscCall(PetscFESetNumComponents(*fem, Nc));
1355:   PetscCall(PetscFESetBasisSpace(*fem, P));
1356:   PetscCall(PetscFESetDualSpace(*fem, Q));
1357:   PetscCall(PetscFESetQuadrature(*fem, q));
1358:   PetscCall(PetscFESetFaceQuadrature(*fem, fq));
1359:   if (prefix) {
1360:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
1361:     PetscCall(PetscFESetFromOptions(*fem));
1362:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL));
1363:   }
1364:   PetscCall(PetscFESetUp(*fem));
1365:   /* Cleanup */
1366:   PetscCall(PetscSpaceDestroy(&P));
1367:   PetscCall(PetscDualSpaceDestroy(&Q));
1368:   PetscCall(PetscQuadratureDestroy(&q));
1369:   PetscCall(PetscQuadratureDestroy(&fq));
1370:   /* Set finite element name */
1371:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k));
1372:   PetscCall(PetscFESetName(*fem, name));
1373:   PetscFunctionReturn(PETSC_SUCCESS);
1374: }

1376: /*@C
1377:   DMPlexCreateGmshFromFile - Create a `DMPLEX` mesh from a Gmsh file

1379:   Input Parameters:
1380: + comm        - The MPI communicator
1381: . filename    - Name of the Gmsh file
1382: - interpolate - Create faces and edges in the mesh

1384:   Output Parameter:
1385: . dm - The `DM` object representing the mesh

1387:   Level: beginner

1389: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
1390: @*/
1391: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1392: {
1393:   PetscViewer     viewer;
1394:   PetscMPIInt     rank;
1395:   int             fileType;
1396:   PetscViewerType vtype;

1398:   PetscFunctionBegin;
1399:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

1401:   /* Determine Gmsh file type (ASCII or binary) from file header */
1402:   if (rank == 0) {
1403:     GmshFile gmsh[1];
1404:     char     line[PETSC_MAX_PATH_LEN];
1405:     int      snum;
1406:     float    version;
1407:     int      fileFormat;

1409:     PetscCall(PetscArrayzero(gmsh, 1));
1410:     PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer));
1411:     PetscCall(PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII));
1412:     PetscCall(PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ));
1413:     PetscCall(PetscViewerFileSetName(gmsh->viewer, filename));
1414:     /* Read only the first two lines of the Gmsh file */
1415:     PetscCall(GmshReadSection(gmsh, line));
1416:     PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1417:     PetscCall(GmshReadString(gmsh, line, 2));
1418:     snum       = sscanf(line, "%f %d", &version, &fileType);
1419:     fileFormat = (int)roundf(version * 10);
1420:     PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1421:     PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1422:     PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1423:     PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1424:     PetscCall(PetscViewerDestroy(&gmsh->viewer));
1425:   }
1426:   PetscCallMPI(MPI_Bcast(&fileType, 1, MPI_INT, 0, comm));
1427:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

1429:   /* Create appropriate viewer and build plex */
1430:   PetscCall(PetscViewerCreate(comm, &viewer));
1431:   PetscCall(PetscViewerSetType(viewer, vtype));
1432:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
1433:   PetscCall(PetscViewerFileSetName(viewer, filename));
1434:   PetscCall(DMPlexCreateGmsh(comm, viewer, interpolate, dm));
1435:   PetscCall(PetscViewerDestroy(&viewer));
1436:   PetscFunctionReturn(PETSC_SUCCESS);
1437: }

1439: /*@
1440:   DMPlexCreateGmsh - Create a `DMPLEX` mesh from a Gmsh file viewer

1442:   Collective

1444:   Input Parameters:
1445: + comm        - The MPI communicator
1446: . viewer      - The `PetscViewer` associated with a Gmsh file
1447: - interpolate - Create faces and edges in the mesh

1449:   Output Parameter:
1450: . dm - The `DM` object representing the mesh

1452:   Options Database Keys:
1453: + -dm_plex_gmsh_hybrid        - Force triangular prisms to use tensor order
1454: . -dm_plex_gmsh_periodic      - Read Gmsh periodic section and construct a periodic Plex
1455: . -dm_plex_gmsh_highorder     - Generate high-order coordinates
1456: . -dm_plex_gmsh_project       - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
1457: . -dm_plex_gmsh_use_regions   - Generate labels with region names
1458: . -dm_plex_gmsh_mark_vertices - Add vertices to generated labels
1459: . -dm_plex_gmsh_multiple_tags - Allow multiple tags for default labels
1460: - -dm_plex_gmsh_spacedim <d>  - Embedding space dimension, if different from topological dimension

1462:   Level: beginner

1464:   Notes:
1465:   The Gmsh file format is described in <http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format>

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

1469: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1470: @*/
1471: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1472: {
1473:   GmshMesh    *mesh          = NULL;
1474:   PetscViewer  parentviewer  = NULL;
1475:   PetscBT      periodicVerts = NULL;
1476:   PetscBT     *periodicCells = NULL;
1477:   DM           cdm, cdmCell = NULL;
1478:   PetscSection cs, csCell   = NULL;
1479:   Vec          coordinates, coordinatesCell;
1480:   DMLabel      cellSets = NULL, faceSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1481:   PetscInt     dim = 0, coordDim = -1, order = 0, maxHeight = 0;
1482:   PetscInt     numNodes = 0, numElems = 0, numVerts = 0, numCells = 0, vStart, vEnd;
1483:   PetscInt     cell, cone[8], e, n, v, d;
1484:   PetscBool    binary, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, multipleTags = PETSC_FALSE;
1485:   PetscBool    hybrid = interpolate, periodic = PETSC_TRUE;
1486:   PetscBool    highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1487:   PetscBool    isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1488:   PetscMPIInt  rank;

1490:   PetscFunctionBegin;
1491:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1492:   PetscObjectOptionsBegin((PetscObject)viewer);
1493:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options");
1494:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL));
1495:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL));
1496:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet));
1497:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL));
1498:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL));
1499:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL));
1500:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL));
1501:   PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE));
1502:   PetscCall(PetscOptionsBoundedInt("-dm_localize_height", "Localize edges and faces in addition to cells", "", maxHeight, &maxHeight, NULL, 0));
1503:   PetscOptionsHeadEnd();
1504:   PetscOptionsEnd();

1506:   PetscCall(GmshCellInfoSetUp());

1508:   PetscCall(DMCreate(comm, dm));
1509:   PetscCall(DMSetType(*dm, DMPLEX));
1510:   PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));

1512:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &binary));

1514:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1515:   if (binary) {
1516:     parentviewer = viewer;
1517:     PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer));
1518:   }

1520:   if (rank == 0) {
1521:     GmshFile  gmsh[1];
1522:     char      line[PETSC_MAX_PATH_LEN];
1523:     PetscBool match;

1525:     PetscCall(PetscArrayzero(gmsh, 1));
1526:     gmsh->viewer = viewer;
1527:     gmsh->binary = binary;

1529:     PetscCall(GmshMeshCreate(&mesh));

1531:     /* Read mesh format */
1532:     PetscCall(GmshReadSection(gmsh, line));
1533:     PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1534:     PetscCall(GmshReadMeshFormat(gmsh));
1535:     PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line));

1537:     /* OPTIONAL Read physical names */
1538:     PetscCall(GmshReadSection(gmsh, line));
1539:     PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match));
1540:     if (match) {
1541:       PetscCall(GmshExpect(gmsh, "$PhysicalNames", line));
1542:       PetscCall(GmshReadPhysicalNames(gmsh, mesh));
1543:       PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line));
1544:       /* Initial read for entity section */
1545:       PetscCall(GmshReadSection(gmsh, line));
1546:     }

1548:     /* Read entities */
1549:     if (gmsh->fileFormat >= 40) {
1550:       PetscCall(GmshExpect(gmsh, "$Entities", line));
1551:       PetscCall(GmshReadEntities(gmsh, mesh));
1552:       PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line));
1553:       /* Initial read for nodes section */
1554:       PetscCall(GmshReadSection(gmsh, line));
1555:     }

1557:     /* Read nodes */
1558:     PetscCall(GmshExpect(gmsh, "$Nodes", line));
1559:     PetscCall(GmshReadNodes(gmsh, mesh));
1560:     PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line));

1562:     /* Read elements */
1563:     PetscCall(GmshReadSection(gmsh, line));
1564:     PetscCall(GmshExpect(gmsh, "$Elements", line));
1565:     PetscCall(GmshReadElements(gmsh, mesh));
1566:     PetscCall(GmshReadEndSection(gmsh, "$EndElements", line));

1568:     /* Read periodic section (OPTIONAL) */
1569:     if (periodic) {
1570:       PetscCall(GmshReadSection(gmsh, line));
1571:       PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic));
1572:     }
1573:     if (periodic) {
1574:       PetscCall(GmshExpect(gmsh, "$Periodic", line));
1575:       PetscCall(GmshReadPeriodic(gmsh, mesh));
1576:       PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line));
1577:     }

1579:     PetscCall(PetscFree(gmsh->wbuf));
1580:     PetscCall(PetscFree(gmsh->sbuf));
1581:     PetscCall(PetscFree(gmsh->nbuf));

1583:     dim      = mesh->dim;
1584:     order    = mesh->order;
1585:     numNodes = mesh->numNodes;
1586:     numElems = mesh->numElems;
1587:     numVerts = mesh->numVerts;
1588:     numCells = mesh->numCells;

1590:     {
1591:       GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1592:       GmshElement *elemB = elemA ? elemA + mesh->numCells - 1 : NULL;
1593:       int          ptA   = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1594:       int          ptB   = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1595:       isSimplex          = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1596:       isHybrid           = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1597:       hasTetra           = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1598:     }
1599:   }

1601:   if (parentviewer) PetscCall(PetscViewerRestoreSubViewer(parentviewer, PETSC_COMM_SELF, &viewer));

1603:   {
1604:     int buf[6];

1606:     buf[0] = (int)dim;
1607:     buf[1] = (int)order;
1608:     buf[2] = periodic;
1609:     buf[3] = isSimplex;
1610:     buf[4] = isHybrid;
1611:     buf[5] = hasTetra;

1613:     PetscCallMPI(MPI_Bcast(buf, 6, MPI_INT, 0, comm));

1615:     dim       = buf[0];
1616:     order     = buf[1];
1617:     periodic  = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1618:     isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1619:     isHybrid  = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1620:     hasTetra  = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1621:   }

1623:   if (!highOrderSet) highOrder = (order > 1) ? PETSC_TRUE : PETSC_FALSE;
1624:   PetscCheck(!highOrder || !isHybrid, comm, PETSC_ERR_SUP, "No support for discretization on hybrid meshes yet");

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

1629:   /* Allocate the cell-vertex mesh */
1630:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numVerts));
1631:   for (cell = 0; cell < numCells; ++cell) {
1632:     GmshElement   *elem  = mesh->elements + cell;
1633:     DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1634:     if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1635:     PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts));
1636:     PetscCall(DMPlexSetCellType(*dm, cell, ctype));
1637:   }
1638:   for (v = numCells; v < numCells + numVerts; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1639:   PetscCall(DMSetUp(*dm));

1641:   /* Add cell-vertex connections */
1642:   for (cell = 0; cell < numCells; ++cell) {
1643:     GmshElement *elem = mesh->elements + cell;
1644:     for (v = 0; v < elem->numVerts; ++v) {
1645:       const PetscInt nn = elem->nodes[v];
1646:       const PetscInt vv = mesh->vertexMap[nn];
1647:       cone[v]           = numCells + vv;
1648:     }
1649:     PetscCall(DMPlexReorderCell(*dm, cell, cone));
1650:     PetscCall(DMPlexSetCone(*dm, cell, cone));
1651:   }

1653:   PetscCall(DMSetDimension(*dm, dim));
1654:   PetscCall(DMPlexSymmetrize(*dm));
1655:   PetscCall(DMPlexStratify(*dm));
1656:   if (interpolate) {
1657:     DM idm;

1659:     PetscCall(DMPlexInterpolate(*dm, &idm));
1660:     PetscCall(DMDestroy(dm));
1661:     *dm = idm;
1662:   }
1663:   PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd));

1665:   if (rank == 0) {
1666:     const PetscInt Nr = useregions ? mesh->numRegions : 0;

1668:     PetscCall(PetscCalloc1(Nr, &regionSets));
1669:     for (cell = 0, e = 0; e < numElems; ++e) {
1670:       GmshElement *elem = mesh->elements + e;

1672:       /* Create cell sets */
1673:       if (elem->dim == dim && dim > 0) {
1674:         if (elem->numTags > 0) {
1675:           PetscInt Nt = elem->numTags, t, r;

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

1681:             if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag));
1682:             for (r = 0; r < Nr; ++r) {
1683:               if (mesh->regionDims[r] != dim) continue;
1684:               if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], cell, tag));
1685:             }
1686:           }
1687:         }
1688:         cell++;
1689:       }

1691:       /* Create face sets */
1692:       if (elem->numTags && interpolate && elem->dim == dim - 1) {
1693:         PetscInt        joinSize;
1694:         const PetscInt *join = NULL;
1695:         PetscInt        Nt   = elem->numTags, t, r;

1697:         /* Find the relevant facet with vertex joins */
1698:         for (v = 0; v < elem->numVerts; ++v) {
1699:           const PetscInt nn = elem->nodes[v];
1700:           const PetscInt vv = mesh->vertexMap[nn];
1701:           cone[v]           = vStart + vv;
1702:         }
1703:         PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1704:         PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex facet for Gmsh element %" PetscInt_FMT " (Plex cell %" PetscInt_FMT ")", elem->id, e);
1705:         for (t = 0; t < Nt; ++t) {
1706:           const PetscInt  tag     = elem->tags[t];
1707:           const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1709:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag));
1710:           for (r = 0; r < Nr; ++r) {
1711:             if (mesh->regionDims[r] != dim - 1) continue;
1712:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], join[0], tag));
1713:           }
1714:         }
1715:         PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1716:       }

1718:       /* Create edge sets */
1719:       if (elem->numTags && interpolate && dim > 2 && elem->dim == 1) {
1720:         PetscInt        joinSize;
1721:         const PetscInt *join = NULL;
1722:         PetscInt        Nt   = elem->numTags, t, r;

1724:         /* Find the relevant edge with vertex joins */
1725:         for (v = 0; v < elem->numVerts; ++v) {
1726:           const PetscInt nn = elem->nodes[v];
1727:           const PetscInt vv = mesh->vertexMap[nn];
1728:           cone[v]           = vStart + vv;
1729:         }
1730:         PetscCall(DMPlexGetJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1731:         PetscCheck(joinSize == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Could not determine Plex edge for Gmsh element %" PetscInt_FMT " (Plex cell %" PetscInt_FMT ")", elem->id, e);
1732:         for (t = 0; t < Nt; ++t) {
1733:           const PetscInt  tag     = elem->tags[t];
1734:           const PetscBool generic = !Nr && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1736:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Edge Sets", join[0], tag));
1737:           for (r = 0; r < Nr; ++r) {
1738:             if (mesh->regionDims[r] != 1) continue;
1739:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], join[0], tag));
1740:           }
1741:         }
1742:         PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1743:       }

1745:       /* Create vertex sets */
1746:       if (elem->numTags && elem->dim == 0 && markvertices) {
1747:         const PetscInt nn  = elem->nodes[0];
1748:         const PetscInt vv  = mesh->vertexMap[nn];
1749:         const PetscInt tag = elem->tags[0];
1750:         PetscInt       r;

1752:         if (!Nr) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1753:         for (r = 0; r < Nr; ++r) {
1754:           if (mesh->regionDims[r] != 0) continue;
1755:           if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag));
1756:         }
1757:       }
1758:     }
1759:     if (markvertices) {
1760:       for (v = 0; v < numNodes; ++v) {
1761:         const PetscInt  vv   = mesh->vertexMap[v];
1762:         const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS];
1763:         PetscInt        r, t;

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

1769:           if (tag == -1) continue;
1770:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1771:           for (r = 0; r < Nr; ++r) {
1772:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag));
1773:           }
1774:         }
1775:       }
1776:     }
1777:     PetscCall(PetscFree(regionSets));
1778:   }

1780:   { /* Create Cell/Face/Vertex Sets labels at all processes */
1781:     enum {
1782:       n = 4
1783:     };
1784:     PetscBool flag[n];

1786:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1787:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1788:     flag[2] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1789:     flag[3] = marker ? PETSC_TRUE : PETSC_FALSE;
1790:     PetscCallMPI(MPI_Bcast(flag, n, MPIU_BOOL, 0, comm));
1791:     if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1792:     if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1793:     if (flag[2]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1794:     if (flag[3]) PetscCall(DMCreateLabel(*dm, "marker"));
1795:   }

1797:   if (periodic) {
1798:     PetscCall(PetscBTCreate(numVerts, &periodicVerts));
1799:     for (n = 0; n < numNodes; ++n) {
1800:       if (mesh->vertexMap[n] >= 0) {
1801:         if (PetscUnlikely(mesh->periodMap[n] != n)) {
1802:           PetscInt m = mesh->periodMap[n];
1803:           PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n]));
1804:           PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m]));
1805:         }
1806:       }
1807:     }
1808:     PetscCall(DMGetCoordinateDM(*dm, &cdm));
1809:     PetscCall(PetscMalloc1(maxHeight + 1, &periodicCells));
1810:     for (PetscInt h = 0; h <= maxHeight; ++h) {
1811:       PetscInt pStart, pEnd;

1813:       PetscCall(DMPlexGetHeightStratum(*dm, h, &pStart, &pEnd));
1814:       PetscCall(PetscBTCreate(pEnd - pStart, &periodicCells[h]));
1815:       for (PetscInt p = pStart; p < pEnd; ++p) {
1816:         PetscInt *closure = NULL;
1817:         PetscInt  Ncl;

1819:         PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1820:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1821:           if (closure[cl] >= vStart && closure[cl] < vEnd) {
1822:             if (PetscUnlikely(PetscBTLookup(periodicVerts, closure[cl] - vStart))) {
1823:               PetscCall(PetscBTSet(periodicCells[h], p - pStart));
1824:               break;
1825:             }
1826:           }
1827:         }
1828:         PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1829:       }
1830:     }
1831:   }

1833:   /* Setup coordinate DM */
1834:   if (coordDim < 0) coordDim = dim;
1835:   PetscCall(DMSetCoordinateDim(*dm, coordDim));
1836:   PetscCall(DMGetCoordinateDM(*dm, &cdm));
1837:   if (highOrder) {
1838:     PetscFE         fe;
1839:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1840:     PetscDTNodeType nodeType   = PETSCDTNODES_EQUISPACED;

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

1844:     PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
1845:     PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view"));
1846:     PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe));
1847:     PetscCall(PetscFEDestroy(&fe));
1848:     PetscCall(DMCreateDS(cdm));
1849:   }
1850:   if (periodic) {
1851:     PetscCall(DMClone(cdm, &cdmCell));
1852:     PetscCall(DMSetCellCoordinateDM(*dm, cdmCell));
1853:   }

1855:   /* Create coordinates */
1856:   if (highOrder) {
1857:     PetscInt     maxDof = GmshNumNodes_HEX(order) * coordDim;
1858:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
1859:     PetscSection section;
1860:     PetscScalar *cellCoords;

1862:     PetscCall(DMSetLocalSection(cdm, NULL));
1863:     PetscCall(DMGetLocalSection(cdm, &cs));
1864:     PetscCall(PetscSectionClone(cs, &section));
1865:     PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */

1867:     PetscCall(DMCreateLocalVector(cdm, &coordinates));
1868:     PetscCall(PetscMalloc1(maxDof, &cellCoords));
1869:     for (cell = 0; cell < numCells; ++cell) {
1870:       GmshElement *elem     = mesh->elements + cell;
1871:       const int   *lexorder = GmshCellMap[elem->cellType].lexorder();
1872:       int          s        = 0;
1873:       for (n = 0; n < elem->numNodes; ++n) {
1874:         while (lexorder[n + s] < 0) ++s;
1875:         const PetscInt node = elem->nodes[lexorder[n + s]];
1876:         for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d];
1877:       }
1878:       if (s) {
1879:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1880:         PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1881:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
1882:         PetscReal   hexBottomWeights[27] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1883:         PetscReal   hexFrontWeights[27]  = {-0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1884:         PetscReal   hexLeftWeights[27]   = {-0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0};
1885:         PetscReal   hexRightWeights[27]  = {0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25, 0.0, 0.0, 0.5, 0.0, 0.0, -0.25};
1886:         PetscReal   hexBackWeights[27]   = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25};
1887:         PetscReal   hexTopWeights[27]    = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
1888:         PetscReal   hexCenterWeights[27] = {-0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, 0.0, 0.0, 0.0, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25, 0.25, 0.0, 0.25, -0.25, 0.25, -0.25};
1889:         PetscReal  *sdWeights2[9]        = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
1890:         PetscReal  *sdWeights3[27]       = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
1891:                                             NULL, NULL, NULL, NULL, hexTopWeights,    NULL, NULL, NULL, NULL};
1892:         PetscReal **sdWeights[4]         = {NULL, NULL, sdWeights2, sdWeights3};

1894:         /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
1895:         for (n = 0; n < elem->numNodes + s; ++n) {
1896:           if (lexorder[n] >= 0) continue;
1897:           for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0;
1898:           for (int bn = 0; bn < elem->numNodes + s; ++bn) {
1899:             if (lexorder[bn] < 0) continue;
1900:             const PetscReal *weights = sdWeights[coordDim][n];
1901:             const PetscInt   bnode   = elem->nodes[lexorder[bn]];
1902:             for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d];
1903:           }
1904:         }
1905:       }
1906:       PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES));
1907:     }
1908:     PetscCall(PetscSectionDestroy(&section));
1909:     PetscCall(PetscFree(cellCoords));
1910:   } else {
1911:     PetscInt    *nodeMap;
1912:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
1913:     PetscScalar *pointCoords;

1915:     PetscCall(DMGetCoordinateSection(*dm, &cs));
1916:     PetscCall(PetscSectionSetNumFields(cs, 1));
1917:     PetscCall(PetscSectionSetFieldComponents(cs, 0, coordDim));
1918:     PetscCall(PetscSectionSetChart(cs, numCells, numCells + numVerts));
1919:     for (v = numCells; v < numCells + numVerts; ++v) {
1920:       PetscCall(PetscSectionSetDof(cs, v, coordDim));
1921:       PetscCall(PetscSectionSetFieldDof(cs, v, 0, coordDim));
1922:     }
1923:     PetscCall(PetscSectionSetUp(cs));

1925:     /* We need to localize coordinates on cells */
1926:     if (periodic) {
1927:       PetscInt newStart = PETSC_MAX_INT, newEnd = -1, pStart, pEnd;

1929:       PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell));
1930:       PetscCall(PetscSectionSetNumFields(csCell, 1));
1931:       PetscCall(PetscSectionSetFieldComponents(csCell, 0, coordDim));
1932:       for (PetscInt h = 0; h <= maxHeight; h++) {
1933:         PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
1934:         newStart = PetscMin(newStart, pStart);
1935:         newEnd   = PetscMax(newEnd, pEnd);
1936:       }
1937:       PetscCall(PetscSectionSetChart(csCell, newStart, newEnd));
1938:       for (PetscInt h = 0; h <= maxHeight; h++) {
1939:         PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
1940:         for (PetscInt p = pStart; p < pEnd; ++p) {
1941:           PetscInt *closure = NULL;
1942:           PetscInt  Ncl, Nv = 0;

1944:           if (PetscUnlikely(PetscBTLookup(periodicCells[h], p - pStart))) {
1945:             PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1946:             for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1947:               if (closure[cl] >= vStart && closure[cl] < vEnd) ++Nv;
1948:             }
1949:             PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1950:             PetscCall(PetscSectionSetDof(csCell, p, Nv * coordDim));
1951:             PetscCall(PetscSectionSetFieldDof(csCell, p, 0, Nv * coordDim));
1952:           }
1953:         }
1954:       }
1955:       PetscCall(PetscSectionSetUp(csCell));
1956:       PetscCall(DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell));
1957:     }

1959:     PetscCall(DMCreateLocalVector(cdm, &coordinates));
1960:     PetscCall(VecGetArray(coordinates, &pointCoords));
1961:     PetscCall(PetscMalloc1(numVerts, &nodeMap));
1962:     for (n = 0; n < numNodes; n++)
1963:       if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n;
1964:     for (v = 0; v < numVerts; ++v) {
1965:       PetscInt off, node = nodeMap[v];

1967:       PetscCall(PetscSectionGetOffset(cs, numCells + v, &off));
1968:       for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
1969:     }
1970:     PetscCall(VecRestoreArray(coordinates, &pointCoords));

1972:     if (periodic) {
1973:       PetscInt cStart, cEnd;

1975:       PetscCall(DMPlexGetHeightStratum(cdmCell, 0, &cStart, &cEnd));
1976:       PetscCall(DMCreateLocalVector(cdmCell, &coordinatesCell));
1977:       PetscCall(VecGetArray(coordinatesCell, &pointCoords));
1978:       for (PetscInt c = cStart; c < cEnd; ++c) {
1979:         GmshElement *elem    = mesh->elements + c - cStart;
1980:         PetscInt    *closure = NULL;
1981:         PetscInt     verts[8];
1982:         PetscInt     Ncl, Nv = 0;

1984:         for (PetscInt v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v];
1985:         PetscCall(DMPlexReorderCell(cdmCell, c, cone));
1986:         PetscCall(DMPlexGetTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
1987:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1988:           if (closure[cl] >= vStart && closure[cl] < vEnd) verts[Nv++] = closure[cl];
1989:         }
1990:         PetscCheck(Nv == elem->numVerts, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of vertices %" PetscInt_FMT " in closure does not match number of vertices %" PetscInt_FMT " in Gmsh cell", Nv, elem->numVerts);
1991:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1992:           const PetscInt point = closure[cl];
1993:           PetscInt       hStart, h;

1995:           PetscCall(DMPlexGetPointHeight(cdmCell, point, &h));
1996:           if (h > maxHeight) continue;
1997:           PetscCall(DMPlexGetHeightStratum(cdmCell, h, &hStart, NULL));
1998:           if (PetscUnlikely(PetscBTLookup(periodicCells[h], point - hStart))) {
1999:             PetscInt *pclosure = NULL;
2000:             PetscInt  Npcl, off, v;

2002:             PetscCall(PetscSectionGetOffset(csCell, point, &off));
2003:             PetscCall(DMPlexGetTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2004:             for (PetscInt pcl = 0; pcl < Npcl * 2; pcl += 2) {
2005:               if (pclosure[pcl] >= vStart && pclosure[pcl] < vEnd) {
2006:                 for (v = 0; v < Nv; ++v)
2007:                   if (verts[v] == pclosure[pcl]) break;
2008:                 PetscCheck(v < Nv, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not find vertex %" PetscInt_FMT " in closure of cell %" PetscInt_FMT, pclosure[pcl], c);
2009:                 for (PetscInt d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[cone[v] * 3 + d];
2010:                 ++v;
2011:               }
2012:             }
2013:             PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2014:           }
2015:         }
2016:         PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
2017:       }
2018:       PetscCall(VecSetBlockSize(coordinatesCell, coordDim));
2019:       PetscCall(VecRestoreArray(coordinatesCell, &pointCoords));
2020:       PetscCall(DMSetCellCoordinatesLocal(*dm, coordinatesCell));
2021:       PetscCall(VecDestroy(&coordinatesCell));
2022:     }
2023:     PetscCall(PetscFree(nodeMap));
2024:     PetscCall(PetscSectionDestroy(&csCell));
2025:     PetscCall(DMDestroy(&cdmCell));
2026:   }

2028:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
2029:   PetscCall(VecSetBlockSize(coordinates, coordDim));
2030:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
2031:   PetscCall(VecDestroy(&coordinates));

2033:   PetscCall(GmshMeshDestroy(&mesh));
2034:   PetscCall(PetscBTDestroy(&periodicVerts));
2035:   if (periodic) {
2036:     for (PetscInt h = 0; h <= maxHeight; ++h) PetscCall(PetscBTDestroy(periodicCells + h));
2037:     PetscCall(PetscFree(periodicCells));
2038:   }

2040:   if (highOrder && project) {
2041:     PetscFE         fe;
2042:     const char      prefix[]   = "dm_plex_gmsh_project_";
2043:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
2044:     PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;

2046:     if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
2047:     PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
2048:     PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view"));
2049:     PetscCall(DMProjectCoordinates(*dm, fe));
2050:     PetscCall(PetscFEDestroy(&fe));
2051:   }

2053:   PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));
2054:   PetscFunctionReturn(PETSC_SUCCESS);
2055: }