Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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  Stéphane Veys <stephane.veys@gmail.com>
7  Date: 2010-08-05
8 
9  Copyright (C) 2010 Université 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 
32 #if defined (FEELPP_HAS_MADLIB_H)
33 #include <MAdLib.h>
34 #endif // FEELPP_HAS_MADLIB_H
35 
36 
37 
38 using namespace Feel;
39 using namespace boost::numeric::ublas;
40 
49 inline
51 makeAbout()
52 {
53  AboutData about( "residualestimator" ,
54  "residualestimator" ,
55  "0.2",
56  "nD(n=1,2,3) Residual Estimator on Laplacian equation",
57  Feel::AboutData::License_GPL,
58  "Copyright (c) 2010 Universite Joseph Fourier" );
59 
60  about.addAuthor( "Christophe Prud'homme", "developer", "christophe.prudhomme@feelpp.org", "" );
61  about.addAuthor( "Stéphane Veys", "developer", "stephane.veys@gmail.com", "" );
62  return about;
63 
64 }
72 template<int Dim, int Order = 1>
74  :
75 public Simget
76 {
77  typedef Simget super;
78 public:
80  typedef double value_type;
81 
85  typedef boost::shared_ptr<backend_type> backend_ptrtype;
86 
87 
91  typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
95  typedef typename backend_type::vector_ptrtype vector_ptrtype;
96 
102  typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
103 
106  typedef boost::shared_ptr<p0_space_type> p0_space_ptrtype;
108  typedef typename p0_space_type::element_type p0_element_type;
109 
110 
112  typedef bases<Lagrange<1,Scalar> > p1_basis_type;
114  typedef boost::shared_ptr<p1_space_type> p1_space_ptrtype;
115  typedef typename p1_space_type::element_type p1_element_type;
116 
117 
119  typedef bases<Lagrange<Order,Scalar> > basis_type;
120 
124  typedef boost::shared_ptr<space_type> space_ptrtype;
126  typedef typename space_type::element_type element_type;
127 
131  typedef boost::shared_ptr<export_type> export_ptrtype;
132 
136  ResidualEstimator( AboutData const& about )
137  :
138  super(),
139  M_backend( backend_type::build() ),
140  M_backendP1( backend_type::build() ),
141  meshSize( 0.1 ),
142  exporter( export_type::New( "gmsh", this->about().appName() ) ),
143  order( 1 ),
144  dim( 1 ),
145  shape( "hypercube" ),
146  fn( 1 ),
147  alpha( 3 ),
148  beta( 10 ),
149  weakdir( 1 ),
150  error_type( 1 ),
151  tol( 1e-2 ),
152  penaldir( 50 )
153 
154  {
155  }
157  :
158  super(),
159  M_backend( backend_type::build( this->vm() ) ),
160  M_backendP1( backend_type::build( this->vm() ) ),
161  meshSize( this->vm()["hsize"].template as<double>() ),
162  exporter( export_type::New( "gmsh", this->about().appName() ) ),
163  order( this->vm()["order"].template as<int>() ),
164  dim( this->vm()["dim"].template as<int>() ),
165  shape( this->vm()["shape"].template as<std::string>() ),
166  fn( this->vm()["fn"].template as<int>() ),
167  alpha( this->vm()["alpha"].template as<double>() ),
168  beta( this->vm()["beta"].template as<double>() ),
169  weakdir( this->vm()["weakdir"].template as<int>() ),
170  error_type( this->vm()["adapt-error-type"].template as<int>() ),
171  tol( this->vm()["adapt-tolerance"].template as<double>() ),
172  penaldir( this->vm()["penaldir"].template as<double>() )
173 
174  {
175  }
176 
177  void run();
178 
179  void run( const double* X, unsigned long P, double* Y, unsigned long N );
180 
181  // this function will move to the mesh library in the mesh class
182  BOOST_PARAMETER_CONST_MEMBER_FUNCTION(
183  ( mesh_ptrtype ), // return type
184  adapt, // 2. function name
185 
186  tag, // 3. namespace of tag types
187 
188  ( required
189  ( h, * ) ) // 4. one required parameter, and
190 
191  ( optional
192  ( maxit, *( boost::is_integral<mpl::_> ), 10 )
193  ( hmin, *( boost::is_arithmetic<mpl::_> ), 1e-2 )
194  ( hmax, *( boost::is_arithmetic<mpl::_> ), 2 )
195  ( model, *, "" )
196  ( statistics, *( boost::is_integral<mpl::_> ), 0 )
197  ( update, *( boost::is_integral<mpl::_> ), MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES|MESH_RENUMBER )
198  ( collapseOnBoundary, *( boost::is_integral<mpl::_> ), true )
199  ( collapseOnBoundaryTolerance, *( boost::is_arithmetic<mpl::_> ), 1e-6 ) ) // 5. optional
200  )
201  {
202 #if defined (FEELPP_HAS_MADLIB_H)
203  saveGMSHMesh( _filename="inputmesh.msh",
204  _parametricnodes=!model.empty(),
205  _mesh=h.mesh() );
206  MAd::pGModel themodel = 0;
207  GM_create( &themodel,"theModel" );
208 
209  if ( !model.empty() )
210  GM_read( themodel, model );
211 
212  else
213  GM_read( themodel, "inputmesh.msh" );
214 
215  MAd::pMesh amesh = MAd::M_new( themodel );
216  MAd::M_load( amesh, "inputmesh.msh" );
217 
218  MAd::PWLSField * sizeField = new MAd::PWLSField( amesh );
219  sizeField->setCurrentSize();
220  auto _elit = h.mesh()->beginElement();
221  auto _elen = h.mesh()->endElement();
222 
223  for ( ; _elit != _elen; ++_elit )
224  {
225  for ( int l = 0; l < _elit->numPoints; ++l )
226  {
227  int dof = h.functionSpace()->dof()->localToGlobal( _elit->id(), l, 0 ).get<0>();
228  int pid = _elit->point( l ).id()+1;
229  sizeField->setSize( pid , h( dof ) );
230  }
231  }
232 
233  MAd::MeshAdapter* ma = new MAd::MeshAdapter( amesh,sizeField );
234 
235  ma->setMaxIterationsNumber( maxit );
236  ma->setEdgeLenSqBounds( 1.0/3.0, 3.0 );
237  ma->setNoSwapQuality( 0.1 );
238  ma->setSliverQuality( 0.02 );
239  ma->setSliverPermissionInESplit( true, 10. );
240  ma->setSliverPermissionInECollapse( true, 0.1 );
241  ma->snapVertices();
242  //ma->setEdgeLenSqBounds( hmin*hmin, hmax*hmax );
243 #if 1
244  ma->setGeoTracking( !model.empty() );
245 #if 0
246  true,
247  0,
248  1.,
249  false,
250  false );
251 #endif
252 #endif
253 
254  if ( statistics )
255  {
256  std::cout << "Statistics before optimization: \n";
257  ma->printStatistics( std::cout );
258  ma->writePos( "meanRatioBefore.pos",MAd::OD_MEANRATIO );
259  }
260 
261  // Optimize
262  // ---------
263  ma->run();
264 
265  if ( statistics )
266  {
267  // Outputs final mesh
268  // -------------------
269  std::cout << "Statistics after optimization: \n";
270  ma->printStatistics( std::cout );
271  ma->writePos( "meanRatioAfter.pos",MAd::OD_MEANRATIO );
272  }
273  MAd::M_writeMsh ( amesh, "result.msh", 2, NULL );
274 
275  return loadGMSHMesh( _mesh=new mesh_type,
276  _filename="result.msh",
277  _update=update );
278 
279 #endif // FEELPP_HAS_MADLIB_H
280 
281  }
282 
283 private:
284 
286  backend_ptrtype M_backend;
287  backend_ptrtype M_backendP1;
288 
290  double meshSize;
291 
293  int dim;
294 
296  int order;
297 
299  std::string shape;
300 
302  int fn;
304  double alpha;
305  double beta;
306  bool weakdir;
307  double penaldir;
308  int error_type;
309  double tol;
310 
311 
313  mesh_ptrtype mesh;
314 
316  export_ptrtype exporter;
317 
319  p0_space_ptrtype P0h;
320  p1_space_ptrtype P1h;
321 
322  //double estimatorH1, estimatorL2, estimator;
323  p1_element_type h_new;
324  std::string msh_name;
325  bool first_time;
326 
327  int tag_Neumann;
328  int tag_Dirichlet;
329 }; // ResidualEstimator
330 
331 template<int Dim, int Order>
332 void
334 {
335  if ( dim && dim != Dim ) return ;
336 
337  if ( order && order != Order ) return ;
338 
339  std::cout << "------------------------------------------------------------\n";
340  std::cout << "Execute ResidualEstimator<" << Dim << ">\n";
341  std::vector<double> X;
342  X.push_back( meshSize );
343 
344  if ( shape == "hypercube" )
345  X.push_back( 1 );
346 
347  else if ( shape == "ellipsoid" )
348  X.push_back( 2 );
349 
350  else // default is simplex
351  X.push_back( 0 );
352 
353  X.push_back( fn );
354  X.push_back( alpha );
355  X.push_back( beta );
356  X.push_back( weakdir );
357  X.push_back( penaldir );
358  first_time=true;
359  X.push_back( first_time );
360  X.push_back( error_type );
361  X.push_back( tol );
362  std::vector<double> Y( 4 );
363  run( X.data(), X.size(), Y.data(), Y.size() );
364 
365 }
366 template<int Dim, int Order>
367 void
368 ResidualEstimator<Dim,Order>::run( const double* X, unsigned long P, double* Y, unsigned long n )
369 {
370  /*
371  * set parameters
372  */
373  meshSize = X[0];
374 
375  if ( X[1] == 0 ) shape = "simplex";
376 
377  if ( X[1] == 1 ) shape = "hypercube";
378 
379  if ( X[1] == 2 ) shape = "ellipsoid";
380 
381  fn = X[2];
382  alpha = X[3];
383  beta = X[4];
384  weakdir = X[5];
385  penaldir = X[6];
386  first_time = X[7];
387  error_type = X[8];
388  tol = X[9];
389 
390  double estimatorH1, estimatorL2, estimator;
391 
392  if ( first_time )
393  {
394  if ( !this->vm().count( "nochdir" ) )
395  Environment::changeRepository( boost::format( "doc/tutorial/%1%/%2%-%3%/P%4%/h_%5%/" )
396  % this->about().appName()
397  % shape
398  % Dim
399  % Order
400  % meshSize );
401 
402  mesh = createGMSHMesh( _mesh=new mesh_type,
403  _parametricnodes=1,
404  _desc=domain( _name=( boost::format( "%1%-%2%" ) % shape % Dim ).str() ,
405  _usenames=true,
406  _shape=shape,
407  _dim=Dim,
408  _h=X[0] ),
409  _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES );
410  tag_Neumann = mesh->markerName( "Neumann" );
411  tag_Dirichlet = mesh->markerName( "Dirichlet" );
412 
413  }//end if(first_time)
414 
415 
420  // [toto1]
421  P0h = p0_space_type::New( mesh );
422  P1h = p1_space_type::New( mesh );
423  space_ptrtype Xh = space_type::New( mesh );
424  element_type u( Xh, "u" );
425  element_type v( Xh, "v" );
426  // [toto1]
427 
428 
429 
435  // [toto2]
436  //# marker1 #
437  value_type pi = M_PI;
439  auto fn1 = ( 1-Px()*Px() )*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() );
440  auto fn2 = sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() );
441  auto g = chi( fn==1 )*fn1 + chi( fn==2 )*fn2;
442  auto grad_g =
443  trans( chi( fn==1 )*( ( -2*Px()*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() )+beta*fn1 )*unitX()+
444  -2*Py()*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*unitY()+
445  -2*Pz()*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*unitZ() )+
446  chi( fn==2 )*( +( alpha*pi*cos( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )+beta*fn2 )*unitX()+
447  -alpha*pi*sin( alpha*pi*Px() )*sin( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )*unitY()+
448  -alpha*pi*sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*sin( alpha*pi*Pz() )*exp( beta*Px() )*unitZ() ) );
450  auto minus_laplacian_g =
451  ( 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 +
452  2*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*chi( Dim >= 2 ) +
453  2*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*chi( Dim == 3 ) ) +
454  chi( fn == 2 )* (
455  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() ) )*
456  ( cos( alpha*pi*Py() )*chi( Dim>=2 ) + chi( Dim==1 ) ) * ( cos( alpha*pi*Pz() )*chi( Dim==3 ) + chi( Dim<=2 ) )
457  )
458 
459  );
460 
461  //# endmarker1 #
462  // [toto2]
463 
464 
465  using namespace Feel::vf;
466 
473  // [toto3]
474  //# marker2 #
475  vector_ptrtype F( M_backend->newVector( Xh ) );
476  form1( _test=Xh, _vector=F, _init=true ) =
477  integrate( elements( mesh ), minus_laplacian_g*id( v ) )+
478  integrate( markedfaces( mesh, tag_Neumann ),
479  grad_g*vf::N()*id( v ) );
480 
481  //# endmarker2 #
482  if ( this->comm().size() != 1 || weakdir )
483  {
484  //# marker41 #
485  form1( _test=Xh, _vector=F ) +=
486  integrate( markedfaces( mesh,tag_Dirichlet ), g*( -grad( v )*vf::N()+penaldir*id( v )/hFace() ) );
487  //# endmarker41 #
488  }
489 
490  F->close();
491 
492  // [toto3]
493 
499  //# marker3 #
500  // [toto4]
501  sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );
502  // [toto4]
503 
508  // [toto5]
509  form2( Xh, Xh, D, _init=true ) =
510  integrate( elements( mesh ), gradt( u )*trans( grad( v ) ) );
511  // [toto5]
512  //# endmarker3 #
513 
514  if ( this->comm().size() != 1 || weakdir )
515  {
522  // [toto6]
523  //# marker10 #
524  form2( Xh, Xh, D ) +=
525  integrate( markedfaces( mesh,tag_Dirichlet ),
526  -( gradt( u )*vf::N() )*id( v )
527  -( grad( v )*vf::N() )*idt( u )
528  +penaldir*id( v )*idt( u )/hFace() );
529  D->close();
530  //# endmarker10 #
531  // [toto6]
532  }
533 
534  else
535  {
541  // [toto7]
542  //# marker5 #
543  D->close();
544  form2( Xh, Xh, D ) +=
545  on( markedfaces( mesh, tag_Dirichlet ), u, F, g );
546  //# endmarker5 #
547  // [toto7]
548 
549  }
550 
551 
556  // [toto8]
558  backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );
559  // [toto8]
560 
565  // [toto9]
566  //# marker7 #
567  double L2exact = math::sqrt( integrate( elements( mesh ),g*g ).evaluate()( 0,0 ) );
568  double L2error2 =integrate( elements( mesh ),( idv( u )-g )*( idv( u )-g ) ).evaluate()( 0,0 );
569  double L2error = math::sqrt( L2error2 );
570  // [toto9]
571  double semiH1error2 = integrate( elements( mesh ),( gradv( u )-grad_g )*( gradv( u )-grad_g ) ).evaluate()( 0,0 );
572  double H1error = math::sqrt( L2error2+semiH1error2 );
573  double H1exact = math::sqrt( integrate( elements( mesh ),g*g+grad_g*trans( grad_g ) ).evaluate()( 0,0 ) );
574 
575 
576  auto RealL2ErrorP0 = integrate( elements( mesh ), ( idv( u )-g )*( idv( u )-g ), _Q<10>() ).broken( P0h );
577  auto RealSemiH1ErrorP0 = integrate( elements( mesh ), ( gradv( u )-grad_g )* trans( gradv( u )-grad_g ), _Q<10>() ).broken( P0h );
578  auto H1RealErrorP0 = RealL2ErrorP0;
579  H1RealErrorP0 += RealSemiH1ErrorP0;
580  H1RealErrorP0 = H1RealErrorP0.sqrt();
581  /*******************residual estimator**********************/
582  //the source terme is given by : minus_laplacian_g
583 
584 
585  auto term1 = vf::h()*( minus_laplacian_g+trace( hessv( u ) ) );
586  auto term2 = jumpv( gradv( u ) );
587  auto term3 = gradv( u )*vf::N()-grad_g*N();
588 
589 
590  auto rho = integrate( elements( mesh ), term1*term1, _Q<10>() ).broken( P0h ).sqrt();
591 
592 
593  rho += integrate( internalfaces( mesh ),0.25*vf::h()*term2*term2 ).broken( P0h ).sqrt();
594 
595  rho += integrate( markedfaces( mesh,tag_Neumann ),
596  vf::h()*term3*term3 ).broken( P0h ).sqrt();
597 
598 
599  auto h=vf::project( P0h, elements( mesh ), vf::h() );
600  auto npen=vf::project( P0h, elements( mesh ), vf::nPEN() );
601  auto H1estimator = rho;
602 
603  //if we use real error for adaptation then H1errorP1 will be the projection of H1RealErrorP0 on P1 space
604  //else H1errorP1 will be the projection of H1estimator on P1 space
605  p1_element_type H1errorP1;
606 
607  if ( error_type==2 )
608  {
609  H1errorP1 = element_div( vf::sum( P1h, idv( H1RealErrorP0 )*meas() ), vf::sum( P1h, meas() ) );
610  }
611 
612  else if ( error_type==1 )
613  {
614  H1errorP1 = element_div( vf::sum( P1h, idv( H1estimator )*meas() ), vf::sum( P1h, meas() ) );
615  }
616 
617  else
618  {
619  std::cout<<"Problem with parameter adapt-error-type, please choice between 1 and 2"<<std::endl;
620  return;
621  }
622 
623  //auto new_hsize=vf::project( P1h, elements(mesh), vf::pow(vf::pow(vf::h(),Order)*(1e-4)/idv(H1estimatorP1),1./Order));
624  //auto new_hsize=vf::project( P1h, elements(mesh), vf::h()*(1e-2)/idv(H1estimatorP1) );
625 
626  estimatorH1=math::sqrt( H1estimator.pow( 2 ).sum() );
627  estimatorL2=math::sqrt( element_product( H1estimator,h ).pow( 2 ).sum() );
628 
629 
630  Y[0] = L2error/L2exact;
631  Y[1] = H1error/H1exact;
632  Y[2] = estimatorL2/L2exact;
633  Y[3] = estimatorH1/H1exact;
634 
635  h_new = P1h->element();
636 
637  h_new = vf::project( P1h, elements( mesh ),
638  vf::max( vf::pow(
639  vf::pow( vf::h(),Order )*( tol )/idv( H1errorP1 ),
640  1./Order ),
641  this->vm()["adapt-hmin"].template as<double>() ) );
642  /**********************end of residual estimaor*************/
643 
644 
647  element_type exact_solution( Xh, "exact_solution" );
648  exact_solution = vf::project( Xh, elements( mesh ), g );
649  element_type u_minus_exact( Xh, "u-g" );
650  u_minus_exact = vf::project( Xh, elements( mesh ), idv( u )-g );
651 
652 
653  if ( exporter->doExport() )
654  {
655  LOG(INFO) << "exportResults starts\n";
656 
657  exporter->step( 0 )->setMesh( mesh );
658  exporter->step( 0 )->add( "unknown", u );
659  exporter->step( 0 )->add( "exact solution", exact_solution );
660  exporter->step( 0 )->add( "u - exact solution", u_minus_exact ) ;
661  exporter->step( 0 )->add( "nPEN" , npen );
662  exporter->step( 0 )->add( "H1 error estimator P0" , H1estimator );
663  exporter->step( 0 )->add( "H1 Real error P0" , H1RealErrorP0 );
664  exporter->step( 0 )->add( "H1 error P1 (used for determination of new hsize)" , H1errorP1 );
665  exporter->step( 0 )->add( "new hsize" , h_new );
666 
667  exporter->save();
668  LOG(INFO) << "exportResults done\n";
669  }
670 
671 
672  LOG(INFO)<< " real L2 error : "<<Y[0]<<"\n";
673  LOG(INFO)<< " estimated L2 error "<<Y[2]<<"\n";
674  LOG(INFO)<< " real H1 error : "<<Y[1]<<"\n";
675  LOG(INFO)<< " estimated H1 error "<<Y[3]<<"\n";
676 
677 
678  std::ostringstream geostr;
679 
680  if ( this->vm()["gmshmodel"].template as<bool>() )
681  {
682  if ( this->vm()["gmshgeo"].template as<bool>() )
683  geostr << shape << "-" << Dim << ".geo";
684 
685  else
686  geostr << shape << "-" << Dim << ".msh";
687  }
688 
689  mesh = adapt( _h=h_new,
690  _model=geostr.str(),
691  _hmin=this->vm()["adapt-hmin"].template as<double>(),
692  _hmax=this->vm()["adapt-hmax"].template as<double>() );
693 } // ResidualEstimator::run
694 
695 
696 
697 
698 
699 
700 

Generated on Fri Oct 25 2013 14:24:24 for Feel++ by doxygen 1.8.4