vsxgrad-test.hpp

Go to the documentation of this file.

#ifndef VSXGRAD_TEST_HPP
#define VSXGRAD_TEST_HPP
 
#include <boost/math/constants/constants.hpp>
#include <vsxgrad.hpp>
 
namespace HArDCore2D
{
  
  static const double PI = boost::math::constants::pi<double>();
  using std::sin;
 
  //------------------------------------------------------------------------------
  
  static VSXGrad::FunctionType
  trigonometric = [](const VectorRd & x) -> VectorRd {
                           return sin(PI * x(0)) * sin(PI * x(1)) * VectorRd(1., 1.);
             };
  
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_trigonometric = [](const Eigen::Vector2d & x) -> MatrixRd {
              MatrixRd G = MatrixRd::Zero();
        G.row(0) << PI * Eigen::Vector2d(
                                cos(PI * x(0)) * sin(PI * x(1)),
                                sin(PI * x(0)) * cos(PI * x(1))
                                ).transpose();
        G.row(1) << PI * Eigen::Vector2d(
                                cos(PI * x(0)) * sin(PI * x(1)),
                                sin(PI * x(0)) * cos(PI * x(1))
                                ).transpose();
 
                return G;
                  };
 
  static std::function<double(const Eigen::Vector2d&)> 
  rot_trigonometric = [](const Eigen::Vector2d& x) -> double {
    double d1f2 = PI * std::cos(PI * x(0)) * std::sin(PI * x(1));
    double d2f1 = PI * std::sin(PI * x(0)) * std::cos(PI * x(1));
    return d1f2 - d2f1;
  };
  //------------------------------------------------------------------------------
 
  static std::function<VectorRd(const VectorRd&)>
  constant = [](const VectorRd & x) -> VectorRd {
              return VectorRd(1., 1.);
            };
 
  static std::function<MatrixRd(const VectorRd&)>
  grad_constant = [](const VectorRd & x) -> MatrixRd {
               return MatrixRd::Zero();
             };
  
static std::function<double(const VectorRd&)>
rot_constant = [](const VectorRd & x) -> double {
    return 0.;
            };
 
  //------------------------------------------------------------------------------
 
  static std::function<VectorRd(const VectorRd&)>
  linear = [](const VectorRd & x) -> VectorRd {
            return (1. + x(0) + 2. * x(1)) * VectorRd(-1., 1.);
          };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_linear = [](const Eigen::Vector2d & x) -> MatrixRd {
       MatrixRd G = MatrixRd::Zero();
       G.row(0) << -VectorRd(1., 2.).transpose();
       G.row(1) << VectorRd(1., 2.).transpose();
       return G;
               };
 
static std::function<double(const VectorRd&)>
rot_linear = [](const VectorRd & x) -> double {
    return 3.0; 
};
 
//------------------------------------------------------------------------------
  
  static std::function<VectorRd(const VectorRd&)>
  quadratic = [](const VectorRd & x) -> VectorRd {
             return linear(x) + (std::pow(x(0), 2) + 2. * std::pow(x(1), 2)) * VectorRd(-1., 0.);
             };
 
  static std::function<MatrixRd(const Eigen::Vector2d&)>
  grad_quadratic = [](const Eigen::Vector2d & x) -> MatrixRd {
        MatrixRd G = grad_linear(x);
        G.row(0) += -1 * VectorRd(2*x(0), 4*x(1)).transpose();      
            return G; 
              };
 
  static std::function<double(const VectorRd&)>
  rot_quadratic = [](const VectorRd & x) -> double {
    return 3.0 + 4.0 * 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