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 */ 9 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) 19 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) 29 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) 38 39 #define GMSH_MAX_ORDER 10 40 41 static inline int GmshLexOrder_VTX(int p, int lex[], int node) 42 { 43 lex[0] = node++; 44 (void)p; 45 return node; 46 } 47 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 } 66 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 } 94 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 } 123 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 } 167 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 } 223 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 } 298 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 } 305