fe_base.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_BASE_H
#define LIBMESH_FE_BASE_H
 
// Local includes
#include "libmesh/libmesh_common.h"
#include "libmesh/compare_types.h"
#include "libmesh/fe_abstract.h"
#include "libmesh/fe_transformation_base.h"
#include "libmesh/point.h"
#include "libmesh/reference_counted_object.h"
#include "libmesh/tensor_tools.h"
#include "libmesh/type_n_tensor.h"
#include "libmesh/vector_value.h"
#include "libmesh/dense_matrix.h"
 
// C++ includes
#include <cstddef>
#include <vector>
#include <memory>
 
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;
template <typename T> class FETransformationBase;
class FEType;
 
#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
 
template <typename OutputType>
class FEGenericBase : public FEAbstract
{
protected:
 
  FEGenericBase (const unsigned int dim,
                 const FEType & fet);
 
public:
 
  virtual ~FEGenericBase();
 
  static std::unique_ptr<FEGenericBase> build (const unsigned int dim,
                                               const FEType & type);
 
  typedef OutputType                                                      OutputShape;
  typedef typename TensorTools::IncrementRank<OutputShape>::type          OutputGradient;
  typedef typename TensorTools::IncrementRank<OutputGradient>::type       OutputTensor;
  typedef typename TensorTools::DecrementRank<OutputShape>::type          OutputDivergence;
  typedef typename TensorTools::MakeNumber<OutputShape>::type             OutputNumber;
  typedef typename TensorTools::IncrementRank<OutputNumber>::type         OutputNumberGradient;
  typedef typename TensorTools::IncrementRank<OutputNumberGradient>::type OutputNumberTensor;
  typedef typename TensorTools::DecrementRank<OutputNumber>::type         OutputNumberDivergence;
 
 
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  static std::unique_ptr<FEGenericBase> build_InfFE (const unsigned int dim,
                                                     const FEType & type);
 
#endif
 
#ifdef LIBMESH_ENABLE_AMR
 
  static void compute_proj_constraints (DofConstraints & constraints,
                                        DofMap & dof_map,
                                        const unsigned int variable_number,
                                        const Elem * elem);
 
  static void coarsened_dof_values(const NumericVector<Number> & global_vector,
                                   const DofMap & dof_map,
                                   const Elem * coarse_elem,
                                   DenseVector<Number> & coarse_dofs,
                                   const unsigned int var,
                                   const bool use_old_dof_indices = false);
 
  static void coarsened_dof_values(const NumericVector<Number> & global_vector,
                                   const DofMap & dof_map,
                                   const Elem * coarse_elem,
                                   DenseVector<Number> & coarse_dofs,
                                   const bool use_old_dof_indices = false);
 
#endif // #ifdef LIBMESH_ENABLE_AMR
 
#ifdef LIBMESH_ENABLE_PERIODIC
 
  static void compute_periodic_constraints (DofConstraints & constraints,
                                            DofMap & dof_map,
                                            const PeriodicBoundaries & boundaries,
                                            const MeshBase & mesh,
                                            const PointLocatorBase * point_locator,
                                            const unsigned int variable_number,
                                            const Elem * elem);
 
#endif // LIBMESH_ENABLE_PERIODIC
 
  const std::vector<std::vector<OutputShape>> & get_phi() const
  { libmesh_assert(!calculations_started || calculate_phi);
    calculate_phi = true; return phi; }
 
  const std::vector<std::vector<OutputShape>> & get_dual_phi() const
  {
    libmesh_assert(!calculations_started || calculate_dual);
    calculate_dual = true;
    // Dual phi computation relies on primal phi computation
    this->request_phi();
    return dual_phi;
  }
 
  virtual void request_phi() const override
  { get_phi(); }
 
  virtual void request_dual_phi() const override
  { get_dual_phi(); }
 
  const std::vector<std::vector<OutputGradient>> & get_dphi() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = calculate_dphiref = true; return dphi; }
 
  const std::vector<std::vector<OutputGradient>> & get_dual_dphi() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = calculate_dual = calculate_dphiref = true; return dual_dphi; }
 
  virtual void request_dphi() const override
  { get_dphi(); }
 
  virtual void request_dual_dphi() const override
  { get_dual_dphi(); }
 
  const DenseMatrix<Real> & get_dual_coeff() const
  { return dual_coeff; }
 
  virtual_for_inffe
  const std::vector<std::vector<OutputShape>> & get_curl_phi() const
  { libmesh_assert(!calculations_started || calculate_curl_phi);
    calculate_curl_phi = calculate_dphiref = true; return curl_phi; }
 
  virtual_for_inffe
  const std::vector<std::vector<OutputDivergence>> & get_div_phi() const
  { libmesh_assert(!calculations_started || calculate_div_phi);
    calculate_div_phi = calculate_dphiref = true; return div_phi; }
 
  const std::vector<std::vector<OutputShape>> & get_dphidx() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = calculate_dphiref = true; return dphidx; }
 
  const std::vector<std::vector<OutputShape>> & get_dphidy() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = calculate_dphiref = true; return dphidy; }
 
  const std::vector<std::vector<OutputShape>> & get_dphidz() const
  { libmesh_assert(!calculations_started || calculate_dphi);
    calculate_dphi = calculate_dphiref = true; return dphidz; }
 
  const std::vector<std::vector<OutputShape>> & get_dphidxi() const
  { libmesh_assert(!calculations_started || calculate_dphiref);
    calculate_dphiref = true; return dphidxi; }
 
  const std::vector<std::vector<OutputShape>> & get_dphideta() const
  { libmesh_assert(!calculations_started || calculate_dphiref);
    calculate_dphiref = true; return dphideta; }
 
  const std::vector<std::vector<OutputShape>> & get_dphidzeta() const
  { libmesh_assert(!calculations_started || calculate_dphiref);
    calculate_dphiref = true; return dphidzeta; }
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  const std::vector<std::vector<OutputTensor>> & get_d2phi() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phi; }
 
  const std::vector<std::vector<OutputTensor>> & get_dual_d2phi() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dual = calculate_dphiref = true; return dual_d2phi; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidx2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidx2; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidxdy() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidxdy; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidxdz() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidxdz; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidy2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi =  calculate_dphiref = true; return d2phidy2; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidydz() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidydz; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidz2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidz2; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidxi2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidxi2; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidxideta() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidxideta; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidxidzeta() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidxidzeta; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phideta2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phideta2; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidetadzeta() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidetadzeta; }
 
  const std::vector<std::vector<OutputShape>> & get_d2phidzeta2() const
  { libmesh_assert(!calculations_started || calculate_d2phi);
    calculate_d2phi = calculate_dphiref = true; return d2phidzeta2; }
 
#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  const std::vector<OutputGradient> & get_dphase() const
  { return dphase; }
 
 
  virtual const std::vector<Real> & get_Sobolev_weight() const
  { return weight; }
 
  virtual const std::vector<RealGradient> & get_Sobolev_dweight() const
  { return dweight; }
 
  virtual const std::vector<Real> & get_Sobolev_weightxR_sq() const
  { return weight; }
 
  virtual const std::vector<RealGradient> & get_Sobolev_dweightxR_sq() const
  { return dweight; }
 
  virtual const std::vector<std::vector<OutputShape>> & get_phi_over_decayxR () const
  { return get_phi();}
 
  virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decayxR () const
  { return get_dphi();}
 
  virtual const std::vector<std::vector<OutputGradient>> & get_dphi_over_decay () const
  { return get_dphi();}
 
#endif
 
  virtual void print_phi(std::ostream & os) const override;
  virtual void print_dual_phi(std::ostream & os) const override;
 
  virtual void print_dphi(std::ostream & os) const override;
  virtual void print_dual_dphi(std::ostream & os) const override;
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  virtual void print_d2phi(std::ostream & os) const override;
  virtual void print_dual_d2phi(std::ostream & os) const override;
 
#endif
 
 
protected:
 
 
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  virtual void init_base_shape_functions(const std::vector<Point> & qp,
                                         const Elem * e) = 0;
 
#endif
 
  virtual_for_inffe
  void determine_calculations();
 
  bool calculating_nothing() const
  {
    return calculate_nothing &&
      !this->calculate_phi && !this->calculate_dphi &&
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
      !this->calculate_d2phi &&
#endif
      !this->calculate_curl_phi && !this->calculate_div_phi &&
      !this->calculate_map;
  }
 
  virtual void compute_shape_functions(const Elem * elem, const std::vector<Point> & qp) override;
 
  void compute_dual_shape_coeffs(const std::vector<Real> & JxW, const std::vector<std::vector<OutputShape>> & phi);
 
  void compute_dual_shape_functions();
 
  std::unique_ptr<FETransformationBase<OutputType>> _fe_trans;
 
  std::vector<std::vector<OutputShape>>   phi;
  std::vector<std::vector<OutputShape>>   dual_phi;
 
  std::vector<std::vector<OutputGradient>>  dphi;
  std::vector<std::vector<OutputGradient>>  dual_dphi;
 
  mutable DenseMatrix<Real> dual_coeff;
 
  std::vector<std::vector<OutputShape>> curl_phi;
 
  std::vector<std::vector<OutputDivergence>> div_phi;
 
  std::vector<std::vector<OutputShape>>   dphidxi;
 
  std::vector<std::vector<OutputShape>>   dphideta;
 
  std::vector<std::vector<OutputShape>>   dphidzeta;
 
  std::vector<std::vector<OutputShape>>   dphidx;
 
  std::vector<std::vector<OutputShape>>   dphidy;
 
  std::vector<std::vector<OutputShape>>   dphidz;
 
 
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
 
  std::vector<std::vector<OutputTensor>>  d2phi;
  std::vector<std::vector<OutputTensor>>  dual_d2phi;
 
  std::vector<std::vector<OutputShape>>   d2phidxi2;
 
  std::vector<std::vector<OutputShape>>   d2phidxideta;
 
  std::vector<std::vector<OutputShape>>   d2phidxidzeta;
 
  std::vector<std::vector<OutputShape>>   d2phideta2;
 
  std::vector<std::vector<OutputShape>>   d2phidetadzeta;
 
  std::vector<std::vector<OutputShape>>   d2phidzeta2;
 
  std::vector<std::vector<OutputShape>>   d2phidx2;
 
  std::vector<std::vector<OutputShape>>   d2phidxdy;
 
  std::vector<std::vector<OutputShape>>   d2phidxdz;
 
  std::vector<std::vector<OutputShape>>   d2phidy2;
 
  std::vector<std::vector<OutputShape>>   d2phidydz;
 
  std::vector<std::vector<OutputShape>>   d2phidz2;
 
#endif
 
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  //--------------------------------------------------------------
  /* protected members for infinite elements, which are accessed
   * from the outside through some inline functions
   */
 
 
  std::vector<OutputGradient> dphase;
 
  std::vector<RealGradient> dweight;
 
  std::vector<Real>  weight;
 
#endif
 
private:
 
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
  template <unsigned int friend_Dim, FEFamily friend_T_radial, InfMapType friend_T_map>
  friend class InfFE;
 
#endif
 
 
};
 
// --------------------------------------------------------------------
// Generic templates. We specialize for OutputType = Real, so these are
// only used for OutputType = RealVectorValue
template <typename OutputType>
void FEGenericBase<OutputType>::compute_dual_shape_functions ()
{
  libmesh_error_msg(
      "Computation of dual shape functions for vector finite element "
      "families is not currently implemented");
}
 
template <typename OutputType>
void FEGenericBase<OutputType>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/)
{
  libmesh_error_msg(
      "Computation of dual shape functions for vector finite element "
      "families is not currently implemented");
}
 
// -----------------------------------------------------------
// Forward declaration of specialization
template <>
void FEGenericBase<Real>::compute_dual_shape_functions();
 
template <>
void FEGenericBase<Real>::compute_dual_shape_coeffs(const std::vector<Real> & /*JxW*/, const std::vector<std::vector<OutputShape>> & /*phi_vals*/);
 
 
// Typedefs for convenience and backwards compatibility
typedef FEGenericBase<Real> FEBase;
typedef FEGenericBase<RealGradient> FEVectorBase;
 
 
 
 
// ------------------------------------------------------------
// FEGenericBase class inline members
template <typename OutputType>
inline
FEGenericBase<OutputType>::FEGenericBase(const unsigned int d,
                                         const FEType & fet) :
  FEAbstract(d,fet),
  _fe_trans( FETransformationBase<OutputType>::build(fet) ),
  phi(),
  dual_phi(),
  dphi(),
  dual_dphi(),
  curl_phi(),
  div_phi(),
  dphidxi(),
  dphideta(),
  dphidzeta(),
  dphidx(),
  dphidy(),
  dphidz()
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
  ,d2phi(),
  dual_d2phi(),
  d2phidxi2(),
  d2phidxideta(),
  d2phidxidzeta(),
  d2phideta2(),
  d2phidetadzeta(),
  d2phidzeta2(),
  d2phidx2(),
  d2phidxdy(),
  d2phidxdz(),
  d2phidy2(),
  d2phidydz(),
  d2phidz2()
#endif
#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
  ,dphase(),
  dweight(),
  weight()
#endif
{
}
 
 
 
template <typename OutputType>
inline
FEGenericBase<OutputType>::~FEGenericBase()
{
}
 
} // namespace libMesh
 
#endif // LIBMESH_FE_BASE_H