vsxgrad.hpp

Go to the documentation of this file.

#ifndef VSXGRAD_HPP
#define VSXGRAD_HPP
 
#include <sxgrad.hpp>
#include <sxcurl.hpp>
#include <unsupported/Eigen/KroneckerProduct>
#include <globaldofspace.hpp>
#include <vsxcurl.hpp>
 
 
namespace HArDCore2D
{
 
  class VSXGrad : public VariableDOFSpace
  {
  public:
    // Type of continuous functions corresponding to the discrete space VSXgrad
    typedef std::function<Eigen::Vector2d(const Eigen::Vector2d &)> FunctionType;
 
    struct LocalOperators
    {
      LocalOperators(
                     const Eigen::MatrixXd & _rot,     
                     const Eigen::MatrixXd & _potential 
                     )
        : rot(_rot),
          potential(_potential)
      {
        // Do nothing
      }
 
      Eigen::MatrixXd rot;
      Eigen::MatrixXd potential;
    };
 
 
    // Types for element bases
    typedef Family<MonomialScalarBasisCell> PolyBasisCellType;
    typedef TensorizedVectorFamily<DDRCore::PolyBasisCellType, dimspace> Poly2BasisCellType;
    typedef MatrixFamily<DDRCore::PolyBasisCellType, dimspace> Poly2x2BasisCellType;
    typedef Family<CurlBasis<ShiftedBasis<MonomialScalarBasisCell> > > RolyBasisCellType;
    typedef Family<RolyComplBasisCell> RolyComplBasisCellType;  
    // Type for edge basis
    typedef Family<MonomialScalarBasisEdge> PolyBasisEdgeType;
    typedef TensorizedVectorFamily<DDRCore::PolyBasisEdgeType, dimspace> Poly2BasisEdgeType;
 
    // Structure to store vectorial/matrix cell bases (the scalar ones are in m_sxgrad)
    struct CellBases
    { 
      std::unique_ptr<Poly2BasisCellType> Polykmo2;
      std::unique_ptr<Poly2BasisCellType> Polykptwo2;
      std::unique_ptr<Poly2x2BasisCellType> Polykpo2x2;
    };
  
    // Structure to store vectorial/matrix edge bases (the scalar ones are in m_sxgrad)
    struct EdgeBases
    {
      std::unique_ptr<Poly2BasisEdgeType> Polykmo2;
      std::unique_ptr<Poly2BasisEdgeType> Polyk2;
      std::unique_ptr<Poly2BasisEdgeType> Polykpo2;
      std::unique_ptr<Poly2BasisEdgeType> Polykptwo2;
    };
 
    VSXGrad(const DDRCore & ddr_core, const SerendipityProblem & sp, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_sxgrad.mesh();
    }
    
    const SXGrad & sXgrad() const
    {
      return m_sxgrad;
    }
    
    const size_t & degree() const
    {
      return m_sxgrad.degree();
    }
    
    Eigen::VectorXd interpolate(
          const FunctionType & q, 
          const int doe_cell = -1 
          ) const;
    
    //---------------------------------------------------------------------//
    //---- Gradients and potential reconstructions, and stabilisation -----//
    
    /* Operators from sXgrad are easily extended into operators on vsXgrad using the Kronecker product
     These extension corresponds to the abovementioned tensorized basis of sXgrad and, for the codomain,
     to the tensorized basis of R^2 \otimes P (for P the polynomial space in which the operator of sXgrad
     lands). We need to do a change of basis of that codomain to represent the operators in the bases
     stored in the structures EdgeBases, FaceBases and CellBases.
    */
    
    //Rename spaces so that their name fit with the fact that we will take degree K+1
 
    
    inline auto PolykE(size_t iE) const -> const PolyBasisEdgeType& {
      return *m_sxgrad.edgeBases(iE).Polykmo;
    }
    inline auto PolykpoE(size_t iE) const -> const PolyBasisEdgeType& {
      return *m_sxgrad.edgeBases(iE).Polyk;
    }
 
    inline auto PolykptwoE(size_t iE) const -> const PolyBasisEdgeType& {
      return *m_sxgrad.edgeBases(iE).Polykpo;
    }
 
 
    inline auto Polykpo(size_t iT) const -> const PolyBasisCellType& {
      return *m_sxgrad.cellBases(iT).Polyk;
    }
 
        
    inline auto Polykmo(size_t iT) const -> const PolyBasisCellType& {
      return *m_sxgrad.cellBases(iT).Polykmtwo;
    }
 
    
    inline auto Polyk(size_t iT) const -> const PolyBasisCellType& {
      return *m_sxgrad.cellBases(iT).Polykmo;
    }
  
    
    inline auto Polykptwo(size_t iT) const -> const PolyBasisCellType& {
      return *m_sxgrad.cellBases(iT).Polykpo;
    }
 
    
    inline auto Polyk2(size_t iT) const -> const Poly2BasisCellType& {
      return *m_sxgrad.cellBases(iT).Polykmo2;
    }
 
    inline auto Rolyk(size_t iT) const -> const RolyBasisCellType& {
      return *m_sxgrad.cellBases(iT).Rolykmo;
    }
 
        inline auto RolyComplk(size_t iT) const -> const RolyComplBasisCellType& {
      return *m_sxgrad.cellBases(iT).RolyComplkmo;
    }
 
 
    
    inline const Eigen::MatrixXd edgeGradient(size_t iE) const
    {
      // The codomain for the edge gradient of sXgrad is Polyk, so we need to change from
      // R^3 \otimes Polyk to Polyk3 (which is, by construction, Polyk \otimes R^3)
      return PermuteTensorization(m_sxgrad.edgeBases(iE).Polyk->dimension(), dimspace) 
                * Eigen::KroneckerProduct(m_sxgrad.edgeGradient(iE), Eigen::Matrix2d::Identity());
    }
    
    inline const Eigen::MatrixXd edgeGradient(const Edge & E) const
    {
      return edgeGradient(E.global_index());
    }
 
    inline const Eigen::MatrixXd edgePotential(size_t iE) const
    {
      // The codomain for the edge potential of sXgrad is Polykpo, so we need to change from
      // R^3 \otimes Polykpo to Polykpo3 which is, by construction, Polykpo \otimes R^3
      return PermuteTensorization(m_sxgrad.edgeBases(iE).Polykpo->dimension(), dimspace) 
                * Eigen::KroneckerProduct(m_sxgrad.edgePotential(iE), Eigen::Matrix2d::Identity());
    }
 
    inline const Eigen::MatrixXd edgePotential(const Edge & E) const
    {
      return edgePotential(E.global_index());
    }
 
 
    inline const Eigen::MatrixXd cellGradient(size_t iT) const
    {
      // The codomain for the cell gradient of sXgrad is Polyk3, so we need to change from
      // R^3 \otimes Polyk3 to Polyk3x3. 
      // Here, Polyk3 = Polyk \otimes R^3, where the R^3=R^3_nabla component represents the derivatives.
      // We first permute the codomain of m_sxgrad.cellGradient(iT) to make it R^3_nabla \otimes Polyk.
      Eigen::MatrixXd GradPermutedBasis 
          = PermuteTensorization(dimspace, m_sxgrad.cellBases(iT).Polyk->dimension()) * m_sxgrad.cellGradient(iT);
 
      // The Kronecker product then gives the cellGradient of vsXgrad in (R^3 \otimes R^3_nabla) \otimes Polyk,
      // where the first copy of R^3 lists the components of the functions in the vectorized vsXgrad.
      // We permute this basis to Polyk \otimes (R^3 \otimes R^3_nabla). In this basis, the first index
      // of R^3 \otimes R^3_nabla goes over the components of the functions, and the second over the coordinates
      // derivatives. 
      // Polyk3x3 is Polyk \otimes M_N(R) with the basis of M_N(R) listing the canonical
      // basis vectors in columns; so identifying the basis of R^3 \otimes R^3_nabla with the basis of M_N(R)
      // organises the functions components in columns and their gradients in rows, which is exactly what
      // we want. Hence, the matrix obtained on Polyk \otimes (R^3 \otimes R^3_nabla) can directly be interpreted
      // as the matrix on Polyk \otimes M_N(R) = Polyk3x3.
      return PermuteTensorization(m_sxgrad.cellBases(iT).Polyk->dimension(), dimspace*dimspace) 
                * Eigen::KroneckerProduct(GradPermutedBasis, Eigen::Matrix2d::Identity());
    }
 
    inline const Eigen::MatrixXd cellGradient(const Cell & T) const
    {
      return cellGradient(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellSymGradient(size_t iT) const
    {
      return cellBases(iT).Polykpo2x2->symmetriseOperator() * cellGradient(iT);
    }
 
    inline const Eigen::MatrixXd cellSymGradient(const Cell & T) const
    {
      return cellSymGradient(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellDivergence(size_t iT) const
    {
      return cellBases(iT).Polykpo2x2->traceOperator() * cellGradient(iT);
    }
 
    inline const Eigen::MatrixXd cellDivergence(const Cell & T) const
    {
      return cellDivergence(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellPotential(size_t iT) const
    {
      // The codomain for the cell potential of sXgrad is Polykpo, so we need to change from
      // R^3 \otimes Polykpo to Polykpo3 which is, by construction, Polykpo \otimes R^3
      return PermuteTensorization(m_sxgrad.cellBases(iT).Polykpo->dimension(), dimspace) 
                * Eigen::KroneckerProduct(m_sxgrad.cellPotential(iT), Eigen::Matrix2d::Identity());
    }
 
    inline const Eigen::MatrixXd cellPotential(const Cell & T) const
    {
      return cellPotential(T.global_index());
    }
    
    Eigen::MatrixXd computeGradientFull(
                                    const size_t iT,
                                    const VSXCurl & vsx_curl
                                    ) const;
 
    Eigen::MatrixXd _permute_bases(size_t iT) const;
 
    Eigen::MatrixXd _injection(size_t iT);
 
    Eigen::MatrixXd computeStabilisation(
                                     const size_t iT, 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const
    {
      // The H^1-stabilisation is the L^2-stabilisation divided by h_T^2
      double hTm2 = std::pow(mesh().cell(iT)->diam(), -2);
      return hTm2 * Eigen::KroneckerProduct(m_sxgrad.computeStabilisation(iT, weight), Eigen::Matrix2d::Identity()); 
    }
 
    std::vector<double> computeH1norms(
                          const std::vector<Eigen::VectorXd> & list_dofs, 
                          const std::string typegrad="full" 
                          ) const;
         
    Eigen::MatrixXd computeL2Stabilisation(
                                     const size_t iT, 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const
    {
      return Eigen::KroneckerProduct(m_sxgrad.computeStabilisation(iT, weight), Eigen::Matrix2d::Identity()); 
    }
 
    Eigen::MatrixXd computeScalarProduct(
                                     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 ;
 
    std::vector<VectorRd> computeVertexValues(
                  const Eigen::VectorXd & u   
                  ) const;
         
    //-----------------------//       
    //---- 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());
    }
    
    inline const EdgeBases & edgeBases(size_t iE) const
    {
      assert( m_edge_bases[iE] );
      return *m_edge_bases[iE].get();
    }
 
    inline const EdgeBases & edgeBases(const Edge & E) const
    {
      return edgeBases(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()];
    }
 
  private:
    EdgeBases _compute_edge_bases(size_t iE);
    CellBases _compute_cell_bases(size_t iT);
    
    const DDRCore & m_ddr_core;
    const SerendipityProblem& m_sp; 
    const SXGrad m_sxgrad; // Not a reference as it is not created outside this class
    const DDRCore m_ddrcore_for_sxcurl;   
    const SerendipityProblem m_sp_for_sxcurl;        
    const SXCurl m_sxcurl; // Not a reference as it is not created outside this class
 
    // Cell bases
    std::vector<std::unique_ptr<CellBases> > m_cell_bases;
    // Edge bases
    std::vector<std::unique_ptr<EdgeBases> > m_edge_bases;
 
    bool m_use_threads;
    std::ostream & m_output;
    
    LocalOperators  _compute_cell_operators(size_t iT);
 
        // Containers for local operators
    std::vector<std::unique_ptr<LocalOperators> > m_cell_operators;
    
  };
 
} // end of namespace HArDCore2D
#endif