fe_abstract.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_ABSTRACT_H
#define LIBMESH_FE_ABSTRACT_H
 
// Local includes
#include "libmesh/reference_counted_object.h"
#include "libmesh/point.h"
#include "libmesh/vector_value.h"
#include "libmesh/fe_type.h"
#include "libmesh/fe_map.h"
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
#include "libmesh/tensor_value.h"
#endif
 
// C++ includes
#include <cstddef>
#include <vector>
#include <memory>
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
#define virtual_for_inffe virtual
#else
#define virtual_for_inffe
#endif
 
 
namespace libMesh
{
 
 
// forward declarations
template <typename T> class DenseMatrix;
template <typename T> class DenseVector;
class BoundaryInfo;
class DofConstraints;
class DofMap;
class Elem;
class MeshBase;
template <typename T> class NumericVector;
class QBase;
enum ElemType : int;
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
class NodeConstraints;
#endif
 
#ifdef LIBMESH_ENABLE_PERIODIC
class PeriodicBoundaries;
class PointLocatorBase;
#endif
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
class InfFE;
#endif
 
 
 
class FEAbstract : public ReferenceCountedObject<FEAbstract>
{
protected:
 
  FEAbstract (const unsigned int dim,
              const FEType & fet);
 
public:
 
  virtual ~FEAbstract();
 
  static std::unique_ptr<FEAbstract> build (const unsigned int dim,
                                            const FEType & type);
 
  virtual void reinit (const Elem * elem,
                       const std::vector<Point> * const pts = nullptr,
                       const std::vector<Real> * const weights = nullptr) = 0;
 
  virtual void reinit_dual_shape_coeffs (const Elem * /*elem*/,
                                         const std::vector<Point> & /*pts*/,
                                         const std::vector<Real> & /*JxW*/)
  {
    libmesh_error_msg("Customized dual shape coefficient calculation has not been implemented for this FE type.");
  }
 
  virtual void reinit_default_dual_shape_coeffs (const Elem * /*elem*/)
  {
    libmesh_error_msg("Default dual shape coefficient calculation has not been implemented for this FE type.");
  }
 
  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) = 0;
 
  virtual void edge_reinit (const Elem * elem,
                            const unsigned int edge,
                            const Real tolerance = TOLERANCE,
                            const std::vector<Point> * pts = nullptr,
                            const std::vector<Real> * weights = nullptr) = 0;
 
  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) = 0;
 
#ifdef LIBMESH_ENABLE_DEPRECATED
  static bool on_reference_element(const Point & p,
                                   const ElemType t,
                                   const Real eps = TOLERANCE);
#endif // LIBMESH_ENABLE_DEPRECATED
 
  static void get_refspace_nodes(const ElemType t,
                                 std::vector<Point> & nodes);
 
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  static void compute_node_constraints (NodeConstraints & constraints,
                                        const Elem * elem);
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
#ifdef LIBMESH_ENABLE_PERIODIC
 
#ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
  static void compute_periodic_node_constraints (NodeConstraints & constraints,
                                                 const PeriodicBoundaries & boundaries,
                                                 const MeshBase & mesh,
                                                 const PointLocatorBase * point_locator,
                                                 const Elem * elem);
#endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
 
#endif // LIBMESH_ENABLE_PERIODIC
 
  unsigned int get_dim() const
  { return dim; }
 
  void get_nothing() const
  { calculate_nothing = true; }
 
  virtual_for_inffe
  const std::vector<Point> & get_xyz() const
  { calculate_map = true; return this->_fe_map->get_xyz(); }
 
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
  virtual const std::vector<Real> & get_JxWxdecay_sq () const
  { return get_JxW();}
#endif
 
  virtual_for_inffe
  const std::vector<Real> & get_JxW() const
  { calculate_map = true; return this->_fe_map->get_JxW(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_dxyzdxi() const
  { calculate_map = true; return this->_fe_map->get_dxyzdxi(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_dxyzdeta() const
  { calculate_map = true; return this->_fe_map->get_dxyzdeta(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_dxyzdzeta() const
  { return _fe_map->get_dxyzdzeta(); }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdxi2() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdxi2(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdeta2() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdeta2(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdzeta2() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdzeta2(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdxideta() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdxideta(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdxidzeta() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdxidzeta(); }
 
  virtual_for_inffe
  const std::vector<RealGradient> & get_d2xyzdetadzeta() const
  { calculate_map = true; return this->_fe_map->get_d2xyzdetadzeta(); }
 
#endif
 
  virtual_for_inffe
  const std::vector<Real> & get_dxidx() const
  { calculate_map = true; return this->_fe_map->get_dxidx(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_dxidy() const
  { calculate_map = true; return this->_fe_map->get_dxidy(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_dxidz() const
  { calculate_map = true; return this->_fe_map->get_dxidz(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_detadx() const
  { calculate_map = true; return this->_fe_map->get_detadx(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_detady() const
  { calculate_map = true; return this->_fe_map->get_detady(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_detadz() const
  { calculate_map = true; return this->_fe_map->get_detadz(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_dzetadx() const
  { calculate_map = true; return this->_fe_map->get_dzetadx(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_dzetady() const
  { calculate_map = true; return this->_fe_map->get_dzetady(); }
 
  virtual_for_inffe
  const std::vector<Real> & get_dzetadz() const
  { calculate_map = true; return this->_fe_map->get_dzetadz(); }
 
  virtual_for_inffe
  const std::vector<std::vector<Point>> & get_tangents() const
  { calculate_map = true; return this->_fe_map->get_tangents(); }
 
  virtual_for_inffe
  const std::vector<Point> & get_normals() const
  { calculate_map = true; return this->_fe_map->get_normals(); }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  virtual_for_inffe
  const std::vector<Real> & get_curvatures() const
  { calculate_map = true; return this->_fe_map->get_curvatures();}
 
#endif
 
  virtual void attach_quadrature_rule (QBase * q) = 0;
 
  virtual unsigned int n_shape_functions () const = 0;
 
  virtual unsigned int n_quadrature_points () const;
 
  const Elem * get_elem()  const { return _elem; }
 
  ElemType get_type() const { return _elem_type; }
 
  unsigned int get_p_level() const { return _p_level; }
 
  FEType get_fe_type()  const { return fe_type; }
 
  Order get_order() const
  { return fe_type.order + _p_level; }
 
  void set_fe_order(int new_order) { fe_type.order = new_order; }
 
  virtual FEContinuity get_continuity() const = 0;
 
  virtual bool is_hierarchic() const = 0;
 
  FEFamily get_family()  const { return fe_type.family; }
 
  const FEMap & get_fe_map() const { return *_fe_map.get(); }
  FEMap & get_fe_map() { return *_fe_map.get(); }
 
  void print_JxW(std::ostream & os) const;
 
  virtual void print_phi(std::ostream & os) const =0;
  virtual void print_dual_phi(std::ostream & os) const =0;
 
  virtual void print_dphi(std::ostream & os) const =0;
  virtual void print_dual_dphi(std::ostream & os) const =0;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  virtual void print_d2phi(std::ostream & os) const =0;
  virtual void print_dual_d2phi(std::ostream & os) const =0;
 
#endif
 
  void print_xyz(std::ostream & os) const;
 
  void print_info(std::ostream & os) const;
 
  friend std::ostream & operator << (std::ostream & os, const FEAbstract & fe);
 
  virtual void request_phi() const = 0;
  virtual void request_dual_phi() const = 0;
 
  virtual void request_dphi() const = 0;
  virtual void request_dual_dphi() const = 0;
 
  void set_calculate_dual(const bool val){calculate_dual = val; }
 
  void set_calculate_default_dual_coeff(const bool val){calculate_default_dual_coeff = val; }
 
  void add_p_level_in_reinit(bool value) { _add_p_level_in_reinit = value; }
 
  bool add_p_level_in_reinit() const { return _add_p_level_in_reinit; }
 
protected:
 
  virtual void compute_shape_functions(const Elem *, const std::vector<Point> & ) =0;
 
  std::unique_ptr<FEMap> _fe_map;
 
 
  const unsigned int dim;
 
  mutable bool calculations_started;
 
  mutable bool calculate_dual;
 
  mutable bool calculate_default_dual_coeff;
 
  mutable bool calculate_nothing;
 
  mutable bool calculate_map;
 
  mutable bool calculate_phi;
 
  mutable bool calculate_dphi;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  mutable bool calculate_d2phi;
#else
  // Otherwise several interfaces need to be redone.
  const bool calculate_d2phi=false;
 
#endif
 
  mutable bool calculate_curl_phi;
 
  mutable bool calculate_div_phi;
 
  mutable bool calculate_dphiref;
 
 
  FEType fe_type;
 
  ElemType _elem_type;
 
  const Elem * _elem;
 
  unsigned int _elem_p_level;
 
  unsigned int _p_level;
 
  QBase * qrule;
 
  bool shapes_on_quadrature;
 
  unsigned int _n_total_qp;
 
  virtual bool shapes_need_reinit() const = 0;
 
  bool _add_p_level_in_reinit;
};
 
} // namespace libMesh
 
#endif // LIBMESH_FE_ABSTRACT_H