fe.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 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/int_range.h"
#include "libmesh/libmesh.h"
 
// C++ includes
#include <cstddef>
 
namespace libMesh
{
 
// forward declarations
class DofConstraints;
class DofMap;
class QGauss;
 
#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<L2_LAGRANGE_VEC>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<HIERARCHIC_VEC>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<L2_HIERARCHIC_VEC>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<NEDELEC_ONE>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<MONOMIAL_VEC>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<RAVIART_THOMAS>
{
  typedef RealVectorValue type;
};
 
template<>
struct FEOutputType<L2_RAVIART_THOMAS>
{
  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(const FEType fet,
                            const Elem * elem,
                            const unsigned int i,
                            const Point & p,
                            const bool add_p_level = true);
 
  static void shapes(const Elem * elem,
                     const Order o,
                     const unsigned int i,
                     const std::vector<Point> & p,
                     std::vector<OutputShape> & v,
                     const bool add_p_level = true);
 
 
  static void all_shapes(const Elem * elem,
                         const Order o,
                         const std::vector<Point> & p,
                         std::vector<std::vector<OutputShape> > & v,
                         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);
 
  static OutputShape shape_deriv(const FEType fet,
                                 const Elem * elem,
                                 const unsigned int i,
                                 const unsigned int j,
                                 const Point & p,
                                 const bool add_p_level = true);
 
  static void shape_derivs(const Elem * elem,
                           const Order o,
                           const unsigned int i,
                           const unsigned int j,
                           const std::vector<Point> & p,
                           std::vector<OutputShape> & v,
                           const bool add_p_level = true);
 
  static void all_shape_derivs(const Elem * elem,
                               const Order o,
                               const std::vector<Point> & p,
                               std::vector<std::vector<OutputShape>> * comps[3],
                               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);
 
  static OutputShape shape_second_deriv(const FEType fet,
                                        const Elem * elem,
                                        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,
                         bool add_p_level = true,
                         const unsigned vdim = 1);
 
  static void side_nodal_soln(const Elem * elem, const Order o,
                              const unsigned int side,
                              const std::vector<Number> & elem_soln,
                              std::vector<Number> & nodal_soln_on_side,
                              bool add_p_level = true,
                              const unsigned vdim = 1);
 
  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(const Elem * e,
                             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_at_node(const Elem & e,
                                     const Order o,
                                     const unsigned int n);
 
  static unsigned int n_dofs_per_elem(const ElemType t,
                                      const Order o);
 
  static unsigned int n_dofs_per_elem(const Elem & e,
                                      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,
                           bool add_p_level=true);
  static void dofs_on_edge(const Elem * const elem,
                           const Order o,
                           unsigned int e,
                           std::vector<unsigned int> & di,
                           bool add_p_level=true);
 
  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, 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, 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_dual_shape_coeffs (const Elem * elem,
                                         const std::vector<Point> & pts,
                                         const std::vector<Real> & JxW) override;
 
   virtual void reinit_default_dual_shape_coeffs (const Elem * elem) 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 edge_map (const Elem * elem,
                         const Elem * edge,
                         const unsigned int e,
                         const std::vector<Point> & reference_edge_points,
                         std::vector<Point> &       reference_points);
 
  virtual void attach_quadrature_rule (QBase * q) 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);
 
  static void default_all_shape_derivs (const Elem * elem,
                                        const Order o,
                                        const std::vector<Point> & p,
                                        std::vector<std::vector<OutputShape>> * comps[3],
                                        const bool add_p_level = true);
 
 
  void init_dual_shape_functions(unsigned int n_shapes, unsigned int n_qp);
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  virtual void init_base_shape_functions(const std::vector<Point> & qp,
                                         const Elem * e) override;
 
#endif
 
  static void default_shapes (const Elem * elem,
                              const Order o,
                              const unsigned int i,
                              const std::vector<Point> & p,
                              std::vector<OutputShape> & v,
                              const bool add_p_level = true)
    {
      libmesh_assert_equal_to(p.size(), v.size());
      for (auto vi : index_range(v))
        v[vi] = FE<Dim,T>::shape (elem, o, i, p[vi], add_p_level);
    }
 
  static void default_all_shapes (const Elem * elem,
                                  const Order o,
                                  const std::vector<Point> & p,
                                  std::vector<std::vector<OutputShape>> & v,
                                  const bool add_p_level = true)
    {
      for (auto i : index_range(v))
        {
          libmesh_assert_equal_to ( p.size(), v[i].size() );
          FE<Dim,T>::shapes (elem, o, i, p, v[i], add_p_level);
        }
    }
 
  static void default_shape_derivs (const Elem * elem,
                                    const Order o,
                                    const unsigned int i,
                                    const unsigned int j,
                                    const std::vector<Point> & p,
                                    std::vector<OutputShape> & v,
                                    const bool add_p_level = true)
    {
      libmesh_assert_equal_to(p.size(), v.size());
      for (auto vi : index_range(v))
        v[vi] = FE<Dim,T>::shape_deriv (elem, o, i, j, p[vi], add_p_level);
    }
 
  static void default_side_nodal_soln(const Elem * elem, const Order o,
                                      const unsigned int side,
                                      const std::vector<Number> & elem_soln,
                                      std::vector<Number> & nodal_soln_on_side,
                                      bool add_p_level = true,
                                      const unsigned vdim = 1);
 
  std::vector<Point> cached_nodes;
  std::vector<bool> cached_edges, cached_faces;
 
  void cache(const Elem * elem);
 
  bool matches_cache(const Elem * elem);
 
  ElemType last_side;
 
  ElemType 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 FEL2LagrangeVec : public FE<Dim,L2_LAGRANGE_VEC>
{
public:
 
  explicit
  FEL2LagrangeVec(const FEType & fet) :
    FE<Dim,L2_LAGRANGE_VEC> (fet)
  {}
};
 
 
 
template <unsigned int Dim>
class FEHierarchicVec : public FE<Dim,HIERARCHIC_VEC>
{
public:
 
  explicit
  FEHierarchicVec(const FEType & fet) :
    FE<Dim,HIERARCHIC_VEC> (fet)
  {}
};
 
 
template <unsigned int Dim>
class FEL2HierarchicVec : public FE<Dim,L2_HIERARCHIC_VEC>
{
public:
 
  explicit
  FEL2HierarchicVec(const FEType & fet) :
    FE<Dim,L2_HIERARCHIC_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)
  {}
};
 
template <unsigned int Dim>
class FERaviartThomas : public FE<Dim,RAVIART_THOMAS>
{
public:
 
  explicit
  FERaviartThomas(const FEType & fet) :
    FE<Dim,RAVIART_THOMAS> (fet)
  {}
};
 
template <unsigned int Dim>
class FEL2RaviartThomas : public FE<Dim,L2_RAVIART_THOMAS>
{
public:
 
  explicit
  FEL2RaviartThomas(const FEType & fet) :
    FE<Dim,L2_RAVIART_THOMAS> (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;
 
}
 
template <typename OutputShape>
OutputShape fe_fdm_deriv(const Elem * elem,
                         const Order order,
                         const unsigned int i,
                         const unsigned int j,
                         const Point & p,
                         const bool add_p_level,
                         OutputShape(*shape_func)
                           (const Elem *, const Order,
                            const unsigned int, const Point &,
                            const bool));
 
template <typename OutputShape>
OutputShape fe_fdm_deriv(const ElemType type,
                         const Order order,
                         const unsigned int i,
                         const unsigned int j,
                         const Point & p,
                         OutputShape(*shape_func)
                           (const ElemType, const Order,
                            const unsigned int, const Point &));
 
template <typename OutputShape>
OutputShape fe_fdm_deriv(const ElemType type,
                         const Order order,
                         const Elem * elem,
                         const unsigned int i,
                         const unsigned int j,
                         const Point & p,
                         OutputShape(*shape_func)
                           (const ElemType type, const Order,
                            const Elem *, const unsigned int,
                            const Point &));
 
 
template <typename OutputShape>
OutputShape
fe_fdm_second_deriv(const Elem * elem,
                    const Order order,
                    const unsigned int i,
                    const unsigned int j,
                    const Point & p,
                    const bool add_p_level,
                    OutputShape(*deriv_func)
                      (const Elem *, const Order,
                       const unsigned int, const unsigned int,
                       const Point &, const bool));
 
template <typename OutputShape>
OutputShape fe_fdm_second_deriv(const ElemType type,
                                const Order order,
                                const unsigned int i,
                                const unsigned int j,
                                const Point & p,
                                OutputShape(*deriv_func)
                                  (const ElemType, const Order,
                                   const unsigned int,
                                   const unsigned int,
                                   const Point &));
 
template <typename OutputShape>
OutputShape fe_fdm_second_deriv(const ElemType type,
                                const Order order,
                                const Elem * elem,
                                const unsigned int i,
                                const unsigned int j,
                                const Point & p,
                                OutputShape(*deriv_func)
                                  (const ElemType, const Order,
                                   const Elem *,
                                   const unsigned int,
                                   const unsigned int,
                                   const Point &));
 
void lagrange_nodal_soln(const Elem * elem,
                         const Order order,
                         const std::vector<Number> & elem_soln,
                         std::vector<Number> &       nodal_soln,
                         bool add_p_level = true);
 
unsigned int monomial_n_dofs(const ElemType t, const Order o);
 
unsigned int monomial_n_dofs(const Elem * e, const Order o);
 
// shapes[i][j] is shape function phi_i at point p[j]
void rational_fe_weighted_shapes(const Elem * elem,
                                 const FEType underlying_fe_type,
                                 std::vector<std::vector<Real>> & shapes,
                                 const std::vector<Point> & p,
                                 const bool add_p_level);
 
// shapes[i][q] is shape function phi_i at point p[q]
// derivs[j][i][q] is dphi_i/dxi_j at p[q]
void rational_fe_weighted_shapes_derivs(const Elem * elem,
                                        const FEType fe_type,
                                        std::vector<std::vector<Real>> & shapes,
                                        std::vector<std::vector<std::vector<Real>>> & derivs,
                                        const std::vector<Point> & p,
                                        const bool add_p_level);
 
Real rational_fe_shape(const Elem & elem,
                       const FEType underlying_fe_type,
                       const unsigned int i,
                       const Point & p,
                       const bool add_p_level);
 
Real rational_fe_shape_deriv(const Elem & elem,
                             const FEType underlying_fe_type,
                             const unsigned int i,
                             const unsigned int j,
                             const Point & p,
                             const bool add_p_level);
 
Real rational_fe_shape_second_deriv(const Elem & elem,
                                    const FEType underlying_fe_type,
                                    const unsigned int i,
                                    const unsigned int j,
                                    const Point & p,
                                    const bool add_p_level);
 
void rational_all_shapes (const Elem & elem,
                          const FEType underlying_fe_type,
                          const std::vector<Point> & p,
                          std::vector<std::vector<Real>> & v,
                          const bool add_p_level);
 
template <typename OutputShape>
void rational_all_shape_derivs (const Elem & elem,
                                const FEType underlying_fe_type,
                                const std::vector<Point> & p,
                                std::vector<std::vector<OutputShape>> * comps[3],
                                const bool add_p_level);
 
} // namespace libMesh
 
 
// Full specialization of all n_dofs type functions, for every
// dimension, with both original ElemType and new Elem signatures
#define LIBMESH_DEFAULT_NDOFS(MyType) \
template <> unsigned int FE<0,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
template <> unsigned int FE<1,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
template <> unsigned int FE<2,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
template <> unsigned int FE<3,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
 \
template <> unsigned int FE<0,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
template <> unsigned int FE<1,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
template <> unsigned int FE<2,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
template <> unsigned int FE<3,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
 \
template <> unsigned int FE<0,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<1,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<2,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<3,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
 \
template <> unsigned int FE<0,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
template <> unsigned int FE<1,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
template <> unsigned int FE<2,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
template <> unsigned int FE<3,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
 \
template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
 \
template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); }
 
 
#define LIBMESH_DEFAULT_VEC_NDOFS(MyType) \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs(t, o); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs(t, o); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs(t, o); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs(t, o); } \
 \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,MyType>::n_dofs(e, o); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,MyType>::n_dofs(e, o); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,MyType>::n_dofs(e, o); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,MyType>::n_dofs(e, o); } \
 \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(t, o, n); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(t, o, n); } \
 \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(e.type(), o, n); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(e.type(), o, n); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(e.type(), o, n); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(e.type(), o, n); } \
 \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs_per_elem(t, o); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs_per_elem(t, o); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(t, o); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(t, o); } \
 \
template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,MyType>::n_dofs_per_elem(e.type(), o); } \
template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,MyType>::n_dofs_per_elem(e.type(), o); } \
template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(e.type(), o); } \
template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(e.type(), o); }
 
 
#define LIBMESH_DEFAULT_VECTORIZED_FE(MyDim, MyType) \
template<>                                           \
void FE<MyDim,MyType>::all_shapes                    \
  (const Elem * elem,                                \
   const Order o,                                    \
   const std::vector<Point> & p,                     \
   std::vector<std::vector<OutputShape>> & v,        \
   const bool add_p_level)                           \
{                                                    \
  FE<MyDim,MyType>::default_all_shapes               \
    (elem,o,p,v,add_p_level);                        \
}                                                    \
                                                     \
template<>                                           \
void FE<MyDim,MyType>::shapes                        \
  (const Elem * elem,                                \
   const Order o,                                    \
   const unsigned int i,                             \
   const std::vector<Point> & p,                     \
   std::vector<OutputShape> & v,                     \
   const bool add_p_level)                           \
{                                                    \
  FE<MyDim,MyType>::default_shapes                   \
    (elem,o,i,p,v,add_p_level);                      \
}                                                    \
                                                     \
template<>                                           \
void FE<MyDim,MyType>::shape_derivs                  \
  (const Elem * elem,                                \
   const Order o,                                    \
   const unsigned int i,                             \
   const unsigned int j,                             \
   const std::vector<Point> & p,                     \
   std::vector<OutputShape> & v,                     \
   const bool add_p_level)                           \
{                                                    \
  FE<MyDim,MyType>::default_shape_derivs             \
    (elem,o,i,j,p,v,add_p_level);                    \
}                                                    \
                                                     \
template<>                                           \
void FE<MyDim,MyType>::all_shape_derivs              \
  (const Elem * elem,                                \
   const Order o,                                    \
   const std::vector<Point> & p,                     \
   std::vector<std::vector<OutputShape>> * comps[3], \
   const bool add_p_level)                           \
{                                                    \
  FE<MyDim,MyType>::default_all_shape_derivs         \
    (elem,o,p,comps,add_p_level);                    \
}
 
 
#endif // LIBMESH_FE_H