xdivdiv-test.hpp

Go to the documentation of this file.

// Test of the xdivdiv class
//
// Author: Daniele Di Pietro (daniele.di-pietro@umontpellier.fr), Jerome Droniou (jerome.droniou@monash.edu)
//
 
 
#ifndef XDIVDIV_TEST_HPP
#define XDIVDIV_TEST_HPP
 
#include <iostream>
 
#include <boost/math/constants/constants.hpp>
 
#include <Eigen/Sparse>
#include <unsupported/Eigen/SparseExtra>
 
#include <mesh.hpp>
 
#include <xdivdiv.hpp>
 
 
namespace HArDCore2D
{
 
   
  struct KirchhoffLove
  {
    typedef Eigen::SparseMatrix<double> SystemMatrixType;
    
    typedef std::function<double(const Eigen::Vector2d &)> ForcingTermType;
    typedef std::function<double(const Eigen::Vector2d &)> DeflectionType;
    typedef std::function<Eigen::Matrix2d(const Eigen::Vector2d &)> MomentTensorType;
    typedef std::function<double(const Eigen::Vector2d &, const Edge &)> MomentTensorEdgeDerivativeType;
 
 
    KirchhoffLove(
                   const PlatesCore & platescore,       
                   bool use_threads,                    
                   std::ostream & output = std::cout    
                   );
 
    inline size_t dimensionSpace() const
    {
      return m_xdivdiv.dimension() + m_platescore.mesh().n_cells() * PolynomialSpaceDimension<Cell>::Poly(m_platescore.degree()-2);
    }
 
    inline const XDivDiv & xDivDiv() const {
      return m_xdivdiv;
    }
    
    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;
    }
    
  private:
    const PlatesCore & m_platescore;
//    bool m_use_threads;
    std::ostream & m_output;
 
    XDivDiv m_xdivdiv;
    
    SystemMatrixType m_A;
    Eigen::VectorXd m_b;   
 
    double m_stab_par;
  };
 
  //------------------------------------------------------------------------------
  // Exact solutions
  //------------------------------------------------------------------------------
 
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
  using std::exp;
  using std::pow;
 
  //------------------------------------------------------------------------------
  // Trigonometric
  
  static KirchhoffLove::DeflectionType
  trigonometric_u = [](const VectorRd & x) -> double {
    return sin(PI*x(0))*sin(PI*x(1));
  };
 
  static KirchhoffLove::ForcingTermType
  trigonometric_f = [](const VectorRd & x) -> double {
    return 4.*pow(PI, 4)*sin(PI*x(0))*sin(PI*x(1));
  };
 
  static KirchhoffLove::MomentTensorType
  trigonometric_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
    Eigen::Matrix2d M;
    M.row(0) << -sin(PI*x(0))*sin(PI*x(1)), cos(PI*x(0))*cos(PI*x(1));
    M.row(1) << cos(PI*x(0))*cos(PI*x(1)), -sin(PI*x(0))*sin(PI*x(1));
 
    return -pow(PI, 2) * M;
  };
 
  static KirchhoffLove::MomentTensorEdgeDerivativeType
  trigonometric_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
 
    Eigen::Vector2d grad_sigma11 {
      -cos(PI*x(0))*sin(PI*x(1)), -sin(PI*x(0))*cos(PI*x(1))
    };
    grad_sigma11 *= -pow(PI, 3);
 
    Eigen::Vector2d grad_sigma12 {
      -sin(PI*x(0))*cos(PI*x(1)), -sin(PI*x(1))*cos(PI*x(0))
    };
    grad_sigma12 *= -pow(PI, 3);
 
    Eigen::Vector2d grad_sigma22 {
      -cos(PI*x(0))*sin(PI*x(1)), -sin(PI*x(0))*cos(PI*x(1))
    };
    grad_sigma22 *= -pow(PI, 3);
 
    Eigen::Vector2d tE = E.tangent();
    Eigen::Vector2d nE = E.normal();
 
    Eigen::Matrix2d dtE_sigma;
    dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
    dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
 
    Eigen::Vector2d div_sigma {
      grad_sigma11(0) + grad_sigma12(1),
      grad_sigma12(0) + grad_sigma22(1)
    };
 
    return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
  };
  
  //------------------------------------------------------------------------------
  // Simple
  
  static KirchhoffLove::DeflectionType
  simple_u = [](const VectorRd & x) -> double {
    return 1.;
  };
 
//  static KirchhoffLove::GradientDeflectionType
//  simple_grad_u = [](const VectorRd & x) -> VectorRd {
//    return VectorRd(0., 0.);
//  };
 
  static KirchhoffLove::ForcingTermType
  simple_f = [](const VectorRd & x) -> double {
    return 0.;
  };
 
  static KirchhoffLove::MomentTensorType
  simple_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
    // Here, M=hess(u) (modulo factor pi^2) and sigma is the opposite
    Eigen::Matrix2d M;
    M.row(0) << 0., 0.;
    M.row(1) << 0., 0.;
 
    return M;
  };
 
  static KirchhoffLove::MomentTensorEdgeDerivativeType
  simple_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
 
    Eigen::Vector2d grad_sigma11 {
      0., 0.
    };
    grad_sigma11 *= -pow(PI, 3);
 
    Eigen::Vector2d grad_sigma12 {
      0., 0.
    };
    grad_sigma12 *= -pow(PI, 3);
 
    Eigen::Vector2d grad_sigma22 {
      0., 0.
    };
    grad_sigma22 *= -pow(PI, 3);
 
    Eigen::Vector2d tE = E.tangent();
    Eigen::Vector2d nE = E.normal();
 
    Eigen::Matrix2d dtE_sigma;
    dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
    dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
 
    Eigen::Vector2d div_sigma {
      grad_sigma11(0) + grad_sigma12(1),
      grad_sigma12(0) + grad_sigma22(1)
    };
 
    return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
  };
 
  //------------------------------------------------------------------------------
  // Quartic
  
  static KirchhoffLove::DeflectionType
  quartic_u = [](const VectorRd & x) -> double {
    return pow(x(0),2)*pow(1.-x(0),2) + pow(x(1),2)*pow(1.-x(1),2);
  };
 
  static KirchhoffLove::ForcingTermType
  quartic_f = [](const VectorRd & x) -> double {
    return 48.;
  };
 
  static KirchhoffLove::MomentTensorType
  quartic_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
 
    Eigen::Matrix2d M;
    M <<
      2.*(pow(x(0),2)+4.*x(0)*(x(0)-1.)+pow(1.-x(0),2)), 0.,
      0., 2.*(pow(x(1),2)+4.*x(1)*(x(1)-1.)+pow(1.-x(1),2));
 
    return -M;
  };
 
  static KirchhoffLove::MomentTensorEdgeDerivativeType
  quartic_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
 
    Eigen::Vector2d grad_sigma11 { -24.*x(0)+12., 0. };
    Eigen::Vector2d grad_sigma12 { 0., 0. };
    Eigen::Vector2d grad_sigma22 { 0., -24.*x(1)+12. };
 
    Eigen::Vector2d tE = E.tangent();
    Eigen::Vector2d nE = E.normal();
 
    Eigen::Matrix2d dtE_sigma;
    dtE_sigma <<
      grad_sigma11.dot(tE), grad_sigma12.dot(tE),
      grad_sigma12.dot(tE), grad_sigma22.dot(tE);
 
    Eigen::Vector2d div_sigma {
      grad_sigma11(0) + grad_sigma12(1),
      grad_sigma12(0) + grad_sigma22(1)
    };
 
    return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
  };
 
  //------------------------------------------------------------------------------
  // Biquartic
  
  static KirchhoffLove::DeflectionType
  biquartic_u = [](const VectorRd & x) -> double {
    return pow(x(0),2)*pow(1.-x(0),2)*pow(x(1),2)*pow(1.-x(1),2);
  };
 
  static KirchhoffLove::ForcingTermType
  biquartic_f = [](const VectorRd & x) -> double {
    return 24*pow(x(0), 2)*pow(x(0) - 1, 2) + 32*x(0)*x(1)*(x(0) - 1)*(x(1) - 1) + 8*x(0)*x(1)*(x(0) - 1)*(4*x(1) - 2) + 8*x(0)*x(1)*(4*x(0) - 2)*(x(1) - 1) + 8*x(0)*x(1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 8*x(0)*(x(0) - 1)*(x(1) - 1)*(4*x(1) - 2) + 8*x(0)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 24*pow(x(1), 2)*pow(x(1) - 1, 2) + 8*x(1)*(x(0) - 1)*(4*x(0) - 2)*(x(1) - 1) + 8*x(1)*(x(0) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 2*(4*x(0) - 4)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1));
  };
 
  static KirchhoffLove::MomentTensorType
  biquartic_sigma = [](const VectorRd & x) -> Eigen::Matrix2d {
    Eigen::Matrix2d M;
      M.row(0) << 2*pow(x(1), 2)*pow(x(1) - 1, 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2)), 4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1));      
      M.row(1) << 4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)), 2*pow(x(0), 2)*pow(x(0) - 1, 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2));
    return -M;
  };
 
  static KirchhoffLove::MomentTensorEdgeDerivativeType
  biquartic_sigma_DE = [](const Eigen::Vector2d & x, const Edge & E) -> double {
 
   Eigen::Vector2d grad_sigma11 {
      2*pow(x(1), 2)*(12*x(0) - 6)*pow(x(1) - 1, 2),
      2*pow(x(1), 2)*(2*x(1) - 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2)) + 4*x(1)*pow(x(1) - 1, 2)*(pow(x(0), 2) + 4*x(0)*(x(0) - 1) + pow(x(0) - 1, 2))
    };
 
    Eigen::Vector2d grad_sigma12 {
      4*x(0)*x(1)*(x(0) - 1)*(x(1) - 1)*(4*x(1) - 2) + 4*x(0)*x(1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 4*x(1)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)), 
      4*x(0)*x(1)*(x(0) - 1)*(4*x(0) - 2)*(x(1) - 1) + 4*x(0)*x(1)*(x(0) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1)) + 4*x(0)*(x(0) - 1)*(x(1) - 1)*(x(0)*x(1) + x(0)*(x(1) - 1) + x(1)*(x(0) - 1) + (x(0) - 1)*(x(1) - 1))
    };
 
    Eigen::Vector2d grad_sigma22 {
     2*pow(x(0), 2)*(2*x(0) - 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2)) + 4*x(0)*pow(x(0) - 1, 2)*(pow(x(1), 2) + 4*x(1)*(x(1) - 1) + pow(x(1) - 1, 2)), 
     2*pow(x(0), 2)*pow(x(0) - 1, 2)*(12*x(1) - 6)
    };
    
    // The previous calculations are actually grad_hess(u)
    grad_sigma11 *= -1.;
    grad_sigma12 *= -1.;
    grad_sigma22 *= -1.;
 
    Eigen::Vector2d tE = E.tangent();
    Eigen::Vector2d nE = E.normal();
 
    Eigen::Matrix2d dtE_sigma;
    dtE_sigma.row(0) << grad_sigma11.dot(tE), grad_sigma12.dot(tE);
    dtE_sigma.row(1) << grad_sigma12.dot(tE), grad_sigma22.dot(tE);
 
    Eigen::Vector2d div_sigma {
      grad_sigma11(0) + grad_sigma12(1),
      grad_sigma12(0) + grad_sigma22(1)
    };
 
    return (dtE_sigma*nE).dot(tE) + div_sigma.dot(nE);
  };
  
} // end of namespace HArDCore2D
 
#endif