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