https://github.com/feelpp/feelpp
Tip revision: c4463293f684bf60433b75d88a9fde449eaf4599 authored by Vincent Huber on 29 August 2014, 14:09:41 UTC
Update CMakeLists.txt
Update CMakeLists.txt
Tip revision: c446329
test_raviartthomas.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=tcl: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: 2010-05-08
Copyright (C) 2010 Université Joseph Fourier (Grenoble I)
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 2.1 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_rt.cpp
\author Christophe Prud'homme <christophe.prudhomme@feelpp.org>
\date 2010-05-08
*/
#define BOOST_TEST_MODULE Raviar-Thomas polynomials test
// Boost.Test
#define USE_TEST 1
#include <boost/test/unit_test.hpp>
using boost::unit_test::test_suite;
#include <feel/feelpoly/raviartthomas.hpp>
#include <polynomial_set.hpp>
BOOST_AUTO_TEST_SUITE( rt_testsuite )
BOOST_AUTO_TEST_CASE( ravth_2d_0 )
{
using namespace Feel;
typedef RaviartThomas<0>::apply<2>::type rt0_type;
rt0_type rt0;
typedef boost::numeric::ublas::matrix<double> matrix_type ;
typedef boost::numeric::ublas::vector<double> Points_type ;
typedef boost::numeric::ublas::vector<Points_type> vectors_type ;
matrix_type m( 6 ,3 ) ;
m( 0 ,0 )=m( 0 ,1 )=m( 1 ,0 )=m( 1 ,2 )=m( 2 ,1 )=m( 3 ,0 )=m( 3 ,2 )=m( 4 ,0 )=m( 4 ,1 )=m( 5 ,2 ) = 1./4.;
m( 0 ,2 )=m( 1 ,1 )=m( 2 ,2 )=m( 3 ,1 )=m( 4 ,2 )=m( 5 ,1 ) = 0. ;
m( 2 ,0 )=m( 5 ,0 )= -1./4. ;
Points_type p1( 2 );
p1( 0 )=-1./3.;
p1( 1 ) =-1./3.;
Points_type p2( 2 );
p2( 0 )=-1.;
p2( 1 ) =-1.;
Points_type p3( 2 );
p3( 0 )= 1.;
p3( 1 ) =-1.;
Points_type p4( 2 );
p4( 0 )=-1.;
p4( 1 ) = 1.;
Points_type p5( 2 );
p5( 0 )=-1.;
p5( 1 ) = 0.;
Points_type p6( 2 );
p6( 0 )=-1./3.;
p6( 1 ) =-1./3.;
Points_type p7( 2 );
p7( 0 )= 0.;
p7( 1 ) = 0.;
Points_type p8( 2 );
p8( 0 )= 0.;
p8( 1 ) =-1.;
Points_type p9( 2 );
p9( 0 )=-1./2.;
p9( 1 ) = 1./2.;
Points_type p10( 2 );
p10( 0 )= 1./2.;
p10( 1 ) =-1./2.;
vectors_type P( 10 ) ;
P( 0 )=p1;
P( 1 )=p2;
P( 2 )=p3;
P( 3 )=p4;
P( 4 )=p5;
P( 5 )=p6;
P( 6 )=p7;
P( 7 )=p8;
P( 8 )=p9;
P( 9 )=p10 ;
Polynomialset<10,6,3,10> Poly( m ) ;
BOOST_TEST_MESSAGE( "**************************************************Analytic method*********************************************" << "\n" );
BOOST_TEST_MESSAGE( "coefficients= " << m << "\n" );
BOOST_TEST_MESSAGE( "Points= " << P << "\n" );
BOOST_TEST_MESSAGE( "evaluate_at_Points= " <<Poly.evaluates_Points ( P ) << "\n" );
BOOST_TEST_MESSAGE( "*****************************************Feel method***********************************************************" << "\n" );
rt0_type::points_type pts( 2,10 );
pts( 0 ,0 ) = -1./3. ;
pts( 1 ,0 ) = -1./3. ;
pts( 0 ,1 ) = -1. ;
pts( 1 ,1 ) = -1. ;
pts( 0 ,2 ) = 1. ;
pts( 1 ,2 ) = -1. ;
pts( 0 ,3 ) = -1. ;
pts( 1 ,3 ) = 1. ;
pts( 0 ,4 ) = -1. ;
pts( 1 ,4 ) = 0. ;
pts( 0 ,5 )=-1./3.;
pts( 1 ,5 ) =-1./3.;
pts( 0 ,6 )= 0.;
pts( 1 ,6 ) = 0.;
pts( 0 ,7 )= 0.;
pts( 1 ,7 ) =-1.;
pts( 0 ,8 )=-1./2.;
pts( 1 ,8 ) = 1./2.;
pts( 0 ,9 )= 1./2.;
pts( 1 ,9 ) =-1./2.;
BOOST_TEST_MESSAGE( "pts= " << pts << "\n" );
auto eval_at_pts = rt0.evaluate( pts );
BOOST_TEST_MESSAGE( "eval_at_pts= "<< eval_at_pts << "\n" );
matrix_type diff( Poly.evaluates_Points( P )-eval_at_pts );
double error = ublas::norm_frobenius( diff );
BOOST_CHECK_SMALL( error, 1e-10 );
}
BOOST_AUTO_TEST_CASE( ravth_3d_0 )
{
using namespace Feel;
typedef RaviartThomas<0>::apply<3>::type rt0_type;
typedef GeoND<3,Simplex<3, 1, 3> >::point_type point_type;
typedef boost::numeric::ublas::matrix<double> matrix_type ;
typedef boost::numeric::ublas::vector<double> Points_type ;
typedef boost::numeric::ublas::vector<Points_type> vectors_type ;
GeoND<3,Simplex<3, 1, 3> > tetra;
point_type V1;
V1( 0 )=-1;
V1( 1 )=-1;
V1( 2 )=-1;
point_type V2;
V2( 0 )=1;
V2( 1 )=-1;
V2( 2 )=-1;
point_type V3;
V3( 0 )=-1;
V3( 1 )=1;
V3( 2 )=-1;
point_type V4;
V4( 0 )=-1;
V4( 1 )=-1;
V4( 2 )=1;
tetra.setPoint( 0, V1 );
tetra.setPoint( 1, V2 );
tetra.setPoint( 2, V3 );
tetra.setPoint( 3, V4 );
tetra.update();
BOOST_TEST_MESSAGE( "[tetra] barycenter = " << tetra.barycenter() << "\n" );
BOOST_TEST_MESSAGE( "[tetra] measure = " << tetra.measure()<< "\n" );
BOOST_TEST_MESSAGE( tetra.faceMeasure( 0 ) << "\n" );
BOOST_TEST_MESSAGE( tetra.faceMeasure( 1 ) << "\n" );
BOOST_TEST_MESSAGE( tetra.faceMeasure( 2 ) << "\n" );
BOOST_TEST_MESSAGE( tetra.faceMeasure( 3 ) << "\n" );
rt0_type rt0;
matrix_type m( 12 ,4 ) ;
m( 0 ,0 )=m( 0 ,1 )=m( 1 ,0 )=m( 1 ,2 )=m( 2 ,0 )=m( 2 ,3 )=m( 3 ,1 )=m( 4 ,0 )
=m( 4 ,2 )=m( 5 ,0 )=m( 5 ,3 )=m( 6,0 ) = m( 6,1 ) =m( 7 ,2 ) =m( 8 , 0 )
=m( 8, 3 )=m( 9,0 ) = m( 9 ,1 ) =m( 10,0 ) =m( 10 ,2 )
=m( 11,3 ) = 1./( 3*tetra.measure() ) ;
m( 0 ,2 )= m( 0 ,3 )= m( 1 ,1 )= m( 1 ,3 )= m( 2 ,1 )=m( 2 ,2 )
=m( 3 ,2 )= m( 3 ,3 )= m( 4 ,1 )= m( 4 ,3 )=m( 5 ,1 )= m( 5 ,2 )
= m( 6 ,2 )= m( 6 ,3 ) = m( 7 ,1 )= m( 7 ,3 )= m( 8 ,1 )= m( 8 ,2 )=
m( 9 ,2 )= m( 9 ,3 )= m( 10 ,1 )= m( 10 ,3 )= m( 11 ,1 )= m( 11 ,2 )= 0. ;
m( 3 ,0 )=m( 7 ,0 )=m( 11 , 0 ) = -1./( 3*tetra.measure() ) ;
Points_type p1( 3 );
p1( 0 )=-1.;
p1( 1 ) =-1.;
p1( 2 ) = -1. ;
Points_type p2( 3 );
p2( 0 )= 1.;
p2( 1 ) =-1.;
p2( 2 ) = -1. ;
Points_type p3( 3 );
p3( 0 )= -1.;
p3( 1 ) =1.;
p3( 2 ) = -1. ;
Points_type p4( 3 );
p4( 0 )=-1.;
p4( 1 ) =-1.;
p4( 2 ) = 1. ;
Points_type p5( 3 );
p5( 0 )=-1.;
p5( 1 ) = 0.;
p5( 2 ) = -1. ;
Points_type p6( 3 );
p6( 0 )=-1.;
p6( 1 ) =-1.;
p6( 2 ) = 0.;
Points_type p7( 3 );
p7( 0 )= -1./3.;
p7( 1 ) =-1./3.;
p7( 2 ) =-1./3.;
Points_type p8( 3 );
p8( 0 )= -1.;
p8( 1 ) =-1./3.;
p8( 2 ) =-1./3.;
Points_type p9( 3 );
p9( 0 )=-1./3.;
p9( 1 ) = -1.;
p9( 2 ) = -1./3.;
Points_type p10( 3 );
p10( 0 )= -1./3.;
p10( 1 ) =-1./3.;
p10( 2 ) =-1.;
vectors_type P( 10 ) ;
P( 0 )=p1;
P( 1 )=p2;
P( 2 )=p3;
P( 3 )=p4;
P( 4 )=p5;
P( 5 )=p6;
P( 6 )=p7;
P( 7 )=p8;
P( 8 )=p9;
P( 9 )=p10 ;
Polynomialset<20,12,4,10> Poly( m ) ;
BOOST_TEST_MESSAGE( "**************************************************Analytic method*********************************************" << "\n" );
BOOST_TEST_MESSAGE( "coefficients= " << m << "\n" );
BOOST_TEST_MESSAGE( "\n" );
BOOST_TEST_MESSAGE( "Points= " << P << "\n" );
BOOST_TEST_MESSAGE( "\n" );
BOOST_TEST_MESSAGE( "evaluate_at_Points= " <<Poly.evaluates_Points ( P ) << "\n" );
BOOST_TEST_MESSAGE( "\n" );
BOOST_TEST_MESSAGE( "*****************************************Feel method***********************************************************" << "\n" );
rt0_type::points_type pts( 3,10 );
pts( 0 ,0 ) = -1. ;
pts( 1 ,0 ) = -1. ;
pts( 2 ,0 ) = -1. ;
pts( 0 ,1 ) = 1. ;
pts( 1 ,1 ) = -1. ;
pts( 2 ,1 ) = -1. ;
pts( 0 ,2 ) = -1. ;
pts( 1 ,2 ) = 1. ;
pts( 2 ,2 ) =-1. ;
pts( 0 ,3 ) = -1. ;
pts( 1 ,3 ) = -1. ;
pts( 2 ,3 ) = 1. ;
pts( 0 ,4 ) = -1.;
pts( 1 ,4 ) = 0. ;
pts( 2 ,4 ) = -1. ;
pts( 0 ,5 ) = -1.;
pts( 1 ,5 ) = -1.;
pts( 2 ,5 ) = 0. ;
pts( 0 ,6 )=-1./3.;
pts( 1 ,6 ) = -1./3.;
pts( 2 ,6 ) = -1./3. ;
pts( 0 ,7 )= -1.;
pts( 1 ,7 ) = -1./3.;
pts( 2 ,7 ) = -1./3. ;
pts( 0 ,8 )= -1./3.;
pts( 1 ,8 ) = -1.;
pts( 2 ,8 ) = -1./3. ;
pts( 0 ,9 )= -1./3.;
pts( 1 ,9 ) =-1./3.;
pts( 2 ,9 ) = -1. ;
BOOST_TEST_MESSAGE( "pts= " << pts << "\n" );
BOOST_TEST_MESSAGE( "\n" );
auto eval_at_pts = rt0.evaluate( pts );
BOOST_TEST_MESSAGE( "eval_at_pts= "<< eval_at_pts << "\n" );
double error = ublas::norm_frobenius( Poly.evaluates_Points( P )-eval_at_pts );
BOOST_CHECK_SMALL ( error, 1e-14 ) ;
}
BOOST_AUTO_TEST_CASE( ravth_2d_1 )
{
Feel::Assert::setLog( "test_ravth.log" );
using namespace Feel;
BOOST_TEST_MESSAGE( "start" );
typedef RaviartThomas<1>::apply<2>::type rt0_type;
typedef boost::numeric::ublas::matrix<double> matrix_type ;
typedef boost::numeric::ublas::vector<double> Points_type ;
typedef boost::numeric::ublas::vector<Points_type> vectors_type ;
rt0_type rt0;
BOOST_TEST_MESSAGE( "RT_1 instantiated" );
matrix_type m( 16 ,6 ) ;
/* for(uint k = 0 ; k < 12 ;++ k )
{
m(k , 3) = m(k , 4) = m(k ,5) =0 ;
}
*/
m( 0 ,0 )=m( 0 ,1 )=m( 1 ,0 )=m( 1 ,2 )=m( 2 ,1 )=m( 3 ,0 )=m( 3 ,2 )
=m( 4 ,0 )=m( 4 ,1 )=m( 5 ,2 ) = 1./4. ;
m( 2 ,0 )=m( 5 ,0 )= -1./4. ;
m( 6 ,0 )=m( 6 ,1 )= m( 13 ,0 )=m( 14 ,0 ) = -3./4. ;
m( 7 ,0 )=m( 7 ,2 )= 3./4. ;
m( 8 ,0 )=m( 8 ,1 )= m( 10 ,0 )=m( 10 ,1 )= m( 12 ,0 )=m( 15 ,0 ) = -6./4. ;
m( 9 ,0 )=m( 9 ,2 )= m( 11 ,0 )=m( 11 ,2 )= 6./4. ;
m( 11 ,1 )= 12./4. ;
m( 12 ,1 )=m( 15 ,2 )= -8./4. ;
m( 12 ,2 )=m( 12 ,3 )= m( 15 ,3 )= -4./4. ;
m( 12 ,4 )=m( 13 ,5 )= m( 14 ,4 )=m( 15 ,5 )= -2./4. ;
m( 13 ,2 ) = -5./4. ;
for ( uint k = 0 ; k < 12 ; ++ k )
{
m( k , 3 ) = m( k , 4 ) = m( k ,5 ) =0 ;
}
BOOST_TEST_MESSAGE( "coeff matrix created " );
// BOOST_TEST_MESSAGE( " coefficients= " << m << "\n" );
/*
m(0 ,0)=m(0 ,1)=m(1 ,0)=m(1 ,2)=m(2 ,1)=m(3 ,0)=m(3 ,2)=m(4 ,0)
=m(4 ,1)=m(5 ,2) = 1./4.;
m(0 ,2)=m(1 ,1)=m(2 ,2)=m(3 ,1)=m(4 ,2)=m(5 ,1) = 0. ;
m(2 ,0)=m(5 ,0)= -1./4. ;
*/
Points_type p1( 2 );
p1( 0 )=-1./3.;
p1( 1 ) =-1./3.;
Points_type p2( 2 );
p2( 0 )=-1.;
p2( 1 ) =-1.;
Points_type p3( 2 );
p3( 0 )= 1.;
p3( 1 ) =-1.;
Points_type p4( 2 );
p4( 0 )=-1.;
p4( 1 ) = 1.;
Points_type p5( 2 );
p5( 0 )=-1.;
p5( 1 ) = 0.;
Points_type p6( 2 );
p6( 0 )=-1./3.;
p6( 1 ) =-1./3.;
Points_type p7( 2 );
p7( 0 )= 0.;
p7( 1 ) = 0.;
Points_type p8( 2 );
p8( 0 )= 0.;
p8( 1 ) =-1.;
Points_type p9( 2 );
p9( 0 )=-1./2.;
p9( 1 ) = 1./2.;
Points_type p10( 2 );
p10( 0 )= 1./2.;
p10( 1 ) =-1./2.;
vectors_type P( 10 ) ;
P( 0 )=p1;
P( 1 )=p2;
P( 2 )=p3;
P( 3 )=p4;
P( 4 )=p5;
P( 5 )=p6;
P( 6 )=p7;
P( 7 )=p8;
P( 8 )=p9;
P( 9 )=p10 ;
BOOST_TEST_MESSAGE( "points created " );
Polynomialset<10,16,6,10> Poly( m ) ;
// Points_type p(2); p(0)= 0.; p(1) = 0.;
// vectors_type P(1 ) ;
// P(0)=p ;
BOOST_TEST_MESSAGE( "polyset created " );
BOOST_TEST_MESSAGE( "**************************************************Analytic method*********************************************" << "\n" );
BOOST_TEST_MESSAGE( "coefficients= " << m << "\n" );
BOOST_TEST_MESSAGE( "Points= " << P << "\n" );
BOOST_TEST_MESSAGE( "evaluate_at_Points= " <<Poly.evaluates_Points ( P ) << "\n" );
BOOST_TEST_MESSAGE( "*****************************************Feel method***********************************************************" << "\n" );
rt0_type::points_type pts( 2,10 );
// pts(0 ,0) = 0. ; pts(1 ,0) = 0. ;
pts( 0 ,0 ) = -1./3. ;
pts( 1 ,0 ) = -1./3. ;
pts( 0 ,1 ) = -1. ;
pts( 1 ,1 ) = -1. ;
pts( 0 ,2 ) = 1. ;
pts( 1 ,2 ) = -1. ;
pts( 0 ,3 ) = -1. ;
pts( 1 ,3 ) = 1. ;
pts( 0 ,4 ) = -1. ;
pts( 1 ,4 ) = 0. ;
pts( 0 ,5 )=-1./3.;
pts( 1 ,5 ) =-1./3.;
pts( 0 ,6 )= 0.;
pts( 1 ,6 ) = 0.;
pts( 0 ,7 )= 0.;
pts( 1 ,7 ) =-1.;
pts( 0 ,8 )=-1./2.;
pts( 1 ,8 ) = 1./2.;
pts( 0 ,9 )= 1./2.;
pts( 1 ,9 ) =-1./2.;
BOOST_TEST_MESSAGE( "pts= " << pts << "\n" );
auto eval_at_pts = rt0.evaluate( pts );
BOOST_TEST_MESSAGE( "eval_at_pts= "<< eval_at_pts << "\n" );
double error = ublas::norm_frobenius( Poly.evaluates_Points( P )-eval_at_pts );
BOOST_CHECK_SMALL( error, 1e-10 );
}
BOOST_AUTO_TEST_SUITE_END()