Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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
+ Inheritance diagram for ResidualEstimator< Dim, Order >:

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)
 
- Public Member Functions inherited from Feel::Simget
 Simget ()
 
virtual ~Simget ()
 destructor
 
Simgetoperator= (Simget const &o)
 copy operator
 
virtual std::string name () const
 return the name of the simget
 
mpi::communicator comm () const
 
po::variables_map const & vm () const
 
AboutData const & about () const
 
double meshSize () const
 return the mesh size
 
double meshSizeInit () const
 return the mesh size
 
int level () const
 return the refinement level
 
ptree::ptree const & stats () const
 return the statistics associated to the simget after calling run
 
ptree::ptree & stats ()
 return the statistics associated to the simget after calling run
 
void setMeshSize (double h)
 set the mesh size
 
void setMeshSizeInit (double h)
 set the initial mesh size
 
void setLevel (int level)
 set the refinment level if applicable
 
void print (std::ostream &out, std::vector< ptree::ptree > &stats)
 

Additional Inherited Members

- Protected Member Functions inherited from Feel::Simget
SimgetchangeRepository (boost::format fmt)
 
- Protected Attributes inherited from Feel::Simget
int M_level
 
double M_meshSize
 
double M_meshSizeInit
 
ptree::ptree M_stats
 

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 ( )
virtual

simply execute the simget

Implements Feel::Simget.

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

models the input/output relation $ Y=F(X) $

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

P0h = p0_space_type::New( mesh );
P1h = p1_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

Reimplemented from Feel::Simget.

References Feel::Backend< T >::build(), Feel::element_div(), Feel::element_product(), Feel::elements(), Feel::integrate(), Feel::internalfaces(), Feel::markedfaces(), and Feel::project().


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

Generated on Sun Dec 22 2013 13:11:16 for Feel++ by doxygen 1.8.5