30 #include <feel/feel.hpp>
32 #if defined (FEELPP_HAS_MADLIB_H)
34 #endif // FEELPP_HAS_MADLIB_H
39 using namespace boost::numeric::ublas;
56 "nD(n=1,2,3) Residual Estimator on Laplacian equation",
57 Feel::AboutData::License_GPL,
58 "Copyright (c) 2010 Universite Joseph Fourier" );
60 about.
addAuthor(
"Christophe Prud'homme",
"developer",
"christophe.prudhomme@feelpp.org",
"" );
61 about.
addAuthor(
"StÃÂÃÂÃÂÃÂÃÂÃÂÃÂéphane Veys",
"developer",
"stephane.veys@gmail.com",
"" );
72 template<
int Dim,
int Order = 1>
106 typedef boost::shared_ptr<p0_space_type> p0_space_ptrtype;
114 typedef boost::shared_ptr<p1_space_type> p1_space_ptrtype;
115 typedef typename p1_space_type::element_type p1_element_type;
142 exporter(
export_type::New(
"gmsh", this->about().appName() ) ),
145 shape(
"hypercube" ),
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>() )
179 void run(
const double* X,
unsigned long P,
double* Y,
unsigned long N );
182 BOOST_PARAMETER_CONST_MEMBER_FUNCTION(
192 ( maxit, *( boost::is_integral<mpl::_> ), 10 )
193 ( hmin, *( boost::is_arithmetic<mpl::_> ), 1e-2 )
194 ( hmax, *( boost::is_arithmetic<mpl::_> ), 2 )
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 ) )
202 #if defined (FEELPP_HAS_MADLIB_H)
203 saveGMSHMesh( _filename=
"inputmesh.msh",
204 _parametricnodes=!model.empty(),
206 MAd::pGModel themodel = 0;
207 GM_create( &themodel,
"theModel" );
209 if ( !model.empty() )
210 GM_read( themodel, model );
213 GM_read( themodel,
"inputmesh.msh" );
215 MAd::pMesh amesh = MAd::M_new( themodel );
216 MAd::M_load( amesh,
"inputmesh.msh" );
218 MAd::PWLSField * sizeField =
new MAd::PWLSField( amesh );
219 sizeField->setCurrentSize();
220 auto _elit = h.mesh()->beginElement();
221 auto _elen = h.mesh()->endElement();
223 for ( ; _elit != _elen; ++_elit )
225 for (
int l = 0; l < _elit->numPoints; ++l )
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 ) );
233 MAd::MeshAdapter* ma =
new MAd::MeshAdapter( amesh,sizeField );
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 );
244 ma->setGeoTracking( !model.empty() );
256 std::cout <<
"Statistics before optimization: \n";
257 ma->printStatistics( std::cout );
258 ma->writePos(
"meanRatioBefore.pos",MAd::OD_MEANRATIO );
269 std::cout <<
"Statistics after optimization: \n";
270 ma->printStatistics( std::cout );
271 ma->writePos(
"meanRatioAfter.pos",MAd::OD_MEANRATIO );
273 MAd::M_writeMsh ( amesh,
"result.msh", 2, NULL );
275 return loadGMSHMesh( _mesh=
new mesh_type,
276 _filename=
"result.msh",
279 #endif // FEELPP_HAS_MADLIB_H
286 backend_ptrtype M_backend;
287 backend_ptrtype M_backendP1;
316 export_ptrtype exporter;
319 p0_space_ptrtype P0h;
320 p1_space_ptrtype P1h;
323 p1_element_type h_new;
324 std::string msh_name;
331 template<
int Dim,
int Order>
335 if ( dim && dim != Dim ) return ;
337 if ( order && order != Order ) return ;
339 std::cout <<
"------------------------------------------------------------\n";
340 std::cout <<
"Execute ResidualEstimator<" << Dim <<
">\n";
341 std::vector<double> X;
342 X.push_back( meshSize );
344 if ( shape ==
"hypercube" )
347 else if ( shape ==
"ellipsoid" )
354 X.push_back( alpha );
356 X.push_back( weakdir );
357 X.push_back( penaldir );
359 X.push_back( first_time );
360 X.push_back( error_type );
362 std::vector<double> Y( 4 );
363 run( X.data(), X.size(), Y.data(), Y.size() );
366 template<
int Dim,
int Order>
375 if ( X[1] == 0 ) shape =
"simplex";
377 if ( X[1] == 1 ) shape =
"hypercube";
379 if ( X[1] == 2 ) shape =
"ellipsoid";
390 double estimatorH1, estimatorL2, estimator;
394 if ( !this->vm().count(
"nochdir" ) )
395 Environment::changeRepository( boost::format(
"doc/tutorial/%1%/%2%-%3%/P%4%/h_%5%/" )
396 % this->about().appName()
402 mesh = createGMSHMesh( _mesh=
new mesh_type,
404 _desc=domain( _name=( boost::format(
"%1%-%2%" ) % shape % Dim ).str() ,
409 _update=MESH_CHECK|MESH_UPDATE_FACES|MESH_UPDATE_EDGES );
410 tag_Neumann = mesh->markerName(
"Neumann" );
411 tag_Dirichlet = mesh->markerName(
"Dirichlet" );
421 P0h = p0_space_type::New( mesh );
422 P1h = p1_space_type::New( mesh );
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;
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 ) ) +
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 ) )
465 using namespace Feel::vf;
476 form1( _test=Xh, _vector=F, _init=
true ) =
479 grad_g*vf::N()*
id( v ) );
482 if ( this->comm().size() != 1 || weakdir )
485 form1( _test=Xh, _vector=F ) +=
486 integrate(
markedfaces( mesh,tag_Dirichlet ), g*( -grad( v )*vf::N()+penaldir*
id( v )/hFace() ) );
509 form2( Xh, Xh, D, _init=
true ) =
514 if ( this->comm().size() != 1 || weakdir )
524 form2( Xh, Xh, D ) +=
526 -( gradt( u )*vf::N() )*
id( v )
527 -( grad( v )*vf::N() )*idt( u )
528 +penaldir*
id( v )*idt( u )/hFace() );
544 form2( Xh, Xh, D ) +=
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 );
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 ) );
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();
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();
590 auto rho =
integrate(
elements( mesh ), term1*term1, _Q<10>() ).broken( P0h ).sqrt();
596 vf::h()*term3*term3 ).broken( P0h ).sqrt();
601 auto H1estimator = rho;
605 p1_element_type H1errorP1;
609 H1errorP1 =
element_div( vf::sum( P1h, idv( H1RealErrorP0 )*meas() ), vf::sum( P1h, meas() ) );
612 else if ( error_type==1 )
614 H1errorP1 =
element_div( vf::sum( P1h, idv( H1estimator )*meas() ), vf::sum( P1h, meas() ) );
619 std::cout<<
"Problem with parameter adapt-error-type, please choice between 1 and 2"<<std::endl;
626 estimatorH1=math::sqrt( H1estimator.pow( 2 ).sum() );
627 estimatorL2=math::sqrt(
element_product( H1estimator,h ).pow( 2 ).sum() );
630 Y[0] = L2error/L2exact;
631 Y[1] = H1error/H1exact;
632 Y[2] = estimatorL2/L2exact;
633 Y[3] = estimatorH1/H1exact;
635 h_new = P1h->element();
639 vf::pow( vf::h(),Order )*( tol )/idv( H1errorP1 ),
641 this->vm()[
"adapt-hmin"].
template as<double>() ) );
653 if ( exporter->doExport() )
655 LOG(INFO) <<
"exportResults starts\n";
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 );
668 LOG(INFO) <<
"exportResults done\n";
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";
678 std::ostringstream geostr;
680 if ( this->vm()[
"gmshmodel"].
template as<bool>() )
682 if ( this->vm()[
"gmshgeo"].
template as<bool>() )
683 geostr << shape <<
"-" << Dim <<
".geo";
686 geostr << shape <<
"-" << Dim <<
".msh";
689 mesh = adapt( _h=h_new,
691 _hmin=this->vm()[
"adapt-hmin"].template as<double>(),
692 _hmax=this->vm()[
"adapt-hmax"].template as<double>() );
Simulation Object.
Definition: simget.hpp:60
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::marker_face_const_iterator, typename MeshTraits< MeshType >::marker_face_const_iterator > markedfaces(MeshType const &mesh)
Definition: filters.hpp:969
Poly::value_type integrate(Polynomial< Poly, PsetType > const &p)
integrate a polynomial over a convex
Definition: operations.hpp:51
elements_type const & elements() const
Definition: elements.hpp:355
Definition: solverlinear.hpp:33
backend_type::vector_type vector_type
vector type associated with backend
Definition: residualestimator.hpp:93
simplex of dimension Dim
Definition: simplex.hpp:73
Polynomial< Pset, Scalar > project(Pset const &pset, Func const &f, IM const &im)
Definition: operations.hpp:436
VectorEigen< T > element_product(VectorEigen< T > const &v1, VectorEigen< T > const &v2)
Definition: vectoreigen.hpp:782
export Feel generated data to some file formatsUse the visitor and factory pattern.
Definition: exporter.hpp:82
Backend< value_type > backend_type
linear algebra backend factory
Definition: residualestimator.hpp:83
bases< Lagrange< Order, Scalar > > basis_type
the basis type of our approximation space
Definition: residualestimator.hpp:119
Exporter< mesh_type, 1 > export_type
the exporter factory type
Definition: residualestimator.hpp:129
backend_type::sparse_matrix_type sparse_matrix_type
sparse matrix type associated with backend
Definition: residualestimator.hpp:89
Definition: residualestimator.hpp:73
boost::shared_ptr< export_type > export_ptrtype
the exporter factory (shared_ptr<> type)
Definition: residualestimator.hpp:131
AboutData makeAbout()
Create the About data of the OpusApp.
Definition: opusappeigs.cpp:49
Holds information needed by the "About" box and other classes.
Definition: about.hpp:173
boost::shared_ptr< space_type > space_ptrtype
the approximation function space type (shared_ptr<> type)
Definition: residualestimator.hpp:124
p0_space_type::element_type p0_element_type
an element type of the discontinuous function space
Definition: residualestimator.hpp:108
base class for all linear algebra backends
Definition: backend.hpp:133
double value_type
numerical type is double
Definition: residualestimator.hpp:80
space_type::element_type element_type
an element type of the approximation function space
Definition: residualestimator.hpp:126
unifying mesh class
Definition: createsubmesh.hpp:41
static backend_ptrtype build(BackendType=BACKEND_GMM, WorldComm const &worldComm=Environment::worldComm())
Definition: backend.cpp:167
ResidualEstimator(AboutData const &about)
Definition: residualestimator.hpp:136
FunctionSpace< mesh_type, basis_type > space_type
the approximation function space type
Definition: residualestimator.hpp:122
Definition: functionspace.hpp:1374
Simplex< Dim > convex_type
geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1 ...
Definition: residualestimator.hpp:98
backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype
sparse matrix type associated with backend (shared_ptr<> type)
Definition: residualestimator.hpp:91
Definition: matrixsparse.hpp:76
Mesh< convex_type > mesh_type
mesh type
Definition: residualestimator.hpp:100
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:105
boost::tuple< mpl::size_t< MESH_FACES >, typename MeshTraits< MeshType >::location_face_const_iterator, typename MeshTraits< MeshType >::location_face_const_iterator > internalfaces(MeshType const &mesh)
Definition: filters.hpp:1175
bases< Lagrange< 1, Scalar > > p1_basis_type
the basis type of our approximation space
Definition: residualestimator.hpp:112
backend_type::vector_ptrtype vector_ptrtype
vector type associated with backend (shared_ptr<> type)
Definition: residualestimator.hpp:95
ElementType element_div(ElementType const &v1, ElementType const &v2)
Definition: functionspace.hpp:6268
boost::shared_ptr< backend_type > backend_ptrtype
linear algebra backend factory shared_ptr<> type
Definition: residualestimator.hpp:85
void addAuthor(std::string const &name, std::string const &task=0, std::string const &emailAddress=0, std::string const &webAddress=0)
Definition: about.cpp:143
void run()
Definition: residualestimator.hpp:333
boost::shared_ptr< mesh_type > mesh_ptrtype
mesh shared_ptr<> type
Definition: residualestimator.hpp:102