bgg-klplate.hpp

Go to the documentation of this file.

// Implementation of the DDR-BGG scheme for the Kirchhoff-Love plate problem.
//
// Author: Arax Leroy (arax.leroy@umontpellier.fr)
//
 
 
#ifndef BGG_KLPLATE_HPP
#define BGG_KLPLATE_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 <vsxgrad.hpp>
#include <xhess.hpp>
 
namespace HArDCore2D
{
 
  struct KLNorms
  {
    KLNorms(double norm_displacement):
      displacement(norm_displacement),
      energy(norm_displacement)
      {
        // do nothing
      };
    double displacement; 
    double energy; 
  };
 
  struct KirchhoffLove
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    typedef std::function<double(const Eigen::Vector2d &)> ForcingTermType; 
    typedef std::function<double(const Eigen::Vector2d &)> SolutionDisplacementType;
    typedef std::function<Eigen::Vector2d(const Eigen::Vector2d&)> GradientDisplacementType;
 
 
    KirchhoffLove(
                   const DDRCore & ddrcore,               
                   const DDRCore & ddrcore_plus,          
                   const SerendipityProblem & sp,
                   const SerendipityProblem & sp_plus,    
                   const BoundaryConditions & BC_u,       
                   bool use_threads,                      
                   std::ostream & output = std::cout      
                   );
 
    void assembleLinearSystem(
                              const ForcingTermType & f,      
                              const SolutionDisplacementType & u,  
                              const GradientDisplacementType & grad_u, 
                              size_t iT_dirac                   
                              );
 
    inline size_t dimensionSpace() const
    {
      return m_xhess.dimension();
    }
 
    inline size_t nb_bdryDOFs() const
    { 
      return  m_BC_u.n_dir_edges() * m_xhess.numLocalDofsEdge()
             + m_BC_u.n_dir_vertices() * m_xhess.numLocalDofsVertex();
    }
 
    inline size_t sizeSystem() const
    {
      return dimensionSpace() - nb_bdryDOFs();
    }
 
    inline const std::vector<std::pair<size_t,size_t>> & locUKN() const {
      return m_locUKN;
    }
 
    std::vector<size_t> globalDOFIndices(const Cell &T) const
    {
      return m_xhess.globalDOFIndices(T);
    }
    
    inline const VSXGrad & vsxGrad() const
    {
      return m_vsxgrad;
    }
    
    inline const XHess & xHess() const
    {
      return m_xhess;
    }
 
    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 & bdryMatrix() const {
      return m_bdryMatrix;
    }
 
    inline const Eigen::VectorXd & bdryValues() const {
      return m_bdryValues;
    }
 
    inline const double & stabilizationParameter() const {
      return m_stab_par;
    }
 
    inline double & stabilizationParameter() {
      return m_stab_par;
    }
 
    KLNorms computeNorms(
                     const Eigen::VectorXd & v 
                    ) const;
 
    template<typename outValue, typename Fct>
    std::function<outValue(const Eigen::Vector2d &)> contractPara(const Fct &F) const 
    {
      std::function<outValue(const Eigen::Vector2d &)> f = [this, &F](const Eigen::Vector2d &x)->outValue { return F(x);}; 
      return f;
    }
 
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                      
                                const ForcingTermType & f,       
                                size_t iT_dirac                   
                                );
 
 
 
    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                  
                                      );
 
    const DDRCore & m_ddrcore;
    bool m_use_threads;
    std::ostream & m_output;
    VSXGrad m_vsxgrad;
    XHess m_xhess;
    BoundaryConditions m_BC_u;
    SystemMatrixType m_bdryMatrix;  // Matrix with columns corresponding to Dirichlet DOFs; multiplied by the boundary values, yields the modification to the system RHS for non-homogeneous BC
    Eigen::VectorXd m_bdryValues;   // Boundary values
    SystemMatrixType m_A;           // Matrix and RHS for system (after removing Dirichlet DOFs)
    Eigen::VectorXd m_b;
    double m_stab_par;
    std::vector<std::pair<size_t,size_t>> m_locUKN;    // Indicates the location, among the DOFs, of the unknowns after the Dirichlet DOFs are removed (the unknowns are between m_locUCK[i].first and m_locUCK[i].first+m_locUKN[i].second for all i)
    Eigen::VectorXi m_DOFtoUKN;       // Maps a dof i to its unknown in the system without BC, or to -1 if the DOF is a Dirichlet DOF
  };
 
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
  using std::exp;
  using std::pow;
 
  //------------------------------------------------------------------------------
  //    Constant
  static KirchhoffLove::SolutionDisplacementType  
  constant_u = [](const VectorRd & x) -> double {
                     return 1.;
                   };
  
  static KirchhoffLove::GradientDisplacementType
  constant_grad_u = [](const VectorRd & /*x*/) -> VectorRd{
                    return VectorRd::Zero();
                };
 
  static KirchhoffLove::ForcingTermType
  constant_f = [](const VectorRd & x) -> double {
                 return 0.;
               };
 
  //------------------------------------------------------------------------------
  //    Polynomial 
 
  static KirchhoffLove::SolutionDisplacementType
  polynomial_u = [](const VectorRd & x) -> double {
    return pow(x(0), 5) + pow(x(1), 5);
  };
 
  static KirchhoffLove::GradientDisplacementType
  polynomial_grad_u = [](const VectorRd & x) -> VectorRd {
      VectorRd grad;
      grad(0) = 5.0 * pow(x(0), 4);
      grad(1) = 5.0 * pow(x(1), 4);
  return grad;
  };
  
  static KirchhoffLove::ForcingTermType
  polynomial_f = [](const VectorRd & x) -> double {
                 return 120. * (x(0) + x(1));
               };
 
//------------------------------------------------------------------------------
//    Trigonometric solution
 
static KirchhoffLove::SolutionDisplacementType
trigo_u = [](const VectorRd & x) -> double {
  return std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
};
 
static KirchhoffLove::GradientDisplacementType
trigo_grad_u = [](const VectorRd & x) -> VectorRd {
  VectorRd grad;
  grad(0) = M_PI * std::cos(M_PI * x(0)) * std::sin(M_PI * x(1));
  grad(1) = M_PI * std::sin(M_PI * x(0)) * std::cos(M_PI * x(1));
  return grad;
};
 
static KirchhoffLove::ForcingTermType
trigo_f = [](const VectorRd & x) -> double {
  double u = std::sin(M_PI * x(0)) * std::sin(M_PI * x(1));
  return 4.0 * std::pow(M_PI, 4) * u;
};
 
//------------------------------------------------------------------------------
//    Dirac function f
 
static KirchhoffLove::SolutionDisplacementType
dirac_u = [](const VectorRd & x) -> double {
  return 1.0;
};
 
static KirchhoffLove::GradientDisplacementType
dirac_grad_u = [](const VectorRd & x) -> VectorRd {
  VectorRd grad;
  grad(0) = 0.0;
  grad(1) = 0.0;
  return grad;
};
 
// static KirchhoffLove::ForcingTermType
// dirac_f = [](const VectorRd & x) -> double {
//   return 0.0;
// };
 
//Dirac-like solution 
static KirchhoffLove::ForcingTermType
dirac_f = [](const VectorRd & x) -> double {
  const double rayon = 1e-2;
  const VectorRd center(0.3, 0.7);
 
  if ((x - center).squaredNorm() < rayon*rayon) {
      return 1.0 / (M_PI*rayon*rayon);  // to normalise //
  } else {
      return 0.0;
  }
};
} // end of namespace HArDCore2D
 
#endif