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

Dim	the 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_type >	bdf_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_type >	export_type
	the exporter factory type
typedef boost::shared_ptr < mesh_type >	mesh_ptrtype
	mesh shared_ptr<> type
typedef Mesh< convex_type >	mesh_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:

• laplacian_parabolic.cpp