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