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
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...
 
- 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
 
virtual void run (const double *X, unsigned long P, double *Y, unsigned long N)
 
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>
Laplacian_parabolic< Dim, Order >::Laplacian_parabolic ( )
inline

Constructor

Member Function Documentation

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

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

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

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]

Implements Feel::Simget.

References Feel::elements(), Feel::integrate(), Feel::markedfaces(), and Feel::project().


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

Generated on Fri Oct 25 2013 14:24:27 for Feel++ by doxygen 1.8.4