MFEM  v4.6.0 Finite element discretization library
fe_h1.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // H1 Finite Element classes
13
14 #include "fe_h1.hpp"
15
16 namespace mfem
17 {
18
19 using namespace std;
20
21 H1_SegmentElement::H1_SegmentElement(const int p, const int btype)
22  : NodalTensorFiniteElement(1, p, VerifyClosed(btype), H1_DOF_MAP)
23 {
24  const double *cp = poly1d.ClosedPoints(p, b_type);
25
27  shape_x.SetSize(p+1);
28  dshape_x.SetSize(p+1);
29  d2shape_x.SetSize(p+1);
30 #endif
31
32  Nodes.IntPoint(0).x = cp[0];
33  Nodes.IntPoint(1).x = cp[p];
34  for (int i = 1; i < p; i++)
35  {
36  Nodes.IntPoint(i+1).x = cp[i];
37  }
38 }
39
41  Vector &shape) const
42 {
43  const int p = order;
44
46  Vector shape_x(p+1);
47 #endif
48
49  basis1d.Eval(ip.x, shape_x);
50
51  shape(0) = shape_x(0);
52  shape(1) = shape_x(p);
53  for (int i = 1; i < p; i++)
54  {
55  shape(i+1) = shape_x(i);
56  }
57 }
58
60  DenseMatrix &dshape) const
61 {
62  const int p = order;
63
65  Vector shape_x(p+1), dshape_x(p+1);
66 #endif
67
68  basis1d.Eval(ip.x, shape_x, dshape_x);
69
70  dshape(0,0) = dshape_x(0);
71  dshape(1,0) = dshape_x(p);
72  for (int i = 1; i < p; i++)
73  {
74  dshape(i+1,0) = dshape_x(i);
75  }
76 }
77
79  DenseMatrix &Hessian) const
80 {
81  const int p = order;
82
84  Vector shape_x(p+1), dshape_x(p+1), d2shape_x(p+1);
85 #endif
86
87  basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
88
89  Hessian(0,0) = d2shape_x(0);
90  Hessian(1,0) = d2shape_x(p);
91  for (int i = 1; i < p; i++)
92  {
93  Hessian(i+1,0) = d2shape_x(i);
94  }
95 }
96
97 void H1_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
98 {
99  const int p = order;
100  const double *cp = poly1d.ClosedPoints(p, b_type);
101
102  switch (vertex)
103  {
104  case 0:
105  dofs(0) = poly1d.CalcDelta(p, (1.0 - cp[0]));
106  dofs(1) = poly1d.CalcDelta(p, (1.0 - cp[p]));
107  for (int i = 1; i < p; i++)
108  {
109  dofs(i+1) = poly1d.CalcDelta(p, (1.0 - cp[i]));
110  }
111  break;
112
113  case 1:
114  dofs(0) = poly1d.CalcDelta(p, cp[0]);
115  dofs(1) = poly1d.CalcDelta(p, cp[p]);
116  for (int i = 1; i < p; i++)
117  {
118  dofs(i+1) = poly1d.CalcDelta(p, cp[i]);
119  }
120  break;
121  }
122 }
123
124
126  : NodalTensorFiniteElement(2, p, VerifyClosed(btype), H1_DOF_MAP)
127 {
128  const double *cp = poly1d.ClosedPoints(p, b_type);
129
131  const int p1 = p + 1;
132
133  shape_x.SetSize(p1);
134  shape_y.SetSize(p1);
135  dshape_x.SetSize(p1);
136  dshape_y.SetSize(p1);
137  d2shape_x.SetSize(p1);
138  d2shape_y.SetSize(p1);
139 #endif
140
141  int o = 0;
142  for (int j = 0; j <= p; j++)
143  {
144  for (int i = 0; i <= p; i++)
145  {
146  Nodes.IntPoint(dof_map[o++]).Set2(cp[i], cp[j]);
147  }
148  }
149 }
150
152  Vector &shape) const
153 {
154  const int p = order;
155
157  Vector shape_x(p+1), shape_y(p+1);
158 #endif
159
160  basis1d.Eval(ip.x, shape_x);
161  basis1d.Eval(ip.y, shape_y);
162
163  for (int o = 0, j = 0; j <= p; j++)
164  for (int i = 0; i <= p; i++)
165  {
166  shape(dof_map[o++]) = shape_x(i)*shape_y(j);
167  }
168 }
169
171  DenseMatrix &dshape) const
172 {
173  const int p = order;
174
176  Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
177 #endif
178
179  basis1d.Eval(ip.x, shape_x, dshape_x);
180  basis1d.Eval(ip.y, shape_y, dshape_y);
181
182  for (int o = 0, j = 0; j <= p; j++)
183  {
184  for (int i = 0; i <= p; i++)
185  {
186  dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
187  dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
188  }
189  }
190 }
191
193  DenseMatrix &Hessian) const
194 {
195  const int p = order;
196
198  Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1),
199  d2shape_x(p+1), d2shape_y(p+1);
200 #endif
201
202  basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
203  basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
204
205  for (int o = 0, j = 0; j <= p; j++)
206  {
207  for (int i = 0; i <= p; i++)
208  {
209  Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j);
210  Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j);
211  Hessian(dof_map[o],2) = shape_x(i)*d2shape_y(j); o++;
212  }
213  }
214 }
215
216 void H1_QuadrilateralElement::ProjectDelta(int vertex, Vector &dofs) const
217 {
218  const int p = order;
219  const double *cp = poly1d.ClosedPoints(p, b_type);
220
222  Vector shape_x(p+1), shape_y(p+1);
223 #endif
224
225  for (int i = 0; i <= p; i++)
226  {
227  shape_x(i) = poly1d.CalcDelta(p, (1.0 - cp[i]));
228  shape_y(i) = poly1d.CalcDelta(p, cp[i]);
229  }
230
231  switch (vertex)
232  {
233  case 0:
234  for (int o = 0, j = 0; j <= p; j++)
235  for (int i = 0; i <= p; i++)
236  {
237  dofs(dof_map[o++]) = shape_x(i)*shape_x(j);
238  }
239  break;
240  case 1:
241  for (int o = 0, j = 0; j <= p; j++)
242  for (int i = 0; i <= p; i++)
243  {
244  dofs(dof_map[o++]) = shape_y(i)*shape_x(j);
245  }
246  break;
247  case 2:
248  for (int o = 0, j = 0; j <= p; j++)
249  for (int i = 0; i <= p; i++)
250  {
251  dofs(dof_map[o++]) = shape_y(i)*shape_y(j);
252  }
253  break;
254  case 3:
255  for (int o = 0, j = 0; j <= p; j++)
256  for (int i = 0; i <= p; i++)
257  {
258  dofs(dof_map[o++]) = shape_x(i)*shape_y(j);
259  }
260  break;
261  }
262 }
263
264
266  : NodalTensorFiniteElement(3, p, VerifyClosed(btype), H1_DOF_MAP)
267 {
268  const double *cp = poly1d.ClosedPoints(p, b_type);
269
271  const int p1 = p + 1;
272
273  shape_x.SetSize(p1);
274  shape_y.SetSize(p1);
275  shape_z.SetSize(p1);
276  dshape_x.SetSize(p1);
277  dshape_y.SetSize(p1);
278  dshape_z.SetSize(p1);
279  d2shape_x.SetSize(p1);
280  d2shape_y.SetSize(p1);
281  d2shape_z.SetSize(p1);
282 #endif
283
284  int o = 0;
285  for (int k = 0; k <= p; k++)
286  for (int j = 0; j <= p; j++)
287  for (int i = 0; i <= p; i++)
288  {
289  Nodes.IntPoint(dof_map[o++]).Set3(cp[i], cp[j], cp[k]);
290  }
291 }
292
294  Vector &shape) const
295 {
296  const int p = order;
297
299  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
300 #endif
301
302  basis1d.Eval(ip.x, shape_x);
303  basis1d.Eval(ip.y, shape_y);
304  basis1d.Eval(ip.z, shape_z);
305
306  for (int o = 0, k = 0; k <= p; k++)
307  for (int j = 0; j <= p; j++)
308  for (int i = 0; i <= p; i++)
309  {
310  shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
311  }
312 }
313
315  DenseMatrix &dshape) const
316 {
317  const int p = order;
318
320  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
321  Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
322 #endif
323
324  basis1d.Eval(ip.x, shape_x, dshape_x);
325  basis1d.Eval(ip.y, shape_y, dshape_y);
326  basis1d.Eval(ip.z, shape_z, dshape_z);
327
328  for (int o = 0, k = 0; k <= p; k++)
329  for (int j = 0; j <= p; j++)
330  for (int i = 0; i <= p; i++)
331  {
332  dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
333  dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
334  dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
335  }
336 }
337
339  DenseMatrix &Hessian) const
340 {
341  const int p = order;
342
344  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
345  Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
346  Vector d2shape_x(p+1), d2shape_y(p+1), d2shape_z(p+1);
347 #endif
348
349  basis1d.Eval(ip.x, shape_x, dshape_x, d2shape_x);
350  basis1d.Eval(ip.y, shape_y, dshape_y, d2shape_y);
351  basis1d.Eval(ip.z, shape_z, dshape_z, d2shape_z);
352
353  for (int o = 0, k = 0; k <= p; k++)
354  for (int j = 0; j <= p; j++)
355  for (int i = 0; i <= p; i++)
356  {
357  Hessian(dof_map[o],0) = d2shape_x(i)* shape_y(j)* shape_z(k);
358  Hessian(dof_map[o],1) = dshape_x(i)* dshape_y(j)* shape_z(k);
359  Hessian(dof_map[o],2) = dshape_x(i)* shape_y(j)* dshape_z(k);
360  Hessian(dof_map[o],3) = shape_x(i)*d2shape_y(j)* shape_z(k);
361  Hessian(dof_map[o],4) = shape_x(i)* dshape_y(j)* dshape_z(k);
362  Hessian(dof_map[o],5) = shape_x(i)* shape_y(j)*d2shape_z(k);
363  o++;
364  }
365 }
366
367 void H1_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
368 {
369  const int p = order;
370  const double *cp = poly1d.ClosedPoints(p,b_type);
371
373  Vector shape_x(p+1), shape_y(p+1);
374 #endif
375
376  for (int i = 0; i <= p; i++)
377  {
378  shape_x(i) = poly1d.CalcDelta(p, (1.0 - cp[i]));
379  shape_y(i) = poly1d.CalcDelta(p, cp[i]);
380  }
381
382  switch (vertex)
383  {
384  case 0:
385  for (int o = 0, k = 0; k <= p; k++)
386  for (int j = 0; j <= p; j++)
387  for (int i = 0; i <= p; i++)
388  {
389  dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_x(k);
390  }
391  break;
392  case 1:
393  for (int o = 0, k = 0; k <= p; k++)
394  for (int j = 0; j <= p; j++)
395  for (int i = 0; i <= p; i++)
396  {
397  dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_x(k);
398  }
399  break;
400  case 2:
401  for (int o = 0, k = 0; k <= p; k++)
402  for (int j = 0; j <= p; j++)
403  for (int i = 0; i <= p; i++)
404  {
405  dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_x(k);
406  }
407  break;
408  case 3:
409  for (int o = 0, k = 0; k <= p; k++)
410  for (int j = 0; j <= p; j++)
411  for (int i = 0; i <= p; i++)
412  {
413  dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_x(k);
414  }
415  break;
416  case 4:
417  for (int o = 0, k = 0; k <= p; k++)
418  for (int j = 0; j <= p; j++)
419  for (int i = 0; i <= p; i++)
420  {
421  dofs(dof_map[o++]) = shape_x(i)*shape_x(j)*shape_y(k);
422  }
423  break;
424  case 5:
425  for (int o = 0, k = 0; k <= p; k++)
426  for (int j = 0; j <= p; j++)
427  for (int i = 0; i <= p; i++)
428  {
429  dofs(dof_map[o++]) = shape_y(i)*shape_x(j)*shape_y(k);
430  }
431  break;
432  case 6:
433  for (int o = 0, k = 0; k <= p; k++)
434  for (int j = 0; j <= p; j++)
435  for (int i = 0; i <= p; i++)
436  {
437  dofs(dof_map[o++]) = shape_y(i)*shape_y(j)*shape_y(k);
438  }
439  break;
440  case 7:
441  for (int o = 0, k = 0; k <= p; k++)
442  for (int j = 0; j <= p; j++)
443  for (int i = 0; i <= p; i++)
444  {
445  dofs(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_y(k);
446  }
447  break;
448  }
449 }
450
451 H1_TriangleElement::H1_TriangleElement(const int p, const int btype)
452  : NodalFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
453  FunctionSpace::Pk)
454 {
455  const double *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
456
458  shape_x.SetSize(p + 1);
459  shape_y.SetSize(p + 1);
460  shape_l.SetSize(p + 1);
461  dshape_x.SetSize(p + 1);
462  dshape_y.SetSize(p + 1);
463  dshape_l.SetSize(p + 1);
464  ddshape_x.SetSize(p + 1);
465  ddshape_y.SetSize(p + 1);
466  ddshape_l.SetSize(p + 1);
467  u.SetSize(dof);
468  du.SetSize(dof, dim);
469  ddu.SetSize(dof, (dim * (dim + 1)) / 2 );
470 #else
471  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
472 #endif
473
474  int p2p3 = 2*p + 3;
475  auto idx = [p2p3](int i, int j) { return ((p2p3-j)*j)/2+i; };
477
478  // vertices
479  lex_ordering[idx(0,0)] = 0;
480  Nodes.IntPoint(0).Set2(cp[0], cp[0]);
481  lex_ordering[idx(p,0)] = 1;
482  Nodes.IntPoint(1).Set2(cp[p], cp[0]);
483  lex_ordering[idx(0,p)] = 2;
484  Nodes.IntPoint(2).Set2(cp[0], cp[p]);
485
486  // edges
487  int o = 3;
488  for (int i = 1; i < p; i++)
489  {
490  lex_ordering[idx(i,0)] = o;
491  Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
492  }
493  for (int i = 1; i < p; i++)
494  {
495  lex_ordering[idx(p-i,i)] = o;
496  Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
497  }
498  for (int i = 1; i < p; i++)
499  {
500  lex_ordering[idx(0,p-i)] = o;
501  Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
502  }
503
504  // interior
505  for (int j = 1; j < p; j++)
506  for (int i = 1; i + j < p; i++)
507  {
508  const double w = cp[i] + cp[j] + cp[p-i-j];
509  lex_ordering[idx(i,j)] = o;
510  Nodes.IntPoint(o++).Set2(cp[i]/w, cp[j]/w);
511  }
512
513  DenseMatrix T(dof);
514  for (int k = 0; k < dof; k++)
515  {
517  poly1d.CalcBasis(p, ip.x, shape_x);
518  poly1d.CalcBasis(p, ip.y, shape_y);
519  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
520
521  o = 0;
522  for (int j = 0; j <= p; j++)
523  for (int i = 0; i + j <= p; i++)
524  {
525  T(o++, k) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
526  }
527  }
528
529  Ti.Factor(T);
530  // mfem::out << "H1_TriangleElement(" << p << ") : "; Ti.TestInversion();
531 }
532
534  Vector &shape) const
535 {
536  const int p = order;
537
539  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1), u(dof);
540 #endif
541
542  poly1d.CalcBasis(p, ip.x, shape_x);
543  poly1d.CalcBasis(p, ip.y, shape_y);
544  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
545
546  for (int o = 0, j = 0; j <= p; j++)
547  for (int i = 0; i + j <= p; i++)
548  {
549  u(o++) = shape_x(i)*shape_y(j)*shape_l(p-i-j);
550  }
551
552  Ti.Mult(u, shape);
553 }
554
556  DenseMatrix &dshape) const
557 {
558  const int p = order;
559
561  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
562  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
563  DenseMatrix du(dof, dim);
564 #endif
565
566  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
567  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
568  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l);
569
570  for (int o = 0, j = 0; j <= p; j++)
571  for (int i = 0; i + j <= p; i++)
572  {
573  int k = p - i - j;
574  du(o,0) = ((dshape_x(i)* shape_l(k)) -
575  ( shape_x(i)*dshape_l(k)))*shape_y(j);
576  du(o,1) = ((dshape_y(j)* shape_l(k)) -
577  ( shape_y(j)*dshape_l(k)))*shape_x(i);
578  o++;
579  }
580
581  Ti.Mult(du, dshape);
582 }
583
585  DenseMatrix &ddshape) const
586 {
587  const int p = order;
589  Vector shape_x(p + 1), shape_y(p + 1), shape_l(p + 1);
590  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_l(p + 1);
591  Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
592  DenseMatrix ddu(dof, dim);
593 #endif
594
595  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
596  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
597  poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l, dshape_l, ddshape_l);
598
599  for (int o = 0, j = 0; j <= p; j++)
600  for (int i = 0; i + j <= p; i++)
601  {
602  int k = p - i - j;
603  // u_xx, u_xy, u_yy
604  ddu(o,0) = ((ddshape_x(i) * shape_l(k)) - 2. * (dshape_x(i) * dshape_l(k)) +
605  (shape_x(i) * ddshape_l(k))) * shape_y(j);
606  ddu(o,1) = (((shape_x(i) * ddshape_l(k)) - dshape_x(i) * dshape_l(k)) * shape_y(
607  j)) + (((dshape_x(i) * shape_l(k)) - (shape_x(i) * dshape_l(k))) * dshape_y(j));
608  ddu(o,2) = ((ddshape_y(j) * shape_l(k)) - 2. * (dshape_y(j) * dshape_l(k)) +
609  (shape_y(j) * ddshape_l(k))) * shape_x(i);
610  o++;
611  }
612
613  Ti.Mult(ddu, ddshape);
614 }
615
616
618  : NodalFiniteElement(3, Geometry::TETRAHEDRON, ((p + 1)*(p + 2)*(p + 3))/6,
619  p, FunctionSpace::Pk)
620 {
621  const double *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
622
624  shape_x.SetSize(p + 1);
625  shape_y.SetSize(p + 1);
626  shape_z.SetSize(p + 1);
627  shape_l.SetSize(p + 1);
628  dshape_x.SetSize(p + 1);
629  dshape_y.SetSize(p + 1);
630  dshape_z.SetSize(p + 1);
631  dshape_l.SetSize(p + 1);
632  ddshape_x.SetSize(p + 1);
633  ddshape_y.SetSize(p + 1);
634  ddshape_z.SetSize(p + 1);
635  ddshape_l.SetSize(p + 1);
636  u.SetSize(dof);
637  du.SetSize(dof, dim);
638  ddu.SetSize(dof, (dim * (dim + 1)) / 2);
639 #else
640  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
641 #endif
642
643  auto tri = [](int k) { return (k*(k + 1))/2; };
644  auto tet = [](int k) { return (k*(k + 1)*(k + 2))/6; };
645  int ndof = tet(p+1);
646  auto idx = [tri, tet, p, ndof](int i, int j, int k)
647  {
648  return ndof - tet(p - k) - tri(p + 1 - k - j) + i;
649  };
650
652
653  // vertices
654  lex_ordering[idx(0,0,0)] = 0;
655  Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
656  lex_ordering[idx(p,0,0)] = 1;
657  Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
658  lex_ordering[idx(0,p,0)] = 2;
659  Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
660  lex_ordering[idx(0,0,p)] = 3;
661  Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
662
663  // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
664  int o = 4;
665  for (int i = 1; i < p; i++) // (0,1)
666  {
667  lex_ordering[idx(i,0,0)] = o;
668  Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
669  }
670  for (int i = 1; i < p; i++) // (0,2)
671  {
672  lex_ordering[idx(0,i,0)] = o;
673  Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
674  }
675  for (int i = 1; i < p; i++) // (0,3)
676  {
677  lex_ordering[idx(0,0,i)] = o;
678  Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
679  }
680  for (int i = 1; i < p; i++) // (1,2)
681  {
682  lex_ordering[idx(p-i,i,0)] = o;
683  Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
684  }
685  for (int i = 1; i < p; i++) // (1,3)
686  {
687  lex_ordering[idx(p-i,0,i)] = o;
688  Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
689  }
690  for (int i = 1; i < p; i++) // (2,3)
691  {
692  lex_ordering[idx(0,p-i,i)] = o;
693  Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
694  }
695
696  // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
697  for (int j = 1; j < p; j++)
698  for (int i = 1; i + j < p; i++) // (1,2,3)
699  {
700  lex_ordering[idx(p-i-j,i,j)] = o;
701  double w = cp[i] + cp[j] + cp[p-i-j];
702  Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
703  }
704  for (int j = 1; j < p; j++)
705  for (int i = 1; i + j < p; i++) // (0,3,2)
706  {
707  lex_ordering[idx(0,j,i)] = o;
708  double w = cp[i] + cp[j] + cp[p-i-j];
709  Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
710  }
711  for (int j = 1; j < p; j++)
712  for (int i = 1; i + j < p; i++) // (0,1,3)
713  {
714  lex_ordering[idx(i,0,j)] = o;
715  double w = cp[i] + cp[j] + cp[p-i-j];
716  Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
717  }
718  for (int j = 1; j < p; j++)
719  for (int i = 1; i + j < p; i++) // (0,2,1)
720  {
721  lex_ordering[idx(j,i,0)] = o;
722  double w = cp[i] + cp[j] + cp[p-i-j];
723  Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
724  }
725
726  // interior
727  for (int k = 1; k < p; k++)
728  for (int j = 1; j + k < p; j++)
729  for (int i = 1; i + j + k < p; i++)
730  {
731  lex_ordering[idx(i,j,k)] = o;
732  double w = cp[i] + cp[j] + cp[k] + cp[p-i-j-k];
733  Nodes.IntPoint(o++).Set3(cp[i]/w, cp[j]/w, cp[k]/w);
734  }
735
736  DenseMatrix T(dof);
737  for (int m = 0; m < dof; m++)
738  {
740  poly1d.CalcBasis(p, ip.x, shape_x);
741  poly1d.CalcBasis(p, ip.y, shape_y);
742  poly1d.CalcBasis(p, ip.z, shape_z);
743  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
744
745  o = 0;
746  for (int k = 0; k <= p; k++)
747  for (int j = 0; j + k <= p; j++)
748  for (int i = 0; i + j + k <= p; i++)
749  {
750  T(o++, m) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
751  }
752  }
753
754  Ti.Factor(T);
755  // mfem::out << "H1_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
756 }
757
759  Vector &shape) const
760 {
761  const int p = order;
762
764  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
765  Vector u(dof);
766 #endif
767
768  poly1d.CalcBasis(p, ip.x, shape_x);
769  poly1d.CalcBasis(p, ip.y, shape_y);
770  poly1d.CalcBasis(p, ip.z, shape_z);
771  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
772
773  for (int o = 0, k = 0; k <= p; k++)
774  for (int j = 0; j + k <= p; j++)
775  for (int i = 0; i + j + k <= p; i++)
776  {
777  u(o++) = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(p-i-j-k);
778  }
779
780  Ti.Mult(u, shape);
781 }
782
784  DenseMatrix &dshape) const
785 {
786  const int p = order;
787
789  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
790  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
791  DenseMatrix du(dof, dim);
792 #endif
793
794  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
795  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
796  poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
797  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
798
799  for (int o = 0, k = 0; k <= p; k++)
800  for (int j = 0; j + k <= p; j++)
801  for (int i = 0; i + j + k <= p; i++)
802  {
803  int l = p - i - j - k;
804  du(o,0) = ((dshape_x(i)* shape_l(l)) -
805  ( shape_x(i)*dshape_l(l)))*shape_y(j)*shape_z(k);
806  du(o,1) = ((dshape_y(j)* shape_l(l)) -
807  ( shape_y(j)*dshape_l(l)))*shape_x(i)*shape_z(k);
808  du(o,2) = ((dshape_z(k)* shape_l(l)) -
809  ( shape_z(k)*dshape_l(l)))*shape_x(i)*shape_y(j);
810  o++;
811  }
812
813  Ti.Mult(du, dshape);
814 }
815
817  DenseMatrix &ddshape) const
818 {
819  const int p = order;
820
822  Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1), shape_l(p + 1);
823  Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1), dshape_l(p + 1);
824  Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_z(p + 1), ddshape_l(p + 1);
825  DenseMatrix ddu(dof, ((dim + 1) * dim) / 2);
826 #endif
827
828  poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
829  poly1d.CalcBasis(p, ip.y, shape_y, dshape_y, ddshape_y);
830  poly1d.CalcBasis(p, ip.z, shape_z, dshape_z, ddshape_z);
831  poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l, ddshape_l);
832
833  for (int o = 0, k = 0; k <= p; k++)
834  for (int j = 0; j + k <= p; j++)
835  for (int i = 0; i + j + k <= p; i++)
836  {
837  // u_xx, u_xy, u_xz, u_yy, u_yz, u_zz
838  int l = p - i - j - k;
839  ddu(o,0) = ((ddshape_x(i) * shape_l(l)) - 2. * (dshape_x(i) * dshape_l(l)) +
840  (shape_x(i) * ddshape_l(l))) * shape_y(j) * shape_z(k);
841  ddu(o,1) = ((dshape_y(j) * ((dshape_x(i) * shape_l(l)) -
842  (shape_x(i) * dshape_l(l)))) +
843  (shape_y(j) * ((ddshape_l(l) * shape_x(i)) -
844  (dshape_x(i) * dshape_l(l)))))* shape_z(k);
845  ddu(o,2) = ((dshape_z(k) * ((dshape_x(i) * shape_l(l)) -
846  (shape_x(i) * dshape_l(l)))) +
847  (shape_z(k) * ((ddshape_l(l) * shape_x(i)) -
848  (dshape_x(i) * dshape_l(l)))))* shape_y(j);
849  ddu(o,3) = ((ddshape_y(j) * shape_l(l)) - 2. * (dshape_y(j) * dshape_l(l)) +
850  (shape_y(j) * ddshape_l(l))) * shape_x(i) * shape_z(k);
851  ddu(o,4) = ((dshape_z(k) * ((dshape_y(j) * shape_l(l)) -
852  (shape_y(j)*dshape_l(l))) ) +
853  (shape_z(k)* ((ddshape_l(l)*shape_y(j)) -
854  (dshape_y(j) * dshape_l(l)) ) ) )* shape_x(i);
855  ddu(o,5) = ((ddshape_z(k) * shape_l(l)) - 2. * (dshape_z(k) * dshape_l(l)) +
856  (shape_z(k) * ddshape_l(l))) * shape_y(j) * shape_x(i);
857  o++;
858  }
859  Ti.Mult(ddu, ddshape);
860 }
861
862 // TODO: use a FunctionSpace specific to wedges instead of Qk.
864  const int btype)
865  : NodalFiniteElement(3, Geometry::PRISM, ((p + 1)*(p + 1)*(p + 2))/2,
866  p, FunctionSpace::Qk),
867  TriangleFE(p, btype),
868  SegmentFE(p, btype)
869 {
871  t_shape.SetSize(TriangleFE.GetDof());
872  s_shape.SetSize(SegmentFE.GetDof());
873  t_dshape.SetSize(TriangleFE.GetDof(), 2);
874  s_dshape.SetSize(SegmentFE.GetDof(), 1);
875 #endif
876
877  t_dof.SetSize(dof);
878  s_dof.SetSize(dof);
879
880  int p2p3 = 2*p + 3, ntri = ((p + 1)*(p + 2))/2;
881  auto idx = [p2p3,ntri](int i, int j, int k)
882  {
883  return k*ntri + ((p2p3-j)*j)/2+i;
884  };
885
887  int o = 0;
888
889  // Nodal DoFs
890  lex_ordering[idx(0,0,0)] = o++;
891  lex_ordering[idx(p,0,0)] = o++;
892  lex_ordering[idx(0,p,0)] = o++;
893  lex_ordering[idx(0,0,p)] = o++;
894  lex_ordering[idx(p,0,p)] = o++;
895  lex_ordering[idx(0,p,p)] = o++;
896  t_dof[0] = 0; s_dof[0] = 0;
897  t_dof[1] = 1; s_dof[1] = 0;
898  t_dof[2] = 2; s_dof[2] = 0;
899  t_dof[3] = 0; s_dof[3] = 1;
900  t_dof[4] = 1; s_dof[4] = 1;
901  t_dof[5] = 2; s_dof[5] = 1;
902
903  // Edge DoFs
904  int k = 0;
905  int ne = p-1;
906  for (int i=1; i<p; i++)
907  {
908  lex_ordering[idx(i,0,0)] = o + 0*ne + k;
909  lex_ordering[idx(p-i,i,0)] = o + 1*ne + k;
910  lex_ordering[idx(0,p-i,0)] = o + 2*ne + k;
911  lex_ordering[idx(i,0,p)] = o + 3*ne + k;
912  lex_ordering[idx(p-i,i,p)] = o + 4*ne + k;
913  lex_ordering[idx(0,p-i,p)] = o + 5*ne + k;
914  lex_ordering[idx(0,0,i)] = o + 6*ne + k;
915  lex_ordering[idx(p,0,i)] = o + 7*ne + k;
916  lex_ordering[idx(0,p,i)] = o + 8*ne + k;
917  t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
918  t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
919  t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
920  t_dof[5 + 3 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 3 * ne + i] = 1;
921  t_dof[5 + 4 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 4 * ne + i] = 1;
922  t_dof[5 + 5 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 5 * ne + i] = 1;
923  t_dof[5 + 6 * ne + i] = 0; s_dof[5 + 6 * ne + i] = i + 1;
924  t_dof[5 + 7 * ne + i] = 1; s_dof[5 + 7 * ne + i] = i + 1;
925  t_dof[5 + 8 * ne + i] = 2; s_dof[5 + 8 * ne + i] = i + 1;
926  ++k;
927  }
928  o += 9*ne;
929
930  // Triangular Face DoFs
931  k=0;
932  int nt = (p-1)*(p-2)/2;
933  for (int j=1; j<p; j++)
934  {
935  for (int i=1; i<p-j; i++)
936  {
937  int l = j - p + (((2 * p - 1) - i) * i) / 2;
938  lex_ordering[idx(i,j,0)] = o+l;
939  lex_ordering[idx(i,j,p)] = o+nt+k;
940  t_dof[6 + 9 * ne + k] = 3 * p + l; s_dof[6 + 9 * ne + k] = 0;
941  t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
942  k++;
943  }
944  }
945  o += 2*nt;
946
948  k=0;
949  int nq = (p-1)*(p-1);
950  for (int j=1; j<p; j++)
951  {
952  for (int i=1; i<p; i++)
953  {
954  lex_ordering[idx(i,0,j)] = o+k;
955  lex_ordering[idx(p-i,i,j)] = o+nq+k;
956  lex_ordering[idx(0,p-i,j)] = o+2*nq+k;
957
958  t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
959  t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
960  t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
961
962  s_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 1 + j;
963  s_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 1 + j;
964  s_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 1 + j;
965
966  k++;
967  }
968  }
969  o += 3*nq;
970
971  // Interior DoFs
972  int m=0;
973  for (k=1; k<p; k++)
974  {
975  int l=0;
976  for (int j=1; j<p; j++)
977  {
978  for (int i=1; i+j<p; i++)
979  {
980  lex_ordering[idx(i,j,k)] = o++;
981  t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
982  s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
983  l++; m++;
984  }
985  }
986  }
987
988  // Define Nodes
989  const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
990  const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
991  for (int i=0; i<dof; i++)
992  {
993  Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
994  Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
995  Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
996  }
997 }
998
1000  Vector &shape) const
1001 {
1003  Vector t_shape(TriangleFE.GetDof());
1004  Vector s_shape(SegmentFE.GetDof());
1005 #endif
1006
1007  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1008
1009  TriangleFE.CalcShape(ip, t_shape);
1010  SegmentFE.CalcShape(ipz, s_shape);
1011
1012  for (int i=0; i<dof; i++)
1013  {
1014  shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
1015  }
1016 }
1017
1019  DenseMatrix &dshape) const
1020 {
1022  Vector t_shape(TriangleFE.GetDof());
1023  DenseMatrix t_dshape(TriangleFE.GetDof(), 2);
1024  Vector s_shape(SegmentFE.GetDof());
1025  DenseMatrix s_dshape(SegmentFE.GetDof(), 1);
1026 #endif
1027
1028  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1029
1030  TriangleFE.CalcShape(ip, t_shape);
1031  TriangleFE.CalcDShape(ip, t_dshape);
1032  SegmentFE.CalcShape(ipz, s_shape);
1033  SegmentFE.CalcDShape(ipz, s_dshape);
1034
1035  for (int i=0; i<dof; i++)
1036  {
1037  dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
1038  dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
1039  dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
1040  }
1041 }
1042
1043 }
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:314
Array< int > lex_ordering
Definition: fe_base.hpp:711
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
static double CalcDelta(const int p, const double x)
Evaluate a representation of a Delta function at point x.
Definition: fe_base.hpp:1127
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_h1.cpp:216
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:783
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3943
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
int dim
Dimension of reference space.
Definition: fe_base.hpp:236
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1675
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3897
STL namespace.
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:170
const double * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition: fe_base.hpp:1071
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:1018
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_h1.cpp:78
H1_HexahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_HexahedronElement of order p and BasisType btype.
Definition: fe_h1.cpp:265
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_h1.cpp:367
static void CalcBasis(const int p, const double x, double *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition: fe_base.hpp:1087
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:40
Class for standard nodal finite elements.
Definition: fe_base.hpp:708
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_h1.cpp:816
Poly_1D::Basis & basis1d
Definition: fe_base.hpp:1188
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:758
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:555
H1_QuadrilateralElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_QuadrilateralElement of order p and BasisType btype.
Definition: fe_h1.cpp:125
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:59
H1_SegmentElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_SegmentElement of order p and BasisType btype.
Definition: fe_h1.cpp:21
void Set2(const double x1, const double x2)
Definition: intrules.hpp:86
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_h1.cpp:192
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:76
IntegrationRule Nodes
Definition: fe_base.hpp:246
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:999
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
Poly_1D poly1d
Definition: fe.cpp:28
H1_TetrahedronElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TetrahedronElement of order p and BasisType btype.
Definition: fe_h1.cpp:617
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &ddshape) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_h1.cpp:584
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:533
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_h1.cpp:338
H1_TriangleElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_TriangleElement of order p and BasisType btype.
Definition: fe_h1.cpp:451
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:293
Class for integration point with weight.
Definition: intrules.hpp:31
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:243
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:151
H1_WedgeElement(const int p, const int btype=BasisType::GaussLobatto)
Construct the H1_WedgeElement of order p and BasisType btype.
Definition: fe_h1.cpp:863
static int VerifyClosed(int b_type)
Ensure that the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition: fe_base.hpp:624
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property). ...
Definition: fe_base.hpp:641
Vector data type.
Definition: vector.hpp:58
Describes the function space on each element.
Definition: fe_base.hpp:216
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_h1.cpp:97
double u(const Vector &xvec)
Definition: lor_mms.hpp:22
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:243