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
Laplacian with homogeneous Dirichlet conditions

Table of Contents

Author
Feel++ Consortium
Date
2013-02-11




Theory

This part explains how to solve the Laplacian equation for homogeneous dirichlet conditions,

\( \left\{ \begin{aligned} -\Delta u & = f & \text{on}\;\Omega \;, \\ u & = 0 & \text{on}\;\partial\Omega \;,\\ \end{aligned} \right. \)


where \(u\in\Omega\) is the unknown "trial" function and \(\Omega\) the domain.

We multiply each part of the first equation by a "test" function \(v\in H_0^1(\Omega)\) and we integrate the resulting equation on the domain \(\Omega\),

\( \begin{aligned} -\int_\Omega \Delta u v = \int_\Omega f v \;. \end{aligned} \)


We can integrate by parts this equation (Green Theorem) to obtain the variationnal formulation,

\[ \begin{aligned} \int_\Omega \nabla u \nabla v -\underbrace{ \int_{\partial\Omega} \frac{\partial u}{\partial n} v }_{= 0} =\int_\Omega f v \; \end{aligned} \]

where \(n\) denotes a unit outward normal vector to the boundary. We can rewrite the problem as find \(u\in H_0^1(\Omega)\) such that for all \(v\in H_0^1(\Omega)\),

\( \begin{aligned} a(u,v)=l(v) \;, \end{aligned} \)


where \(a\) is a bilinear form, continuous, coercive and \(l\) a linear form.

Implementation

Let's take a look at the Feel++ code (source "doc/manual/tutorial/mylaplacian.cpp").
We consider for this example \(f=1\) constant.

int main(int argc, char**argv )
{
// initialize feel++
using namespace Feel;
Environment env( _argc=argc, _argv=argv,
_about=about(_name="mylaplacian",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
// create mesh
auto mesh = unitSquare();
// function space
auto Vh = Pch<1>( mesh );
auto u = Vh->element();
auto v = Vh->element();
// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
_expr=gradt(u)*trans(grad(v)) );
// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=id(v));
// apply the boundary condition
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=constant(0.) );
// solve the equation a(u,v) = l(v)
a.solve(_rhs=l,_solution=u);
// export results
auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->save();
}

As you can see, the program looks very close to the mathematical formulation.
We use the form2() function to define the bilinear form and form1() for the linear one (see Forms and Solver ).
The gradient for the trial functions is declared with the gradt() expression where as grad() is used for the test functions (see Keywords). Note that we need to transpose the second vector to perform the scalar product.

To introduce the homogeneous dirichlet conditions on the boundary, we use the function on(). Once the variationnal formulation and the boundary conditions are set, we call the solver with solve().