sxcurl.hpp

Go to the documentation of this file.

#ifndef SXCURL_HPP
#define SXCURL_HPP
 
#include <variabledofspace.hpp>
#include <ddrcore.hpp>
#include <integralweight.hpp>
#include <xcurl.hpp>
#include <sxgrad.hpp>
#include <serendipity_problem.hpp>
 
namespace HArDCore2D
{
 
  class SXCurl : public VariableDOFSpace
  {
  public:
    typedef std::function<Eigen::Vector2d(const Eigen::Vector2d &)> FunctionType;
 
 
    struct TransferOperators
    {
      TransferOperators(
                     const Eigen::MatrixXd & _serendipity,    
                     const Eigen::MatrixXd & _extension,     
                     const Eigen::MatrixXd & _reduction      
                     )
        : serendipity(_serendipity),
          extension(_extension),
          reduction(_reduction)
      {
        // Do nothing
      }
 
      Eigen::MatrixXd serendipity;
      Eigen::MatrixXd extension;
      Eigen::MatrixXd reduction;
    };
 
    SXCurl(const DDRCore & ddr_core, const SerendipityProblem & ser_pro, 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;
    //----------------------------------//
    //---- Cell transfer operators -----//
    inline const Eigen::MatrixXd & ScurlCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).serendipity;
    }
 
    inline const Eigen::MatrixXd & EcurlCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).extension;
    }
 
    inline const Eigen::MatrixXd & RcurlCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).reduction;
    }
 
    inline const Eigen::MatrixXd & ScurlCell(const Cell & T) const
    {
      return ScurlCell(T.global_index());
    }
 
    inline const Eigen::MatrixXd & EcurlCell(const Cell & T) const
    {
      return EcurlCell(T.global_index());
    }
 
    inline const Eigen::MatrixXd & RcurlCell(const Cell & T) const
    {
      return RcurlCell(T.global_index());
    }
 
    //-----------------------------------------------------------------//
    //---- Full curl and potential reconstructions, and L2 product ----//
 
    inline const Eigen::MatrixXd cellCurl(size_t iT) const
    {
      return m_xcurl.cellOperators(iT).curl * EcurlCell(iT);
    }
 
    inline const Eigen::MatrixXd cellCurl(const Cell & T) const
    {
      return cellCurl(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellPotential(size_t iT) const
    {
      return m_xcurl.cellOperators(iT).potential * EcurlCell(iT);
    }
 
    inline const Eigen::MatrixXd cellPotential(const Cell & T) const
    {
      return cellPotential(T.global_index());
    }
 
    Eigen::MatrixXd computeStabilisation(
                                 const size_t iT,
                                 const double & penalty_factor = 1.,
                                 const IntegralWeight & weight = IntegralWeight(1.)
                                 ) const;
                                 
    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
    {
      return EcurlCell(iT).transpose() 
                * m_xcurl.computeL2Product(iT, penalty_factor, mass_Pk2_T, weight)
                  *EcurlCell(iT);  
    }
 
         
    Eigen::MatrixXd computeL2ProductGradient(
                                     const size_t iT, 
                                     const SXGrad & sx_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;
 
    //-----------------------//       
    //---- Getters ----------//
 
    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());
    }
 
  private:
    TransferOperators _compute_cell_transfer_operators(size_t iT);
    
    const DDRCore & m_ddr_core;
    const SerendipityProblem & m_ser_pro;
    const XCurl m_xcurl;
       
    // Containers for serendipity, extension and reduction operators
    std::vector<std::unique_ptr<TransferOperators> > m_cell_transfer_operators;
 
    bool m_use_threads;
    std::ostream & m_output;
 
  };
 
} // end of namespace HArDCore2D
#endif