Actual source code: gmshlex.h

      1: #pragma once
      2: /*
      3:  S: simplex  B: box
      4:  N: size     I: index  L: loop
      5:  p: degree (aka order in Gmsh)
      6:  1,2,3: topological dimension
      7:  i,j,k: coordinate indices
      8: */

     10: #define SN1(p)          ((p) + 1)
     11: #define SN2(p)          (SN1(p) * SN1((p) + 1) / 2)
     12: #define SN3(p)          (SN2(p) * SN1((p) + 2) / 3)
     13: #define SI1(p, i)       (i)
     14: #define SI2(p, i, j)    ((i) + (SN2(p) - SN2((p) - (j))))
     15: #define SI3(p, i, j, k) (SI2((p) - (k), i, j) + (SN3(p) - SN3((p) - (k))))
     16: #define SL1(p, i)       for ((i) = 1; (i) < (p); ++(i))
     17: #define SL2(p, i, j)    SL1((p) - 1, i) SL1((p) - (i), j)
     18: #define SL3(p, i, j, k) SL1((p) - 2, i) SL1((p) - (i), j) SL1((p) - (i) - (j), k)

     20: #define BN1(p)          ((p) + 1)
     21: #define BN2(p)          (BN1(p) * BN1(p))
     22: #define BN3(p)          (BN2(p) * BN1(p))
     23: #define BI1(p, i)       (i)
     24: #define BI2(p, i, j)    ((i) + (j) * BN1(p))
     25: #define BI3(p, i, j, k) ((i) + BI2(p, j, k) * BN1(p))
     26: #define BL1(p, i)       for ((i) = 1; (i) < (p); ++(i))
     27: #define BL2(p, i, j)    BL1(p, i) BL1(p, j)
     28: #define BL3(p, i, j, k) BL1(p, i) BL1(p, j) BL1(p, k)

     30: #define GmshNumNodes_VTX(p) (1)
     31: #define GmshNumNodes_SEG(p) SN1(p)
     32: #define GmshNumNodes_TRI(p) SN2(p)
     33: #define GmshNumNodes_QUA(p) BN2(p)
     34: #define GmshNumNodes_TET(p) SN3(p)
     35: #define GmshNumNodes_HEX(p) BN3(p)
     36: #define GmshNumNodes_PRI(p) (SN2(p) * BN1(p))
     37: #define GmshNumNodes_PYR(p) (((p) + 1) * ((p) + 2) * (2 * (p) + 3) / 6)

     39: #define GMSH_MAX_ORDER 10

     41: static inline int GmshLexOrder_VTX(int p, int lex[], int node)
     42: {
     43:   lex[0] = node++;
     44:   (void)p;
     45:   return node;
     46: }

     48: static inline int GmshLexOrder_SEG(int p, int lex[], int node)
     49: {
     50: #define loop1(i) SL1(p, i)
     51: #define index(i) SI1(p, i)
     52:   int i;
     53:   /* trivial case */
     54:   if (p == 0) lex[0] = node++;
     55:   if (p == 0) return node;
     56:   /* vertex nodes */
     57:   lex[index(0)] = node++;
     58:   lex[index(p)] = node++;
     59:   if (p == 1) return node;
     60:   /* internal cell nodes */
     61:   loop1(i) lex[index(i)] = node++;
     62:   return node;
     63: #undef loop1
     64: #undef index
     65: }

     67: static inline int GmshLexOrder_TRI(int p, int lex[], int node)
     68: {
     69: #define loop1(i)    SL1(p, i)
     70: #define loop2(i, j) SL2(p, i, j)
     71: #define index(i, j) SI2(p, i, j)
     72:   int i, j, *sub, buf[SN2(GMSH_MAX_ORDER)];
     73:   /* trivial case */
     74:   if (p == 0) lex[0] = node++;
     75:   if (p == 0) return node;
     76:   /* vertex nodes */
     77:   lex[index(0, 0)] = node++;
     78:   lex[index(p, 0)] = node++;
     79:   lex[index(0, p)] = node++;
     80:   if (p == 1) return node;
     81:   /* internal edge nodes */
     82:   loop1(i) lex[index(i, 0)]     = node++;
     83:   loop1(j) lex[index(p - j, j)] = node++;
     84:   loop1(j) lex[index(0, p - j)] = node++;
     85:   if (p == 2) return node;
     86:   /* internal cell nodes */
     87:   node                         = GmshLexOrder_TRI(p - 3, sub = buf, node);
     88:   loop2(j, i) lex[index(i, j)] = *sub++;
     89:   return node;
     90: #undef loop1
     91: #undef loop2
     92: #undef index
     93: }

     95: static inline int GmshLexOrder_QUA(int p, int lex[], int node)
     96: {
     97: #define loop1(i)    BL1(p, i)
     98: #define loop2(i, j) BL2(p, i, j)
     99: #define index(i, j) BI2(p, i, j)
    100:   int i, j, *sub, buf[BN2(GMSH_MAX_ORDER)];
    101:   /* trivial case */
    102:   if (p == 0) lex[0] = node++;
    103:   if (p == 0) return node;
    104:   /* vertex nodes */
    105:   lex[index(0, 0)] = node++;
    106:   lex[index(p, 0)] = node++;
    107:   lex[index(p, p)] = node++;
    108:   lex[index(0, p)] = node++;
    109:   if (p == 1) return node;
    110:   /* internal edge nodes */
    111:   loop1(i) lex[index(i, 0)]     = node++;
    112:   loop1(j) lex[index(p, j)]     = node++;
    113:   loop1(i) lex[index(p - i, p)] = node++;
    114:   loop1(j) lex[index(0, p - j)] = node++;
    115:   /* internal cell nodes */
    116:   node                         = GmshLexOrder_QUA(p - 2, sub = buf, node);
    117:   loop2(j, i) lex[index(i, j)] = *sub++;
    118:   return node;
    119: #undef loop1
    120: #undef loop2
    121: #undef index
    122: }

    124: static inline int GmshLexOrder_TET(int p, int lex[], int node)
    125: {
    126: #define loop1(i)       SL1(p, i)
    127: #define loop2(i, j)    SL2(p, i, j)
    128: #define loop3(i, j, k) SL3(p, i, j, k)
    129: #define index(i, j, k) SI3(p, i, j, k)
    130:   int i, j, k, *sub, buf[SN3(GMSH_MAX_ORDER)];
    131:   /* trivial case */
    132:   if (p == 0) lex[0] = node++;
    133:   if (p == 0) return node;
    134:   /* vertex nodes */
    135:   lex[index(0, 0, 0)] = node++;
    136:   lex[index(p, 0, 0)] = node++;
    137:   lex[index(0, p, 0)] = node++;
    138:   lex[index(0, 0, p)] = node++;
    139:   if (p == 1) return node;
    140:   /* internal edge nodes */
    141:   loop1(i) lex[index(i, 0, 0)]     = node++;
    142:   loop1(j) lex[index(p - j, j, 0)] = node++;
    143:   loop1(j) lex[index(0, p - j, 0)] = node++;
    144:   loop1(k) lex[index(0, 0, p - k)] = node++;
    145:   loop1(j) lex[index(0, j, p - j)] = node++;
    146:   loop1(i) lex[index(i, 0, p - i)] = node++;
    147:   if (p == 2) return node;
    148:   /* internal face nodes */
    149:   node                                    = GmshLexOrder_TRI(p - 3, sub = buf, node);
    150:   loop2(i, j) lex[index(i, j, 0)]         = *sub++;
    151:   node                                    = GmshLexOrder_TRI(p - 3, sub = buf, node);
    152:   loop2(k, i) lex[index(i, 0, k)]         = *sub++;
    153:   node                                    = GmshLexOrder_TRI(p - 3, sub = buf, node);
    154:   loop2(j, k) lex[index(0, j, k)]         = *sub++;
    155:   node                                    = GmshLexOrder_TRI(p - 3, sub = buf, node);
    156:   loop2(j, i) lex[index(i, j, p - i - j)] = *sub++;
    157:   if (p == 3) return node;
    158:   /* internal cell nodes */
    159:   node                               = GmshLexOrder_TET(p - 4, sub = buf, node);
    160:   loop3(k, j, i) lex[index(i, j, k)] = *sub++;
    161:   return node;
    162: #undef loop1
    163: #undef loop2
    164: #undef loop3
    165: #undef index
    166: }

    168: static inline int GmshLexOrder_HEX(int p, int lex[], int node)
    169: {
    170: #define loop1(i)       BL1(p, i)
    171: #define loop2(i, j)    BL2(p, i, j)
    172: #define loop3(i, j, k) BL3(p, i, j, k)
    173: #define index(i, j, k) BI3(p, i, j, k)
    174:   int i, j, k, *sub, buf[BN3(GMSH_MAX_ORDER)];
    175:   /* trivial case */
    176:   if (p == 0) lex[0] = node++;
    177:   if (p == 0) return node;
    178:   /* vertex nodes */
    179:   lex[index(0, 0, 0)] = node++;
    180:   lex[index(p, 0, 0)] = node++;
    181:   lex[index(p, p, 0)] = node++;
    182:   lex[index(0, p, 0)] = node++;
    183:   lex[index(0, 0, p)] = node++;
    184:   lex[index(p, 0, p)] = node++;
    185:   lex[index(p, p, p)] = node++;
    186:   lex[index(0, p, p)] = node++;
    187:   if (p == 1) return node;
    188:   /* internal edge nodes */
    189:   loop1(i) lex[index(i, 0, 0)]     = node++;
    190:   loop1(j) lex[index(0, j, 0)]     = node++;
    191:   loop1(k) lex[index(0, 0, k)]     = node++;
    192:   loop1(j) lex[index(p, j, 0)]     = node++;
    193:   loop1(k) lex[index(p, 0, k)]     = node++;
    194:   loop1(i) lex[index(p - i, p, 0)] = node++;
    195:   loop1(k) lex[index(p, p, k)]     = node++;
    196:   loop1(k) lex[index(0, p, k)]     = node++;
    197:   loop1(i) lex[index(i, 0, p)]     = node++;
    198:   loop1(j) lex[index(0, j, p)]     = node++;
    199:   loop1(j) lex[index(p, j, p)]     = node++;
    200:   loop1(i) lex[index(p - i, p, p)] = node++;
    201:   /* internal face nodes */
    202:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    203:   loop2(i, j) lex[index(i, j, 0)]     = *sub++;
    204:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    205:   loop2(k, i) lex[index(i, 0, k)]     = *sub++;
    206:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    207:   loop2(j, k) lex[index(0, j, k)]     = *sub++;
    208:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    209:   loop2(k, j) lex[index(p, j, k)]     = *sub++;
    210:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    211:   loop2(k, i) lex[index(p - i, p, k)] = *sub++;
    212:   node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    213:   loop2(j, i) lex[index(i, j, p)]     = *sub++;
    214:   /* internal cell nodes */
    215:   node                               = GmshLexOrder_HEX(p - 2, sub = buf, node);
    216:   loop3(k, j, i) lex[index(i, j, k)] = *sub++;
    217:   return node;
    218: #undef loop1
    219: #undef loop2
    220: #undef loop3
    221: #undef index
    222: }

    224: static inline int GmshLexOrder_PRI(int p, int lex[], int node)
    225: {
    226: #define loop1(i)       BL1(p, i)
    227: #define loops(i, j)    SL2(p, i, j)
    228: #define loopb(i, j)    BL2(p, i, j)
    229: #define index(i, j, k) (SI2(p, i, j) + BI1(p, k) * SN2(p))
    230:   int i, j, k, *sub, buf[BN2(GMSH_MAX_ORDER)];
    231:   /* trivial case */
    232:   if (p == 0) lex[0] = node++;
    233:   if (p == 0) return node;
    234:   /* vertex nodes */
    235:   lex[index(0, 0, 0)] = node++;
    236:   lex[index(p, 0, 0)] = node++;
    237:   lex[index(0, p, 0)] = node++;
    238:   lex[index(0, 0, p)] = node++;
    239:   lex[index(p, 0, p)] = node++;
    240:   lex[index(0, p, p)] = node++;
    241:   if (p == 1) return node;
    242:   /* internal edge nodes */
    243:   loop1(i) lex[index(i, 0, 0)]     = node++;
    244:   loop1(j) lex[index(0, j, 0)]     = node++;
    245:   loop1(k) lex[index(0, 0, k)]     = node++;
    246:   loop1(j) lex[index(p - j, j, 0)] = node++;
    247:   loop1(k) lex[index(p, 0, k)]     = node++;
    248:   loop1(k) lex[index(0, p, k)]     = node++;
    249:   loop1(i) lex[index(i, 0, p)]     = node++;
    250:   loop1(j) lex[index(0, j, p)]     = node++;
    251:   loop1(j) lex[index(p - j, j, p)] = node++;
    252:   if (p >= 3) {
    253:     /* internal bottom face nodes */
    254:     node                            = GmshLexOrder_TRI(p - 3, sub = buf, node);
    255:     loops(i, j) lex[index(i, j, 0)] = *sub++;
    256:     /* internal top face nodes */
    257:     node                            = GmshLexOrder_TRI(p - 3, sub = buf, node);
    258:     loops(j, i) lex[index(i, j, p)] = *sub++;
    259:   }
    260:   if (p >= 2) {
    261:     /* internal front face nodes */
    262:     node                            = GmshLexOrder_QUA(p - 2, sub = buf, node);
    263:     loopb(k, i) lex[index(i, 0, k)] = *sub++;
    264:     /* internal left face nodes */
    265:     node                            = GmshLexOrder_QUA(p - 2, sub = buf, node);
    266:     loopb(j, k) lex[index(0, j, k)] = *sub++;
    267:     /* internal back face nodes */
    268:     node                                = GmshLexOrder_QUA(p - 2, sub = buf, node);
    269:     loopb(k, j) lex[index(p - j, j, k)] = *sub++;
    270:   }
    271:   if (p >= 3) {
    272:     /* internal cell nodes */
    273:     typedef struct {
    274:       int i, j;
    275:     } pair;
    276:     pair ij[SN2(GMSH_MAX_ORDER)], tmp[SN2(GMSH_MAX_ORDER)];
    277:     int  m = GmshLexOrder_TRI(p - 3, sub = buf, 0), l = 0;
    278:     loops(j, i)
    279:     {
    280:       tmp[l].i = i;
    281:       tmp[l].j = j;
    282:       l++;
    283:     }
    284:     for (l = 0; l < m; ++l) ij[sub[l]] = tmp[l];
    285:     for (l = 0; l < m; ++l) {
    286:       i                            = ij[l].i;
    287:       j                            = ij[l].j;
    288:       node                         = GmshLexOrder_SEG(p - 2, sub = buf, node);
    289:       loop1(k) lex[index(i, j, k)] = *sub++;
    290:     }
    291:   }
    292:   return node;
    293: #undef loop1
    294: #undef loops
    295: #undef loopb
    296: #undef index
    297: }

    299: static inline int GmshLexOrder_PYR(int p, int lex[], int node)
    300: {
    301:   int i, m = GmshNumNodes_PYR(p);
    302:   for (i = 0; i < m; ++i) lex[i] = node++; /* TODO */
    303:   return node;
    304: }