#include "epot_bicgstabsolver.hpp"
#include "geometry.hpp"
#include "func_solid.hpp"
#include "epot_efield.hpp"
#include "meshvectorfield.hpp"
#include "particledatabase.hpp"
#include "geomplotter.hpp"
#include "ibsimu.hpp"
#include "error.hpp"


bool solid1( double x, double y, double z )
{
    double r = sqrt(x*x+y*y);
    return( z <= 0.02 && r >= 0.018 );
}


bool solid2( double x, double y, double z )
{
    double r = sqrt(x*x+y*y);
    return( z >= 0.03 && z <= 0.04 && r >= 0.02 );
}


bool solid3( double x, double y, double z )
{
    double r = sqrt(x*x+y*y);
    return( z >= 0.06 && r >= 0.03 && r >= 0.07 - 0.5*z );
}


void simu( void )
{
    Geometry geom( MODE_3D, Int3D(101,101,241), Vec3D(0,0,0), 5e-4 );

    Solid *s1 = new FuncSolid( solid1 );
    geom.set_solid( 7, s1 );
    Solid *s2 = new FuncSolid( solid2 );
    geom.set_solid( 8, s2 );
    Solid *s3 = new FuncSolid( solid3 );
    geom.set_solid( 9, s3 );

    geom.set_boundary( 1, Bound(BOUND_NEUMANN,     0.0  ) );
    geom.set_boundary( 2, Bound(BOUND_NEUMANN,     0.0  ) );
    geom.set_boundary( 3, Bound(BOUND_NEUMANN,     0.0  ) );
    geom.set_boundary( 4, Bound(BOUND_NEUMANN,     0.0  ) );
    geom.set_boundary( 5, Bound(BOUND_DIRICHLET,  -3.0e3) );
    geom.set_boundary( 6, Bound(BOUND_DIRICHLET,  -1.0e3) );

    geom.set_boundary( 7, Bound(BOUND_DIRICHLET,  -3.0e3) );
    geom.set_boundary( 8, Bound(BOUND_DIRICHLET, -14.0e3) );
    geom.set_boundary( 9, Bound(BOUND_DIRICHLET,  -1.0e3) );

    geom.build_mesh();
    geom.save( "geom.dat" );
}


int main( int argc, char **argv )
{
    try {
	ibsimu.set_message_threshold( MSG_VERBOSE, 1 );
	ibsimu.set_thread_count( 4 );
        simu();
    } catch ( Error e ) {
        e.print_error_message( ibsimu.message( 0 ) );
        exit( 1 );
    }

    return( 0 );
}