30 #include <feel/feel.hpp>
31 #include <feel/feeldiscr/elementdiv.hpp>
32 #include <feel/feelfilters/creategmshmesh.hpp>
49 AboutData about(
"residualestimator" ,
52 "nD(n=1,2,3) Residual Estimator on Laplacian equation",
53 Feel::AboutData::License_GPL,
54 "Copyright (c) 2010 Universite Joseph Fourier" );
56 about.addAuthor(
"Christophe Prud'homme",
"developer",
"christophe.prudhomme@feelpp.org",
"" );
57 about.addAuthor(
"StÃÂÃÂÃÂÃÂÃÂÃÂÃÂéphane Veys",
"developer",
"stephane.veys@gmail.com",
"" );
68 template<
int Dim,
int Order = 1>
101 typedef FunctionSpace<mesh_type, bases<Lagrange<0,Scalar,Discontinuous> > >
p0_space_type;
102 typedef boost::shared_ptr<p0_space_type> p0_space_ptrtype;
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;
135 M_backend( backend_type::build() ),
136 M_backendP1( backend_type::build() ),
138 exporter( export_type::New(
"gmsh", this->about().appName() ) ),
141 shape(
"hypercube" ),
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>() )
175 void run(
const double* X,
unsigned long P,
double* Y,
unsigned long N );
178 BOOST_PARAMETER_CONST_MEMBER_FUNCTION(
188 ( maxit, *( boost::is_integral<mpl::_> ), 10 )
189 ( hmin, *( boost::is_arithmetic<mpl::_> ), 1e-2 )
190 ( hmax, *( boost::is_arithmetic<mpl::_> ), 2 )
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 ) )
204 backend_ptrtype M_backend;
205 backend_ptrtype M_backendP1;
234 export_ptrtype exporter;
237 p0_space_ptrtype P0h;
238 p1_space_ptrtype P1h;
241 p1_element_type h_new;
242 std::string msh_name;
249 template<
int Dim,
int Order>
253 if ( dim && dim != Dim ) return ;
255 if ( order && order != Order ) return ;
257 std::cout <<
"------------------------------------------------------------\n";
258 std::cout <<
"Execute ResidualEstimator<" << Dim <<
">\n";
259 std::vector<double> X;
260 X.push_back( meshSize );
262 if ( shape ==
"hypercube" )
265 else if ( shape ==
"ellipsoid" )
272 X.push_back( alpha );
274 X.push_back( weakdir );
275 X.push_back( penaldir );
277 X.push_back( first_time );
278 X.push_back( error_type );
280 std::vector<double> Y( 4 );
281 run( X.data(), X.size(), Y.data(), Y.size() );
284 template<
int Dim,
int Order>
293 if ( X[1] == 0 ) shape =
"simplex";
295 if ( X[1] == 1 ) shape =
"hypercube";
297 if ( X[1] == 2 ) shape =
"ellipsoid";
308 double estimatorH1, estimatorL2, estimator;
312 if ( !this->vm().count(
"nochdir" ) )
313 Environment::changeRepository( boost::format(
"doc/tutorial/%1%/%2%-%3%/P%4%/h_%5%/" )
314 % this->about().appName()
320 mesh = createGMSHMesh( _mesh=
new mesh_type,
322 _desc=domain( _name=( boost::format(
"%1%-%2%" ) % shape % Dim ).str() ,
327 _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES );
328 tag_Neumann = mesh->markerName(
"Neumann" );
329 tag_Dirichlet = mesh->markerName(
"Dirichlet" );
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" );
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;
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 ) ) +
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 ) )
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 ) );
400 if ( this->comm().size() != 1 || weakdir )
403 form1( _test=Xh, _vector=F ) +=
404 integrate( markedfaces( mesh,tag_Dirichlet ), g*( -grad( v )*vf::N()+penaldir*
id( v )/hFace() ) );
419 sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );
427 form2( Xh, Xh, D, _init=
true ) =
428 integrate( elements( mesh ), gradt( u )*trans( grad( v ) ) );
432 if ( this->comm().size() != 1 || weakdir )
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() );
462 form2( Xh, Xh, D ) +=
463 on( markedfaces( mesh, tag_Dirichlet ), u, F, g );
476 backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );
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 );
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 ) );
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();
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();
508 auto rho = integrate( elements( mesh ), term1*term1, _Q<10>() ).broken( P0h ).sqrt();
511 rho += integrate( internalfaces( mesh ),0.25*vf::h()*term2*term2 ).broken( P0h ).sqrt();
513 rho += integrate( markedfaces( mesh,tag_Neumann ),
514 vf::h()*term3*term3 ).broken( P0h ).sqrt();
517 auto h=vf::project( P0h, elements( mesh ), vf::h() );
518 auto npen=vf::project( P0h, elements( mesh ), vf::nPEN() );
519 auto H1estimator = rho;
523 p1_element_type H1errorP1;
527 H1errorP1 = div( vf::sum( P1h, idv( H1RealErrorP0 )*meas() ), vf::sum( P1h, meas() ) );
530 else if ( error_type==1 )
532 H1errorP1 = div( vf::sum( P1h, idv( H1estimator )*meas() ), vf::sum( P1h, meas() ) );
537 std::cout<<
"Problem with parameter adapt-error-type, please choice between 1 and 2"<<std::endl;
544 estimatorH1=math::sqrt( H1estimator.pow( 2 ).sum() );
545 estimatorL2=math::sqrt( element_product( H1estimator,h ).pow( 2 ).sum() );
548 Y[0] = L2error/L2exact;
549 Y[1] = H1error/H1exact;
550 Y[2] = estimatorL2/L2exact;
551 Y[3] = estimatorH1/H1exact;
553 h_new = P1h->element();
555 h_new = vf::project( P1h, elements( mesh ),
557 vf::pow( vf::h(),Order )*( tol )/idv( H1errorP1 ),
559 this->vm()[
"adapt-hmin"].
template as<double>() ) );
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 );
571 if ( exporter->doExport() )
573 LOG(INFO) <<
"exportResults starts\n";
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 );
586 LOG(INFO) <<
"exportResults done\n";
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";
596 std::ostringstream geostr;
598 if ( this->vm()[
"gmshmodel"].
template as<bool>() )
600 if ( this->vm()[
"gmshgeo"].
template as<bool>() )
601 geostr << shape <<
"-" << Dim <<
".geo";
604 geostr << shape <<
"-" << Dim <<
".msh";
607 mesh = adapt( _h=h_new,
609 _hmin=this->vm()[
"adapt-hmin"].template as<double>(),
610 _hmax=this->vm()[
"adapt-hmax"].template as<double>() );
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
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