sxcurl-test.hpp

Go to the documentation of this file.

#ifndef XCURL_TEST_HPP
#define XCURL_TEST_HPP
 
#include <boost/math/constants/constants.hpp>
 
namespace HArDCore2D
{
  
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
 
  //------------------------------------------------------------------------------
    static std::function<double(const Eigen::Vector2d&)>
  trigonometric_scalar = [](const Eigen::Vector2d & x) -> double {
               return sin(PI * x(0)) * sin(PI * x(1));
             };
  
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  grad_trigonometric_scalar = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
                return PI * Eigen::Vector2d(
                                cos(PI * x(0)) * sin(PI * x(1)),
                                sin(PI * x(0)) * cos(PI * x(1))
                                );
                  };
  //------------------------------------------------------------------------------
 
  static std::function<double(const Eigen::Vector2d&)>
  constant_scalar = [](const Eigen::Vector2d & x) -> double {
              return 1.;
            };
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  grad_constant_scalar = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
               return Eigen::Vector2d::Zero();
             };
 
  //------------------------------------------------------------------------------
 
  static std::function<double(const Eigen::Vector2d&)>
  linear_scalar = [](const Eigen::Vector2d & x) -> double {
            return 1. + x(0) + 2. * x(1);
          };
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  grad_linear_scalar = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
             return Eigen::Vector2d(1., 2.);
               };
 
  //------------------------------------------------------------------------------
  
  static std::function<double(const Eigen::Vector2d&)>
  quadratic_scalar = [](const Eigen::Vector2d & x) -> double {
               return linear_scalar(x) + std::pow(x(0), 2) + 2. * std::pow(x(1), 2);
             };
 
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  grad_quadratic_scalar = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
                return grad_linear_scalar(x) + Eigen::Vector2d(2. * x(0), 4. * x(1));
              };
 
 
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------  
  static std::function<Eigen::Vector2d(const Eigen::Vector2d&)>
  constant_vector = [](const Eigen::Vector2d & x) -> Eigen::Vector2d {
              return Eigen::Vector2d(1., -1.);
            };
 
  static std::function<double(const Eigen::Vector2d&)>
  rot_constant_vector = [](const Eigen::Vector2d & x) -> double {
               return 0.;
             };
 
  //------------------------------------------------------------------------------
 
  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<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));
    };
 
  //------------------------------------------------------------------------------
 
  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
    //std::cout<<"fx_basis 0: "<<fX_basis_quad.shape()[0] <<", fxsize: "<<fX.size()<<", fx_basis 1: "<<fX_basis_quad.shape()[1]<<", quad size: "<< quad_X.size()<<std::endl;
    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