/* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=cpp:et:sw=4:ts=4:sts=4

  This file is part of the Feel library

  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
  Author(s): Cecile Daversin  <cecile.daversin@lncmi.cnrs.fr>
       Date: 2011-12-07

  Copyright (C) 2011 UJF
  Copyright (C) 2011 CNRS

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 3.0 of the License, or (at your option) any later version.

  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
*/
/**
   \file test_hcurl.cpp
   \author Christophe Prud'homme <christophe.prudhomme@feelpp.org>
   \author Cecile Daversin <cecile.daversin@lncmi.cnrs.fr>
   \date 2011-12-07
 */
// make sure that the init_unit_test function is defined by UTF
//#define BOOST_TEST_MAIN
// give a name to the testsuite
#define BOOST_TEST_MODULE H_curl approximation
// disable the main function creation, use our own
//#define BOOST_TEST_NO_MAIN

#include <feel/feelcore/testsuite.hpp>

#include <feel/feelalg/backend.hpp>
#include <feel/feelfilters/creategmshmesh.hpp>
#include <feel/feelfilters/domain.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/form.hpp>
#include <feel/feelvf/operators.hpp>
#include <feel/feelvf/operations.hpp>
#include <feel/feelvf/on.hpp>
#include <feel/feelpoly/nedelec.hpp>
#include <feel/feelvf/vf.hpp>
#include <feel/feelvf/print.hpp>
#include <feel/feeldiscr/projector.hpp>
#include <feel/feeldiscr/ned1h.hpp>
#include <feel/feeldiscr/pchv.hpp>

#include <feel/feelfilters/loadmesh.hpp>

namespace Feel
{

/**
 * This routine returns the list of options using the
 * boost::program_options library. The data returned is typically used
 * as an argument of a Feel::Application subclass.
 *
 * \return the list of options
 */
inline
po::options_description
makeOptions()
{
    po::options_description testhcurloptions( "test h_curl options" );
    testhcurloptions.add_options()
        //( "hsize", po::value<double>()->default_value( 0.1 ), "mesh size" )
        //( "xmin", po::value<double>()->default_value( -1 ), "xmin of the reference element" )
        //( "ymin", po::value<double>()->default_value( -1 ), "ymin of the reference element" )
        //( "zmin", po::value<double>()->default_value( -1 ), "zmin of the reference element" )
        ( "tol", po::value<double>()->default_value( 2.5e-1 ), "tolerance for the solution" )
        ( "u_exact_2d", po::value<std::string>()->default_value("{1-y^2,1-x^2}:x:y"), "exact solution" )
        ( "f_exact_2d", po::value<std::string>()->default_value("{3-y^2,3-x^2}:x:y"), "exact rhs" )
        ( "u_exact_3d", po::value<std::string>()->default_value("{(1-y^2)*(1-z^2),(1-x^2)*(1-z^2),(1-x^2)*(1-y^2)}:x:y:z"), "exact solution" )
        ( "f_exact_3d", po::value<std::string>()->default_value("{y^2*z^2-3*y^2-3*z^2+5,x^2*z^2-3*x^2-3*z^2+5,x^2*y^2-3*x^2-3*y^2+5}:x:y:z"), "exact rhs" )
    ;
    return testhcurloptions.add( Feel::feel_options() );
}
#if 0
inline
AboutData
makeAbout()
{
    AboutData about( "test_hcurl" ,
                     "test_hcurl" ,
                     "0.1",
                     "Test for h_curl space (Dim=2 Order=1)",
                     AboutData::License_GPL,
                     "Copyright (c) 2009 Universite Joseph Fourier" );
    about.addAuthor( "Cecile Daversin", "developer", "cecile.daversin@lncmi.cnrs.fr", "" );
    about.addAuthor( "Christophe Prud'homme", "developer", "christophe.prudhomme@feelpp.org", "" );
    return about;

}
#endif
using namespace Feel;

template<int Dim>
class TestHCurl
    :
public Application
{
    typedef Application super;

public:

    //! numerical type is double
    typedef double value_type;

    //! linear algebra backend factory
    typedef Backend<value_type> backend_type;
    //! linear algebra backend factory shared_ptr<> type
    typedef typename std::shared_ptr<backend_type> backend_ptrtype ;
    //! sparse matrix type associated with backend
    typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
    //! sparse matrix type associated with backend (shared_ptr<> type)
    typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    //! vector type associated with backend
    typedef typename backend_type::vector_type vector_type;
    //! vector type associated with backend (shared_ptr<> type)
    typedef typename backend_type::vector_ptrtype vector_ptrtype;

    //! geometry entities type composing the mesh, here Simplex in Dimension Dim of Order G_order
    typedef Simplex<Dim,1> convex_type;
    //! mesh type
    typedef Mesh<convex_type> mesh_type;
    //! mesh shared_ptr<> type
    typedef std::shared_ptr<mesh_type> mesh_ptrtype;

    //! the basis type of our approximation space
    //typedef bases<Nedelec<0> > basis_type;
    typedef bases<Nedelec<0,NedelecKind::NED1> > basis_type;
    typedef bases<Lagrange<1,Scalar> > lagrange_basis_s_type; //P1 scalar space
    typedef bases<Lagrange<1,Vectorial> > lagrange_basis_v_type; //P1 vectorial space

    //! the approximation function space type
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    typedef FunctionSpace<mesh_type, lagrange_basis_s_type> lagrange_space_s_type;
    typedef FunctionSpace<mesh_type, lagrange_basis_v_type> lagrange_space_v_type;
    //! the approximation function space type (shared_ptr<> type)
    typedef std::shared_ptr<space_type> space_ptrtype;
    typedef std::shared_ptr<lagrange_space_s_type> lagrange_space_s_ptrtype;
    typedef std::shared_ptr<lagrange_space_v_type> lagrange_space_v_ptrtype;
    //! an element type of the approximation function space
    typedef typename space_type::element_type element_type;

    //! the exporter factory type
    typedef Exporter<mesh_type> export_type;
    //! the exporter factory (shared_ptr<> type)
    typedef std::shared_ptr<export_type> export_ptrtype;

    /**
     * Constructor
     */
    TestHCurl()
        :
        super(),
        M_backend( backend_type::build( soption( _name="backend" ) ) ),
        meshSize( doption("gmsh.hsize") ),
        M_penaldir(1e5)
    {
        Feel::cout << "[TestHCurl]\n";
        Feel::cout << "hsize = " << meshSize << std::endl;

        this->changeRepository( boost::format( "%1%/h_%2%/" )
                                % this->about().appName()
                                % doption("gmsh.hsize") );

        mesh = loadMesh(_mesh = new mesh_type);
        M_Xh = space_type::New(mesh);

        M_a = M_backend->newMatrix(_test=M_Xh,_trial=M_Xh);
        M_f = M_backend->newVector(_test=M_Xh);
    }

    /**
     * run the application
     */
    //void twoElementsMesh();
    //void eightElementsMesh();
    void testProjector();
    void exampleProblem1();

private:
    mesh_ptrtype mesh;
    space_ptrtype M_Xh;

    //! linear algebra backend
    backend_ptrtype M_backend;
    vector_ptrtype M_f;
    sparse_matrix_ptrtype M_a;

    //! mesh characteristic size
    double meshSize;

    double M_penaldir;

}; //TestHCurl

// Resolve problem curl(curl(u)) + u = f with cross_prod(u,n) = 0 on boundary
template<int Dim>
void
TestHCurl<Dim>::exampleProblem1()
{
    // M_Xh : space build with Nedelec elements
    auto Vh = Pchv<1>( mesh );
    auto u = M_Xh->element();
    auto phi = M_Xh->element();

    //auto u_exact = expr<2,1>( "{1-y*y,1-x*x}:x:y" );
    //auto v = Vh->element();

    auto a = form2( _test=M_Xh, _trial=M_Xh, _matrix=M_a );
    auto l = form1( _test=M_Xh, _vector=M_f );
    //variationnal formulation : curl(curl(u)) + u = f

    std::string u_opt = "u_exact_"+std::to_string(Dim)+"d";
    std::string f_opt = "f_exact_"+std::to_string(Dim)+"d";
    auto u_exact = expr<Dim,1>( soption(u_opt) );
    auto f = expr<Dim,1>( soption(f_opt) );

    l = integrate( _range=elements(mesh), _expr=inner(f,id(phi)) );
    //Dirichlet bc on weak form
    l += integrate( _range=boundaryfaces(mesh),
                    _expr=- inner(curl(phi),cross(u_exact,N()))
                    + M_penaldir*trans( cross(u_exact,N()) )*cross(id(phi),N())/hFace() );

    a = integrate(_range=elements(mesh),
                  _expr=inner(curlt(u),curl(phi)) + inner(idt(u),id(phi)) );
    //a += on( _range=boundaryfaces( mesh ),_element=u, _rhs=l, _expr=cst(0.) );

    //Dirichlet bc on weak form
    a += integrate(_range=boundaryfaces(mesh),
                   _expr=-inner(curlt(u),cross(id(phi),N()))
                   - inner(curl(phi),cross(idt(u),N()))
                   + M_penaldir*trans( cross(idt(u),N()) )*cross(id(phi),N())/hFace() );

    a.solve( _solution=u, _rhs=l, _rebuild=true);

    // L2 projection of solution u
    auto u_L2proj = M_Xh->element();
    auto phi_L2proj = M_Xh->element();

    auto a_L2proj = form2( _test=M_Xh, _trial=M_Xh );
    a_L2proj = integrate( _range=elements(mesh), _expr = trans(idt(u_L2proj))*id(phi_L2proj) );
    auto f_L2proj = form1( _test=M_Xh );
    f_L2proj = integrate( _range=elements(mesh), _expr = trans(idv(u))*id(phi_L2proj) );
    a_L2proj.solve( _solution=u_L2proj, _rhs=f_L2proj, _rebuild=true);

    double errorU = normL2(_range=elements(mesh), _expr= idv(u)-u_exact);
    auto errorProj = normL2(_range=elements(mesh), _expr=idv(u)-idv(u_L2proj) );

    BOOST_CHECK_SMALL( errorProj, 1e-13 );
    BOOST_CHECK_SMALL( errorU, doption("tol") );

#if 0
    auto l2proj = opProjection( _domainSpace=M_Xh, _imageSpace=M_Xh, _type=L2 );
    auto error_hcurl = l2proj->project( _expr=trans( u_exact - idv(u) ) );
#endif

    std::string pro1_name = "problem1";
    std::string exporterName = ( boost::format( "%1%-%2%-%3%" )
                                 % this->about().appName()
                                 % ( boost::format( "%1%-%2%-%3%" ) % "hypercube" % 2 % 1 ).str()
                                 % pro1_name ).str();
    auto exporter_pro1 = exporter(_mesh=mesh,_name=exporterName );
    exporter_pro1->addRegions();
    exporter_pro1->add( "u", u );
    exporter_pro1->add( "proj_L2_u", u_L2proj );
    exporter_pro1->add( "u_exact", u_exact );
    exporter_pro1->save();

}

#if 0
void
TestHCurl::testProjector()
{
    mesh_ptrtype mesh = loadMesh(_mesh = new mesh_type);
    auto Nh = Ned1h<0>( mesh );
    lagrange_space_v_ptrtype Yh_v = lagrange_space_v_type::New( mesh ); //lagrange vectorial space
    lagrange_space_s_ptrtype Yh_s = lagrange_space_s_type::New( mesh ); //lagrange scalar space

    auto E = Py()*unitX() + Px()*unitY();
    auto f = cst(0.); //curld[2d](E) = f

    // curl2D(u1,u2) = d(u2)/dx1 - d(u1)/dx2
    // 2D case : only curlx is initialized to curl2D(u1,u2) (curl is a vector with only one component initialized)

    // L2 projection (Lagrange)
    auto l2_lagV = opProjection( _domainSpace=Yh_v, _imageSpace=Yh_v, _type=L2 ); //l2 vectorial proj
    auto l2_lagS = opProjection( _domainSpace=Yh_s, _imageSpace=Yh_s, _type=L2 ); //l2 scalar proj
    auto E_pL2_lag = l2_lagV->project( _expr= trans(E) );
    auto error_pL2_lag = l2_lagS->project( _expr=curlxv(E_pL2_lag) - f );

    // L2 projection (Nedelec)
    auto l2_ned = opProjection( _domainSpace=Nh, _imageSpace=Nh, _type=L2 );
    auto E_pL2_ned = l2_ned->project( _expr= trans(E) );
    auto error_pL2_ned = l2_lagS->project( _expr=curlxv(E_pL2_lag) - f );

    // H1 projection (Lagrange)
    auto h1_lagV = opProjection( _domainSpace=Yh_v, _imageSpace=Yh_v, _type=H1 ); //h1 vectorial proj
    auto h1_lagS = opProjection( _domainSpace=Yh_s, _imageSpace=Yh_s, _type=H1 ); //h1 scalar proj
    auto E_pH1_lag = h1_lagV->project( _expr= trans(E), _grad_expr=mat<2,2>(cst(0.),cst(1.),cst(1.),cst(0.)) );
    auto error_pH1_lag = l2_lagS->project( _expr=curlxv(E_pH1_lag) - f );

    // H1 projection (Nedelec)
    auto h1_ned = opProjection( _domainSpace=Nh, _imageSpace=Nh, _type=H1 ); //h1 vectorial proj
    auto E_pH1_ned = h1_ned->project( _expr= trans(E), _grad_expr=mat<2,2>(cst(0.),cst(1.),cst(1.),cst(0.)) );
    auto error_pH1_ned = l2_lagS->project( _expr=curlxv(E_pH1_ned) - f );

    // HCURL projection (Lagrange)
    auto hcurl_lagV = opProjection( _domainSpace=Yh_v, _imageSpace=Yh_v, _type=HCURL );
    auto hcurl_lagS = opProjection( _domainSpace=Yh_s, _imageSpace=Yh_s, _type=HCURL );
    auto E_pHCURL_lag = hcurl_lagV->project( _expr= trans(E), _div_expr=cst(0.) );
    auto error_pHCURL_lag = l2_lagS->project( _expr=curlxv(E_pHCURL_lag) - f );

    // HCURL projection (Nedelec)
    auto hcurl = opProjection( _domainSpace=Nh, _imageSpace=Nh, _type=HCURL ); //hdiv proj (RT elts)
    auto E_pHCURL_ned = hcurl->project( _expr= trans(E), _div_expr=cst(0.) );
    auto error_pHCURL_ned = l2_lagS->project( _expr=curlxv(E_pHCURL_ned) - f );

    BOOST_TEST_MESSAGE("L2 projection [Lagrange]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pL2_lag, error_pL2_lag ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pL2_lag, error_pL2_lag ) ), 1e-13 );
    BOOST_TEST_MESSAGE("L2 projection [NED]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pL2_ned, error_pL2_ned ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pL2_ned, error_pL2_ned ) ), 1e-13 );
    BOOST_TEST_MESSAGE("H1 projection [Lagrange]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pH1_lag, error_pH1_lag ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pH1_lag, error_pH1_lag ) ), 1e-13 );
    BOOST_TEST_MESSAGE("H1 projection [NED]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pH1_ned, error_pH1_ned ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pH1_ned, error_pH1_ned ) ), 1e-13 );
    BOOST_TEST_MESSAGE("HDIV projection [Lagrange]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pHCURL_lag, error_pHCURL_lag ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pHCURL_lag, error_pHCURL_lag ) ), 1e-13 );
    BOOST_TEST_MESSAGE("HDIV projection [NED]: error[div(E)-f]");
    std::cout << "error L2: " << math::sqrt( l2_lagS->energy( error_pHCURL_ned, error_pHCURL_ned ) ) << "\n";
    BOOST_CHECK_SMALL( math::sqrt( l2_lagS->energy( error_pHCURL_ned, error_pHCURL_ned ) ), 1e-13 );

    std::string proj_name = "projection";
    std::string exporterName = ( boost::format( "%1%-%2%-%3%" )
                                 % this->about().appName()
                                 % ( boost::format( "%1%-%2%-%3%" ) % "hypercube" % 2 % 1 ).str()
                                 % proj_name ).str();
    auto exporter_proj = exporter(_mesh=mesh,_name=exporterName );
    exporter_proj->step( 0 )->add( "proj_L2_E[Lagrange]", E_pL2_lag );
    exporter_proj->step( 0 )->add( "proj_L2_E[NED]", E_pL2_ned );
    exporter_proj->step( 0 )->add( "proj_H1_E[Lagrange]", E_pH1_lag );
    exporter_proj->step( 0 )->add( "proj_H1_E[NED]", E_pH1_ned );
    exporter_proj->step( 0 )->add( "proj_HDiv_E[Lagrange]", E_pHCURL_lag );
    exporter_proj->step( 0 )->add( "proj_HDiv_E[NED]", E_pHCURL_ned );
    exporter_proj->save();
}
#endif
}


FEELPP_ENVIRONMENT_WITH_OPTIONS( Feel::makeAboutDefault("test_hcurl"), Feel::makeOptions() )

BOOST_AUTO_TEST_SUITE( space )

#if 0
BOOST_AUTO_TEST_CASE( test_hcurl_projection )
{
    BOOST_TEST_MESSAGE( "test_hcurl_N0 : projection on one real element" );
    Feel::TestHCurl t;
    t.testProjector();
    BOOST_TEST_MESSAGE( "test_hcurl_N0 : projection on one real element done" );
}
#endif

BOOST_AUTO_TEST_CASE( test_hcurl_example_2d )
{
    BOOST_TEST_MESSAGE( "test_hcurl on example 1" );
    Feel::TestHCurl<2> t;
    t.exampleProblem1();
    BOOST_TEST_MESSAGE( "test_hcurl_N0 on example 1 done" );
}

BOOST_AUTO_TEST_CASE( test_hcurl_example_3d )
{
    BOOST_TEST_MESSAGE( "test_hcurl on example 1" );
    Feel::TestHCurl<3> t;
    t.exampleProblem1();
    BOOST_TEST_MESSAGE( "test_hcurl_N0 on example 1 done" );
}

BOOST_AUTO_TEST_SUITE_END()