/* -*- 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.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::worldCommPtr();
    auto worldCommDomain2 = Environment::worldCommPtr();

    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->worldCommPtr());
    auto backend2 = backend(_rebuild=true,_worldcomm=Xh2->worldCommPtr());
    //---------------------------------------------------------------------------------------//
    // 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) );
    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) );
    form2( _test=Xh2, _trial=Xh2, _matrix=A2 ) +=
        on( _range=boundaryfaces(mesh2),
            _element=u2,_rhs=F2,
            _expr=cst(0.) );

    //---------------------------------------------------------------------------------------//
    // solve
    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