fe_base.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_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