xhess_full.hpp

Go to the documentation of this file.

#ifndef XHESS_HPP
#define XHESS_HPP
 
#include <globaldofspace.hpp>
#include <ddrcore.hpp>
#include <integralweight.hpp>
#include "serendipity_problem.hpp"
#include "xgrad.hpp"
#include "sxgrad.hpp"
#include "vsxgrad.hpp"
 
namespace HArDCore2D
{
  class XHess : public GlobalDOFSpace
  {
  public:
    typedef std::function<double(const Eigen::Vector2d &)> FunctionType;
    typedef std::function<Eigen::Vector2d(const Eigen::Vector2d &)> GradFunctionType;
 
 
    struct LocalOperators
    {
      LocalOperators(
                     const Eigen::MatrixXd & _gradient, 
                     const Eigen::MatrixXd & _potential 
                     )
        : gradient(_gradient),
          potential(_potential)
      {
        // Do nothing
      }
      
      Eigen::MatrixXd gradient;
      Eigen::MatrixXd potential;
    };
    struct TransferOperators
    {
      TransferOperators(
                     const Eigen::MatrixXd & _serendipity,   
                     const Eigen::MatrixXd & _extension   
                     )
        : serendipity(_serendipity),
          extension(_extension)
      {
        // Do nothing
      }
 
      Eigen::MatrixXd serendipity;
      Eigen::MatrixXd extension;
    };
    Eigen::MatrixXd SerExtension(size_t iT) const;
    XHess(const DDRCore & ddr_core,  const SerendipityProblem & ser_pro, bool use_threads = true, std::ostream & output = std::cout);//ajout senredipity ser_pro
 
    const Mesh & mesh() const
    {
      return m_ddr_core.mesh();
    }
    
    const size_t & degree() const
    {
      return m_ddr_core.degree();
    }
    
    Eigen::VectorXd interpolate(
                                const FunctionType & q, 
                                const GradFunctionType & Dq, 
                                const int deg_quad = -1 
                                ) const;
                                
        //---------------------------------//
    //---- Cell transfer operators -----//
    inline const Eigen::MatrixXd & SgradCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).serendipity;
    }
 
    inline const Eigen::MatrixXd & SgradCell(const Cell & T) const
    {
      return SgradCell(T.global_index());
    }
 
    inline const Eigen::MatrixXd & EgradCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).extension;
    }
 
    inline const Eigen::MatrixXd & EgradCell(const Cell & T) const
    {
      return EgradCell(T.global_index());
    }
 
    inline const TransferOperators & TcellOperators(size_t iT) const
    {
      return *m_cell_transfer_operators[iT];
    }
 
    inline const TransferOperators & TcellOperators(const Cell & T) const
    {
      return *m_cell_transfer_operators[T.global_index()];
    }
 
    inline const LocalOperators & edgeOperators(size_t iE) const
    {
      return *m_edge_operators[iE];
    }
 
    inline const LocalOperators & edgeOperators(const Edge & E) const
    {
      return *m_edge_operators[E.global_index()];
    }
    
    inline const LocalOperators & cellOperators(size_t iT) const
    {
      return *m_cell_operators[iT];
    }
 
    inline const LocalOperators & cellOperators(const Cell & T) const
    {
      return *m_cell_operators[T.global_index()];
    }
 
    inline const DDRCore::CellBases & cellBases(size_t iT) const
    {
      return m_ddr_core.cellBases(iT);
    }
 
    inline const DDRCore::CellBases & cellBases(const Cell & T) const
    {
      return m_ddr_core.cellBases(T.global_index());
    }
    
    inline const DDRCore::EdgeBases & edgeBases(size_t iE) const
    {
      return m_ddr_core.edgeBases(iE);
    }
 
    inline const DDRCore::EdgeBases & edgeBases(const Edge & E) const
    {
      return m_ddr_core.edgeBases(E.global_index());
    }
 
    // The mass matrix of P^{k+1}(T) is the most expensive mass matrix in the calculation of this norm, which
    // is why there's the option of passing it as parameter if it's been already pre-computed when the norm is called.
 
    Eigen::MatrixXd computeStabilisation(
                                     const size_t iT, 
                                     const double & penalty_factor = 1., 
                                     const Eigen::MatrixXd & mass_Pkpo_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const;
 
    Eigen::MatrixXd computeL2Product(
                                     const size_t iT, 
                                     const double & penalty_factor = 1., 
                                     const Eigen::MatrixXd & mass_Pkpo_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const;
 
 
 
    Eigen::MatrixXd computeGradientFull(
                                     const size_t iT, 
                                     const VSXGrad & vsx_grad
    ) const;
 
    Eigen::MatrixXd potential2(
                                     const size_t iT, 
                                     const VSXGrad & vsx_grad
    ) const;
 
    double evaluatePotential(
                             const size_t iT, 
                             const Eigen::VectorXd & vT, 
                             const VectorRd & x 
                             ) const;
 
    double computeL2Norm(const Eigen::VectorXd & v) const;    
    
  private:    
    Eigen::MatrixXd _compute_cell_serendipity_operator(size_t iT); 
    Eigen::MatrixXd _compute_cell_extension_operator(size_t iT) const; 
 
    LocalOperators _compute_edge_gradient_potential(size_t iE);
    LocalOperators _compute_cell_gradient_potential(size_t iT);
    double _compute_squared_l2_norm(size_t iT, const Eigen::VectorXd & vT) const;
    
    const DDRCore & m_ddr_core;
    const SerendipityProblem & m_ser_pro;  
    const SXGrad m_sxgrad;                  
    
    bool m_use_threads;
    std::ostream & m_output;
 
    // Containers for local operators
    std::vector<std::unique_ptr<LocalOperators> > m_edge_operators;
    std::vector<std::unique_ptr<LocalOperators> > m_cell_operators;
    
    // Containers for serendipity
    std::vector<std::shared_ptr<TransferOperators> > m_cell_transfer_operators; 
  };
 
} // end of namespace HArDCore2D
 
#endif