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
GmshLexOrder_VTX(int p,int lex[],int node)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
GmshLexOrder_SEG(int p,int lex[],int node)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
GmshLexOrder_TRI(int p,int lex[],int node)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
GmshLexOrder_QUA(int p,int lex[],int node)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
GmshLexOrder_TET(int p,int lex[],int node)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
GmshLexOrder_HEX(int p,int lex[],int node)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
GmshLexOrder_PRI(int p,int lex[],int node)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
GmshLexOrder_PYR(int p,int lex[],int node)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