bgg-rmplate.hpp

Go to the documentation of this file.

// Implementation of the discrete de Rham sequence for the Reissner-Mindlin plate problem.
//
// Author: Arax Leroy (arax.leroy@umontpellier.fr)
//
 
#ifndef BGG_RMPLATE_HPP
#define BGG_RMPLATE_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 <vsxcurl.hpp>
#include <vsxgrad.hpp>
#include <xhess_full.hpp>
 
namespace HArDCore2D
{
 
  struct RMParameters
  {
    RMParameters( const double thickness, 
                  const double young_modulus,
                  const double poisson_ratio
                ):
             t(thickness),
             E(young_modulus),
             nu(poisson_ratio)
             {
                // Compute derived quantities
                beta0 = E/(12*(1.+nu));
                beta1 = E*nu/(12*(1.-std::pow(nu,2)));
                kappa = 5*E/(12*(1.+nu));
             }
       
    double t; // plate thickness
    double E; // Young modulus
    double nu; // Poisson ratio
    double beta0;
    double beta1;
    double kappa;
  };
 
  struct RMNorms
  {
    RMNorms(double norm_rotation, double norm_displacement, double norm_kirchoff):
      rotation(norm_rotation),
      displacement(norm_displacement),
      kirchoff(norm_kirchoff),
      energy(std::sqrt( std::pow(norm_rotation, 2)+std::pow(norm_displacement, 2)+std::pow(norm_kirchoff, 2) ) )
      {
        // do nothing
      };
    
    double rotation; 
    double displacement; 
    double kirchoff; 
    double energy; 
  
  };
 
  struct ReissnerMindlin
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<double(const RMParameters &, const Eigen::Vector2d &)> ForcingTermType;
    typedef std::function<Eigen::Vector2d(const RMParameters &, const Eigen::Vector2d &)> SolutionRotationType;
    typedef std::function<Eigen::Matrix2d(const RMParameters &, const Eigen::Vector2d &)> GradientRotationType;
    typedef std::function<double(const RMParameters &, const Eigen::Vector2d &)> SolutionDisplacementType;
    typedef std::function<Eigen::Vector2d(const RMParameters&, const Eigen::Vector2d&)> GradientDisplacementType;
 
 
    ReissnerMindlin(
                   const DDRCore & ddrcore,               
                   const DDRCore & ddrcore_plus,          
                   const SerendipityProblem & sp,
                   const SerendipityProblem & sp_plus,    
                   const RMParameters & para,             
                   const BoundaryConditions & BC_theta,   
                   const BoundaryConditions & BC_u,       
                   bool use_threads,                      
                   std::ostream & output = std::cout      
                   );
 
    void assembleLinearSystem(
                              const ForcingTermType & f,      
                              const SolutionRotationType & theta,  
                              const GradientRotationType & grad_theta,  
                              const SolutionDisplacementType & u,  
                              const GradientDisplacementType & grad_u 
                              );
 
 
    inline size_t dimensionSpace() const
    {
      return m_vsxgrad.dimension() + m_xhess.dimension();
    }
 
    inline size_t nb_bdryDOFs() const
    { 
      return m_BC_theta.n_dir_edges() * m_vsxgrad.numLocalDofsEdge(0)
             + m_BC_theta.n_dir_vertices() * m_vsxgrad.numLocalDofsVertex(0)
             + m_BC_u.n_dir_edges() * (m_xhess.numLocalDofsEdge() - m_ddrcore.edgeBases(0).Polyk->dimension())
             + m_BC_u.n_dir_vertices() * (m_xhess.numLocalDofsVertex()-2); 
    }
 
    inline size_t sizeSystem() const
    {
      return dimensionSpace() - nb_bdryDOFs();
    }
 
inline const std::vector<std::pair<size_t,size_t>> & locUKN() const {
  return m_locUKN;
}
 
inline const Eigen::VectorXi & DOFtoUKN() const {
  return m_DOFtoUKN;
}
 
inline const BoundaryConditions & BC_u() const {
  return m_BC_u;
}
 
std::vector<size_t> globalDOFIndices(const Cell &T) const
{
      std::vector<size_t> I_vsxgrad_T = m_vsxgrad.globalDOFIndices(T);
      std::vector<size_t> I_xhess_T = m_xhess.globalDOFIndices(T);
      size_t dim_T = m_vsxgrad.dimensionCell(T.global_index()) + m_xhess.dimensionCell(T.global_index());
      std::vector<size_t> I_T(dim_T);
      size_t dim_vsxgrad = m_vsxgrad.dimension();
      auto it_I_T = std::copy(I_vsxgrad_T.begin(), I_vsxgrad_T.end(), I_T.begin());
      std::transform(I_xhess_T.begin(), I_xhess_T.end(), it_I_T, [&dim_vsxgrad](const size_t & index) { return index + dim_vsxgrad; });
      return I_T;
    }
 
    inline const RMParameters & para() const
    {
      return m_para;
    }
    
    inline const VSXGrad & vsxGrad() const
    {
      return m_vsxgrad;
    }
 
    inline const VSXCurl & vsxCurl() const
    {
      return m_vsxcurl;
    }
 
 
    inline const XHessFull & 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;
    }
 
    RMNorms 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(m_para,x);};
      return f;
    }
    
 
  private:
    std::pair<Eigen::MatrixXd, Eigen::VectorXd>
    _compute_local_contribution(
                                size_t iT,                      
                                const GradientRotationType & grad_theta,      
                                const ForcingTermType & f       
                                );
 
 
 
    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;
    const RMParameters m_para;
    bool m_use_threads;
    std::ostream & m_output;
    VSXGrad m_vsxgrad;
    VSXCurl m_vsxcurl;
    XHessFull m_xhess;
    BoundaryConditions m_BC_theta;
    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
//------------------------------------------------------------------------------
//    Test complex: u = x^2 / 2, grad_u = (x,0)
//------------------------------------------------------------------------------
static ReissnerMindlin::SolutionRotationType
constant_theta = [](const RMParameters &, const VectorRd &) -> VectorRd {
    return VectorRd(1, 0.);
};
 
static ReissnerMindlin::GradientRotationType
constant_grad_theta = [](const RMParameters &, const VectorRd &) -> Eigen::Matrix2d {
    return Eigen::Matrix2d::Zero();
};
 
static ReissnerMindlin::SolutionDisplacementType
constant_u = [](const RMParameters &, const VectorRd & x) -> double {
    return x(0);
};
 
static ReissnerMindlin::GradientDisplacementType
constant_grad_u = [](const RMParameters &, const VectorRd &) -> VectorRd {
    VectorRd g;
    g.setZero();
    g(0) = 1.0;
    return g;
};
 
static ReissnerMindlin::ForcingTermType
constant_f = [](const RMParameters &, const VectorRd &) -> double {
    return 0.0;
};
  //------------------------------------------------------------------------------
  // Polynomial
  
static ReissnerMindlin::SolutionRotationType
polynomial_theta = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
                 return VectorRd(
                          pow(x(1),3)*pow(x(1)-1,3)*pow(x(0),2)*pow(x(0)-1,2)*(2*x(0)-1),
                          pow(x(0),3)*pow(x(0)-1,3)*pow(x(1),2)*pow(x(1)-1,2)*(2*x(1)-1)
                          );
               };
               
static ReissnerMindlin::GradientRotationType
polynomial_grad_theta = [](const RMParameters & para, const VectorRd & x) -> Eigen::Matrix2d {
    const double x0 = x(0);
    const double x1 = x(1);
 
    // theta_0 = A(x0) B(x1)
    const double A =
        std::pow(x0, 2) * std::pow(x0 - 1.0, 2) * (2.0 * x0 - 1.0);
 
    const double dA =
        2.0 * x0 * std::pow(x0 - 1.0, 2) * (2.0 * x0 - 1.0)
      + std::pow(x0, 2) * 2.0 * (x0 - 1.0) * (2.0 * x0 - 1.0)
      + std::pow(x0, 2) * std::pow(x0 - 1.0, 2) * 2.0;
 
    const double B =
        std::pow(x1, 3) * std::pow(x1 - 1.0, 3);
 
    const double dB =
        3.0 * std::pow(x1, 2) * std::pow(x1 - 1.0, 3)
      + 3.0 * std::pow(x1, 3) * std::pow(x1 - 1.0, 2);
 
    // theta_1 = C(x0) D(x1)
    const double C =
        std::pow(x0, 3) * std::pow(x0 - 1.0, 3);
 
    const double dC =
        3.0 * std::pow(x0, 2) * std::pow(x0 - 1.0, 3)
      + 3.0 * std::pow(x0, 3) * std::pow(x0 - 1.0, 2);
 
    const double D =
        std::pow(x1, 2) * std::pow(x1 - 1.0, 2) * (2.0 * x1 - 1.0);
 
    const double dD =
        2.0 * x1 * std::pow(x1 - 1.0, 2) * (2.0 * x1 - 1.0)
      + std::pow(x1, 2) * 2.0 * (x1 - 1.0) * (2.0 * x1 - 1.0)
      + std::pow(x1, 2) * std::pow(x1 - 1.0, 2) * 2.0;
 
    Eigen::Matrix2d G = Eigen::Matrix2d::Zero();
 
    G(0,0) = dA * B;  // d theta_0 / dx
    G(0,1) = A * dB;  // d theta_0 / dy
    G(1,0) = dC * D;  // d theta_1 / dx
    G(1,1) = C * dD;  // d theta_1 / dy
 
    return G;
};        
 
  static ReissnerMindlin::SolutionDisplacementType  
  polynomial_u = [](const RMParameters & para, const VectorRd & x) -> double {
                 double val = 0;
                 
                 val += (1./3.)*pow(x(0),3)*pow(x(0)-1,3)*pow(x(1),3)*pow(x(1)-1,3);
                 
                 val -= ( 2*pow(para.t,2) / (5*(1.-para.nu)) ) *
                          ( pow(x(1),3)*pow(x(1)-1,3)*x(0)*(x(0)-1)*(5*x(0)*x(0)-5*x(0)+1)
                            + pow(x(0),3)*pow(x(0)-1,3)*x(1)*(x(1)-1)*(5*x(1)*x(1)-5*x(1)+1) );
                 
                 return val;
                 };
 
  static ReissnerMindlin::ForcingTermType
  polynomial_f = [](const RMParameters & para, const VectorRd & x) -> double {
                double val = 0;
                
                val += 12*x(1)*(x(1)-1)*(5*x(0)*x(0)-5*x(0)+1) * 
                        ( 2*x(1)*x(1)*(x(1)-1)*(x(1)-1) + x(0)*(x(0)-1)*(5*x(1)*x(1)-5*x(1)+1) );
                val += 12*x(0)*(x(0)-1)*(5*x(1)*x(1)-5*x(1)+1) * 
                        ( 2*x(0)*x(0)*(x(0)-1)*(x(0)-1) + x(1)*(x(1)-1)*(5*x(0)*x(0)-5*x(0)+1) );
                
                return val * para.E / (12* (1.-para.nu*para.nu) );
               };
 
static ReissnerMindlin::GradientDisplacementType 
polynomial_grad_u = [](const RMParameters & para, const VectorRd & x) -> VectorRd
{
    const double t2 = para.t * para.t;
    const double c  = 2.0 * t2 / (5.0 * (1.0 - para.nu));
 
    const double x0 = x(0);
    const double x1 = x(1);
 
    const double P0  = pow(x0,3) * pow(x0 - 1.0,3);
    const double P1  = pow(x1,3) * pow(x1 - 1.0,3);
 
    const double dP0 =
        3.0 * pow(x0,2) * pow(x0 - 1.0,3)
      + 3.0 * pow(x0,3) * pow(x0 - 1.0,2);
 
    const double dP1 =
        3.0 * pow(x1,2) * pow(x1 - 1.0,3)
      + 3.0 * pow(x1,3) * pow(x1 - 1.0,2);
 
    const double Q0 = x0 * (x0 - 1.0) * (5.0*x0*x0 - 5.0*x0 + 1.0);
    const double Q1 = x1 * (x1 - 1.0) * (5.0*x1*x1 - 5.0*x1 + 1.0);
 
    const double dQ0 =
        (x0 - 1.0) * (5.0*x0*x0 - 5.0*x0 + 1.0)
      + x0 * (5.0*x0*x0 - 5.0*x0 + 1.0)
      + x0 * (x0 - 1.0) * (10.0*x0 - 5.0);
 
    const double dQ1 =
        (x1 - 1.0) * (5.0*x1*x1 - 5.0*x1 + 1.0)
      + x1 * (5.0*x1*x1 - 5.0*x1 + 1.0)
      + x1 * (x1 - 1.0) * (10.0*x1 - 5.0);
 
    VectorRd g;
    g.setZero();
 
    // d/dx0
    g(0) =
        (1.0 / 3.0) * dP0 * P1
      - c * ( P1 * dQ0 + dP0 * Q1 );
 
    // d/dx1
    g(1) =
        (1.0 / 3.0) * P0 * dP1
      - c * ( dP1 * Q0 + P0 * dQ1 );
 
    return g;
};
  //------------------------------------------------------------------------------
  // Analytical solution (see paper)
  //  (BC are not zero here)  
 
  static std::function<double(const VectorRd &)> 
          an_g = [](const VectorRd &x)->double { return sin(PI*x(0))*sin(PI*x(1)); };
          
  static std::function<VectorRd(const VectorRd &)>
          an_GRAD_g = [](const VectorRd & x) -> VectorRd {
                 return VectorRd( PI*cos(PI*x(0))*sin(PI*x(1)), PI*sin(PI*x(0))*cos(PI*x(1)) ); 
               };
               
  static std::function<MatrixRd(const VectorRd &)>
          an_HESS_g = [](const VectorRd & x) -> MatrixRd {
                MatrixRd H = MatrixRd::Zero();
                H.row(0) << -PI*PI*sin(PI*x(0))*sin(PI*x(1)), PI*PI*cos(PI*x(0))*cos(PI*x(1));
                H.row(1) << PI*PI*cos(PI*x(0))*cos(PI*x(1)), -PI*PI*sin(PI*x(0))*sin(PI*x(1));
                return H;
              };              
 
  static std::function<double(const VectorRd &)> 
          an_LAPL_g = [](const VectorRd &x)->double { return -2*pow(PI,2)*sin(PI*x(0))*sin(PI*x(1)); };
 
  static std::function<double(const VectorRd &)> 
          an_V = [](const VectorRd &x)->double { return x(0) * exp(-x(0))*cos(x(1)); };
 
  static std::function<VectorRd(const VectorRd &)>
          an_GRAD_V = [](const VectorRd & x) -> VectorRd {
                 return VectorRd((1.-x(0))*exp(-x(0))*cos(x(1)), -x(0)*exp(-x(0))*sin(x(1))); 
               };
            
  static std::function<MatrixRd(const VectorRd &)>
          an_HESS_V = [](const VectorRd & x) -> MatrixRd {
                MatrixRd H = MatrixRd::Zero();
                H.row(0) << (x(0)-2.)*exp(-x(0))*cos(x(1)), -(1.-x(0))*exp(-x(0))*sin(x(1));
                H.row(1) << -(1.-x(0))*exp(-x(0))*sin(x(1)), -x(0)*exp(-x(0))*cos(x(1));
                return H;
              };              
 
  static std::function<double(const VectorRd &)> 
          an_LAPL_V = [](const VectorRd &x)->double { return -2.* exp(-x(0))*cos(x(1)); };
 
 
  static ReissnerMindlin::SolutionDisplacementType  
          an_v = [](const RMParameters & para, const VectorRd & x) -> double {
                     return pow(para.t,3) * an_V(x/para.t) + an_g(x);
                   };
 
  static ReissnerMindlin::SolutionRotationType
          an_GRAD_v = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
                 return pow(para.t, 2) * an_GRAD_V(x/para.t) + an_GRAD_g(x);
               };
  
  static ReissnerMindlin::GradientRotationType
          an_HESS_v = [](const RMParameters & para, const VectorRd & x) -> MatrixRd {
                  return para.t * an_HESS_V(x/para.t) + an_HESS_g(x);
              };
          
  static ReissnerMindlin::SolutionDisplacementType  
          an_LAPL_v = [](const RMParameters & para, const VectorRd & x) -> double {
                     return para.t * (an_LAPL_V)(x/para.t) + an_LAPL_g(x);
                   };
 
  static ReissnerMindlin::SolutionRotationType analytical_theta = an_GRAD_v;
  
  static ReissnerMindlin::GradientRotationType analytical_grad_theta = an_HESS_v;
 
  static ReissnerMindlin::SolutionDisplacementType  
  analytical_u = [](const RMParameters & para, const VectorRd & x) -> double {
                     return an_v(para, x) - pow(para.t,2) * ( (para.beta0+para.beta1)/para.kappa ) * an_LAPL_v(para, x);
                   };
 
  static ReissnerMindlin::ForcingTermType
  analytical_f = [](const RMParameters & para, const VectorRd & x) -> double {
                 return (para.beta0+para.beta1) * 4 * pow(PI,4)*sin(PI*x(0))*sin(PI*x(1));
               };
 
 // ajout 
 static std::function<VectorRd(const VectorRd &)>
an_GRAD_LAPL_g = [](const VectorRd & x) -> VectorRd {
  return VectorRd(
    -2.0 * pow(PI,3) * cos(PI*x(0)) * sin(PI*x(1)),
    -2.0 * pow(PI,3) * sin(PI*x(0)) * cos(PI*x(1))
  );
};
 
static std::function<VectorRd(const VectorRd &)>
an_GRAD_LAPL_V = [](const VectorRd & x) -> VectorRd {
  return VectorRd(
    2.0 * exp(-x(0)) * cos(x(1)),
    2.0 * exp(-x(0)) * sin(x(1))
  );
};
 
static ReissnerMindlin::SolutionRotationType
an_GRAD_LAPL_v = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
  return an_GRAD_LAPL_V(x/para.t) + an_GRAD_LAPL_g(x);
};
 
 
static ReissnerMindlin::GradientDisplacementType
analytical_grad_u = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
  return an_GRAD_v(para, x)
       - pow(para.t, 2) * ((para.beta0 + para.beta1) / para.kappa)
         * an_GRAD_LAPL_v(para, x);
};
 
 
 
  //------------------------------------------------------------------------------
  // Unkown exact solution, just useful to fix BC and load and plot the approximation solution
  
  static ReissnerMindlin::SolutionRotationType
  ukn_theta = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
                 return VectorRd(0., 0.);
               };
 
  static ReissnerMindlin::GradientRotationType
  ukn_grad_theta = [](const RMParameters & para, const VectorRd & x) -> Eigen::Matrix2d {
                 return Eigen::Matrix2d::Zero();
               };               
 
  static ReissnerMindlin::SolutionDisplacementType  
  ukn_u = [](const RMParameters & para, const VectorRd & x) -> double {
                     return x(0);
                   };
  
  static ReissnerMindlin::GradientDisplacementType
  ukn_grad_u = [](const RMParameters & /*para*/, const VectorRd & /*x*/) -> VectorRd {
  VectorRd grad = VectorRd::Zero();
  grad(0) = 1.0;
                    return grad;
                  };
 
  static ReissnerMindlin::ForcingTermType
  ukn_f = [](const RMParameters & para, const VectorRd & x) -> double {
                double val = 1.;
                if (x(0)<.5){
                  val = -1.;
                }
                return val;
               };
 
  //------------------------------------------------------------------------------
  // Kirchoff limit: solution of E/(12(1-nu^2)) Delta^2 u = g
  //  Useful to evaluate convergence as t->0
  
  static ReissnerMindlin::SolutionRotationType
  kir_theta = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
                double coef = 12.*(1-pow(para.nu,2))/ (para.E * 4.*pow(PI, 4));
                return coef * PI * VectorRd(cos(PI*x(0))*sin(PI*x(1)), sin(PI*x(0))*cos(PI*x(1)));
               };
 
  static ReissnerMindlin::GradientRotationType
  kir_grad_theta = [](const RMParameters & para, const VectorRd & x) -> Eigen::Matrix2d {
                 double coef = 12.*(1-pow(para.nu,2))/ (para.E * 4.*pow(PI, 4));
                 Eigen::Matrix2d H = Eigen::Matrix2d::Zero();
                 H.row(0) << -PI*PI*sin(PI*x(0))*sin(PI*x(1)), PI*PI*cos(PI*x(0))*cos(PI*x(1));
                 H.row(1) <<  PI*PI*cos(PI*x(0))*cos(PI*x(1)), -PI*PI*sin(PI*x(0))*sin(PI*x(1));
                 return coef * H;
               };               
 
  static ReissnerMindlin::SolutionDisplacementType  
  kir_u = [](const RMParameters & para, const VectorRd & x) -> double {
                  double coef = 12.*(1-pow(para.nu,2))/ (para.E * 4.*pow(PI, 4));
                  return coef  * sin(PI*x(0))*sin(PI*x(1));
               };
static ReissnerMindlin::GradientDisplacementType
kir_grad_u = [](const RMParameters & para, const VectorRd & x) -> VectorRd {
                double coef = 12. * (1 - pow(para.nu, 2)) / (para.E * 4. * pow(PI, 4));
                VectorRd grad = VectorRd::Zero();
                grad(0) = coef * PI * cos(PI * x(0)) * sin(PI * x(1));
                grad(1) = coef * PI * sin(PI * x(0)) * cos(PI * x(1));
                return grad;
              };
 
  static ReissnerMindlin::ForcingTermType
  kir_f = [](const RMParameters & para, const VectorRd & x) -> double {
                return sin(PI*x(0))*sin(PI*x(1));
               };
 
 
} // end of namespace HArDCore2D
 
#endif