Actual source code: plexgmsh.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    167: static GmshCellInfo GmshCellMap[150];

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    392: #define GMSH_MAX_TAGS 16

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    937:   PetscFunctionBegin;
    938:   PetscCall(GmshReadPetscCount(gmsh, count, 4));
    939:   PetscCall(GmshEntitiesCreate(count, &mesh->entities));
    940:   for (dim = 0; dim < 4; ++dim) {
    941:     for (index = 0; index < count[dim]; ++index) {
    942:       PetscCall(GmshReadInt(gmsh, &eid, 1));
    943:       PetscCall(GmshEntitiesAdd(mesh->entities, (PetscInt)index, dim, eid, &entity));
    944:       PetscCall(GmshReadDouble(gmsh, entity->bbox, (dim == 0) ? 3 : 6));
    945:       PetscCall(GmshReadPetscCount(gmsh, &numTags, 1));
    946:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
    947:       PetscCall(GmshReadInt(gmsh, tags, numTags));
    948:       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);
    949:       PetscCall(PetscIntCast(numTags, &entity->numTags));
    950:       for (PetscInt i = 0; i < entity->numTags; ++i) entity->tags[i] = tags[i];
    951:       if (dim == 0) continue;
    952:       PetscCall(GmshReadPetscCount(gmsh, &numTags, 1));
    953:       PetscCall(GmshBufferGet(gmsh, numTags, sizeof(int), &tags));
    954:       PetscCall(GmshReadInt(gmsh, tags, numTags));
    955:       /* Currently, we do not save the ids for the bounding entities */
    956:     }
    957:   }
    958:   PetscFunctionReturn(PETSC_SUCCESS);
    959: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1485:   Level: beginner

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

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

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

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

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

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

    1540:   Collective

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

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

    1550:   Options Database Keys:
    1551: + -dm_plex_gmsh_hybrid        - Force triangular prisms to use tensor order
    1552: . -dm_plex_gmsh_periodic      - Read Gmsh periodic section and construct a periodic Plex
    1553: . -dm_plex_gmsh_highorder     - Generate high-order coordinates
    1554: . -dm_plex_gmsh_project       - Project high-order coordinates to a different space, use the prefix dm_plex_gmsh_project_ to define the space
    1555: . -dm_plex_gmsh_use_generic   - Generate generic labels, i.e. Cell Sets, Face Sets, etc.
    1556: . -dm_plex_gmsh_use_regions   - Generate labels with region names
    1557: . -dm_plex_gmsh_mark_vertices - Add vertices 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, 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_multiple_tags", "Allow multiple tags for default labels", "DMPlexCreateGmsh", multipleTags, &multipleTags, NULL));
    1602:   PetscCall(PetscOptionsBoundedInt("-dm_plex_gmsh_spacedim", "Embedding space dimension", "DMPlexCreateGmsh", coordDim, &coordDim, NULL, PETSC_DECIDE));
    1603:   PetscCall(PetscOptionsBoundedInt("-dm_localize_height", "Localize edges and faces in addition to cells", "", maxHeight, &maxHeight, NULL, 0));
    1604:   PetscOptionsHeadEnd();
    1605:   PetscOptionsEnd();

    1607:   PetscCall(GmshCellInfoSetUp());

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

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

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

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

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

    1630:     PetscCall(GmshMeshCreate(&mesh));

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

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

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

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

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

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

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

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

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

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

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

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

    1727:   {
    1728:     int buf[6];

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1892:         if (vv < 0) continue;
    1893:         for (t = 0; t < GMSH_MAX_TAGS; ++t) {
    1894:           const PetscInt  tag     = tags[t];
    1895:           const PetscBool generic = usegeneric && (!t || multipleTags) ? PETSC_TRUE : PETSC_FALSE;

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

    1908:   { /* Create Cell/Face/Edge/Vertex Sets labels at all processes */
    1909:     enum {
    1910:       n = 5
    1911:     };
    1912:     PetscBool flag[n];

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

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

    1943:       PetscCall(DMPlexGetHeightStratum(*dm, h, &pStart, &pEnd));
    1944:       PetscCall(PetscBTCreate(pEnd - pStart, &periodicCells[h]));
    1945:       for (PetscInt p = pStart; p < pEnd; ++p) {
    1946:         PetscInt *closure = NULL;
    1947:         PetscInt  Ncl;

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

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

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

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

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

    1992:     PetscCall(DMSetLocalSection(cdm, NULL));
    1993:     PetscCall(DMGetLocalSection(cdm, &cs));
    1994:     PetscCall(PetscSectionClone(cs, &section));
    1995:     PetscCall(DMPlexSetClosurePermutationTensor(cdm, 0, section)); /* XXX Implement DMPlexSetClosurePermutationLexicographic() */

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

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

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

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

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

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

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

    2097:       PetscCall(PetscSectionGetOffset(cs, numCells + v, &off));
    2098:       for (d = 0; d < coordDim; ++d) pointCoords[off + d] = (PetscReal)coords[node * 3 + d];
    2099:     }
    2100:     PetscCall(VecRestoreArray(coordinates, &pointCoords));

    2102:     if (periodic) {
    2103:       PetscInt cStart, cEnd;

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

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

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

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

    2158:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
    2159:   PetscCall(VecSetBlockSize(coordinates, coordDim));
    2160:   PetscCall(DMSetCoordinatesLocal(*dm, coordinates));
    2161:   PetscCall(VecDestroy(&coordinates));

    2163:   PetscCall(GmshMeshDestroy(&mesh));
    2164:   PetscCall(PetscBTDestroy(&periodicVerts));
    2165:   if (periodic) {
    2166:     for (PetscInt h = 0; h <= maxHeight; ++h) PetscCall(PetscBTDestroy(periodicCells + h));
    2167:     PetscCall(PetscFree(periodicCells));
    2168:   }

    2170:   if (highOrder && project) {
    2171:     PetscFE         fe;
    2172:     const char      prefix[]   = "dm_plex_gmsh_project_";
    2173:     PetscBool       continuity = periodic ? PETSC_FALSE : PETSC_TRUE;
    2174:     PetscDTNodeType nodeType   = PETSCDTNODES_GAUSSJACOBI;

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

    2183:   PetscCall(PetscLogEventEnd(DMPLEX_CreateGmsh, *dm, NULL, NULL, NULL));
    2184:   PetscFunctionReturn(PETSC_SUCCESS);
    2185: }