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 Namespaces Files Functions Variables Typedefs Pages
Feel::Laplacian< Dim, Order, Cont, Entity, FType > Class Template Reference

Detailed Description

template<int Dim, int Order, typename Cont, template< uint16_type, uint16_type, uint16_type > class Entity, template< uint16_type > class FType>
class Feel::Laplacian< Dim, Order, Cont, Entity, FType >

Laplacian Solver using discontinous approximation spaces

solve \( -\Delta u = f\) on \(\Omega\) and \(u= g\) on \(\Gamma\)

Inherits Application.

Public Types

typedef boost::shared_ptr
< backend_type > 
backend_ptrtype
 
typedef Backend< value_type > backend_type
 
typedef Entity< Dim, 1, Dim > entity_type
 
typedef boost::shared_ptr
< export_type > 
export_ptrtype
 
typedef Exporter< mesh_type > export_type
 
typedef boost::shared_ptr
< mesh_type > 
mesh_ptrtype
 
typedef Mesh< entity_type > mesh_type
 
typedef p0_space_type::element_type p0_element_type
 
typedef FunctionSpace
< mesh_type, bases< Lagrange
< 0, Scalar > >, Discontinuous > 
p0_space_type
 
typedef
backend_type::sparse_matrix_ptrtype 
sparse_matrix_ptrtype
 
typedef
backend_type::sparse_matrix_type 
sparse_matrix_type
 
typedef double value_type
 
typedef
backend_type::vector_ptrtype 
vector_ptrtype
 
typedef backend_type::vector_type vector_type
 

Public Member Functions

void operator() ()
 
void run ()
 

Static Public Attributes

static const uint16_type imOrder = 2*Order
 

Member Function Documentation

template<int Dim, int Order, typename Cont , template< uint16_type, uint16_type, uint16_type > class Entity, template< uint16_type > class FType>
void Feel::Laplacian< Dim, Order, Cont, Entity, FType >::operator() ( )
inline

alias for run()

template<int Dim, int Order, typename Cont , template< uint16_type, uint16_type, uint16_type > class Entity, template< uint16_type > class FType>
void Laplacian< Dim, Order, Cont, Entity, FType >::run ( )

run the convergence test

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 gproj( Xh, "v" );

define \(g\) the expression of the exact solution and \(f\) the expression of the right hand side such that \(g\) is the exact solution

*/
//# marker1 #
bool weak_dirichlet = this->vm()["weakdir"].template as<int>();
value_type penaldir = this->vm()["penaldir"].template as<double>();
value_type nu = this->vm()["nu"].template as<double>();
std::string exact =
this->vm()[(boost::format("exact%1%D")%Dim).str()].template
as<std::string>();
std::string rhs =
this->vm()[(boost::format("rhs%1%D")%Dim).str()].template
as<std::string>();
LOG(INFO) << "exact = " << exact << "\n";
auto vars=symbols<Dim>();
ex ff;
/*
* if the right hand side is an empty string then use the string exact and
* compute its laplacian to set the right hand side
*/
if ( rhs.empty() )
ff=-nu*laplacian(exact,vars);
else
ff = parse(rhs,vars);
LOG(INFO) << "rhs="<< ff << "\n";
auto g = expr(exact,vars);
auto f = expr(ff,vars);
auto gradg = expr<1,Dim,2>(grad(exact,vars), vars );
LOG(INFO) << "grad(g)="<< grad(exact,vars) << "\n";
// build Xh-interpolant of g
gproj = project( _space=Xh, _range=elements( mesh ), _expr=g );
//# endmarker1 #

Construction of the right hand side. F is the vector that holds the algebraic representation of the right habd side of the problem

*/
//# marker2 #
auto F = backend()->newVector( Xh );
auto l = form1( _test=Xh, _vector=F );
l = integrate( _range=elements( mesh ), _expr=f*id( v ) )+
integrate( _range=markedfaces( mesh, "Neumann" ),
_expr=nu*gradg*vf::N()*id( v )
);
//# endmarker2 #
if ( weak_dirichlet )
{
//# marker41 #
l += integrate( _range=markedfaces( mesh,"Dirichlet" ),
_expr=nu*g*( -grad( v )*vf::N()
+ penaldir*id( v )/hFace() ) );
//# endmarker41 #
}

create the matrix that will hold the algebraic representation of the left hand side

*/
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 ) ) );

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\)
*/
//# marker10 #
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() ) );
//# endmarker10 #

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 #
a += on( _range=markedfaces( mesh, "Dirichlet" ),
_element=u, _rhs=F, _expr=g );
//# endmarker5 #

solve the system

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

compute the

*/
//# marker7 #
double L2error =normL2( _range=elements( mesh ),_expr=( idv( u )-g ) );
LOG(INFO) << "||error||_L2=" << L2error << "\n";
//# endmarker7 #

save the results

*/
element_type e( Xh, "e" );
e = project( _space=Xh, _range=elements( mesh ), _expr=g );
export_ptrtype exporter( export_type::New() );
if ( exporter->doExport() )
{
LOG(INFO) << "exportResults starts\n";
exporter->step( 0 )->setMesh( mesh );
exporter->step( 0 )->add( "solution", u );
exporter->step( 0 )->add( "exact", e );
exporter->save();
LOG(INFO) << "exportResults done\n";
}

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