https://github.com/feelpp/feelpp
Tip revision: 756e99c230d0621fd98af36e5c6c90575c5e3eac authored by Christophe Prud'homme on 12 April 2016, 15:51:58 UTC
basic doc
basic doc
Tip revision: 756e99c
test_twolaplaciansdistributed.cpp
/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-*/
#if 1
#define BOOST_TEST_MODULE test_twolaplaciansdistributed
#include <testsuite/testsuite.hpp>
#else
#define USE_BOOST_TEST 0
#endif
#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/pch.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelvf/vf.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelfilters/geotool.hpp>
#include <feel/feeldiscr/operatorinterpolation.hpp>
#include <feel/feeldiscr/operatorlagrangep1.hpp>
namespace test_twolaplaciansdistributed
{
using namespace Feel;
using namespace Feel::vf;
/*_________________________________________________*
* Options
*_________________________________________________*/
inline
po::options_description
makeOptions()
{
po::options_description desc_options( "test_twolaplaciansdistributed options" );
desc_options.add_options()
( "hsize", po::value<double>()->default_value( 0.02 ), "mesh size" )
( "1d-hsize", po::value<double>()->default_value( 0.02 ), "mesh size 1d" )
;
return desc_options.add( Feel::feel_options() ).add( backend_options( "backend1" ) );
}
/*_________________________________________________*
* About
*_________________________________________________*/
inline
AboutData
makeAbout()
{
AboutData about( "Test_Twolaplaciansdistributed" ,
"Test_Twolaplaciansdistributed" ,
"0.1",
"test two laplacian distributed",
Feel::AboutData::License_GPL,
"Copyright (c) 2011 Universite Joseph Fourier" );
about.addAuthor( "Vincent Chabannes", "developer", "vincent.chabannes@feelpp.org", "" );
return about;
}
/*_________________________________________________*
* About
*_________________________________________________*/
void run()
{
typedef Mesh< Simplex<2,1,2> > mesh_type;
double meshSize = doption(_name="hsize");
const int domainColor1=1;
const int domainColor2=2;
const int nTotalProc = Environment::worldComm().size();
auto worldCommDomain1 = Environment::worldComm();
auto worldCommDomain2 = Environment::worldComm();
if ( nTotalProc>1)
{
std::vector<int> MapWorld(nTotalProc);
for (int proc = 0 ; proc < nTotalProc; ++proc)
{
if (proc < nTotalProc/2 )
MapWorld[proc] = domainColor1;
else
MapWorld[proc] = domainColor2;
}
auto worldCommColored = WorldComm(MapWorld);
//worldCommColored.showMe();
worldCommDomain1 = worldCommColored.subWorldComm(domainColor1);
worldCommDomain2 = worldCommColored.subWorldComm(domainColor2);
}
//worldCommDomain1.showMe();
//worldCommDomain2.showMe();
//Environment::worldComm().showMe();
//---------------------------------------------------------------------------------------//
// meshes
GeoTool::Node x11( 0,0 );
GeoTool::Node x12( 1,1 );
GeoTool::Rectangle Omega1( meshSize,"Omega",x11,x12 );
Omega1.setMarker( _type="line",_name="Paroi",_markerAll=true );
Omega1.setMarker( _type="surface",_name="Omega",_markerAll=true );
auto mesh1 = Omega1.createMesh(_mesh=new mesh_type,
_name ="omega1_"+ mesh_type::shape_type::name(),
_partitions=worldCommDomain1.localSize(),
_worldcomm=worldCommDomain1 );
GeoTool::Node x21( -1,0.3 );
GeoTool::Node x22( -1+0.5,0.3 );
GeoTool::Circle Omega2( meshSize,"Omega",x21,x22 );
Omega2.setMarker( _type="line",_name="Paroi",_markerAll=true );
Omega2.setMarker( _type="surface",_name="Omega",_markerAll=true );
auto mesh2 = Omega2.createMesh(_mesh=new mesh_type,
_name ="omega2_"+ mesh_type::shape_type::name(),
_partitions=worldCommDomain2.localSize(),
_worldcomm=worldCommDomain2 );
//---------------------------------------------------------------------------------------//
// functionspaces
auto Xh1 = Pch<1>( mesh1 );
auto Xh2 = Pch<1>( mesh2 );
auto u1 = Xh1->element();
auto u2 = Xh2->element();
//---------------------------------------------------------------------------------------//
// backends
auto backend1 = backend(_rebuild=true,_worldcomm=Xh1->worldComm());
auto backend2 = backend(_rebuild=true,_worldcomm=Xh2->worldComm());
//---------------------------------------------------------------------------------------//
// init matrix and vectors
auto A1 = backend1->newMatrix(_test=Xh1,_trial=Xh1);
auto A2 = backend2->newMatrix(_test=Xh2,_trial=Xh2);
auto F1 = backend1->newVector(Xh1);
auto F2 = backend2->newVector(Xh2);
//---------------------------------------------------------------------------------------//
// assembly
form2( _test=Xh1, _trial=Xh1, _matrix=A1 ) +=
integrate( _range=elements(mesh1), _expr=gradt(u1)*trans(grad(u1)) );
form1( _test=Xh1, _vector=F1 ) +=
integrate( _range=elements(mesh1), _expr=id(u1) );
if (Xh1->worldComm().isActive()) //(todo vincent : simplifier)
{
//A1->close();F1->close();
form2( _test=Xh1, _trial=Xh1, _matrix=A1 ) +=
on( _range=boundaryfaces(mesh1),
_element=u1,_rhs=F1,
_expr=cst(0.) );
}
form2( _test=Xh2, _trial=Xh2, _matrix=A2 ) +=
integrate( _range=elements(mesh2), _expr=gradt(u2)*trans(grad(u2)) );
form1( _test=Xh2, _vector=F2 ) +=
integrate( _range=elements(mesh2), _expr=id(u2) );
if (Xh2->worldComm().isActive()) //(todo vincent : simplifier)
{
//A2->close();F2->close();
form2( _test=Xh2, _trial=Xh2, _matrix=A2 ) +=
on( _range=boundaryfaces(mesh2),
_element=u2,_rhs=F2,
_expr=cst(0.) );
}
//---------------------------------------------------------------------------------------//
// solve (todo vincent : simplifier (+rajouter le worldcomm automatiquement dans prec defaut de la fonction solve et nlsolve )
if (Xh1->worldComm().isActive())
{
backend1->solve(_matrix=A1,_solution=u1,_rhs=F1/*,_prec=prec1*/);
double l2NormU1 = u1.l2Norm();
#if USE_BOOST_TEST
BOOST_CHECK( l2NormU1 > 1e-5 );
#else
CHECK( l2NormU1 > 1e-5 ) << "l2Norm of u1 must be non null";
#endif
}
if (Xh2->worldComm().isActive())
{
backend2->solve(_matrix=A2,_solution=u2,_rhs=F2/*,_prec=prec2*/);
double l2NormU2 = u2.l2Norm();
#if USE_BOOST_TEST
BOOST_CHECK( l2NormU2 > 1e-5 );
#else
CHECK( l2NormU2 > 1e-5 ) << "l2Norm of u2 must be non null";
#endif
}
//---------------------------------------------------------------------------------------//
// exports
auto myexporter1 = exporter( _mesh=mesh1, _name="MyExportDomain1" );
myexporter1->step(0)->add( "u1", u1 );
myexporter1->save();
auto myexporter2 = exporter( _mesh=mesh2, _name="MyExportDomain2" );
myexporter2->step(0)->add( "u2", u2 );
myexporter2->save();
}
} //namespace test_twolaplaciansdistributed
/*_________________________________________________*
* Main
*_________________________________________________*/
#if USE_BOOST_TEST
FEELPP_ENVIRONMENT_WITH_OPTIONS( test_twolaplaciansdistributed::makeAbout(),
test_twolaplaciansdistributed::makeOptions() )
BOOST_AUTO_TEST_SUITE( twolaplaciansdistributed )
BOOST_AUTO_TEST_CASE( twolaplaciansdistributed1 )
{
test_twolaplaciansdistributed::run();
}
BOOST_AUTO_TEST_SUITE_END()
#else
int
main( int argc, char** argv )
{
Feel::Environment env( Feel::_argc=argc,Feel::_argv=argv,
Feel::_about=test_twolaplaciansdistributed::makeAbout(),
Feel::_desc=test_twolaplaciansdistributed::makeOptions() );
test_twolaplaciansdistributed::run();
}
#endif