ddr-quadrot.hpp

Go to the documentation of this file.

// 
// Solver for a quad-rot problem
// Author: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr)
//
 
 
#ifndef QUADROT_HPP
#define QUADROT_HPP
 
#include <iostream>
 
#include <boost/math/constants/constants.hpp>
 
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
 
#include <mesh.hpp>
#include <BoundaryConditions/BoundaryConditions.hpp>   // To re-order the boundary edges and vertices
#include <BoundaryConditions/BChandlers.hpp>   
 
#include <xgrad.hpp>
#include <xrotrot.hpp>
#include <xrot.hpp>
 
namespace HArDCore2D
{
 
  struct QuadRot {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<VectorRd(const VectorRd &)> PotentialType;
    typedef std::function<double(const VectorRd &)> RotorType;
    typedef std::function<VectorRd(const VectorRd &)> RotRotType;
    typedef std::function<double(const VectorRd &)> LagrangeMultiplierType;
    typedef std::function<VectorRd(const VectorRd &)> ForcingTermType;
    
    QuadRot(
        const DDRCore & ddrcore,          
            const BoundaryConditions & bc_u,  
            const BoundaryConditions & bc_p,  
        bool use_threads = true,          
        std::ostream & output = std::cout 
        );
 
    inline const XGrad & xGrad() const
    {
      return m_xgrad;
    }
 
    inline const XRotRot & xRotRot() const
    {
      return m_xrotrot;
    }
 
    inline const XRot & xRot() const
    {
      return m_xrot;
    }
    
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    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;
    }
 
    std::vector<size_t> globalDOFIndices(const Cell & T) const;
 
    inline size_t dimensionLocalSpace(size_t iT) const
    {
      return m_xrotrot.dimensionCell(iT) + m_xgrad.dimensionCell(iT);
    }
 
    inline size_t dimensionSpace() const
    {
      return m_xrotrot.dimension() + m_xgrad.dimension();
    }
 
    inline size_t numBoundaryDofs() const
    {
      return m_bc_u.n_dir_vertices() * m_xrotrot.numLocalDofsVertex()
        + m_bc_u.n_dir_edges() * m_xrotrot.numLocalDofsEdge()
        + m_bc_p.n_dir_vertices() * m_xgrad.numLocalDofsVertex()
        + m_bc_p.n_dir_edges() * m_xgrad.numLocalDofsEdge();
    }
    
    inline size_t dimensionLinearSystem() const
    {
      return dimensionSpace() - numBoundaryDofs();
    }
 
    inline int dofToUnknown(size_t i) const
    {
      return m_dof_to_unknown(i);
    }
 
    Eigen::VectorXd assembleLinearSystem(
                                         const ForcingTermType & f,       
                                         const PotentialType & u,         
                                         const RotorType & rot_u,         
                                         const LagrangeMultiplierType & p 
                                         );
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                size_t iT,                
                const ForcingTermType & f 
                );
 
    void _assemble_local_contribution(
                                      size_t iT,                                               
                                      const std::pair<Eigen::MatrixXd, Eigen::VectorXd> & lsT, 
                                      std::list<Eigen::Triplet<double> > & triplets_A,         
                                      Eigen::VectorXd & rsh_b,                                 
                                      std::list<Eigen::Triplet<double> > & triplets_boundary_A 
                                      );
    
    const DDRCore & m_ddrcore;
    BoundaryConditions m_bc_u;
    BoundaryConditions m_bc_p;
    bool m_use_threads;
    std::ostream & m_output;
    double m_stab_par;
    XGrad m_xgrad;
    XRotRot m_xrotrot;
    XRot m_xrot;
    SystemMatrixType m_A;
    SystemMatrixType m_B;
    Eigen::VectorXd m_b;
    // Indicates the location, among the DOFs, of the unknowns after the Dirichlet DOFs are removed
    // (m_unknowns_stride[i].first = starting index, m_unknowns_stride[i].second = number of unknowns)
    std::vector<std::pair<size_t, size_t>> m_unknowns_stride;
    Eigen::VectorXi m_dof_to_unknown;
  };
 
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
 
  //------------------------------------------------------------------------------
  // Linear
 
  static QuadRot::PotentialType
  linear_u = [](const VectorRd & x) -> VectorRd {
    VectorRd u_x;
    u_x << -(x(1) - 0.5), (x(0) - 0.5);
    return u_x;
  };
 
  static QuadRot::RotorType
  linear_rot_u = [](const VectorRd & x) -> double {
    return 2.;
  };
 
  static QuadRot::RotRotType
  linear_rotrot_u = [](const VectorRd & x)->VectorRd {
    return VectorRd::Zero();
  };
  
  static QuadRot::LagrangeMultiplierType
  linear_p = [](const VectorRd & x) -> double {
    return 0.;
  };
 
  static QuadRot::ForcingTermType
  linear_f = [](const VectorRd & x) -> VectorRd {
    VectorRd f_x;
    f_x << 0, 0;
    return f_x;
  };
 
  //------------------------------------------------------------------------------
  // Quadratic
 
  static QuadRot::PotentialType
  quadratic_u = [](const VectorRd & x) -> VectorRd {
    VectorRd u_x;
    u_x << -pow(x(1) - 0.5, 2), pow(x(0) - 0.5, 2);
    return u_x;
  };
 
  static QuadRot::RotorType
  quadratic_rot_u = [](const VectorRd & x) -> double {
    return 2. * ( x(0) + x(1) - 1 );
  };
 
  static QuadRot::RotRotType
  quadratic_rotrot_u = [](const VectorRd & x)->VectorRd {
    VectorRd rotrot_u_x;
    rotrot_u_x << 2., -2.;
    return rotrot_u_x;
  };
  
  static QuadRot::LagrangeMultiplierType
  quadratic_p = [](const VectorRd & x) -> double {
    return 0.;
  };
 
  static QuadRot::ForcingTermType
  quadratic_f = [](const VectorRd & x) -> VectorRd {
    VectorRd f_x;
    f_x << 0., 0.;
    return f_x;
  };
 
  //------------------------------------------------------------------------------
  // Cubic
  
  static QuadRot::PotentialType
  cubic_u = [](const VectorRd & x) -> VectorRd {
    VectorRd u_x;
    u_x << -pow(x(1) - 0.5, 3), pow(x(0) - 0.5, 3);
    return u_x;
  };
 
  static QuadRot::RotorType
  cubic_rot_u = [](const VectorRd & x) -> double {
    return 3. * ( pow(x(0) - 0.5, 2) + pow(x(1) - 0.5, 2) );
  };
 
  static QuadRot::RotRotType
  cubic_rotrot_u = [](const VectorRd & x)->VectorRd {
    VectorRd rotrot_u_x;
    rotrot_u_x << 6. * (x(1) - 0.5), -6. * (x(0) - 0.5);
    return rotrot_u_x;
  };
  
  static QuadRot::LagrangeMultiplierType
  cubic_p = [](const VectorRd & x) -> double {
    return 0.;
  };
 
  static QuadRot::ForcingTermType
  cubic_f = [](const VectorRd & x) -> VectorRd {
    VectorRd f_x;
    f_x << 0., 0.;
    return f_x;
  };
  
  //------------------------------------------------------------------------------
  // Quartic
 
  static QuadRot::PotentialType
  quartic_u = [](const VectorRd & x) -> VectorRd {
    VectorRd u_x;
    u_x << -pow(x(1) - 0.5, 4), pow(x(0) - 0.5, 4);
    return u_x;
  };
 
  static QuadRot::RotorType
  quartic_rot_u = [](const VectorRd & x) -> double {
    return 4. * ( pow(x(0)- 0.5, 3) + pow(x(1) - 0.5, 3) );
  };
 
  static QuadRot::RotRotType
  quartic_rotrot_u = [](const VectorRd & x)->VectorRd {
    VectorRd rotrot_u_x;
    rotrot_u_x << 12. * pow(x(1) - 0.5, 2),
      -12. * pow(x(0) - 0.5 , 2);
    return rotrot_u_x;
  };
 
  // rot rot rot (u) = -24. * ( x(0) + x(1) - 1 )
  
  static QuadRot::LagrangeMultiplierType
  quartic_p = [](const VectorRd & x) -> double {
    return x(0) * (1. - x(0)) * x(1) * (1. - x(1));
  };
 
  static QuadRot::ForcingTermType
  quartic_f = [](const VectorRd & x) -> VectorRd {
    VectorRd f_x;
    f_x << -24. + (1. - 2. * x(0)) * x(1) * (1. - x(1)),
      24. + x(0) * (1. - x(0)) * (1. - 2. * x(1));
    return f_x;
  }; 
  
  //------------------------------------------------------------------------------
  // Trigonometric
 
  static QuadRot::PotentialType
  trigonometric_u = [](const VectorRd & x) -> VectorRd {
    VectorRd u_x;
    u_x << -sin(PI * (x(0) + x(1))), sin(PI * (x(0) + x(1)));
    return u_x;
  };
 
  static QuadRot::RotorType
  trigonometric_rot_u = [](const VectorRd & x) -> double {
    return 2. * PI * cos(PI * (x(0) + x(1)));
  };
 
  static QuadRot::RotRotType
  trigonometric_rotrot_u = [](const VectorRd & x) -> VectorRd {
    VectorRd rotrot_u_x;
    rotrot_u_x << -sin(PI * (x(0) + x(1))), sin(PI * (x(0) + x(1)));
    return pow(PI, 2) * rotrot_u_x;
  };
 
  static QuadRot::RotorType
  trigonometric_rotrotrot_u = [](const VectorRd & x) -> double {
    return 4. * pow(PI, 3) * cos(PI * (x(0) + x(1)));
  };
  
  static QuadRot::LagrangeMultiplierType
  trigonometric_p = [](const VectorRd & x) -> double {
    return sin(PI * x(0)) * sin(PI * x(1));
  };
 
  static QuadRot::ForcingTermType
  trigonometric_f = [](const VectorRd & x) -> VectorRd {
    VectorRd f_x;
    f_x << -4. * pow(PI, 4) * sin(PI * (x(0) + x(1))) + PI * cos(PI * x(0)) * sin(PI * x(1)),
      4. * pow(PI, 4) * sin(PI * (x(0) + x(1))) + PI * sin(PI * x(0)) * cos(PI * x(1));
    return f_x;
  };  
 
} // end of namespace HArDCore2D
 
#endif