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

#include <residualestimator.hpp>

Detailed Description

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

Laplacian Solver using continuous approximation spaces solve \( -\Delta u = f \) on \(\Omega\) and \( u= g \) on \(\Gamma \)

Template Parameters
Dimthe geometric dimension of the problem (e.g. Dim=1, 2 or 3)
Orderthe approximation order

Inherits Simget.

Public Types

typedef boost::shared_ptr
< backend_type
backend_ptrtype
 linear algebra backend factory shared_ptr<> type
 
typedef Backend< value_typebackend_type
 linear algebra backend factory
 
typedef bases< Lagrange< Order,
Scalar > > 
basis_type
 the basis type of our approximation space
 
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
< export_type
export_ptrtype
 the exporter factory (shared_ptr<> type)
 
typedef Exporter< mesh_type, 1 > export_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 p0_space_type::element_type p0_element_type
 an element type of the \(P_0\) discontinuous function space
 
typedef boost::shared_ptr
< p0_space_type
p0_space_ptrtype
 
typedef FunctionSpace
< mesh_type, bases< Lagrange
< 0, Scalar, Discontinuous > > > 
p0_space_type
 function space that holds piecewise constant ( \(P_0\)) functions (e.g. to store material properties or partitioning
 
typedef bases< Lagrange
< 1, Scalar > > 
p1_basis_type
 the basis type of our approximation space
 
typedef p1_space_type::element_type p1_element_type
 
typedef boost::shared_ptr
< p1_space_type > 
p1_space_ptrtype
 
typedef FunctionSpace
< mesh_type, p1_basis_type
p1_space_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
backend_type::sparse_matrix_ptrtype 
sparse_matrix_ptrtype
 sparse matrix type associated with backend (shared_ptr<> type)
 
typedef
backend_type::sparse_matrix_type 
sparse_matrix_type
 sparse matrix type associated with backend
 
typedef double value_type
 numerical type is double
 
typedef
backend_type::vector_ptrtype 
vector_ptrtype
 vector type associated with backend (shared_ptr<> type)
 
typedef backend_type::vector_type vector_type
 vector type associated with backend
 

Public Member Functions

 ResidualEstimator (AboutData const &about)
 
void run ()
 
void run (const double *X, unsigned long P, double *Y, unsigned long N)
 

Constructor & Destructor Documentation

template<int Dim, int Order = 1>
ResidualEstimator< Dim, Order >::ResidualEstimator ( AboutData const &  about)
inline

Constructor

Member Function Documentation

template<int Dim, int Order>
void ResidualEstimator< Dim, Order >::run ( const double *  X,
unsigned long  P,
double *  Y,
unsigned long  N 
)

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

P0h = p0_space_type::New( mesh );
P1h = p1_space_type::New( mesh );
space_ptrtype Xh = space_type::New( mesh );
element_type u( Xh, "u" );
element_type v( 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 #
value_type pi = M_PI;
auto fn1 = ( 1-Px()*Px() )*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() );
auto fn2 = sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() );
auto g = chi( fn==1 )*fn1 + chi( fn==2 )*fn2;
auto grad_g =
trans( chi( fn==1 )*( ( -2*Px()*( 1-Py()*Py() )*( 1-Pz()*Pz() )*exp( beta*Px() )+beta*fn1 )*unitX()+
-2*Py()*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*unitY()+
-2*Pz()*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*unitZ() )+
chi( fn==2 )*( +( alpha*pi*cos( alpha*pi*Px() )*cos( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )+beta*fn2 )*unitX()+
-alpha*pi*sin( alpha*pi*Px() )*sin( alpha*pi*Py() )*cos( alpha*pi*Pz() )*exp( beta*Px() )*unitY()+
-alpha*pi*sin( alpha*pi*Px() )*cos( alpha*pi*Py() )*sin( alpha*pi*Pz() )*exp( beta*Px() )*unitZ() ) );
auto minus_laplacian_g =
( 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 +
2*( 1-Px()*Px() )*( 1-Pz()*Pz() )*exp( beta*Px() )*chi( Dim >= 2 ) +
2*( 1-Px()*Px() )*( 1-Py()*Py() )*exp( beta*Px() )*chi( Dim == 3 ) ) +
chi( fn == 2 )* (
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() ) )*
( cos( alpha*pi*Py() )*chi( Dim>=2 ) + chi( Dim==1 ) ) * ( cos( alpha*pi*Pz() )*chi( Dim==3 ) + chi( Dim<=2 ) )
)
);
//# endmarker1 #

deduce from expression the type of g (thanks to keyword 'auto')

deduce from expression the type of laplacian (thanks to keyword 'auto')

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

//# marker2 #
vector_ptrtype F( M_backend->newVector( Xh ) );
form1( _test=Xh, _vector=F, _init=true ) =
integrate( elements( mesh ), minus_laplacian_g*id( v ) )+
integrate( markedfaces( mesh, tag_Neumann ),
grad_g*vf::N()*id( v ) );
//# endmarker2 #
if ( this->comm().size() != 1 || weakdir )
{
//# marker41 #
form1( _test=Xh, _vector=F ) +=
integrate( markedfaces( mesh,tag_Dirichlet ), g*( -grad( v )*vf::N()+penaldir*id( v )/hFace() ) );
//# endmarker41 #
}
F->close();

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

sparse_matrix_ptrtype D( M_backend->newMatrix( Xh, Xh ) );

assemble \(\int_\Omega \nu \nabla u \cdot \nabla v\)

form2( Xh, Xh, D, _init=true ) =
integrate( elements( mesh ), 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 #
    form2( Xh, Xh, D ) +=
    integrate( markedfaces( mesh,tag_Dirichlet ),
    -( gradt( u )*vf::N() )*id( v )
    -( grad( v )*vf::N() )*idt( u )
    +penaldir*id( v )*idt( u )/hFace() );
    D->close();
    //# endmarker10 #

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

  4. first close the matrix (the matrix must be closed first before any manipulation )
  5. modify the matrix by cancelling out the rows and columns of D that are associated with the Dirichlet dof
    //# marker5 #
    D->close();
    form2( Xh, Xh, D ) +=
    on( markedfaces( mesh, tag_Dirichlet ), u, F, g );
    //# endmarker5 #

solve the system

backend_type::build()->solve( _matrix=D, _solution=u, _rhs=F );

solve \( D u = F \)

compute the \(L_2\) norm of the error

//# marker7 #
double L2exact = math::sqrt( integrate( elements( mesh ),g*g ).evaluate()( 0,0 ) );
double L2error2 =integrate( elements( mesh ),( idv( u )-g )*( idv( u )-g ) ).evaluate()( 0,0 );
double L2error = math::sqrt( L2error2 );

save the results project the exact solution


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