MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tfe.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
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 #ifndef MFEM_TEMPLATE_FINITE_ELEMENTS
13 #define MFEM_TEMPLATE_FINITE_ELEMENTS
14 
15 #include "../config/tconfig.hpp"
16 #include "fe_coll.hpp"
17 
18 namespace mfem
19 {
20 
21 // Templated finite element classes, cf. fe.?pp
22 
23 template <typename real_t>
24 void CalcShapeMatrix(const FiniteElement &fe, const IntegrationRule &ir,
25  real_t *B, const Array<int> *dof_map = NULL)
26 {
27  // - B must be (nip x dof) with column major storage
28  // - The inverse of dof_map is applied to reorder the local dofs.
29  int nip = ir.GetNPoints();
30  int dof = fe.GetDof();
31  Vector shape(dof);
32 
33  for (int ip = 0; ip < nip; ip++)
34  {
35  fe.CalcShape(ir.IntPoint(ip), shape);
36  for (int id = 0; id < dof; id++)
37  {
38  int orig_id = dof_map ? (*dof_map)[id] : id;
39  B[ip+nip*id] = shape(orig_id);
40  }
41  }
42 }
43 
44 template <typename real_t>
45 void CalcGradTensor(const FiniteElement &fe, const IntegrationRule &ir,
46  real_t *G, const Array<int> *dof_map = NULL)
47 {
48  // - G must be (nip x dim x dof) with column major storage
49  // - The inverse of dof_map is applied to reorder the local dofs.
50  int dim = fe.GetDim();
51  int nip = ir.GetNPoints();
52  int dof = fe.GetDof();
53  DenseMatrix dshape(dof, dim);
54 
55  for (int ip = 0; ip < nip; ip++)
56  {
57  fe.CalcDShape(ir.IntPoint(ip), dshape);
58  for (int id = 0; id < dof; id++)
59  {
60  int orig_id = dof_map ? (*dof_map)[id] : id;
61  for (int d = 0; d < dim; d++)
62  {
63  G[ip+nip*(d+dim*id)] = dshape(orig_id, d);
64  }
65  }
66  }
67 }
68 
69 template <typename real_t>
70 void CalcShapes(const FiniteElement &fe, const IntegrationRule &ir,
71  real_t *B, real_t *G, const Array<int> *dof_map)
72 {
73  if (B) { mfem::CalcShapeMatrix(fe, ir, B, dof_map); }
74  if (G) { mfem::CalcGradTensor(fe, ir, G, dof_map); }
75 }
76 
77 // H1 finite elements
78 
79 template <Geometry::Type G, int P>
81 
82 template <int P>
83 class H1_FiniteElement<Geometry::SEGMENT, P>
84 {
85 public:
87  static const int dim = 1;
88  static const int degree = P;
89  static const int dofs = P+1;
90 
91  static const bool tensor_prod = true;
92  static const int dofs_1d = P+1;
93 
94  // Type for run-time parameter for the constructor
95  typedef int parameter_type;
96 
97 protected:
100  parameter_type type; // run-time specified basis type
101  void Init(const parameter_type type_)
102  {
103  type = type_;
104  if (type == BasisType::Positive)
105  {
107  my_fe = fe;
108  my_dof_map = &fe->GetDofMap();
109  }
110  else
111  {
112  int pt_type = BasisType::GetQuadrature1D(type);
113  H1_SegmentElement *fe = new H1_SegmentElement(P, pt_type);
114  my_fe = fe;
115  my_dof_map = &fe->GetDofMap();
116  }
117  }
118 
119 public:
121  {
122  Init(type_);
123  }
125  {
126  const H1_FECollection *h1_fec =
127  dynamic_cast<const H1_FECollection *>(&fec);
128  MFEM_ASSERT(h1_fec, "invalid FiniteElementCollection");
129  Init(h1_fec->GetBasisType());
130  }
131  ~H1_FiniteElement() { delete my_fe; }
132 
133  template <typename real_t>
134  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
135  {
136  mfem::CalcShapes(*my_fe, ir, B, G, my_dof_map);
137  }
138  template <typename real_t>
139  void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
140  {
141  CalcShapes(ir, B, G);
142  }
143  const Array<int> *GetDofMap() const { return my_dof_map; }
144 };
145 
146 template <int P>
147 class H1_FiniteElement<Geometry::TRIANGLE, P>
148 {
149 public:
151  static const int dim = 2;
152  static const int degree = P;
153  static const int dofs = ((P + 1)*(P + 2))/2;
154 
155  static const bool tensor_prod = false;
156 
157  // Type for run-time parameter for the constructor
158  typedef int parameter_type;
159 
160 protected:
162  parameter_type type; // run-time specified basis type
163  void Init(const parameter_type type_)
164  {
165  type = type_;
166  if (type == BasisType::Positive)
167  {
168  my_fe = new H1Pos_TriangleElement(P);
169  }
170  else
171  {
172  int pt_type = BasisType::GetQuadrature1D(type);
173  my_fe = new H1_TriangleElement(P, pt_type);
174  }
175  }
176 
177 public:
179  {
180  Init(type_);
181  }
183  {
184  const H1_FECollection *h1_fec =
185  dynamic_cast<const H1_FECollection *>(&fec);
186  MFEM_ASSERT(h1_fec, "invalid FiniteElementCollection");
187  Init(h1_fec->GetBasisType());
188  }
189  ~H1_FiniteElement() { delete my_fe; }
190 
191  template <typename real_t>
192  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
193  {
194  mfem::CalcShapes(*my_fe, ir, B, G, NULL);
195  }
196  const Array<int> *GetDofMap() const { return NULL; }
197 };
198 
199 template <int P>
200 class H1_FiniteElement<Geometry::SQUARE, P>
201 {
202 public:
204  static const int dim = 2;
205  static const int degree = P;
206  static const int dofs = (P+1)*(P+1);
207 
208  static const bool tensor_prod = true;
209  static const int dofs_1d = P+1;
210 
211  // Type for run-time parameter for the constructor
212  typedef int parameter_type;
213 
214 protected:
215  const FiniteElement *my_fe, *my_fe_1d;
217  parameter_type type; // run-time specified basis type
218  void Init(const parameter_type type_)
219  {
220  type = type_;
221  if (type == BasisType::Positive)
222  {
224  my_fe = fe;
225  my_dof_map = &fe->GetDofMap();
226  my_fe_1d = new L2Pos_SegmentElement(P);
227  }
228  else
229  {
230  int pt_type = BasisType::GetQuadrature1D(type);
231  H1_QuadrilateralElement *fe = new H1_QuadrilateralElement(P, pt_type);
232  my_fe = fe;
233  my_dof_map = &fe->GetDofMap();
234  my_fe_1d = new L2_SegmentElement(P, pt_type);
235  }
236  }
237 
238 public:
240  {
241  Init(type_);
242  }
244  {
245  const H1_FECollection *h1_fec =
246  dynamic_cast<const H1_FECollection *>(&fec);
247  MFEM_ASSERT(h1_fec, "invalid FiniteElementCollection");
248  Init(h1_fec->GetBasisType());
249  }
250  ~H1_FiniteElement() { delete my_fe; delete my_fe_1d; }
251 
252  template <typename real_t>
253  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
254  {
255  mfem::CalcShapes(*my_fe, ir, B, G, my_dof_map);
256  }
257  template <typename real_t>
258  void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
259  {
260  mfem::CalcShapes(*my_fe_1d, ir, B, G, NULL);
261  }
262  const Array<int> *GetDofMap() const { return my_dof_map; }
263 };
264 
265 template <int P>
266 class H1_FiniteElement<Geometry::TETRAHEDRON, P>
267 {
268 public:
270  static const int dim = 3;
271  static const int degree = P;
272  static const int dofs = ((P + 1)*(P + 2)*(P + 3))/6;
273 
274  static const bool tensor_prod = false;
275 
276  // Type for run-time parameter for the constructor
277  typedef int parameter_type;
278 
279 protected:
281  parameter_type type; // run-time specified basis type
282  void Init(const parameter_type type_)
283  {
284  type = type_;
285  if (type == BasisType::Positive)
286  {
287  my_fe = new H1Pos_TetrahedronElement(P);
288  }
289  else
290  {
291  int pt_type = BasisType::GetQuadrature1D(type);
292  my_fe = new H1_TetrahedronElement(P, pt_type);
293  }
294  }
295 
296 public:
298  {
299  Init(type_);
300  }
302  {
303  const H1_FECollection *h1_fec =
304  dynamic_cast<const H1_FECollection *>(&fec);
305  MFEM_ASSERT(h1_fec, "invalid FiniteElementCollection");
306  Init(h1_fec->GetBasisType());
307  }
308  ~H1_FiniteElement() { delete my_fe; }
309 
310  template <typename real_t>
311  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
312  {
313  mfem::CalcShapes(*my_fe, ir, B, G, NULL);
314  }
315  const Array<int> *GetDofMap() const { return NULL; }
316 };
317 
318 template <int P>
319 class H1_FiniteElement<Geometry::CUBE, P>
320 {
321 public:
323  static const int dim = 3;
324  static const int degree = P;
325  static const int dofs = (P+1)*(P+1)*(P+1);
326 
327  static const bool tensor_prod = true;
328  static const int dofs_1d = P+1;
329 
330  // Type for run-time parameter for the constructor
331  typedef int parameter_type;
332 
333 protected:
334  const FiniteElement *my_fe, *my_fe_1d;
336  parameter_type type; // run-time specified basis type
337 
338  void Init(const parameter_type type_)
339  {
340  type = type_;
341  if (type == BasisType::Positive)
342  {
344  my_fe = fe;
345  my_dof_map = &fe->GetDofMap();
346  my_fe_1d = new L2Pos_SegmentElement(P);
347  }
348  else
349  {
350  int pt_type = BasisType::GetQuadrature1D(type);
351  H1_HexahedronElement *fe = new H1_HexahedronElement(P, pt_type);
352  my_fe = fe;
353  my_dof_map = &fe->GetDofMap();
354  my_fe_1d = new L2_SegmentElement(P, pt_type);
355  }
356  }
357 
358 public:
360  {
361  Init(type_);
362  }
364  {
365  const H1_FECollection *h1_fec =
366  dynamic_cast<const H1_FECollection *>(&fec);
367  MFEM_ASSERT(h1_fec, "invalid FiniteElementCollection");
368  Init(h1_fec->GetBasisType());
369  }
370  ~H1_FiniteElement() { delete my_fe; delete my_fe_1d; }
371 
372  template <typename real_t>
373  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
374  {
375  mfem::CalcShapes(*my_fe, ir, B, G, my_dof_map);
376  }
377  template <typename real_t>
378  void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
379  {
380  mfem::CalcShapes(*my_fe_1d, ir, B, G, NULL);
381  }
382  const Array<int> *GetDofMap() const { return my_dof_map; }
383 };
384 
385 
386 // L2 finite elements
387 
388 template <Geometry::Type G, int P, typename L2_FE_type, typename L2Pos_FE_type,
389  int DOFS, bool TP>
391 {
392 public:
393  static const Geometry::Type geom = G;
395  static const int degree = P;
396  static const int dofs = DOFS;
397 
398  static const bool tensor_prod = TP;
399  static const int dofs_1d = P+1;
400 
401  // Type for run-time parameter for the constructor
402  typedef int parameter_type;
403 
404 protected:
406  parameter_type type; // run-time specified basis type
407 
408  void Init(const parameter_type type_)
409  {
410  type = type_;
411  switch (type)
412  {
413  case BasisType::Positive:
414  my_fe = new L2Pos_FE_type(P);
415  my_fe_1d = (TP && dim != 1) ? new L2Pos_SegmentElement(P) : NULL;
416  break;
417 
418  default:
419  int pt_type = BasisType::GetQuadrature1D(type);
420  my_fe = new L2_FE_type(P, pt_type);
421  my_fe_1d = (TP && dim != 1) ? new L2_SegmentElement(P, pt_type) :
422  NULL;
423  }
424  }
425 
427  { Init(type); }
428 
430  {
431  const L2_FECollection *l2_fec =
432  dynamic_cast<const L2_FECollection *>(&fec);
433  MFEM_ASSERT(l2_fec, "invalid FiniteElementCollection");
434  Init(l2_fec->GetBasisType());
435  }
436 
437  ~L2_FiniteElement_base() { delete my_fe; delete my_fe_1d; }
438 
439 public:
440  template <typename real_t>
441  void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *Grad) const
442  {
443  mfem::CalcShapes(*my_fe, ir, B, Grad, NULL);
444  }
445  template <typename real_t>
446  void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *Grad) const
447  {
448  mfem::CalcShapes(dim == 1 ? *my_fe : *my_fe_1d, ir, B, Grad, NULL);
449  }
450  const Array<int> *GetDofMap() const { return NULL; }
451 };
452 
453 
454 template <Geometry::Type G, int P>
456 
457 
458 template <int P>
459 class L2_FiniteElement<Geometry::SEGMENT, P>
460  : public L2_FiniteElement_base<
461  Geometry::SEGMENT,P,L2_SegmentElement,L2Pos_SegmentElement,P+1,true>
462 {
463 protected:
466 public:
469  : base_class(type_) { }
471  : base_class(fec) { }
472 };
473 
474 
475 template <int P>
476 class L2_FiniteElement<Geometry::TRIANGLE, P>
477  : public L2_FiniteElement_base<Geometry::TRIANGLE,P,L2_TriangleElement,
478  L2Pos_TriangleElement,((P+1)*(P+2))/2,false>
479 {
480 protected:
482  L2Pos_TriangleElement,((P+1)*(P+2))/2,false> base_class;
483 public:
486  : base_class(type_) { }
488  : base_class(fec) { }
489 };
490 
491 
492 template <int P>
493 class L2_FiniteElement<Geometry::SQUARE, P>
494  : public L2_FiniteElement_base<Geometry::SQUARE,P,L2_QuadrilateralElement,
495  L2Pos_QuadrilateralElement,(P+1)*(P+1),true>
496 {
497 protected:
500 public:
503  : base_class(type_) { }
505  : base_class(fec) { }
506 };
507 
508 
509 template <int P>
510 class L2_FiniteElement<Geometry::TETRAHEDRON, P>
511  : public L2_FiniteElement_base<Geometry::TETRAHEDRON,P,L2_TetrahedronElement,
512  L2Pos_TetrahedronElement,((P+1)*(P+2)*(P+3))/6,false>
513 {
514 protected:
516  L2Pos_TetrahedronElement,((P+1)*(P+2)*(P+3))/6,false> base_class;
517 public:
520  : base_class(type_) { }
522  : base_class(fec) { }
523 };
524 
525 
526 template <int P>
527 class L2_FiniteElement<Geometry::CUBE, P>
528  : public L2_FiniteElement_base<Geometry::CUBE,P,L2_HexahedronElement,
529  L2Pos_HexahedronElement,(P+1)*(P+1)*(P+1),true>
530 {
531 protected:
533  L2Pos_HexahedronElement,(P+1)*(P+1)*(P+1),true> base_class;
534 public:
537  : base_class(type_) { }
539  : base_class(fec) { }
540 };
541 
542 } // namespace mfem
543 
544 #endif // MFEM_TEMPLATE_FINITE_ELEMENTS
const Array< int > * GetDofMap() const
Definition: tfe.hpp:315
void Init(const parameter_type type_)
Definition: tfe.hpp:282
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
Abstract class for Finite Elements.
Definition: fe.hpp:232
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:308
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
L2_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:504
static const Geometry::Type geom
Definition: tfe.hpp:393
static const int dim
Definition: tfe.hpp:394
L2_FiniteElement_base< Geometry::SEGMENT, P, L2_SegmentElement, L2Pos_SegmentElement, P+1, true > base_class
Definition: tfe.hpp:465
const Array< int > * GetDofMap() const
Definition: tfe.hpp:262
L2_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:470
L2_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:538
const Geometry::Type geom
Definition: ex1.cpp:40
void Init(const parameter_type type_)
Definition: tfe.hpp:408
H1_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:301
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void CalcShapeMatrix(const FiniteElement &fe, const IntegrationRule &ir, real_t *B, const Array< int > *dof_map=NULL)
Definition: tfe.hpp:24
static int GetQuadrature1D(int b_type)
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition: fe.hpp:61
const Array< int > * GetDofMap() const
Definition: tfe.hpp:196
H1_FiniteElement(const parameter_type type_=BasisType::GaussLobatto)
Definition: tfe.hpp:297
H1_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:182
static const int degree
Definition: tfe.hpp:395
const Array< int > * GetDofMap() const
Definition: tfe.hpp:382
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:311
base_class::parameter_type parameter_type
Definition: tfe.hpp:484
L2_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:521
H1_FiniteElement(const parameter_type type_=BasisType::GaussLobatto)
Definition: tfe.hpp:120
void Init(const parameter_type type_)
Definition: tfe.hpp:163
L2_FiniteElement_base(const FiniteElementCollection &fec)
Definition: tfe.hpp:429
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
H1_FiniteElement(const parameter_type type_=BasisType::GaussLobatto)
Definition: tfe.hpp:178
int GetBasisType() const
Definition: fe_coll.hpp:182
int GetBasisType() const
Definition: fe_coll.hpp:107
void CalcGradTensor(const FiniteElement &fe, const IntegrationRule &ir, real_t *G, const Array< int > *dof_map=NULL)
Definition: tfe.hpp:45
H1_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:363
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
static const bool tensor_prod
Definition: tfe.hpp:398
L2_FiniteElement_base< Geometry::TETRAHEDRON, P, L2_TetrahedronElement, L2Pos_TetrahedronElement,((P+1)*(P+2)*(P+3))/6, false > base_class
Definition: tfe.hpp:516
base_class::parameter_type parameter_type
Definition: tfe.hpp:518
H1_FiniteElement(const parameter_type type_=BasisType::GaussLobatto)
Definition: tfe.hpp:239
base_class::parameter_type parameter_type
Definition: tfe.hpp:501
L2_FiniteElement(const parameter_type type_=BasisType::GaussLegendre)
Definition: tfe.hpp:536
base_class::parameter_type parameter_type
Definition: tfe.hpp:535
static const int dofs_1d
Definition: tfe.hpp:399
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:134
void Init(const parameter_type type_)
Definition: tfe.hpp:218
L2_FiniteElement_base(const parameter_type type)
Definition: tfe.hpp:426
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *Grad) const
Definition: tfe.hpp:441
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
parameter_type type
Definition: tfe.hpp:406
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:192
A class to initialize the size of a Tensor.
Definition: dtensor.hpp:54
const Array< int > & GetDofMap() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe.hpp:1850
L2_FiniteElement(const parameter_type type_=BasisType::GaussLegendre)
Definition: tfe.hpp:485
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:253
L2_FiniteElement(const parameter_type type_=BasisType::GaussLegendre)
Definition: tfe.hpp:468
H1_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:243
void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:258
H1_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:124
void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:378
int dim
Definition: ex24.cpp:43
void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:139
L2_FiniteElement(const parameter_type type_=BasisType::GaussLegendre)
Definition: tfe.hpp:519
H1_FiniteElement(const parameter_type type_=BasisType::GaussLobatto)
Definition: tfe.hpp:359
const FiniteElement * my_fe_1d
Definition: tfe.hpp:405
Vector data type.
Definition: vector.hpp:48
const FiniteElement * my_fe
Definition: tfe.hpp:405
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
static const int dofs
Definition: tfe.hpp:396
Bernstein polynomials.
Definition: fe.hpp:35
base_class::parameter_type parameter_type
Definition: tfe.hpp:467
void CalcShapes(const IntegrationRule &ir, real_t *B, real_t *G) const
Definition: tfe.hpp:373
L2_FiniteElement(const parameter_type type_=BasisType::GaussLegendre)
Definition: tfe.hpp:502
void Init(const parameter_type type_)
Definition: tfe.hpp:338
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
L2_FiniteElement_base< Geometry::SQUARE, P, L2_QuadrilateralElement, L2Pos_QuadrilateralElement,(P+1)*(P+1), true > base_class
Definition: tfe.hpp:499
void CalcShapes(const FiniteElement &fe, const IntegrationRule &ir, real_t *B, real_t *G, const Array< int > *dof_map)
Definition: tfe.hpp:70
void Calc1DShapes(const IntegrationRule &ir, real_t *B, real_t *Grad) const
Definition: tfe.hpp:446
const Array< int > * GetDofMap() const
Definition: tfe.hpp:143
void Init(const parameter_type type_)
Definition: tfe.hpp:101
L2_FiniteElement(const FiniteElementCollection &fec)
Definition: tfe.hpp:487
const Array< int > * GetDofMap() const
Definition: tfe.hpp:450
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143