fe.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2018 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 RealGradient type;
};

template<>
struct FEOutputType<NEDELEC_ONE>
{
  typedef RealGradient 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);

  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);

  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);

  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);

  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);

  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);

  static Point map_xi (const Elem * elem,
                       const Point & reference_point);

  static Point map_eta (const Elem * elem,
                        const Point & reference_point);

  static Point map_zeta (const Elem * elem,
                         const Point & 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)
  {}

  static Real hermite_raw_shape_second_deriv(const unsigned int basis_num,
                                             const Real xi);
  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);

template <>
Real FE<2,SUBDIVISION>::shape_deriv(const Elem * elem,
                                    const Order order,
                                    const unsigned int i,
                                    const unsigned int j,
                                    const Point & p);

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);


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);

  static Real regular_shape_second_deriv(const unsigned int i,
                                         const unsigned int j,
                                         const Real v,
                                         const Real w);


  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)
  {}
};




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