vsxcurl-test.hpp

Go to the documentation of this file.

#ifndef VSXCURL_TEST_HPP
#define VSXCURL_TEST_HPP
 
#include <boost/math/constants/constants.hpp>
#include <vsxcurl.hpp>
 
namespace HArDCore2D
{
static const double PI = boost::math::constants::pi<double>();
using std::sin;
 
//------------------------------------------------------------------------------
 
static std::function<Eigen::Matrix2d(const Eigen::Vector2d&)> constant_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Matrix2d {
  return Eigen::Matrix2d::Constant(1.);
};
 
static std::function<Eigen::Vector2d(const Eigen::Vector2d&)> rot_constant_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Vector2d {
  return Eigen::Vector2d::Zero();
};
 
//------------------------------------------------------------------------------
 
static std::function<Eigen::Matrix2d(const Eigen::Vector2d&)> linear_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Matrix2d {
  Eigen::Matrix2d F;
  F << 1. + x(0), x(1),
       2. * x(1), 3. - x(0);
  return F;
};
 
static std::function<Eigen::Vector2d(const Eigen::Vector2d&)> rot_linear_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Vector2d {
  return Eigen::Vector2d(0., -3.);
};
 
//------------------------------------------------------------------------------
 
static std::function<Eigen::Matrix2d(const Eigen::Vector2d&)> trigonometric_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Matrix2d {
  Eigen::Matrix2d F;
  double s = sin(PI * x(0)) * sin(PI * x(1));
  double c = cos(PI * x(0)) * cos(PI * x(1));
  F << s,  c,
       -c, s;
  return F;
};
 
static std::function<Eigen::Vector2d(const Eigen::Vector2d&)> rot_trigonometric_tensor = 
[](const Eigen::Vector2d& x) -> Eigen::Vector2d {
  double px = PI * x(0);
  double py = PI * x(1);
  double comp1 = -2.0 * PI * sin(px) * cos(py);
  double comp2 =  0.0;
  return Eigen::Vector2d(comp1, comp2);
};
//------------------------------------------------------------------------------
 
template<typename T>
double squared_l2_error(
  const std::function<T(const Eigen::Vector2d &)> & f,
  const Eigen::VectorXd & fX,
  const boost::multi_array<T, 2> & fX_basis_quad,
  const QuadratureRule & quad_X
)
{
  assert(fX_basis_quad.shape()[0] == (size_t)fX.size() &&
         fX_basis_quad.shape()[1] == quad_X.size());
  
  double err = 0.;
  for (size_t iqn = 0; iqn < quad_X.size(); iqn++) {
    T f_iqn = f(quad_X[iqn].vector());
    
    T fX_iqn = fX(0) * fX_basis_quad[0][iqn];
    for (size_t i = 1; i < fX_basis_quad.shape()[0]; i++) {
      fX_iqn += fX(i) * fX_basis_quad[i][iqn];
    }
 
    T diff_iqn = f_iqn - fX_iqn;
    err += quad_X[iqn].w * diff_iqn.squaredNorm(); 
  }
  return err;
}
}
 
#endif