1 /* 2 S: simplex B: box 3 N: size I: index L: loop 4 p: degree (aka order in Gmsh) 5 1,2,3: topological dimension 6 i,j,k: coordinate indices 7 */ 8 9 #define SN1(p) ((p) + 1) 10 #define SN2(p) (SN1(p) * SN1((p) + 1) / 2) 11 #define SN3(p) (SN2(p) * SN1((p) + 2) / 3) 12 #define SI1(p, i) ((i)) 13 #define SI2(p, i, j) ((i) + (SN2(p) - SN2((p) - (j)))) 14 #define SI3(p, i, j, k) (SI2((p) - (k), i, j) + (SN3(p) - SN3((p) - (k)))) 15 #define SL1(p, i) for ((i) = 1; (i) < (p); ++(i)) 16 #define SL2(p, i, j) SL1((p)-1, i) SL1((p) - (i), j) 17 #define SL3(p, i, j, k) SL1((p)-2, i) SL1((p) - (i), j) SL1((p) - (i) - (j), k) 18 19 #define BN1(p) ((p) + 1) 20 #define BN2(p) (BN1(p) * BN1(p)) 21 #define BN3(p) (BN2(p) * BN1(p)) 22 #define BI1(p, i) ((i)) 23 #define BI2(p, i, j) ((i) + (j)*BN1(p)) 24 #define BI3(p, i, j, k) ((i) + BI2(p, j, k) * BN1(p)) 25 #define BL1(p, i) for ((i) = 1; (i) < (p); ++(i)) 26 #define BL2(p, i, j) BL1(p, i) BL1(p, j) 27 #define BL3(p, i, j, k) BL1(p, i) BL1(p, j) BL1(p, k) 28 29 #define GmshNumNodes_VTX(p) (1) 30 #define GmshNumNodes_SEG(p) SN1(p) 31 #define GmshNumNodes_TRI(p) SN2(p) 32 #define GmshNumNodes_QUA(p) BN2(p) 33 #define GmshNumNodes_TET(p) SN3(p) 34 #define GmshNumNodes_HEX(p) BN3(p) 35 #define GmshNumNodes_PRI(p) (SN2(p) * BN1(p)) 36 #define GmshNumNodes_PYR(p) (((p) + 1) * ((p) + 2) * (2 * (p) + 3) / 6) 37 38 #define GMSH_MAX_ORDER 10 39 40 static inline int GmshLexOrder_VTX(int p, int lex[], int node) { 41 lex[0] = node++; 42 (void)p; 43 return node; 44 } 45 46 static inline int GmshLexOrder_SEG(int p, int lex[], int node) { 47 #define loop1(i) SL1(p, i) 48 #define index(i) SI1(p, i) 49 int i; 50 /* trivial case */ 51 if (p == 0) lex[0] = node++; 52 if (p == 0) return node; 53 /* vertex nodes */ 54 lex[index(0)] = node++; 55 lex[index(p)] = node++; 56 if (p == 1) return node; 57 /* internal cell nodes */ 58 loop1(i) lex[index(i)] = node++; 59 return node; 60 #undef loop1 61 #undef index 62 } 63 64 static inline int GmshLexOrder_TRI(int p, int lex[], int node) { 65 #define loop1(i) SL1(p, i) 66 #define loop2(i, j) SL2(p, i, j) 67 #define index(i, j) SI2(p, i, j) 68 int i, j, *sub, buf[SN2(GMSH_MAX_ORDER)]; 69 /* trivial case */ 70 if (p == 0) lex[0] = node++; 71 if (p == 0) return node; 72 /* vertex nodes */ 73 lex[index(0, 0)] = node++; 74 lex[index(p, 0)] = node++; 75 lex[index(0, p)] = node++; 76 if (p == 1) return node; 77 /* internal edge nodes */ 78 loop1(i) lex[index(i, 0)] = node++; 79 loop1(j) lex[index(p - j, j)] = node++; 80 loop1(j) lex[index(0, p - j)] = node++; 81 if (p == 2) return node; 82 /* internal cell nodes */ 83 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 84 loop2(j, i) lex[index(i, j)] = *sub++; 85 return node; 86 #undef loop1 87 #undef loop2 88 #undef index 89 } 90 91 static inline int GmshLexOrder_QUA(int p, int lex[], int node) { 92 #define loop1(i) BL1(p, i) 93 #define loop2(i, j) BL2(p, i, j) 94 #define index(i, j) BI2(p, i, j) 95 int i, j, *sub, buf[BN2(GMSH_MAX_ORDER)]; 96 /* trivial case */ 97 if (p == 0) lex[0] = node++; 98 if (p == 0) return node; 99 /* vertex nodes */ 100 lex[index(0, 0)] = node++; 101 lex[index(p, 0)] = node++; 102 lex[index(p, p)] = node++; 103 lex[index(0, p)] = node++; 104 if (p == 1) return node; 105 /* internal edge nodes */ 106 loop1(i) lex[index(i, 0)] = node++; 107 loop1(j) lex[index(p, j)] = node++; 108 loop1(i) lex[index(p - i, p)] = node++; 109 loop1(j) lex[index(0, p - j)] = node++; 110 /* internal cell nodes */ 111 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 112 loop2(j, i) lex[index(i, j)] = *sub++; 113 return node; 114 #undef loop1 115 #undef loop2 116 #undef index 117 } 118 119 static inline int GmshLexOrder_TET(int p, int lex[], int node) { 120 #define loop1(i) SL1(p, i) 121 #define loop2(i, j) SL2(p, i, j) 122 #define loop3(i, j, k) SL3(p, i, j, k) 123 #define index(i, j, k) SI3(p, i, j, k) 124 int i, j, k, *sub, buf[SN3(GMSH_MAX_ORDER)]; 125 /* trivial case */ 126 if (p == 0) lex[0] = node++; 127 if (p == 0) return node; 128 /* vertex nodes */ 129 lex[index(0, 0, 0)] = node++; 130 lex[index(p, 0, 0)] = node++; 131 lex[index(0, p, 0)] = node++; 132 lex[index(0, 0, p)] = node++; 133 if (p == 1) return node; 134 /* internal edge nodes */ 135 loop1(i) lex[index(i, 0, 0)] = node++; 136 loop1(j) lex[index(p - j, j, 0)] = node++; 137 loop1(j) lex[index(0, p - j, 0)] = node++; 138 loop1(k) lex[index(0, 0, p - k)] = node++; 139 loop1(j) lex[index(0, j, p - j)] = node++; 140 loop1(i) lex[index(i, 0, p - i)] = node++; 141 if (p == 2) return node; 142 /* internal face nodes */ 143 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 144 loop2(i, j) lex[index(i, j, 0)] = *sub++; 145 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 146 loop2(k, i) lex[index(i, 0, k)] = *sub++; 147 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 148 loop2(j, k) lex[index(0, j, k)] = *sub++; 149 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 150 loop2(j, i) lex[index(i, j, p - i - j)] = *sub++; 151 if (p == 3) return node; 152 /* internal cell nodes */ 153 node = GmshLexOrder_TET(p - 4, sub = buf, node); 154 loop3(k, j, i) lex[index(i, j, k)] = *sub++; 155 return node; 156 #undef loop1 157 #undef loop2 158 #undef loop3 159 #undef index 160 } 161 162 static inline int GmshLexOrder_HEX(int p, int lex[], int node) { 163 #define loop1(i) BL1(p, i) 164 #define loop2(i, j) BL2(p, i, j) 165 #define loop3(i, j, k) BL3(p, i, j, k) 166 #define index(i, j, k) BI3(p, i, j, k) 167 int i, j, k, *sub, buf[BN3(GMSH_MAX_ORDER)]; 168 /* trivial case */ 169 if (p == 0) lex[0] = node++; 170 if (p == 0) return node; 171 /* vertex nodes */ 172 lex[index(0, 0, 0)] = node++; 173 lex[index(p, 0, 0)] = node++; 174 lex[index(p, p, 0)] = node++; 175 lex[index(0, p, 0)] = node++; 176 lex[index(0, 0, p)] = node++; 177 lex[index(p, 0, p)] = node++; 178 lex[index(p, p, p)] = node++; 179 lex[index(0, p, p)] = node++; 180 if (p == 1) return node; 181 /* internal edge nodes */ 182 loop1(i) lex[index(i, 0, 0)] = node++; 183 loop1(j) lex[index(0, j, 0)] = node++; 184 loop1(k) lex[index(0, 0, k)] = node++; 185 loop1(j) lex[index(p, j, 0)] = node++; 186 loop1(k) lex[index(p, 0, k)] = node++; 187 loop1(i) lex[index(p - i, p, 0)] = node++; 188 loop1(k) lex[index(p, p, k)] = node++; 189 loop1(k) lex[index(0, p, k)] = node++; 190 loop1(i) lex[index(i, 0, p)] = node++; 191 loop1(j) lex[index(0, j, p)] = node++; 192 loop1(j) lex[index(p, j, p)] = node++; 193 loop1(i) lex[index(p - i, p, p)] = node++; 194 /* internal face nodes */ 195 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 196 loop2(i, j) lex[index(i, j, 0)] = *sub++; 197 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 198 loop2(k, i) lex[index(i, 0, k)] = *sub++; 199 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 200 loop2(j, k) lex[index(0, j, k)] = *sub++; 201 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 202 loop2(k, j) lex[index(p, j, k)] = *sub++; 203 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 204 loop2(k, i) lex[index(p - i, p, k)] = *sub++; 205 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 206 loop2(j, i) lex[index(i, j, p)] = *sub++; 207 /* internal cell nodes */ 208 node = GmshLexOrder_HEX(p - 2, sub = buf, node); 209 loop3(k, j, i) lex[index(i, j, k)] = *sub++; 210 return node; 211 #undef loop1 212 #undef loop2 213 #undef loop3 214 #undef index 215 } 216 217 static inline int GmshLexOrder_PRI(int p, int lex[], int node) { 218 #define loop1(i) BL1(p, i) 219 #define loops(i, j) SL2(p, i, j) 220 #define loopb(i, j) BL2(p, i, j) 221 #define index(i, j, k) (SI2(p, i, j) + BI1(p, k) * SN2(p)) 222 int i, j, k, *sub, buf[BN2(GMSH_MAX_ORDER)]; 223 /* trivial case */ 224 if (p == 0) lex[0] = node++; 225 if (p == 0) return node; 226 /* vertex nodes */ 227 lex[index(0, 0, 0)] = node++; 228 lex[index(p, 0, 0)] = node++; 229 lex[index(0, p, 0)] = node++; 230 lex[index(0, 0, p)] = node++; 231 lex[index(p, 0, p)] = node++; 232 lex[index(0, p, p)] = node++; 233 if (p == 1) return node; 234 /* internal edge nodes */ 235 loop1(i) lex[index(i, 0, 0)] = node++; 236 loop1(j) lex[index(0, j, 0)] = node++; 237 loop1(k) lex[index(0, 0, k)] = node++; 238 loop1(j) lex[index(p - j, j, 0)] = node++; 239 loop1(k) lex[index(p, 0, k)] = node++; 240 loop1(k) lex[index(0, p, k)] = node++; 241 loop1(i) lex[index(i, 0, p)] = node++; 242 loop1(j) lex[index(0, j, p)] = node++; 243 loop1(j) lex[index(p - j, j, p)] = node++; 244 if (p >= 3) { 245 /* internal bottom face nodes */ 246 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 247 loops(i, j) lex[index(i, j, 0)] = *sub++; 248 /* internal top face nodes */ 249 node = GmshLexOrder_TRI(p - 3, sub = buf, node); 250 loops(j, i) lex[index(i, j, p)] = *sub++; 251 } 252 if (p >= 2) { 253 /* internal front face nodes */ 254 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 255 loopb(k, i) lex[index(i, 0, k)] = *sub++; 256 /* internal left face nodes */ 257 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 258 loopb(j, k) lex[index(0, j, k)] = *sub++; 259 /* internal back face nodes */ 260 node = GmshLexOrder_QUA(p - 2, sub = buf, node); 261 loopb(k, j) lex[index(p - j, j, k)] = *sub++; 262 } 263 if (p >= 3) { 264 /* internal cell nodes */ 265 typedef struct { 266 int i, j; 267 } pair; 268 pair ij[SN2(GMSH_MAX_ORDER)], tmp[SN2(GMSH_MAX_ORDER)]; 269 int m = GmshLexOrder_TRI(p - 3, sub = buf, 0), l = 0; 270 loops(j, i) { 271 tmp[l].i = i; 272 tmp[l].j = j; 273 l++; 274 } 275 for (l = 0; l < m; ++l) ij[sub[l]] = tmp[l]; 276 for (l = 0; l < m; ++l) { 277 i = ij[l].i; 278 j = ij[l].j; 279 node = GmshLexOrder_SEG(p - 2, sub = buf, node); 280 loop1(k) lex[index(i, j, k)] = *sub++; 281 } 282 } 283 return node; 284 #undef loop1 285 #undef loops 286 #undef loopb 287 #undef index 288 } 289 290 static inline int GmshLexOrder_PYR(int p, int lex[], int node) { 291 int i, m = GmshNumNodes_PYR(p); 292 for (i = 0; i < m; ++i) { lex[i] = node++; } /* TODO */ 293 return node; 294 } 295