excurl.hpp

Go to the documentation of this file.

#ifndef EXCURL_HPP
#define EXCURL_HPP
 
#include <xcurl.hpp>
 
namespace HArDCore2D
{
 
  class EXCurl : public GlobalDOFSpace
  {
  public:
    typedef std::function<Eigen::Vector2d(const Eigen::Vector2d &)> FunctionType;
    typedef MatrixFamily<DDRCore::PolyBasisCellType, dimspace> Polyk2x2Type;
    typedef TensorizedVectorFamily<DDRCore::PolyBasisCellType, dimspace> Polykpo2Type;
    
    struct hhoLocalOperators
    {
      hhoLocalOperators(
                     const Eigen::MatrixXd & _gradient,     
                     const Eigen::MatrixXd & _symmetric_gradient,     
                     const Eigen::MatrixXd & _potential, 
                     const Eigen::MatrixXd & _stabilisation 
                     )
        : gradient(_gradient),
          symmetric_gradient(_symmetric_gradient),
          potential(_potential),
          stabilisation(_stabilisation)
      {
        // Do nothing
      }
 
      Eigen::MatrixXd gradient;
      Eigen::MatrixXd symmetric_gradient;
      Eigen::MatrixXd potential;
      Eigen::MatrixXd stabilisation;
    };
 
    EXCurl(const DDRCore & ddr_core, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_ddr_core.mesh();
    }
    
    const size_t & degree() const
    {
      return m_ddr_core.degree();
    }
 
    Eigen::VectorXd interpolate(
                const FunctionType & v, 
                const int deg_quad = -1 
                ) const;
 
    // We don't store the potential and curl, as they are very cheap to construct from the (stored) operators in m_xcurl
    inline const Eigen::MatrixXd ddrPotential(size_t iT) const
    {
      return m_xcurl.cellOperators(iT).potential * m_reduction[iT];
    };
 
    inline const Eigen::MatrixXd ddrCurl(size_t iT) const
    {
      return m_xcurl.cellOperators(iT).curl * m_reduction[iT];
    };
 
    inline const Eigen::MatrixXd ddrPotential(const Cell & T) const
    {
      return ddrPotential(T.global_index());
    };
 
    inline const Eigen::MatrixXd ddrCurl(const Cell & T) const
    {
      return ddrCurl(T.global_index());
    };
 
    inline const hhoLocalOperators & hhoOperators(size_t iT) const
    {
      assert( m_hho_operators[iT] );
      return *m_hho_operators[iT];
    }
 
    inline const hhoLocalOperators & hhoOperators(const Cell & T) const
    {
      return * m_hho_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 Polyk2x2Type & Polyk2x2(size_t iT) const
    {
      // Make sure that the basis has been created
      assert( m_Polyk2x2[iT] );
      return *m_Polyk2x2[iT].get();
    }
 
    inline const Polykpo2Type & Polykpo2(size_t iT) const
    {
      // Make sure that the basis has been created
      assert( m_Polykpo2[iT] );
      return *m_Polykpo2[iT].get();
    }
 
    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(T)^2 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 computeL2Product(
                                     const size_t iT, 
                                     const double & penalty_factor = 1., 
                                     const Eigen::MatrixXd & mass_Pk2_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const;
                                     
    Eigen::MatrixXd computeL2ProductGradient(
                                     const size_t iT, 
                                     const XGrad & x_grad, 
                                     const std::string & side, 
                                     const double & penalty_factor = 1., 
                                     const Eigen::MatrixXd & mass_Pk2_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const;
 
    Eigen::MatrixXd computeL2Product_with_Ops(
                                     const size_t iT, 
                                     const std::vector<Eigen::MatrixXd> & leftOp, 
                                     const std::vector<Eigen::MatrixXd> & rightOp, 
                                     const double & penalty_factor, 
                                     const Eigen::MatrixXd & w_mass_Pk2_T, 
                                     const IntegralWeight & weight 
                                     ) const;
                                
 
  private:
    hhoLocalOperators _compute_hho_operators(size_t iT);
      
    const DDRCore & m_ddr_core;
    const XCurl m_xcurl; // Underlying standard XCurl space
    std::vector<std::unique_ptr<Polyk2x2Type> > m_Polyk2x2;
    std::vector<std::unique_ptr<Polykpo2Type> > m_Polykpo2;
    std::vector<Eigen::MatrixXd> m_reduction; // reduction matrix to go from EXcurl DOFs to XCurl DOFs 
    bool m_use_threads;
    std::ostream & m_output;
    
    // Containers for local hho operators
    std::vector<std::unique_ptr<hhoLocalOperators> > m_hho_operators;
  };
 
} // end of namespace HArDCore2D
#endif