https://github.com/feelpp/feelpp
Tip revision: 95faa71eb659d1a1ddd5ae72a3d6cd879bc7ba02 authored by vdoyeux on 26 July 2016, 15:48:57 UTC
start working on self label
start working on self label
Tip revision: 95faa71
test_im.cpp
/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=cpp:et:sw=4:ts=4:sts=4
This file is part of the Feel library
Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
Date: 2005-07-29
Copyright (C) 2005,2006 EPFL
Copyright (C) 2007, 2008 Université Joseph Fourier Grenoble 1
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 3.0 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
\file test_im.cpp
\author Christophe Prud'homme <christophe.prudhomme@feelpp.org>
\date 2005-07-29
*/
#define BOOST_TEST_MODULE integration methods test
// Boost.Test
#include <boost/test/unit_test.hpp>
using boost::unit_test::test_suite;
#include <feel/feelcore/feel.hpp>
#include <feel/feelpoly/im.hpp>
#include <feel/feelpoly/imsimplex.hpp>
using namespace Feel;
struct F
{
F() : i( 0 )
{
BOOST_TEST_MESSAGE( "setup fixture" );
}
~F()
{
BOOST_TEST_MESSAGE( "teardown fixture" );
}
int i;
};
template<typename T>
T two( typename Feel::node<T>::type const& /*t*/ )
{
return 2.0;
}
template<typename T>
T one( typename Feel::node<T>::type const& /*t*/ )
{
return 1.0;
}
template<int P, typename T>
T xp( typename Feel::node<T>::type const& t )
{
if ( P ==0 ) return 1.0;
T res = t[0];
for ( int i = 1; i < P; ++ i )
res *= t[0];
//std::cout << "res xp <" << P << ">: " << res << "\n";
return res;
//return pow( t, T(P) );
}
template<typename T>
struct P2N
{
typedef T value_type;
typedef typename Feel::node<T>::type node_type;
P2N( int alpha1, int alpha2 )
:
M_alpha1( alpha1 ),
M_alpha2( alpha2 )
{}
value_type operator()( node_type const& t )
{
value_type res = Feel::math::pow( t[0], M_alpha1 )* Feel::math::pow( t[1], M_alpha2 );
return res;
}
value_type integral() const
{
value_type a1 = Feel::math::pow( -1.0, M_alpha2+1 )/( M_alpha2 +1 );
value_type a2 = 1./( M_alpha1+M_alpha2+2 )-1./( M_alpha1+1 );
value_type a3 = Feel::math::pow( -1.0, M_alpha1+M_alpha2+2 )/( M_alpha1 + M_alpha2+2 );
value_type a4 = Feel::math::pow( -1.0, M_alpha1+1 )/( M_alpha1 + 1 );
return a1*( a2-( a3-a4 ) );
}
int M_alpha1;
int M_alpha2;
};
template<typename T>
struct P3N
{
typedef T value_type;
typedef typename Feel::node<T>::type node_type;
P3N( int alpha1, int alpha2 = 1, int alpha3 = 1 )
:
M_alpha1( alpha1 ),
M_alpha2( alpha2 ),
M_alpha3( alpha3 )
{}
value_type operator()( node_type const& t )
{
value_type res = Feel::math::pow( t[0], M_alpha1 )* Feel::math::pow( t[1], M_alpha2 )* Feel::math::pow( t[2], M_alpha3 );
return res;
}
value_type integral( int face )
{
value_type a = M_alpha1;
value_type b = M_alpha2;
value_type c = M_alpha3;
switch ( face )
{
default:
case FACE_0:
return 0;
case FACE_1:
return 0;
case FACE_2:
return 0;
case FACE_3:
return -( Feel::math::pow( -1.,a )+1 )*Feel::math::pow( -1.,c )/( ( a+1 )*( a+3 ) );
}
}
value_type integral() const
{
value_type a = M_alpha1;
return ( 2*Feel::math::pow( a,3. )*Feel::math::pow( ( -1. ),a )+21*Feel::math::pow( a,2. )*Feel::math::pow( ( -1. ),a )+
67*a*Feel::math::pow( -1.,a )+63*Feel::math::pow( -1.,a )+3*Feel::math::pow( a,2. )+21*a+33 )/( 3*( a+1 )*( a+2 )*( a+3 )*( a+4 )*( a+5 ) );
}
int M_alpha1;
int M_alpha2;
int M_alpha3;
};
template<int P, typename T>
T xpdim( typename Feel::node<T>::type const& t )
{
if ( P ==0 ) return 1.0;
T res = t[0];
for ( int i = 1; i < P; ++ i )
res *= t[i];
//std::cout << "res xp <" << P << ">: " << res << "\n";
return res;
//return pow( t, T(P) );
}
template<typename T>
T sint( typename Feel::node<T>::type const& t )
{
return sin( t[0] );
}
template<typename T>
T cost( typename Feel::node<T>::type const& t )
{
return cos( t[0] );
}
template<typename T>
T xp1( typename Feel::node<T>::type const& t )
{
// return pow(t[0]*t[1],T(P));
return t[0];
}
template<typename T>
T xp2( typename Feel::node<T>::type const& t )
{
// return pow(t[0]*t[1],T(P));
return t[0]*t[0];
}
template<typename T>
T xp3( typename Feel::node<T>::type const& t )
{
// return pow(t[0]*t[1],T(P));
return t[0]*t[0]*t[0];
}
template<typename T>
T xp4( typename Feel::node<T>::type const& t )
{
// return pow(t[0]*t[1],T(P));
return t[0]*t[0]*t[0]*t[0];
}
template<int P, typename T>
T exp2( typename Feel::node<T>::type const& t )
{
// return pow(t[0]*t[1],T(P));
return exp( t[0]+t[1] );
}
template<typename T>
T sint2( typename Feel::node<T>::type const& t )
{
return sin( t[0] )*sin( t[1] );
}
template<typename T>
T coscos( typename Feel::node<T>::type const& t )
{
return cos( t[0] )*cos( t[1] );
}
template<Feel::uint16_type D,
Feel::uint16_type N,
typename T,
template<class Convex, Feel::uint16_type O, typename T2> class QPS = Feel::Gauss>
class TestImPK
{
public:
typedef T value_type;
typedef typename Feel::node<value_type>::type node_type;
typedef typename mpl::if_<mpl::equal_to<mpl::int_<D>,mpl::int_<1> >,
mpl::identity<node_type>,
mpl::identity<node_type> >::type::type arg_type;
//IM<D, N, T, Simplex, QPS> im;
template<int IMORDER> struct MyIm
:
public IM<D,N,T,Simplex>
{};
TestImPK( value_type res,
boost::function<value_type( arg_type const& )> const& func,
value_type tol = Feel::type_traits<T>::epsilon() )
:
M_res( res ),
M_tol( tol ),
M_face( -2 ),
M_func( func )
{
}
TestImPK( int __face,
value_type res,
boost::function<value_type( arg_type const& )> const& func,
value_type tol = Feel::type_traits<T>::epsilon() )
:
M_res( res ),
M_tol( tol ),
M_face( __face ),
M_func( func )
{
}
value_type operator()() const
{
MyIm<N> im;
value_type res = 0.0;
value_type integral = 0.0;
if ( M_face == -2 )
{
integral = im.integrateAtPoints( M_func );
}
if ( M_face != -2 )
{
integral = im.integrateAtPoints( Feel::IntegrationFaceEnum( M_face ), M_func );
}
res = math::abs( integral - M_res );
if ( ! ( res < M_tol ) )
{
std::cout << "IM: " << im << "\n";
std::cout << "TestIm<"<<D << "," << typeid( T ).name() << ">:: integral = " << integral << "\n";
std::cout << "TestIm<"<<D << "," << typeid( T ).name()<< ">:: exact value = " << M_res << "\n";
std::cout << "TestIm<"<<D << "," << typeid( T ).name()<< ">:: error = " << res << "\n";
std::cout << "TestIm<"<<D << "," << typeid( T ).name()<< ">:: tolerance = " << M_tol << "\n";
}
BOOST_CHECK( res < M_tol );
return res;
}
value_type M_res;
value_type M_tol;
int M_face;
boost::function<value_type( arg_type const& )> M_func;
};
#define TESTS_XP( T ) \
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, T>( 0.0 , xp<1,T> ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,10,T>( T(2.0)/11.0 , xp<10,T> ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,2, T>( T(2.0)/3.0 , xp<2,T> ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,3, T>( 0.0 , xp<3,T> ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,20,T>( T(2.0)/21.0 , xp<20,T>) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,21,T>( 0.0 , xp<21,T>) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<1,15,T>( 0.0 , xp<21,T>) ) ) ); \
/**/
#define TESTS_SIN( T ) \
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, T>( 0.0 , sint<T> ) ) ) );\
test->add( BOOST_TEST_CASE( ( TestImPK<1,2, T>( 0.0 , sint<T> ) ) ) );\
test->add( BOOST_TEST_CASE( ( TestImPK<1,3, T>( 0.0 , sint<T> ) ) ) );\
test->add( BOOST_TEST_CASE( ( TestImPK<1,20,T>( 0.0 , sint<T> ) ) ) );\
/**/
#define TESTS_XP2( T ) \
test->add( BOOST_TEST_CASE( ( TestImPK<2,5, T>( -T(2.0)/3.0 , xp2<1,T> ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,17, T>( 1.0+T(1.0)/(exp(T(1.0))*exp(T(1.0))) , exp2<1,T>, 1e-15 ) ) ) ); \
/**/
#define TESTS_SIN2( T ) \
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, T>( 0.0 , sint2<T>, 1e-1 ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,2, T>( 0.0 , sint2<T>, 1.5e-2 ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,3, T>( 0.0 , sint2<T>, 5e-4 ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,4, T>( 0.0 , sint2<T>, 5e-4 ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,6, T>( 0.0 , sint2<T>, 5e-4 ) ) ) ); \
test->add( BOOST_TEST_CASE( ( TestImPK<2,20,T>( 0.0 , sint2<T>, 1.5e-16 ) ) ) ); \
/**/
#define TESTS_COSCOS( T ) \
test->add( BOOST_TEST_CASE( ( TestImPK<2,40,T>( 2.0*sin(T(1))*sin(T(1)) , coscos<T> ) ) ) );\
/**/
template<int D, int N, typename T>
class TestImQK
{
public:
typedef T value_type;
typedef typename Feel::node<T>::type node_type;
TestImQK( value_type res,
boost::function<value_type( node_type const& )> const& func,
value_type tol = Feel::type_traits<T>::epsilon() )
:
M_res( res ),
M_tol( tol ),
M_func( func )
{}
value_type operator()() const
{
using namespace Feel;
IM<D, N, T, Hypercube> im;
value_type res = math::abs( im.integrateAtPoints( M_func ) - M_res );
if ( !( res < M_tol ) )
{
std::cout << "TestImQK:: int = " << im.integrateAtPoints( M_func ) << "\n";
std::cout << "TestImQK:: exact = " << M_res << "\n";
std::cout << "TestImQK:: res = " << res << "\n";
std::cout << "TestImQK:: tol = " << M_tol << "\n";
}
BOOST_CHECK( res < M_tol );
return res;
}
value_type M_res;
value_type M_tol;
boost::function<value_type( node_type const& )> M_func;
};
#define TESTS_QK_COSCOS( T ) \
test->add( BOOST_TEST_CASE( ( TestImQK2D<40, T>( 4.0*sin(T(1.0))*sin(T(1.0)) , coscos<T>, 1.0E-14 ) ) ) );\
/**/
#if 1
// automatically registered test cases could be organized in test suites
BOOST_FIXTURE_TEST_SUITE( im1d_double_suite, F )
BOOST_AUTO_TEST_CASE( im1d_test1 )
{
TestImPK<1,1, double> t1( 2.0, one<double> );
t1();
}
BOOST_AUTO_TEST_CASE( im1d_test2 )
{
TestImPK<1,1, double> t2( 2.0*sin( double( 1.0 ) ) , cost<double>, 1.1E-2 );
t2();
}
BOOST_AUTO_TEST_CASE( im1d_test3 )
{
TestImPK<1,2, double> t3( 2.0*sin( double( 1.0 ) ) , cost<double>, 1.1E-2 );
t3();
}
BOOST_AUTO_TEST_CASE( im1d_test4 )
{
TestImPK<1,3, double> t4( 2.0*sin( double( 1.0 ) ) , cost<double>, 1.0E-4 );
t4();
}
BOOST_AUTO_TEST_CASE( im1d_test5 )
{
TestImPK<1,50,double> t5( 2.0*sin( double( 1.0 ) ) , cost<double> );
t5();
}
BOOST_AUTO_TEST_CASE( im1d_test6 )
{
TestImPK<1,1, double> t6( Feel::FACE_0, -1.0 , xp<1,double> );
t6();
}
BOOST_AUTO_TEST_CASE( im1d_test7 )
{
TestImPK<1,1, double> t7( Feel::FACE_1, 1.0 , xp<1,double> );
t7();
}
BOOST_AUTO_TEST_CASE( im1d_face_test1 )
{
TestImPK<1,1, double> t0( Feel::FACE_0, -1.0 , xp<1,double> );
t0();
TestImPK<1,1, double> t1( Feel::FACE_1, 1.0 , xp<1,double> );
t1();
}
BOOST_AUTO_TEST_CASE( im1d_face_test2 )
{
TestImPK<1,2, double> t0( Feel::FACE_0, -1.0 , xp<1,double> );
t0();
TestImPK<1,2, double> t1( Feel::FACE_1, 1.0 , xp<1,double> );
t1();
}
BOOST_AUTO_TEST_CASE( im1d_face_test3 )
{
TestImPK<1,3, double> t0( Feel::FACE_0, -1.0 , xp<1,double> );
t0();
TestImPK<1,3, double> t1( Feel::FACE_1, 1.0 , xp<1,double> );
t1();
}
BOOST_AUTO_TEST_CASE( im2d_test0 )
{
TestImPK<2,1, double> t1( P2N<double>( 0, 0 ).integral(), P2N<double>( 0, 0 ) );
t1();
}
BOOST_AUTO_TEST_CASE( im2d_test1 )
{
const int N = 1;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_test2 )
{
const int N = 2;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_test3 )
{
const int N = 3;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_test4 )
{
const int N = 4;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_test5 )
{
const int N = 5;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_test6 )
{
const int N = 6;
for ( int alpha1 = 0; alpha1 < N; ++ alpha1 )
{
int alpha2 = N-alpha1;
TestImPK<2,N, double> t10( P2N<double>( alpha1, alpha2 ).integral(), P2N<double>( alpha1, alpha2 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im2d_face_test1 )
{
TestImPK<2,1, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,1, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,1, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im2d_face_test2 )
{
TestImPK<2,2, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,2, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,2, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im2d_face_test3 )
{
TestImPK<2,3, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,3, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,3, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im2d_face_test4 )
{
TestImPK<2,4, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,4, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,4, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im2d_face_test5 )
{
TestImPK<2,5, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,5, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,5, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im2d_face_test6 )
{
TestImPK<2,6, double> t0( Feel::FACE_1, 2.0 , one<double> );
t0();
TestImPK<2,6, double> t1( Feel::FACE_0, Feel::math::sqrt( double( 8.0 ) ) , one<double> );
t1();
TestImPK<2,6, double> t2( Feel::FACE_2, 2.0 , one<double> );
t2();
}
BOOST_AUTO_TEST_CASE( im3d_test1 )
{
const int N = 2;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test2 )
{
const int N = 3;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test3 )
{
const int N = 4;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test5 )
{
const int N = 5;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test6 )
{
const int N = 6;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test7 )
{
const int N = 7;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test11 )
{
const int N = 11;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test13 )
{
const int N = 13;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test15 )
{
const int N = 15;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test17 )
{
const int N = 17;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test19 )
{
const int N = 19;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_test20 )
{
const int N = 20;
for ( int alpha1 = 0; alpha1 <= N-2; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = 1;
TestImPK<3,N, double> t10( P3N<double>( alpha1, alpha2, alpha3 ).integral(), P3N<double>( alpha1, alpha2, alpha3 ), 5e-13 );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test1 )
{
TestImPK<3,1, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,1, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,1, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,1, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
const int N = 1;
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test2 )
{
const int N = 2;
TestImPK<3,N, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,N, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,N, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,N, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test3 )
{
const int N = 3;
TestImPK<3,N, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,N, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,N, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,N, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test4 )
{
const int N = 4;
TestImPK<3,N, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,N, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,N, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,N, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test5 )
{
const int N = 5;
TestImPK<3,N, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,N, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,N, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,N, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_CASE( im3d_face_test6 )
{
const int N = 6;
TestImPK<3,N, double> t1( Feel::FACE_2, 2.0 , one<double> );
t1();
TestImPK<3,N, double> t2( Feel::FACE_1, 2.0 , one<double> );
t2();
TestImPK<3,N, double> t3( Feel::FACE_0, 2.0*Feel::math::sqrt( double( 3.0 ) ), one<double> );
t3();
TestImPK<3,N, double> t4( Feel::FACE_3, 2.0 , one<double> );
t4();
for ( int alpha1 = 0; alpha1 < N-1; ++ alpha1 )
{
int alpha2 = 1;
int alpha3 = N-1-alpha1;
TestImPK<3,N, double> t10( Feel::FACE_3, P3N<double>( alpha1, alpha2, alpha3 ).integral( FACE_3 ), P3N<double>( alpha1, alpha2, alpha3 ) );
t10();
}
}
BOOST_AUTO_TEST_SUITE_END()
#else
template<typename T>
void add_tests( test_suite* test )
{
using namespace Feel;
typedef T double;
// 1D quadrature on segment [-1;1]
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, value_type>( 2.0, one<value_type> ) ) ) );
TESTS_XP( value_type );
TESTS_SIN( value_type );
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, value_type>( 2.0*sin( value_type( 1.0 ) ) , cost<value_type>, 1.1E-2 ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<1,2, value_type>( 2.0*sin( value_type( 1.0 ) ) , cost<value_type>, 1.1E-2 ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<1,3, value_type>( 2.0*sin( value_type( 1.0 ) ) , cost<value_type>, 1.0E-4 ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<1,50,value_type>( 2.0*sin( value_type( 1.0 ) ) , cost<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, value_type>( Feel::FACE_0, -1.0 , xp<1,value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<1,1, value_type>( Feel::FACE_1, 1.0 , xp<1,value_type> ) ) ) );
// 2D quadrature on triangle [-1;1]x[-1;-x]
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, value_type>( 2.0 , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, value_type>( Feel::FACE_1, 2.0 , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, value_type>( Feel::FACE_0, Feel::math::sqrt( value_type( 8.0 ) ) , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, value_type>( Feel::FACE_2, 2.0 , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<2,5, value_type>( -2.0/value_type( 3.0 ) , xp2<1,value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<2,50, value_type>( value_type( 1.0 )+1.0/( exp( value_type( 1.0 ) )*exp( value_type( 1.0 ) ) ) , exp2<1,value_type> ) ) ) );
// 2D Quadratures on faces of the reference triangle
test->add( BOOST_TEST_CASE( ( TestImPK<2,1, value_type>( Feel::ALL_FACES, value_type( 4.0 )+sqrt( value_type( 8.0 ) ), one<value_type> ) ) ) );
// 3D quadrature [-1;1]x-[-1;-x]x[-1;-x]
test->add( BOOST_TEST_CASE( ( TestImPK<3,1, value_type>( 8.0/value_type( 6.0 ) , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<3,1, value_type>( Feel::FACE_2, 2.0 , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<3,1, value_type>( Feel::FACE_1, 2.0 , one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<3,1, value_type>( Feel::FACE_0, 2.0*Feel::math::sqrt( value_type( 3.0 ) ), one<value_type> ) ) ) );
test->add( BOOST_TEST_CASE( ( TestImPK<3,1, value_type>( Feel::FACE_3, 2.0 , one<value_type> ) ) ) );
TESTS_XP2( value_type );
TESTS_SIN2( value_type );
TESTS_COSCOS( value_type );
// 2D quadrature
test->add( BOOST_TEST_CASE( ( TestImQK<2,50, value_type>( 4.0*sin( value_type( 1.0 ) )*sin( value_type( 1.0 ) ) , coscos<value_type>, 1.8e-15 ) ) ) );
// 3D quadrature
test->add( BOOST_TEST_CASE( ( TestImQK<3,1, value_type>( 0.0, xpdim<3,value_type>, 1.8e-15 ) ) ) );
// 4D quadrature
test->add( BOOST_TEST_CASE( ( TestImQK<4,1, value_type>( 0.0, xpdim<4,value_type>, 1.8e-15 ) ) ) );
}
#if 1
test_suite*
init_unit_test_suite( int /*argc*/, char** /*argv*/ )
{
test_suite* test = BOOST_TEST_SUITE( "Integration methods test suite" );
add_tests<double>( test );
#if defined( FEELPP_HAS_QD_QD_H )
unsigned int old_cw;
fpu_fix_start( &old_cw );
add_tests<dd_real>( test );
add_tests<qd_real>( test );
#endif /* FEELPP_HAS_QD_QD_H */
#if defined( FEELPP_HAS_MPFR )
std::cout << "[feelpoly::testsuite::test_im] test with arbitrary precision\n";
mp::mp_init( 200 );
add_tests<mp_type>( test );
#endif
return test;
}
#else
int main()
{
//unsigned int old_cw;
//fpu_fix_start(&old_cw);
TestImPK<1,1, double> im0( Feel::FACE_0, 1.0 , one<double> );
im0();
TestImPK<1,1, double> im1( Feel::FACE_1, 1.0 , one<double> );
im1();
TestImPK<1,1, double> im11( Feel::FACE_0,-1.0 , xp<1,double> );
im11();
TestImPK<1,1, double> im2( Feel::FACE_0, 1.0 , xp<2,double> );
im2();
TestImPK<1,1, double> im3( Feel::FACE_1, 1.0 , xp<2,double> );
im3();
#if 0
#if 0
//test_suite* test = BOOST_TEST_SUITE( "Integration methods test suite" );
//add_tests<double>( test );
TestImPK<2,1, double> im( 2.0 , one<double> );
im();
TestImPK<2,1, dd_real> imdd( 2.0 , one<dd_real> );
imdd();
#endif
TestImQK<2,1, double> im1( 0.0 , xp<2,dd_real> );
im1();
TestImQK<3,1, double> im2( 0.0 , xp<3,dd_real> );
im2();
Feel::Tesseract tess;
std::cout << "Tess.Dim = " << tess.nDim << "\n";
TestImQK<4,1, double> im3( 0.0 , xp<4,dd_real> );
im3();
#endif
}
#endif
#endif // 1