Actual source code: gmshlex.h
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: */
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)
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)
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)
39: #define GMSH_MAX_ORDER 10
41: static inline int GmshLexOrder_VTX(int p, int lex[], int node)
42: {
43: lex[0] = node++;
44: (void)p;
45: return node;
46: }
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: }
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: }
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: }
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: }
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: }
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: }
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: }