fe.h

Go to the documentation of this file.

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


//-------------------------------------------------------------
// FE class definition
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;

  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;

  virtual bool is_hierarchic() const;

  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 = NULL,
                       const std::vector<Real>* const weights = NULL);

  virtual void reinit (const Elem* elem,
                       const unsigned int side,
                       const Real tolerance = TOLERANCE,
                       const std::vector<Point>* const pts = NULL,
                       const std::vector<Real>* const weights = NULL);

  virtual void edge_reinit (const Elem* elem,
                            const unsigned int edge,
                            const Real tolerance = TOLERANCE,
                            const std::vector<Point>* const pts = NULL,
                            const std::vector<Real>* const weights = NULL);

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

  virtual void attach_quadrature_rule (QBase* q);

  virtual unsigned int n_quadrature_points () const;

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

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

#endif

  std::vector<Point> cached_nodes;

  ElemType last_side;

  unsigned int last_edge;
};



//-------------------------------------------------------------
// FEHierarchic class definition
template <unsigned int Dim>
class FEClough : public FE<Dim,CLOUGH>
{
public:

  explicit
  FEClough(const FEType& fet);
};



//-------------------------------------------------------------
// FEHierarchic class definition
template <unsigned int Dim>
class FEHermite : public FE<Dim,HERMITE>
{
public:

  explicit
  FEHermite(const FEType& 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);
};

//-------------------------------------------------------------
// FESubdivision class definition


// template specialization prototypes, needed for being able to
// call them from inside FESubdivision::init_shape_functions

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 = NULL,
                       const std::vector<Real>* const weights = NULL);

  virtual void reinit (const Elem*,
                       const unsigned int,
                       const Real = TOLERANCE,
                       const std::vector<Point>* const = NULL,
                       const std::vector<Real>* const = NULL)
  { libmesh_not_implemented(); }

  virtual void attach_quadrature_rule (QBase* q);

  virtual void init_shape_functions(const std::vector<Point>& qp,
                                    const Elem* elem);

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



//-------------------------------------------------------------
// FEHierarchic class definition
template <unsigned int Dim>
class FEHierarchic : public FE<Dim,HIERARCHIC>
{
public:

  explicit
  FEHierarchic(const FEType& fet);
};



//-------------------------------------------------------------
// FEL2Hierarchic class definition
template <unsigned int Dim>
class FEL2Hierarchic : public FE<Dim,L2_HIERARCHIC>
{
public:

  explicit
  FEL2Hierarchic(const FEType& fet);
};



//-------------------------------------------------------------
// FELagrange class definition
template <unsigned int Dim>
class FELagrange : public FE<Dim,LAGRANGE>
{
public:

  explicit
  FELagrange(const FEType& fet);
};


//-------------------------------------------------------------
// FEL2Lagrange class definition
template <unsigned int Dim>
class FEL2Lagrange : public FE<Dim,L2_LAGRANGE>
{
public:

  explicit
  FEL2Lagrange(const FEType& fet);
};


//-------------------------------------------------------------
// FEMonomial class definition
template <unsigned int Dim>
class FEMonomial : public FE<Dim,MONOMIAL>
{
public:

  explicit
  FEMonomial(const FEType& fet);
};


//-------------------------------------------------------------
// FEScalar class definition
template <unsigned int Dim>
class FEScalar : public FE<Dim,SCALAR>
{
public:

  explicit
  FEScalar(const FEType& fet);
};


//-------------------------------------------------------------
// FEXYZ class definition
template <unsigned int Dim>
class FEXYZ : public FE<Dim,XYZ>
{
public:

  explicit
  FEXYZ(const FEType& fet);

  virtual void reinit (const Elem* elem,
                       const std::vector<Point>* const pts = NULL,
                       const std::vector<Real>* const weights = NULL)
  { 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 = NULL,
                       const std::vector<Real>* const weights = NULL);


protected:

  virtual void init_shape_functions(const std::vector<Point>& qp,
                                    const Elem* e);

  virtual void compute_shape_functions(const Elem* elem, const std::vector<Point>& qp);

  void compute_face_values (const Elem* elem,
                            const Elem* side,
                            const std::vector<Real>& weights);
};

//-------------------------------------------------------------
// FELagrangeVec class definition
template <unsigned int Dim>
class FELagrangeVec : public FE<Dim,LAGRANGE_VEC>
{
public:

  explicit
  FELagrangeVec(const FEType& fet);

};

//-------------------------------------------------------------
// FENedelecOne class definition
template <unsigned int Dim>
class FENedelecOne : public FE<Dim,NEDELEC_ONE>
{
public:

  explicit
  FENedelecOne(const FEType& 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;

}




// ------------------------------------------------------------
// FE class inline members
template <unsigned int Dim, FEFamily T>
inline
FE<Dim,T>::FE (const FEType& fet) :
  FEGenericBase<typename FEOutputType<T>::type> (Dim,fet),
  last_side(INVALID_ELEM),
  last_edge(libMesh::invalid_uint)
{
  // Sanity check.  Make sure the
  // Family specified in the template instantiation
  // matches the one in the FEType object
  libmesh_assert_equal_to (T, this->get_family());
}



// ------------------------------------------------------------
// FEClough class inline members
template <unsigned int Dim>
inline
FEClough<Dim>::FEClough (const FEType& fet) :
  FE<Dim,CLOUGH> (fet)
{
}



// ------------------------------------------------------------
// FEHermite class inline members
template <unsigned int Dim>
inline
FEHermite<Dim>::FEHermite (const FEType& fet) :
  FE<Dim,HERMITE> (fet)
{
}



// ------------------------------------------------------------
// FEHierarchic class inline members
template <unsigned int Dim>
inline
FEHierarchic<Dim>::FEHierarchic (const FEType& fet) :
  FE<Dim,HIERARCHIC> (fet)
{
}



// ------------------------------------------------------------
// FEL2Hierarchic class inline members
template <unsigned int Dim>
inline
FEL2Hierarchic<Dim>::FEL2Hierarchic (const FEType& fet) :
  FE<Dim,L2_HIERARCHIC> (fet)
{
}



// ------------------------------------------------------------
// FELagrange class inline members
template <unsigned int Dim>
inline
FELagrange<Dim>::FELagrange (const FEType& fet) :
  FE<Dim,LAGRANGE> (fet)
{
}

// ------------------------------------------------------------
// FELagrangeVec class inline members
template <unsigned int Dim>
inline
FELagrangeVec<Dim>::FELagrangeVec (const FEType& fet) :
  FE<Dim,LAGRANGE_VEC> (fet)
{
}

// ------------------------------------------------------------
// FEL2Lagrange class inline members
template <unsigned int Dim>
inline
FEL2Lagrange<Dim>::FEL2Lagrange (const FEType& fet) :
  FE<Dim,L2_LAGRANGE> (fet)
{
}



// ------------------------------------------------------------
// FEMonomial class inline members
template <unsigned int Dim>
inline
FEMonomial<Dim>::FEMonomial (const FEType& fet) :
  FE<Dim,MONOMIAL> (fet)
{
}




// ------------------------------------------------------------
// FEXYZ class inline members
template <unsigned int Dim>
inline
FEXYZ<Dim>::FEXYZ (const FEType& fet) :
  FE<Dim,XYZ> (fet)
{
}

// ------------------------------------------------------------
// FEScalar class inline members
template <unsigned int Dim>
inline
FEScalar<Dim>::FEScalar (const FEType& fet) :
  FE<Dim,SCALAR> (fet)
{
}

// ------------------------------------------------------------
// FENedelecOne class inline members
template <unsigned int Dim>
inline
FENedelecOne<Dim>::FENedelecOne (const FEType& fet) :
  FE<Dim,NEDELEC_ONE> (fet)
{
}

} // namespace libMesh

#endif // LIBMESH_FE_H