MFEM  v4.4.0 Finite element discretization library
fe_base.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #ifndef MFEM_FE_BASE
13 #define MFEM_FE_BASE
14
15 #include "../intrules.hpp"
16 #include "../geom.hpp"
17
18 #include <map>
19
20 namespace mfem
21 {
22
23 /// Possible basis types. Note that not all elements can use all BasisType(s).
24 class BasisType
25 {
26 public:
27  enum
28  {
29  Invalid = -1,
30  GaussLegendre = 0, ///< Open type
31  GaussLobatto = 1, ///< Closed type
32  Positive = 2, ///< Bernstein polynomials
33  OpenUniform = 3, ///< Nodes: x_i = (i+1)/(n+1), i=0,...,n-1
34  ClosedUniform = 4, ///< Nodes: x_i = i/(n-1), i=0,...,n-1
35  OpenHalfUniform = 5, ///< Nodes: x_i = (i+1/2)/n, i=0,...,n-1
36  Serendipity = 6, ///< Serendipity basis (squares / cubes)
37  ClosedGL = 7, ///< Closed GaussLegendre
38  IntegratedGLL = 8, ///< Integrated GLL indicator functions
39  NumBasisTypes = 9 /**< Keep track of maximum types to prevent
40  hard-coding */
41  };
42  /** @brief If the input does not represents a valid BasisType, abort with an
43  error; otherwise return the input. */
44  static int Check(int b_type)
45  {
46  MFEM_VERIFY(0 <= b_type && b_type < NumBasisTypes,
47  "unknown BasisType: " << b_type);
48  return b_type;
49  }
50  /** @brief If the input does not represents a valid nodal BasisType, abort
51  with an error; otherwise return the input. */
52  static int CheckNodal(int b_type)
53  {
54  MFEM_VERIFY(Check(b_type) != Positive && b_type != IntegratedGLL,
55  "invalid nodal BasisType: " << Name(b_type));
56  return b_type;
57  }
58  /** @brief Get the corresponding Quadrature1D constant, when that makes
59  sense; otherwise return Quadrature1D::Invalid. */
61  {
62  switch (b_type)
63  {
66  case Positive: return Quadrature1D::ClosedUniform; // <-----
73  }
75  }
76  /// Return the nodal BasisType corresponding to the Quadrature1D type.
77  static int GetNodalBasis(int qpt_type)
78  {
79  switch (qpt_type)
80  {
87  }
88  return Invalid;
89  }
90  /// Check and convert a BasisType constant to a string identifier.
91  static const char *Name(int b_type)
92  {
93  static const char *name[] =
94  {
95  "Gauss-Legendre", "Gauss-Lobatto", "Positive (Bernstein)",
96  "Open uniform", "Closed uniform", "Open half uniform",
97  "Serendipity", "Closed Gauss-Legendre",
98  "Integrated Gauss-Lobatto indicator"
99  };
100  return name[Check(b_type)];
101  }
102  /// Check and convert a BasisType constant to a char basis identifier.
103  static char GetChar(int b_type)
104  {
105  static const char ident[]
106  = { 'g', 'G', 'P', 'u', 'U', 'o', 'S', 'c', 'i' };
107  return ident[Check(b_type)];
108  }
109  /// Convert char basis identifier to a BasisType constant.
110  static int GetType(char b_ident)
111  {
112  switch (b_ident)
113  {
114  case 'g': return GaussLegendre;
115  case 'G': return GaussLobatto;
116  case 's': return GaussLobatto;
117  case 'P': return Positive;
118  case 'u': return OpenUniform;
119  case 'U': return ClosedUniform;
120  case 'o': return OpenHalfUniform;
121  case 'S': return Serendipity;
122  case 'c': return ClosedGL;
123  case 'i': return IntegratedGLL;
124  }
125  MFEM_ABORT("unknown BasisType identifier");
126  return -1;
127  }
128 };
129
130
131 /** @brief Structure representing the matrices/tensors needed to evaluate (in
132  reference space) the values, gradients, divergences, or curls of a
133  FiniteElement at a the quadrature points of a given IntegrationRule. */
134 /** Object of this type are typically created and owned by the respective
135  FiniteElement object. */
137 {
138 public:
139  /// The FiniteElement that created and owns this object.
140  /** This pointer is not owned. */
141  const class FiniteElement *FE;
142
143  /** @brief IntegrationRule that defines the quadrature points at which the
144  basis functions of the #FE are evaluated. */
145  /** This pointer is not owned. */
147
148  /// Type of data stored in the arrays #B, #Bt, #G, and #Gt.
149  enum Mode
150  {
151  /** @brief Full multidimensional representation which does not use tensor
152  product structure. The ordering of the degrees of freedom is as
153  defined by #FE */
155
156  /** @brief Tensor product representation using 1D matrices/tensors with
157  dimensions using 1D number of quadrature points and degrees of
158  freedom. */
159  /** When representing a vector-valued FiniteElement, two DofToQuad objects
160  are used to describe the "closed" and "open" 1D basis functions
161  (TODO). */
163  };
164
165  /// Describes the contents of the #B, #Bt, #G, and #Gt arrays, see #Mode.
167
168  /** @brief Number of degrees of freedom = number of basis functions. When
169  #mode is TENSOR, this is the 1D number. */
170  int ndof;
171
172  /** @brief Number of quadrature points. When #mode is TENSOR, this is the 1D
173  number. */
174  int nqpt;
175
176  /// Basis functions evaluated at quadrature points.
177  /** The storage layout is column-major with dimensions:
178  - #nqpt x #ndof, for scalar elements, or
179  - #nqpt x dim x #ndof, for vector elements, (TODO)
180
181  where
182
183  - dim = dimension of the finite element reference space when #mode is
184  FULL, and dim = 1 when #mode is TENSOR. */
186
187  /// Transpose of #B.
188  /** The storage layout is column-major with dimensions:
189  - #ndof x #nqpt, for scalar elements, or
190  - #ndof x #nqpt x dim, for vector elements (TODO). */
192
193  /** @brief Gradients/divergences/curls of basis functions evaluated at
195  /** The storage layout is column-major with dimensions:
196  - #nqpt x dim x #ndof, for scalar elements, or
197  - #nqpt x #ndof, for H(div) vector elements (TODO), or
198  - #nqpt x cdim x #ndof, for H(curl) vector elements (TODO),
199
200  where
201
202  - dim = dimension of the finite element reference space when #mode is
203  FULL, and 1 when #mode is TENSOR,
204  - cdim = 1/1/3 in 1D/2D/3D, respectively, when #mode is FULL, and cdim =
205  1 when #mode is TENSOR. */
207
208  /// Transpose of #G.
209  /** The storage layout is column-major with dimensions:
210  - #ndof x #nqpt x dim, for scalar elements, or
211  - #ndof x #nqpt, for H(div) vector elements (TODO), or
212  - #ndof x #nqpt x cdim, for H(curl) vector elements (TODO). */
214 };
215
216
217 /// Describes the function space on each element
219 {
220 public:
221  enum
222  {
223  Pk, ///< Polynomials of order k
224  Qk, ///< Tensor products of polynomials of order k
225  rQk ///< Refined tensor products of polynomials of order k
226  };
227 };
228
229 class ElementTransformation;
230 class Coefficient;
231 class VectorCoefficient;
232 class MatrixCoefficient;
233
234 /// Abstract class for all finite elements.
236 {
237 protected:
238  int dim; ///< Dimension of reference space
239  int vdim; ///< Vector dimension of vector-valued basis functions
240  int cdim; ///< Dimension of curl for vector-valued basis functions
241  Geometry::Type geom_type; ///< Geometry::Type of the reference element
244  mutable
245  int dof, ///< Number of degrees of freedom
246  order; ///< Order/degree of the shape functions
247  mutable int orders[Geometry::MaxDim]; ///< Anisotropic orders
250  mutable DenseMatrix vshape; // Dof x VDim
251 #endif
252  /// Container for all DofToQuad objects created by the FiniteElement.
253  /** Multiple DofToQuad objects may be needed when different quadrature rules
254  or different DofToQuad::Mode are used. */
256
257 public:
258  /// Enumeration for range_type and deriv_range_type
260
261  /** @brief Enumeration for MapType: defines how reference functions are
262  mapped to physical space.
263
264  A reference function \f$\hat u(\hat x) \f$ can be mapped to a function
265  \f$u(x) \f$ on a general physical element in following ways:
266  - \f$x = T(\hat x) \f$ is the image of the reference point \f$\hat x \f$
267  - \f$J = J(\hat x) \f$ is the Jacobian matrix of the transformation T
268  - \f$w = w(\hat x) = det(J) \f$ is the transformation weight factor for square J
269  - \f$w = w(\hat x) = det(J^t J)^{1/2} \f$ is the transformation weight factor in general
270  */
271  enum MapType
272  {
273  VALUE, /**< For scalar fields; preserves point values
274  \f$u(x) = \hat u(\hat x) \f$ */
275  INTEGRAL, /**< For scalar fields; preserves volume integrals
276  \f$u(x) = (1/w) \hat u(\hat x) \f$ */
277  H_DIV, /**< For vector fields; preserves surface integrals of the
278  normal component \f$u(x) = (J/w) \hat u(\hat x) \f$ */
279  H_CURL /**< For vector fields; preserves line integrals of the
280  tangential component
281  \f$u(x) = J^{-t} \hat u(\hat x) \f$ (square J),
282  \f$u(x) = J(J^t J)^{-1} \hat u(\hat x) \f$ (general J) */
283  };
284
285  /** @brief Enumeration for DerivType: defines which derivative method
286  is implemented.
287
288  Each FiniteElement class implements up to one type of derivative. The
289  value returned by GetDerivType() indicates which derivative method is
290  implemented.
291  */
293  {
294  NONE, ///< No derivatives implemented
295  GRAD, ///< Implements CalcDShape methods
296  DIV, ///< Implements CalcDivShape methods
297  CURL ///< Implements CalcCurlShape methods
298  };
299
300  /** @brief Construct FiniteElement with given
301  @param D Reference space dimension
302  @param G Geometry type (of type Geometry::Type)
303  @param Do Number of degrees of freedom in the FiniteElement
304  @param O Order/degree of the FiniteElement
305  @param F FunctionSpace type of the FiniteElement
306  */
307  FiniteElement(int D, Geometry::Type G, int Do, int O,
308  int F = FunctionSpace::Pk);
309
310  /// Returns the reference space dimension for the finite element
311  int GetDim() const { return dim; }
312
313  /// Returns the vector dimension for vector-valued finite elements
314  int GetVDim() const { return vdim; }
315
316  /// Returns the dimension of the curl for vector-valued finite elements
317  int GetCurlDim() const { return cdim; }
318
319  /// Returns the Geometry::Type of the reference element
321
322  /// Returns the number of degrees of freedom in the finite element
323  int GetDof() const { return dof; }
324
325  /** @brief Returns the order of the finite element. In the case of
326  anisotropic orders, returns the maximum order. */
327  int GetOrder() const { return order; }
328
329  /** @brief Returns true if the FiniteElement basis *may be using* different
330  orders/degrees in different spatial directions. */
331  bool HasAnisotropicOrders() const { return orders[0] != -1; }
332
333  /// Returns an array containing the anisotropic orders/degrees.
334  const int *GetAnisotropicOrders() const { return orders; }
335
336  /// Returns the type of FunctionSpace on the element.
337  int Space() const { return func_space; }
338
339  /// Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
340  int GetRangeType() const { return range_type; }
341
342  /** @brief Returns the FiniteElement::RangeType of the element derivative, either
343  SCALAR or VECTOR. */
344  int GetDerivRangeType() const { return deriv_range_type; }
345
346  /** @brief Returns the FiniteElement::MapType of the element describing how reference
347  functions are mapped to physical space, one of {VALUE, INTEGRAL
348  H_DIV, H_CURL}. */
349  int GetMapType() const { return map_type; }
350
351
352  /** @brief Returns the FiniteElement::DerivType of the element describing the
353  spatial derivative method implemented, one of {NONE, GRAD,
354  DIV, CURL}. */
355  int GetDerivType() const { return deriv_type; }
356
357  /** @brief Returns the FiniteElement::DerivType of the element describing how
358  reference function derivatives are mapped to physical space, one of {VALUE,
359  INTEGRAL, H_DIV, H_CURL}. */
360  int GetDerivMapType() const { return deriv_map_type; }
361
362  /** @brief Evaluate the values of all shape functions of a scalar finite
363  element in reference space at the given point @a ip. */
364  /** The size (#dof) of the result Vector @a shape must be set in advance. */
365  virtual void CalcShape(const IntegrationPoint &ip,
366  Vector &shape) const = 0;
367
368  /** @brief Evaluate the values of all shape functions of a scalar finite
369  element in physical space at the point described by @a Trans. */
370  /** The size (#dof) of the result Vector @a shape must be set in advance. */
371  void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const;
372
373  /** @brief Evaluate the gradients of all shape functions of a scalar finite
374  element in reference space at the given point @a ip. */
375  /** Each row of the result DenseMatrix @a dshape contains the derivatives of
376  one shape function. The size (#dof x #dim) of @a dshape must be set in
378  virtual void CalcDShape(const IntegrationPoint &ip,
379  DenseMatrix &dshape) const = 0;
380
381  /** @brief Evaluate the gradients of all shape functions of a scalar finite
382  element in physical space at the point described by @a Trans. */
383  /** Each row of the result DenseMatrix @a dshape contains the derivatives of
384  one shape function. The size (#dof x SDim) of @a dshape must be set in
385  advance, where SDim >= #dim is the physical space dimension as described
386  by @a Trans. */
388
389  /// Get a const reference to the nodes of the element
390  const IntegrationRule & GetNodes() const { return Nodes; }
391
392  // virtual functions for finite elements on vector spaces
393
394  /** @brief Evaluate the values of all shape functions of a *vector* finite
395  element in reference space at the given point @a ip. */
396  /** Each row of the result DenseMatrix @a shape contains the components of
397  one vector shape function. The size (#dof x #dim) of @a shape must be set
399  virtual void CalcVShape(const IntegrationPoint &ip,
400  DenseMatrix &shape) const;
401
402  /** @brief Evaluate the values of all shape functions of a *vector* finite
403  element in physical space at the point described by @a Trans. */
404  /** Each row of the result DenseMatrix @a shape contains the components of
405  one vector shape function. The size (#dof x SDim) of @a shape must be set
406  in advance, where SDim >= #dim is the physical space dimension as
407  described by @a Trans. */
408  virtual void CalcVShape(ElementTransformation &Trans,
409  DenseMatrix &shape) const;
410
411  /// Equivalent to the CalcVShape() method with the same arguments.
413  { CalcVShape(Trans, shape); }
414
415  /** @brief Evaluate the divergence of all shape functions of a *vector*
416  finite element in reference space at the given point @a ip. */
417  /** The size (#dof) of the result Vector @a divshape must be set in advance.
418  */
419  virtual void CalcDivShape(const IntegrationPoint &ip,
420  Vector &divshape) const;
421
422  /** @brief Evaluate the divergence of all shape functions of a *vector*
423  finite element in physical space at the point described by @a Trans. */
424  /** The size (#dof) of the result Vector @a divshape must be set in advance.
425  */
426  void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const;
427
428  /** @brief Evaluate the curl of all shape functions of a *vector* finite
429  element in reference space at the given point @a ip. */
430  /** Each row of the result DenseMatrix @a curl_shape contains the components
431  of the curl of one vector shape function. The size (#dof x CDim) of
432  @a curl_shape must be set in advance, where CDim = 3 for #dim = 3 and
433  CDim = 1 for #dim = 2. */
434  virtual void CalcCurlShape(const IntegrationPoint &ip,
435  DenseMatrix &curl_shape) const;
436
437  /** @brief Evaluate the curl of all shape functions of a *vector* finite
438  element in physical space at the point described by @a Trans. */
439  /** Each row of the result DenseMatrix @a curl_shape contains the components
440  of the curl of one vector shape function. The size (#dof x CDim) of
441  @a curl_shape must be set in advance, where CDim = 3 for #dim = 3 and
442  CDim = 1 for #dim = 2. */
444  DenseMatrix &curl_shape) const;
445
446  /** @brief Get the dofs associated with the given @a face.
447  @a *dofs is set to an internal array of the local dofc on the
448  face, while *ndofs is set to the number of dofs on that face.
449  */
450  virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
451
452  /** @brief Evaluate the Hessians of all shape functions of a scalar finite
453  element in reference space at the given point @a ip. */
454  /** Each row of the result DenseMatrix @a Hessian contains upper triangular
455  part of the Hessian of one shape function.
456  The order in 2D is {u_xx, u_xy, u_yy}.
457  The size (#dof x (#dim (#dim+1)/2) of @a Hessian must be set in advance.*/
458  virtual void CalcHessian (const IntegrationPoint &ip,
459  DenseMatrix &Hessian) const;
460
461  /** @brief Evaluate the Hessian of all shape functions of a scalar finite
462  element in reference space at the given point @a ip. */
463  /** The size (#dof, #dim*(#dim+1)/2) of @a Hessian must be set in advance. */
465  DenseMatrix& Hessian) const;
466
467  /** @brief Evaluate the Laplacian of all shape functions of a scalar finite
468  element in reference space at the given point @a ip. */
469  /** The size (#dof) of @a Laplacian must be set in advance. */
471  Vector& Laplacian) const;
472
474  Vector& Laplacian) const;
475
476  /** @brief Return the local interpolation matrix @a I (Dof x Dof) where the
477  fine element is the image of the base geometry under the given
478  transformation. */
480  DenseMatrix &I) const;
481
482  /** @brief Return a local restriction matrix @a R (Dof x Dof) mapping fine
483  dofs to coarse dofs.
484
485  The fine element is the image of the base geometry under the given
486  transformation, @a Trans.
487
488  The assumption in this method is that a subset of the coarse dofs can be
489  expressed only in terms of the dofs of the given fine element.
490
491  Rows in @a R corresponding to coarse dofs that cannot be expressed in
492  terms of the fine dofs will be marked as invalid by setting the first
493  entry (column 0) in the row to infinity().
494
495  This method assumes that the dimensions of @a R are set before it is
496  called. */
498  DenseMatrix &R) const;
499
500  /** @brief Return interpolation matrix, @a I, which maps dofs from a coarse
501  element, @a fe, to the fine dofs on @a this finite element. */
502  /** @a Trans represents the mapping from the reference element of @a this
503  element into a subset of the reference space of the element @a fe, thus
504  allowing the "coarse" FiniteElement to be different from the "fine"
505  FiniteElement as when h-refinement is combined with p-refinement or
506  p-derefinement. It is assumed that both finite elements use the same
507  FiniteElement::MapType. */
508  virtual void GetTransferMatrix(const FiniteElement &fe,
510  DenseMatrix &I) const;
511
512  /** @brief Given a coefficient and a transformation, compute its projection
513  (approximation) in the local finite dimensional space in terms
514  of the degrees of freedom. */
515  /** The approximation used to project is usually local interpolation of
516  degrees of freedom. The derived class could use other methods not
517  implemented yet, e.g. local L2 projection. */
518  virtual void Project(Coefficient &coeff,
519  ElementTransformation &Trans, Vector &dofs) const;
520
521  /** @brief Given a vector coefficient and a transformation, compute its
522  projection (approximation) in the local finite dimensional space
523  in terms of the degrees of freedom. (VectorFiniteElements) */
524  /** The approximation used to project is usually local interpolation of
525  degrees of freedom. The derived class could use other methods not
526  implemented yet, e.g. local L2 projection. */
527  virtual void Project(VectorCoefficient &vc,
528  ElementTransformation &Trans, Vector &dofs) const;
529
530  /** @brief Given a vector of values at the finite element nodes and a
531  transformation, compute its projection (approximation) in the local
532  finite dimensional space in terms of the degrees of freedom. Valid for
533  VectorFiniteElements. */
535  Vector &dofs) const;
536
537  /** @brief Given a matrix coefficient and a transformation, compute an
538  approximation ("projection") in the local finite dimensional space in
539  terms of the degrees of freedom. For VectorFiniteElements, the rows of
540  the coefficient are projected in the vector space. */
541  virtual void ProjectMatrixCoefficient(
542  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
543
544  /** @brief Project a delta function centered on the given @a vertex in
545  the local finite dimensional space represented by the @a dofs. */
546  virtual void ProjectDelta(int vertex, Vector &dofs) const;
547
548  /** @brief Compute the embedding/projection matrix from the given
549  FiniteElement onto 'this' FiniteElement. The ElementTransformation is
550  included to support cases when the projection depends on it. */
551  virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
552  DenseMatrix &I) const;
553
554  /** @brief Compute the discrete gradient matrix from the given FiniteElement
555  onto 'this' FiniteElement. The ElementTransformation is included to
556  support cases when the matrix depends on it. */
557  virtual void ProjectGrad(const FiniteElement &fe,
560
561  /** @brief Compute the discrete curl matrix from the given FiniteElement onto
562  'this' FiniteElement. The ElementTransformation is included to support
563  cases when the matrix depends on it. */
564  virtual void ProjectCurl(const FiniteElement &fe,
566  DenseMatrix &curl) const;
567
568  /** @brief Compute the discrete divergence matrix from the given
569  FiniteElement onto 'this' FiniteElement. The ElementTransformation is
570  included to support cases when the matrix depends on it. */
571  virtual void ProjectDiv(const FiniteElement &fe,
573  DenseMatrix &div) const;
574
575  /** @brief Return a DofToQuad structure corresponding to the given
576  IntegrationRule using the given DofToQuad::Mode. */
577  /** See the documentation for DofToQuad for more details. */
580  /// Deconstruct the FiniteElement
581  virtual ~FiniteElement();
582
583  /** @brief Return true if the BasisType of @a b_type is closed
584  (has Quadrature1D points on the boundary). */
585  static bool IsClosedType(int b_type)
586  {
587  const int q_type = BasisType::GetQuadrature1D(b_type);
588  return ((q_type != Quadrature1D::Invalid) &&
590  }
591
592  /** @brief Return true if the BasisType of @a b_type is open
593  (doesn't have Quadrature1D points on the boundary). */
594  static bool IsOpenType(int b_type)
595  {
596  const int q_type = BasisType::GetQuadrature1D(b_type);
597  return ((q_type != Quadrature1D::Invalid) &&
599  }
600
601  /** @brief Ensure that the BasisType of @a b_type is closed
602  (has Quadrature1D points on the boundary). */
603  static int VerifyClosed(int b_type)
604  {
605  MFEM_VERIFY(IsClosedType(b_type),
606  "invalid closed basis type: " << b_type);
607  return b_type;
608  }
609
610  /** @brief Ensure that the BasisType of @a b_type is open
611  (doesn't have Quadrature1D points on the boundary). */
612  static int VerifyOpen(int b_type)
613  {
614  MFEM_VERIFY(IsOpenType(b_type), "invalid open basis type: " << b_type);
615  return b_type;
616  }
617
618  /** @brief Ensure that the BasisType of @a b_type nodal
619  (satisfies the interpolation property). */
620  static int VerifyNodal(int b_type)
621  {
622  return BasisType::CheckNodal(b_type);
623  }
624 };
625
626
627 /** @brief Class for finite elements with basis functions
628  that return scalar values. */
630 {
631 protected:
633  mutable Vector c_shape;
634 #endif
635
637  {
638  MFEM_VERIFY(fe.GetRangeType() == SCALAR,
639  "'fe' must be a ScalarFiniteElement");
640  return static_cast<const ScalarFiniteElement &>(fe);
641  }
642
644  const IntegrationRule &ir,
646
647 public:
648  /** @brief Construct ScalarFiniteElement with given
649  @param D Reference space dimension
650  @param G Geometry type (of type Geometry::Type)
651  @param Do Number of degrees of freedom in the FiniteElement
652  @param O Order/degree of the FiniteElement
653  @param F FunctionSpace type of the FiniteElement
654  */
656  int F = FunctionSpace::Pk)
658  : FiniteElement(D, G, Do, O, F)
660 #else
661  : FiniteElement(D, G, Do, O, F), c_shape(dof)
663 #endif
664
665  /** @brief Set the FiniteElement::MapType of the element to either VALUE or
666  INTEGRAL. Also sets the FiniteElement::DerivType to GRAD if the
667  FiniteElement::MapType is VALUE. */
668  virtual void SetMapType(int M)
669  {
670  MFEM_VERIFY(M == VALUE || M == INTEGRAL, "unknown MapType");
671  map_type = M;
672  deriv_type = (M == VALUE) ? GRAD : NONE;
673  }
674
675
676  /** @brief Get the matrix @a I that defines nodal interpolation
677  @a between this element and the refined element @a fine_fe. */
679  DenseMatrix &I,
680  const ScalarFiniteElement &fine_fe) const;
681
682  /** @brief Get matrix @a I "Interpolation" defined through local
683  L2-projection in the space defined by the @a fine_fe. */
684  /** If the "fine" elements cannot represent all basis functions of the
685  "coarse" element, then boundary values from different sub-elements are
686  generally different. */
688  DenseMatrix &I,
689  const ScalarFiniteElement &fine_fe) const;
690
691  /** @brief Get restriction matrix @a R defined through local L2-projection
692  in the space defined by the @a coarse_fe. */
693  /** If the "fine" elements cannot represent all basis functions of the
694  "coarse" element, then boundary values from different sub-elements are
695  generally different. */
697  DenseMatrix &R,
698  const ScalarFiniteElement &coarse_fe) const;
699
702 };
703
704
705 /// Class for standard nodal finite elements.
707 {
708 protected:
710  void ProjectCurl_2D(const FiniteElement &fe,
712  DenseMatrix &curl) const;
713
714 public:
715  /** @brief Construct NodalFiniteElement with given
716  @param D Reference space dimension
717  @param G Geometry type (of type Geometry::Type)
718  @param Do Number of degrees of freedom in the FiniteElement
719  @param O Order/degree of the FiniteElement
720  @param F FunctionSpace type of the FiniteElement
721  */
723  int F = FunctionSpace::Pk)
724  : ScalarFiniteElement(D, G, Do, O, F) { }
725
727  DenseMatrix &I) const
728  { NodalLocalInterpolation(Trans, I, *this); }
729
731  DenseMatrix &R) const;
732
733  virtual void GetTransferMatrix(const FiniteElement &fe,
735  DenseMatrix &I) const
736  { CheckScalarFE(fe).NodalLocalInterpolation(Trans, I, *this); }
737
738  virtual void Project (Coefficient &coeff,
739  ElementTransformation &Trans, Vector &dofs) const;
740
741  virtual void Project (VectorCoefficient &vc,
742  ElementTransformation &Trans, Vector &dofs) const;
743
744  // (mc.height x mc.width) @ DOFs -> (Dof x mc.width x mc.height) in dofs
745  virtual void ProjectMatrixCoefficient(
746  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
747
748  virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
749  DenseMatrix &I) const;
750
751  virtual void ProjectGrad(const FiniteElement &fe,
754
755  virtual void ProjectDiv(const FiniteElement &fe,
757  DenseMatrix &div) const;
758
759  /** @brief Get an Array<int> that maps lexicographically ordered indices to
760  the indices of the respective nodes/dofs/basis functions. Lexicographic
761  ordering of nodes is defined in terms of reference-space coordinates
762  (x,y,z). Lexicographically ordered nodes are listed first in order of
763  increasing x-coordinate, and then in order of increasing y-coordinate,
764  and finally in order of increasing z-coordinate.
765
766  For example, the six nodes of a quadratic triangle are lexicographically
767  ordered as follows:
768
769  5
770  |\
771  3 4
772  | \
773  0-1-2
774
775  The resulting array may be empty if the DOFs are already ordered
776  lexicographically, or if the finite element does not support creating
777  this permutation. The array returned is the same as the array given by
778  TensorBasisElement::GetDofMap, but it is also available for non-tensor
779  elements. */
781 };
782
783 /** @brief Intermediate class for finite elements whose basis functions return
784  vector values. */
786 {
787  // Hide the scalar functions CalcShape and CalcDShape.
788 private:
789  /// Overrides the scalar CalcShape function to print an error.
790  virtual void CalcShape(const IntegrationPoint &ip,
791  Vector &shape) const;
792
793  /// Overrides the scalar CalcDShape function to print an error.
794  virtual void CalcDShape(const IntegrationPoint &ip,
795  DenseMatrix &dshape) const;
796
797 protected:
798  bool is_nodal;
800  mutable DenseMatrix JtJ;
802 #endif
803  void SetDerivMembers();
804
806  DenseMatrix &shape) const;
807
809  DenseMatrix &shape) const;
810
811  /** @brief Project a vector coefficient onto the RT basis functions
812  @param nk Face normal vectors for this element type
813  @param d2n Offset into nk for each degree of freedom
814  @param vc Vector coefficient to be projected
815  @param Trans Transformation from reference to physical coordinates
816  @param dofs Expansion coefficients for the approximation of vc
817  */
818  void Project_RT(const double *nk, const Array<int> &d2n,
820  Vector &dofs) const;
821
822  /// Projects the vector of values given at FE nodes to RT space
823  /** Project vector values onto the RT basis functions
824  @param nk Face normal vectors for this element type
825  @param d2n Offset into nk for each degree of freedom
826  @param vc Vector values at each interpolation point
827  @param Trans Transformation from reference to physical coordinates
828  @param dofs Expansion coefficients for the approximation of vc
829  */
830  void Project_RT(const double *nk, const Array<int> &d2n,
832  Vector &dofs) const;
833
834  /// Project the rows of the matrix coefficient in an RT space
836  const double *nk, const Array<int> &d2n,
837  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
838
839  /** @brief Project vector-valued basis functions onto the RT basis functions
840  @param nk Face normal vectors for this element type
841  @param d2n Offset into nk for each degree of freedom
842  @param fe Vector-valued finite element basis
843  @param Trans Transformation from reference to physical coordinates
844  @param I Expansion coefficients for the approximation of each basis
845  function
846
847  Note: If the FiniteElement, fe, is scalar-valued the projection will
848  assume that a FiniteElementSpace is being used to define a vector
849  field using the scalar basis functions for each component of the
850  vector field.
851  */
852  void Project_RT(const double *nk, const Array<int> &d2n,
854  DenseMatrix &I) const;
855
856  // rotated gradient in 2D
857  void ProjectGrad_RT(const double *nk, const Array<int> &d2n,
860
861  // Compute the curl as a discrete operator from ND FE (fe) to ND FE (this).
862  // The natural FE for the range is RT, so this is an approximation.
863  void ProjectCurl_ND(const double *tk, const Array<int> &d2t,
865  DenseMatrix &curl) const;
866
867  void ProjectCurl_RT(const double *nk, const Array<int> &d2n,
869  DenseMatrix &curl) const;
870
871  /** @brief Project a vector coefficient onto the ND basis functions
872  @param tk Edge tangent vectors for this element type
873  @param d2t Offset into tk for each degree of freedom
874  @param vc Vector coefficient to be projected
875  @param Trans Transformation from reference to physical coordinates
876  @param dofs Expansion coefficients for the approximation of vc
877  */
878  void Project_ND(const double *tk, const Array<int> &d2t,
880  Vector &dofs) const;
881
882  /// Projects the vector of values given at FE nodes to ND space
883  /** Project vector values onto the ND basis functions
884  @param tk Edge tangent vectors for this element type
885  @param d2t Offset into tk for each degree of freedom
886  @param vc Vector values at each interpolation point
887  @param Trans Transformation from reference to physical coordinates
888  @param dofs Expansion coefficients for the approximation of vc
889  */
890  void Project_ND(const double *tk, const Array<int> &d2t,
892  Vector &dofs) const;
893
894  /// Project the rows of the matrix coefficient in an ND space
896  const double *tk, const Array<int> &d2t,
897  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
898
899  /** @brief Project vector-valued basis functions onto the ND basis functions
900  @param tk Edge tangent vectors for this element type
901  @param d2t Offset into tk for each degree of freedom
902  @param fe Vector-valued finite element basis
903  @param Trans Transformation from reference to physical coordinates
904  @param I Expansion coefficients for the approximation of each basis
905  function
906
907  Note: If the FiniteElement, fe, is scalar-valued the projection will
908  assume that a FiniteElementSpace is being used to define a vector
909  field using the scalar basis functions for each component of the
910  vector field.
911  */
912  void Project_ND(const double *tk, const Array<int> &d2t,
914  DenseMatrix &I) const;
915
916  void ProjectGrad_ND(const double *tk, const Array<int> &d2t,
919
922  DenseMatrix &I) const;
923
925  const double *nk, const Array<int> &d2n,
927  DenseMatrix &I) const;
928
931  DenseMatrix &I) const;
932
934  const double *tk, const Array<int> &d2t,
936  DenseMatrix &I) const;
937
938  void LocalRestriction_RT(const double *nk, const Array<int> &d2n,
940  DenseMatrix &R) const;
941
942  void LocalRestriction_ND(const double *tk, const Array<int> &d2t,
944  DenseMatrix &R) const;
945
947  {
948  if (fe.GetRangeType() != VECTOR)
949  { mfem_error("'fe' must be a VectorFiniteElement"); }
950  return static_cast<const VectorFiniteElement &>(fe);
951  }
952
953 public:
954  VectorFiniteElement (int D, Geometry::Type G, int Do, int O, int M,
955  int F = FunctionSpace::Pk);
956 };
957
958
959 /// @brief Class for computing 1D special polynomials and their associated basis
960 /// functions
961 class Poly_1D
962 {
963 public:
964  /// One-dimensional basis evaluation type
965  enum EvalType
966  {
967  ChangeOfBasis = 0, ///< Use change of basis, O(p^2) Evals
968  Barycentric = 1, ///< Use barycentric Lagrangian interpolation, O(p) Evals
969  Positive = 2, ///< Fast evaluation of Bernstein polynomials
970  Integrated = 3, ///< Integrated indicator functions (cf. Gerritsma)
971  NumEvalTypes = 4 ///< Keep count of the number of eval types
972  };
973
974  /// @brief Class for evaluating 1D nodal, positive (Bernstein), or integrated
975  /// (Gerritsma) bases.
976  class Basis
977  {
978  private:
979  EvalType etype; ///< Determines how the basis functions should be evaluated.
981  mutable Vector x, w;
982  /// The following data members are used for "integrated basis type", which
983  /// is defined in terms of nodal basis of one degree higher.
984  ///@{
985  mutable Vector u_aux, d_aux, d2_aux;
986  ///@}
987  /// @brief An auxiliary nodal basis used to evaluate the integrated basis.
988  /// This member variable is NULL whenever etype != Integrated.
989  Basis *auxiliary_basis;
990  /// Should the integrated basis functions be scaled? See ScaleIntegrated.
991  bool scale_integrated;
992
993  public:
994  /// Create a nodal or positive (Bernstein) basis of degree @a p
995  Basis(const int p, const double *nodes, EvalType etype = Barycentric);
996  /// Evaluate the basis functions at point @a x in [0,1]
997  void Eval(const double x, Vector &u) const;
998  /// @brief Evaluate the basis functions and their derivatives at point @a
999  /// x in [0,1]
1000  void Eval(const double x, Vector &u, Vector &d) const;
1001  /// @brief Evaluate the basis functions and their first two derivatives at
1002  /// point @a x in [0,1]
1003  void Eval(const double x, Vector &u, Vector &d, Vector &d2) const;
1004  /// @brief Evaluate the "integrated" basis type using pre-computed closed
1005  /// basis derivatives.
1006  ///
1007  /// This basis is given by the negative partial sum of the corresponding
1008  /// closed basis derivatives. The closed basis derivatives are given by @a
1009  /// d, and the result is stored in @a i.
1010  void EvalIntegrated(const Vector &d, Vector &i) const;
1011  /// @brief Set whether the "integrated" basis should be scaled by the
1012  /// subcell sizes. Has no effect for non-integrated bases.
1013  ///
1014  /// Generally, this should be true for mfem::FiniteElement::MapType VALUE
1015  /// and false for all other map types. If this option is enabled, the
1016  /// basis functions will be scaled by the widths of the subintervals, so
1017  /// that the basis functions represent mean values. Otherwise, the basis
1018  /// functions represent integrated values.
1019  void ScaleIntegrated(bool scale_integrated_);
1020  /// Returns true if the basis is "integrated", false otherwise.
1021  bool IsIntegratedType() const { return etype == Integrated; }
1022  ~Basis();
1023  };
1024
1025 private:
1026  typedef std::map< int, Array<double*>* > PointsMap;
1027  typedef std::map< int, Array<Basis*>* > BasisMap;
1028
1029  MemoryType h_mt;
1030  PointsMap points_container;
1031  BasisMap bases_container;
1032
1033  static Array2D<int> binom;
1034
1035  static void CalcMono(const int p, const double x, double *u);
1036  static void CalcMono(const int p, const double x, double *u, double *d);
1037
1038  static void CalcChebyshev(const int p, const double x, double *u);
1039  static void CalcChebyshev(const int p, const double x, double *u, double *d);
1040  static void CalcChebyshev(const int p, const double x, double *u, double *d,
1041  double *dd);
1042
1044
1045 public:
1046  Poly_1D(): h_mt(MemoryType::HOST) { }
1047
1048  /** @brief Get a pointer to an array containing the binomial coefficients "p
1049  choose k" for k=0,...,p for the given p. */
1050  static const int *Binom(const int p);
1051
1052  /** @brief Get the coordinates of the points of the given BasisType,
1053  @a btype.
1054
1055  @param[in] p The polynomial degree; the number of points is p+1.
1056  @param[in] btype The BasisType.
1057
1058  @return A pointer to an array containing the p+1 coordinates of the
1059  points. Returns NULL if the BasisType has no associated set of
1060  points. */
1061  const double *GetPoints(const int p, const int btype);
1062
1063  /// Get coordinates of an open (GaussLegendre) set of points if degree @a p
1064  const double *OpenPoints(const int p,
1065  const int btype = BasisType::GaussLegendre)
1066  { return GetPoints(p, btype); }
1067
1068  /// Get coordinates of a closed (GaussLegendre) set of points if degree @a p
1069  const double *ClosedPoints(const int p,
1070  const int btype = BasisType::GaussLobatto)
1071  { return GetPoints(p, btype); }
1072
1073  /** @brief Get a Poly_1D::Basis object of the given degree and BasisType,
1074  @a btype.
1075
1076  @param[in] p The polynomial degree of the basis.
1077  @param[in] btype The BasisType.
1078
1079  @return A reference to an object of type Poly_1D::Basis that represents
1080  the requested basis type. */
1081  Basis &GetBasis(const int p, const int btype);
1082
1083  /** @brief Evaluate the values of a hierarchical 1D basis at point x
1084  hierarchical = k-th basis function is degree k polynomial */
1085  static void CalcBasis(const int p, const double x, double *u)
1086  // { CalcMono(p, x, u); }
1087  // Bernstein basis is not hierarchical --> does not work for triangles
1088  // and tetrahedra
1089  // { CalcBernstein(p, x, u); }
1090  // { CalcLegendre(p, x, u); }
1091  { CalcChebyshev(p, x, u); }
1092
1093  /// Evaluate the values and derivatives of a hierarchical 1D basis at point @a x
1094  static void CalcBasis(const int p, const double x, double *u, double *d)
1095  // { CalcMono(p, x, u, d); }
1096  // { CalcBernstein(p, x, u, d); }
1097  // { CalcLegendre(p, x, u, d); }
1098  { CalcChebyshev(p, x, u, d); }
1099
1100  /// Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x
1101  static void CalcBasis(const int p, const double x, double *u, double *d,
1102  double *dd)
1103  // { CalcMono(p, x, u, d); }
1104  // { CalcBernstein(p, x, u, d); }
1105  // { CalcLegendre(p, x, u, d); }
1106  { CalcChebyshev(p, x, u, d, dd); }
1107
1108  /// Evaluate a representation of a Delta function at point x
1109  static double CalcDelta(const int p, const double x)
1110  { return pow(x, (double) p); }
1111
1112  /** @brief Compute the points for the Chebyshev polynomials of order @a p
1113  and place them in the already allocated @a x array. */
1114  static void ChebyshevPoints(const int p, double *x);
1115
1116  /** @brief Compute the @a p terms in the expansion of the binomial (x + y)^p
1117  and store them in the already allocated @a u array. */
1118  static void CalcBinomTerms(const int p, const double x, const double y,
1119  double *u);
1120  /** @brief Compute the terms in the expansion of the binomial (x + y)^p and
1121  their derivatives with respect to x assuming that dy/dx = -1. Store the
1122  results in the already allocated @a u and @a d arrays.*/
1123  static void CalcBinomTerms(const int p, const double x, const double y,
1124  double *u, double *d);
1125  /** @brief Compute the derivatives (w.r.t. x) of the terms in the expansion
1126  of the binomial (x + y)^p assuming that dy/dx = -1. Store the results
1127  in the already allocated @a d array.*/
1128  static void CalcDBinomTerms(const int p, const double x, const double y,
1129  double *d);
1130
1131  /** @brief Compute the values of the Bernstein basis functions of order
1132  @a p at coordinate @a x and store the results in the already allocated
1133  @a u array. */
1134  static void CalcBernstein(const int p, const double x, double *u)
1135  { CalcBinomTerms(p, x, 1. - x, u); }
1136
1137  /** @brief Compute the values and derivatives of the Bernstein basis functions
1138  of order @a p at coordinate @a x and store the results in the already allocated
1139  @a u and @a d arrays. */
1140  static void CalcBernstein(const int p, const double x, double *u, double *d)
1141  { CalcBinomTerms(p, x, 1. - x, u, d); }
1142
1143  static void CalcLegendre(const int p, const double x, double *u);
1144  static void CalcLegendre(const int p, const double x, double *u, double *d);
1145
1146  ~Poly_1D();
1147 };
1148
1149 extern Poly_1D poly1d;
1150
1151
1152 /// An element defined as an ND tensor product of 1D elements on a segment,
1153 /// square, or cube
1155 {
1156 protected:
1157  int b_type;
1161
1162 public:
1164  {
1167  Sr_DOF_MAP = 2, // Sr = Serendipity
1168  };
1169
1170  TensorBasisElement(const int dims, const int p, const int btype,
1171  const DofMapType dmtype);
1172
1173  int GetBasisType() const { return b_type; }
1174
1175  const Poly_1D::Basis& GetBasis1D() const { return basis1d; }
1176
1177  /** @brief Get an Array<int> that maps lexicographically ordered indices to
1178  the indices of the respective nodes/dofs/basis functions. If the dofs are
1179  ordered lexicographically, i.e. the mapping is identity, the returned
1180  Array will be empty. */
1181  const Array<int> &GetDofMap() const { return dof_map; }
1182
1184  {
1185  switch (dim)
1186  {
1187  case 1: return Geometry::SEGMENT;
1188  case 2: return Geometry::SQUARE;
1189  case 3: return Geometry::CUBE;
1190  default:
1191  MFEM_ABORT("invalid dimension: " << dim);
1192  return Geometry::INVALID;
1193  }
1194  }
1195
1196  /// Return @a base raised to the power @a dim.
1197  static int Pow(int base, int dim)
1198  {
1199  switch (dim)
1200  {
1201  case 1: return base;
1202  case 2: return base*base;
1203  case 3: return base*base*base;
1204  default: MFEM_ABORT("invalid dimension: " << dim); return -1;
1205  }
1206  }
1207 };
1208
1210  public TensorBasisElement
1211 {
1212 public:
1213  NodalTensorFiniteElement(const int dims, const int p, const int btype,
1214  const DofMapType dmtype);
1215
1218  {
1219  return (mode == DofToQuad::FULL) ?
1222  }
1223
1224  virtual void SetMapType(const int map_type_);
1225
1226  virtual void GetTransferMatrix(const FiniteElement &fe,
1228  DenseMatrix &I) const
1229  {
1230  if (basis1d.IsIntegratedType())
1231  {
1232  CheckScalarFE(fe).ScalarLocalInterpolation(Trans, I, *this);
1233  }
1234  else
1235  {
1237  }
1238  }
1239 };
1240
1242  public TensorBasisElement
1243 {
1244 private:
1246
1247 protected:
1249
1250 public:
1251  VectorTensorFiniteElement(const int dims, const int d, const int p,
1252  const int cbtype, const int obtype,
1253  const int M, const DofMapType dmtype);
1254
1255  // For 1D elements: there is only an "open basis", no "closed basis"
1256  VectorTensorFiniteElement(const int dims, const int d, const int p,
1257  const int obtype, const int M,
1258  const DofMapType dmtype);
1259
1262
1265
1268  const bool closed) const;
1269
1271 };
1272
1274  const IntegrationPoint &pt, Vector &x);
1275
1276 } // namespace mfem
1277
1278 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:235
Integrated GLL indicator functions.
Definition: fe_base.hpp:38
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
ScalarFiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct ScalarFiniteElement with given.
Definition: fe_base.hpp:655
static void CalcBernstein(const int p, const double x, double *u, double *d)
Compute the values and derivatives of the Bernstein basis functions of order p at coordinate x and st...
Definition: fe_base.hpp:1140
Array< int > lex_ordering
Definition: fe_base.hpp:709
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input...
Definition: fe_base.hpp:44
virtual void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Definition: fe_base.cpp:244
static double CalcDelta(const int p, const double x)
Evaluate a representation of a Delta function at point x.
Definition: fe_base.hpp:1109
void ProjectMatrixCoefficient_RT(const double *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an RT space.
Definition: fe_base.cpp:963
EvalType
One-dimensional basis evaluation type.
Definition: fe_base.hpp:965
void ProjectMatrixCoefficient_ND(const double *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an ND space.
Definition: fe_base.cpp:1180
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_base.cpp:39
void ScaleIntegrated(bool scale_integrated_)
Set whether the &quot;integrated&quot; basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition: fe_base.cpp:1879
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:2464
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition: fe_base.cpp:2407
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1347
Basis(const int p, const double *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis of degree p.
Definition: fe_base.cpp:1558
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:162
bool IsIntegratedType() const
Returns true if the basis is &quot;integrated&quot;, false otherwise.
Definition: fe_base.hpp:1021
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:174
Base class for vector Coefficients that optionally depend on time and space.
static int GetNodalBasis(int qpt_type)
Return the nodal BasisType corresponding to the Quadrature1D type.
Definition: fe_base.hpp:77
virtual void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in reference space at the give...
Definition: fe_base.cpp:289
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:125
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
Definition: fe_base.cpp:57
int GetDerivMapType() const
Returns the FiniteElement::DerivType of the element describing how reference function derivatives are...
Definition: fe_base.hpp:360
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation (&quot;projection&quot;) in the local...
Definition: fe_base.cpp:695
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
Definition: fe_base.cpp:903
virtual void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_base.cpp:202
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:656
FiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct FiniteElement with given.
Definition: fe_base.cpp:22
int dim
Dimension of reference space.
Definition: fe_base.hpp:238
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition: fe_base.cpp:2206
aka closed Newton-Cotes
Definition: intrules.hpp:298
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:174
Array< double > Gt
Transpose of G.
Definition: fe_base.hpp:213
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:776
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:170
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition: fe_base.hpp:60
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
Definition: fe_base.cpp:837
MapType
Enumeration for MapType: defines how reference functions are mapped to physical space.
Definition: fe_base.hpp:271
void ProjectCurl_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:1093
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:337
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition: fe_base.cpp:106
virtual void SetMapType(const int map_type_)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.cpp:2418
DenseMatrix vshape
Definition: fe_base.hpp:250
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
Nodes: x_i = (i+1/2)/n, i=0,...,n-1.
Definition: fe_base.hpp:35
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition: fe_base.hpp:241
aka open Newton-Cotes
Definition: intrules.hpp:297
DerivType
Enumeration for DerivType: defines which derivative method is implemented.
Definition: fe_base.hpp:292
Intermediate class for finite elements whose basis functions return vector values.
Definition: fe_base.hpp:785
static void CalcLegendre(const int p, const double x, double *u)
Definition: fe_base.cpp:2032
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
Definition: fe_base.cpp:382
static const char * Name(int b_type)
Check and convert a BasisType constant to a string identifier.
Definition: fe_base.hpp:91
void ProjectCurl_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:1131
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:160
Use change of basis, O(p^2) Evals.
Definition: fe_base.hpp:967
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.hpp:733
Definition: fe_base.cpp:2482
static int GetType(char b_ident)
Convert char basis identifier to a BasisType constant.
Definition: fe_base.hpp:110
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:1069
static int CheckNodal(int b_type)
If the input does not represents a valid nodal BasisType, abort with an error; otherwise return the i...
Definition: fe_base.hpp:52
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition: fe_base.hpp:317
Tensor products of polynomials of order k.
Definition: fe_base.hpp:224
Possible basis types. Note that not all elements can use all BasisType(s).
Definition: fe_base.hpp:24
Serendipity basis (squares / cubes)
Definition: fe_base.hpp:36
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition: fe_base.cpp:70
int cdim
Dimension of curl for vector-valued basis functions.
Definition: fe_base.hpp:240
int vdim
Vector dimension of vector-valued basis functions.
Definition: fe_base.hpp:239
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:1085
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
Definition: fe_base.cpp:137
static void CalcBinomTerms(const int p, const double x, const double y, double *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
Definition: fe_base.cpp:1938
const double * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition: fe_base.cpp:2131
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:167
Class for standard nodal finite elements.
Definition: fe_base.hpp:706
void ScalarLocalRestriction(ElementTransformation &Trans, DenseMatrix &R, const ScalarFiniteElement &coarse_fe) const
Get restriction matrix R defined through local L2-projection in the space defined by the coarse_fe...
Definition: fe_base.cpp:459
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
VectorTensorFiniteElement(const int dims, const int d, const int p, const int cbtype, const int obtype, const int M, const DofMapType dmtype)
Definition: fe_base.cpp:2431
Class for finite elements with basis functions that return scalar values.
Definition: fe_base.hpp:629
static const ScalarFiniteElement & CheckScalarFE(const FiniteElement &fe)
Definition: fe_base.hpp:636
int GetBasisType() const
Definition: fe_base.hpp:1173
virtual void SetMapType(int M)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.hpp:668
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.hpp:1216
Poly_1D::Basis & basis1d
Definition: fe_base.hpp:1159
Fast evaluation of Bernstein polynomials.
Definition: fe_base.hpp:969
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
static void ChebyshevPoints(const int p, double *x)
Compute the points for the Chebyshev polynomials of order p and place them in the already allocated x...
Definition: fe_base.cpp:1906
Class for computing 1D special polynomials and their associated basis functions.
Definition: fe_base.hpp:961
const int * GetAnisotropicOrders() const
Returns an array containing the anisotropic orders/degrees.
Definition: fe_base.hpp:334
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition: fe_base.hpp:166
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...
class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition: fe_base.hpp:141
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
RangeType
Enumeration for range_type and deriv_range_type.
Definition: fe_base.hpp:259
static int VerifyOpen(int b_type)
Ensure that the BasisType of b_type is open (doesn&#39;t have Quadrature1D points on the boundary)...
Definition: fe_base.hpp:612
int GetDerivRangeType() const
Returns the FiniteElement::RangeType of the element derivative, either SCALAR or VECTOR.
Definition: fe_base.hpp:344
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the &quot;integrated&quot; basis type using pre-computed closed basis derivatives. ...
Definition: fe_base.cpp:1854
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition: fe_base.hpp:247
virtual ~FiniteElement()
Deconstruct the FiniteElement.
Definition: fe_base.cpp:373
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe_base.cpp:191
Use barycentric Lagrangian interpolation, O(p) Evals.
Definition: fe_base.hpp:968
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition: fe_base.cpp:2155
const double * OpenPoints(const int p, const int btype=BasisType::GaussLegendre)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition: fe_base.hpp:1064
static void CalcBasis(const int p, const double x, double *u, double *d, double *dd)
Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x...
Definition: fe_base.hpp:1101
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I &quot;Interpolation&quot; defined through local L2-projection in the space defined by the fine_fe...
Definition: fe_base.cpp:418
void ProjectGrad_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition: fe_base.cpp:1279
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition: fe_base.hpp:412
static const int MaxDim
Definition: geom.hpp:43
bool HasAnisotropicOrders() const
Returns true if the FiniteElement basis may be using different orders/degrees in different spatial di...
Definition: fe_base.hpp:331
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:314
IntegrationRule Nodes
Definition: fe_base.hpp:248
Array< double > Bt
Transpose of B.
Definition: fe_base.hpp:191
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1434
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:854
Integrated indicator functions (cf. Gerritsma)
Definition: fe_base.hpp:970
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
Definition: fe_base.cpp:915
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:31
static bool IsClosedType(int b_type)
Return true if the BasisType of b_type is closed (has Quadrature1D points on the boundary).
Definition: fe_base.hpp:585
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_base.cpp:100
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition: fe_base.cpp:181
Base class for Matrix Coefficients that optionally depend on time and space.
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients &quot;p choose k&quot; for k=0,...,p for the given p.
Definition: fe_base.cpp:1889
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition: fe_base.cpp:51
Closed GaussLegendre.
Definition: fe_base.hpp:37
Array< int > inv_dof_map
Definition: fe_base.hpp:1160
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:136
const Array< int > & GetLexicographicOrdering() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:780
Polynomials of order k.
Definition: fe_base.hpp:223
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:866
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:185
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_base.hpp:1181
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
Definition: fe_base.cpp:95
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:366
static char GetChar(int b_type)
Check and convert a BasisType constant to a char basis identifier.
Definition: fe_base.hpp:103
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_base.cpp:64
void LocalRestriction_ND(const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
Definition: fe_base.cpp:1516
Class for integration point with weight.
Definition: intrules.hpp:25
Host memory; using new[] and delete[].
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition: fe_base.hpp:149
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation (&quot;projection&quot;) in the local...
Definition: fe_base.cpp:143
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition: fe_base.hpp:154
const Poly_1D::Basis & GetBasis1D() const
Definition: fe_base.hpp:1175
static bool IsOpenType(int b_type)
Return true if the BasisType of b_type is open (doesn&#39;t have Quadrature1D points on the boundary)...
Definition: fe_base.hpp:594
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition: fe_base.hpp:726
Definition: fe_base.cpp:2473
void Project_RT(const double *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the RT basis functions.
Definition: fe_base.cpp:926
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:245
Container for all DofToQuad objects created by the FiniteElement.
Definition: fe_base.hpp:255
Nodes: x_i = (i+1)/(n+1), i=0,...,n-1.
Definition: fe_base.hpp:33
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:603
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:805
Implements CalcDivShape methods.
Definition: fe_base.hpp:296
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Definition: fe_base.cpp:112
int dim
Definition: ex24.cpp:53
NodalFiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct NodalFiniteElement with given.
Definition: fe_base.hpp:722
void ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:587
aka &quot;open half&quot; Newton-Cotes
Definition: intrules.hpp:299
Array< double > G
Definition: fe_base.hpp:206
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition: fe_base.hpp:976
static int VerifyNodal(int b_type)
Ensure that the BasisType of b_type nodal (satisfies the interpolation property). ...
Definition: fe_base.hpp:620
Bernstein polynomials.
Definition: fe_base.hpp:32
static int Pow(int base, int dim)
Return base raised to the power dim.
Definition: fe_base.hpp:1197
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1622
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1388
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition: fe_base.cpp:611
Vector data type.
Definition: vector.hpp:60
static void CalcBernstein(const int p, const double x, double *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition: fe_base.hpp:1134
int GetDerivType() const
Returns the FiniteElement::DerivType of the element describing the spatial derivative method implemen...
Definition: fe_base.hpp:355
static Geometry::Type GetTensorProductGeometry(int dim)
Definition: fe_base.hpp:1183
Refined tensor products of polynomials of order k.
Definition: fe_base.hpp:225
Nodes: x_i = i/(n-1), i=0,...,n-1.
Definition: fe_base.hpp:34
void LocalRestriction_RT(const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
Definition: fe_base.cpp:1473
Describes the function space on each element.
Definition: fe_base.hpp:218
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.cpp:118
static const VectorFiniteElement & CheckVectorFE(const FiniteElement &fe)
Definition: fe_base.hpp:946
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Implements CalcDShape methods.
Definition: fe_base.hpp:295
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:497
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...
void Project_ND(const double *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the ND basis functions.
Definition: fe_base.cpp:1150
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:340
static void CalcBasis(const int p, const double x, double *u, double *d)
Evaluate the values and derivatives of a hierarchical 1D basis at point x.
Definition: fe_base.hpp:1094
A Class that defines 1-D numerical quadrature rules on [0,1].
Definition: intrules.hpp:265
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Definition: fe_base.cpp:627
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.hpp:1226
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
Definition: fe_base.cpp:2002
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Definition: fe_base.hpp:146
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_base.cpp:149
Definition: fe_base.cpp:545
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:245
aka closed Gauss Legendre
Definition: intrules.hpp:300
Keep count of the number of eval types.
Definition: fe_base.hpp:971
Poly_1D poly1d
Definition: fe.cpp:28
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1300
No derivatives implemented.
Definition: fe_base.hpp:294
Implements CalcCurlShape methods.
Definition: fe_base.hpp:297
void ProjectGrad_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition: fe_base.cpp:1066