Actual source code: plexgmsh.c

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

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

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

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

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

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

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

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

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

119: #define GmshCellEntry(cellType, polytope, dim, order) {cellType, GMSH_##polytope, dim, order, GmshNumNodes_##polytope(1), GmshNumNodes_##polytope(order), GmshLexOrder_##polytope##_##order}

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

124:   GmshCellEntry(1, SEG, 1, 1),   GmshCellEntry(8, SEG, 1, 2),   GmshCellEntry(26, SEG, 1, 3),
125:   GmshCellEntry(27, SEG, 1, 4),  GmshCellEntry(28, SEG, 1, 5),  GmshCellEntry(62, SEG, 1, 6),
126:   GmshCellEntry(63, SEG, 1, 7),  GmshCellEntry(64, SEG, 1, 8),  GmshCellEntry(65, SEG, 1, 9),
127:   GmshCellEntry(66, SEG, 1, 10),

129:   GmshCellEntry(2, TRI, 2, 1),   GmshCellEntry(9, TRI, 2, 2),   GmshCellEntry(21, TRI, 2, 3),
130:   GmshCellEntry(23, TRI, 2, 4),  GmshCellEntry(25, TRI, 2, 5),  GmshCellEntry(42, TRI, 2, 6),
131:   GmshCellEntry(43, TRI, 2, 7),  GmshCellEntry(44, TRI, 2, 8),  GmshCellEntry(45, TRI, 2, 9),
132:   GmshCellEntry(46, TRI, 2, 10),

134:   GmshCellEntry(3, QUA, 2, 1),   GmshCellEntry(10, QUA, 2, 2),  {16, GMSH_QUA, 2, 2, 4, 8,  GmshLexOrder_QUA_2_Serendipity},
135:   GmshCellEntry(36, QUA, 2, 3),  GmshCellEntry(37, QUA, 2, 4),  GmshCellEntry(38, QUA, 2, 5),
136:   GmshCellEntry(47, QUA, 2, 6),  GmshCellEntry(48, QUA, 2, 7),  GmshCellEntry(49, QUA, 2, 8),
137:   GmshCellEntry(50, QUA, 2, 9),  GmshCellEntry(51, QUA, 2, 10),

139:   GmshCellEntry(4, TET, 3, 1),   GmshCellEntry(11, TET, 3, 2),  GmshCellEntry(29, TET, 3, 3),
140:   GmshCellEntry(30, TET, 3, 4),  GmshCellEntry(31, TET, 3, 5),  GmshCellEntry(71, TET, 3, 6),
141:   GmshCellEntry(72, TET, 3, 7),  GmshCellEntry(73, TET, 3, 8),  GmshCellEntry(74, TET, 3, 9),
142:   GmshCellEntry(75, TET, 3, 10),

144:   GmshCellEntry(5, HEX, 3, 1),   GmshCellEntry(12, HEX, 3, 2),  {17, GMSH_HEX, 3, 2, 8, 20, GmshLexOrder_HEX_2_Serendipity},
145:   GmshCellEntry(92, HEX, 3, 3),  GmshCellEntry(93, HEX, 3, 4),  GmshCellEntry(94, HEX, 3, 5),
146:   GmshCellEntry(95, HEX, 3, 6),  GmshCellEntry(96, HEX, 3, 7),  GmshCellEntry(97, HEX, 3, 8),
147:   GmshCellEntry(98, HEX, 3, 9),  GmshCellEntry(-1, HEX, 3, 10),

149:   GmshCellEntry(6, PRI, 3, 1),   GmshCellEntry(13, PRI, 3, 2),  GmshCellEntry(90, PRI, 3, 3),
150:   GmshCellEntry(91, PRI, 3, 4),  GmshCellEntry(106, PRI, 3, 5), GmshCellEntry(107, PRI, 3, 6),
151:   GmshCellEntry(108, PRI, 3, 7), GmshCellEntry(109, PRI, 3, 8), GmshCellEntry(110, PRI, 3, 9),
152:   GmshCellEntry(-1, PRI, 3, 10),

154:   GmshCellEntry(7, PYR, 3, 1),   GmshCellEntry(14, PYR, 3, 2),  GmshCellEntry(118, PYR, 3, 3),
155:   GmshCellEntry(119, PYR, 3, 4), GmshCellEntry(120, PYR, 3, 5), GmshCellEntry(121, PYR, 3, 6),
156:   GmshCellEntry(122, PYR, 3, 7), GmshCellEntry(123, PYR, 3, 8), GmshCellEntry(124, PYR, 3, 9),
157:   GmshCellEntry(-1, PYR, 3, 10)

159: #if 0
160:   {20, GMSH_TRI, 2, 3, 3,  9, NULL},
161:   {18, GMSH_PRI, 3, 2, 6, 15, NULL},
162:   {19, GMSH_PYR, 3, 2, 5, 13, NULL},
163: #endif
164: };

166: static GmshCellInfo GmshCellMap[150];

168: static PetscErrorCode GmshCellInfoSetUp(void)
169: {
170:   size_t           i, n;
171:   static PetscBool called = PETSC_FALSE;

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

189: #define GmshCellTypeCheck(ct) \
190:   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_); \
191:                            PetscCheck(GmshCellMap[_ct_].polytope != -1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Unsupported Gmsh element type %d", _ct_);)

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

209: /*
210:    Returns an array of count items each with a sizeof(eltsize)
211: */
212: static PetscErrorCode GmshBufferGet(GmshFile *gmsh, PetscCount count, size_t eltsize, void *buf)
213: {
214:   size_t size = count * eltsize;

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

226: /*
227:     Returns an array of count items each with the size determined by the GmshFile
228: */
229: static PetscErrorCode GmshBufferSizeGet(GmshFile *gmsh, PetscCount count, void *buf)
230: {
231:   size_t dataSize = (size_t)gmsh->dataSize;
232:   size_t size     = count * dataSize;

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

244: /*
245:     Reads an array of count items each with the size determined by the PetscDataType
246: */
247: static PetscErrorCode GmshRead(GmshFile *gmsh, void *buf, PetscCount count, PetscDataType dtype)
248: {
249:   PetscInt icount;

251:   PetscFunctionBegin;
252:   PetscCall(PetscIntCast(count, &icount));
253:   PetscCall(PetscViewerRead(gmsh->viewer, buf, icount, NULL, dtype));
254:   if (gmsh->byteSwap) PetscCall(PetscByteSwap(buf, dtype, icount));
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }

258: static PetscErrorCode GmshReadString(GmshFile *gmsh, char *buf, PetscCount count)
259: {
260:   PetscInt icount;

262:   PetscFunctionBegin;
263:   PetscCall(PetscIntCast(count, &icount));
264:   PetscCall(PetscViewerRead(gmsh->viewer, buf, icount, NULL, PETSC_STRING));
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: static PetscErrorCode GmshMatch(PETSC_UNUSED GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN], PetscBool *match)
269: {
270:   PetscFunctionBegin;
271:   PetscCall(PetscStrcmp(line, Section, match));
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: static PetscErrorCode GmshExpect(GmshFile *gmsh, const char Section[], char line[PETSC_MAX_PATH_LEN])
276: {
277:   PetscBool match;

279:   PetscFunctionBegin;
280:   PetscCall(GmshMatch(gmsh, Section, line, &match));
281:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file, expecting %s, not %s", Section, line);
282:   PetscFunctionReturn(PETSC_SUCCESS);
283: }

285: static PetscErrorCode GmshReadSection(GmshFile *gmsh, char line[PETSC_MAX_PATH_LEN])
286: {
287:   PetscBool match;

289:   PetscFunctionBegin;
290:   while (PETSC_TRUE) {
291:     PetscCall(GmshReadString(gmsh, line, 1));
292:     PetscCall(GmshMatch(gmsh, "$Comments", line, &match));
293:     if (!match) break;
294:     while (PETSC_TRUE) {
295:       PetscCall(GmshReadString(gmsh, line, 1));
296:       PetscCall(GmshMatch(gmsh, "$EndComments", line, &match));
297:       if (match) break;
298:     }
299:   }
300:   PetscFunctionReturn(PETSC_SUCCESS);
301: }

303: static PetscErrorCode GmshReadEndSection(GmshFile *gmsh, const char EndSection[], char line[PETSC_MAX_PATH_LEN])
304: {
305:   PetscFunctionBegin;
306:   PetscCall(GmshReadString(gmsh, line, 1));
307:   PetscCall(GmshExpect(gmsh, EndSection, line));
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

311: /*
312:      Read into buf[] count number of PetscInt integers (the file storage size may be different than PetscInt)
313: */
314: static PetscErrorCode GmshReadPetscInt(GmshFile *gmsh, PetscInt *buf, PetscCount count)
315: {
316:   PetscCount i;
317:   size_t     dataSize = (size_t)gmsh->dataSize;

319:   PetscFunctionBegin;
320:   if (dataSize == sizeof(PetscInt)) {
321:     PetscCall(GmshRead(gmsh, buf, count, PETSC_INT));
322:   } else if (dataSize == sizeof(int)) {
323:     int *ibuf = NULL;
324:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
325:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_ENUM));
326:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
327:   } else if (dataSize == sizeof(long)) {
328:     long *ibuf = NULL;
329:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
330:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_LONG));
331:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
332:   } else if (dataSize == sizeof(PetscInt64)) {
333:     PetscInt64 *ibuf = NULL;
334:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
335:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_INT64));
336:     for (i = 0; i < count; ++i) buf[i] = (PetscInt)ibuf[i];
337:   }
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: /*
342:      Read into buf[] count number of PetscCount integers  (the file storage size may be different than PetscCount)
343: */
344: static PetscErrorCode GmshReadPetscCount(GmshFile *gmsh, PetscCount *buf, PetscCount count)
345: {
346:   PetscCount i;
347:   size_t     dataSize = (size_t)gmsh->dataSize;

349:   PetscFunctionBegin;
350:   if (dataSize == sizeof(PetscCount)) {
351:     PetscCall(GmshRead(gmsh, buf, count, PETSC_COUNT));
352:   } else if (dataSize == sizeof(int)) {
353:     int *ibuf = NULL;
354:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
355:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_ENUM));
356:     for (i = 0; i < count; ++i) buf[i] = (PetscCount)ibuf[i];
357:   } else if (dataSize == sizeof(long)) {
358:     long *ibuf = NULL;
359:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
360:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_LONG));
361:     for (i = 0; i < count; ++i) buf[i] = (PetscCount)ibuf[i];
362:   } else if (dataSize == sizeof(PetscInt64)) {
363:     PetscInt64 *ibuf = NULL;
364:     PetscCall(GmshBufferSizeGet(gmsh, count, &ibuf));
365:     PetscCall(GmshRead(gmsh, ibuf, count, PETSC_INT64));
366:     for (i = 0; i < count; ++i) buf[i] = (PetscCount)ibuf[i];
367:   }
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }

371: /*
372:      Read into buf[] count number of PetscEnum integers
373: */
374: static PetscErrorCode GmshReadInt(GmshFile *gmsh, int *buf, PetscCount count)
375: {
376:   PetscFunctionBegin;
377:   PetscCall(GmshRead(gmsh, buf, count, PETSC_ENUM));
378:   PetscFunctionReturn(PETSC_SUCCESS);
379: }

381: /*
382:      Read into buf[] count number of double
383: */
384: static PetscErrorCode GmshReadDouble(GmshFile *gmsh, double *buf, PetscCount count)
385: {
386:   PetscFunctionBegin;
387:   PetscCall(GmshRead(gmsh, buf, count, PETSC_DOUBLE));
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: #define GMSH_MAX_TAGS 16

393: typedef struct {
394:   PetscInt id;                  /* Entity ID */
395:   PetscInt dim;                 /* Dimension */
396:   double   bbox[6];             /* Bounding box */
397:   PetscInt numTags;             /* Size of tag array */
398:   int      tags[GMSH_MAX_TAGS]; /* Tag array */
399: } GmshEntity;

401: typedef struct {
402:   GmshEntity *entity[4];
403:   PetscHMapI  entityMap[4];
404: } GmshEntities;

406: static PetscErrorCode GmshEntitiesCreate(PetscCount count[4], GmshEntities **entities)
407: {
408:   PetscFunctionBegin;
409:   PetscCall(PetscNew(entities));
410:   for (PetscInt dim = 0; dim < 4; ++dim) {
411:     PetscCall(PetscCalloc1(count[dim], &(*entities)->entity[dim]));
412:     PetscCall(PetscHMapICreate(&(*entities)->entityMap[dim]));
413:   }
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

417: static PetscErrorCode GmshEntitiesDestroy(GmshEntities **entities)
418: {
419:   PetscInt dim;

421:   PetscFunctionBegin;
422:   if (!*entities) PetscFunctionReturn(PETSC_SUCCESS);
423:   for (dim = 0; dim < 4; ++dim) {
424:     PetscCall(PetscFree((*entities)->entity[dim]));
425:     PetscCall(PetscHMapIDestroy(&(*entities)->entityMap[dim]));
426:   }
427:   PetscCall(PetscFree(*entities));
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: static PetscErrorCode GmshEntitiesAdd(GmshEntities *entities, PetscInt index, PetscInt dim, PetscInt eid, GmshEntity **entity)
432: {
433:   PetscFunctionBegin;
434:   PetscCall(PetscHMapISet(entities->entityMap[dim], eid, index));
435:   entities->entity[dim][index].dim = dim;
436:   entities->entity[dim][index].id  = eid;
437:   if (entity) *entity = &entities->entity[dim][index];
438:   PetscFunctionReturn(PETSC_SUCCESS);
439: }

441: static PetscErrorCode GmshEntitiesGet(GmshEntities *entities, PetscInt dim, PetscInt eid, GmshEntity **entity)
442: {
443:   PetscInt index;

445:   PetscFunctionBegin;
446:   PetscCall(PetscHMapIGet(entities->entityMap[dim], eid, &index));
447:   *entity = &entities->entity[dim][index];
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

451: typedef struct {
452:   PetscInt *id;  /* Node IDs */
453:   double   *xyz; /* Coordinates */
454:   PetscInt *tag; /* Physical tag */
455: } GmshNodes;

457: static PetscErrorCode GmshNodesCreate(PetscCount count, GmshNodes **nodes)
458: {
459:   PetscFunctionBegin;
460:   PetscCall(PetscNew(nodes));
461:   PetscCall(PetscMalloc1(count * 1, &(*nodes)->id));
462:   PetscCall(PetscMalloc1(count * 3, &(*nodes)->xyz));
463:   PetscCall(PetscMalloc1(count * GMSH_MAX_TAGS, &(*nodes)->tag));
464:   PetscFunctionReturn(PETSC_SUCCESS);
465: }

467: static PetscErrorCode GmshNodesDestroy(GmshNodes **nodes)
468: {
469:   PetscFunctionBegin;
470:   if (!*nodes) PetscFunctionReturn(PETSC_SUCCESS);
471:   PetscCall(PetscFree((*nodes)->id));
472:   PetscCall(PetscFree((*nodes)->xyz));
473:   PetscCall(PetscFree((*nodes)->tag));
474:   PetscCall(PetscFree(*nodes));
475:   PetscFunctionReturn(PETSC_SUCCESS);
476: }

478: typedef struct {
479:   PetscInt  id;                  /* Element ID */
480:   PetscInt  dim;                 /* Dimension */
481:   PetscInt  cellType;            /* Cell type */
482:   PetscInt  numVerts;            /* Size of vertex array */
483:   PetscInt  numNodes;            /* Size of node array */
484:   PetscInt *nodes;               /* Vertex/Node array */
485:   PetscInt  numTags;             /* Size of physical tag array */
486:   int       tags[GMSH_MAX_TAGS]; /* Physical tag array */
487: } GmshElement;

489: static PetscErrorCode GmshElementsCreate(PetscCount count, GmshElement **elements)
490: {
491:   PetscFunctionBegin;
492:   PetscCall(PetscCalloc1(count, elements));
493:   PetscFunctionReturn(PETSC_SUCCESS);
494: }

496: static PetscErrorCode GmshElementsDestroy(GmshElement **elements)
497: {
498:   PetscFunctionBegin;
499:   if (!*elements) PetscFunctionReturn(PETSC_SUCCESS);
500:   PetscCall(PetscFree(*elements));
501:   PetscFunctionReturn(PETSC_SUCCESS);
502: }

504: typedef struct {
505:   PetscInt       dim;
506:   PetscInt       order;
507:   GmshEntities  *entities;
508:   PetscInt       numNodes;
509:   GmshNodes     *nodelist;
510:   PetscInt       numElems;
511:   GmshElement   *elements;
512:   PetscInt       numVerts;
513:   PetscInt       numCells;
514:   PetscInt      *periodMap;
515:   PetscInt      *vertexMap;
516:   PetscSegBuffer segbuf;
517:   PetscInt       numRegions;
518:   PetscInt      *regionDims;
519:   PetscInt      *regionTags;
520:   char         **regionNames;
521: } GmshMesh;

523: static PetscErrorCode GmshMeshCreate(GmshMesh **mesh)
524: {
525:   PetscFunctionBegin;
526:   PetscCall(PetscNew(mesh));
527:   PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 0, &(*mesh)->segbuf));
528:   PetscFunctionReturn(PETSC_SUCCESS);
529: }

531: static PetscErrorCode GmshMeshDestroy(GmshMesh **mesh)
532: {
533:   PetscInt r;

535:   PetscFunctionBegin;
536:   if (!*mesh) PetscFunctionReturn(PETSC_SUCCESS);
537:   PetscCall(GmshEntitiesDestroy(&(*mesh)->entities));
538:   PetscCall(GmshNodesDestroy(&(*mesh)->nodelist));
539:   PetscCall(GmshElementsDestroy(&(*mesh)->elements));
540:   PetscCall(PetscFree((*mesh)->periodMap));
541:   PetscCall(PetscFree((*mesh)->vertexMap));
542:   PetscCall(PetscSegBufferDestroy(&(*mesh)->segbuf));
543:   for (r = 0; r < (*mesh)->numRegions; ++r) PetscCall(PetscFree((*mesh)->regionNames[r]));
544:   PetscCall(PetscFree3((*mesh)->regionDims, (*mesh)->regionTags, (*mesh)->regionNames));
545:   PetscCall(PetscFree(*mesh));
546:   PetscFunctionReturn(PETSC_SUCCESS);
547: }

549: static PetscErrorCode GmshReadNodes_v22(GmshFile *gmsh, GmshMesh *mesh)
550: {
551:   PetscViewer viewer   = gmsh->viewer;
552:   PetscBool   byteSwap = gmsh->byteSwap;
553:   char        line[PETSC_MAX_PATH_LEN];
554:   int         n, t, num, nid, snum;
555:   GmshNodes  *nodes;

557:   PetscFunctionBegin;
558:   PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
559:   snum = sscanf(line, "%d", &num);
560:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
561:   PetscCall(GmshNodesCreate(num, &nodes));
562:   mesh->numNodes = num;
563:   mesh->nodelist = nodes;
564:   for (n = 0; n < num; ++n) {
565:     double *xyz = nodes->xyz + n * 3;
566:     PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
567:     PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
568:     if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
569:     if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
570:     nodes->id[n] = nid;
571:     for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
572:   }
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: /* Gmsh elements can be of any dimension/co-dimension, so we need to traverse the
577:    file contents multiple times to figure out the true number of cells and facets
578:    in the given mesh. To make this more efficient we read the file contents only
579:    once and store them in memory, while determining the true number of cells. */
580: static PetscErrorCode GmshReadElements_v22(GmshFile *gmsh, GmshMesh *mesh)
581: {
582:   PetscViewer  viewer   = gmsh->viewer;
583:   PetscBool    binary   = gmsh->binary;
584:   PetscBool    byteSwap = gmsh->byteSwap;
585:   char         line[PETSC_MAX_PATH_LEN];
586:   int          i, c, p, num, ibuf[1 + 4 + 1000], snum;
587:   int          cellType, numElem, numVerts, numNodes, numTags;
588:   GmshElement *elements;
589:   PetscInt    *nodeMap = gmsh->nodeMap;

591:   PetscFunctionBegin;
592:   PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
593:   snum = sscanf(line, "%d", &num);
594:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
595:   PetscCall(GmshElementsCreate(num, &elements));
596:   mesh->numElems = num;
597:   mesh->elements = elements;
598:   for (c = 0; c < num;) {
599:     PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
600:     if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));

602:     cellType = binary ? ibuf[0] : ibuf[1];
603:     numElem  = binary ? ibuf[1] : 1;
604:     numTags  = ibuf[2];

606:     PetscCall(GmshCellTypeCheck(cellType));
607:     numVerts = GmshCellMap[cellType].numVerts;
608:     numNodes = GmshCellMap[cellType].numNodes;

610:     for (i = 0; i < numElem; ++i, ++c) {
611:       GmshElement *element = elements + c;
612:       const int    off = binary ? 0 : 1, nint = 1 + numTags + numNodes - off;
613:       const int   *id = ibuf, *nodes = ibuf + 1 + numTags, *tags = ibuf + 1;
614:       PetscCall(PetscViewerRead(viewer, ibuf + off, nint, NULL, PETSC_ENUM));
615:       if (byteSwap) PetscCall(PetscByteSwap(ibuf + off, PETSC_ENUM, nint));
616:       element->id       = id[0];
617:       element->dim      = GmshCellMap[cellType].dim;
618:       element->cellType = cellType;
619:       element->numVerts = numVerts;
620:       element->numNodes = numNodes;
621:       element->numTags  = PetscMin(numTags, GMSH_MAX_TAGS);
622:       PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
623:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
624:       for (p = 0; p < element->numTags; p++) element->tags[p] = tags[p];
625:     }
626:   }
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: /*
631: $Entities
632:   numPoints(unsigned long) numCurves(unsigned long)
633:     numSurfaces(unsigned long) numVolumes(unsigned long)
634:   // points
635:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
636:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
637:     numPhysicals(unsigned long) phyisicalTag[...](int)
638:   ...
639:   // curves
640:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
641:      boxMaxX(double) boxMaxY(double) boxMaxZ(double)
642:      numPhysicals(unsigned long) physicalTag[...](int)
643:      numBREPVert(unsigned long) tagBREPVert[...](int)
644:   ...
645:   // surfaces
646:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
647:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
648:     numPhysicals(unsigned long) physicalTag[...](int)
649:     numBREPCurve(unsigned long) tagBREPCurve[...](int)
650:   ...
651:   // volumes
652:   tag(int) boxMinX(double) boxMinY(double) boxMinZ(double)
653:     boxMaxX(double) boxMaxY(double) boxMaxZ(double)
654:     numPhysicals(unsigned long) physicalTag[...](int)
655:     numBREPSurfaces(unsigned long) tagBREPSurfaces[...](int)
656:   ...
657: $EndEntities
658: */
659: static PetscErrorCode GmshReadEntities_v40(GmshFile *gmsh, GmshMesh *mesh)
660: {
661:   PetscViewer viewer   = gmsh->viewer;
662:   PetscBool   byteSwap = gmsh->byteSwap;
663:   long        lnum, lbuf[4];
664:   int         dim, eid, numTags, *ibuf, t;
665:   PetscCount  index, count[4];
666:   PetscInt    i, num;
667:   GmshEntity *entity = NULL;

669:   PetscFunctionBegin;
670:   PetscCall(PetscViewerRead(viewer, lbuf, 4, NULL, PETSC_LONG));
671:   if (byteSwap) PetscCall(PetscByteSwap(lbuf, PETSC_LONG, 4));
672:   for (i = 0; i < 4; ++i) count[i] = (PetscCount)lbuf[i];
673:   PetscCall(GmshEntitiesCreate(count, &mesh->entities));
674:   for (dim = 0; dim < 4; ++dim) {
675:     for (index = 0; index < count[dim]; ++index) {
676:       PetscCall(PetscViewerRead(viewer, &eid, 1, NULL, PETSC_ENUM));
677:       if (byteSwap) PetscCall(PetscByteSwap(&eid, PETSC_ENUM, 1));
678:       PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
679:       PetscCall(PetscViewerRead(viewer, entity->bbox, 6, NULL, PETSC_DOUBLE));
680:       if (byteSwap) PetscCall(PetscByteSwap(entity->bbox, PETSC_DOUBLE, 6));
681:       PetscCall(PetscViewerRead(viewer, &lnum, 1, NULL, PETSC_LONG));
682:       if (byteSwap) PetscCall(PetscByteSwap(&lnum, PETSC_LONG, 1));
683:       PetscCall(PetscIntCast(lnum, &num));
684:       PetscCall(GmshBufferGet(gmsh, num, sizeof(int), &ibuf));
685:       PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
686:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
687:       entity->numTags = numTags = (int)PetscMin(num, GMSH_MAX_TAGS);
688:       for (t = 0; t < numTags; ++t) entity->tags[t] = ibuf[t];
689:       if (dim == 0) continue;
690:       PetscCall(PetscViewerRead(viewer, &lnum, 1, NULL, PETSC_LONG));
691:       if (byteSwap) PetscCall(PetscByteSwap(&lnum, PETSC_LONG, 1));
692:       PetscCall(GmshBufferGet(gmsh, lnum, sizeof(int), &ibuf));
693:       PetscCall(PetscIntCast(lnum, &num));
694:       PetscCall(PetscViewerRead(viewer, ibuf, num, NULL, PETSC_ENUM));
695:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, num));
696:     }
697:   }
698:   PetscFunctionReturn(PETSC_SUCCESS);
699: }

701: /*
702: $Nodes
703:   numEntityBlocks(unsigned long) numNodes(unsigned long)
704:   tagEntity(int) dimEntity(int) typeNode(int) numNodes(unsigned long)
705:     tag(int) x(double) y(double) z(double)
706:     ...
707:   ...
708: $EndNodes
709: */
710: static PetscErrorCode GmshReadNodes_v40(GmshFile *gmsh, GmshMesh *mesh)
711: {
712:   PetscViewer viewer   = gmsh->viewer;
713:   PetscBool   byteSwap = gmsh->byteSwap;
714:   long        block, node, n, t, numEntityBlocks, numTotalNodes, numNodes;
715:   int         info[3], nid;
716:   GmshNodes  *nodes;

718:   PetscFunctionBegin;
719:   PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
720:   if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
721:   PetscCall(PetscViewerRead(viewer, &numTotalNodes, 1, NULL, PETSC_LONG));
722:   if (byteSwap) PetscCall(PetscByteSwap(&numTotalNodes, PETSC_LONG, 1));
723:   PetscCall(GmshNodesCreate(numTotalNodes, &nodes));
724:   PetscCall(PetscIntCast(numTotalNodes, &mesh->numNodes));
725:   mesh->nodelist = nodes;
726:   for (n = 0, block = 0; block < numEntityBlocks; ++block) {
727:     PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
728:     PetscCall(PetscViewerRead(viewer, &numNodes, 1, NULL, PETSC_LONG));
729:     if (byteSwap) PetscCall(PetscByteSwap(&numNodes, PETSC_LONG, 1));
730:     if (gmsh->binary) {
731:       size_t   nbytes = sizeof(int) + 3 * sizeof(double);
732:       char    *cbuf   = NULL; /* dummy value to prevent warning from compiler about possible uninitialized value */
733:       PetscInt icnt;

735:       PetscCall(GmshBufferGet(gmsh, numNodes, nbytes, &cbuf));
736:       PetscCall(PetscIntCast(numNodes * nbytes, &icnt));
737:       PetscCall(PetscViewerRead(viewer, cbuf, icnt, NULL, PETSC_CHAR));
738:       for (node = 0; node < numNodes; ++node, ++n) {
739:         char   *cnid = cbuf + node * nbytes, *cxyz = cnid + sizeof(int);
740:         double *xyz = nodes->xyz + n * 3;

742:         if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cnid, PETSC_ENUM, 1));
743:         if (!PetscBinaryBigEndian()) PetscCall(PetscByteSwap(cxyz, PETSC_DOUBLE, 3));
744:         PetscCall(PetscMemcpy(&nid, cnid, sizeof(int)));
745:         PetscCall(PetscMemcpy(xyz, cxyz, 3 * sizeof(double)));
746:         if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
747:         if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
748:         nodes->id[n] = nid;
749:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
750:       }
751:     } else {
752:       for (node = 0; node < numNodes; ++node, ++n) {
753:         double *xyz = nodes->xyz + n * 3;

755:         PetscCall(PetscViewerRead(viewer, &nid, 1, NULL, PETSC_ENUM));
756:         PetscCall(PetscViewerRead(viewer, xyz, 3, NULL, PETSC_DOUBLE));
757:         if (byteSwap) PetscCall(PetscByteSwap(&nid, PETSC_ENUM, 1));
758:         if (byteSwap) PetscCall(PetscByteSwap(xyz, PETSC_DOUBLE, 3));
759:         nodes->id[n] = nid;
760:         for (t = 0; t < GMSH_MAX_TAGS; ++t) nodes->tag[n * GMSH_MAX_TAGS + t] = -1;
761:       }
762:     }
763:   }
764:   PetscFunctionReturn(PETSC_SUCCESS);
765: }

767: /*
768: $Elements
769:   numEntityBlocks(unsigned long) numElements(unsigned long)
770:   tagEntity(int) dimEntity(int) typeEle(int) numElements(unsigned long)
771:     tag(int) numVert[...](int)
772:     ...
773:   ...
774: $EndElements
775: */
776: static PetscErrorCode GmshReadElements_v40(GmshFile *gmsh, GmshMesh *mesh)
777: {
778:   PetscViewer  viewer   = gmsh->viewer;
779:   PetscBool    byteSwap = gmsh->byteSwap;
780:   long         c, block, numEntityBlocks, numTotalElements, elem, numElements;
781:   int          p, info[3], *ibuf = NULL;
782:   int          eid, dim, cellType, numVerts, numNodes, numTags;
783:   GmshEntity  *entity = NULL;
784:   GmshElement *elements;
785:   PetscInt    *nodeMap = gmsh->nodeMap, icnt;

787:   PetscFunctionBegin;
788:   PetscCall(PetscViewerRead(viewer, &numEntityBlocks, 1, NULL, PETSC_LONG));
789:   if (byteSwap) PetscCall(PetscByteSwap(&numEntityBlocks, PETSC_LONG, 1));
790:   PetscCall(PetscViewerRead(viewer, &numTotalElements, 1, NULL, PETSC_LONG));
791:   if (byteSwap) PetscCall(PetscByteSwap(&numTotalElements, PETSC_LONG, 1));
792:   PetscCall(GmshElementsCreate(numTotalElements, &elements));
793:   PetscCall(PetscIntCast(numTotalElements, &mesh->numElems));
794:   mesh->elements = elements;
795:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
796:     PetscCall(PetscViewerRead(viewer, info, 3, NULL, PETSC_ENUM));
797:     if (byteSwap) PetscCall(PetscByteSwap(info, PETSC_ENUM, 3));
798:     eid      = info[0];
799:     dim      = info[1];
800:     cellType = info[2];
801:     PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
802:     PetscCall(GmshCellTypeCheck(cellType));
803:     numVerts = GmshCellMap[cellType].numVerts;
804:     numNodes = GmshCellMap[cellType].numNodes;
805:     numTags  = (int)entity->numTags;
806:     PetscCall(PetscViewerRead(viewer, &numElements, 1, NULL, PETSC_LONG));
807:     if (byteSwap) PetscCall(PetscByteSwap(&numElements, PETSC_LONG, 1));
808:     PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numElements, sizeof(int), &ibuf));
809:     PetscCall(PetscIntCast((1 + numNodes) * numElements, &icnt));
810:     PetscCall(PetscViewerRead(viewer, ibuf, icnt, NULL, PETSC_ENUM));
811:     if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, (1 + numNodes) * numElements));
812:     for (elem = 0; elem < numElements; ++elem, ++c) {
813:       GmshElement *element = elements + c;
814:       const int   *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;
815:       element->id       = id[0];
816:       element->dim      = dim;
817:       element->cellType = cellType;
818:       element->numVerts = numVerts;
819:       element->numNodes = numNodes;
820:       element->numTags  = numTags;
821:       PetscCall(PetscSegBufferGet(mesh->segbuf, (size_t)element->numNodes, &element->nodes));
822:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
823:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
824:     }
825:   }
826:   PetscFunctionReturn(PETSC_SUCCESS);
827: }

829: static PetscErrorCode GmshReadPeriodic_v40(GmshFile *gmsh, PetscInt periodicMap[])
830: {
831:   PetscViewer viewer     = gmsh->viewer;
832:   int         fileFormat = gmsh->fileFormat;
833:   PetscBool   binary     = gmsh->binary;
834:   PetscBool   byteSwap   = gmsh->byteSwap;
835:   int         numPeriodic, snum, i;
836:   char        line[PETSC_MAX_PATH_LEN];
837:   PetscInt   *nodeMap = gmsh->nodeMap;

839:   PetscFunctionBegin;
840:   if (fileFormat == 22 || !binary) {
841:     PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
842:     snum = sscanf(line, "%d", &numPeriodic);
843:     PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
844:   } else {
845:     PetscCall(PetscViewerRead(viewer, &numPeriodic, 1, NULL, PETSC_ENUM));
846:     if (byteSwap) PetscCall(PetscByteSwap(&numPeriodic, PETSC_ENUM, 1));
847:   }
848:   for (i = 0; i < numPeriodic; i++) {
849:     int    ibuf[3], correspondingDim = -1, correspondingTag = -1, primaryTag = -1, correspondingNode, primaryNode;
850:     long   j, nNodes;
851:     double affine[16];

853:     if (fileFormat == 22 || !binary) {
854:       PetscCall(PetscViewerRead(viewer, line, 3, NULL, PETSC_STRING));
855:       snum = sscanf(line, "%d %d %d", &correspondingDim, &correspondingTag, &primaryTag);
856:       PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
857:     } else {
858:       PetscCall(PetscViewerRead(viewer, ibuf, 3, NULL, PETSC_ENUM));
859:       if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 3));
860:       correspondingDim = ibuf[0];
861:       correspondingTag = ibuf[1];
862:       primaryTag       = ibuf[2];
863:     }
864:     (void)correspondingDim;
865:     (void)correspondingTag;
866:     (void)primaryTag; /* unused */

868:     if (fileFormat == 22 || !binary) {
869:       PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
870:       snum = sscanf(line, "%ld", &nNodes);
871:       if (snum != 1) { /* discard transformation and try again */
872:         PetscCall(PetscViewerRead(viewer, line, -PETSC_MAX_PATH_LEN, NULL, PETSC_STRING));
873:         PetscCall(PetscViewerRead(viewer, line, 1, NULL, PETSC_STRING));
874:         snum = sscanf(line, "%ld", &nNodes);
875:         PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
876:       }
877:     } else {
878:       PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
879:       if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
880:       if (nNodes == -1) { /* discard transformation and try again */
881:         PetscCall(PetscViewerRead(viewer, affine, 16, NULL, PETSC_DOUBLE));
882:         PetscCall(PetscViewerRead(viewer, &nNodes, 1, NULL, PETSC_LONG));
883:         if (byteSwap) PetscCall(PetscByteSwap(&nNodes, PETSC_LONG, 1));
884:       }
885:     }

887:     for (j = 0; j < nNodes; j++) {
888:       if (fileFormat == 22 || !binary) {
889:         PetscCall(PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING));
890:         snum = sscanf(line, "%d %d", &correspondingNode, &primaryNode);
891:         PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
892:       } else {
893:         PetscCall(PetscViewerRead(viewer, ibuf, 2, NULL, PETSC_ENUM));
894:         if (byteSwap) PetscCall(PetscByteSwap(ibuf, PETSC_ENUM, 2));
895:         correspondingNode = ibuf[0];
896:         primaryNode       = ibuf[1];
897:       }
898:       correspondingNode              = (int)nodeMap[correspondingNode];
899:       primaryNode                    = (int)nodeMap[primaryNode];
900:       periodicMap[correspondingNode] = primaryNode;
901:     }
902:   }
903:   PetscFunctionReturn(PETSC_SUCCESS);
904: }

906: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
907: $Entities
908:   numPoints(size_t) numCurves(size_t)
909:     numSurfaces(size_t) numVolumes(size_t)
910:   pointTag(int) X(double) Y(double) Z(double)
911:     numPhysicalTags(size_t) physicalTag(int) ...
912:   ...
913:   curveTag(int) minX(double) minY(double) minZ(double)
914:     maxX(double) maxY(double) maxZ(double)
915:     numPhysicalTags(size_t) physicalTag(int) ...
916:     numBoundingPoints(size_t) pointTag(int) ...
917:   ...
918:   surfaceTag(int) minX(double) minY(double) minZ(double)
919:     maxX(double) maxY(double) maxZ(double)
920:     numPhysicalTags(size_t) physicalTag(int) ...
921:     numBoundingCurves(size_t) curveTag(int) ...
922:   ...
923:   volumeTag(int) minX(double) minY(double) minZ(double)
924:     maxX(double) maxY(double) maxZ(double)
925:     numPhysicalTags(size_t) physicalTag(int) ...
926:     numBoundngSurfaces(size_t) surfaceTag(int) ...
927:   ...
928: $EndEntities
929: */
930: static PetscErrorCode GmshReadEntities_v41(GmshFile *gmsh, GmshMesh *mesh)
931: {
932:   PetscCount  count[4], index, numTags;
933:   int         dim, eid, *tags = NULL;
934:   GmshEntity *entity = NULL;

936:   PetscFunctionBegin;
937:   PetscCall(GmshReadPetscCount(gmsh, count, 4));
938:   PetscCall(GmshEntitiesCreate(count, &mesh->entities));
939:   for (dim = 0; dim < 4; ++dim) {
940:     for (index = 0; index < count[dim]; ++index) {
941:       PetscCall(GmshReadInt(gmsh, &eid, 1));
942:       PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
943:       PetscCall(GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6));
944:       PetscCall(GmshReadPetscCount(gmsh, &numTags, 1));
945:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
946:       PetscCall(GmshReadInt(gmsh, tags, numTags));
947:       PetscCheck(numTags <= GMSH_MAX_TAGS, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "PETSc currently supports up to %" PetscInt_FMT " tags per entity, not %" PetscCount_FMT, (PetscInt)GMSH_MAX_TAGS, numTags);
948:       PetscCall(PetscIntCast(numTags, &entity->numTags));
949:       for (PetscInt i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
950:       if (dim == 0) continue;
951:       PetscCall(GmshReadPetscCount(gmsh, &numTags, 1));
952:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
953:       PetscCall(GmshReadInt(gmsh, tags, numTags));
954:       /* Currently, we do not save the ids for the bounding entities */
955:     }
956:   }
957:   PetscFunctionReturn(PETSC_SUCCESS);
958: }

960: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
961: $Nodes
962:   numEntityBlocks(size_t) numNodes(size_t)
963:     minNodeTag(size_t) maxNodeTag(size_t)
964:   entityDim(int) entityTag(int) parametric(int; 0 or 1) numNodesBlock(size_t)
965:     nodeTag(size_t)
966:     ...
967:     x(double) y(double) z(double)
968:        < u(double; if parametric and entityDim = 1 or entityDim = 2) >
969:        < v(double; if parametric and entityDim = 2) >
970:     ...
971:   ...
972: $EndNodes
973: */
974: static PetscErrorCode GmshReadNodes_v41(GmshFile *gmsh, GmshMesh *mesh)
975: {
976:   int         info[3], dim, eid, parametric;
977:   PetscCount  sizes[4], numEntityBlocks, numNodes, numNodesBlock = 0;
978:   PetscInt    numTags;
979:   GmshEntity *entity = NULL;
980:   GmshNodes  *nodes;

982:   PetscFunctionBegin;
983:   PetscCall(GmshReadPetscCount(gmsh, sizes, 4));
984:   numEntityBlocks = sizes[0];
985:   numNodes        = sizes[1];
986:   PetscCall(GmshNodesCreate(numNodes, &nodes));
987:   PetscCall(PetscIntCast(numNodes, &mesh->numNodes));
988:   mesh->nodelist = nodes;
989:   if (numEntityBlocks && !mesh->entities) PetscCall(PetscInfo(NULL, "File specifies %" PetscCount_FMT " entity blocks, but was missing the $Entities section\n", numEntityBlocks));
990:   for (PetscCount block = 0, node = 0; block < numEntityBlocks; ++block, node += numNodesBlock) {
991:     PetscCall(GmshReadInt(gmsh, info, 3));
992:     dim        = info[0];
993:     eid        = info[1];
994:     parametric = info[2];
995:     if (mesh->entities) PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
996:     numTags = entity ? entity->numTags : 0;
997:     PetscCheck(!parametric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Parametric coordinates not supported");
998:     PetscCall(GmshReadPetscCount(gmsh, &numNodesBlock, 1));
999:     PetscCall(GmshReadPetscInt(gmsh, nodes->id + node, numNodesBlock));
1000:     PetscCall(GmshReadDouble(gmsh, nodes->xyz + node * 3, numNodesBlock * 3));
1001:     for (PetscCount n = 0; n < numNodesBlock; ++n) {
1002:       PetscInt *tags = &nodes->tag[node * GMSH_MAX_TAGS];

1004:       for (PetscInt t = 0; t < numTags; ++t) tags[n * GMSH_MAX_TAGS + t] = entity->tags[t];
1005:       for (PetscInt t = numTags; t < GMSH_MAX_TAGS; ++t) tags[n * GMSH_MAX_TAGS + t] = -1;
1006:     }
1007:   }
1008:   PetscCall(PetscIntCast(sizes[2], &gmsh->nodeStart));
1009:   PetscCall(PetscIntCast(sizes[3] + 1, &gmsh->nodeEnd));
1010:   PetscFunctionReturn(PETSC_SUCCESS);
1011: }

1013: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1014: $Elements
1015:   numEntityBlocks(size_t) numElements(size_t)
1016:     minElementTag(size_t) maxElementTag(size_t)
1017:   entityDim(int) entityTag(int) elementType(int; see below) numElementsBlock(size_t)
1018:     elementTag(size_t) nodeTag(size_t) ...
1019:     ...
1020:   ...
1021: $EndElements
1022: */
1023: static PetscErrorCode GmshReadElements_v41(GmshFile *gmsh, GmshMesh *mesh)
1024: {
1025:   int          info[3], eid, dim, cellType;
1026:   PetscCount   sizes[4], numEntityBlocks, numElements, numBlockElements, numVerts, numNodes, numTags, block, elem, c, p;
1027:   GmshEntity  *entity = NULL;
1028:   GmshElement *elements;
1029:   PetscInt    *nodeMap = gmsh->nodeMap, *ibuf = NULL;

1031:   PetscFunctionBegin;
1032:   PetscCall(GmshReadPetscCount(gmsh, sizes, 4));
1033:   numEntityBlocks = sizes[0];
1034:   numElements     = sizes[1];
1035:   PetscCall(GmshElementsCreate(numElements, &elements));
1036:   PetscCall(PetscIntCast(numElements, &mesh->numElems));
1037:   mesh->elements = elements;
1038:   if (numEntityBlocks && !mesh->entities) PetscCall(PetscInfo(NULL, "File specifies %" PetscCount_FMT " entity blocks, but was missing the $Entities section\n", numEntityBlocks));
1039:   for (c = 0, block = 0; block < numEntityBlocks; ++block) {
1040:     PetscCall(GmshReadInt(gmsh, info, 3));
1041:     dim      = info[0];
1042:     eid      = info[1];
1043:     cellType = info[2];
1044:     if (mesh->entities) PetscCall(GmshEntitiesGet(mesh->entities, dim, eid, &entity));
1045:     PetscCall(GmshCellTypeCheck(cellType));
1046:     numVerts = GmshCellMap[cellType].numVerts;
1047:     numNodes = GmshCellMap[cellType].numNodes;
1048:     numTags  = entity ? entity->numTags : 0;
1049:     PetscCall(GmshReadPetscCount(gmsh, &numBlockElements, 1));
1050:     PetscCall(GmshBufferGet(gmsh, (1 + numNodes) * numBlockElements, sizeof(PetscInt), &ibuf));
1051:     PetscCall(GmshReadPetscInt(gmsh, ibuf, (1 + numNodes) * numBlockElements));
1052:     for (elem = 0; elem < numBlockElements; ++elem, ++c) {
1053:       GmshElement    *element = elements + c;
1054:       const PetscInt *id = ibuf + elem * (1 + numNodes), *nodes = id + 1;

1056:       element->id       = id[0];
1057:       element->dim      = dim;
1058:       element->cellType = cellType;
1059:       PetscCall(PetscIntCast(numVerts, &element->numVerts));
1060:       PetscCall(PetscIntCast(numNodes, &element->numNodes));
1061:       PetscCall(PetscIntCast(numTags, &element->numTags));
1062:       PetscCall(PetscSegBufferGet(mesh->segbuf, element->numNodes, &element->nodes));
1063:       for (p = 0; p < element->numNodes; p++) element->nodes[p] = nodeMap[nodes[p]];
1064:       for (p = 0; p < element->numTags; p++) element->tags[p] = entity->tags[p];
1065:     }
1066:   }
1067:   PetscFunctionReturn(PETSC_SUCCESS);
1068: }

1070: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1071: $Periodic
1072:   numPeriodicLinks(size_t)
1073:   entityDim(int) entityTag(int) entityTagPrimary(int)
1074:   numAffine(size_t) value(double) ...
1075:   numCorrespondingNodes(size_t)
1076:     nodeTag(size_t) nodeTagPrimary(size_t)
1077:     ...
1078:   ...
1079: $EndPeriodic
1080: */
1081: static PetscErrorCode GmshReadPeriodic_v41(GmshFile *gmsh, PetscInt periodicMap[])
1082: {
1083:   int        info[3];
1084:   double     dbuf[16];
1085:   PetscCount numPeriodicLinks, numAffine, numCorrespondingNodes;
1086:   PetscInt  *nodeMap = gmsh->nodeMap, *nodeTags = NULL;

1088:   PetscFunctionBegin;
1089:   PetscCall(GmshReadPetscCount(gmsh, &numPeriodicLinks, 1));
1090:   for (PetscCount link = 0; link < numPeriodicLinks; ++link) {
1091:     PetscCall(GmshReadInt(gmsh, info, 3));
1092:     PetscCall(GmshReadPetscCount(gmsh, &numAffine, 1));
1093:     PetscCall(GmshReadDouble(gmsh, dbuf, numAffine));
1094:     PetscCall(GmshReadPetscCount(gmsh, &numCorrespondingNodes, 1));
1095:     PetscCall(GmshBufferGet(gmsh, numCorrespondingNodes, sizeof(PetscInt), &nodeTags));
1096:     PetscCall(GmshReadPetscInt(gmsh, nodeTags, numCorrespondingNodes * 2));
1097:     for (PetscCount node = 0; node < numCorrespondingNodes; ++node) {
1098:       PetscInt correspondingNode = nodeMap[nodeTags[node * 2 + 0]];
1099:       PetscInt primaryNode       = nodeMap[nodeTags[node * 2 + 1]];

1101:       periodicMap[correspondingNode] = primaryNode;
1102:     }
1103:   }
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1108: $MeshFormat // same as MSH version 2
1109:   version(ASCII double; currently 4.1)
1110:   file-type(ASCII int; 0 for ASCII mode, 1 for binary mode)
1111:   data-size(ASCII int; sizeof(size_t))
1112:   < int with value one; only in binary mode, to detect endianness >
1113: $EndMeshFormat
1114: */
1115: static PetscErrorCode GmshReadMeshFormat(GmshFile *gmsh)
1116: {
1117:   char  line[PETSC_MAX_PATH_LEN];
1118:   int   snum, fileType, fileFormat, dataSize, checkEndian;
1119:   float version;

1121:   PetscFunctionBegin;
1122:   PetscCall(GmshReadString(gmsh, line, 3));
1123:   snum       = sscanf(line, "%f %d %d", &version, &fileType, &dataSize);
1124:   fileFormat = (int)roundf(version * 10);
1125:   PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1126:   PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1127:   PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1128:   PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1129:   PetscCheck(!gmsh->binary || fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is binary but Gmsh file is ASCII");
1130:   PetscCheck(gmsh->binary || !fileType, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Viewer is ASCII but Gmsh file is binary");
1131:   PetscCheck(fileFormat > 40 || dataSize == sizeof(double), PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Data size %d is not valid for a Gmsh file", dataSize);
1132:   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);
1133:   gmsh->fileFormat = fileFormat;
1134:   gmsh->dataSize   = dataSize;
1135:   gmsh->byteSwap   = PETSC_FALSE;
1136:   if (gmsh->binary) {
1137:     PetscCall(GmshReadInt(gmsh, &checkEndian, 1));
1138:     if (checkEndian != 1) {
1139:       PetscCall(PetscByteSwap(&checkEndian, PETSC_ENUM, 1));
1140:       PetscCheck(checkEndian == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to detect endianness in Gmsh file header: %s", line);
1141:       gmsh->byteSwap = PETSC_TRUE;
1142:     }
1143:   }
1144:   PetscFunctionReturn(PETSC_SUCCESS);
1145: }

1147: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1148: Neper: https://neper.info/ adds this section

1150: $MeshVersion
1151:   <major>.<minor>,<patch>
1152: $EndMeshVersion
1153: */
1154: static PetscErrorCode GmshReadMeshVersion(GmshFile *gmsh)
1155: {
1156:   char line[PETSC_MAX_PATH_LEN];
1157:   int  snum, major, minor, patch;

1159:   PetscFunctionBegin;
1160:   PetscCall(GmshReadString(gmsh, line, 1));
1161:   snum = sscanf(line, "%d.%d.%d", &major, &minor, &patch);
1162:   PetscCheck(snum == 3, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1163:   PetscFunctionReturn(PETSC_SUCCESS);
1164: }

1166: /* http://gmsh.info/dev/doc/texinfo/gmsh.html#MSH-file-format
1167: Neper: https://neper.info/ adds this section

1169: $Domain
1170:   <shape>
1171: $EndDomain
1172: */
1173: static PetscErrorCode GmshReadMeshDomain(GmshFile *gmsh)
1174: {
1175:   char line[PETSC_MAX_PATH_LEN];

1177:   PetscFunctionBegin;
1178:   PetscCall(GmshReadString(gmsh, line, 1));
1179:   PetscFunctionReturn(PETSC_SUCCESS);
1180: }

1182: /*
1183: PhysicalNames
1184:   numPhysicalNames(ASCII int)
1185:   dimension(ASCII int) physicalTag(ASCII int) "name"(127 characters max)
1186:   ...
1187: $EndPhysicalNames
1188: */
1189: static PetscErrorCode GmshReadPhysicalNames(GmshFile *gmsh, GmshMesh *mesh)
1190: {
1191:   char line[PETSC_MAX_PATH_LEN], name[128 + 2], *p = NULL, *q = NULL, *r = NULL;
1192:   int  snum, region, dim, tag;

1194:   PetscFunctionBegin;
1195:   PetscCall(GmshReadString(gmsh, line, 1));
1196:   snum             = sscanf(line, "%d", &region);
1197:   mesh->numRegions = region;
1198:   PetscCheck(snum == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1199:   PetscCall(PetscMalloc3(mesh->numRegions, &mesh->regionDims, mesh->numRegions, &mesh->regionTags, mesh->numRegions, &mesh->regionNames));
1200:   for (region = 0; region < mesh->numRegions; ++region) {
1201:     PetscCall(GmshReadString(gmsh, line, 2));
1202:     snum = sscanf(line, "%d %d", &dim, &tag);
1203:     PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1204:     PetscCall(GmshReadString(gmsh, line, -(PetscInt)sizeof(line)));
1205:     PetscCall(PetscStrchr(line, '"', &p));
1206:     PetscCheck(p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1207:     PetscCall(PetscStrrchr(line, '"', &q));
1208:     PetscCheck(q != p, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "File is not a valid Gmsh file");
1209:     PetscCall(PetscStrrchr(line, ':', &r));
1210:     if (p != r) q = r;
1211:     PetscCall(PetscStrncpy(name, p + 1, (size_t)(q - p - 1)));
1212:     mesh->regionDims[region] = dim;
1213:     mesh->regionTags[region] = tag;
1214:     PetscCall(PetscStrallocpy(name, &mesh->regionNames[region]));
1215:   }
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

1219: static PetscErrorCode GmshReadEntities(GmshFile *gmsh, GmshMesh *mesh)
1220: {
1221:   PetscFunctionBegin;
1222:   switch (gmsh->fileFormat) {
1223:   case 41:
1224:     PetscCall(GmshReadEntities_v41(gmsh, mesh));
1225:     break;
1226:   default:
1227:     PetscCall(GmshReadEntities_v40(gmsh, mesh));
1228:     break;
1229:   }
1230:   PetscFunctionReturn(PETSC_SUCCESS);
1231: }

1233: static PetscErrorCode GmshReadNodes(GmshFile *gmsh, GmshMesh *mesh)
1234: {
1235:   PetscFunctionBegin;
1236:   switch (gmsh->fileFormat) {
1237:   case 41:
1238:     PetscCall(GmshReadNodes_v41(gmsh, mesh));
1239:     break;
1240:   case 40:
1241:     PetscCall(GmshReadNodes_v40(gmsh, mesh));
1242:     break;
1243:   default:
1244:     PetscCall(GmshReadNodes_v22(gmsh, mesh));
1245:     break;
1246:   }

1248:   { /* Gmsh v2.2/v4.0 does not provide min/max node tags */
1249:     if (mesh->numNodes > 0 && gmsh->nodeEnd >= gmsh->nodeStart) {
1250:       PetscInt   tagMin = PETSC_INT_MAX, tagMax = PETSC_INT_MIN, n;
1251:       GmshNodes *nodes = mesh->nodelist;
1252:       for (n = 0; n < mesh->numNodes; ++n) {
1253:         const PetscInt tag = nodes->id[n];
1254:         tagMin             = PetscMin(tag, tagMin);
1255:         tagMax             = PetscMax(tag, tagMax);
1256:       }
1257:       gmsh->nodeStart = tagMin;
1258:       gmsh->nodeEnd   = tagMax + 1;
1259:     }
1260:   }

1262:   { /* Support for sparse node tags */
1263:     PetscInt   n, t;
1264:     GmshNodes *nodes = mesh->nodelist;
1265:     PetscCall(PetscMalloc1(gmsh->nodeEnd - gmsh->nodeStart, &gmsh->nbuf));
1266:     for (t = 0; t < gmsh->nodeEnd - gmsh->nodeStart; ++t) gmsh->nbuf[t] = PETSC_INT_MIN;
1267:     gmsh->nodeMap = gmsh->nbuf - gmsh->nodeStart;
1268:     for (n = 0; n < mesh->numNodes; ++n) {
1269:       const PetscInt tag = nodes->id[n];
1270:       PetscCheck(gmsh->nodeMap[tag] < 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Repeated node tag %" PetscInt_FMT, tag);
1271:       gmsh->nodeMap[tag] = n;
1272:     }
1273:   }
1274:   PetscFunctionReturn(PETSC_SUCCESS);
1275: }

1277: static PetscErrorCode GmshReadElements(GmshFile *gmsh, GmshMesh *mesh)
1278: {
1279:   PetscFunctionBegin;
1280:   switch (gmsh->fileFormat) {
1281:   case 41:
1282:     PetscCall(GmshReadElements_v41(gmsh, mesh));
1283:     break;
1284:   case 40:
1285:     PetscCall(GmshReadElements_v40(gmsh, mesh));
1286:     break;
1287:   default:
1288:     PetscCall(GmshReadElements_v22(gmsh, mesh));
1289:     break;
1290:   }

1292:   { /* Reorder elements by codimension and polytope type */
1293:     PetscInt     ne       = mesh->numElems;
1294:     GmshElement *elements = mesh->elements;
1295:     PetscInt     keymap[GMSH_NUM_POLYTOPES], nk = 0;
1296:     PetscInt     offset[GMSH_NUM_POLYTOPES + 1], e, k;

1298:     for (k = 0; k < GMSH_NUM_POLYTOPES; ++k) keymap[k] = PETSC_INT_MIN;
1299:     PetscCall(PetscMemzero(offset, sizeof(offset)));

1301:     keymap[GMSH_TET] = nk++;
1302:     keymap[GMSH_HEX] = nk++;
1303:     keymap[GMSH_PRI] = nk++;
1304:     keymap[GMSH_PYR] = nk++;
1305:     keymap[GMSH_TRI] = nk++;
1306:     keymap[GMSH_QUA] = nk++;
1307:     keymap[GMSH_SEG] = nk++;
1308:     keymap[GMSH_VTX] = nk++;

1310:     PetscCall(GmshElementsCreate(mesh->numElems, &mesh->elements));
1311: #define key(eid) keymap[GmshCellMap[elements[eid].cellType].polytope]
1312:     for (e = 0; e < ne; ++e) offset[1 + key(e)]++;
1313:     for (k = 1; k < nk; ++k) offset[k] += offset[k - 1];
1314:     for (e = 0; e < ne; ++e) mesh->elements[offset[key(e)]++] = elements[e];
1315: #undef key
1316:     PetscCall(GmshElementsDestroy(&elements));
1317:   }

1319:   { /* Mesh dimension and order */
1320:     GmshElement *elem = mesh->numElems ? mesh->elements : NULL;
1321:     mesh->dim         = elem ? GmshCellMap[elem->cellType].dim : 0;
1322:     mesh->order       = elem ? GmshCellMap[elem->cellType].order : 0;
1323:   }

1325:   {
1326:     PetscBT  vtx;
1327:     PetscInt dim = mesh->dim, e, n, v;

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

1331:     /* Compute number of cells and set of vertices */
1332:     mesh->numCells = 0;
1333:     for (e = 0; e < mesh->numElems; ++e) {
1334:       GmshElement *elem = mesh->elements + e;
1335:       if (elem->dim == dim && dim > 0) mesh->numCells++;
1336:       for (v = 0; v < elem->numVerts; v++) PetscCall(PetscBTSet(vtx, elem->nodes[v]));
1337:     }

1339:     /* Compute numbering for vertices */
1340:     mesh->numVerts = 0;
1341:     PetscCall(PetscMalloc1(mesh->numNodes, &mesh->vertexMap));
1342:     for (n = 0; n < mesh->numNodes; ++n) mesh->vertexMap[n] = PetscBTLookup(vtx, n) ? mesh->numVerts++ : PETSC_INT_MIN;

1344:     PetscCall(PetscBTDestroy(&vtx));
1345:   }
1346:   PetscFunctionReturn(PETSC_SUCCESS);
1347: }

1349: static PetscErrorCode GmshReadPeriodic(GmshFile *gmsh, GmshMesh *mesh)
1350: {
1351:   PetscInt n;

1353:   PetscFunctionBegin;
1354:   PetscCall(PetscMalloc1(mesh->numNodes, &mesh->periodMap));
1355:   for (n = 0; n < mesh->numNodes; ++n) mesh->periodMap[n] = n;
1356:   switch (gmsh->fileFormat) {
1357:   case 41:
1358:     PetscCall(GmshReadPeriodic_v41(gmsh, mesh->periodMap));
1359:     break;
1360:   default:
1361:     PetscCall(GmshReadPeriodic_v40(gmsh, mesh->periodMap));
1362:     break;
1363:   }

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

1369:   /* Renumber vertices (filter out correspondings) */
1370:   mesh->numVerts = 0;
1371:   for (n = 0; n < mesh->numNodes; ++n)
1372:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1373:       if (mesh->periodMap[n] == n) /* is primary */
1374:         mesh->vertexMap[n] = mesh->numVerts++;
1375:   for (n = 0; n < mesh->numNodes; ++n)
1376:     if (mesh->vertexMap[n] >= 0)   /* is vertex */
1377:       if (mesh->periodMap[n] != n) /* is corresponding  */
1378:         mesh->vertexMap[n] = mesh->vertexMap[mesh->periodMap[n]];
1379:   PetscFunctionReturn(PETSC_SUCCESS);
1380: }

1382: #define DM_POLYTOPE_VERTEX DM_POLYTOPE_POINT
1383: static const DMPolytopeType DMPolytopeMap[] = {
1384:   /* GMSH_VTX */ DM_POLYTOPE_VERTEX,
1385:   /* GMSH_SEG */ DM_POLYTOPE_SEGMENT,
1386:   /* GMSH_TRI */ DM_POLYTOPE_TRIANGLE,
1387:   /* GMSH_QUA */ DM_POLYTOPE_QUADRILATERAL,
1388:   /* GMSH_TET */ DM_POLYTOPE_TETRAHEDRON,
1389:   /* GMSH_HEX */ DM_POLYTOPE_HEXAHEDRON,
1390:   /* GMSH_PRI */ DM_POLYTOPE_TRI_PRISM,
1391:   /* GMSH_PYR */ DM_POLYTOPE_PYRAMID,       DM_POLYTOPE_UNKNOWN};

1393: static inline DMPolytopeType DMPolytopeTypeFromGmsh(PetscInt cellType)
1394: {
1395:   return DMPolytopeMap[GmshCellMap[cellType].polytope];
1396: }

1398: static PetscErrorCode GmshCreateFE(MPI_Comm comm, const char prefix[], PetscBool isSimplex, PetscBool continuity, PetscDTNodeType nodeType, PetscInt dim, PetscInt Nc, PetscInt k, PetscFE *fem)
1399: {
1400:   DM              K;
1401:   PetscSpace      P;
1402:   PetscDualSpace  Q;
1403:   PetscQuadrature q, fq;
1404:   PetscBool       isTensor = isSimplex ? PETSC_FALSE : PETSC_TRUE;
1405:   PetscBool       endpoint = PETSC_TRUE;
1406:   char            name[32];

1408:   PetscFunctionBegin;
1409:   /* Create space */
1410:   PetscCall(PetscSpaceCreate(comm, &P));
1411:   PetscCall(PetscSpaceSetType(P, PETSCSPACEPOLYNOMIAL));
1412:   PetscCall(PetscSpacePolynomialSetTensor(P, isTensor));
1413:   PetscCall(PetscSpaceSetNumComponents(P, Nc));
1414:   PetscCall(PetscSpaceSetNumVariables(P, dim));
1415:   PetscCall(PetscSpaceSetDegree(P, k, PETSC_DETERMINE));
1416:   if (prefix) {
1417:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, prefix));
1418:     PetscCall(PetscSpaceSetFromOptions(P));
1419:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)P, NULL));
1420:     PetscCall(PetscSpaceGetDegree(P, &k, NULL));
1421:   }
1422:   PetscCall(PetscSpaceSetUp(P));
1423:   /* Create dual space */
1424:   PetscCall(PetscDualSpaceCreate(comm, &Q));
1425:   PetscCall(PetscDualSpaceSetType(Q, PETSCDUALSPACELAGRANGE));
1426:   PetscCall(PetscDualSpaceLagrangeSetTensor(Q, isTensor));
1427:   PetscCall(PetscDualSpaceLagrangeSetContinuity(Q, continuity));
1428:   PetscCall(PetscDualSpaceLagrangeSetNodeType(Q, nodeType, endpoint, 0));
1429:   PetscCall(PetscDualSpaceSetNumComponents(Q, Nc));
1430:   PetscCall(PetscDualSpaceSetOrder(Q, k));
1431:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, DMPolytopeTypeSimpleShape(dim, isSimplex), &K));
1432:   PetscCall(PetscDualSpaceSetDM(Q, K));
1433:   PetscCall(DMDestroy(&K));
1434:   if (prefix) {
1435:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, prefix));
1436:     PetscCall(PetscDualSpaceSetFromOptions(Q));
1437:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Q, NULL));
1438:   }
1439:   PetscCall(PetscDualSpaceSetUp(Q));
1440:   /* Create quadrature */
1441:   if (isSimplex) {
1442:     PetscCall(PetscDTStroudConicalQuadrature(dim, 1, k + 1, -1, +1, &q));
1443:     PetscCall(PetscDTStroudConicalQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1444:   } else {
1445:     PetscCall(PetscDTGaussTensorQuadrature(dim, 1, k + 1, -1, +1, &q));
1446:     PetscCall(PetscDTGaussTensorQuadrature(dim - 1, 1, k + 1, -1, +1, &fq));
1447:   }
1448:   /* Create finite element */
1449:   PetscCall(PetscFECreate(comm, fem));
1450:   PetscCall(PetscFESetType(*fem, PETSCFEBASIC));
1451:   PetscCall(PetscFESetNumComponents(*fem, Nc));
1452:   PetscCall(PetscFESetBasisSpace(*fem, P));
1453:   PetscCall(PetscFESetDualSpace(*fem, Q));
1454:   PetscCall(PetscFESetQuadrature(*fem, q));
1455:   PetscCall(PetscFESetFaceQuadrature(*fem, fq));
1456:   if (prefix) {
1457:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, prefix));
1458:     PetscCall(PetscFESetFromOptions(*fem));
1459:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*fem, NULL));
1460:   }
1461:   PetscCall(PetscFESetUp(*fem));
1462:   /* Cleanup */
1463:   PetscCall(PetscSpaceDestroy(&P));
1464:   PetscCall(PetscDualSpaceDestroy(&Q));
1465:   PetscCall(PetscQuadratureDestroy(&q));
1466:   PetscCall(PetscQuadratureDestroy(&fq));
1467:   /* Set finite element name */
1468:   PetscCall(PetscSNPrintf(name, sizeof(name), "%s%" PetscInt_FMT, isSimplex ? "P" : "Q", k));
1469:   PetscCall(PetscFESetName(*fem, name));
1470:   PetscFunctionReturn(PETSC_SUCCESS);
1471: }

1473: /*@
1474:   DMPlexCreateGmshFromFile - Create a `DMPLEX` mesh from a Gmsh file

1476:   Input Parameters:
1477: + comm        - The MPI communicator
1478: . filename    - Name of the Gmsh file
1479: - interpolate - Create faces and edges in the mesh

1481:   Output Parameter:
1482: . dm - The `DM` object representing the mesh

1484:   Level: beginner

1486: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
1487: @*/
1488: PetscErrorCode DMPlexCreateGmshFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
1489: {
1490:   PetscViewer     viewer;
1491:   PetscMPIInt     rank;
1492:   int             fileType;
1493:   PetscViewerType vtype;

1495:   PetscFunctionBegin;
1496:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

1498:   /* Determine Gmsh file type (ASCII or binary) from file header */
1499:   if (rank == 0) {
1500:     GmshFile gmsh[1];
1501:     char     line[PETSC_MAX_PATH_LEN];
1502:     int      snum;
1503:     float    version;
1504:     int      fileFormat;

1506:     PetscCall(PetscArrayzero(gmsh, 1));
1507:     PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &gmsh->viewer));
1508:     PetscCall(PetscViewerSetType(gmsh->viewer, PETSCVIEWERASCII));
1509:     PetscCall(PetscViewerFileSetMode(gmsh->viewer, FILE_MODE_READ));
1510:     PetscCall(PetscViewerFileSetName(gmsh->viewer, filename));
1511:     /* Read only the first two lines of the Gmsh file */
1512:     PetscCall(GmshReadSection(gmsh, line));
1513:     PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1514:     PetscCall(GmshReadString(gmsh, line, 2));
1515:     snum       = sscanf(line, "%f %d", &version, &fileType);
1516:     fileFormat = (int)roundf(version * 10);
1517:     PetscCheck(snum == 2, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Unable to parse Gmsh file header: %s", line);
1518:     PetscCheck(fileFormat >= 22, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at least 2.2", (double)version);
1519:     PetscCheck((int)version != 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f not supported", (double)version);
1520:     PetscCheck(fileFormat <= 41, PETSC_COMM_SELF, PETSC_ERR_SUP, "Gmsh file version %3.1f must be at most 4.1", (double)version);
1521:     PetscCall(PetscViewerDestroy(&gmsh->viewer));
1522:   }
1523:   PetscCallMPI(MPI_Bcast(&fileType, 1, MPI_INT, 0, comm));
1524:   vtype = (fileType == 0) ? PETSCVIEWERASCII : PETSCVIEWERBINARY;

1526:   /* Create appropriate viewer and build plex */
1527:   PetscCall(PetscViewerCreate(comm, &viewer));
1528:   PetscCall(PetscViewerSetType(viewer, vtype));
1529:   PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
1530:   PetscCall(PetscViewerFileSetName(viewer, filename));
1531:   PetscCall(DMPlexCreateGmsh(comm, viewer, interpolate, dm));
1532:   PetscCall(PetscViewerDestroy(&viewer));
1533:   PetscFunctionReturn(PETSC_SUCCESS);
1534: }

1536: /*@
1537:   DMPlexCreateGmsh - Create a `DMPLEX` mesh from a Gmsh file viewer

1539:   Collective

1541:   Input Parameters:
1542: + comm        - The MPI communicator
1543: . viewer      - The `PetscViewer` associated with a Gmsh file
1544: - interpolate - Create faces and edges in the mesh

1546:   Output Parameter:
1547: . dm - The `DM` object representing the mesh

1549:   Options Database Keys:
1550: + -dm_plex_gmsh_hybrid               - Force triangular prisms to use tensor order
1551: . -dm_plex_gmsh_periodic             - Read Gmsh periodic section and construct a periodic Plex
1552: . -dm_plex_gmsh_highorder            - Generate high-order coordinates
1553: . -dm_plex_gmsh_project              - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
1554: . -dm_plex_gmsh_use_generic          - Generate generic labels, i.e. Cell Sets, Face Sets, etc.
1555: . -dm_plex_gmsh_use_regions          - Generate labels with region names
1556: . -dm_plex_gmsh_mark_vertices        - Add vertices to generated labels
1557: . -dm_plex_gmsh_mark_vertices_strict - Add vertices included in a region to generated labels
1558: . -dm_plex_gmsh_multiple_tags        - Allow multiple tags for default labels
1559: - -dm_plex_gmsh_spacedim <d>         - Embedding space dimension, if different from topological dimension

1561:   Level: beginner

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

1566:   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. Alternatively, `-dm_plex_gmsh_use_regions` creates labels based on the region names from the `PhysicalNames` section, and all tags are used, but you can retain the generic labels using `-dm_plex_gmsh_use_generic`.

1568: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMCreate()`
1569: @*/
1570: PetscErrorCode DMPlexCreateGmsh(MPI_Comm comm, PetscViewer viewer, PetscBool interpolate, DM *dm)
1571: {
1572:   GmshMesh    *mesh          = NULL;
1573:   PetscViewer  parentviewer  = NULL;
1574:   PetscBT      periodicVerts = NULL;
1575:   PetscBT     *periodicCells = NULL;
1576:   DM           cdm, cdmCell = NULL;
1577:   PetscSection cs, csCell   = NULL;
1578:   Vec          coordinates, coordinatesCell;
1579:   DMLabel      cellSets = NULL, faceSets = NULL, edgeSets = NULL, vertSets = NULL, marker = NULL, *regionSets;
1580:   PetscInt     dim = 0, coordDim = -1, order = 0, maxHeight = 0;
1581:   PetscInt     numNodes = 0, numElems = 0, numVerts = 0, numCells = 0, vStart, vEnd;
1582:   PetscInt     cell, cone[8], e, n, v, d;
1583:   PetscBool    usegeneric = PETSC_TRUE, useregions = PETSC_FALSE, markvertices = PETSC_FALSE, markverticesstrict = PETSC_FALSE, multipleTags = PETSC_FALSE;
1584:   PetscBool    flg, binary, hybrid = interpolate, periodic = PETSC_TRUE;
1585:   PetscBool    highOrder = PETSC_TRUE, highOrderSet, project = PETSC_FALSE;
1586:   PetscBool    isSimplex = PETSC_FALSE, isHybrid = PETSC_FALSE, hasTetra = PETSC_FALSE;
1587:   PetscMPIInt  rank;

1589:   PetscFunctionBegin;
1590:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
1591:   PetscObjectOptionsBegin((PetscObject)viewer);
1592:   PetscOptionsHeadBegin(PetscOptionsObject, "DMPlex Gmsh options");
1593:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_hybrid", "Force triangular prisms to use tensor order", "DMPlexCreateGmsh", hybrid, &hybrid, NULL));
1594:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_periodic", "Read Gmsh periodic section", "DMPlexCreateGmsh", periodic, &periodic, NULL));
1595:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_highorder", "Generate high-order coordinates", "DMPlexCreateGmsh", highOrder, &highOrder, &highOrderSet));
1596:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_project", "Project high-order coordinates to a different space", "DMPlexCreateGmsh", project, &project, NULL));
1597:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_regions", "Generate labels with region names", "DMPlexCreateGmsh", useregions, &useregions, NULL));
1598:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_use_generic", "Generate generic labels, i.e. Cell Sets, Face Sets, etc", "DMPlexCreateGmsh", usegeneric, &usegeneric, &flg));
1599:   if (!flg && useregions) usegeneric = PETSC_FALSE;
1600:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices", "Add vertices to generated labels", "DMPlexCreateGmsh", markvertices, &markvertices, NULL));
1601:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_mark_vertices_strict", "Add only directly tagged vertices to generated labels", "DMPlexCreateGmsh", markverticesstrict, &markverticesstrict, NULL));
1602:   PetscCall(PetscOptionsBool("-dm_plex_gmsh_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL));
1603:   PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE));
1604:   PetscCall(PetscOptionsBoundedInt("-dm_localize_height", "Localize edges and faces in addition to cells", "", maxHeight, &maxHeight, NULL, 0));
1605:   PetscOptionsHeadEnd();
1606:   PetscOptionsEnd();

1608:   PetscCall(GmshCellInfoSetUp());

1610:   PetscCall(DMCreate(comm, dm));
1611:   PetscCall(DMSetType(*dm, DMPLEX));
1612:   PetscCall(PetscLogEventBegin(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));

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

1616:   /* Binary viewers read on all ranks, get subviewer to read only in rank 0 */
1617:   if (binary) {
1618:     parentviewer = viewer;
1619:     PetscCall(PetscViewerGetSubViewer(parentviewer, PETSC_COMM_SELF, &viewer));
1620:   }

1622:   if (rank == 0) {
1623:     GmshFile  gmsh[1];
1624:     char      line[PETSC_MAX_PATH_LEN];
1625:     PetscBool match;

1627:     PetscCall(PetscArrayzero(gmsh, 1));
1628:     gmsh->viewer = viewer;
1629:     gmsh->binary = binary;

1631:     PetscCall(GmshMeshCreate(&mesh));

1633:     /* Read mesh format */
1634:     PetscCall(GmshReadSection(gmsh, line));
1635:     PetscCall(GmshExpect(gmsh, "$MeshFormat", line));
1636:     PetscCall(GmshReadMeshFormat(gmsh));
1637:     PetscCall(GmshReadEndSection(gmsh, "$EndMeshFormat", line));

1639:     /* OPTIONAL Read mesh version (Neper only) */
1640:     PetscCall(GmshReadSection(gmsh, line));
1641:     PetscCall(GmshMatch(gmsh, "$MeshVersion", line, &match));
1642:     if (match) {
1643:       PetscCall(GmshExpect(gmsh, "$MeshVersion", line));
1644:       PetscCall(GmshReadMeshVersion(gmsh));
1645:       PetscCall(GmshReadEndSection(gmsh, "$EndMeshVersion", line));
1646:       /* Initial read for entity section */
1647:       PetscCall(GmshReadSection(gmsh, line));
1648:     }

1650:     /* OPTIONAL Read mesh domain (Neper only) */
1651:     PetscCall(GmshMatch(gmsh, "$Domain", line, &match));
1652:     if (match) {
1653:       PetscCall(GmshExpect(gmsh, "$Domain", line));
1654:       PetscCall(GmshReadMeshDomain(gmsh));
1655:       PetscCall(GmshReadEndSection(gmsh, "$EndDomain", line));
1656:       /* Initial read for entity section */
1657:       PetscCall(GmshReadSection(gmsh, line));
1658:     }

1660:     /* OPTIONAL Read physical names */
1661:     PetscCall(GmshMatch(gmsh, "$PhysicalNames", line, &match));
1662:     if (match) {
1663:       PetscCall(GmshExpect(gmsh, "$PhysicalNames", line));
1664:       PetscCall(GmshReadPhysicalNames(gmsh, mesh));
1665:       PetscCall(GmshReadEndSection(gmsh, "$EndPhysicalNames", line));
1666:       /* Initial read for entity section */
1667:       PetscCall(GmshReadSection(gmsh, line));
1668:     }

1670:     /* OPTIONAL Read entities */
1671:     if (gmsh->fileFormat >= 40) {
1672:       PetscCall(GmshMatch(gmsh, "$Entities", line, &match));
1673:       if (match) {
1674:         PetscCall(GmshExpect(gmsh, "$Entities", line));
1675:         PetscCall(GmshReadEntities(gmsh, mesh));
1676:         PetscCall(GmshReadEndSection(gmsh, "$EndEntities", line));
1677:         /* Initial read for nodes section */
1678:         PetscCall(GmshReadSection(gmsh, line));
1679:       }
1680:     }

1682:     /* Read nodes */
1683:     PetscCall(GmshExpect(gmsh, "$Nodes", line));
1684:     PetscCall(GmshReadNodes(gmsh, mesh));
1685:     PetscCall(GmshReadEndSection(gmsh, "$EndNodes", line));

1687:     /* Read elements */
1688:     PetscCall(GmshReadSection(gmsh, line));
1689:     PetscCall(GmshExpect(gmsh, "$Elements", line));
1690:     PetscCall(GmshReadElements(gmsh, mesh));
1691:     PetscCall(GmshReadEndSection(gmsh, "$EndElements", line));

1693:     /* Read periodic section (OPTIONAL) */
1694:     if (periodic) {
1695:       PetscCall(GmshReadSection(gmsh, line));
1696:       PetscCall(GmshMatch(gmsh, "$Periodic", line, &periodic));
1697:     }
1698:     if (periodic) {
1699:       PetscCall(GmshExpect(gmsh, "$Periodic", line));
1700:       PetscCall(GmshReadPeriodic(gmsh, mesh));
1701:       PetscCall(GmshReadEndSection(gmsh, "$EndPeriodic", line));
1702:     }

1704:     PetscCall(PetscFree(gmsh->wbuf));
1705:     PetscCall(PetscFree(gmsh->sbuf));
1706:     PetscCall(PetscFree(gmsh->nbuf));

1708:     dim      = mesh->dim;
1709:     order    = mesh->order;
1710:     numNodes = mesh->numNodes;
1711:     numElems = mesh->numElems;
1712:     numVerts = mesh->numVerts;
1713:     numCells = mesh->numCells;

1715:     {
1716:       GmshElement *elemA = mesh->numCells > 0 ? mesh->elements : NULL;
1717:       GmshElement *elemB = PetscSafePointerPlusOffset(elemA, mesh->numCells - 1);
1718:       int          ptA   = elemA ? GmshCellMap[elemA->cellType].polytope : -1;
1719:       int          ptB   = elemB ? GmshCellMap[elemB->cellType].polytope : -1;
1720:       isSimplex          = (ptA == GMSH_QUA || ptA == GMSH_HEX) ? PETSC_FALSE : PETSC_TRUE;
1721:       isHybrid           = (ptA == ptB) ? PETSC_FALSE : PETSC_TRUE;
1722:       hasTetra           = (ptA == GMSH_TET) ? PETSC_TRUE : PETSC_FALSE;
1723:     }
1724:   }

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

1728:   {
1729:     int buf[6];

1731:     buf[0] = (int)dim;
1732:     buf[1] = (int)order;
1733:     buf[2] = periodic;
1734:     buf[3] = isSimplex;
1735:     buf[4] = isHybrid;
1736:     buf[5] = hasTetra;

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

1740:     dim       = buf[0];
1741:     order     = buf[1];
1742:     periodic  = buf[2] ? PETSC_TRUE : PETSC_FALSE;
1743:     isSimplex = buf[3] ? PETSC_TRUE : PETSC_FALSE;
1744:     isHybrid  = buf[4] ? PETSC_TRUE : PETSC_FALSE;
1745:     hasTetra  = buf[5] ? PETSC_TRUE : PETSC_FALSE;
1746:   }

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

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

1754:   /* Allocate the cell-vertex mesh */
1755:   PetscCall(DMPlexSetChart(*dm, 0, numCells + numVerts));
1756:   for (cell = 0; cell < numCells; ++cell) {
1757:     GmshElement   *elem  = mesh->elements + cell;
1758:     DMPolytopeType ctype = DMPolytopeTypeFromGmsh(elem->cellType);
1759:     if (hybrid && hasTetra && ctype == DM_POLYTOPE_TRI_PRISM) ctype = DM_POLYTOPE_TRI_PRISM_TENSOR;
1760:     PetscCall(DMPlexSetConeSize(*dm, cell, elem->numVerts));
1761:     PetscCall(DMPlexSetCellType(*dm, cell, ctype));
1762:   }
1763:   for (v = numCells; v < numCells + numVerts; ++v) PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT));
1764:   PetscCall(DMSetUp(*dm));

1766:   /* Add cell-vertex connections */
1767:   for (cell = 0; cell < numCells; ++cell) {
1768:     GmshElement *elem = mesh->elements + cell;
1769:     for (v = 0; v < elem->numVerts; ++v) {
1770:       const PetscInt nn = elem->nodes[v];
1771:       const PetscInt vv = mesh->vertexMap[nn];
1772:       cone[v]           = numCells + vv;
1773:     }
1774:     PetscCall(DMPlexReorderCell(*dm, cell, cone));
1775:     PetscCall(DMPlexSetCone(*dm, cell, cone));
1776:   }

1778:   PetscCall(DMSetDimension(*dm, dim));
1779:   PetscCall(DMPlexSymmetrize(*dm));
1780:   PetscCall(DMPlexStratify(*dm));
1781:   if (interpolate) {
1782:     DM idm;

1784:     PetscCall(DMPlexInterpolate(*dm, &idm));
1785:     PetscCall(DMDestroy(dm));
1786:     *dm = idm;
1787:   }
1788:   PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd));

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

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

1797:       /* Create cell sets */
1798:       if (elem->dim == dim && dim > 0) {
1799:         if (elem->numTags > 0) {
1800:           PetscInt Nt = elem->numTags, t, r;

1802:           for (t = 0; t < Nt; ++t) {
1803:             const PetscInt  tag     = elem->tags[t];
1804:             const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1806:             if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &cellSets, "Cell Sets", cell, tag));
1807:             for (r = 0; r < Nr; ++r) {
1808:               if (mesh->regionDims[r] != dim) continue;
1809:               if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], cell, tag));
1810:             }
1811:           }
1812:         }
1813:         cell++;
1814:       }

1816:       /* Create face sets */
1817:       if (elem->numTags && interpolate && elem->dim == dim - 1) {
1818:         PetscInt        joinSize;
1819:         const PetscInt *join = NULL;
1820:         PetscInt        Nt   = elem->numTags, pdepth;

1822:         /* Find the relevant facet with vertex joins */
1823:         for (v = 0; v < elem->numVerts; ++v) {
1824:           const PetscInt nn = elem->nodes[v];
1825:           const PetscInt vv = mesh->vertexMap[nn];
1826:           cone[v]           = vStart + vv;
1827:         }
1828:         PetscCall(DMPlexGetFullJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1829:         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);
1830:         PetscCall(DMPlexGetPointDepth(*dm, join[0], &pdepth));
1831:         PetscCheck(pdepth == dim - 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Plex facet %" PetscInt_FMT " for Gmsh element %" PetscInt_FMT " had depth %" PetscInt_FMT " != %" PetscInt_FMT, join[0], elem->id, pdepth, dim - 1);
1832:         for (PetscInt t = 0; t < Nt; ++t) {
1833:           const PetscInt  tag     = elem->tags[t];
1834:           const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1836:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &faceSets, "Face Sets", join[0], tag));
1837:           for (PetscInt r = 0; r < Nr; ++r) {
1838:             if (mesh->regionDims[r] != dim - 1) continue;
1839:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], join[0], tag));
1840:           }
1841:         }
1842:         PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1843:       }

1845:       /* Create edge sets */
1846:       if (elem->numTags && interpolate && dim > 2 && elem->dim == 1) {
1847:         PetscInt        joinSize;
1848:         const PetscInt *join = NULL;
1849:         PetscInt        Nt   = elem->numTags, t, r;

1851:         /* Find the relevant edge with vertex joins */
1852:         for (v = 0; v < elem->numVerts; ++v) {
1853:           const PetscInt nn = elem->nodes[v];
1854:           const PetscInt vv = mesh->vertexMap[nn];
1855:           cone[v]           = vStart + vv;
1856:         }
1857:         PetscCall(DMPlexGetJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1858:         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);
1859:         for (t = 0; t < Nt; ++t) {
1860:           const PetscInt  tag     = elem->tags[t];
1861:           const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1863:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &edgeSets, "Edge Sets", join[0], tag));
1864:           for (r = 0; r < Nr; ++r) {
1865:             if (mesh->regionDims[r] != 1) continue;
1866:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], join[0], tag));
1867:           }
1868:         }
1869:         PetscCall(DMPlexRestoreJoin(*dm, elem->numVerts, cone, &joinSize, &join));
1870:       }

1872:       /* Create vertex sets */
1873:       if (elem->numTags && elem->dim == 0 && (markverticesstrict || markvertices)) {
1874:         const PetscInt nn = elem->nodes[0];
1875:         const PetscInt vv = mesh->vertexMap[nn];
1876:         PetscInt       Nt = elem->numTags;

1878:         for (PetscInt t = 0; t < Nt; ++t) {
1879:           const PetscInt tag = elem->tags[t];

1881:           if (vv < 0) continue;
1882:           if (usegeneric) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1883:           for (PetscInt r = 0; r < Nr; ++r) {
1884:             if (mesh->regionDims[r] != 0) continue;
1885:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag));
1886:           }
1887:         }
1888:       }
1889:     }
1890:     if (markvertices) {
1891:       for (v = 0; v < numNodes; ++v) {
1892:         const PetscInt  vv   = mesh->vertexMap[v];
1893:         const PetscInt *tags = &mesh->nodelist->tag[v * GMSH_MAX_TAGS];
1894:         PetscInt        r, t;

1896:         if (vv < 0) continue;
1897:         for (t = 0; t < GMSH_MAX_TAGS; ++t) {
1898:           const PetscInt  tag     = tags[t];
1899:           const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

1901:           if (tag == -1) continue;
1902:           if (generic) PetscCall(DMSetLabelValue_Fast(*dm, &vertSets, "Vertex Sets", vStart + vv, tag));
1903:           for (r = 0; r < Nr; ++r) {
1904:             if (mesh->regionTags[r] == tag) PetscCall(DMSetLabelValue_Fast(*dm, &regionSets[r], mesh->regionNames[r], vStart + vv, tag));
1905:           }
1906:         }
1907:       }
1908:     }
1909:     PetscCall(PetscFree(regionSets));
1910:   }

1912:   { /* Create Cell/Face/Edge/Vertex Sets labels at all processes */
1913:     enum {
1914:       n = 5
1915:     };
1916:     PetscBool flag[n];

1918:     flag[0] = cellSets ? PETSC_TRUE : PETSC_FALSE;
1919:     flag[1] = faceSets ? PETSC_TRUE : PETSC_FALSE;
1920:     flag[2] = edgeSets ? PETSC_TRUE : PETSC_FALSE;
1921:     flag[3] = vertSets ? PETSC_TRUE : PETSC_FALSE;
1922:     flag[4] = marker ? PETSC_TRUE : PETSC_FALSE;
1923:     PetscCallMPI(MPI_Bcast(flag, n, MPI_C_BOOL, 0, comm));
1924:     if (flag[0]) PetscCall(DMCreateLabel(*dm, "Cell Sets"));
1925:     if (flag[1]) PetscCall(DMCreateLabel(*dm, "Face Sets"));
1926:     if (flag[2]) PetscCall(DMCreateLabel(*dm, "Edge Sets"));
1927:     if (flag[3]) PetscCall(DMCreateLabel(*dm, "Vertex Sets"));
1928:     if (flag[4]) PetscCall(DMCreateLabel(*dm, "marker"));
1929:   }

1931:   if (periodic) {
1932:     PetscCall(PetscBTCreate(numVerts, &periodicVerts));
1933:     for (n = 0; n < numNodes; ++n) {
1934:       if (mesh->vertexMap[n] >= 0) {
1935:         if (PetscUnlikely(mesh->periodMap[n] != n)) {
1936:           PetscInt m = mesh->periodMap[n];
1937:           PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[n]));
1938:           PetscCall(PetscBTSet(periodicVerts, mesh->vertexMap[m]));
1939:         }
1940:       }
1941:     }
1942:     PetscCall(DMGetCoordinateDM(*dm, &cdm));
1943:     PetscCall(PetscMalloc1(maxHeight + 1, &periodicCells));
1944:     for (PetscInt h = 0; h <= maxHeight; ++h) {
1945:       PetscInt pStart, pEnd;

1947:       PetscCall(DMPlexGetHeightStratum(*dm, h, &pStart, &pEnd));
1948:       PetscCall(PetscBTCreate(pEnd - pStart, &periodicCells[h]));
1949:       for (PetscInt p = pStart; p < pEnd; ++p) {
1950:         PetscInt *closure = NULL;
1951:         PetscInt  Ncl;

1953:         PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1954:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
1955:           if (closure[cl] >= vStart && closure[cl] < vEnd) {
1956:             if (PetscUnlikely(PetscBTLookup(periodicVerts, closure[cl] - vStart))) {
1957:               PetscCall(PetscBTSet(periodicCells[h], p - pStart));
1958:               break;
1959:             }
1960:           }
1961:         }
1962:         PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
1963:       }
1964:     }
1965:   }

1967:   /* Setup coordinate DM */
1968:   if (coordDim < 0) coordDim = dim;
1969:   PetscCall(DMSetCoordinateDim(*dm, coordDim));
1970:   PetscCall(DMGetCoordinateDM(*dm, &cdm));
1971:   if (highOrder) {
1972:     PetscFE         fe;
1973:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
1974:     PetscDTNodeType nodeType   = PETSCDTNODES_EQUISPACED;

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

1978:     PetscCall(GmshCreateFE(comm, NULL, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
1979:     PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_fe_view"));
1980:     PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)fe));
1981:     PetscCall(PetscFEDestroy(&fe));
1982:     PetscCall(DMCreateDS(cdm));
1983:   }
1984:   if (periodic) {
1985:     PetscCall(DMClone(cdm, &cdmCell));
1986:     PetscCall(DMSetCellCoordinateDM(*dm, cdmCell));
1987:   }

1989:   /* Create coordinates */
1990:   if (highOrder) {
1991:     PetscInt     maxDof = GmshNumNodes_HEX(order) * coordDim;
1992:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
1993:     PetscSection section;
1994:     PetscScalar *cellCoords;

1996:     PetscCall(DMSetLocalSection(cdm, NULL));
1997:     PetscCall(DMGetLocalSection(cdm, &cs));
1998:     PetscCall(PetscSectionClone(cs, &section));
1999:     PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */

2001:     PetscCall(DMCreateLocalVector(cdm, &coordinates));
2002:     PetscCall(PetscMalloc1(maxDof, &cellCoords));
2003:     for (cell = 0; cell < numCells; ++cell) {
2004:       GmshElement *elem     = mesh->elements + cell;
2005:       const int   *lexorder = GmshCellMap[elem->cellType].lexorder();
2006:       int          s        = 0;
2007:       for (n = 0; n < elem->numNodes; ++n) {
2008:         while (lexorder[n + s] < 0) ++s;
2009:         const PetscInt node = elem->nodes[lexorder[n + s]];
2010:         for (d = 0; d < coordDim; ++d) cellCoords[(n + s) * coordDim + d] = (PetscReal)coords[node * 3 + d];
2011:       }
2012:       if (s) {
2013:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
2014:         PetscReal quaCenterWeights[9] = {-0.25, 0.5, -0.25, 0.5, 0.0, 0.5, -0.25, 0.5, -0.25};
2015:         /* For the coordinate mapping we weight vertices by -1/4 and edges by 1/2, which we get from Q_2 interpolation */
2016:         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};
2017:         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};
2018:         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};
2019:         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};
2020:         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};
2021:         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};
2022:         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};
2023:         PetscReal  *sdWeights2[9]        = {NULL, NULL, NULL, NULL, quaCenterWeights, NULL, NULL, NULL, NULL};
2024:         PetscReal  *sdWeights3[27]       = {NULL, NULL, NULL, NULL, hexBottomWeights, NULL, NULL, NULL, NULL, NULL, hexFrontWeights, NULL, hexLeftWeights, hexCenterWeights, hexRightWeights, NULL, hexBackWeights, NULL,
2025:                                             NULL, NULL, NULL, NULL, hexTopWeights,    NULL, NULL, NULL, NULL};
2026:         PetscReal **sdWeights[4]         = {NULL, NULL, sdWeights2, sdWeights3};

2028:         /* Missing entries in serendipity cell, only works for 8-node quad and 20-node hex */
2029:         for (n = 0; n < elem->numNodes + s; ++n) {
2030:           if (lexorder[n] >= 0) continue;
2031:           for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] = 0.0;
2032:           for (int bn = 0; bn < elem->numNodes + s; ++bn) {
2033:             if (lexorder[bn] < 0) continue;
2034:             const PetscReal *weights = sdWeights[coordDim][n];
2035:             const PetscInt   bnode   = elem->nodes[lexorder[bn]];
2036:             for (d = 0; d < coordDim; ++d) cellCoords[n * coordDim + d] += weights[bn] * (PetscReal)coords[bnode * 3 + d];
2037:           }
2038:         }
2039:       }
2040:       PetscCall(DMPlexVecSetClosure(cdm, section, coordinates, cell, cellCoords, INSERT_VALUES));
2041:     }
2042:     PetscCall(PetscSectionDestroy(&section));
2043:     PetscCall(PetscFree(cellCoords));
2044:   } else {
2045:     PetscInt    *nodeMap;
2046:     double      *coords = mesh ? mesh->nodelist->xyz : NULL;
2047:     PetscScalar *pointCoords;

2049:     PetscCall(DMGetCoordinateSection(*dm, &cs));
2050:     PetscCall(PetscSectionSetNumFields(cs, 1));
2051:     PetscCall(PetscSectionSetFieldComponents(cs, 0, coordDim));
2052:     PetscCall(PetscSectionSetChart(cs, numCells, numCells + numVerts));
2053:     for (v = numCells; v < numCells + numVerts; ++v) {
2054:       PetscCall(PetscSectionSetDof(cs, v, coordDim));
2055:       PetscCall(PetscSectionSetFieldDof(cs, v, 0, coordDim));
2056:     }
2057:     PetscCall(PetscSectionSetUp(cs));

2059:     /* We need to localize coordinates on cells */
2060:     if (periodic) {
2061:       PetscInt newStart = PETSC_INT_MAX, newEnd = -1, pStart, pEnd;

2063:       PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)cdmCell), &csCell));
2064:       PetscCall(PetscSectionSetNumFields(csCell, 1));
2065:       PetscCall(PetscSectionSetFieldComponents(csCell, 0, coordDim));
2066:       for (PetscInt h = 0; h <= maxHeight; h++) {
2067:         PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
2068:         newStart = PetscMin(newStart, pStart);
2069:         newEnd   = PetscMax(newEnd, pEnd);
2070:       }
2071:       PetscCall(PetscSectionSetChart(csCell, newStart, newEnd));
2072:       for (PetscInt h = 0; h <= maxHeight; h++) {
2073:         PetscCall(DMPlexGetHeightStratum(cdmCell, h, &pStart, &pEnd));
2074:         for (PetscInt p = pStart; p < pEnd; ++p) {
2075:           PetscInt *closure = NULL;
2076:           PetscInt  Ncl, Nv = 0;

2078:           if (PetscUnlikely(PetscBTLookup(periodicCells[h], p - pStart))) {
2079:             PetscCall(DMPlexGetTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
2080:             for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2081:               if (closure[cl] >= vStart && closure[cl] < vEnd) ++Nv;
2082:             }
2083:             PetscCall(DMPlexRestoreTransitiveClosure(*dm, p, PETSC_TRUE, &Ncl, &closure));
2084:             PetscCall(PetscSectionSetDof(csCell, p, Nv * coordDim));
2085:             PetscCall(PetscSectionSetFieldDof(csCell, p, 0, Nv * coordDim));
2086:           }
2087:         }
2088:       }
2089:       PetscCall(PetscSectionSetUp(csCell));
2090:       PetscCall(DMSetCellCoordinateSection(*dm, PETSC_DETERMINE, csCell));
2091:     }

2093:     PetscCall(DMCreateLocalVector(cdm, &coordinates));
2094:     PetscCall(VecGetArray(coordinates, &pointCoords));
2095:     PetscCall(PetscMalloc1(numVerts, &nodeMap));
2096:     for (n = 0; n < numNodes; n++)
2097:       if (mesh->vertexMap[n] >= 0) nodeMap[mesh->vertexMap[n]] = n;
2098:     for (v = 0; v < numVerts; ++v) {
2099:       PetscInt off, node = nodeMap[v];

2101:       PetscCall(PetscSectionGetOffset(cs, numCells + v, &off));
2102:       for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
2103:     }
2104:     PetscCall(VecRestoreArray(coordinates, &pointCoords));

2106:     if (periodic) {
2107:       PetscInt cStart, cEnd;

2109:       PetscCall(DMPlexGetHeightStratum(cdmCell, 0, &cStart, &cEnd));
2110:       PetscCall(DMCreateLocalVector(cdmCell, &coordinatesCell));
2111:       PetscCall(VecGetArray(coordinatesCell, &pointCoords));
2112:       for (PetscInt c = cStart; c < cEnd; ++c) {
2113:         GmshElement *elem    = mesh->elements + c - cStart;
2114:         PetscInt    *closure = NULL;
2115:         PetscInt     verts[8];
2116:         PetscInt     Ncl, Nv = 0;

2118:         for (PetscInt v = 0; v < elem->numVerts; ++v) cone[v] = elem->nodes[v];
2119:         PetscCall(DMPlexReorderCell(cdmCell, c, cone));
2120:         PetscCall(DMPlexGetTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
2121:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2122:           if (closure[cl] >= vStart && closure[cl] < vEnd) verts[Nv++] = closure[cl];
2123:         }
2124:         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);
2125:         for (PetscInt cl = 0; cl < Ncl * 2; cl += 2) {
2126:           const PetscInt point = closure[cl];
2127:           PetscInt       hStart, h;

2129:           PetscCall(DMPlexGetPointHeight(cdmCell, point, &h));
2130:           if (h > maxHeight) continue;
2131:           PetscCall(DMPlexGetHeightStratum(cdmCell, h, &hStart, NULL));
2132:           if (PetscUnlikely(PetscBTLookup(periodicCells[h], point - hStart))) {
2133:             PetscInt *pclosure = NULL;
2134:             PetscInt  Npcl, off, v;

2136:             PetscCall(PetscSectionGetOffset(csCell, point, &off));
2137:             PetscCall(DMPlexGetTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2138:             for (PetscInt pcl = 0; pcl < Npcl * 2; pcl += 2) {
2139:               if (pclosure[pcl] >= vStart && pclosure[pcl] < vEnd) {
2140:                 for (v = 0; v < Nv; ++v)
2141:                   if (verts[v] == pclosure[pcl]) break;
2142:                 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);
2143:                 for (PetscInt d = 0; d < coordDim; ++d) pointCoords[off++] = (PetscReal)coords[cone[v] * 3 + d];
2144:                 ++v;
2145:               }
2146:             }
2147:             PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, point, PETSC_TRUE, &Npcl, &pclosure));
2148:           }
2149:         }
2150:         PetscCall(DMPlexRestoreTransitiveClosure(cdmCell, c, PETSC_TRUE, &Ncl, &closure));
2151:       }
2152:       PetscCall(VecSetBlockSize(coordinatesCell, coordDim));
2153:       PetscCall(VecRestoreArray(coordinatesCell, &pointCoords));
2154:       PetscCall(DMSetCellCoordinatesLocal(*dm, coordinatesCell));
2155:       PetscCall(VecDestroy(&coordinatesCell));
2156:     }
2157:     PetscCall(PetscFree(nodeMap));
2158:     PetscCall(PetscSectionDestroy(&csCell));
2159:     PetscCall(DMDestroy(&cdmCell));
2160:   }

2162:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
2163:   PetscCall(VecSetBlockSize(coordinates, coordDim));
2164:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
2165:   PetscCall(VecDestroy(&coordinates));

2167:   PetscCall(GmshMeshDestroy(&mesh));
2168:   PetscCall(PetscBTDestroy(&periodicVerts));
2169:   if (periodic) {
2170:     for (PetscInt h = 0; h <= maxHeight; ++h) PetscCall(PetscBTDestroy(periodicCells + h));
2171:     PetscCall(PetscFree(periodicCells));
2172:   }

2174:   if (highOrder && project) {
2175:     PetscFE         fe;
2176:     const char      prefix[]   = "dm_plex_gmsh_project_";
2177:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
2178:     PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;

2180:     if (isSimplex) continuity = PETSC_FALSE; /* XXX FIXME Requires DMPlexSetClosurePermutationLexicographic() */
2181:     PetscCall(GmshCreateFE(comm, prefix, isSimplex, continuity, nodeType, dim, coordDim, order, &fe));
2182:     PetscCall(PetscFEViewFromOptions(fe, NULL, "-dm_plex_gmsh_project_fe_view"));
2183:     PetscCall(DMSetCoordinateDisc(*dm, fe, PETSC_FALSE, PETSC_TRUE));
2184:     PetscCall(PetscFEDestroy(&fe));
2185:   }

2187:   PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));
2188:   PetscFunctionReturn(PETSC_SUCCESS);
2189: }