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