lasxgrad.hpp

Go to the documentation of this file.

#ifndef LASXGRAD_HPP
#define LASXGRAD_HPP
 
#include <sxgrad.hpp>
#include <liealgebra.hpp>
#include <unsupported/Eigen/KroneckerProduct>
 
namespace HArDCore3D
{
 
  class LASXGrad : public VariableDOFSpace
  {
  public:
    typedef std::function<double(const Eigen::Vector3d &)> FunctionType;
    typedef std::vector<FunctionType> LAFunctionType;
 
 
    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;
    };
 
    LASXGrad(const LieAlgebra & lie_algebra, const SXGrad & sxgrad, bool use_threads = true, std::ostream & output = std::cout);
 
    const Mesh & mesh() const
    {
      return m_sxgrad.mesh();
    }
    
    const size_t & degree() const
    {
      return m_sxgrad.degree();
    }
 
    inline const LieAlgebra & lieAlg() const
    {
      return m_lie_algebra;
    }
 
    const SerendipityProblem & serPro() const
    {
      return m_sxgrad.serPro();
    }
 
    Eigen::VectorXd interpolate(
          const LAFunctionType & q, 
          const int doe_cell = -1, 
          const int doe_face = -1, 
          const int doe_edge = -1 
          ) const;
 
    //---------------------------------//
    //---- Face transfer operators-----//
    
    inline const Eigen::MatrixXd & SgradFace(size_t iF) const
    {
      return (*m_face_transfer_operators[iF]).serendipity;
    }
 
    inline const Eigen::MatrixXd & EgradFace(size_t iF) const
    {
      return (*m_face_transfer_operators[iF]).extension;
    }
 
    inline const Eigen::MatrixXd & RgradFace(size_t iF) const
    {
      return (*m_face_transfer_operators[iF]).reduction;
    }
 
    inline const Eigen::MatrixXd & SgradFace(const Face & F) const
    {
      return SgradFace(F.global_index());
    }
 
    inline const Eigen::MatrixXd & EgradFace(const Face & F) const
    {
      return EgradFace(F.global_index());
    }
 
    inline const Eigen::MatrixXd & RgradFace(const Face & F) const
    {
      return RgradFace(F.global_index());
    }
    
    //---------------------------------//
    //---- Cell transfer operators -----//
    inline const Eigen::MatrixXd & SgradCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).serendipity;
    }
 
    inline const Eigen::MatrixXd & EgradCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).extension;
    }
 
    inline const Eigen::MatrixXd & RgradCell(size_t iT) const
    {
      return (*m_cell_transfer_operators[iT]).reduction;
    }
 
    inline const Eigen::MatrixXd & SgradCell(const Cell & T) const
    {
      return SgradCell(T.global_index());
    }
 
    inline const Eigen::MatrixXd & EgradCell(const Cell & T) const
    {
      return EgradCell(T.global_index());
    }
 
    inline const Eigen::MatrixXd & RgradCell(const Cell & T) const
    {
      return RgradCell(T.global_index());
    }
 
    //---------------------------------------------------------------------//
    //---- Full gradient and potential reconstructions, and L2 product ----//
    
    inline const Eigen::MatrixXd edgeGradient(size_t iE) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.edgeGradient(iE), id);
    }
    
    inline const Eigen::MatrixXd edgeGradient(const Edge & E) const
    {
      return edgeGradient(E.global_index());
    }
 
    inline const Eigen::MatrixXd edgePotential(size_t iE) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.edgePotential(iE), id);
    }
 
    inline const Eigen::MatrixXd edgePotential(const Edge & E) const
    {
      return edgePotential(E.global_index());
    }
 
    inline const Eigen::MatrixXd faceGradient(size_t iF) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.faceGradient(iF), id);
    }
    
    inline const Eigen::MatrixXd faceGradient(const Face & F) const
    {
      return faceGradient(F.global_index());
    }
 
    inline const Eigen::MatrixXd facePotential(size_t iF) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.facePotential(iF), id);
    }
 
    inline const Eigen::MatrixXd facePotential(const Face & F) const
    {
      return facePotential(F.global_index());
    }
 
    inline const Eigen::MatrixXd cellGradient(size_t iT) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.cellGradient(iT), id);
    }
 
    inline const Eigen::MatrixXd cellGradient(const Cell & T) const
    {
      return cellGradient(T.global_index());
    }
 
    inline const Eigen::MatrixXd cellPotential(size_t iT) const
    {
      Eigen::MatrixXd id = Eigen::MatrixXd::Identity(m_lie_algebra.dimension(), m_lie_algebra.dimension());
      return Eigen::KroneckerProduct(m_sxgrad.cellPotential(iT), id);
    }
 
    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_Pk3_T = Eigen::MatrixXd::Zero(1,1), 
                                     const IntegralWeight & weight = IntegralWeight(1.) 
                                     ) const
    {
      return Eigen::KroneckerProduct(m_sxgrad.computeL2Product(iT,penalty_factor,mass_Pk3_T, weight), m_lie_algebra.massMatrix()); 
    }
 
         
    //-----------------------//       
    //---- Getters ----------//
 
    inline const DDRCore::CellBases & cellBases(size_t iT) const
    {
      return m_sxgrad.cellBases(iT);
    }
 
    inline const DDRCore::CellBases & cellBases(const Cell & T) const
    {
      return m_sxgrad.cellBases(T.global_index());
    }
    
    inline const DDRCore::FaceBases & faceBases(size_t iF) const
    {
      return m_sxgrad.faceBases(iF);
    }
 
    inline const DDRCore::FaceBases & faceBases(const Face & F) const
    {
      return m_sxgrad.faceBases(F.global_index());
    }
    
    inline const DDRCore::EdgeBases & edgeBases(size_t iE) const
    {
      return m_sxgrad.edgeBases(iE);
    }
 
    inline const DDRCore::EdgeBases & edgeBases(const Edge & E) const
    {
      return m_sxgrad.edgeBases(E.global_index());
    }
 
  private:
    TransferOperators _compute_face_transfer_operators(size_t iF);
    TransferOperators _compute_cell_transfer_operators(size_t iT);
    
    const LieAlgebra & m_lie_algebra;
    const SXGrad & m_sxgrad;
       
    // Containers for serendipity, extension and reduction operators
    std::vector<std::unique_ptr<TransferOperators> > m_face_transfer_operators;
    std::vector<std::unique_ptr<TransferOperators> > m_cell_transfer_operators;
 
    bool m_use_threads;
    std::ostream & m_output;
 
  };
 
} // end of namespace HArDCore3D
#endif