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