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