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: }