diff --git a/docs/modules/ROOT/examples/02-environment.cpp b/docs/modules/ROOT/examples/02-environment.cpp index 278f5a6..43884fd 100644 --- a/docs/modules/ROOT/examples/02-environment.cpp +++ b/docs/modules/ROOT/examples/02-environment.cpp @@ -1,6 +1,12 @@ #include #include "ut.hpp" +// Just for information "ut.hpp" file is missing. Do not forget to install the lib boots +// +// https://github.com/boost-ext/ut/tree/master + + + int main( int argc, char* argv[] ) { using namespace Feel; diff --git a/docs/modules/ROOT/examples/07-myfunctionspace.cpp b/docs/modules/ROOT/examples/07-myfunctionspace.cpp index 7bf3def..2144fd3 100644 --- a/docs/modules/ROOT/examples/07-myfunctionspace.cpp +++ b/docs/modules/ROOT/examples/07-myfunctionspace.cpp @@ -39,6 +39,8 @@ #include #include #include +#include "ut.hpp" + using namespace Feel; @@ -51,6 +53,8 @@ int main( int argc, char** argv ) _author="Feel++ Consortium", _email="feelpp-devel@feelpp.org" ) ); + + // tag::mesh[] // create the mesh auto mesh = loadMesh(_mesh=new Mesh>); @@ -61,11 +65,14 @@ int main( int argc, char** argv ) auto Xh = Pch<2>( mesh ); // end::space[] + // tag::expression[] auto g = expr<4>( soption(_name="functions.g")); auto gradg = grad<3>(g); // end::expression[] + + // tag::interpolant[] // tag::element[] // elements of \( u,w \in X_h \) @@ -79,13 +86,30 @@ int main( int argc, char** argv ) // build the interpolant of the interpolation error w.on( _range=elements( mesh ), _expr=idv( u )-g ); + + // compute L2 norms \(||\cdot||_{L^2}\) - double L2g = normL2( elements( mesh ), g ); - double H1g = normL2( elements( mesh ), _expr=g,_grad_expr=gradg ); - double L2uerror = normL2( elements( mesh ), ( idv( u )-g ) ); - double H1uerror = normH1( elements( mesh ), _expr=( idv( u )-g ), - _grad_expr=( gradv( u )-gradg ) ); + double L2g = normL2(_range=elements( mesh ), _expr=g ); + double H1g = normL2(_range=elements( mesh ), _expr=g,_grad_expr=gradg ); + + double L2uerror = normL2(_range=elements( mesh ), _expr=( idv( u )-g ) ); + + + auto range=elements(mesh); //OK + auto expr=idv(u)-g; //OK + auto grad_expr=gradv( u )-gradg; //OK + double H1uerror = normH1(_range=range, + _expr=expr, + _grad_expr=grad_expr + ); + + +/* + + + + Feel::cout << "||u-g||_0 = " << L2uerror/L2g << std::endl; Feel::cout << "||u-g||_1 = " << H1uerror/H1g << std::endl; // end::interpolant[] @@ -103,5 +127,7 @@ int main( int argc, char** argv ) e->save(); // end::export[] +*/ + } // end::all[] diff --git a/docs/modules/ROOT/examples/10-model2.cpp b/docs/modules/ROOT/examples/10-model2.cpp index d9b652a..5764774 100644 --- a/docs/modules/ROOT/examples/10-model2.cpp +++ b/docs/modules/ROOT/examples/10-model2.cpp @@ -1,19 +1,28 @@ #include #include + +const int MODEL_DIM = 3; + + int main(int argc, char**argv ) { using namespace Feel; po::options_description laplacianoptions( "Laplacian options" ); laplacianoptions.add_options() - ("myVerbose", po::value< bool >()-> default_value( true ), "Display information during execution") - ; + ("myVerbose", po::value< bool >()-> default_value( true ), "Display information during execution"); + Environment env( _argc=argc, _argv=argv, _desc=laplacianoptions, _about=about(_name="aniso_laplacian", _author="Feel++ Consortium", _email="feelpp-devel@feelpp.org")); + ModelProperties model; // Will load --mod-file + map_scalar_field<2> bc_u { model.boundaryConditions().getScalarFields<2>("heat","dirichlet") }; + + + ModelMaterials materials = model.materials(); if(boption("myVerbose") && Environment::isMasterRank() ) std::cout << "Model " << Environment::expand( soption("mod-file")) << " loaded." << std::endl; @@ -25,6 +34,7 @@ int main(int argc, char**argv ) auto k11 = Vh->element(); auto k12 = Vh->element(); auto k22 = Vh->element(); + #if MODEL_DIM == 3 auto k13 = Vh->element(); auto k23 = Vh->element(); @@ -35,27 +45,31 @@ int main(int argc, char**argv ) l = integrate(_range=elements(mesh),_expr=f*id(v)); for(auto it : materials) { - auto mat = material(it); - if(boption("myVerbose") && Environment::isMasterRank() ) - std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first << " with diffusion coef [" -#if MODEL_DIM == 3 - << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "]," - << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "]," - << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]" -#else - << "[" << it.second.k11() << "," << it.second.k12() << "]," - << "[" << it.second.k12() << "," << it.second.k22() << "]]" -#endif - << std::endl; - k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11())); - k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12())); - k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22())); -#if MODEL_DIM == 3 - k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13()); - k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23()); - k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33()); -#endif + + + auto mat = material(it); + if(boption("myVerbose") && Environment::isMasterRank() ) + std::cout << "[Materials] - Laoding data for " << it.second.name() << " that apply on marker " << it.first << " with diffusion coef [" + #if MODEL_DIM == 3 + << "[" << it.second.k11() << "," << it.second.k12() << "," << it.second.k13() << "]," + << "[" << it.second.k12() << "," << it.second.k22() << "," << it.second.k23() << "]," + << "[" << it.second.k13() << "," << it.second.k23() << "," << it.second.k33() << "]]" + #else + << "[" << it.second.k11() << "," << it.second.k12() << "]," + << "[" << it.second.k12() << "," << it.second.k22() << "]]" + #endif + << std::endl; + k11.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k11())); + k12.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k12())); + k22.on(_range=markedelements(mesh,it.first),_expr=cst(it.second.k22())); + #if MODEL_DIM == 3 + k13 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k13()); + k23 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k23()); + k33 += vf::project(_space=Vh,_range=markedelements(mesh,marker(it)),_expr=mat.k33()); + #endif } + + #if MODEL_DIM == 2 a += integrate(_range=elements(mesh),_expr=inner(mat(idv(k11), idv(k12), idv(k12), idv(k22) )*trans(gradt(u)),trans(grad(v))) ); #else @@ -66,6 +80,7 @@ int main(int argc, char**argv ) std::cout << "[BC] - Applying " << it.second << " on " << it.first << std::endl; a+=on(_range=markedfaces(mesh,it.first), _rhs=l, _element=u, _expr=it.second ); } + a.solve(_rhs=l,_solution=u); auto e = exporter( _mesh=mesh ); for(int i = 0; i < 3; i ++){ diff --git a/src2/Cases/Circle/Data2d.cfg b/src2/Cases/Circle/Data2d.cfg new file mode 100644 index 0000000..8055510 --- /dev/null +++ b/src2/Cases/Circle/Data2d.cfg @@ -0,0 +1,29 @@ + + + +[functions] +# Dirichlet +g=x^2+y^2:x:y +# right hand side +f=-4 +# Robin left hand side +a=1 +# Robin right hand side +b=2*(x*nx+y*ny)+x^2+y^2:x:y:nx:ny +# Neumann +c=2*(x*nx+y*ny):x:y:nx:ny +# mu: diffusion term (laplacian) +mu=1 + +[checker] +check=true +solution=x**2+y**2 +filename=$cfgdir/checker.json +script=$cfgdir/../python/laplacian.py + +[gmsh] +filename=$cfgdir/Data2d.geo + +[exporter] +format=ensightgold +geometry=static diff --git a/src2/Cases/Circle/Data2d.geo b/src2/Cases/Circle/Data2d.geo new file mode 100644 index 0000000..cca75ac --- /dev/null +++ b/src2/Cases/Circle/Data2d.geo @@ -0,0 +1,17 @@ +h=0.1; + +// Gmsh project created +Point(1) = {0, 0, 0, h}; +Point(2) = {1, 0, 0, h}; +Point(3) = {-1, 0, 0, h}; +Point(4) = {0, 1, 0, h}; +Point(5) = {0, -1, 0, h}; +Circle(1) = {2, 1, 4}; +Circle(2) = {4, 1, 3}; +Circle(3) = {3, 1, 5}; +Circle(5) = {5, 1, 2}; +Line Loop(6) = {1, 2, 3, 5}; +Plane Surface(7) = {6}; +Physical Line("Dirichlet") = {1, 2, 3, 5}; +Physical Surface(8) = {7}; + diff --git a/src2/Laplacian.cpp b/src2/Laplacian.cpp new file mode 100644 index 0000000..a6cc20d --- /dev/null +++ b/src2/Laplacian.cpp @@ -0,0 +1,144 @@ +#include + + +const int FEELPP_DIM=2; + + +int main(int argc, char**argv ) +{ + using namespace Feel; + using Feel::cout; + po::options_description laplacianoptions( "Laplacian options" ); + laplacianoptions.add_options() + ( "no-solve", po::value()->default_value( false ), "No solve" ); + + Environment env( _argc=argc, _argv=argv, + _desc=laplacianoptions, + _about=about(_name="qs_laplacian", + _author="Feel++ Consortium", + _email="feelpp-devel@feelpp.org")); + + + Feel::cout << "proc " << Environment::rank() + <<" of "<< Environment::numberOfProcessors() + << std::endl; + + tic(); + auto mesh = loadMesh(_mesh=new Mesh>); + toc("loadMesh"); + + tic(); + auto Vh = Pch<2>( mesh ); + auto u = Vh->element("u"); + auto mu = expr(soption(_name="functions.mu")); // diffusion term + auto f = expr( soption(_name="functions.f"), "f" ); + auto r_1 = expr( soption(_name="functions.a"), "a" ); // Robin left hand side expression + auto r_2 = expr( soption(_name="functions.b"), "b" ); // Robin right hand side expression + auto n = expr( soption(_name="functions.c"), "c" ); // Neumann expression + + std::map inputs{{"dim",std::to_string(dimension(mesh))}, + {"k",soption("k")},{"r_1",soption("r_1")},{"u",""},{"un",soption("un")},{"f",soption("f")},{"g",soption("g")},{"r_2",soption("r_2")}}; + + + auto thechecker = checker( _name= "L1/H1 convergence", + _solution_key="p", + _gradient_key="grad_p", + _inputs=inputs); + + + auto solution = expr(thechecker.solution(), "solution" ); + auto g = thechecker.check()?solution:expr( soption(_name="functions.g"), "g" ); + + auto v = Vh->element( g, "g" ); + toc("Vh"); + + + + tic(); + auto l = form1( _test=Vh ); + l = integrate(_range=elements(mesh), + _expr=f*id(v)); + l+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_2*id(v)); + l+=integrate(_range=markedfaces(mesh,"Neumann"), _expr=n*id(v)); + toc("l"); + + //Until here OK + + + tic(); + auto a = form2( _trial=Vh, _test=Vh); + tic(); + a = integrate(_range=elements(mesh), + _expr=mu*inner(gradt(u),grad(v)) ); + toc("a"); + a+=integrate(_range=markedfaces(mesh,"Robin"), _expr=r_1*idt(u)*id(v)); + a+=on(_range=markedfaces(mesh,"Dirichlet"), _rhs=l, _element=u, _expr=g ); + //! if no markers Robin Neumann or Dirichlet are present in the mesh then + //! impose Dirichlet boundary conditions over the entire boundary + if ( !mesh->hasAnyMarker({"Robin", "Neumann","Dirichlet"}) ) + a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u, _expr=g ); + toc("a"); + + + tic(); + //! solve the linear system, find u s.t. a(u,v)=l(v) for all v + if ( !boption( "no-solve" ) ) + a.solve(_rhs=l,_solution=u); + toc("a.solve"); + + + + tic(); + auto e = exporter( _mesh=mesh ); + e->addRegions(); + e->add( "uh", u ); + + + + if (thechecker.check()) + { + v.on(_range=elements(mesh), _expr=solution ); + e->add( "solution", v ); + } + + e->save(); + toc("Exporter"); + + + // compute l2 and h1 norm of u-u_h where u=solution + auto norms = [=]( std::string const& solution ) ->std::map + { + tic(); + double l2 = normL2(_range=elements(mesh), _expr=idv(u)-expr(solution) ); + toc("L2 error norm"); + tic(); + + auto range=elements(mesh); + auto expr=idv(u)-Feel::expr(solution); + auto grad_expr=gradv(u)-grad<2>(Feel::expr(solution)); + + + double h1 = normH1(_range=range, + _expr=expr, + _grad_expr=grad_expr + ); + + + + toc("H1 error norm"); + return { { "L2", l2 }, { "H1", h1 } }; + }; + + + + int status = thechecker.runOnce( norms, rate::hp( mesh->hMax(), Vh->fe()->order() ) ); + + // exit status = 0 means no error + return !status; + + //return 0; + +} + + +