excurl-test.hpp

Go to the documentation of this file.

#ifndef EXCURL_TEST_HPP
#define EXCURL_TEST_HPP
 
#include <boost/math/constants/constants.hpp>
 
namespace HArDCore2D
{
  
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
  using std::cos;
 
  //------------------------------------------------------------------------------
  
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  constant_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
              return Eigen::Vector2d(1., 2.);
            };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_constant_vector = [](const Eigen::Vector2d & x) -> double {
               return 0.;
             };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_constant_vector = [](const Eigen::Vector2d & x) -> MatrixRd {
               return MatrixRd::Zero();
             };
  //------------------------------------------------------------------------------
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  linear_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
            return Eigen::Vector2d(
                       1. + x(0) - x(1),
                       2. + 2. * x(1) + x(0)
                       );
          };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_linear_vector = [](const Eigen::Vector2d & x) -> double {
             return 2;
               };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_linear_vector = [](const Eigen::Vector2d & x) -> MatrixRd {
         MatrixRd G = MatrixRd::Zero();
         G.row(0) << 1, -1;
         G.row(1) << 1, 2;         
               return G;
             };
 
 
  //------------------------------------------------------------------------------
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  quadratic_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
            return Eigen::Vector2d(
                       1. + x(0)*x(0) - x(1)*x(0),
                       2. + 2. * x(1)*x(1) + x(0)*x(1)
                       );
          };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_quadratic_vector = [](const Eigen::Vector2d & x) -> double {
             return x(0)+x(1);
               };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_quadratic_vector = [](const Eigen::Vector2d & x) -> MatrixRd {
         MatrixRd G = MatrixRd::Zero();
         G.row(0) << 2.*x(0)-x(1), -x(0);
         G.row(1) << x(1), 4.*x(1)+x(0);         
               return G;
             };
 
  //------------------------------------------------------------------------------
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  cubic_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
            return Eigen::Vector2d(
                       1. + std::pow(x(0),3) - x(1)*std::pow(x(0),2),
                       2. + 2. * std::pow(x(1), 3) + std::pow(x(1),2)*x(0)
                       );
          };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_cubic_vector = [](const Eigen::Vector2d & x) -> double {
             return x(0)*x(0)+x(1)*x(1);
               };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_cubic_vector = [](const Eigen::Vector2d & x) -> MatrixRd {
         MatrixRd G = MatrixRd::Zero();
         G.row(0) << 3*x(0)*x(0)-2*x(1)*x(0), -x(0)*x(0);
         G.row(1) << x(1)*x(1), 6*x(1)*x(1)+2*x(1)*x(0);         
               return G;
             };
 
  //------------------------------------------------------------------------------
  
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  trigonometric_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
               return Eigen::Vector2d(
                          sin(PI * x(0)) * sin(PI * x(1)),
                          sin(PI * x(0)) * sin(PI * x(1))
                          );
             };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_trigonometric_vector
  = [](const Eigen::Vector2d & x) -> double {
      return PI*cos(PI*x(0))*sin(PI*x(1)) - PI*sin(PI*x(0))*cos(PI*x(1));
    };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_trigonometric_vector = [](const Eigen::Vector2d & x) -> MatrixRd {
         MatrixRd G = MatrixRd::Zero();
         G.row(0) << PI*cos(PI*x(0))*sin(PI*x(1)), PI*sin(PI*x(0))*cos(PI*x(1));
         G.row(1) << PI*cos(PI*x(0))*sin(PI*x(1)), PI*sin(PI*x(0))*cos(PI*x(1));         
               return G;
             };
 
 
  //------------------------------------------------------------------------------
  // Scalar function for gradient
  static std::function<double(const Eigen::Vector2d&)>
  trigonometric_q = [](const Eigen::Vector2d & x) -> double {
               return sin(PI * x(0)) * sin(PI * x(1));
             };
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  grad_trigonometric_q
  = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
      return Eigen::Vector2d(PI*cos(PI*x(0))*sin(PI*x(1)), PI*sin(PI*x(0))*cos(PI*x(1)));
    };
 
  //------------------------------------------------------------------------------
 
  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
              )
  {
    // Check that the dimensions are consistent    
    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];
      } // for i
 
      T diff_iqn = f_iqn - fX_iqn;
 
      err += quad_X[iqn].w * scalar_product(diff_iqn, diff_iqn);
    } // for iqn
    return err;
  }
 
} // end of namespace HArDCore2D
 
#endif