Revision af03cba596d2a9ea242804f26c247846d1ee0e34 authored by Vincent Chabannes on 09 February 2021, 14:07:24 UTC, committed by Vincent Chabannes on 09 February 2021, 14:07:24 UTC
1 parent 18b1aa2
test_twospacesmap.cpp
#define BOOST_TEST_MODULE test_twospaces_map
#include <feel/feelcore/testsuite.hpp>
#include <feel/feelfilters/unitcircle.hpp>
#include <feel/feelvf/vf.hpp>
#include <feel/feeldiscr/pchv.hpp>
#include <feel/feeldiscr/thch.hpp>
#include <feel/feelvf/form.hpp>
#include <feel/feelfilters/savegmshmesh.hpp>
#include <feel/feelfilters/loadmesh.hpp>
#include <feel/feeldiscr/twospacesmap.hpp>
#include <feel/feelfilters/exporter.hpp>
using namespace Feel;
FEELPP_ENVIRONMENT_NO_OPTIONS
BOOST_AUTO_TEST_SUITE( test_twospacesmap )
BOOST_AUTO_TEST_CASE( onespace )
{
typedef Mesh<Simplex<2,1>> mesh_type;
typedef Pchv_type<mesh_type,2> space_type;
auto mesh = unitCircle();
auto Xh = Pchv<2>( mesh );
auto u = Xh->element();
u.on( _range=elements(mesh), _expr=20*abs(sin(2*Pi*Px())*exp(Px()*Py()+Py()*Py()))*oneX()
+ 20*abs(cos(2*Pi*Px())*exp(Px()*Py()+Py()*Py()))*oneY() );
int size = u.size();
std::vector<int> randoms;
std::set<int> element_ids;
std::srand (std::time(NULL));
if ( Environment::isMasterRank() )
for ( int i=0; i<10; i++ )
randoms.push_back( std::rand() % size );
boost::mpi::broadcast( Environment::worldComm(), randoms, Environment::worldComm().masterRank() );
for ( auto index : randoms )
{
auto searchGpDof = Xh->dof()->searchGlobalProcessDof( index );
if ( boost::get<0>( searchGpDof ) )
{
size_type gpdof = boost::get<1>( searchGpDof );
for ( auto const& dof : Xh->dof()->globalDof( gpdof ) )
{
size_type eltId = dof.second.elementId();
if ( Xh->mesh()->element( eltId ).isGhostCell() )
continue;
element_ids.insert( eltId );
}
}
}
auto newmesh = createSubmesh( _mesh=mesh, _range=idelements(mesh,element_ids.begin(), element_ids.end()) );
saveGMSHMesh(_mesh=newmesh, _filename="submesh.msh" );
Environment::worldComm().barrier();
auto seqmesh = loadMesh( _mesh=new mesh_type,
_filename="submesh.msh", _worldcomm= Environment::worldCommSeqPtr() );
auto Rh = space_type::New(_mesh=seqmesh );
auto my_map = TwoSpacesMap<space_type>( Rh, Xh );
auto ur = Rh->element();
ur.on( _range=elements(seqmesh), _expr=20*abs(sin(2*Pi*Px())*exp(Px()*Py()+Py()*Py()))*oneX()
+ 20*abs(cos(2*Pi*Px())*exp(Px()*Py()+Py()*Py())
)*oneY() );
for ( auto index : randoms )
{
double value = 0;
int proc_number = Xh->map().procOnGlobalCluster(index);
if ( Environment::worldComm().globalRank()==proc_number )
{
auto searchGpDof = Xh->dof()->searchGlobalProcessDof( index );
CHECK( boost::get<0>( searchGpDof ) ) << "GPDof not found\n";
size_type gpdof = boost::get<1>( searchGpDof );
value = u( gpdof );
}
boost::mpi::broadcast( Environment::worldComm(), value, proc_number );
int index_r = my_map.clusterToSequential( index );
double value_r = ur( index_r );
BOOST_CHECK_CLOSE( value, value_r, 1e-9 );
}
my_map.project( ur, u );
double diff = normL2( elements(newmesh), idv(u)-idv(ur) );
BOOST_CHECK_SMALL( diff, 1e-9 );
}
BOOST_AUTO_TEST_CASE( composite )
{
typedef Mesh<Simplex<2,1>> mesh_type;
typedef THch_type<1,mesh_type> space_type;
auto mesh = unitCircle();
auto Xh = space_type::New( mesh );
auto Xh0 = Xh->functionSpace<0>();
auto Xh1 = Xh->functionSpace<1>();
auto U = Xh->element();
auto u = U.template element<0>();
auto p = U.template element<1>();
auto expr_u = 20*sin(2*Pi*Px())*exp(Px()*Py()+Py()*Py())*oneX() + 20*cos(2*Pi*Px())*exp(Px()*Py()+Py()*Py())*oneY();
auto expr_p = cos(Pi*Px()*pow(Py(),3))*exp( Px()-Py()*(Px() - Py() ) );
u.on( _range=elements(mesh), _expr=expr_u );
p.on( _range=elements(mesh), _expr=expr_p );
int size = U.size();
std::vector<int> randoms;
std::set<int> element_ids;
// build a random set of index in the range of the vector
std::srand (std::time(NULL));
if ( Environment::isMasterRank() )
for ( int i=0; i<10; i++ )
randoms.push_back( std::rand() % size );
boost::mpi::broadcast( Environment::worldComm(), randoms, Environment::worldComm().masterRank() );
for ( auto index : randoms )
{
auto searchGpDof = Xh->dof()->searchGlobalProcessDof( index );
if ( boost::get<0>( searchGpDof ) )
{
size_type gpdof = boost::get<1>( searchGpDof );
int space = Xh->dof()->databaseIndexFromContainerId( gpdof );
gpdof = Xh->dof()->containerIdToDofId( space, gpdof );
if ( space==0 )
{
for ( auto const& dof : Xh0->dof()->globalDof( gpdof ) )
{
size_type eltId = dof.second.elementId();
if ( Xh0->mesh()->element( eltId ).isGhostCell() )
continue;
element_ids.insert( eltId );
}
}
else if ( space==1 )
{
for ( auto const& dof : Xh1->dof()->globalDof( gpdof ) )
{
size_type eltId = dof.second.elementId();
if ( Xh1->mesh()->element( eltId ).isGhostCell() )
continue;
element_ids.insert( eltId );
}
}
}
}
auto newmesh = createSubmesh( _mesh=mesh, _range=idelements(mesh,element_ids.begin(), element_ids.end()) );
saveGMSHMesh(_mesh=newmesh, _filename="submesh.msh" );
Environment::worldComm().barrier();
auto seqmesh = loadMesh( _mesh=new mesh_type,
_filename="submesh.msh", _worldcomm= Environment::worldCommSeqPtr() );
auto Rh = space_type::New(_mesh=seqmesh );
auto my_map = TwoSpacesMap<space_type>( Rh, Xh );
auto Ur = Rh->element();
auto ur = Ur.template element<0>();
auto pr = Ur.template element<1>();
ur.on( _range=elements(seqmesh), _expr=expr_u );
pr.on( _range=elements(seqmesh), _expr=expr_p );
for ( auto index : randoms )
{
double value = 0;
int proc_number = Xh->map().procOnGlobalCluster(index);
if ( Environment::worldComm().globalRank()==proc_number )
{
auto searchGpDof = Xh->dof()->searchGlobalProcessDof( index );
CHECK( boost::get<0>( searchGpDof ) ) << "GPDof not found\n";
size_type gpdof = boost::get<1>( searchGpDof );
value = U( gpdof );
}
boost::mpi::broadcast( Environment::worldComm(), value, proc_number );
int index_r = my_map.clusterToSequential( index );
CHECK( index_r>=0 )<<"index "<< index <<" was not found in p_to_s map\n";
double value_r = Ur( index_r );
BOOST_CHECK_CLOSE( value, value_r, 1e-9 );
}
my_map.project( Ur, U );
double diff = normL2( elements(newmesh), idv(u)-idv(ur) );
BOOST_CHECK_SMALL( diff, 1e-9 );
diff = normL2( elements(newmesh), idv(p)-idv(pr) );
BOOST_CHECK_SMALL( diff, 1e-9 );
auto Rh0 = Rh->functionSpace<0>();
auto Rh1 = Rh->functionSpace<1>();
auto v = Xh0->element();
auto q = Xh1->element();
auto vr = Rh0->element();
auto qr = Rh1->element();
v.on( _range=elements(mesh), _expr=expr_u );
q.on( _range=elements(mesh), _expr=expr_p );
my_map.template project<0>( vr, v );
my_map.template project<1>( qr, q );
diff = normL2( elements(newmesh), idv(v)-idv(vr) );
BOOST_CHECK_SMALL( diff, 1e-9 );
diff = normL2( elements(newmesh), idv(q)-idv(qr) );
BOOST_CHECK_SMALL( diff, 1e-9 );
}
BOOST_AUTO_TEST_SUITE_END()
Computing file changes ...