vhhospace.hpp

Go to the documentation of this file.

// Core data structures and methods required to implement the Hybrid High-Order in 2D, for vector-valued functions
//
// Provides:
//  - Polynomial spaces on the element and edges
//  - Interpolator of smooth functions
//  - Full gradient, potential and stabilisation bilinear form in the elements
//
// Author: Alessandra Crippa (alessandra.crippa@ifpen.fr)
//
 
/*
 *
 *      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 VHHOSPACE_HPP
#define VHHOSPACE_HPP
 
#include <iostream>
#include <globaldofspace.hpp>
#include <basis.hpp>
#include <polynomialspacedimension.hpp>
 
 
namespace HArDCore2D
{
 
  //------------------------------------------------------------------------------
 
  // Type for selecting some cells (return 0 or 1 for each cell T)
  typedef std::function<bool(const Cell &)> CellSelection;
  static const CellSelection allcells = [](const Cell &)->bool {return true;};
 
  class VHHOSpace : public GlobalDOFSpace
  {
  public:
    // Types for element bases
    typedef Family<MonomialScalarBasisCell> PolyBasisCellType;
    typedef TensorizedVectorFamily<PolyBasisCellType, dimspace> PolydBasisCellType;
    typedef MatrixFamily<PolyBasisCellType, dimspace> PolydxdBasisCellType;
    // added AC
    typedef Family<PolydxdBasisCellType> PolySymdxdBasisCellType;
 
    // Types for edge bases
    typedef Family<MonomialScalarBasisEdge> PolyBasisEdgeType;
    typedef Family<TensorizedVectorFamily<PolyBasisEdgeType, dimspace>> PolydBasisEdgeType;
 
    // Types for functions to interpolate
    typedef std::function<VectorRd(const VectorRd &)> FunctionType;
    // AC - type needed only for tests
    typedef std::function<MatrixRd(const VectorRd &)> GradientType;
    
 
    struct CellBases
    {
      typedef Cell GeometricSupport;
 
      std::unique_ptr<PolyBasisCellType> Polykpo;
      std::unique_ptr<PolyBasisCellType> Polyk;
      std::unique_ptr<PolydBasisCellType> Polykpod;
      std::unique_ptr<PolydBasisCellType> Polykd;
      std::unique_ptr<PolydxdBasisCellType> Polykdxd;
      // added AC
      std::unique_ptr<PolySymdxdBasisCellType> PolySymkdxd;
    };
 
 
    struct EdgeBases
    {
      typedef Edge GeometricSupport;
 
      std::unique_ptr<PolyBasisEdgeType> Polyk;
      std::unique_ptr<PolydBasisEdgeType> Polykd;
    };
    
    struct LocalOperators
    {
      LocalOperators(
                     const Eigen::MatrixXd & _gradient, 
                       // Added AC - symmetric gradient
                     const Eigen::MatrixXd & _symmetric_gradient, 
                     const Eigen::MatrixXd & _potential, 
                     // added AC - potential from symmetrix gradient definition
                     const Eigen::MatrixXd & _potential_sym,
                     const Eigen::MatrixXd & _stabilisation, 
                     // added AC to remove diff operator
                     const Eigen::MatrixXd & _difference
                     )
        : gradient(_gradient),
          symmetric_gradient(_symmetric_gradient), // added AC
          potential(_potential),
          potential_sym(_potential_sym), // added AC
          stabilisation(_stabilisation),
          difference(_difference) //AC remove
      {
        // Do nothing
      }
      
      Eigen::MatrixXd gradient;
      Eigen::MatrixXd symmetric_gradient; //added AC
      Eigen::MatrixXd potential;
      Eigen::MatrixXd potential_sym; // added AC
      Eigen::MatrixXd stabilisation;
      Eigen::MatrixXd difference; // AC remove
    };
 
    VHHOSpace(const Mesh & mesh, size_t K, const CellSelection & BoundaryStab, bool use_threads = true, std::ostream & output = std::cout);
 
    VHHOSpace(const Mesh & mesh, size_t K, bool use_threads = true, std::ostream & output = std::cout)
            : VHHOSpace(mesh, K, allcells, use_threads, output) {};
    
    const Mesh & mesh() const
    {
      return m_mesh;
    }
 
    const size_t & degree() const
    {
      return m_K;
    }
 
    const CellSelection & boundaryStab() const
    {
      return m_boundary_stab;
    }
 
    Eigen::VectorXd interpolate(
          const FunctionType & q, 
          const int doe_cell = -1, 
          const int doe_face = -1 
          ) const;
 
    // AC
    Eigen::VectorXd local_interpolate(
            size_t iT, 
            const FunctionType & q, 
            const int doe_cell = -1, 
            const int doe_face = -1 
    ) const;
 
      // AC
      Eigen::VectorXd components(
              size_t iT,
              const FunctionType & q, 
              const int doe = -1 
      ) const;
 
    // AC
    Eigen::VectorXd components(
            size_t iT,
            const GradientType & q, 
            const int doe = -1 
    ) const;
 
    inline const CellBases & cellBases(size_t iT) const
    {
      // Make sure that the basis has been created
      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
    {
      // Make sure that the basis has been created
      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 & operators(size_t iT) const
    {
      assert( m_operators[iT] );
      return *m_operators[iT];
    }
 
    inline const LocalOperators & operators(const Cell & T) const
    {
      return operators(T.global_index());
    }
    
    std::vector<std::pair<double,double>> computeNorms(
                   const std::vector<Eigen::VectorXd> & list_dofs   
                  ) const;
 
    std::vector<Eigen::VectorXd> computeVertexValues(
                                                     const Eigen::VectorXd & u 
                                                     ) const;
   
  private:
    CellBases _construct_cell_bases(size_t iT);
 
    EdgeBases _construct_edge_bases(size_t iE);
    
    LocalOperators _compute_operators(size_t iT);
    
    // Pointer to the mesh
    const Mesh & m_mesh;
    // Degrees
    const size_t m_K;
    // Choice of boundary stabilisation
    CellSelection m_boundary_stab;
    // Parallel or not
    bool m_use_threads;
    // Output stream
    std::ostream & m_output;    
    
    // Cell bases
    std::vector<std::unique_ptr<CellBases> > m_cell_bases;
    // Edge bases
    std::vector<std::unique_ptr<EdgeBases> > m_edge_bases;
 
    // Local operators
    std::vector<std::unique_ptr<LocalOperators> > m_operators;
        
  };
 
  
} // end of namespace HArDCore2D
 
#endif // VHHOSPACE_HPP