https://github.com/feelpp/feelpp
Tip revision: 94701fa62fba97dbfafab22e917f6b3bcfcf08d2 authored by Christophe Prud'homme on 03 March 2013, 10:33:36 UTC
Merge branch 'release/version-0.92'
Merge branch 'release/version-0.92'
Tip revision: 94701fa
test_twolaplaciansdistributed.cpp
/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*-*/
#define TWOLAP_USE_BOOST_TEST 0
#if TWOLAP_USE_BOOST_TEST
#define BOOST_TEST_MODULE test_twolaplaciansdistributed
#include <testsuite/testsuite.hpp>
#endif
#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.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;
typedef Application Application_type;
typedef boost::shared_ptr<Application_type> Application_ptrtype;
/*_________________________________________________*
* 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( "laplacian" ) );
}
/*_________________________________________________*
* About
*_________________________________________________*/
inline
AboutData
makeAbout()
{
AboutData about( "Test_Twolaplaciansdistributed" ,
"Test_Twolaplaciansdistributed" ,
"0.1",
"test bilinear and linear form for several trial/test mesh ",
Feel::AboutData::License_GPL,
"Copyright (c) 2011 Universite Joseph Fourier" );
about.addAuthor( "Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "" );
return about;
}
/*_________________________________________________*
* About
*_________________________________________________*/
void run( Application_ptrtype & theApp )
{
typedef Mesh< Simplex<2,1,2> > mesh_type;
double meshSize = theApp->vm()["hsize"].as<double>();
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);
if (nTotalProc==1) // a particular test case
{
MapWorld[0]=domainColor2;
MapWorld[1]=domainColor1;
MapWorld[2]=domainColor2;
MapWorld[3]=domainColor2;
MapWorld[4]=domainColor1;
}
else
{
for (int proc = 0 ; proc < nTotalProc; ++proc)
{
if (proc < nTotalProc/2 )
MapWorld[proc] = domainColor1;
else
MapWorld[proc] = domainColor2;
}
}
auto worldCommColored = WorldComm(MapWorld);
worldCommDomain1 = worldCommColored.subWorldComm(domainColor1);
worldCommDomain2 = worldCommColored.subWorldComm(domainColor2);
}
//worldCommColored.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
typedef Lagrange<1, Scalar,Continuous,PointSetFekete> basis_type;
typedef FunctionSpace<mesh_type, bases<basis_type> > space_type;
auto Xh1 = space_type::New(_mesh=mesh1,
_worldscomm=std::vector<WorldComm>(1,worldCommDomain1) );
auto Xh2 = space_type::New(_mesh=mesh2,
_worldscomm=std::vector<WorldComm>(1,worldCommDomain2) );
auto u1 = Xh1->element();
auto u2 = Xh2->element();
//---------------------------------------------------------------------------------------//
// backends
auto backend1 = Backend<double>::build(theApp->vm(),"",Xh1->worldComm());
auto backend2 = Backend<double>::build(theApp->vm(),"",Xh2->worldComm());
//---------------------------------------------------------------------------------------//
// init matrix and vectors (todo vincent : simplifier (tres facile) )
typedef Backend<double> backend_type;
typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
typedef backend_type::vector_ptrtype vector_ptrtype;
sparse_matrix_ptrtype A1,A2;
vector_ptrtype F1,F2;
if (Xh1->worldComm().isActive())
{
A1 = backend1->newMatrix(_test=Xh1,_trial=Xh1);
F1 = backend1->newVector(Xh1);
}
if (Xh2->worldComm().isActive())
{
A2 = backend2->newMatrix(_test=Xh2,_trial=Xh2);
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())
{
auto prec1 = preconditioner(_pc=(PreconditionerType) backend1->pcEnumType() /*LU_PRECOND*/,
_matrix=A1,
_backend= backend1,
_pcfactormatsolverpackage=(MatSolverPackageType) backend1->matSolverPackageEnumType(),
_worldcomm=backend1->comm(),
_prefix=backend1->prefix() );
backend1->solve(_matrix=A1,_solution=u1,_rhs=F1,_prec=prec1);
}
if (Xh2->worldComm().isActive())
{
auto prec2 = preconditioner(_pc=(PreconditionerType) backend2->pcEnumType() /*LU_PRECOND*/,
_matrix=A2,
_backend= backend2,
_pcfactormatsolverpackage=(MatSolverPackageType) backend2->matSolverPackageEnumType(),
_worldcomm=backend2->comm(),
_prefix=backend2->prefix() );
backend2->solve(_matrix=A2,_solution=u2,_rhs=F2,_prec=prec2);
}
//---------------------------------------------------------------------------------------//
// exports
#if 1
auto myexporter1 = Exporter<mesh_type>::New( theApp->vm(), "MyExportDomain1", Xh1->worldComm() );
myexporter1->step(0)->setMesh( mesh1 );
myexporter1->step(0)->add( "u1", u1 );
myexporter1->save();
auto myexporter2 = Exporter<mesh_type>::New( theApp->vm(), "MyExportDomain2", Xh2->worldComm() );
myexporter2->step(0)->setMesh( mesh2 );
myexporter2->step(0)->add( "u2", u2 );
myexporter2->save();
#endif
}
} //namespace test_twolaplaciansdistributed
/*_________________________________________________*
* Main
*_________________________________________________*/
#if TWOLAP_USE_BOOST_TEST
FEELPP_ENVIRONMENT_WITH_OPTIONS( test_twolaplaciansdistributed::makeAbout(),
test_twolaplaciansdistributed::makeOptions() )
BOOST_AUTO_TEST_SUITE( twolaplaciansdistributed )
typedef Feel::Application Application_type;
typedef boost::shared_ptr<Application_type> Application_ptrtype;
BOOST_AUTO_TEST_CASE( twolaplaciansdistributed1 )
{
auto theApp = Application_ptrtype( new Application_type );
if ( theApp->vm().count( "help" ) )
{
std::cout << theApp->optionsDescription() << "\n";
exit( 0 );
}
theApp->changeRepository( boost::format( "/testsuite/feeldiscr/%1%/" )
% theApp->about().appName() );
test_twolaplaciansdistributed::run( theApp );
}
BOOST_AUTO_TEST_SUITE_END()
#else
int
main( int argc, char** argv )
{
typedef Feel::Application Application_type;
typedef boost::shared_ptr<Application_type> Application_ptrtype;
Feel::Environment env( Feel::_argc=argc,Feel::_argv=argv,
Feel::_about=test_twolaplaciansdistributed::makeAbout(),
Feel::_desc=test_twolaplaciansdistributed::makeOptions() );
auto theApp = Application_ptrtype( new Application_type );
if ( theApp->vm().count( "help" ) )
{
std::cout << theApp->optionsDescription() << "\n";
exit( 0 );
}
theApp->changeRepository( boost::format( "/testsuite/feeldiscr/%1%/" )
% theApp->about().appName() );
test_twolaplaciansdistributed::run( theApp );
}
#endif