fem_context.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_FEM_CONTEXT_H
#define LIBMESH_FEM_CONTEXT_H

// Local Includes
#include "libmesh/diff_context.h"
#include "libmesh/id_types.h"
#include "libmesh/fe_type.h"
#include "libmesh/fe_base.h"
#include "libmesh/vector_value.h"

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
#include "libmesh/tensor_value.h"
#endif

// C++ includes
#include <map>
#include <set>

namespace libMesh
{

// Forward Declarations
class BoundaryInfo;
class Elem;
template <typename T> class FEGenericBase;
typedef FEGenericBase<Real> FEBase;
class QBase;
class Point;
template <typename T> class NumericVector;

class FEMContext : public DiffContext
{
public:

  explicit
  FEMContext (const System & sys,
              const std::vector<unsigned int> * active_vars = nullptr,
              bool allocate_local_matrices = true);

  explicit
  FEMContext (const System & sys,
              int extra_quadrature_order,
              const std::vector<unsigned int> * active_vars = nullptr,
              bool allocate_local_matrices = true);


  virtual ~FEMContext ();

  void use_default_quadrature_rules(int extra_quadrature_order=0);

  void use_unweighted_quadrature_rules(int extra_quadrature_order=0);

  bool has_side_boundary_id(boundary_id_type id) const;

  void side_boundary_ids(std::vector<boundary_id_type> & vec_to_fill) const;

  Number interior_value(unsigned int var, unsigned int qp) const;

  Number side_value(unsigned int var, unsigned int qp) const;

  Number point_value(unsigned int var, const Point & p) const;

  Gradient interior_gradient(unsigned int var, unsigned int qp) const;

  Gradient side_gradient(unsigned int var, unsigned int qp) const;

  Gradient point_gradient(unsigned int var, const Point & p) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  Tensor interior_hessian(unsigned int var, unsigned int qp) const;

  Tensor side_hessian(unsigned int var, unsigned int qp) const;

  Tensor point_hessian(unsigned int var, const Point & p) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  Number fixed_interior_value(unsigned int var, unsigned int qp) const;

  Number fixed_side_value(unsigned int var, unsigned int qp) const;

  Number fixed_point_value(unsigned int var, const Point & p) const;

  Gradient fixed_interior_gradient(unsigned int var, unsigned int qp) const;

  Gradient fixed_side_gradient(unsigned int var, unsigned int qp) const;

  Gradient fixed_point_gradient(unsigned int var, const Point & p) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  Tensor fixed_interior_hessian(unsigned int var, unsigned int qp) const;

  Tensor fixed_side_hessian(unsigned int var, unsigned int qp) const;

  Tensor fixed_point_hessian (unsigned int var, const Point & p) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputShape>
  void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
  { this->get_element_fe<OutputShape>(var,fe,this->get_elem_dim()); }

  FEBase * get_element_fe( unsigned int var ) const
  { return this->get_element_fe(var,this->get_elem_dim()); }

  template<typename OutputShape>
  void get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                       unsigned short dim ) const;

  void get_element_fe( unsigned int var, FEAbstract *& fe,
                       unsigned short dim ) const;

  FEBase * get_element_fe( unsigned int var, unsigned short dim ) const;

  template<typename OutputShape>
  void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
  { this->get_side_fe<OutputShape>(var,fe,this->get_elem_dim()); }

  FEBase * get_side_fe( unsigned int var ) const
  { return this->get_side_fe(var,this->get_elem_dim()); }

  template<typename OutputShape>
  void get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                    unsigned short dim ) const;

  void get_side_fe( unsigned int var, FEAbstract *& fe,
                    unsigned short dim ) const;

  FEBase * get_side_fe( unsigned int var, unsigned short dim ) const;

  template<typename OutputShape>
  void get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const;

  void get_edge_fe( unsigned int var, FEAbstract *& fe ) const;

  FEBase * get_edge_fe( unsigned int var ) const;

  template<typename OutputType>
  void interior_value(unsigned int var,
                      unsigned int qp,
                      OutputType & u) const;

  template<typename OutputType>
  void interior_values(unsigned int var,
                       const NumericVector<Number> & _system_vector,
                       std::vector<OutputType> & interior_values_vector) const;

  template<typename OutputType>
  void side_value(unsigned int var,
                  unsigned int qp,
                  OutputType & u) const;

  template<typename OutputType>
  void side_values(unsigned int var,
                   const NumericVector<Number> & _system_vector,
                   std::vector<OutputType> & side_values_vector) const;

  template<typename OutputType>
  void point_value(unsigned int var,
                   const Point & p,
                   OutputType & u,
                   const Real tolerance = TOLERANCE) const;

  template<typename OutputType>
  void interior_gradient(unsigned int var, unsigned int qp,
                         OutputType & du) const;

  template<typename OutputType>
  void interior_gradients(unsigned int var,
                          const NumericVector<Number> & _system_vector,
                          std::vector<OutputType> & interior_gradients_vector) const;

  template<typename OutputType>
  void side_gradient(unsigned int var,
                     unsigned int qp,
                     OutputType & du) const;

  template<typename OutputType>
  void side_gradients(unsigned int var,
                      const NumericVector<Number> & _system_vector,
                      std::vector<OutputType> & side_gradients_vector) const;

  template<typename OutputType>
  void point_gradient(unsigned int var,
                      const Point & p,
                      OutputType & grad_u,
                      const Real tolerance = TOLERANCE) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_hessian(unsigned int var,
                        unsigned int qp,
                        OutputType & d2u) const;

  template<typename OutputType>
  void interior_hessians(unsigned int var,
                         const NumericVector<Number> & _system_vector,
                         std::vector<OutputType> & d2u_vals) const;

  template<typename OutputType>
  void side_hessian(unsigned int var,
                    unsigned int qp,
                    OutputType & d2u) const;

  template<typename OutputType>
  void side_hessians(unsigned int var,
                     const NumericVector<Number> & _system_vector,
                     std::vector<OutputType> & d2u_vals) const;

  template<typename OutputType>
  void point_hessian(unsigned int var,
                     const Point & p,
                     OutputType & hess_u,
                     const Real tolerance = TOLERANCE) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_rate(unsigned int var,
                     unsigned int qp,
                     OutputType & u) const;


  template<typename OutputType>
  void interior_rate_gradient(unsigned int var,
                              unsigned int qp,
                              OutputType & u) const;


  template<typename OutputType>
  void side_rate(unsigned int var,
                 unsigned int qp,
                 OutputType & u) const;

  template<typename OutputType>
  void point_rate(unsigned int var,
                  const Point & p,
                  OutputType & u) const;

  template<typename OutputType>
  void interior_accel(unsigned int var,
                      unsigned int qp,
                      OutputType & u) const;

  template<typename OutputType>
  void side_accel(unsigned int var,
                  unsigned int qp,
                  OutputType & u) const;

  template<typename OutputType>
  void point_accel(unsigned int var,
                   const Point & p,
                   OutputType & u) const;

  template<typename OutputType>
  void fixed_interior_value(unsigned int var,
                            unsigned int qp,
                            OutputType & u) const;

  template<typename OutputType>
  void fixed_side_value(unsigned int var,
                        unsigned int qp,
                        OutputType & u) const;

  template<typename OutputType>
  void fixed_point_value(unsigned int var,
                         const Point & p,
                         OutputType & u,
                         const Real tolerance = TOLERANCE) const;

  template<typename OutputType>
  void fixed_interior_gradient(unsigned int var,
                               unsigned int qp,
                               OutputType & grad_u) const;

  template<typename OutputType>
  void fixed_side_gradient(unsigned int var,
                           unsigned int qp,
                           OutputType & grad_u) const;

  template<typename OutputType>
  void fixed_point_gradient(unsigned int var,
                            const Point & p,
                            OutputType & grad_u,
                            const Real tolerance = TOLERANCE) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void fixed_interior_hessian(unsigned int var,
                              unsigned int qp,
                              OutputType & hess_u) const;

  template<typename OutputType>
  void fixed_side_hessian(unsigned int var,
                          unsigned int qp,
                          OutputType & hess_u) const;

  template<typename OutputType>
  void fixed_point_hessian(unsigned int var,
                           const Point & p,
                           OutputType & hess_u,
                           const Real tolerance = TOLERANCE) const;

#endif // LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType>
  void interior_curl(unsigned int var,
                     unsigned int qp,
                     OutputType & curl_u) const;

  template<typename OutputType>
  void point_curl(unsigned int var,
                  const Point & p,
                  OutputType & curl_u,
                  const Real tolerance = TOLERANCE) const;

  template<typename OutputType>
  void interior_div(unsigned int var,
                    unsigned int qp,
                    OutputType & div_u) const;

  // should be protected:
  virtual void elem_reinit(Real theta) override;

  virtual void elem_side_reinit(Real theta) override;

  virtual void elem_edge_reinit(Real theta) override;

  virtual void nonlocal_reinit(Real theta) override;

  virtual void pre_fe_reinit(const System &,
                             const Elem * e);

  virtual void elem_fe_reinit(const std::vector<Point> * const pts = nullptr);

  virtual void side_fe_reinit();

  virtual void edge_fe_reinit();

  const QBase & get_element_qrule() const
  { return this->get_element_qrule(this->get_elem_dim()); }

  const QBase & get_side_qrule() const
  { return this->get_side_qrule(this->get_elem_dim()); }

  const QBase & get_element_qrule( unsigned short dim ) const
  { libmesh_assert(_element_qrule[dim]);
    return *(this->_element_qrule[dim]); }

  const QBase & get_side_qrule( unsigned short dim ) const
  {
    libmesh_assert(_side_qrule[dim]);
    return *(this->_side_qrule[dim]);
  }

  const QBase & get_edge_qrule() const
  { return *(this->_edge_qrule); }

  virtual void set_mesh_system(System * sys)
  { this->_mesh_sys = sys; }

  const System * get_mesh_system() const
  { return this->_mesh_sys; }

  System * get_mesh_system()
  { return this->_mesh_sys; }

  unsigned int get_mesh_x_var() const
  { return _mesh_x_var; }

  void set_mesh_x_var(unsigned int x_var)
  { _mesh_x_var = x_var; }

  unsigned int get_mesh_y_var() const
  { return _mesh_y_var; }

  void set_mesh_y_var(unsigned int y_var)
  { _mesh_y_var = y_var; }

  unsigned int get_mesh_z_var() const
  { return _mesh_z_var; }

  void set_mesh_z_var(unsigned int z_var)
  { _mesh_z_var = z_var; }

  bool has_elem() const
  { return (this->_elem != nullptr); }

  const Elem & get_elem() const
  { libmesh_assert(this->_elem);
    return *(this->_elem); }

  Elem & get_elem()
  { libmesh_assert(this->_elem);
    return *(const_cast<Elem *>(this->_elem)); }

  unsigned char get_side() const
  { return side; }

  unsigned char get_edge() const
  { return edge; }

  unsigned char get_dim() const
  { return this->_dim; }

  unsigned char get_elem_dim() const
  { return this->_elem ? this->_elem_dim : this->_dim; }

  const std::set<unsigned char> & elem_dimensions() const
  { return _elem_dims; }

  void elem_position_set(Real theta);

  void elem_position_get();

  enum AlgebraicType { NONE = 0,  // Do not reinitialize dof_indices
                       DOFS_ONLY, // Reinitialize dof_indices, not
                                  // algebraic structures
                       CURRENT,   // Use dof_indices, current solution
                       OLD,       // Use old_dof_indices, custom solution
                       OLD_DOFS_ONLY}; // Reinitialize old_dof_indices, not
                                       // algebraic structures

  void set_algebraic_type(const AlgebraicType atype)
  { _atype = atype; }

  /*
   * Get the current AlgebraicType setting
   */
  AlgebraicType algebraic_type() const { return _atype; }

  void set_custom_solution(const NumericVector<Number> * custom_sol)
  { _custom_solution = custom_sol; }

  void set_jacobian_tolerance(Real tol);

  unsigned char side;

  unsigned char edge;

  FEType find_hardest_fe_type();

  void attach_quadrature_rules();

  const std::vector<unsigned int> * active_vars() const { return _active_vars.get(); }

  template<typename OutputShape>
  FEGenericBase<OutputShape> * build_new_fe(const FEGenericBase<OutputShape> * fe,
                                            const Point & p,
                                            const Real tolerance = TOLERANCE,
                                            const int get_derivative_level = -1) const;

protected:

  std::unique_ptr<const std::vector<unsigned int>> _active_vars;

  System * _mesh_sys;

  unsigned int _mesh_x_var, _mesh_y_var, _mesh_z_var;

  AlgebraicType _atype;

  const NumericVector<Number> * _custom_solution;

  mutable std::unique_ptr<FEGenericBase<Real>>         _real_fe;
  mutable std::unique_ptr<FEGenericBase<RealGradient>> _real_grad_fe;
  mutable int _real_fe_derivative_level;
  mutable int _real_grad_fe_derivative_level;

#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
  mutable bool _real_fe_is_inf;
  mutable bool _real_grad_fe_is_inf;
#endif

  template<typename OutputShape>
  FEGenericBase<OutputShape> * cached_fe( const unsigned int elem_dim,
                                          const FEType fe_type,
                                          const int get_derivative_level ) const;

  void set_elem( const Elem * e );

  // gcc-3.4, oracle 12.3 require this typedef to be public
  // in order to use it in a return type
public:

  typedef const DenseSubVector<Number> & (DiffContext::*diff_subsolution_getter)(unsigned int) const;

protected:
  template <typename OutputType>
  struct FENeeded
  {
    // Rank decrementer helper types
    typedef typename TensorTools::DecrementRank<OutputType>::type Rank1Decrement;
    typedef typename TensorTools::DecrementRank<Rank1Decrement>::type Rank2Decrement;

    // Typedefs for "Value getter" function pointer
    typedef typename TensorTools::MakeReal<OutputType>::type value_shape;
    typedef FEGenericBase<value_shape> value_base;
    typedef void (FEMContext::*value_getter) (unsigned int, value_base *&, unsigned short) const;

    // Typedefs for "Grad getter" function pointer
    typedef typename TensorTools::MakeReal<Rank1Decrement>::type grad_shape;
    typedef FEGenericBase<grad_shape> grad_base;
    typedef void (FEMContext::*grad_getter) (unsigned int, grad_base *&, unsigned short) const;

    // Typedefs for "Hessian getter" function pointer
    typedef typename TensorTools::MakeReal<Rank2Decrement>::type hess_shape;
    typedef FEGenericBase<hess_shape> hess_base;
    typedef void (FEMContext::*hess_getter) (unsigned int, hess_base *&, unsigned short) const;
  };



  template<typename OutputType,
           typename FENeeded<OutputType>::value_getter fe_getter,
           diff_subsolution_getter subsolution_getter>
  void some_value(unsigned int var, unsigned int qp, OutputType & u) const;

  template<typename OutputType,
           typename FENeeded<OutputType>::grad_getter fe_getter,
           diff_subsolution_getter subsolution_getter>
  void some_gradient(unsigned int var, unsigned int qp, OutputType & u) const;

#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES

  template<typename OutputType,
           typename FENeeded<OutputType>::hess_getter fe_getter,
           diff_subsolution_getter subsolution_getter>
  void some_hessian(unsigned int var, unsigned int qp, OutputType & u) const;
#endif

  std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>> _element_fe;
  std::vector<std::map<FEType, std::unique_ptr<FEAbstract>>> _side_fe;
  std::map<FEType, std::unique_ptr<FEAbstract>> _edge_fe;


  std::vector<std::vector<FEAbstract *>> _element_fe_var;
  std::vector<std::vector<FEAbstract *>> _side_fe_var;
  std::vector<FEAbstract *> _edge_fe_var;

  const BoundaryInfo & _boundary_info;

  const Elem * _elem;

  unsigned char _dim;

  unsigned char _elem_dim;

  std::set<unsigned char> _elem_dims;

  std::vector<std::unique_ptr<QBase>> _element_qrule;

  std::vector<std::unique_ptr<QBase>> _side_qrule;

  std::unique_ptr<QBase> _edge_qrule;

  int _extra_quadrature_order;

private:
  void init_internal_data(const System & sys);

  void _do_elem_position_set(Real theta);

  void _update_time_from_system(Real theta);
};



// ------------------------------------------------------------
// FEMContext inline methods

inline
void FEMContext::elem_position_set(Real theta)
{
  if (_mesh_sys)
    this->_do_elem_position_set(theta);
}

template<typename OutputShape>
inline
void FEMContext::get_element_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                                 unsigned short dim ) const
{
  libmesh_assert( !_element_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
  fe = cast_ptr<FEGenericBase<OutputShape> *>( (_element_fe_var[dim][var] ) );
}

inline
void FEMContext::get_element_fe( unsigned int var, FEAbstract *& fe,
                                 unsigned short dim ) const
{
  libmesh_assert( !_element_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
  fe = _element_fe_var[dim][var];
}

inline
FEBase * FEMContext::get_element_fe( unsigned int var, unsigned short dim ) const
{
  libmesh_assert( !_element_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_element_fe_var[dim].size() ) );
  return cast_ptr<FEBase *>( (_element_fe_var[dim][var] ) );
}

template<typename OutputShape>
inline
void FEMContext::get_side_fe( unsigned int var, FEGenericBase<OutputShape> *& fe,
                              unsigned short dim ) const
{
  libmesh_assert( !_side_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
  fe = cast_ptr<FEGenericBase<OutputShape> *>( (_side_fe_var[dim][var] ) );
}

inline
void FEMContext::get_side_fe( unsigned int var, FEAbstract *& fe,
                              unsigned short dim ) const
{
  libmesh_assert( !_side_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
  fe = _side_fe_var[dim][var];
}

inline
FEBase * FEMContext::get_side_fe( unsigned int var, unsigned short dim ) const
{
  libmesh_assert( !_side_fe_var[dim].empty() );
  libmesh_assert_less ( var, (_side_fe_var[dim].size() ) );
  return cast_ptr<FEBase *>( (_side_fe_var[dim][var] ) );
}

template<typename OutputShape>
inline
void FEMContext::get_edge_fe( unsigned int var, FEGenericBase<OutputShape> *& fe ) const
{
  libmesh_assert_less ( var, _edge_fe_var.size() );
  fe = cast_ptr<FEGenericBase<OutputShape> *>( _edge_fe_var[var] );
}

inline
void FEMContext::get_edge_fe( unsigned int var, FEAbstract *& fe ) const
{
  libmesh_assert_less ( var, _edge_fe_var.size() );
  fe = _edge_fe_var[var];
}

inline
FEBase * FEMContext::get_edge_fe( unsigned int var ) const
{
  libmesh_assert_less ( var, _edge_fe_var.size() );
  return cast_ptr<FEBase *>( _edge_fe_var[var] );
}


} // namespace libMesh

#endif // LIBMESH_FEM_CONTEXT_H