hho-stokes.hpp

Go to the documentation of this file.

// Implementation of the HHO scheme, with full gradient, for the Stokes equations -div(mu grad u)+grad p=f, div u = g.
//
// Author: Jerome Droniou (jerome.droniou@umontpellier.fr), Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
//
 
/*
 *
 *  This implementation of HHO was developped following the principles described in 
 * Appendix B of the book
 *
 * 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.
 *
 * If you use this code or part of it for a scientific publication, please cite the book
 * above as a reference for the implementation.
 *
 */
 
 
#ifndef HHOSTOKES_HPP
#define HHOSTOKES_HPP
 
#include <iostream>
 
#include <boost/math/constants/constants.hpp>
 
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
 
#include <mesh.hpp>
// #include <mesh_builder.hpp>
#include <local_static_condensation.hpp>
#include <linearsolver.hpp>
 
#include <vhhospace.hpp>
#include <BoundaryConditions/BoundaryConditions.hpp>
#include <integralweight.hpp>
 
 
namespace HArDCore2D
{
  
  //------------------------------------------------------------------------------
  // Stokes model
  //------------------------------------------------------------------------------
 
 
  struct Stokes
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<VectorRd(const VectorRd &)> MomentumForcingTermType;
    typedef std::function<double(const VectorRd &)> CompressibilityForcingTermType;
    typedef std::function<VectorRd(const VectorRd &)> VelocityType;
    typedef std::function<MatrixRd(const VectorRd &)> VelocityGradientType;
    typedef std::function<double(const VectorRd &)> PressureType;
    typedef std::function<VectorRd(const VectorRd &)> PressureGradientType;
    typedef IntegralWeight ViscosityType;
 
    Stokes(
           const VHHOSpace & vhho_space,      
           const GlobalDOFSpace & p_space,    
           const BoundaryConditions & BC,     
           bool use_threads,                  
           std::ostream & output = std::cout  
           );
 
    void assembleLinearSystem(
                              const MomentumForcingTermType & f,      
                              const CompressibilityForcingTermType & g,      
                              const ViscosityType & mu,       
                              const VelocityType & u,         
                              Eigen::VectorXd & UDir          
                              );
                    
    inline size_t nloc_sc_u() const
    {
      return m_nloc_sc_u;
    }
 
    inline size_t nloc_sc_p() const
    {
      return m_nloc_sc_p;
    }
 
    inline size_t numSCDOFs_u() const
    {
      return m_vhhospace.mesh().n_cells() * m_nloc_sc_u;
    }
 
    inline size_t numSCDOFs_p() const
    {
      return m_vhhospace.mesh().n_cells() * m_nloc_sc_p;
    }
 
    inline size_t numSCDOFs() const
    {
      return numSCDOFs_u() + numSCDOFs_p();
    }
 
    inline size_t numDirDOFs() const
    {
      return m_BC.n_dir_edges()*m_vhhospace.numLocalDofsEdge();
    }
    
    inline size_t dimVelocity() const
    {
      return m_vhhospace.dimension();
    }
 
    inline size_t dimPressure() const
    {
      return m_pspace.dimension();
    }
    
    inline size_t numNonSCDOFs() const
    {
      return dimVelocity() + dimPressure() - numSCDOFs() + 1;
    }
 
    inline size_t sizeSystem() const
    {
      return numNonSCDOFs() - numDirDOFs();
    }
 
    inline const VHHOSpace & vhhospace() const
    {
      return m_vhhospace;
    }
    
    inline const GlobalDOFSpace & pspace() const
    {
      return m_pspace;
    }
 
    inline const Mesh & mesh() const
    {
      return m_vhhospace.mesh();
    }
 
    inline const SystemMatrixType & systemMatrix() const {
      return m_A;
    }
 
    inline SystemMatrixType & systemMatrix() {
      return m_A;
    }
 
    inline const Eigen::VectorXd & systemVector() const {
      return m_b;
    }
 
    inline Eigen::VectorXd & systemVector() {
      return m_b;
    }
 
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    inline const SystemMatrixType & scMatrix() const {
      return m_sc_A;
    }
 
    inline Eigen::VectorXd & scVector() {
      return m_sc_b;
    }
 
    Eigen::VectorXd interpolate(
                                const VelocityType & u, 
                                const PressureType & p, 
                                const int doe_cell = -1, 
                                const int doe_face = -1 
                                ) const;    
 
    std::vector<double> computeEnergyNorms(
                                           const std::vector<Eigen::VectorXd> & list_dofs   
                                           ) const;
 
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                            
                                const MomentumForcingTermType & f,      
                                const CompressibilityForcingTermType & g,    
                                const ViscosityType & mu              
                                );
 
    LocalStaticCondensation _compute_static_condensation(const size_t & iT) const;
 
    void _assemble_local_contribution(
                                      size_t iT,                                               
                                      const std::pair<Eigen::MatrixXd, Eigen::VectorXd> & lsT, 
                                      std::list<Eigen::Triplet<double> > & A1,                 
                                      Eigen::VectorXd & b1,                                    
                                      std::list<Eigen::Triplet<double> > & A2,                 
                                      Eigen::VectorXd & b2                                     
                                      );
    
    const VHHOSpace & m_vhhospace;
    const GlobalDOFSpace & m_pspace; 
 
    const BoundaryConditions & m_BC;
    bool m_use_threads;
    std::ostream & m_output;
    const size_t m_nloc_sc_u; // Number of statically condensed velocity DOFs in each cell
    const size_t m_nloc_sc_p; // Number of statically condensed pressure DOFs in each cell
    SystemMatrixType m_A;     // Matrix and RHS of statically condensed system
    Eigen::VectorXd m_b;
    SystemMatrixType m_sc_A;  // Static condensation operator and RHS (to recover statically condensed DOFs)
    Eigen::VectorXd m_sc_b;
    double m_stab_par;
  };
 
  
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
 
  //------------------------------------------------------------------------------
  // Linear
 
  double pressure_scaling = 1.;
  
  static Stokes::VelocityType
  linear_u = [](const VectorRd & x) -> VectorRd {
    return VectorRd(x(0), -x(1));
  };
 
  static Stokes::VelocityGradientType
  linear_gradu = [](const VectorRd & x) -> MatrixRd {
    MatrixRd M = MatrixRd::Zero();
    M << 1., 0., 0., -1.;      
    return M;
  };
 
  static Stokes::PressureType
  linear_p = [](const VectorRd & x) -> double {
    return 0.;
  };
 
  static Stokes::PressureGradientType
  linear_gradp = [](const VectorRd & x) -> VectorRd {
    return VectorRd::Zero();
  };
 
  static Stokes::MomentumForcingTermType
  linear_f = [](const VectorRd & x) -> VectorRd {
    return VectorRd::Zero();
  };
 
  static Stokes::CompressibilityForcingTermType 
  linear_g = [](const VectorRd & x)->double { return linear_gradu(x).trace();};
 
  static Stokes::ViscosityType
  linear_mu = Stokes::ViscosityType(1.);
 
  //------------------------------------------------------------------------------
  // Trigonometric
 
  static Stokes::VelocityType
  trigonometric_u = [](const VectorRd & x) -> VectorRd {
    return exp(x(0)) * VectorRd(
                    -(x(1) * cos(x(1)) + sin(x(1))),
                    x(1) * sin(x(1))
                    );
  };
 
  static Stokes::VelocityGradientType
  trigonometric_gradu = [](const VectorRd & x) -> MatrixRd {
    MatrixRd M = MatrixRd::Zero(); 
    M <<
      -x(1) * cos(x(1)) - sin(x(1)),  -2. * cos(x(1)) + x(1) * sin(x(1)),
      x(1) * sin(x(1)), sin(x(1)) + x(1) * cos(x(1));
    M *= exp(x(0));
    return M;
  };
 
  static Stokes::PressureType
  trigonometric_p = [](const VectorRd & x) -> double {
    return 2. * exp(x(0)) * sin(x(1)) - 2. * (exp(1) - 1.) * (1. - cos(1.));
  };
 
  static Stokes::PressureGradientType
  trigonometric_gradp = [](const VectorRd & x) -> VectorRd {
    return 2. * exp(x(0)) * VectorRd(sin(x(1)), cos(x(1)));
  };
 
  static Stokes::MomentumForcingTermType
  trigonometric_f = [](const VectorRd & x) -> VectorRd {
    return -2. * exp(x(0)) * VectorRd(sin(x(1)), cos(x(1)))
      + trigonometric_gradp(x);
  };
 
  static Stokes::CompressibilityForcingTermType 
  trigonometric_g = [](const VectorRd & x)->double { return trigonometric_gradu(x).trace();};
 
  static Stokes::ViscosityType
  trigonometric_mu = Stokes::ViscosityType(1.);
 
 
} // end of namespace HArDCore3D
 
#endif