fe_abstract.h

Go to the documentation of this file.

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