serendipity_problem.hpp

Go to the documentation of this file.

 
#ifndef SERENDIPITY_PROBLEM_HPP
#define SERENDIPITY_PROBLEM_HPP
 
#include <ddrcore.hpp>
 
namespace HArDCore2D
{
 
  //------------------------------------------------------------------------------
 
  class SerendipityProblem
  {
  public:
    typedef RestrictedBasis<DDRCore::PolyBasisCellType > PolylBasisCellType;
    typedef RestrictedBasis<DDRCore::RolyComplBasisCellType > RolyCompllpoBasisCellType;
    
    /* PartialPivLU does not seem to always work... */
    typedef Eigen::FullPivLU<Eigen::MatrixXd> InverseProblem;
    
    SerendipityProblem(const DDRCore & ddrcore, bool use_threads = true, std::ostream & output = std::cout);
    
    //---------- Serendipity edges and faces, and degrees -------------//
    inline const std::vector<size_t> & serendipityEdges(size_t iT) const
    {
      return m_serendipity_edges[iT];
    }
    
    inline const int n_serendipityEdges(size_t iT) const
    {
      return m_serendipity_edges[iT].size();
    }
 
    inline const int serDegreeCell(size_t iT) const
    {
      return m_ddrcore.degree() + 1 - n_serendipityEdges(iT);
    }
 
 
    //---------- Serendipity dimensions and bases for grad -------------//
    inline size_t dimCellPolyl(size_t iT) const
    {
      return PolynomialSpaceDimension<Cell>::Poly(serDegreeCell(iT));
    }
 
    inline size_t dimCellPolyl(const Cell & T) const
    {
      return dimCellPolyl(T.global_index());
    }
 
    inline const PolylBasisCellType & cellBasisPolyl(size_t iT) const
    {
      // Make sure that the basis has been created
      assert( m_cell_bases_Polyl[iT] );
      return *m_cell_bases_Polyl[iT].get();
 
    }
 
    inline const PolylBasisCellType & cellBasisPolyl(const Cell & T) const
    {
      return cellBasisPolyl(T.global_index()); 
    }
   
    Eigen::VectorXd nDOFs_cells_SXGrad() const;
 
    //---------- Serendipity dimensions and bases for curl -------------//
 
    inline size_t dimCellRolyCompllpo(size_t iT) const
    {
      return PolynomialSpaceDimension<Cell>::RolyCompl(serDegreeCell(iT)+1);
    }
 
    inline size_t dimCellRolyCompllpo(const Cell & T) const
    {
      return dimCellRolyCompllpo(T.global_index());
    }
 
    inline const RolyCompllpoBasisCellType & cellBasisRolyCompllpo(size_t iT) const
    {
      // Make sure that the basis has been created
      assert( m_cell_bases_RolyCompllpo[iT] );
      return *m_cell_bases_RolyCompllpo[iT].get();
 
    }
 
    inline const RolyCompllpoBasisCellType & cellBasisRolyCompllpo(const Cell & T) const
    {
      return cellBasisRolyCompllpo(T.global_index()); 
    }
    
    Eigen::VectorXd nDOFs_cells_SXCurl() const;
    
    
    
    //---------- Serendipity dimensions and bases for rotrot -------------//
    
 
    
    Eigen::VectorXd nDOFs_cells_SXRotRot() const;
 
 
 
    //------------ Serendipity operators ---------------//
 
    const Eigen::MatrixXd SerendipityOperatorCell(const size_t iT, const Eigen::MatrixXd & LT) const;
    
    inline const Eigen::MatrixXd SerendipityOperatorCell(const Cell & T, const Eigen::MatrixXd & LT) const
    {
      return SerendipityOperatorCell(T.global_index(), LT);
    };
 
 
    inline const Mesh & mesh() const
    {
      return m_ddrcore.mesh();
    }
 
    inline const DDRCore & ddrCore() const
    {
      return m_ddrcore;
    }
 
  private:
 
    std::vector<size_t> _compute_serendipity_edges(size_t iT);
 
    InverseProblem _compute_inverse_problem_cells(size_t iT);
 
    // Pointer do DDRCore
    const DDRCore & m_ddrcore;
    // Output stream
    std::ostream & m_output;
 
    // cell basis P^l
    std::vector<std::unique_ptr<PolylBasisCellType> > m_cell_bases_Polyl;
 
    // cell basis R^{c,l+1}
    std::vector<std::unique_ptr<RolyCompllpoBasisCellType> > m_cell_bases_RolyCompllpo;
 
    // List the local indices of edges/faces used for the serendipity on each face/cell
    std::vector<std::vector<size_t>> m_serendipity_edges;  
  
    // Store the inverse of the matrix to compute serendipity operators on cells
    std::vector<InverseProblem> m_inverse_problem_cells;
    
  };
  
} // end of namespace HArDCore2D
 
#endif // SERENDIPITY_PROBLEM_HPP