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
residualestimator.hpp
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 Prud'homme <christophe.prudhomme@feelpp.org>
6  Stephane Veys <stephane.veys@gmail.com>
7  Date: 2010-08-05
8 
9  Copyright (C) 2010 Universite Joseph Fourier (Grenoble I)
10 
11  This library is free software; you can redistribute it and/or
12  modify it under the terms of the GNU Lesser General Public
13  License as published by the Free Software Foundation; either
14  version 2.1 of the License, or (at your option) any later version.
15 
16  This library is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19  Lesser General Public License for more details.
20 
21  You should have received a copy of the GNU Lesser General Public
22  License along with this library; if not, write to the Free Software
23  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24 */
30 #include <feel/feel.hpp>
31 #include <feel/feeldiscr/elementdiv.hpp>
32 #include <feel/feelfilters/creategmshmesh.hpp>
33 
34 using namespace Feel;
35 using namespace boost::numeric::ublas;
36 
45 inline
46 AboutData
47 makeAbout()
48 {
49  AboutData about( "residualestimator" ,
50  "residualestimator" ,
51  "0.2",
52  "nD(n=1,2,3) Residual Estimator on Laplacian equation",
53  Feel::AboutData::License_GPL,
54  "Copyright (c) 2010 Universite Joseph Fourier" );
55 
56  about.addAuthor( "Christophe Prud'homme", "developer", "christophe.prudhomme@feelpp.org", "" );
57  about.addAuthor( "Stéphane Veys", "developer", "stephane.veys@gmail.com", "" );
58  return about;
59 
60 }
68 template<int Dim, int Order = 1>
70  :
71 public Simget
72 {
73  typedef Simget super;
74 public:
76  typedef double value_type;
77 
79  typedef Backend<value_type> backend_type;
81  typedef boost::shared_ptr<backend_type> backend_ptrtype;
82 
83 
85  typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
87  typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
89  typedef typename backend_type::vector_type vector_type;
91  typedef typename backend_type::vector_ptrtype vector_ptrtype;
92 
94  typedef Simplex<Dim> convex_type;
96  typedef Mesh<convex_type> mesh_type;
98  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
99 
101  typedef FunctionSpace<mesh_type, bases<Lagrange<0,Scalar,Discontinuous> > > p0_space_type;
102  typedef boost::shared_ptr<p0_space_type> p0_space_ptrtype;
104  typedef typename p0_space_type::element_type p0_element_type;
105 
106 
108  typedef bases<Lagrange<1,Scalar> > p1_basis_type;
109  typedef FunctionSpace<mesh_type, p1_basis_type> p1_space_type;
110  typedef boost::shared_ptr<p1_space_type> p1_space_ptrtype;
111  typedef typename p1_space_type::element_type p1_element_type;
112 
113 
115  typedef bases<Lagrange<Order,Scalar> > basis_type;
116 
118  typedef FunctionSpace<mesh_type, basis_type> space_type;
120  typedef boost::shared_ptr<space_type> space_ptrtype;
122  typedef typename space_type::element_type element_type;
123 
125  typedef Exporter<mesh_type,1> export_type;
127  typedef boost::shared_ptr<export_type> export_ptrtype;
128 
132  ResidualEstimator( AboutData const& about )
133  :
134  super(),
135  M_backend( backend_type::build() ),
136  M_backendP1( backend_type::build() ),
137  meshSize( 0.1 ),
138  exporter( export_type::New( "gmsh", this->about().appName() ) ),
139  order( 1 ),
140  dim( 1 ),
141  shape( "hypercube" ),
142  fn( 1 ),
143  alpha( 3 ),
144  beta( 10 ),
145  weakdir( 1 ),
146  error_type( 1 ),
147  tol( 1e-2 ),
148  penaldir( 50 )
149 
150  {
151  }
153  :
154  super(),
155  M_backend( backend_type::build( this->vm() ) ),
156  M_backendP1( backend_type::build( this->vm() ) ),
157  meshSize( this->vm()["hsize"].template as<double>() ),
158  exporter( export_type::New( "gmsh", this->about().appName() ) ),
159  order( this->vm()["order"].template as<int>() ),
160  dim( this->vm()["dim"].template as<int>() ),
161  shape( this->vm()["shape"].template as<std::string>() ),
162  fn( this->vm()["fn"].template as<int>() ),
163  alpha( this->vm()["alpha"].template as<double>() ),
164  beta( this->vm()["beta"].template as<double>() ),
165  weakdir( this->vm()["weakdir"].template as<int>() ),
166  error_type( this->vm()["adapt-error-type"].template as<int>() ),
167  tol( this->vm()["adapt-tolerance"].template as<double>() ),
168  penaldir( this->vm()["penaldir"].template as<double>() )
169 
170  {
171  }
172 
173  void run();
174 
175  void run( const double* X, unsigned long P, double* Y, unsigned long N );
176 
177  // this function will move to the mesh library in the mesh class
178  BOOST_PARAMETER_CONST_MEMBER_FUNCTION(
179  ( mesh_ptrtype ), // return type
180  adapt, // 2. function name
181 
182  tag, // 3. namespace of tag types
183 
184  ( required
185  ( h, * ) ) // 4. one required parameter, and
186 
187  ( optional
188  ( maxit, *( boost::is_integral<mpl::_> ), 10 )
189  ( hmin, *( boost::is_arithmetic<mpl::_> ), 1e-2 )
190  ( hmax, *( boost::is_arithmetic<mpl::_> ), 2 )
191  ( model, *, "" )
192  ( statistics, *( boost::is_integral<mpl::_> ), 0 )
193  ( update, *( boost::is_integral<mpl::_> ), MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER )
194  ( collapseOnBoundary, *( boost::is_integral<mpl::_> ), true )
195  ( collapseOnBoundaryTolerance, *( boost::is_arithmetic<mpl::_> ), 1e-6 ) ) // 5. optional
196  )
197  {
198 
199  }
200 
201 private:
202 
204  backend_ptrtype M_backend;
205  backend_ptrtype M_backendP1;
206 
208  double meshSize;
209 
211  int dim;
212 
214  int order;
215 
217  std::string shape;
218 
220  int fn;
222  double alpha;
223  double beta;
224  bool weakdir;
225  double penaldir;
226  int error_type;
227  double tol;
228 
229 
231  mesh_ptrtype mesh;
232 
234  export_ptrtype exporter;
235 
237  p0_space_ptrtype P0h;
238  p1_space_ptrtype P1h;
239 
240  //double estimatorH1, estimatorL2, estimator;
241  p1_element_type h_new;
242  std::string msh_name;
243  bool first_time;
244 
245  int tag_Neumann;
246  int tag_Dirichlet;
247 }; // ResidualEstimator
248 
249 template<int Dim, int Order>
250 void
252 {
253  if ( dim && dim != Dim ) return ;
254 
255  if ( order && order != Order ) return ;
256 
257  std::cout << "------------------------------------------------------------\n";
258  std::cout << "Execute ResidualEstimator<" << Dim << ">\n";
259  std::vector<double> X;
260  X.push_back( meshSize );
261 
262  if ( shape == "hypercube" )
263  X.push_back( 1 );
264 
265  else if ( shape == "ellipsoid" )
266  X.push_back( 2 );
267 
268  else // default is simplex
269  X.push_back( 0 );
270 
271  X.push_back( fn );
272  X.push_back( alpha );
273  X.push_back( beta );
274  X.push_back( weakdir );
275  X.push_back( penaldir );
276  first_time=true;
277  X.push_back( first_time );
278  X.push_back( error_type );
279  X.push_back( tol );
280  std::vector<double> Y( 4 );
281  run( X.data(), X.size(), Y.data(), Y.size() );
282 
283 }
284 template<int Dim, int Order>
285 void
286 ResidualEstimator<Dim,Order>::run( const double* X, unsigned long P, double* Y, unsigned long n )
287 {
288  /*
289  * set parameters
290  */
291  meshSize = X[0];
292 
293  if ( X[1] == 0 ) shape = "simplex";
294 
295  if ( X[1] == 1 ) shape = "hypercube";
296 
297  if ( X[1] == 2 ) shape = "ellipsoid";
298 
299  fn = X[2];
300  alpha = X[3];
301  beta = X[4];
302  weakdir = X[5];
303  penaldir = X[6];
304  first_time = X[7];
305  error_type = X[8];
306  tol = X[9];
307 
308  double estimatorH1, estimatorL2, estimator;
309 
310  if ( first_time )
311  {
312  if ( !this->vm().count( "nochdir" ) )
313  Environment::changeRepository( boost::format( "doc/tutorial/%1%/%2%-%3%/P%4%/h_%5%/" )
314  % this->about().appName()
315  % shape
316  % Dim
317  % Order
318  % meshSize );
319 
320  mesh = createGMSHMesh( _mesh=new mesh_type,
321  _parametricnodes=1,
322  _desc=domain( _name=( boost::format( "%1%-%2%" ) % shape % Dim ).str() ,
323  _usenames=true,
324  _shape=shape,
325  _dim=Dim,
326  _h=X[0] ),
327  _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES );
328  tag_Neumann = mesh->markerName( "Neumann" );
329  tag_Dirichlet = mesh->markerName( "Dirichlet" );
330 
331  }//end if(first_time)
332 
333 
338  // [toto1]
339  P0h = p0_space_type::New( mesh );
340  P1h = p1_space_type::New( mesh );
341  space_ptrtype Xh = space_type::New( mesh );
342  element_type u( Xh, "u" );
343  element_type v( Xh, "v" );
344  // [toto1]
345 
346 
347 
353  // [toto2]
354  //# marker1 #
355  value_type pi = M_PI;
357  auto fn1 = ( 1-Px()*Px() )*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() );
358  auto fn2 = sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() );
359  auto g = chi( fn==1 )*fn1 + chi( fn==2 )*fn2;
360  auto grad_g =
361  trans( chi( fn==1 )*( ( -2*Px()*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() )+beta*fn1 )*unitX()+
362  -2*Py()*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*unitY()+
363  -2*Pz()*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*unitZ() )+
364  chi( fn==2 )*( +( alpha*pi*cos( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )+beta*fn2 )*unitX()+
365  -alpha*pi*sin( alpha*pi*Px() )*sin( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )*unitY()+
366  -alpha*pi*sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*sin( alpha*pi*Pz() )*exp( beta*Px() )*unitZ() ) );
368  auto minus_laplacian_g =
369  ( chi( fn == 1 )*( 2*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() ) + 4*beta*Px()*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() ) - beta*beta *fn1 +
370  2*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*chi( Dim >= 2 ) +
371  2*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*chi( Dim == 3 ) ) +
372  chi( fn == 2 )* (
373  exp( beta*Px() )*( Dim*alpha*alpha*pi*pi*sin( alpha*pi*Px() )-beta*beta*sin( alpha*pi*Px() )-2*beta*alpha*pi*cos( alpha*pi*Px() ) )*
374  ( cos( alpha*pi*Py() )*chi( Dim>=2 ) + chi( Dim==1 ) ) * ( cos( alpha*pi*Pz() )*chi( Dim==3 ) + chi( Dim<=2 ) )
375  )
376 
377  );
378 
379  //# endmarker1 #
380  // [toto2]
381 
382 
383  using namespace Feel::vf;
384 
391  // [toto3]
392  //# marker2 #
393  vector_ptrtype F( M_backend->newVector( Xh ) );
394  form1( _test=Xh, _vector=F, _init=true ) =
395  integrate( elements( mesh ), minus_laplacian_g*id( v ) )+
396  integrate( markedfaces( mesh, tag_Neumann ),
397  grad_g*vf::N()*id( v ) );
398 
399  //# endmarker2 #
400  if ( this->comm().size() != 1 || weakdir )
401  {
402  //# marker41 #
403  form1( _test=Xh, _vector=F ) +=
404  integrate( markedfaces( mesh,tag_Dirichlet ), g*( -grad( v )*vf::N()+penaldir*id( v )/hFace() ) );
405  //# endmarker41 #
406  }
407 
408  F->close();
409 
410  // [toto3]
411 
417  //# marker3 #
418  // [toto4]
419  sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );
420  // [toto4]
421 
426  // [toto5]
427  form2( Xh, Xh, D, _init=true ) =
428  integrate( elements( mesh ), gradt( u )*trans( grad( v ) ) );
429  // [toto5]
430  //# endmarker3 #
431 
432  if ( this->comm().size() != 1 || weakdir )
433  {
440  // [toto6]
441  //# marker10 #
442  form2( Xh, Xh, D ) +=
443  integrate( markedfaces( mesh,tag_Dirichlet ),
444  -( gradt( u )*vf::N() )*id( v )
445  -( grad( v )*vf::N() )*idt( u )
446  +penaldir*id( v )*idt( u )/hFace() );
447  D->close();
448  //# endmarker10 #
449  // [toto6]
450  }
451 
452  else
453  {
459  // [toto7]
460  //# marker5 #
461  D->close();
462  form2( Xh, Xh, D ) +=
463  on( markedfaces( mesh, tag_Dirichlet ), u, F, g );
464  //# endmarker5 #
465  // [toto7]
466 
467  }
468 
469 
474  // [toto8]
476  backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );
477  // [toto8]
478 
483  // [toto9]
484  //# marker7 #
485  double L2exact = math::sqrt( integrate( elements( mesh ),g*g ).evaluate()( 0,0 ) );
486  double L2error2 =integrate( elements( mesh ),( idv( u )-g )*( idv( u )-g ) ).evaluate()( 0,0 );
487  double L2error = math::sqrt( L2error2 );
488  // [toto9]
489  double semiH1error2 = integrate( elements( mesh ),( gradv( u )-grad_g )*( gradv( u )-grad_g ) ).evaluate()( 0,0 );
490  double H1error = math::sqrt( L2error2+semiH1error2 );
491  double H1exact = math::sqrt( integrate( elements( mesh ),g*g+grad_g*trans( grad_g ) ).evaluate()( 0,0 ) );
492 
493 
494  auto RealL2ErrorP0 = integrate( elements( mesh ), ( idv( u )-g )*( idv( u )-g ), _Q<10>() ).broken( P0h );
495  auto RealSemiH1ErrorP0 = integrate( elements( mesh ), ( gradv( u )-grad_g )* trans( gradv( u )-grad_g ), _Q<10>() ).broken( P0h );
496  auto H1RealErrorP0 = RealL2ErrorP0;
497  H1RealErrorP0 += RealSemiH1ErrorP0;
498  H1RealErrorP0 = H1RealErrorP0.sqrt();
499  /*******************residual estimator**********************/
500  //the source terme is given by : minus_laplacian_g
501 
502 
503  auto term1 = vf::h()*( minus_laplacian_g+trace( hessv( u ) ) );
504  auto term2 = jumpv( gradv( u ) );
505  auto term3 = gradv( u )*vf::N()-grad_g*N();
506 
507 
508  auto rho = integrate( elements( mesh ), term1*term1, _Q<10>() ).broken( P0h ).sqrt();
509 
510 
511  rho += integrate( internalfaces( mesh ),0.25*vf::h()*term2*term2 ).broken( P0h ).sqrt();
512 
513  rho += integrate( markedfaces( mesh,tag_Neumann ),
514  vf::h()*term3*term3 ).broken( P0h ).sqrt();
515 
516 
517  auto h=vf::project( P0h, elements( mesh ), vf::h() );
518  auto npen=vf::project( P0h, elements( mesh ), vf::nPEN() );
519  auto H1estimator = rho;
520 
521  //if we use real error for adaptation then H1errorP1 will be the projection of H1RealErrorP0 on P1 space
522  //else H1errorP1 will be the projection of H1estimator on P1 space
523  p1_element_type H1errorP1;
524 
525  if ( error_type==2 )
526  {
527  H1errorP1 = div( vf::sum( P1h, idv( H1RealErrorP0 )*meas() ), vf::sum( P1h, meas() ) );
528  }
529 
530  else if ( error_type==1 )
531  {
532  H1errorP1 = div( vf::sum( P1h, idv( H1estimator )*meas() ), vf::sum( P1h, meas() ) );
533  }
534 
535  else
536  {
537  std::cout<<"Problem with parameter adapt-error-type, please choice between 1 and 2"<<std::endl;
538  return;
539  }
540 
541  //auto new_hsize=vf::project( P1h, elements(mesh), vf::pow(vf::pow(vf::h(),Order)*(1e-4)/idv(H1estimatorP1),1./Order));
542  //auto new_hsize=vf::project( P1h, elements(mesh), vf::h()*(1e-2)/idv(H1estimatorP1) );
543 
544  estimatorH1=math::sqrt( H1estimator.pow( 2 ).sum() );
545  estimatorL2=math::sqrt( element_product( H1estimator,h ).pow( 2 ).sum() );
546 
547 
548  Y[0] = L2error/L2exact;
549  Y[1] = H1error/H1exact;
550  Y[2] = estimatorL2/L2exact;
551  Y[3] = estimatorH1/H1exact;
552 
553  h_new = P1h->element();
554 
555  h_new = vf::project( P1h, elements( mesh ),
556  vf::max( vf::pow(
557  vf::pow( vf::h(),Order )*( tol )/idv( H1errorP1 ),
558  1./Order ),
559  this->vm()["adapt-hmin"].template as<double>() ) );
560  /**********************end of residual estimaor*************/
561 
562 
565  element_type exact_solution( Xh, "exact_solution" );
566  exact_solution = vf::project( Xh, elements( mesh ), g );
567  element_type u_minus_exact( Xh, "u-g" );
568  u_minus_exact = vf::project( Xh, elements( mesh ), idv( u )-g );
569 
570 
571  if ( exporter->doExport() )
572  {
573  LOG(INFO) << "exportResults starts\n";
574 
575  exporter->step( 0 )->setMesh( mesh );
576  exporter->step( 0 )->add( "unknown", u );
577  exporter->step( 0 )->add( "exact solution", exact_solution );
578  exporter->step( 0 )->add( "u - exact solution", u_minus_exact ) ;
579  exporter->step( 0 )->add( "nPEN" , npen );
580  exporter->step( 0 )->add( "H1 error estimator P0" , H1estimator );
581  exporter->step( 0 )->add( "H1 Real error P0" , H1RealErrorP0 );
582  exporter->step( 0 )->add( "H1 error P1 (used for determination of new hsize)" , H1errorP1 );
583  exporter->step( 0 )->add( "new hsize" , h_new );
584 
585  exporter->save();
586  LOG(INFO) << "exportResults done\n";
587  }
588 
589 
590  LOG(INFO)<< " real L2 error : "<<Y[0]<<"\n";
591  LOG(INFO)<< " estimated L2 error "<<Y[2]<<"\n";
592  LOG(INFO)<< " real H1 error : "<<Y[1]<<"\n";
593  LOG(INFO)<< " estimated H1 error "<<Y[3]<<"\n";
594 
595 
596  std::ostringstream geostr;
597 
598  if ( this->vm()["gmshmodel"].template as<bool>() )
599  {
600  if ( this->vm()["gmshgeo"].template as<bool>() )
601  geostr << shape << "-" << Dim << ".geo";
602 
603  else
604  geostr << shape << "-" << Dim << ".msh";
605  }
606 
607  mesh = adapt( _h=h_new,
608  _model=geostr.str(),
609  _hmin=this->vm()["adapt-hmin"].template as<double>(),
610  _hmax=this->vm()["adapt-hmax"].template as<double>() );
611 } // ResidualEstimator::run
612 
613 
614 
615 
616 
617 
618 
backend_type::vector_type vector_type
vector type associated with backend
Definition: residualestimator.hpp:89
Backend< value_type > backend_type
linear algebra backend factory
Definition: residualestimator.hpp:79
bases< Lagrange< Order, Scalar > > basis_type
the basis type of our approximation space
Definition: residualestimator.hpp:115
Exporter< mesh_type, 1 > export_type
the exporter factory type
Definition: residualestimator.hpp:125
backend_type::sparse_matrix_type sparse_matrix_type
sparse matrix type associated with backend
Definition: residualestimator.hpp:85
Definition: residualestimator.hpp:69
boost::shared_ptr< export_type > export_ptrtype
the exporter factory (shared_ptr<> type)
Definition: residualestimator.hpp:127
STL namespace.
boost::shared_ptr< space_type > space_ptrtype
the approximation function space type (shared_ptr<> type)
Definition: residualestimator.hpp:120
p0_space_type::element_type p0_element_type
an element type of the discontinuous function space
Definition: residualestimator.hpp:104
double value_type
numerical type is double
Definition: residualestimator.hpp:76
space_type::element_type element_type
an element type of the approximation function space
Definition: residualestimator.hpp:122
ResidualEstimator(AboutData const &about)
Definition: residualestimator.hpp:132
FunctionSpace< mesh_type, basis_type > space_type
the approximation function space type
Definition: residualestimator.hpp:118
Simplex< Dim > convex_type
geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1 ...
Definition: residualestimator.hpp:94
backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype
sparse matrix type associated with backend (shared_ptr<> type)
Definition: residualestimator.hpp:87
Mesh< convex_type > mesh_type
mesh type
Definition: residualestimator.hpp:96
FunctionSpace< mesh_type, bases< Lagrange< 0, Scalar, Discontinuous > > > p0_space_type
function space that holds piecewise constant ( ) functions (e.g. to store material properties or part...
Definition: residualestimator.hpp:101
bases< Lagrange< 1, Scalar > > p1_basis_type
the basis type of our approximation space
Definition: residualestimator.hpp:108
backend_type::vector_ptrtype vector_ptrtype
vector type associated with backend (shared_ptr<> type)
Definition: residualestimator.hpp:91
Definition: qs_laplacian.doc:2
boost::shared_ptr< backend_type > backend_ptrtype
linear algebra backend factory shared_ptr<> type
Definition: residualestimator.hpp:81
boost::shared_ptr< mesh_type > mesh_ptrtype
mesh shared_ptr<> type
Definition: residualestimator.hpp:98