Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Files Functions Variables Typedefs Pages
error.hpp
Go to the documentation of this file.
1 /* -*- 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=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Trophime <christophe.trophime@lncmi.cnrs.fr>
6  Date: 2012-06-02
7 
8  Copyright (C) 2008-2012 Universite Joseph Fourier (Grenoble I)
9 
10  This program is free software: you can redistribute it and/or modify
11  it under the terms of the GNU General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  This program is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  GNU General Public License for more details.
19 
20  You should have received a copy of the GNU General Public License
21  along with this program. If not, see <http://www.gnu.org/licenses/>.
22 */
29 #ifndef __ERROR_HPP
30 #define __ERROR_HPP 1
31 
32 #include <string>
33 #include <sstream>
34 
35 #include <boost/range/algorithm/for_each.hpp>
36 #include <boost/algorithm/string/split.hpp>
37 
38 #include <boost/parameter.hpp>
39 #include <feel/feelcore/parameter.hpp>
40 
41 #include <feel/feel.hpp>
42 
43 // cf bdf class for instance :
44 
45 // should add default option
46 // params
47 // exact
48 // rhs
49 // convergence_study
50 
51 // should
52 // compute error(s) (both norm and "error")
53 // run cvg study
54 // display result (gnuplot_iostream?)
55 
56 // should
57 // run external code
58 // read output from code
59 // display some error indicator
60 
61 namespace Feel
62 {
63  using GiNaC::symbol;
64  using GiNaC::lst;
65  using GiNaC::ex;
66 
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)
73 
81  template<int Dim, int Order>
82  class ErrorBase
83  {
84  public:
85 
87  typedef Simplex<Dim> convex_type;
89  typedef Mesh<convex_type> mesh_type;
91  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
92 
94  typedef bases<Lagrange<Order,Scalar> > basis_type;
96  typedef FunctionSpace<mesh_type, basis_type> space_type;
98  typedef boost::shared_ptr<space_type> space_ptrtype;
100  typedef typename space_type::element_type element_type;
101 
105  ErrorBase( po::variables_map const& vm, std::string const& prefix)
106  :
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>() )
113  {
114  }
115 
116  ErrorBase( std::string const& name, std::string const& params)
117  :
118  M_exact( name ),
119  M_exact_params(params),
120  M_rhs(""),
121  M_rhs_computed( false ),
122  M_convergence(0),
123  M_convergence_max(0)
124  {
125  }
126 
127  ErrorBase( ErrorBase const &e )
128  :
129  M_exact( e.M_exact ),
130  M_exact_params( e.M_exact_params ),
131  M_rhs( e.M_rhs ),
132  M_rhs_computed( e.M_rhs_computed ),
133  M_convergence( e.M_convergence ),
134  M_convergence_max( e.M_convergence_max )
135  {
136  }
137 
138  ErrorBase& operator=( ErrorBase const& e )
139  {
140  if ( this != &e )
141  {
142  M_exact = e.M_exact ;
143  M_exact_params = e.M_exact_params;
144  M_rhs = e.M_rhs;
145  M_rhs_computed = e.M_rhs_computed;
146  M_convergence = e.M_convergence;
147  M_convergence_max = e.M_convergence_max;
148  }
149 
150  return *this;
151  }
152 
153  virtual ~ErrorBase() {}
155 
156 
160  void setParams ( std::string params )
162  {
163  M_exact_params=params;
164  }
165 
167  void setComputedRhs ( bool doComputeRhs )
168  {
169  M_rhs_computed=doComputeRhs;
170  }
172  void setConvergence ( bool doConvergence )
173  {
174  M_convergence=doConvergence;
175  }
176 
179  {
180  M_convergence_max=n;
181  }
182 
183  void setSolution(std::string const& expression ="", std::string const& p ="")
184  {
185  if ( expression.empty() )
186  FEELPP_ASSERT( M_rhs.size() ).error( "undefined rhs exact" );
187  else
188  M_exact = expression;
189 
190  if ( !p.empty() )
191  M_exact_params = p;
192 
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);
200 
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";
205  }
206 
207  void setRhs(std::string const& expression ="")
208  {
209  if ( expression.empty() )
210  FEELPP_ASSERT( M_rhs.size() ).error( "undefined rhs exact" );
211  else
212  M_rhs = expression;
213 
214  LOG(INFO) << "Loading rhs function : " << M_rhs << std::endl;
215  rhs = parse(M_rhs, vars, parameters);
216  LOG(INFO) << "rhs is : " << rhs << "\n";
217  }
218 
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;
221 
222  void setRhs(t_edp_type * edp)
223  {
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";
227 
228  std::ostringstream rhs_expression;
229  rhs_expression << rhs;
230  M_rhs = rhs_expression.str();
231  }
232 
233  void setRhs(edp_type * edp)
234  {
235  rhs = (*edp)(exact, vars);
236  LOG(INFO) << "computed rhs is : " << rhs << "\n";
237 
238  std::ostringstream rhs_expression;
239  rhs_expression << rhs;
240  M_rhs = rhs_expression.str();
241  }
243 
247 
251  bool computedrhs() const
252  {
253  return M_rhs_computed;
254  }
255 
259  bool convergence() const
260  {
261  return M_convergence;
262  }
263 
268  {
269  return M_convergence_max;
270  }
271 
275  std::vector<symbol> getVars() const
276  {
277  return vars;
278  }
279 
283  std::vector<symbol> getParams() const
284  {
285  return parameters;
286  }
287 
288 
292  std::string getExactSolution() const
293  {
294  return M_exact;
295  }
296 
300  std::string getExactRhs() const
301  {
302  return M_rhs;
303  }
304 
308  ex getRhs() const
309  {
310  // check if exact is defined
311  FEELPP_ASSERT( M_rhs.size() ).error( "undefined rhs" );
312  return rhs;
313  }
314 
318  ex getRhs(const double & value) const
319  {
320  // check if exact is defined
321  FEELPP_ASSERT( M_rhs.size() ).error( "undefined rhs" );
322 
323  // check if values is consistant with params
324  FEELPP_ASSERT( parameters.size() == 1 ).error( "inconsistant values and parameters size" );
325 
326  ex tmp_sol = rhs;
327  tmp_sol = substitute(tmp_sol, parameters[0], value);
328 
329  return tmp_sol;
330  }
331 
335  ex getRhs(const std::vector<double> & values) const
336  {
337  // check if exact is defined
338  FEELPP_ASSERT( M_rhs.size() ).error( "undefined exact solution" );
339 
340  // check if values is consistant with params
341  FEELPP_ASSERT( values.size() == parameters.size() ).error( "inconsistant values and parameters size" );
342 
343  ex tmp_sol = rhs;
344  for (unsigned int i=0; i<values.size(); i++) // may use boost zip iterator
345  tmp_sol = substitute(tmp_sol, parameters[i], values[i]);
346 
347  return tmp_sol;
348  }
349 
353  ex getSolution() const
354  {
355  // check if exact is defined
356  FEELPP_ASSERT( M_exact.size() ).error( "undefined exact solution" );
357  return exact;
358  }
359 
363  ex getSolution(const double & value) const
364  {
365  // check if exact is defined
366  FEELPP_ASSERT( M_exact.size() ).error( "undefined exact solution" );
367 
368  // check if values is consistant with params
369  FEELPP_ASSERT( parameters.size() == 1 ).error( "inconsistant values and parameters size" );
370 
371  ex tmp_sol = exact;
372  tmp_sol = substitute(tmp_sol, parameters[0], value);
373 
374  return tmp_sol;
375  }
376 
380  ex getSolution(const std::vector<double> & values) const
381  {
382  // check if exact is defined
383  FEELPP_ASSERT( M_exact.size() ).error( "undefined exact solution" );
384 
385  // check if values is consistant with params
386  FEELPP_ASSERT( values.size() == parameters.size() ).error( "inconsistant values and parameters size" );
387 
388  ex tmp_sol = exact;
389  for (unsigned int i=0; i<values.size(); i++) // may use boost zip iterator
390  tmp_sol = substitute(tmp_sol, parameters[i], values[i]);
391 
392  return tmp_sol;
393  }
395 
396  void print() const
397  {
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() << " ";});
403  LOG(INFO) << "\n";
404  if ( !M_rhs.empty() )
405  {
406  LOG(INFO) << " rhs : " << getRhs() << "\n";
407  }
408  LOG(INFO) << " convergence : " << convergence() << "\n";
409  LOG(INFO) << " convergence steps : " << numberOfConvergenceSteps() << "\n";
410  }
411 
412  /****************/
413  element_type rhs_project(const space_ptrtype & Xh)
414  {
415  ex solution = getRhs();
416 
417  auto mesh = Xh->mesh();
418  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(solution,vars) );
419  return gproj;
420  }
421 
422  element_type exact_project(const space_ptrtype & Xh)
423  {
424  ex solution = getSolution();
425 
426  auto mesh = Xh->mesh();
427  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(solution,vars) );
428  return gproj;
429  }
430 
431  element_type grad_exact_project(const space_ptrtype & Xh)
432  {
433  ex solution = getSolution();
434  auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
435 
436  auto mesh = Xh->mesh();
437  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(gradg,vars) );
438  return gproj;
439  }
440 
441  element_type error_project(const space_ptrtype & Xh, const element_type & T)
442  {
443  auto gproj = exact_project(Xh);
444  return (gproj - T);
445  }
446 
447  element_type grad_error_project(const space_ptrtype & Xh, const element_type & T)
448  {
449  auto gproj = grad_exact_project(Xh);
450  return (gproj - gradv(T));
451  }
452 
453  element_type rhs_project(const space_ptrtype & Xh, const double & val)
454  {
455  ex solution = getRhs(val);
456 
457  auto mesh = Xh->mesh();
458  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(solution,vars) );
459  return gproj;
460  }
461 
462  element_type exact_project(const space_ptrtype & Xh, const double & val)
463  {
464  ex solution = getSolution(val);
465 
466  auto mesh = Xh->mesh();
467  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(solution,vars) );
468  return gproj;
469  }
470 
471  element_type grad_exact_project(const space_ptrtype & Xh, const double & val)
472  {
473  ex solution = getSolution(val);
474  auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
475 
476  auto mesh = Xh->mesh();
477  auto gproj = vf::project( _space=Xh, _range=elements( mesh ), _expr=expr(gradg,vars) );
478  return gproj;
479  }
480 
481  element_type error_project(const space_ptrtype & Xh, const element_type & T, const double & val)
482  {
483  auto gproj = exact_project(Xh, val);
484  return (gproj - T);
485  }
486 
487  element_type grad_error_project(const space_ptrtype & Xh, const element_type & T, const double & val)
488  {
489  auto gproj = grad_exact_project(Xh, val);
490  return (gproj - gradv(T));
491  }
492  /****************/
493 
494  // may add optionnal weight for AXI???
495  double L2_error(const space_ptrtype & Xh, const element_type & T) const
496  {
497  // replace params by specified values
498  ex solution = getSolution();
499 
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) );
503  }
504 
505  double H1_error(const space_ptrtype & Xh, const element_type & T) const
506  {
507  // replace params by specified values
508  ex solution = getSolution();
509 
510  auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
511  auto mesh = Xh->mesh();
512 
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);
516  }
517 
518  double L2_error(const space_ptrtype & Xh, const element_type & T, const double & values) const
519  {
520  // replace params by specified values
521  ex solution = getSolution(values);
522 
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) );
526  }
527 
528  double H1_error(const space_ptrtype & Xh, const element_type & T, const double & values) const
529  {
530  // replace params by specified values
531  ex solution = getSolution(values);
532 
533  auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
534  auto mesh = Xh->mesh();
535 
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);
539  }
540 
541 
542  protected:
543 
545  std::string M_exact;
546 
548  std::string M_exact_params;
549 
551  std::string M_rhs;
552  bool M_rhs_computed;
553 
554  // convergence study
555  bool M_convergence;
556  int M_convergence_max;
557 
558  ex exact;
559  ex rhs;
560  std::vector<symbol> vars;
561  std::vector<symbol> parameters;
562 
563 
564  };
565 
566 
567 }
568 #endif /* __ERRROR_HPP */
space_type::element_type element_type
an element type of the approximation function space
Definition: error.hpp:100
Mesh< convex_type > mesh_type
mesh type
Definition: error.hpp:89
Definition: error.hpp:82
std::vector< symbol > getVars() const
Definition: error.hpp:275
ex getRhs(const std::vector< double > &values) const
Definition: error.hpp:335
bool convergence() const
Definition: error.hpp:259
ex getSolution() const
Definition: error.hpp:353
boost::shared_ptr< mesh_type > mesh_ptrtype
mesh shared_ptr<> type
Definition: error.hpp:91
void setConvergence(bool doConvergence)
set the convergence flag
Definition: error.hpp:172
std::string getExactRhs() const
Definition: error.hpp:300
FunctionSpace< mesh_type, basis_type > space_type
the approximation function space type
Definition: error.hpp:96
bases< Lagrange< Order, Scalar > > basis_type
the basis type of our approximation space
Definition: error.hpp:94
std::string M_rhs
name of the rhs
Definition: error.hpp:551
std::string M_exact_params
name of the exact solution parameters
Definition: error.hpp:548
void setComputedRhs(bool doComputeRhs)
set the rhs_computed flag
Definition: error.hpp:167
std::string M_exact
name of the exact solution
Definition: error.hpp:545
ex getSolution(const double &value) const
Definition: error.hpp:363
std::string getExactSolution() const
Definition: error.hpp:292
Simplex< Dim > convex_type
geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1 ...
Definition: error.hpp:87
int numberOfConvergenceSteps() const
Definition: error.hpp:267
std::vector< symbol > getParams() const
Definition: error.hpp:283
boost::shared_ptr< space_type > space_ptrtype
the approximation function space type (shared_ptr<> type)
Definition: error.hpp:98
ex getRhs(const double &value) const
Definition: error.hpp:318
bool computedrhs() const
Definition: error.hpp:251
void setnumberOfConvergenceSteps(int n)
set the convergence interations
Definition: error.hpp:178
ex getSolution(const std::vector< double > &values) const
Definition: error.hpp:380
ex getRhs() const
Definition: error.hpp:308