ddr-stokes.hpp

Go to the documentation of this file.

// Implementation of the discrete de Rham sequence for the Stokes problem in curl-curl formulation.
//
// Authors: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
//          Jerome Droniou (jerome.droniou@monash.edu)
//
 
 
#ifndef STOKES_HPP
#define STOKES_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 <xgrad.hpp>
#include <xcurl.hpp>
#include <xdiv.hpp>
 
namespace HArDCore3D
{
 
  struct StokesNorms
  {
    StokesNorms(double norm_u, double norm_curl_u, double norm_p, double norm_grad_p):
      u(norm_u),
      curl_u(norm_curl_u),
      p(norm_p),
      grad_p(norm_grad_p),
      hcurl_u( std::sqrt( std::pow(norm_u, 2) + std::pow(norm_curl_u, 2) ) ),
      hgrad_p( std::sqrt( std::pow(norm_p, 2) + std::pow(norm_grad_p, 2) ) )
      {
        // do nothing
      };
    
    // Constructor without arguments
    StokesNorms() : StokesNorms(0., 0., 0., 0.) {};
    
    double u; 
    double curl_u; 
    double p; 
    double grad_p; 
    double hcurl_u; 
    double hgrad_p; 
  
  };
 
  struct Stokes
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> ForcingTermType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> VelocityType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> VorticityType;
    typedef std::function<double(const Eigen::Vector3d &)> PressureType;
    typedef std::function<Eigen::Vector3d(const Eigen::Vector3d &)> PressureGradientType;
    typedef IntegralWeight ViscosityType;
 
    Stokes(
           const DDRCore & ddrcore,           
           bool use_threads,                  
           std::ostream & output = std::cout  
         );
 
    void assembleLinearSystem(
                              const ForcingTermType & f,      
                              const VelocityType & u, 
                              const VorticityType & omega, 
                              const ViscosityType & nu    
                              );
 
    inline size_t dimension() const
    {
      return m_xcurl.dimension() + m_xgrad.dimension();
    }
 
    inline size_t nbSCDOFs_u() const
    {
      return m_ddrcore.mesh().n_cells() * m_nloc_sc_u;
    }
 
    inline size_t nbSCDOFs_p() const
    {
      return m_ddrcore.mesh().n_cells() * m_nloc_sc_p;
    }
 
    inline size_t nbSCDOFs() const
    {
      return m_ddrcore.mesh().n_cells() * (m_nloc_sc_u + m_nloc_sc_p);
    }
 
    inline size_t sizeSystem() const
    {
      return dimension() + 1 - nbSCDOFs();
    }
 
    inline const XGrad & xGrad() const
    {
      return m_xgrad;
    }
 
    inline const XCurl & xCurl() const
    {
      return m_xcurl;
    }
 
    inline const XDiv & xDiv() const
    {
      return m_xdiv;
    }
 
    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 SystemMatrixType & scMatrix() const {
      return m_sc_A;
    }
 
    inline Eigen::VectorXd & scVector() {
      return m_sc_b;
    }
 
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    std::vector<StokesNorms> computeStokesNorms(
                const std::vector<Eigen::VectorXd> & list_dofs 
                ) const;
 
    std::pair<StokesNorms, StokesNorms> computeContinuousErrorsNorms(
                const Eigen::VectorXd & v, 
                const VelocityType & u, 
                const VorticityType & curl_u, 
                const PressureType & p, 
                const PressureGradientType & grad_p 
                ) const;
 
    double evaluateDefectCommutationGradInterp(
                       const PressureType & p, 
                       const PressureGradientType & grad_p,  
                       const size_t deg_quad_interpolate    
                      ) const;                     
 
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                      
                                const Eigen::VectorXd & interp_f,  
                                const ViscosityType & nu    
                                );
 
    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 DDRCore & m_ddrcore;
    bool m_use_threads;
    std::ostream & m_output;
    const XGrad m_xgrad;
    const XCurl m_xcurl;
    const XDiv m_xdiv;
    const size_t m_nloc_sc_u; // Nb of statically condensed DOFs for velocity in each cell (cell unknowns)
    const size_t m_nloc_sc_p; // Nb of statically condensed DOFs for pressure in each cell (cell unknowns)
    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;
  using std::pow;
 
  //------------------------------------------------------------------------------
  // Trigonometric solution
  //------------------------------------------------------------------------------
 
  double pressure_scaling = 1.;
  static Stokes::VelocityType
  trigonometric_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return Eigen::Vector3d(
                                             0.5 * sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                             0.5 * cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                             -cos(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2))
                                             );
                    };
 
  static Stokes::VorticityType
  trigonometric_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return 3. * PI * Eigen::Vector3d(
                                                      cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      -sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      0.
                                                      );
                         };
 
  static Stokes::PressureType
  trigonometric_p = [](const Eigen::Vector3d & x) -> double {
                      return pressure_scaling * sin(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2));
                    };
 
  static Stokes::PressureGradientType
  trigonometric_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return pressure_scaling * 2. * PI * Eigen::Vector3d(
                                                      cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2)),
                                                      sin(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2))
                                                      );
                         };
 
  static Stokes::ForcingTermType
  trigonometric_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                      return 6. * std::pow(PI, 2) * Eigen::Vector3d(
                                                                    sin(2. * PI * x(0)) * cos(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                                                    cos(2. * PI * x(0)) * sin(2. * PI * x(1)) * cos(2. * PI * x(2)),
                                                                    -2. * cos(2. * PI * x(0)) * cos(2. * PI * x(1)) * sin(2. * PI * x(2))
                                                                    )
                        + trigonometric_grad_p(x);
                    };
 
  static Stokes::ViscosityType
  trigonometric_nu = Stokes::ViscosityType(1.);
 
  //------------------------------------------------------------------------------
  // Linear solution
  //------------------------------------------------------------------------------
 
  static Stokes::VelocityType
  linear_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return Eigen::Vector3d(x(0), -x(1), 0.);
             };
 
  static Stokes::VorticityType
  linear_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d::Zero();
                  };
 
  static Stokes::PressureType
  linear_p = [](const Eigen::Vector3d & x) -> double {
               return pressure_scaling * (x(0) - x(1));
             };
 
  static Stokes::PressureGradientType
  linear_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return pressure_scaling * Eigen::Vector3d(1., -1., 0.);
                  };
 
  static Stokes::ForcingTermType
  linear_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               return linear_grad_p(x);
             };
 
  static Stokes::ViscosityType
  linear_nu = Stokes::ViscosityType(1.);
 
 
  //------------------------------------------------------------------------------
  // Solution derived from a field
  //------------------------------------------------------------------------------
  // Setting f(x,y,z)= (sin(pi x) sin(pi y) sin(pi z))^3, we take 
  //    U=(-d_y f, d_x f, 0), p=0
  // Calculations are made symbolically in matlab
 
  static Stokes::VelocityType
  field_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
               double ux = PI*cos(PI*x(1))*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(1)),2.0)*pow(sin(PI*x(2)),3.0)*3.0;
               double uy = PI*cos(PI*x(0))*pow(sin(PI*x(0)),2.0)*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),3.0)*(-3.0);
               return Eigen::Vector3d(ux, uy, 0.);
             };
 
  static Stokes::VorticityType
  field_curl_u = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                double cux = (PI*PI)*cos(PI*x(0))*cos(PI*x(2))*pow(sin(PI*x(0)),2.0)*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),2.0)*9.0;
                double cuy = (PI*PI)*cos(PI*x(1))*cos(PI*x(2))*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(1)),2.0)*pow(sin(PI*x(2)),2.0)*9.0;
                double cuz = (PI*PI)*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),3.0)*6.0-(PI*PI)*pow(cos(PI*x(0)),2.0)*sin(PI*x(0))*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),3.0)*6.0-(PI*PI)*pow(cos(PI*x(1)),2.0)*pow(sin(PI*x(0)),3.0)*sin(PI*x(1))*pow(sin(PI*x(2)),3.0)*6.0;
 
                return Eigen::Vector3d(cux, cuy, cuz);
                };
 
  static Stokes::PressureType
  field_p = [](const Eigen::Vector3d & x) -> double {
               return 0.;
             };
 
  static Stokes::PressureGradientType
  field_grad_p = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
                    return Eigen::Vector3d(0., 0., 0.);
                  };
 
  static Stokes::ForcingTermType
  field_f = [](const Eigen::Vector3d & x) -> Eigen::Vector3d {
              double fx = (PI*PI*PI)*pow(cos(PI*x(1)),3.0)*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(2)),3.0)*-6.0+(PI*PI*PI)*cos(PI*x(1))*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(1)),2.0)*pow(sin(PI*x(2)),3.0)*3.9E+1-(PI*PI*PI)*pow(cos(PI*x(0)),2.0)*cos(PI*x(1))*sin(PI*x(0))*pow(sin(PI*x(1)),2.0)*pow(sin(PI*x(2)),3.0)*1.8E+1-(PI*PI*PI)*cos(PI*x(1))*pow(cos(PI*x(2)),2.0)*pow(sin(PI*x(0)),3.0)*pow(sin(PI*x(1)),2.0)*sin(PI*x(2))*1.8E+1;
 
              double fy = (PI*PI*PI)*pow(cos(PI*x(0)),3.0)*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),3.0)*6.0-(PI*PI*PI)*cos(PI*x(0))*pow(sin(PI*x(0)),2.0)*pow(sin(PI*x(1)),3.0)*pow(sin(PI*x(2)),3.0)*3.9E+1+(PI*PI*PI)*cos(PI*x(0))*pow(cos(PI*x(1)),2.0)*pow(sin(PI*x(0)),2.0)*sin(PI*x(1))*pow(sin(PI*x(2)),3.0)*1.8E+1+(PI*PI*PI)*cos(PI*x(0))*pow(cos(PI*x(2)),2.0)*pow(sin(PI*x(0)),2.0)*pow(sin(PI*x(1)),3.0)*sin(PI*x(2))*1.8E+1;
 
              double fz = 0.;
               return Eigen::Vector3d(fx, fy, fz);
             };
 
  static Stokes::ViscosityType
  field_nu = Stokes::ViscosityType(1.);
 
  //------------------------------------------------------------------------------
  
 
} // end of namespace HArDCore3D
 
#endif