serendipity_problem.hpp

Go to the documentation of this file.

// Class to select edges/faces to create serendipity spaces and compute serendipity operator
//
//
// Author: Jerome Droniou (jerome.droniou@monash.edu)
//
 
/*
 *
 *      This library was developed around HHO methods, although some parts of it have a more
 * general purpose. If you use this code or part of it in a scientific publication, 
 * please mention the following book as a reference for the underlying principles
 * of HHO schemes:
 *
 * The Hybrid High-Order Method for Polytopal Meshes: Design, Analysis, and Applications. 
 *  D. A. Di Pietro and J. Droniou. Modeling, Simulation and Applications, vol. 19. 
 *  Springer International Publishing, 2020, xxxi + 525p. doi: 10.1007/978-3-030-37203-3. 
 *  url: https://hal.archives-ouvertes.fr/hal-02151813.
 *
 */
 
 
#ifndef SERENDIPITY_PROBLEM_HPP
#define SERENDIPITY_PROBLEM_HPP
 
#include <ddrcore.hpp>
 
namespace HArDCore3D
{
 
  //------------------------------------------------------------------------------
 
  class SerendipityProblem
  {
  public:
    typedef RestrictedBasis<DDRCore::PolyBasisFaceType> PolylBasisFaceType;
    typedef RestrictedBasis<DDRCore::PolyBasisCellType> PolylBasisCellType;
    typedef RestrictedBasis<DDRCore::RolyComplBasisFaceType> RolyCompllpoBasisFaceType;
    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,   
                       bool use_serendipity = true  
                       );
    
    //---------- Serendipity edges and faces, and degrees -------------//
    inline const std::vector<size_t> & serendipityEdges(size_t iF) const
    {
      return m_serendipity_edges[iF];
    }
    
    inline const int n_serendipityEdges(size_t iF) const
    {
      return m_serendipity_edges[iF].size();
    }
 
    inline const int serDegreeFace(size_t iF) const
    {
      return m_ddrcore.degree() + 1 - n_serendipityEdges(iF);
    }
 
    inline const std::vector<size_t> & serendipityFaces(size_t iT) const
    {
      return m_serendipity_faces[iT];
    }
    
    inline const int n_serendipityFaces(size_t iT) const
    {
      return m_serendipity_faces[iT].size();
    }
 
    inline const int serDegreeCell(size_t iT) const
    {
      return m_ddrcore.degree() + 1 - n_serendipityFaces(iT);
    }
 
    //---------- Serendipity dimensions and bases for grad -------------//
    inline size_t dimFacePolyl(size_t iF) const
    {
      return PolynomialSpaceDimension<Face>::Poly(serDegreeFace(iF));
    }
 
    inline size_t dimFacePolyl(const Face & F) const
    {
      return dimFacePolyl(F.global_index());
    }
 
    inline const PolylBasisFaceType & faceBasisPolyl(size_t iF) const
    {
      // Make sure that the basis has been created
      assert( m_face_bases_Polyl[iF] );
      return *m_face_bases_Polyl[iF].get();
 
    }
 
    inline const PolylBasisFaceType & faceBasisPolyl(const Face & F) const
    {
      return faceBasisPolyl(F.global_index()); 
    }
    
    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::VectorXi nDOFs_faces_SXGrad() const;
    Eigen::VectorXi nDOFs_cells_SXGrad() const;
 
    //---------- Serendipity dimensions and bases for curl -------------//
    inline size_t dimFaceRolyCompllpo(size_t iF) const
    {
      return PolynomialSpaceDimension<Face>::RolyCompl(serDegreeFace(iF)+1);
    }
 
    inline size_t dimFaceRolyCompllpo(const Face & F) const
    {
      return dimFaceRolyCompllpo(F.global_index());
    }
 
    inline const RolyCompllpoBasisFaceType & faceBasisRolyCompllpo(size_t iF) const
    {
      // Make sure that the basis has been created
      assert( m_face_bases_RolyCompllpo[iF] );
      return *m_face_bases_RolyCompllpo[iF].get();
 
    }
 
    inline const RolyCompllpoBasisFaceType & faceBasisRolyCompllpo(const Face & F) const
    {
      return faceBasisRolyCompllpo(F.global_index());
    }
    
    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::VectorXi nDOFs_faces_SXCurl() const;
    Eigen::VectorXi nDOFs_cells_SXCurl() const;
 
    //------------ Serendipity operators ---------------//
    const Eigen::MatrixXd SerendipityOperatorFace(const size_t iF, const Eigen::MatrixXd & LF) const;
 
    inline const Eigen::MatrixXd SerendipityOperatorFace(const Face & F, const Eigen::MatrixXd & LF) const
    {
      return SerendipityOperatorFace(F.global_index(), LF);
    };
 
    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 iF);
 
    InverseProblem _compute_inverse_problem_faces(size_t iF);
 
    std::vector<size_t> _compute_serendipity_faces(size_t iT);
 
    InverseProblem _compute_inverse_problem_cells(size_t iT);
 
    // Pointer do DDRCore
    const DDRCore & m_ddrcore;
    // Output stream
    std::ostream & m_output;
    
    // To determine if we want to actually use serendipity or not. False gives back the standard DDR
    bool m_use_serendipity;
 
    // Face and cell basis P^l
    std::vector<std::unique_ptr<PolylBasisFaceType> > m_face_bases_Polyl;
    std::vector<std::unique_ptr<PolylBasisCellType> > m_cell_bases_Polyl;
 
    // Face and cell basis R^{c,l+1}
    std::vector<std::unique_ptr<RolyCompllpoBasisFaceType> > m_face_bases_RolyCompllpo;
    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;  
    std::vector<std::vector<size_t>> m_serendipity_faces;
  
    // Store the inverse of the matrix to compute serendipity operators on faces/cells
    std::vector<InverseProblem> m_inverse_problem_faces;
    std::vector<InverseProblem> m_inverse_problem_cells;
    
  };
  
} // end of namespace HArDCore3D
 
#endif // SERENDIPITY_PROBLEM_HPP