hho-fullgradientdiff.hpp

Go to the documentation of this file.

// Implementation of the HHO scheme, with full gradient, for -div(K nabla u)=f with K scalar piecewise function.
//
// Author: Jerome Droniou (jerome.droniou@monash.edu)
//
 
/*
 *
 *  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 HHOFULLGRADIENTDIFF_HPP
#define HHOFULLGRADIENTDIFF_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 <hhospace.hpp>
#include <BoundaryConditions/BoundaryConditions.hpp>
#include <integralweight.hpp>
 
namespace HArDCore3D
{
 
  struct FullGradientDiffusion
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<double(const VectorRd &)> ForcingTermType;
    typedef std::function<double(const VectorRd &)> SolutionType;
    typedef std::function<VectorRd(const VectorRd &)> SolutionGradientType;
    typedef IntegralWeight PermeabilityType;
 
    FullGradientDiffusion(
                   const HHOSpace & hho_space,           
                   const BoundaryConditions & BC,     
                   bool use_threads,                  
                   std::ostream & output = std::cout  
                   );
 
    void assembleLinearSystem(
                              const ForcingTermType & f,      
                              const PermeabilityType & kappa,    
                              const SolutionType & u,          
                              Eigen::VectorXd & UDir          
                              );
 
    inline size_t numSCDOFs() const
    {
      return m_hhospace.mesh().n_cells() * m_nloc_sc;
    }
 
    inline size_t numDirDOFs() const
    {
      return m_BC.n_dir_faces()*m_hhospace.numLocalDofsFace();
    }
 
    inline size_t numSkeletalDOFs() const
    {
      return m_hhospace.dimension() - numSCDOFs();
    }
 
    inline size_t sizeSystem() const
    {
      return m_hhospace.dimension() - numSCDOFs() - numDirDOFs();
    }
 
    inline const HHOSpace & hhospace() const
    {
      return m_hhospace;
    }
    
    inline const Mesh & mesh() const
    {
      return m_hhospace.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;
    }
 
    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 ForcingTermType & f,      
                                const PermeabilityType & kappa     
                                );
 
    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 HHOSpace & m_hhospace;
    const BoundaryConditions & m_BC;
    bool m_use_threads;
    std::ostream & m_output;
    const size_t m_nloc_sc;   // Number of statically condensed 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;
 
  //------------------------------------------------------------------------------
  
  static const VectorRd vec_a = VectorRd(1.,2.,3.);
  
  static FullGradientDiffusion::SolutionType
  linear_u = [](const VectorRd & x) -> double {
                 return vec_a.dot(x) + 4.;
               };
 
  static FullGradientDiffusion::SolutionGradientType  
  linear_gradu = [](const VectorRd & x) -> VectorRd {
                     return vec_a;
                   };
 
  static FullGradientDiffusion::ForcingTermType
  linear_f = [](const VectorRd & x) -> double {
                 return 0.;
               };
 
  static FullGradientDiffusion::PermeabilityType
  linear_kappa = FullGradientDiffusion::PermeabilityType(1.);
  
 
  //------------------------------------------------------------------------------
 
  static FullGradientDiffusion::SolutionType
  trigonometric_u = [](const VectorRd & x) -> double {
                      return sin(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2));
                    };
 
  static FullGradientDiffusion::SolutionGradientType  
  trigonometric_gradu = [](const VectorRd & x) -> VectorRd {
                          return PI * VectorRd(
                                                     cos(PI*x(0)) * sin(PI*x(1)) * sin(PI*x(2)),
                                                     sin(PI*x(0)) * cos(PI*x(1)) * sin(PI*x(2)),
                                                     sin(PI*x(0)) * sin(PI*x(1)) * cos(PI*x(2))
                                                     );
                        };
 
  static FullGradientDiffusion::ForcingTermType
  trigonometric_f = [](const VectorRd & x) -> double {
                      return 3. * std::pow(PI, 2) * trigonometric_u(x);
                    };
 
  static FullGradientDiffusion::PermeabilityType
  trigonometric_kappa = FullGradientDiffusion::PermeabilityType(1.);
 
 
} // end of namespace HArDCore3D
 
#endif