35 #include <boost/range/algorithm/for_each.hpp>
36 #include <boost/algorithm/string/split.hpp>
38 #include <boost/parameter.hpp>
39 #include <feel/feelcore/parameter.hpp>
41 #include <feel/feel.hpp>
67 BOOST_PARAMETER_NAME(Xh)
68 BOOST_PARAMETER_NAME(T)
69 BOOST_PARAMETER_NAME(exact_solution)
70 BOOST_PARAMETER_NAME(parameters)
71 BOOST_PARAMETER_NAME(convergence)
72 BOOST_PARAMETER_NAME(convergence_steps)
81 template<
int Dim,
int Order>
105 ErrorBase( po::variables_map
const& vm, std::string
const& prefix)
107 M_exact( vm[prefixvm( prefix,
"error.exact" )].as<std::string>() ),
108 M_exact_params( vm[prefixvm( prefix,
"error.params" )].as<std::string>() ),
109 M_rhs( vm[prefixvm( prefix,
"error.rhs" )].as<std::string>() ),
110 M_rhs_computed(vm[prefixvm( prefix,
"error.rhs.computed" )].as<bool>() ),
111 M_convergence(vm[prefixvm( prefix,
"error.convergence" )].as<bool>() ),
112 M_convergence_max(vm[prefixvm( prefix,
"error.convergence.steps" )].as<int>() )
116 ErrorBase( std::string
const& name, std::string
const& params)
119 M_exact_params(params),
121 M_rhs_computed( false ),
127 ErrorBase( ErrorBase
const &e )
129 M_exact( e.M_exact ),
130 M_exact_params( e.M_exact_params ),
132 M_rhs_computed( e.M_rhs_computed ),
133 M_convergence( e.M_convergence ),
134 M_convergence_max( e.M_convergence_max )
138 ErrorBase&
operator=( ErrorBase
const& e )
142 M_exact = e.M_exact ;
143 M_exact_params = e.M_exact_params;
145 M_rhs_computed = e.M_rhs_computed;
146 M_convergence = e.M_convergence;
147 M_convergence_max = e.M_convergence_max;
153 virtual ~ErrorBase() {}
160 void setParams ( std::string params )
163 M_exact_params=params;
169 M_rhs_computed=doComputeRhs;
174 M_convergence=doConvergence;
183 void setSolution(std::string
const& expression =
"", std::string
const& p =
"")
185 if ( expression.empty() )
186 FEELPP_ASSERT( M_rhs.size() ).error(
"undefined rhs exact" );
188 M_exact = expression;
193 vars = symbols<Dim>();
194 std::vector<std::string> lst_params;
195 boost::split(lst_params, M_exact_params, boost::is_any_of(
";"));
196 LOG(INFO) <<
"Loading additionnal parameters : ";
197 boost::for_each(lst_params, [](
const std::string &s){LOG(INFO) << s <<
", ";});
198 LOG(INFO) << std::endl;
199 parameters = symbols(lst_params);
201 LOG(INFO) <<
"Loading function : " << M_exact << std::endl;
202 boost::for_each( parameters, [](symbol
const& s ) {LOG(INFO) <<
"Symbol " << s.get_name() <<
" found\n";});
203 exact = parse(M_exact, vars, parameters);
204 LOG(INFO) <<
"exact solution is : " << exact <<
"\n";
207 void setRhs(std::string
const& expression =
"")
209 if ( expression.empty() )
210 FEELPP_ASSERT( M_rhs.size() ).error(
"undefined rhs exact" );
214 LOG(INFO) <<
"Loading rhs function : " << M_rhs << std::endl;
215 rhs = parse(M_rhs, vars, parameters);
216 LOG(INFO) <<
"rhs is : " << rhs <<
"\n";
219 typedef typename boost::function<ex (ex u, std::vector<symbol> vars, std::vector<symbol> p)> t_edp_type;
220 typedef typename boost::function<ex (ex u, std::vector<symbol> vars)> edp_type;
222 void setRhs(t_edp_type * edp)
224 FEELPP_ASSERT( parameters.size() ).error(
"setRhs: inconsistant numbers of parameters" );
225 rhs = (*edp)(exact, vars, parameters);
226 LOG(INFO) <<
"computed rhs is : " << rhs <<
"\n";
228 std::ostringstream rhs_expression;
229 rhs_expression << rhs;
230 M_rhs = rhs_expression.str();
233 void setRhs(edp_type * edp)
235 rhs = (*edp)(exact, vars);
236 LOG(INFO) <<
"computed rhs is : " << rhs <<
"\n";
238 std::ostringstream rhs_expression;
239 rhs_expression << rhs;
240 M_rhs = rhs_expression.str();
253 return M_rhs_computed;
261 return M_convergence;
269 return M_convergence_max;
311 FEELPP_ASSERT( M_rhs.size() ).error(
"undefined rhs" );
321 FEELPP_ASSERT( M_rhs.size() ).error(
"undefined rhs" );
324 FEELPP_ASSERT( parameters.size() == 1 ).error(
"inconsistant values and parameters size" );
327 tmp_sol = substitute(tmp_sol, parameters[0], value);
335 ex
getRhs(
const std::vector<double> & values)
const
338 FEELPP_ASSERT( M_rhs.size() ).error(
"undefined exact solution" );
341 FEELPP_ASSERT( values.size() == parameters.size() ).error(
"inconsistant values and parameters size" );
344 for (
unsigned int i=0; i<values.size(); i++)
345 tmp_sol = substitute(tmp_sol, parameters[i], values[i]);
356 FEELPP_ASSERT( M_exact.size() ).error(
"undefined exact solution" );
366 FEELPP_ASSERT( M_exact.size() ).error(
"undefined exact solution" );
369 FEELPP_ASSERT( parameters.size() == 1 ).error(
"inconsistant values and parameters size" );
372 tmp_sol = substitute(tmp_sol, parameters[0], value);
383 FEELPP_ASSERT( M_exact.size() ).error(
"undefined exact solution" );
386 FEELPP_ASSERT( values.size() == parameters.size() ).error(
"inconsistant values and parameters size" );
389 for (
unsigned int i=0; i<values.size(); i++)
390 tmp_sol = substitute(tmp_sol, parameters[i], values[i]);
398 LOG(INFO) <<
"============================================================\n";
399 LOG(INFO) <<
"Error Information\n";
400 LOG(INFO) <<
" exact : " << getSolution() <<
"\n";
401 LOG(INFO) <<
" params : ";
402 boost::for_each( parameters, [](symbol
const& s ) {LOG(INFO) << s.get_name() <<
" ";});
404 if ( !M_rhs.empty() )
406 LOG(INFO) <<
" rhs : " << getRhs() <<
"\n";
408 LOG(INFO) <<
" convergence : " << convergence() <<
"\n";
409 LOG(INFO) <<
" convergence steps : " << numberOfConvergenceSteps() <<
"\n";
413 element_type rhs_project(
const space_ptrtype & Xh)
415 ex solution = getRhs();
417 auto mesh = Xh->mesh();
422 element_type exact_project(
const space_ptrtype & Xh)
424 ex solution = getSolution();
426 auto mesh = Xh->mesh();
431 element_type grad_exact_project(
const space_ptrtype & Xh)
433 ex solution = getSolution();
434 auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
436 auto mesh = Xh->mesh();
441 element_type error_project(
const space_ptrtype & Xh,
const element_type & T)
443 auto gproj = exact_project(Xh);
447 element_type grad_error_project(
const space_ptrtype & Xh,
const element_type & T)
449 auto gproj = grad_exact_project(Xh);
450 return (gproj - gradv(T));
453 element_type rhs_project(
const space_ptrtype & Xh,
const double &
val)
455 ex solution = getRhs(val);
457 auto mesh = Xh->mesh();
462 element_type exact_project(
const space_ptrtype & Xh,
const double & val)
464 ex solution = getSolution(val);
466 auto mesh = Xh->mesh();
471 element_type grad_exact_project(
const space_ptrtype & Xh,
const double & val)
473 ex solution = getSolution(val);
474 auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
476 auto mesh = Xh->mesh();
481 element_type error_project(
const space_ptrtype & Xh,
const element_type & T,
const double & val)
483 auto gproj = exact_project(Xh, val);
487 element_type grad_error_project(
const space_ptrtype & Xh,
const element_type & T,
const double & val)
489 auto gproj = grad_exact_project(Xh, val);
490 return (gproj - gradv(T));
495 double L2_error(
const space_ptrtype & Xh,
const element_type & T)
const
498 ex solution = getSolution();
500 auto g = expr(solution,vars);
501 auto mesh = Xh->mesh();
502 return math::sqrt(
integrate(
elements(mesh), Px()*(idv(T) - g)*(idv(T) - g) ).evaluate()(0,0) );
505 double H1_error(
const space_ptrtype & Xh,
const element_type & T)
const
508 ex solution = getSolution();
510 auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
511 auto mesh = Xh->mesh();
513 double L2error = L2_error(Xh, T);
514 double H1seminorm = math::sqrt(
integrate(
elements(mesh), Px()*(gradv(T) - gradg)*trans(gradv(T) - gradg) ).evaluate()(0,0) );
515 return math::sqrt( L2error*L2error + H1seminorm*H1seminorm);
518 double L2_error(
const space_ptrtype & Xh,
const element_type & T,
const double & values)
const
521 ex solution = getSolution(values);
523 auto g = expr(solution,vars);
524 auto mesh = Xh->mesh();
525 return math::sqrt(
integrate(
elements(mesh), Px()*(idv(T) - g)*(idv(T) - g) ).evaluate()(0,0) );
528 double H1_error(
const space_ptrtype & Xh,
const element_type & T,
const double & values)
const
531 ex solution = getSolution(values);
533 auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
534 auto mesh = Xh->mesh();
536 double L2error = L2_error(Xh, T, values);
537 double H1seminorm = math::sqrt(
integrate(
elements(mesh), Px()*(gradv(T) - gradg)*trans(gradv(T) - gradg) ).evaluate()(0,0) );
538 return math::sqrt( L2error*L2error + H1seminorm*H1seminorm);
556 int M_convergence_max;
560 std::vector<symbol> vars;
561 std::vector<symbol> parameters;