vsxcurl.hpp

Go to the documentation of this file.

#ifndef VSXCURL_HPP
#define VSXCURL_HPP
 
#include <sxcurl.hpp>
#include <unsupported/Eigen/KroneckerProduct>
 
 
namespace HArDCore2D
{
     class VSXCurl : public VariableDOFSpace
  {
  public:
    // Type of continuous functions corresponding to the discrete space VScurl
    typedef std::function<Eigen::Matrix2d(const Eigen::Vector2d &)> FunctionType;
    // Types for element bases
    typedef Family<MonomialScalarBasisCell> PolyBasisCellType;
    typedef TensorizedVectorFamily<DDRCore::PolyBasisCellType, dimspace> Poly2BasisCellType;
    typedef MatrixFamily<DDRCore::PolyBasisCellType, dimspace> Poly2x2BasisCellType;
  
    // Type for edge basis
    typedef Family<MonomialScalarBasisEdge> PolyBasisEdgeType;
    typedef TensorizedVectorFamily<DDRCore::PolyBasisEdgeType, dimspace> Poly2BasisEdgeType;
 
 
        struct CellBases
    { 
      std::unique_ptr<Poly2x2BasisCellType> Polyk2x2;
    };
 
    VSXCurl(const SerendipityProblem & sp, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_sxcurl.mesh();
    }
    
    const SXCurl & sXcurl() const
    {
      return m_sxcurl;
    }
    
    const size_t & degree() const
    {
      return m_sxcurl.degree();
    }
 
    Eigen::VectorXd interpolate(
          const FunctionType & q, 
          const int doe_cell = -1 
          ) const;
 
 
    inline auto Polyk(size_t iT) const -> const PolyBasisCellType& {
      return *m_sxcurl.cellBases(iT).Polyk;
    }
    
    //---------------------------------------------------------------------//
    //---- Curl and potential reconstructions, and stabilisation -----//
       inline const Eigen::MatrixXd cellCurl(size_t iT) const
    {
      return PermuteTensorization(m_sxcurl.cellBases(iT).Polyk->dimension(), dimspace) 
                * Eigen::KroneckerProduct(m_sxcurl.cellCurl(iT), Eigen::Matrix2d::Identity());
    }
 
    inline const Eigen::MatrixXd cellCurl(const Cell & T) const
    {
      return cellCurl(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellPotential(size_t iT) const
    {
      Eigen::MatrixXd PotPermutedBasis 
          = PermuteTensorization(dimspace, m_sxcurl.cellBases(iT).Polyk->dimension()) * m_sxcurl.cellPotential(iT);
      return PermuteTensorization( m_sxcurl.cellBases(iT).Polyk->dimension(), dimspace*dimspace) 
                * Eigen::KroneckerProduct(PotPermutedBasis , Eigen::Matrix2d::Identity());
    }
 
    inline const Eigen::MatrixXd cellPotential(const Cell & T) const
    {
      return cellPotential(T.global_index());
    }
 
    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 Eigen::KroneckerProduct(m_sxcurl.computeL2Product(iT, penalty_factor, mass_Pk2_T, weight), Eigen::Matrix2d::Identity()); 
    }
    
 
        //-----------------------//       
    //---- Getters ----------//
 
    inline const CellBases & cellBases(size_t iT) const
    {
      assert( m_cell_bases[iT] );
      return *m_cell_bases[iT].get();
    }
 
    inline const CellBases & cellBases(const Cell & T) const
    {
      return cellBases(T.global_index());
    }
    private:
/*     const DDRCore & m_ddr_core;
    const SerendipityProblem& m_sp; // Warning */
    const SXCurl m_sxcurl; // Not a reference as it is not created outside this class
    CellBases _compute_cell_bases(size_t iT);
 
        // Cell bases
    std::vector<std::unique_ptr<CellBases> > m_cell_bases;
 
 
    bool m_use_threads;
    std::ostream & m_output;
    
  };
 
} // end of namespace HArDCore2D
#endif