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
Laplacian_parabolic< Dim, Order > Class Template Reference

Detailed Description

template<int Dim, int Order = 1>
class Laplacian_parabolic< Dim, Order >

Laplacian equation with instationnary term solver using continuous approximation spaces solve \( \dfrac{\partial u}{\partial t} - nu*\Delta u = f\) on \(\Omega\) and \(u= g\) on \(\Gamma\)

Template Parameters
Dimthe geometric dimension of the problem = 2
+ Inheritance diagram for Laplacian_parabolic< Dim, Order >:

Public Types

typedef bases< Lagrange< Order,
Scalar > > 
basis_type
 the basis type of our approximation space
 
typedef boost::shared_ptr
< bdf_type > 
bdf_ptrtype
 
typedef Bdf< space_typebdf_type
 
typedef Simplex< Dim > convex_type
 geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1
 
typedef space_type::element_type element_type
 an element type of the approximation function space
 
typedef boost::shared_ptr
< error_type
error_ptrtype
 
typedef ErrorBase< Dim, Order > error_type
 
typedef boost::shared_ptr
< export_type
export_ptrtype
 the exporter factory (shared_ptr<> type)
 
typedef Exporter< mesh_typeexport_type
 the exporter factory type
 
typedef boost::shared_ptr
< mesh_type
mesh_ptrtype
 mesh shared_ptr<> type
 
typedef Mesh< convex_typemesh_type
 mesh type
 
typedef boost::shared_ptr
< space_type
space_ptrtype
 the approximation function space type (shared_ptr<> type)
 
typedef FunctionSpace
< mesh_type, basis_type
space_type
 the approximation function space type
 
typedef double value_type
 numerical type is double
 

Public Member Functions

 Laplacian_parabolic ()
 
void run ()
 Function to compute the equation and find the unknown. More...
 

Constructor & Destructor Documentation

template<int Dim, int Order = 1>
Laplacian_parabolic< Dim, Order >::Laplacian_parabolic ( )
inline

Constructor

Member Function Documentation

template<int Dim, int Order>
void Laplacian_parabolic< Dim, Order >::run ( )

Function to compute the equation and find the unknown.

Loading variables from cfg file

Creation of a new mesh depending on the information of the geofile

*/
//std::string access_geofile = ( boost::format( "%1%.geo" ) % geofile ).str();
//gmsh_ptrtype desc_geo;
//desc_geo=geo(_filename = access_geofile, _dim = Dim, _order = 1, _h = meshSize);
//mesh_ptrtype mesh = createGMSHMesh( _mesh=new mesh_type,
// _desc=desc_geo,
// _update=MESH_UPDATE_EDGES|MESH_UPDATE_FACES );
mesh_ptrtype mesh = loadMesh(_mesh=new mesh_type,_filename=this->vm()["geofile"].template as<std::string>());

The function space and some associated elements (functions) are then defined

*/
space_ptrtype Xh = space_type::New( mesh );
// print some information (number of local/global dof in logfile)
Xh->printInfo();
element_type u( Xh, "u" );
element_type v( Xh, "v" );
element_type Rhs( Xh, "rhs" );
element_type gproj(Xh, "exact_proj");

[marker11]

[marker11]

[marker12]

[marker12]

Add extra parameters ( t for example )

*/
auto parameters = cvg->getParams(); //symbols<Dim>();

[marker13]

[marker13]

Initializing u, g and f from initial temperature expression

*/
initial_u_expr = parse(initial_u, vars, parameters);

BDF implementation

*/
// create the BDF structure
bdf_ptrtype M_bdf = bdf( _space=Xh, _name="u" );
// start the time at time-initial
M_bdf->start();
// create the finite difference polynome of the unknown
M_bdf->initialize(u);

create the matrix that will hold the algebraic representation of the left hand side (only stationnary terms)

[marker3]

*/
auto D = backend()->newMatrix( _test=Xh, _trial=Xh );

assemble $ u v$

*/
auto a = form2( _test=Xh, _trial=Xh, _matrix=D );
a = integrate( _range=elements( mesh ), _expr=nu*gradt( u )*trans( grad( v ) ) );

[marker3]

weak dirichlet conditions treatment for the boundaries marked 1 and 3

  1. assemble \(\int_{\partial \Omega} -\nabla u \cdot \mathbf{n} v\)
  2. assemble \(\int_{\partial \Omega} -\nabla v \cdot \mathbf{n} u\)
  3. assemble \(\int_{\partial \Omega} \frac{\gamma}{h} u v\)
*/
a += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr= nu * ( -( gradt( u )*vf::N() )*id( v )
-( grad( v )*vf::N() )*idt( u )
+penaldir*id( v )*idt( u )/hFace() ) );

assemble $ u^{n+1} v$

*/
if( !steady )
a += integrate( _range=elements( mesh ), _expr = cst(M_bdf->polyDerivCoefficient( 0 ))*idt(u)*id(v) );
// time depending term
do{ // temporal loop start
if(cvg->computedrhs())
if ( parameters.size() )
Rhs = cvg->rhs_project(Xh, t0);
else
Rhs = cvg->rhs_project(Xh);
else
Rhs = vf::project( Xh, elements(mesh), cst(0.0) );
if ( !steady && parameters.size() && !cvg->getExactSolution().empty() )
gproj = cvg->exact_project(Xh, M_bdf->time() );
else
gproj = cvg->exact_project(Xh);
//# marker2 #
auto F = backend()->newVector( Xh );
auto l = form1( _test=Xh, _vector=F );
l = integrate( _range=elements( mesh ), _expr=idv(Rhs)*id( v ) );
//# endmarker2 #
if ( weak_dirichlet )
{
//# marker41 #
l += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr=nu*idv(gproj)*( -grad( v )*vf::N()
+ penaldir*id( v )/hFace() ) );
//# endmarker41 #
}

add temporal term to the lhs and the rhs

*/
if( !steady )
{
auto bdf_poly = M_bdf->polyDeriv();
l += integrate( _range = elements(mesh), _expr = idv( bdf_poly )*id(v));
}

add time depending terms for the left hand side

*/
auto Dt = backend()->newMatrix( _test=Xh, _trial=Xh );
auto at = form2( _test=Xh, _trial=Xh, _matrix=Dt );

strong(algebraic) dirichlet conditions treatment for the boundaries marked 1 and 3

  1. first close the matrix (the matrix must be closed first before any manipulation )
  2. modify the matrix by cancelling out the rows and columns of D that are associated with the Dirichlet dof
*/
//# marker5 #
at = on( _range=markedfaces( mesh, "Dirichlet" ),
_element=u, _rhs=F, _expr=idv(gproj) );
//# endmarker5 #
*/
at += a;

solve the system

*/
//# marker6 #
backend( _rebuild=true )->solve( _matrix=Dt, _solution=u, _rhs=F );
//# endmarker6 #

Computing L2 and H1 error

*/
//# marker7 #
if ( !cvg->getExactSolution().empty() )
{
ex solution;
if( !steady && parameters.size() )
solution = cvg->getSolution(M_bdf->time() );
else
solution = cvg->getSolution();
auto g = expr(solution,vars);
auto gradg = expr<1,Dim,2>(grad(solution,vars), vars );
double L2error = normL2( _range=elements( mesh ),_expr=( idv( u )-idv(gproj) ) );
double H1seminorm = math::sqrt( integrate( elements(mesh), (gradv(u) - gradg)*trans(gradv(u) - gradg) ).evaluate()(0,0) );
double H1error = math::sqrt( L2error*L2error + H1seminorm*H1seminorm);
LOG(INFO) << "||error||_L2=" << L2error << "\n";
LOG(INFO) << "||error||_H1=" << H1error << "\n";
std::cout << M_bdf->time() <<"\t" << L2error << "\t" << H1error << "\n";
L2Time_error = L2Time_error + L2error*L2error;
}
//# endmarker7 #

save the results

*/
if ( exp->doExport() )
{
LOG(INFO) << "exportResults starts at " << M_bdf->time() << "\n";
exp->step( M_bdf->time() )->add( "solution", u );
exp->step( M_bdf->time() )->add( "exact", gproj );
exp->step( M_bdf->time() )->add( "rhs", Rhs );
exp->save();
LOG(INFO) << "exportResults done\n";
}

[marker15]

[marker15]

[marker16]

[marker16]


The documentation for this class was generated from the following file: