fe.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2019 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_FE_H
#define LIBMESH_FE_H
 
// Local includes
#include "libmesh/fe_base.h"
#include "libmesh/libmesh.h"
 
// C++ includes
#include <cstddef>
 
namespace libMesh
{
 
// forward declarations
class DofConstraints;
class DofMap;
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
class InfFE;
 
#endif
 
 
template <FEFamily T>
struct FEOutputType
{
  typedef Real type;
};
 
 
template<>
struct FEOutputType<LAGRANGE_VEC>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<NEDELEC_ONE>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<MONOMIAL_VEC>
{
  typedef RealVectorValue type;
};
 
 
template <unsigned int Dim, FEFamily T>
class FE : public FEGenericBase<typename FEOutputType<T>::type>
{
public:
 
  explicit
  FE(const FEType & fet);
 
  typedef typename
  FEGenericBase<typename FEOutputType<T>::type>::OutputShape
  OutputShape;
 
  static OutputShape shape(const ElemType t,
                           const Order o,
                           const unsigned int i,
                           const Point & p);
 
  static OutputShape shape(const Elem * elem,
                           const Order o,
                           const unsigned int i,
                           const Point & p,
                           const bool add_p_level = true);
 
  static OutputShape shape_deriv(const ElemType t,
                                 const Order o,
                                 const unsigned int i,
                                 const unsigned int j,
                                 const Point & p);
 
  static OutputShape shape_deriv(const Elem * elem,
                                 const Order o,
                                 const unsigned int i,
                                 const unsigned int j,
                                 const Point & p,
                                 const bool add_p_level = true);
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  static OutputShape shape_second_deriv(const ElemType t,
                                        const Order o,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p);
 
  static OutputShape shape_second_deriv(const Elem * elem,
                                        const Order o,
                                        const unsigned int i,
                                        const unsigned int j,
                                        const Point & p,
                                        const bool add_p_level = true);
 
#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  static void nodal_soln(const Elem * elem, const Order o,
                         const std::vector<Number> & elem_soln,
                         std::vector<Number> & nodal_soln);
 
  virtual unsigned int n_shape_functions () const override;
 
  static unsigned int n_shape_functions (const ElemType t,
                                         const Order o)
  { return FE<Dim,T>::n_dofs (t,o); }
 
  static unsigned int n_dofs(const ElemType t,
                             const Order o);
 
  static unsigned int n_dofs_at_node(const ElemType t,
                                     const Order o,
                                     const unsigned int n);
 
  static unsigned int n_dofs_per_elem(const ElemType t,
                                      const Order o);
 
  virtual FEContinuity get_continuity() const override;
 
  virtual bool is_hierarchic() const override;
 
  static void dofs_on_side(const Elem * const elem,
                           const Order o,
                           unsigned int s,
                           std::vector<unsigned int> & di);
  static void dofs_on_edge(const Elem * const elem,
                           const Order o,
                           unsigned int e,
                           std::vector<unsigned int> & di);
 
  static Point inverse_map (const Elem * elem,
                            const Point & p,
                            const Real tolerance = TOLERANCE,
                            const bool secure = true) {
    // libmesh_deprecated(); // soon
    return FEMap::inverse_map(Dim, elem, p, tolerance, secure);
  }
 
  static void inverse_map (const Elem * elem,
                           const std::vector<Point> & physical_points,
                           std::vector<Point> &       reference_points,
                           const Real tolerance = TOLERANCE,
                           const bool secure = true) {
    // libmesh_deprecated(); // soon
    FEMap::inverse_map(Dim, elem, physical_points, reference_points,
                       tolerance, secure);
  }
 
  virtual void reinit (const Elem * elem,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) override;
 
  virtual void reinit (const Elem * elem,
                       const unsigned int side,
                       const Real tolerance = TOLERANCE,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) override;
 
  virtual void edge_reinit (const Elem * elem,
                            const unsigned int edge,
                            const Real tolerance = TOLERANCE,
                            const std::vector<Point> * const pts = nullptr,
                            const std::vector<Real> * const weights = nullptr) override;
 
  virtual void side_map (const Elem * elem,
                         const Elem * side,
                         const unsigned int s,
                         const std::vector<Point> & reference_side_points,
                         std::vector<Point> &       reference_points) override;
 
  virtual void attach_quadrature_rule (QBase * q) override;
 
  virtual unsigned int n_quadrature_points () const override;
 
#ifdef LIBMESH_ENABLE_AMR
 
  static void compute_constraints (DofConstraints & constraints,
                                   DofMap & dof_map,
                                   const unsigned int variable_number,
                                   const Elem * elem);
#endif // #ifdef LIBMESH_ENABLE_AMR
 
  virtual bool shapes_need_reinit() const override;
 
  static Point map (const Elem * elem,
                    const Point & reference_point) {
    // libmesh_deprecated(); // soon
    return FEMap::map(Dim, elem, reference_point);
  }
 
  static Point map_xi (const Elem * elem,
                       const Point & reference_point) {
    // libmesh_deprecated(); // soon
    return FEMap::map_deriv(Dim, elem, 0, reference_point);
  }
 
  static Point map_eta (const Elem * elem,
                        const Point & reference_point) {
    // libmesh_deprecated(); // soon
    return FEMap::map_deriv(Dim, elem, 1, reference_point);
  }
 
  static Point map_zeta (const Elem * elem,
                         const Point & reference_point) {
    // libmesh_deprecated(); // soon
    return FEMap::map_deriv(Dim, elem, 2, reference_point);
  }
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
  friend class InfFE;
#endif
 
protected:
 
  virtual void init_shape_functions(const std::vector<Point> & qp,
                                    const Elem * e);
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  virtual void init_base_shape_functions(const std::vector<Point> & qp,
                                         const Elem * e) override;
 
#endif
 
  std::vector<Point> cached_nodes;
 
  ElemType last_side;
 
  unsigned int last_edge;
};
 
 
 
template <unsigned int Dim>
class FEClough : public FE<Dim,CLOUGH>
{
public:
 
  explicit
  FEClough(const FEType & fet) :
    FE<Dim,CLOUGH> (fet)
  {}
};
 
 
 
template <unsigned int Dim>
class FEHermite : public FE<Dim,HERMITE>
{
public:
 
  explicit
  FEHermite(const FEType & fet) :
    FE<Dim,HERMITE> (fet)
  {}
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  static Real hermite_raw_shape_second_deriv(const unsigned int basis_num,
                                             const Real xi);
#endif
  static Real hermite_raw_shape_deriv(const unsigned int basis_num,
                                      const Real xi);
  static Real hermite_raw_shape(const unsigned int basis_num,
                                const Real xi);
};
 
 
 
template <>
Real FE<2,SUBDIVISION>::shape(const Elem * elem,
                              const Order order,
                              const unsigned int i,
                              const Point & p,
                              const bool add_p_level);
 
template <>
Real FE<2,SUBDIVISION>::shape_deriv(const Elem * elem,
                                    const Order order,
                                    const unsigned int i,
                                    const unsigned int j,
                                    const Point & p,
                                    const bool add_p_level);
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
template <>
Real FE<2,SUBDIVISION>::shape_second_deriv(const Elem * elem,
                                           const Order order,
                                           const unsigned int i,
                                           const unsigned int j,
                                           const Point & p,
                                           const bool add_p_level);
 
#endif
 
class FESubdivision : public FE<2,SUBDIVISION>
{
public:
 
  FESubdivision(const FEType & fet);
 
  virtual void reinit (const Elem * elem,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) override;
 
  virtual void reinit (const Elem *,
                       const unsigned int,
                       const Real = TOLERANCE,
                       const std::vector<Point> * const = nullptr,
                       const std::vector<Real> * const = nullptr) override
  { libmesh_not_implemented(); }
 
  virtual void attach_quadrature_rule (QBase * q) override;
 
  virtual void init_shape_functions(const std::vector<Point> & qp,
                                    const Elem * elem) override;
 
  static Real regular_shape(const unsigned int i,
                            const Real v,
                            const Real w);
 
  static Real regular_shape_deriv(const unsigned int i,
                                  const unsigned int j,
                                  const Real v,
                                  const Real w);
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  static Real regular_shape_second_deriv(const unsigned int i,
                                         const unsigned int j,
                                         const Real v,
                                         const Real w);
 
 
#endif // LIBMESH_ENABLE_SECOND_DERIVATIVE
 
  static void loop_subdivision_mask(std::vector<Real> & weights,
                                    const unsigned int valence);
 
 
  static void init_subdivision_matrix(DenseMatrix<Real> & A,
                                      unsigned int valence);
};
 
 
 
template <unsigned int Dim>
class FEHierarchic : public FE<Dim,HIERARCHIC>
{
public:
 
  explicit
  FEHierarchic(const FEType & fet) :
    FE<Dim,HIERARCHIC> (fet)
  {}
};
 
 
 
template <unsigned int Dim>
class FEL2Hierarchic : public FE<Dim,L2_HIERARCHIC>
{
public:
 
  explicit
  FEL2Hierarchic(const FEType & fet) :
    FE<Dim,L2_HIERARCHIC> (fet)
  {}
};
 
 
 
template <unsigned int Dim>
class FELagrange : public FE<Dim,LAGRANGE>
{
public:
 
  explicit
  FELagrange(const FEType & fet) :
    FE<Dim,LAGRANGE> (fet)
  {}
};
 
 
template <unsigned int Dim>
class FEL2Lagrange : public FE<Dim,L2_LAGRANGE>
{
public:
 
  explicit
  FEL2Lagrange(const FEType & fet) :
    FE<Dim,L2_LAGRANGE> (fet)
  {}
};
 
 
template <unsigned int Dim>
class FEMonomial : public FE<Dim,MONOMIAL>
{
public:
 
  explicit
  FEMonomial(const FEType & fet) :
    FE<Dim,MONOMIAL> (fet)
  {}
};
 
 
template <unsigned int Dim>
class FEScalar : public FE<Dim,SCALAR>
{
public:
 
  explicit
  FEScalar(const FEType & fet) :
    FE<Dim,SCALAR> (fet)
  {}
};
 
 
template <unsigned int Dim>
class FEXYZ : public FE<Dim,XYZ>
{
public:
 
  explicit
  FEXYZ(const FEType & fet) :
    FE<Dim,XYZ> (fet)
  {}
 
  virtual void reinit (const Elem * elem,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) override
  { FE<Dim,XYZ>::reinit (elem, pts, weights); }
 
  virtual void reinit (const Elem * elem,
                       const unsigned int side,
                       const Real tolerance = TOLERANCE,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) override;
 
 
protected:
 
  virtual void init_shape_functions(const std::vector<Point> & qp,
                                    const Elem * e) override;
 
  virtual void compute_shape_functions(const Elem * elem, const std::vector<Point> & qp) override;
 
  void compute_face_values (const Elem * elem,
                            const Elem * side,
                            const std::vector<Real> & weights);
};
 
 
 
template <unsigned int Dim>
class FELagrangeVec : public FE<Dim,LAGRANGE_VEC>
{
public:
 
  explicit
  FELagrangeVec(const FEType & fet) :
    FE<Dim,LAGRANGE_VEC> (fet)
  {}
};
 
 
 
template <unsigned int Dim>
class FENedelecOne : public FE<Dim,NEDELEC_ONE>
{
public:
  explicit
  FENedelecOne(const FEType & fet) :
    FE<Dim,NEDELEC_ONE> (fet)
  {}
};
 
template <unsigned int Dim>
class FEMonomialVec : public FE<Dim,MONOMIAL_VEC>
{
public:
 
  explicit
  FEMonomialVec(const FEType & fet) :
    FE<Dim,MONOMIAL_VEC> (fet)
  {}
};
 
 
 
namespace FiniteElements
{
typedef FEClough<2> FEClough2D;
 
typedef FE<1,HIERARCHIC> FEHierarchic1D;
 
typedef FE<2,HIERARCHIC> FEHierarchic2D;
 
typedef FE<3,HIERARCHIC> FEHierarchic3D;
 
 
typedef FE<1,L2_HIERARCHIC> FEL2Hierarchic1D;
 
typedef FE<2,L2_HIERARCHIC> FEL2Hierarchic2D;
 
typedef FE<3,L2_HIERARCHIC> FEL2Hierarchic3D;
 
 
typedef FE<1,LAGRANGE> FELagrange1D;
 
typedef FE<2,LAGRANGE> FELagrange2D;
 
typedef FE<3,LAGRANGE> FELagrange3D;
 
 
typedef FE<1,L2_LAGRANGE> FEL2Lagrange1D;
 
typedef FE<2,L2_LAGRANGE> FEL2Lagrange2D;
 
typedef FE<3,L2_LAGRANGE> FEL2Lagrange3D;
 
 
typedef FE<1,MONOMIAL> FEMonomial1D;
 
typedef FE<2,MONOMIAL> FEMonomial2D;
 
typedef FE<3,MONOMIAL> FEMonomial3D;
 
}
 
} // namespace libMesh
 
#endif // LIBMESH_FE_H