1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include "pointsGenerators.h"
7 #include "GmshDefines.h"
8 #include "MTriangle.h"
9 #include "MQuadrangle.h"
10 #include "MTetrahedron.h"
11 #include "MHexahedron.h"
12 #include "MPrism.h"
13 #include "MPyramid.h"
14 #include "FuncSpaceData.h"
15 
16 // Points Generators
17 
gmshGeneratePoints(FuncSpaceData data,fullMatrix<double> & points)18 void gmshGeneratePoints(FuncSpaceData data, fullMatrix<double> &points)
19 {
20   const int order = data.getSpaceOrder();
21   const bool serendip = data.getSerendipity();
22   switch(data.getType()) {
23   case TYPE_PNT: points = gmshGeneratePointsLine(0); return;
24   case TYPE_LIN: points = gmshGeneratePointsLine(order); return;
25   case TYPE_TRI: points = gmshGeneratePointsTriangle(order, serendip); return;
26   case TYPE_QUA: points = gmshGeneratePointsQuadrangle(order, serendip); return;
27   case TYPE_TET:
28     points = gmshGeneratePointsTetrahedron(order, serendip);
29     return;
30   case TYPE_PRI: points = gmshGeneratePointsPrism(order, serendip); return;
31   case TYPE_HEX: points = gmshGeneratePointsHexahedron(order, serendip); return;
32   case TYPE_PYR:
33     points = gmshGeneratePointsPyramidGeneral(
34       data.getPyramidalSpace(), data.getNij(), data.getNk(), serendip);
35     return;
36   default:
37     Msg::Error("Unknown element type %d for points generation", data.getType());
38     return;
39   }
40 }
41 
gmshGeneratePointsLine(int order)42 fullMatrix<double> gmshGeneratePointsLine(int order)
43 {
44   fullMatrix<double> points = gmshGenerateMonomialsLine(order);
45   if(order == 0) return points;
46   points.scale(2. / order);
47   points.add(-1.);
48   return points;
49 }
50 
gmshGeneratePointsTriangle(int order,bool serendip)51 fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
52 {
53   fullMatrix<double> points = gmshGenerateMonomialsTriangle(order, serendip);
54   if(order == 0) return points;
55   points.scale(1. / order);
56   return points;
57 }
58 
gmshGeneratePointsQuadrangle(int order,bool serendip)59 fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip)
60 {
61   fullMatrix<double> points = gmshGenerateMonomialsQuadrangle(order, serendip);
62   if(order == 0) return points;
63   points.scale(2. / order);
64   points.add(-1.);
65   return points;
66 }
67 
gmshGeneratePointsTetrahedron(int order,bool serendip)68 fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
69 {
70   fullMatrix<double> points = gmshGenerateMonomialsTetrahedron(order, serendip);
71   if(order == 0) return points;
72   points.scale(1. / order);
73   return points;
74 }
75 
gmshGeneratePointsHexahedron(int order,bool serendip)76 fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip)
77 {
78   fullMatrix<double> points = gmshGenerateMonomialsHexahedron(order, serendip);
79   if(order == 0) return points;
80   points.scale(2. / order);
81   points.add(-1.);
82   return points;
83 }
84 
gmshGeneratePointsPrism(int order,bool serendip)85 fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
86 {
87   fullMatrix<double> points = gmshGenerateMonomialsPrism(order, serendip);
88   if(order == 0) return points;
89   fullMatrix<double> tmp;
90   tmp.setAsProxy(points, 0, 2);
91   tmp.scale(1. / order);
92   tmp.setAsProxy(points, 2, 1);
93   tmp.scale(2. / order);
94   tmp.add(-1.);
95   return points;
96 }
97 
gmshGeneratePointsPyramid(int order,bool serendip)98 fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
99 {
100   fullMatrix<double> points = gmshGenerateMonomialsPyramid(order, serendip);
101   if(order == 0) return points;
102   for(int i = 0; i < points.size1(); ++i) {
103     points(i, 2) = 1 - points(i, 2) / order;
104     const double duv = -1. + points(i, 2);
105     points(i, 0) = duv + points(i, 0) * 2. / order;
106     points(i, 1) = duv + points(i, 1) * 2. / order;
107   }
108   return points;
109 }
110 
gmshGeneratePointsPyramidGeneral(bool pyr,int nij,int nk,bool forSerendipPoints)111 fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk,
112                                                     bool forSerendipPoints)
113 {
114   // if pyr:
115   //   div = nk + nij
116   //   monomial(i, j, k) -> (-(1-k')+2*i/div, -(1-k')+2*j/div, (nk-k)/div)
117   // else:
118   //   div = max(nij, nk)
119   //   monomial(i, j, k) -> (-1+2*i/nij)*(1-k'), (-1+2*j/nij)*(1-k'),
120   //   (nk-k)/div)
121   fullMatrix<double> points =
122     gmshGenerateMonomialsPyramidGeneral(pyr, nij, nk, forSerendipPoints);
123   if(points.size1() == 1) return points;
124   const int div = pyr ? nk + nij : std::max(nij, nk);
125   double scale = 2. / div;
126   for(int i = 0; i < points.size1(); ++i) {
127     points(i, 2) = (nk - points(i, 2)) / div;
128     const double duv = 1. - points(i, 2);
129     if(!pyr) scale = 2. / nij * duv;
130     points(i, 0) = -duv + points(i, 0) * scale;
131     points(i, 1) = -duv + points(i, 1) * scale;
132   }
133   return points;
134 }
135 
136 // Monomials Generators
137 
gmshGenerateMonomials(FuncSpaceData data,fullMatrix<double> & monomials)138 void gmshGenerateMonomials(FuncSpaceData data, fullMatrix<double> &monomials)
139 {
140   const int order = data.getSpaceOrder();
141   const bool serendip = data.getSerendipity();
142   switch(data.getType()) {
143   case TYPE_PNT: monomials = gmshGenerateMonomialsLine(0); return;
144   case TYPE_LIN: monomials = gmshGenerateMonomialsLine(order); return;
145   case TYPE_TRI:
146     monomials = gmshGenerateMonomialsTriangle(order, serendip);
147     return;
148   case TYPE_QUA:
149     monomials = gmshGenerateMonomialsQuadrangle(order, serendip);
150     return;
151   case TYPE_TET:
152     monomials = gmshGenerateMonomialsTetrahedron(order, serendip);
153     return;
154   case TYPE_PRI:
155     monomials = gmshGenerateMonomialsPrism(order, serendip);
156     return;
157   case TYPE_HEX:
158     monomials = gmshGenerateMonomialsHexahedron(order, serendip);
159     return;
160   case TYPE_PYR:
161     monomials = gmshGenerateMonomialsPyramidGeneral(
162       data.getPyramidalSpace(), data.getNij(), data.getNk(), serendip);
163     return;
164   default:
165     Msg::Error("Unknown element type %d for monomials generation",
166                data.getType());
167     return;
168   }
169 }
170 
gmshGenerateMonomialsLine(int order,bool serendip)171 fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip)
172 {
173   fullMatrix<double> monomials(order + 1, 1);
174   monomials(0, 0) = 0;
175   if(order > 0) {
176     monomials(1, 0) = order;
177     for(int i = 2; i < order + 1; i++) monomials(i, 0) = i - 1;
178   }
179   return monomials;
180 }
181 
gmshGenerateMonomialsTriangle(int order,bool serendip)182 fullMatrix<double> gmshGenerateMonomialsTriangle(int order, bool serendip)
183 {
184   int nbMonomials = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
185   if(serendip && !order) nbMonomials = 1;
186   fullMatrix<double> monomials(nbMonomials, 2);
187 
188   monomials(0, 0) = 0;
189   monomials(0, 1) = 0;
190 
191   if(order > 0) {
192     monomials(1, 0) = order;
193     monomials(1, 1) = 0;
194 
195     monomials(2, 0) = 0;
196     monomials(2, 1) = order;
197 
198     if(order > 1) {
199       int index = 3;
200       for(int iedge = 0; iedge < 3; ++iedge) {
201         int i0 = MTriangle::edges_tri(iedge, 0);
202         int i1 = MTriangle::edges_tri(iedge, 1);
203 
204         int u_0 = (monomials(i1, 0) - monomials(i0, 0)) / order;
205         int u_1 = (monomials(i1, 1) - monomials(i0, 1)) / order;
206 
207         for(int i = 1; i < order; ++i, ++index) {
208           monomials(index, 0) = monomials(i0, 0) + u_0 * i;
209           monomials(index, 1) = monomials(i0, 1) + u_1 * i;
210         }
211       }
212       if(!serendip && order > 2) {
213         fullMatrix<double> inner = gmshGenerateMonomialsTriangle(order - 3);
214         inner.add(1);
215         monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
216       }
217     }
218   }
219   return monomials;
220 }
221 
gmshGenerateMonomialsQuadrangle(int order,bool forSerendipPoints)222 fullMatrix<double> gmshGenerateMonomialsQuadrangle(int order,
223                                                    bool forSerendipPoints)
224 {
225   int nbMonomials = forSerendipPoints ? order * 4 : (order + 1) * (order + 1);
226   if(forSerendipPoints && !order) nbMonomials = 1;
227   fullMatrix<double> monomials(nbMonomials, 2);
228 
229   monomials(0, 0) = 0;
230   monomials(0, 1) = 0;
231 
232   if(order > 0) {
233     monomials(1, 0) = order;
234     monomials(1, 1) = 0;
235 
236     monomials(2, 0) = order;
237     monomials(2, 1) = order;
238 
239     monomials(3, 0) = 0;
240     monomials(3, 1) = order;
241 
242     if(order > 1) {
243       int index = 4;
244       for(int iedge = 0; iedge < 4; ++iedge) {
245         int i0 = MQuadrangle::edges_quad(iedge, 0);
246         int i1 = MQuadrangle::edges_quad(iedge, 1);
247 
248         int u_0 = (monomials(i1, 0) - monomials(i0, 0)) / order;
249         int u_1 = (monomials(i1, 1) - monomials(i0, 1)) / order;
250 
251         for(int i = 1; i < order; ++i, ++index) {
252           monomials(index, 0) = monomials(i0, 0) + u_0 * i;
253           monomials(index, 1) = monomials(i0, 1) + u_1 * i;
254         }
255       }
256 
257       if(!forSerendipPoints) {
258         fullMatrix<double> inner = gmshGenerateMonomialsQuadrangle(order - 2);
259         inner.add(1);
260         monomials.copy(inner, 0, nbMonomials - index, 0, 2, index, 0);
261       }
262     }
263   }
264   return monomials;
265 }
266 
267 /*
268 00 10 20 30 40 ...
269 01 11 21 31 41 ...
270 02 12
271 03 13
272 04 14
273 ...
274 */
275 
gmshGenerateMonomialsQuadSerendipity(int order)276 fullMatrix<double> gmshGenerateMonomialsQuadSerendipity(int order)
277 {
278   int nbMonomials = order ? order * 4 : 1;
279   fullMatrix<double> monomials(nbMonomials, 2);
280 
281   monomials(0, 0) = 0;
282   monomials(0, 1) = 0;
283 
284   if(order > 0) {
285     monomials(1, 0) = 1;
286     monomials(1, 1) = 0;
287 
288     monomials(2, 0) = 1;
289     monomials(2, 1) = 1;
290 
291     monomials(3, 0) = 0;
292     monomials(3, 1) = 1;
293 
294     if(order > 1) {
295       int index = 3;
296       for(int p = 2; p <= order; ++p) {
297         monomials(++index, 0) = p;
298         monomials(index, 1) = 0;
299 
300         monomials(++index, 0) = p;
301         monomials(index, 1) = 1;
302 
303         monomials(++index, 0) = 1;
304         monomials(index, 1) = p;
305 
306         monomials(++index, 0) = 0;
307         monomials(index, 1) = p;
308       }
309     }
310   }
311   return monomials;
312 }
313 
314 // KH : caveat : node coordinates are not yet coherent with node numbering
315 // associated
316 //              to numbering of principal vertices of face !!!!
gmshGenerateMonomialsTetrahedron(int order,bool serendip)317 fullMatrix<double> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
318 {
319   int nbMonomials = serendip ? 4 + (order - 1) * 6 :
320                                (order + 1) * (order + 2) * (order + 3) / 6;
321   if(serendip && !order) nbMonomials = 1;
322   fullMatrix<double> monomials(nbMonomials, 3);
323 
324   monomials(0, 0) = 0;
325   monomials(0, 1) = 0;
326   monomials(0, 2) = 0;
327 
328   if(order > 0) {
329     monomials(1, 0) = order;
330     monomials(1, 1) = 0;
331     monomials(1, 2) = 0;
332 
333     monomials(2, 0) = 0;
334     monomials(2, 1) = order;
335     monomials(2, 2) = 0;
336 
337     monomials(3, 0) = 0;
338     monomials(3, 1) = 0;
339     monomials(3, 2) = order;
340 
341     // the template has been defined in table edges_tetra and faces_tetra
342     // (MElement.h)
343 
344     if(order > 1) {
345       int index = 4;
346       for(int iedge = 0; iedge < 6; ++iedge) {
347         int i0 = MTetrahedron::edges_tetra(iedge, 0);
348         int i1 = MTetrahedron::edges_tetra(iedge, 1);
349 
350         int u[3];
351         u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
352         u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
353         u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
354 
355         for(int i = 1; i < order; ++i, ++index) {
356           monomials(index, 0) = monomials(i0, 0) + u[0] * i;
357           monomials(index, 1) = monomials(i0, 1) + u[1] * i;
358           monomials(index, 2) = monomials(i0, 2) + u[2] * i;
359         }
360       }
361 
362       if(!serendip && order > 2) {
363         fullMatrix<double> dudv = gmshGenerateMonomialsTriangle(order - 3);
364         dudv.add(1);
365 
366         for(int iface = 0; iface < 4; ++iface) {
367           int i0 = MTetrahedron::faces_tetra(iface, 0);
368           int i1 = MTetrahedron::faces_tetra(iface, 1);
369           int i2 = MTetrahedron::faces_tetra(iface, 2);
370 
371           int u[3];
372           u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
373           u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
374           u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
375           int v[3];
376           v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
377           v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
378           v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
379 
380           for(int i = 0; i < dudv.size1(); ++i, ++index) {
381             monomials(index, 0) =
382               monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
383             monomials(index, 1) =
384               monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
385             monomials(index, 2) =
386               monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
387           }
388         }
389 
390         if(order > 3) {
391           fullMatrix<double> inner =
392             gmshGenerateMonomialsTetrahedron(order - 4);
393           inner.add(1);
394           monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
395         }
396       }
397     }
398   }
399   return monomials;
400 }
401 
gmshGenerateMonomialsPrism(int order,bool forSerendipPoints)402 fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
403 {
404   int nbMonomials = forSerendipPoints ?
405                       6 + (order - 1) * 9 :
406                       (order + 1) * (order + 1) * (order + 2) / 2;
407   if(forSerendipPoints && !order) nbMonomials = 1;
408   fullMatrix<double> monomials(nbMonomials, 3);
409 
410   monomials(0, 0) = 0;
411   monomials(0, 1) = 0;
412   monomials(0, 2) = 0;
413 
414   if(order > 0) {
415     monomials(1, 0) = order;
416     monomials(1, 1) = 0;
417     monomials(1, 2) = 0;
418 
419     monomials(2, 0) = 0;
420     monomials(2, 1) = order;
421     monomials(2, 2) = 0;
422 
423     monomials(3, 0) = 0;
424     monomials(3, 1) = 0;
425     monomials(3, 2) = order;
426 
427     monomials(4, 0) = order;
428     monomials(4, 1) = 0;
429     monomials(4, 2) = order;
430 
431     monomials(5, 0) = 0;
432     monomials(5, 1) = order;
433     monomials(5, 2) = order;
434 
435     if(order > 1) {
436       int index = 6;
437       for(int iedge = 0; iedge < 9; ++iedge) {
438         int i0 = MPrism::edges_prism(iedge, 0);
439         int i1 = MPrism::edges_prism(iedge, 1);
440 
441         int u_1 = (monomials(i1, 0) - monomials(i0, 0)) / order;
442         int u_2 = (monomials(i1, 1) - monomials(i0, 1)) / order;
443         int u_3 = (monomials(i1, 2) - monomials(i0, 2)) / order;
444 
445         for(int i = 1; i < order; ++i, ++index) {
446           monomials(index, 0) = monomials(i0, 0) + i * u_1;
447           monomials(index, 1) = monomials(i0, 1) + i * u_2;
448           monomials(index, 2) = monomials(i0, 2) + i * u_3;
449         }
450       }
451 
452       if(!forSerendipPoints) {
453         fullMatrix<double> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
454         dudvQ.add(1);
455 
456         fullMatrix<double> dudvT;
457         if(order > 2) {
458           dudvT = gmshGenerateMonomialsTriangle(order - 3);
459           dudvT.add(1);
460         }
461 
462         for(int iface = 0; iface < 5; ++iface) {
463           int i0, i1, i2;
464           i0 = MPrism::faces_prism(iface, 0);
465           i1 = MPrism::faces_prism(iface, 1);
466           fullMatrix<double> dudv;
467           if(MPrism::faces_prism(iface, 3) != -1) {
468             i2 = MPrism::faces_prism(iface, 3);
469             dudv.setAsProxy(dudvQ);
470           }
471           else if(order > 2) {
472             i2 = MPrism::faces_prism(iface, 2);
473             dudv.setAsProxy(dudvT);
474           }
475           else
476             continue;
477 
478           int u[3];
479           u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
480           u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
481           u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
482           int v[3];
483           v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
484           v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
485           v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
486 
487           for(int i = 0; i < dudv.size1(); ++i, ++index) {
488             monomials(index, 0) =
489               monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
490             monomials(index, 1) =
491               monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
492             monomials(index, 2) =
493               monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
494           }
495         }
496 
497         if(order > 2) {
498           fullMatrix<double> triMonomials =
499             gmshGenerateMonomialsTriangle(order - 3);
500           fullMatrix<double> lineMonomials =
501             gmshGenerateMonomialsLine(order - 2);
502 
503           for(int i = 0; i < triMonomials.size1(); ++i) {
504             for(int j = 0; j < lineMonomials.size1(); ++j, ++index) {
505               monomials(index, 0) = 1 + triMonomials(i, 0);
506               monomials(index, 1) = 1 + triMonomials(i, 1);
507               monomials(index, 2) = 1 + lineMonomials(j, 0);
508             }
509           }
510         }
511       }
512     }
513   }
514   return monomials;
515 }
516 
gmshGenerateMonomialsPrismSerendipity(int order)517 fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order)
518 {
519   int nbMonomials = order ? 6 + (order - 1) * 9 : 1;
520   fullMatrix<double> monomials(nbMonomials, 3);
521 
522   monomials(0, 0) = 0;
523   monomials(0, 1) = 0;
524   monomials(0, 2) = 0;
525 
526   if(order > 0) {
527     monomials(1, 0) = 1;
528     monomials(1, 1) = 0;
529     monomials(1, 2) = 0;
530 
531     monomials(2, 0) = 0;
532     monomials(2, 1) = 1;
533     monomials(2, 2) = 0;
534 
535     monomials(3, 0) = 0;
536     monomials(3, 1) = 0;
537     monomials(3, 2) = 1;
538 
539     monomials(4, 0) = 1;
540     monomials(4, 1) = 0;
541     monomials(4, 2) = 1;
542 
543     monomials(5, 0) = 0;
544     monomials(5, 1) = 1;
545     monomials(5, 2) = 1;
546 
547     if(order > 1) {
548       const int ind[7][3] = {{2, 0, 0}, {2, 0, 1},
549 
550                              {0, 2, 0}, {0, 2, 1},
551 
552                              {0, 0, 2}, {1, 0, 2}, {0, 1, 2}};
553       int val[3] = {0, 1, -1};
554       int index = 5;
555       for(int p = 2; p <= order; ++p) {
556         val[2] = p;
557         for(int i = 0; i < 7; ++i) {
558           monomials(++index, 0) = val[ind[i][0]];
559           monomials(index, 1) = val[ind[i][1]];
560           monomials(index, 2) = val[ind[i][2]];
561         }
562       }
563 
564       int val0 = 1;
565       int val1 = order - 1;
566       for(int p = 2; p <= order; ++p) {
567         monomials(++index, 0) = val0;
568         monomials(index, 1) = val1;
569         monomials(index, 2) = 0;
570 
571         monomials(++index, 0) = val0;
572         monomials(index, 1) = val1;
573         monomials(index, 2) = 1;
574 
575         ++val0;
576         --val1;
577       }
578     }
579   }
580   return monomials;
581 }
582 
gmshGenerateMonomialsHexahedron(int order,bool forSerendipPoints)583 fullMatrix<double> gmshGenerateMonomialsHexahedron(int order,
584                                                    bool forSerendipPoints)
585 {
586   int nbMonomials = forSerendipPoints ? 8 + (order - 1) * 12 :
587                                         (order + 1) * (order + 1) * (order + 1);
588   if(forSerendipPoints && !order) nbMonomials = 1;
589   fullMatrix<double> monomials(nbMonomials, 3);
590 
591   monomials(0, 0) = 0;
592   monomials(0, 1) = 0;
593   monomials(0, 2) = 0;
594 
595   if(order > 0) {
596     monomials(1, 0) = order;
597     monomials(1, 1) = 0;
598     monomials(1, 2) = 0;
599 
600     monomials(2, 0) = order;
601     monomials(2, 1) = order;
602     monomials(2, 2) = 0;
603 
604     monomials(3, 0) = 0;
605     monomials(3, 1) = order;
606     monomials(3, 2) = 0;
607 
608     monomials(4, 0) = 0;
609     monomials(4, 1) = 0;
610     monomials(4, 2) = order;
611 
612     monomials(5, 0) = order;
613     monomials(5, 1) = 0;
614     monomials(5, 2) = order;
615 
616     monomials(6, 0) = order;
617     monomials(6, 1) = order;
618     monomials(6, 2) = order;
619 
620     monomials(7, 0) = 0;
621     monomials(7, 1) = order;
622     monomials(7, 2) = order;
623 
624     if(order > 1) {
625       int index = 8;
626       for(int iedge = 0; iedge < 12; ++iedge) {
627         int i0 = MHexahedron::edges_hexa(iedge, 0);
628         int i1 = MHexahedron::edges_hexa(iedge, 1);
629 
630         int u_1 = (monomials(i1, 0) - monomials(i0, 0)) / order;
631         int u_2 = (monomials(i1, 1) - monomials(i0, 1)) / order;
632         int u_3 = (monomials(i1, 2) - monomials(i0, 2)) / order;
633 
634         for(int i = 1; i < order; ++i, ++index) {
635           monomials(index, 0) = monomials(i0, 0) + i * u_1;
636           monomials(index, 1) = monomials(i0, 1) + i * u_2;
637           monomials(index, 2) = monomials(i0, 2) + i * u_3;
638         }
639       }
640 
641       if(!forSerendipPoints) {
642         fullMatrix<double> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
643         dudv.add(1);
644 
645         for(int iface = 0; iface < 6; ++iface) {
646           int i0 = MHexahedron::faces_hexa(iface, 0);
647           int i1 = MHexahedron::faces_hexa(iface, 1);
648           int i3 = MHexahedron::faces_hexa(iface, 3);
649 
650           int u[3];
651           u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
652           u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
653           u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
654           int v[3];
655           v[0] = (monomials(i3, 0) - monomials(i0, 0)) / order;
656           v[1] = (monomials(i3, 1) - monomials(i0, 1)) / order;
657           v[2] = (monomials(i3, 2) - monomials(i0, 2)) / order;
658 
659           for(int i = 0; i < dudv.size1(); ++i, ++index) {
660             monomials(index, 0) =
661               monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
662             monomials(index, 1) =
663               monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
664             monomials(index, 2) =
665               monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
666           }
667         }
668 
669         fullMatrix<double> inner = gmshGenerateMonomialsHexahedron(order - 2);
670         inner.add(1);
671         monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
672       }
673     }
674   }
675   return monomials;
676 }
677 
gmshGenerateMonomialsHexaSerendipity(int order)678 fullMatrix<double> gmshGenerateMonomialsHexaSerendipity(int order)
679 {
680   int nbMonomials = order ? 8 + (order - 1) * 12 : 1;
681   fullMatrix<double> monomials(nbMonomials, 3);
682 
683   monomials(0, 0) = 0;
684   monomials(0, 1) = 0;
685   monomials(0, 2) = 0;
686 
687   if(order > 0) {
688     monomials(1, 0) = 1;
689     monomials(1, 1) = 0;
690     monomials(1, 2) = 0;
691 
692     monomials(2, 0) = 1;
693     monomials(2, 1) = 1;
694     monomials(2, 2) = 0;
695 
696     monomials(3, 0) = 0;
697     monomials(3, 1) = 1;
698     monomials(3, 2) = 0;
699 
700     monomials(4, 0) = 0;
701     monomials(4, 1) = 0;
702     monomials(4, 2) = 1;
703 
704     monomials(5, 0) = 1;
705     monomials(5, 1) = 0;
706     monomials(5, 2) = 1;
707 
708     monomials(6, 0) = 1;
709     monomials(6, 1) = 1;
710     monomials(6, 2) = 1;
711 
712     monomials(7, 0) = 0;
713     monomials(7, 1) = 1;
714     monomials(7, 2) = 1;
715 
716     if(order > 1) {
717       const int ind[12][3] = {{2, 0, 0}, {2, 0, 1}, {2, 1, 1}, {2, 1, 0},
718 
719                               {0, 2, 0}, {0, 2, 1}, {1, 2, 1}, {1, 2, 0},
720 
721                               {0, 0, 2}, {0, 1, 2}, {1, 1, 2}, {1, 0, 2}};
722       int val[3] = {0, 1, -1};
723       int index = 7;
724       for(int p = 2; p <= order; ++p) {
725         val[2] = p;
726         for(int i = 0; i < 12; ++i) {
727           monomials(++index, 0) = val[ind[i][0]];
728           monomials(index, 1) = val[ind[i][1]];
729           monomials(index, 2) = val[ind[i][2]];
730         }
731       }
732     }
733   }
734   return monomials;
735 }
736 
gmshGenerateMonomialsPyramid(int order,bool forSerendipPoints)737 fullMatrix<double> gmshGenerateMonomialsPyramid(int order,
738                                                 bool forSerendipPoints)
739 {
740   int nbMonomials = forSerendipPoints ? 5 + (order - 1) * 8 :
741                                         (order + 1) * ((order + 1) + 1) *
742                                           (2 * (order + 1) + 1) / 6;
743   if(forSerendipPoints && !order) nbMonomials = 1;
744 
745   fullMatrix<double> monomials(nbMonomials, 3);
746 
747   monomials(0, 0) = 0;
748   monomials(0, 1) = 0;
749   monomials(0, 2) = order;
750 
751   if(order > 0) {
752     monomials(1, 0) = order;
753     monomials(1, 1) = 0;
754     monomials(1, 2) = order;
755 
756     monomials(2, 0) = order;
757     monomials(2, 1) = order;
758     monomials(2, 2) = order;
759 
760     monomials(3, 0) = 0;
761     monomials(3, 1) = order;
762     monomials(3, 2) = order;
763 
764     monomials(4, 0) = 0;
765     monomials(4, 1) = 0;
766     monomials(4, 2) = 0;
767 
768     if(order > 1) {
769       int index = 5;
770       for(int iedge = 0; iedge < 8; ++iedge) {
771         int i0 = MPyramid::edges_pyramid(iedge, 0);
772         int i1 = MPyramid::edges_pyramid(iedge, 1);
773 
774         int u_1 = (monomials(i1, 0) - monomials(i0, 0)) / order;
775         int u_2 = (monomials(i1, 1) - monomials(i0, 1)) / order;
776         int u_3 = (monomials(i1, 2) - monomials(i0, 2)) / order;
777 
778         for(int i = 1; i < order; ++i, ++index) {
779           monomials(index, 0) = monomials(i0, 0) + i * u_1;
780           monomials(index, 1) = monomials(i0, 1) + i * u_2;
781           monomials(index, 2) = monomials(i0, 2) + i * u_3;
782         }
783       }
784 
785       if(!forSerendipPoints) {
786         fullMatrix<double> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
787         dudvQ.add(1);
788 
789         fullMatrix<double> dudvT;
790         if(order > 2) {
791           dudvT = gmshGenerateMonomialsTriangle(order - 3);
792           dudvT.add(1);
793         }
794 
795         for(int iface = 0; iface < 5; ++iface) {
796           int i0, i1, i2;
797           i0 = MPyramid::faces_pyramid(iface, 0);
798           i1 = MPyramid::faces_pyramid(iface, 1);
799           fullMatrix<double> dudv;
800           if(MPyramid::faces_pyramid(iface, 3) != -1) {
801             i2 = MPyramid::faces_pyramid(iface, 3);
802             dudv.setAsProxy(dudvQ);
803           }
804           else if(order > 2) {
805             i2 = MPyramid::faces_pyramid(iface, 2);
806             dudv.setAsProxy(dudvT);
807           }
808           else
809             continue;
810 
811           int u[3];
812           u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
813           u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
814           u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
815           int v[3];
816           v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
817           v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
818           v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
819 
820           for(int i = 0; i < dudv.size1(); ++i, ++index) {
821             monomials(index, 0) =
822               monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
823             monomials(index, 1) =
824               monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
825             monomials(index, 2) =
826               monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
827           }
828         }
829 
830         if(order > 2) {
831           fullMatrix<double> inner = gmshGenerateMonomialsPyramid(order - 3);
832           fullMatrix<double> prox;
833           prox.setAsProxy(inner, 0, 2);
834           prox.add(1);
835           prox.setAsProxy(inner, 2, 1);
836           prox.add(2);
837           monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
838         }
839       }
840     }
841   }
842   return monomials;
843 }
844 
gmshGenerateMonomialsPyramidSerendipity(int order)845 fullMatrix<double> gmshGenerateMonomialsPyramidSerendipity(int order)
846 {
847   // WARNING: Is it correct?
848   int nbMonomials = order ? 5 + (order - 1) * 8 : 1;
849 
850   fullMatrix<double> monomials(nbMonomials, 3);
851 
852   monomials(0, 0) = 0;
853   monomials(0, 1) = 0;
854   monomials(0, 2) = 0;
855 
856   if(order > 0) {
857     monomials(1, 0) = 0;
858     monomials(1, 1) = 0;
859     monomials(1, 2) = 1;
860 
861     monomials(2, 0) = 1;
862     monomials(2, 1) = 0;
863     monomials(2, 2) = 1;
864 
865     monomials(3, 0) = 0;
866     monomials(3, 1) = 1;
867     monomials(3, 2) = 1;
868 
869     monomials(4, 0) = 1;
870     monomials(4, 1) = 1;
871     monomials(4, 2) = 1;
872 
873     // monomials of tetrahedra plus a bit more
874     if(order > 1) {
875       int idx = 5;
876       for(int i = 2; i <= order; ++i, ++idx) {
877         monomials(idx, 0) = 0;
878         monomials(idx, 1) = 0;
879         monomials(idx, 2) = i;
880       }
881       for(int i = 2; i <= order; ++i, ++idx) {
882         monomials(idx, 0) = i;
883         monomials(idx, 1) = 0;
884         monomials(idx, 2) = i;
885       }
886       for(int i = 2; i <= order; ++i, ++idx) {
887         monomials(idx, 0) = 0;
888         monomials(idx, 1) = i;
889         monomials(idx, 2) = i;
890       }
891       for(int i = 1; i < order; ++i, ++idx) {
892         monomials(idx, 0) = i;
893         monomials(idx, 1) = 0;
894         monomials(idx, 2) = order;
895       }
896       for(int i = 1; i < order; ++i, ++idx) {
897         monomials(idx, 0) = 0;
898         monomials(idx, 1) = i;
899         monomials(idx, 2) = order;
900       }
901       for(int i = 1; i < order; ++i, ++idx) {
902         monomials(idx, 0) = i;
903         monomials(idx, 1) = order - i;
904         monomials(idx, 2) = order;
905       }
906       for(int i = 2; i <= order; ++i, ++idx) {
907         monomials(idx, 0) = i;
908         monomials(idx, 1) = 1;
909         monomials(idx, 2) = i;
910       }
911       for(int i = 2; i <= order; ++i, ++idx) {
912         monomials(idx, 0) = 1;
913         monomials(idx, 1) = i;
914         monomials(idx, 2) = i;
915       }
916     }
917   }
918   return monomials;
919 }
920 
gmshGenerateMonomialsPyramidGeneral(bool pyr,int nij,int nk,bool forSerendipPoints)921 fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij,
922                                                        int nk,
923                                                        bool forSerendipPoints)
924 {
925   if(nij < 0 || nk < 0) {
926     Msg::Error("Wrong arguments for pyramid's monomials generation ! (%d & %d)",
927                nij, nk);
928     nij = nk = 1;
929   }
930   if(!pyr && nk > 0 && nij == 0) {
931     Msg::Error("Wrong argument association for pyramid's monomials generation! "
932                "Setting nij to 1");
933     nij = 1;
934   }
935 
936   // If monomials for pyramidal space, generate them in gmsh convention.
937   if(forSerendipPoints || (pyr && nij == 0))
938     return gmshGenerateMonomialsPyramid(nk, forSerendipPoints);
939 
940   // Otherwise, just put corners at first places,
941   // order of others having no importance.
942 
943   int nbMonomials;
944   if(pyr) {
945     nbMonomials = 0;
946     for(int k = 0; k <= nk; ++k) {
947       nbMonomials += (k + nij + 1) * (k + nij + 1);
948     }
949   }
950   else
951     nbMonomials = (nij + 1) * (nij + 1) * (nk + 1);
952 
953   fullMatrix<double> monomials(nbMonomials, 3);
954 
955   monomials(0, 0) = 0;
956   monomials(0, 1) = 0;
957   monomials(0, 2) = nk;
958 
959   if(nk == 0 && nij == 0) return monomials;
960 
961   // Here, nij > 0
962   int nijBase = pyr ? nk + nij : nij;
963 
964   // Corners
965   monomials(1, 0) = nijBase;
966   monomials(1, 1) = 0;
967   monomials(1, 2) = nk;
968 
969   monomials(2, 0) = nijBase;
970   monomials(2, 1) = nijBase;
971   monomials(2, 2) = nk;
972 
973   monomials(3, 0) = 0;
974   monomials(3, 1) = nijBase;
975   monomials(3, 2) = nk;
976 
977   int index = 4;
978 
979   if(nk > 0) {
980     monomials(4, 0) = 0;
981     monomials(4, 1) = 0;
982     monomials(4, 2) = 0;
983 
984     monomials(5, 0) = nij;
985     monomials(5, 1) = 0;
986     monomials(5, 2) = 0;
987 
988     monomials(6, 0) = nij;
989     monomials(6, 1) = nij;
990     monomials(6, 2) = 0;
991 
992     monomials(7, 0) = 0;
993     monomials(7, 1) = nij;
994     monomials(7, 2) = 0;
995 
996     index = 8;
997   }
998 
999   // Base
1000   if(nijBase > 1) {
1001     for(int iedge = 0; iedge < 4; ++iedge) {
1002       int i0 = iedge;
1003       int i1 = (iedge + 1) % 4;
1004 
1005       int u0 =
1006         static_cast<int>((monomials(i1, 0) - monomials(i0, 0)) / nijBase);
1007       int u1 =
1008         static_cast<int>((monomials(i1, 1) - monomials(i0, 1)) / nijBase);
1009 
1010       for(int i = 1; i < nijBase; ++i, ++index) {
1011         monomials(index, 0) = monomials(i0, 0) + i * u0;
1012         monomials(index, 1) = monomials(i0, 1) + i * u1;
1013         monomials(index, 2) = nk;
1014       }
1015     }
1016 
1017     for(int i = 1; i < nijBase; ++i) {
1018       for(int j = 1; j < nijBase; ++j, ++index) {
1019         monomials(index, 0) = i;
1020         monomials(index, 1) = j;
1021         monomials(index, 2) = nk;
1022       }
1023     }
1024   }
1025 
1026   // Above base
1027   if(nk > 0) {
1028     // Top
1029     if(nij > 1) {
1030       for(int iedge = 0; iedge < 4; ++iedge) {
1031         int i0 = 4 + iedge;
1032         int i1 = 4 + (iedge + 1) % 4;
1033 
1034         int u0 = static_cast<int>((monomials(i1, 0) - monomials(i0, 0)) / nij);
1035         int u1 = static_cast<int>((monomials(i1, 1) - monomials(i0, 1)) / nij);
1036 
1037         for(int i = 1; i < nij; ++i, ++index) {
1038           monomials(index, 0) = monomials(i0, 0) + i * u0;
1039           monomials(index, 1) = monomials(i0, 1) + i * u1;
1040           monomials(index, 2) = 0;
1041         }
1042       }
1043 
1044       for(int i = 1; i < nij; ++i) {
1045         for(int j = 1; j < nij; ++j, ++index) {
1046           monomials(index, 0) = i;
1047           monomials(index, 1) = j;
1048           monomials(index, 2) = 0;
1049         }
1050       }
1051     }
1052 
1053     // Between bottom & top
1054     for(int k = nk - 1; k > 0; --k) {
1055       int currentnij = pyr ? k + nij : nij;
1056       for(int i = 0; i <= currentnij; ++i) {
1057         for(int j = 0; j <= currentnij; ++j, ++index) {
1058           monomials(index, 0) = i;
1059           monomials(index, 1) = j;
1060           monomials(index, 2) = k;
1061         }
1062       }
1063     }
1064   }
1065 
1066   return monomials;
1067 }
1068 
1069 // Ordered points (and monomials)
1070 
gmshGenerateOrderedPointsLine(int order,fullVector<double> & points)1071 void gmshGenerateOrderedPointsLine(int order, fullVector<double> &points)
1072 {
1073   points.resize(order + 1);
1074   for(int i = 0; i < order + 1; ++i) {
1075     points(i) = i / static_cast<double>(order);
1076   }
1077   return;
1078 }
1079 
gmshGenerateOrderedPoints(FuncSpaceData data,fullMatrix<double> & points,bool onBezierDomain)1080 void gmshGenerateOrderedPoints(FuncSpaceData data, fullMatrix<double> &points,
1081                                bool onBezierDomain)
1082 {
1083   gmshGenerateOrderedMonomials(data, points);
1084   if(points.size1() == 1) return;
1085 
1086   const int type = data.getType();
1087   const int order = data.getSpaceOrder();
1088   const bool pyr = data.getPyramidalSpace();
1089 
1090   if(onBezierDomain) {
1091     if(type != TYPE_PYR) {
1092       points.scale(1. / order);
1093       return;
1094     }
1095     else {
1096       // This is used for creating the matrices bez2lag and lag2bez and we want
1097       // the points of the "unshrinked" pyramid. In other words, we want them
1098       // mapped from the pyramid to the cube.
1099       // if pyr:
1100       //   div = nk + nij
1101       //   monomial(i, j, k) -> (i/div/(1-k'), j/div/(1-k'), (nk-k)/div)
1102       // else:
1103       //   div = max(nij, nk)
1104       //   monomial(i, j, k) -> (i/nij, j/nij, (nk-k)/div)
1105       const int nij = data.getNij();
1106       const int nk = data.getNk();
1107       const int div = pyr ? nij + nk : std::max(nij, nk);
1108       double scale = 1. / nij;
1109       for(int i = 0; i < points.size1(); ++i) {
1110         points(i, 2) = (nk - points(i, 2)) / div;
1111         if(pyr) scale = 1. / div / (1 - points(i, 2));
1112         points(i, 0) = points(i, 0) * scale;
1113         points(i, 1) = points(i, 1) * scale;
1114       }
1115     }
1116     return;
1117   }
1118 
1119   switch(type) {
1120   case TYPE_PNT: return;
1121   case TYPE_LIN:
1122   case TYPE_QUA:
1123   case TYPE_HEX:
1124     points.scale(2. / order);
1125     points.add(-1.);
1126     return;
1127   case TYPE_TRI:
1128   case TYPE_TET: points.scale(1. / order); return;
1129   case TYPE_PRI: {
1130     fullMatrix<double> tmp;
1131     tmp.setAsProxy(points, 0, 2);
1132     tmp.scale(1. / order);
1133     tmp.setAsProxy(points, 2, 1);
1134     tmp.scale(2. / order);
1135     tmp.add(-1.);
1136     return;
1137   }
1138   case TYPE_PYR: {
1139     // This is used e.g. for creating sampling points. We want the points to be
1140     // inside the reference pyramids.
1141     // if pyr:
1142     //   div = nk + nij
1143     //   monomial(i, j, k) -> (-(1-k')+2*i/div, -(1-k')+2*j/div, (nk-k)/div)
1144     // else:
1145     //   div = max(nij, nk)
1146     //   monomial(i, j, k) -> (-1+2*i/nij)*(1-k'), (-1+2*j/nij)*(1-k'),
1147     //   (nk-k)/div)
1148     const int nij = data.getNij();
1149     const int nk = data.getNk();
1150     const int div = pyr ? nij + nk : std::max(nij, nk);
1151     double scale = 2. / div;
1152     for(int i = 0; i < points.size1(); ++i) {
1153       points(i, 2) = (nk - points(i, 2)) / div;
1154       const double duv = 1. - points(i, 2);
1155       if(!pyr) scale = 2. / nij * duv;
1156       points(i, 0) = -duv + points(i, 0) * scale;
1157       points(i, 1) = -duv + points(i, 1) * scale;
1158     }
1159     return;
1160   }
1161   default: return;
1162   }
1163 }
1164 
gmshGenerateOrderedMonomials(FuncSpaceData data,fullMatrix<double> & monomials)1165 void gmshGenerateOrderedMonomials(FuncSpaceData data,
1166                                   fullMatrix<double> &monomials)
1167 {
1168   if(data.getSerendipity())
1169     Msg::Warning("Ordered monomials for serendipity elements not implemented");
1170 
1171   int idx, order = data.getSpaceOrder();
1172 
1173   switch(data.getType()) {
1174   case TYPE_LIN:
1175     monomials.resize(order + 1, 1);
1176     idx = 0;
1177     for(int i = 0; i < order + 1; ++i) {
1178       monomials(idx, 0) = i;
1179       ++idx;
1180     }
1181     return;
1182   case TYPE_TRI:
1183     monomials.resize((order + 1) * (order + 2) / 2, 2);
1184     idx = 0;
1185     for(int j = 0; j < order + 1; ++j) {
1186       for(int i = 0; i < order + 1 - j; ++i) {
1187         monomials(idx, 0) = i;
1188         monomials(idx, 1) = j;
1189         ++idx;
1190       }
1191     }
1192     return;
1193   case TYPE_QUA:
1194     monomials.resize((order + 1) * (order + 1), 2);
1195     idx = 0;
1196     for(int j = 0; j < order + 1; ++j) {
1197       for(int i = 0; i < order + 1; ++i) {
1198         monomials(idx, 0) = i;
1199         monomials(idx, 1) = j;
1200         ++idx;
1201       }
1202     }
1203     return;
1204   case TYPE_TET:
1205     monomials.resize((order + 1) * (order + 2) * (order + 3) / 6, 3);
1206     idx = 0;
1207     for(int k = 0; k < order + 1; ++k) {
1208       for(int j = 0; j < order + 1 - k; ++j) {
1209         for(int i = 0; i < order + 1 - j - k; ++i) {
1210           monomials(idx, 0) = i;
1211           monomials(idx, 1) = j;
1212           monomials(idx, 2) = k;
1213           ++idx;
1214         }
1215       }
1216     }
1217     return;
1218   case TYPE_PRI:
1219     monomials.resize((order + 1) * (order + 1) * (order + 2) / 2, 3);
1220     idx = 0;
1221     for(int k = 0; k < order + 1; ++k) {
1222       for(int j = 0; j < order + 1; ++j) {
1223         for(int i = 0; i < order + 1 - j; ++i) {
1224           monomials(idx, 0) = i;
1225           monomials(idx, 1) = j;
1226           monomials(idx, 2) = k;
1227           ++idx;
1228         }
1229       }
1230     }
1231     return;
1232   case TYPE_HEX:
1233     monomials.resize((order + 1) * (order + 1) * (order + 1), 3);
1234     idx = 0;
1235     for(int k = 0; k < order + 1; ++k) {
1236       for(int j = 0; j < order + 1; ++j) {
1237         for(int i = 0; i < order + 1; ++i) {
1238           monomials(idx, 0) = i;
1239           monomials(idx, 1) = j;
1240           monomials(idx, 2) = k;
1241           ++idx;
1242         }
1243       }
1244     }
1245     return;
1246   case TYPE_PYR: {
1247     const int nij = data.getNij();
1248     const int nk = data.getNk();
1249     if(data.getPyramidalSpace()) {
1250       int n = nk + nij;
1251       int numMonomials = (n + 1) * (n + 2) * (2 * n + 3) / 6;
1252       numMonomials -= nij * (nij + 1) * (2 * nij + 1) / 6;
1253       monomials.resize(numMonomials, 3);
1254       idx = 0;
1255       for(int k = nk; k >= 0; --k) {
1256         for(int j = 0; j < k + nij + 1; ++j) {
1257           for(int i = 0; i < k + nij + 1; ++i) {
1258             monomials(idx, 0) = i;
1259             monomials(idx, 1) = j;
1260             monomials(idx, 2) = k;
1261             ++idx;
1262           }
1263         }
1264       }
1265     }
1266     else {
1267       monomials.resize((nij + 1) * (nij + 1) * (nk + 1), 3);
1268       idx = 0;
1269       for(int k = nk; k >= 0; --k) {
1270         for(int j = 0; j < nij + 1; ++j) {
1271           for(int i = 0; i < nij + 1; ++i) {
1272             monomials(idx, 0) = i;
1273             monomials(idx, 1) = j;
1274             monomials(idx, 2) = k;
1275             ++idx;
1276           }
1277         }
1278       }
1279     }
1280     return;
1281   }
1282   default:
1283     Msg::Error("Unknown element type for ordered monomials: %d",
1284                data.getType());
1285     monomials.resize(1, 1);
1286     return;
1287   }
1288 }
1289