hybridcore.hpp

Go to the documentation of this file.

// Data structures and methods to implement hybrid schemes in 2D, with polynomial unknowns
// in the cells and on the edges, such as Hybrid High-order (HHO) schemes.
//
//
// 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 HYBRIDCORE_HPP
#define HYBRIDCORE_HPP
 
#include <cassert>
#include <cmath>
 
#include <functional>
#include <memory>
#include <string>
#include <vector>
#include <mesh.hpp>
//#include <cell.hpp>
//#include <edge.hpp>
#include <quadraturerule.hpp>
#include <basis.hpp>
#include <iostream>
#include <parallel_for.hpp>
 
namespace HArDCore2D {
 
 
  // -----------------------------------------
  //          Free Functions
  // -----------------------------------------
  
  // Dimension of polynomial spaces -- only specialisations are relevant
  template<typename GeometricSupport>
  inline size_t const DimPoly(int m);
 
  template<>
  inline const size_t DimPoly<Cell>(const int m)   
    {
      return (m >= 0 ? (m + 1) * (m + 2) / 2 : 0);
    };
 
  template<>
  inline const size_t DimPoly<Edge>(const int m)   
    {
      return (m >= 0 ? m + 1 : 0);
    }
 
  //-----------------------------------------------------------------------------
  //                          UVector class definition
  //-----------------------------------------------------------------------------
  class UVector
  {
    public:
      UVector(
        const Eigen::VectorXd & values,   
        const Mesh & mesh,               
        const int cell_deg,             
        const size_t edge_deg           
      );
 
    inline Eigen::VectorXd & asVectorXd() const {
      return m_values;
    }
 
    inline const int get_cell_deg() const{
      return m_cell_deg;
    }
 
    inline const size_t get_edge_deg() const{
      return m_edge_deg;
    }
  
    inline const size_t n_cell_dofs() const{
      return DimPoly<Cell>( std::max(m_cell_deg,0) );
    }
 
    inline const size_t n_edge_dofs() const{
      return DimPoly<Edge>( m_edge_deg );
    }
 
    inline const size_t n_total_cell_dofs() const{
      return m_mesh.n_cells() * n_cell_dofs();
    }
 
    Eigen::VectorXd restr(size_t iT) const;
 
    UVector operator+(const UVector& b){
      assert(m_cell_deg == b.get_cell_deg() || m_edge_deg == b.get_edge_deg() );
      return UVector(m_values + b.asVectorXd(), m_mesh, m_cell_deg, m_edge_deg);
    }
 
    UVector operator-(const UVector& b){
      assert(m_cell_deg == b.get_cell_deg() || m_edge_deg == b.get_edge_deg() );
      return UVector(m_values - b.asVectorXd(), m_mesh, m_cell_deg, m_edge_deg);
    }
 
    void operator+=(const UVector& b){
      assert(m_cell_deg == b.get_cell_deg() || m_edge_deg == b.get_edge_deg() );
      this->asVectorXd() = m_values + b.asVectorXd();
    }
 
    void operator/=(const double & r){
      this->asVectorXd() = m_values / r;
    }
 
    double operator()(size_t index) const{
      return m_values(index);
    }
 
    private:
      mutable Eigen::VectorXd m_values;
      const Mesh& m_mesh;
      const int m_cell_deg;
      const size_t m_edge_deg;
 
  };
 
 
  // ----------------------------------------------------------------------------
  //                            HybridCore class definition
  // ----------------------------------------------------------------------------
 
  class HybridCore 
  {
 
  public:
 
    HybridCore(
           const Mesh* mesh_ptr, 
           const int cell_deg, 
           const size_t edge_deg, 
           const bool use_threads = true, 
           std::ostream & output = std::cout, 
           const bool ortho = true 
           ); 
 
    //---- Basis functions types -----//
    typedef Family<MonomialScalarBasisCell> PolyCellBasisType; 
    typedef Family<MonomialScalarBasisEdge> PolyEdgeBasisType; 
 
    //-------- Getters --------//
    inline const Mesh* get_mesh() const {return m_mesh;}
 
    inline const int CellDegree() const { return m_cell_deg; };
    inline const int CellDegreePos() const { return m_cell_deg_pos; };
    inline const size_t EdgeDegree() const { return m_edge_deg; };
 
    inline const PolyCellBasisType & CellBasis(size_t iT) const
    {
      // Make sure that the basis has been created
      assert( m_cell_basis[iT] );
      return *m_cell_basis[iT].get();
    }
 
    inline const PolyEdgeBasisType & EdgeBasis(size_t iE) const
    {
      // Make sure that the basis has been created
      assert( m_edge_basis[iE] );
      return *m_edge_basis[iE].get();
    }
 
    //-------- Functions ------------//
    double L2norm(
      const UVector &Xh     
    ) const;
 
    double H1norm(
      const UVector &Xh     
      ) const;
 
    template<typename ContinuousFunction>
    UVector interpolate(
                const ContinuousFunction& f,    
        const int deg_cell,    
        const size_t deg_edge,    
                size_t doe  
                ) const;  
    Eigen::VectorXd compute_weights(size_t iT) const;
    
    double evaluate_in_cell(const UVector Xh, size_t iT, VectorRd x) const; 
    double evaluate_in_edge(const UVector Xh, size_t iE, VectorRd x) const;
 
    Eigen::VectorXd VertexValues(
                 const UVector Xh,   
                 const std::string from_dofs 
                 );
 
  private:
    // Mesh
    const Mesh* m_mesh;  // Pointer to mesh data
    // Degree of the cell polynomials (can be -1).
    const int m_cell_deg;
    const size_t m_cell_deg_pos;
    // Degree of the edge polynomials
    const int m_edge_deg;
    // Using threads
    const bool m_use_threads;
    // Output stream
    std::ostream & m_output;
    // Orthonormalise or not the basis functions
    const bool m_ortho;
 
    // Cell and edges bases
    std::vector<std::unique_ptr<PolyCellBasisType>> m_cell_basis;
    std::vector<std::unique_ptr<PolyEdgeBasisType>> m_edge_basis;
 
    // Creates the cell and edge bases
    PolyCellBasisType _construct_cell_basis(size_t iT);
    PolyEdgeBasisType _construct_edge_basis(size_t iF);
 
  };
 
 
 
  // -------------------------------------------------------------
  // ------- Interpolates continuous function on discrete space
 
  template<typename ContinuousFunction>
  UVector HybridCore::interpolate(const ContinuousFunction& f, int cell_deg, size_t edge_deg, size_t doe) const {
        
    size_t cell_deg_pos = std::max(cell_deg, 0);
    assert(cell_deg <= CellDegree() || edge_deg <= EdgeDegree() || cell_deg > -1);
    
    // Sizes of local interpolation
    size_t n_cell_dofs = DimPoly<Cell>(cell_deg_pos);
    size_t n_total_cell_dofs = m_mesh->n_cells() * n_cell_dofs;
    size_t n_edge_dofs = DimPoly<Edge>(edge_deg);
    size_t n_total_edge_dofs = m_mesh->n_edges() * n_edge_dofs;
    Eigen::VectorXd XTF = Eigen::VectorXd::Zero( n_total_cell_dofs + n_total_edge_dofs );
 
    // Edge projections
    std::function<void(size_t, size_t)> construct_all_edge_projections
      = [&](size_t start, size_t end)->void
      {
        for (size_t iE = start; iE < end; iE++) {
          Edge E = *m_mesh->edge(iE);
 
          // Compute the L2 projection on edge
          QuadratureRule quadE = generate_quadrature_rule(E, doe);
          boost::multi_array<double, 2> phiE_quadE = evaluate_quad<Function>::compute(EdgeBasis(iE), quadE);
          Eigen::VectorXd UF = l2_projection<PolyEdgeBasisType>(f, EdgeBasis(iE), quadE, phiE_quadE);
 
          // Fill in the complete vector of DOFs
          XTF.segment(n_total_cell_dofs + iE * n_edge_dofs, n_edge_dofs) = UF;
        }
      };
    parallel_for(m_mesh->n_edges(), construct_all_edge_projections, m_use_threads);
 
    // Cell projections
    std::function<void(size_t, size_t)> construct_all_cell_projections
      = [&](size_t start, size_t end)->void
      {
        for (size_t iT = start; iT < end; iT++) {
          Cell* cell = m_mesh->cell(iT);
 
          // L2 projection in cell
          // Mass matrix in cell
          QuadratureRule quadT = generate_quadrature_rule(*cell, doe);
          boost::multi_array<double, 2> phiT_quadT = evaluate_quad<Function>::compute(CellBasis(iT), quadT);
          Eigen::MatrixXd MT = compute_gram_matrix(phiT_quadT, phiT_quadT, quadT, n_cell_dofs, n_cell_dofs, "sym");
 
          // Vector of integrals of f against basis functions
          Eigen::VectorXd bT = Eigen::VectorXd::Zero(n_cell_dofs);
          for (size_t i = 0; i < n_cell_dofs; i++) {
            for (size_t iqn = 0; iqn < quadT.size(); iqn++){
              bT(i) += quadT[iqn].w * phiT_quadT[i][iqn] * f(quadT[iqn].vector());
            }
          }
 
          // Vector of coefficients (on cell basis functions) of the L2(T) projection of f
          Eigen::VectorXd UT = MT.ldlt().solve(bT);
 
          // Fill in the complete vector of DOFs
          size_t offset_T = iT * n_cell_dofs;
          XTF.segment(offset_T, n_cell_dofs) = UT;
 
          // Special case of L=-1, we replace the previously computed cell value with the proper average of edge values
          if ( cell_deg==-1 ){
            size_t nedges = cell->n_edges();
              // Weights corresponding to the values in cells/on edges
              Eigen::VectorXd barycoefT = compute_weights(iT);
              // We transform these weights so that they apply to the coefficients on the basis functions
              VectorRd xT = cell->center_mass();
              double phiT_cst = CellBasis(iT).function(0, xT);
              for (size_t ilE = 0; ilE < nedges; ilE++){
                VectorRd xE = cell->edge(ilE)->center_mass();
                size_t iE = cell->edge(ilE)->global_index();
                double phiE_cst = EdgeBasis(iE).function(0, xE);
                barycoefT(ilE) *= phiE_cst / phiT_cst;
              }
 
              XTF(iT) = 0;
              for (size_t ilE = 0; ilE < nedges; ilE++) {
                size_t iE = cell->edge(ilE)->global_index();
                XTF(iT) += barycoefT(ilE) * XTF(n_total_cell_dofs + iE);
              }
          }
 
        }
      };
    parallel_for(m_mesh->n_cells(), construct_all_cell_projections, m_use_threads);
 
    return UVector(XTF, *get_mesh(), cell_deg, edge_deg);
  }
 
 
 
}  // end of namespace HArDCore2D
 
#endif /* HYBRIDCORE_HPP */