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 in a perforated domain
Author
Feel++ Consortium



Problem statement

We solve for the laplacian with homogeneous Dirichlet conditions in a domain with a hole

\[ \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.

The variational formulation reads, find $u \in H^1_0(\Omega)$ such that $\forall v \in H^1_0(\Omega)$

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

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)$,

\[ a(u,v)&=l(v) \;, \]

where $a$ is a bilinear form, continuous, coercive and $l$ a linear form.

Implementation

We defined $\Omega$ as the unit square with a circle inside of radius $0.25$

typedef Mesh<Simplex<2> > mesh_type;
GeoTool::Rectangle R1( 0.05,"R1",GeoTool::Node( 0,0 ),GeoTool::Node( 1,1 ) );
GeoTool::Circle C1( 0.05,"C1",GeoTool::Node( 0.5,0.5 ),GeoTool::Node( 0.75,0.75 ) );
auto R1mesh = R1.createMesh(_mesh=new mesh_type,_name="R1" );
auto C1mesh = C1.createMesh(_mesh=new mesh_type,_name="C1" );
auto R1mC1mesh = ( R1-C1 ).createMesh(_mesh=new mesh_type,_name="R1-C1" );

We consider for this example $f=1$ constant.

auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=id(v));

The complete example is here

using namespace Feel;
Environment env( _argc=argc, _argv=argv,
_desc=feel_options(),
_directory=".",
_about=about(_name="laplacian-with_holes",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
typedef Mesh<Simplex<2> > mesh_type;
GeoTool::Rectangle R1( 0.05,"R1",GeoTool::Node( 0,0 ),GeoTool::Node( 1,1 ) );
GeoTool::Circle C1( 0.05,"C1",GeoTool::Node( 0.5,0.5 ),GeoTool::Node( 0.75,0.75 ) );
auto R1mesh = R1.createMesh(_mesh=new mesh_type,_name="R1" );
auto C1mesh = C1.createMesh(_mesh=new mesh_type,_name="C1" );
auto R1mC1mesh = ( R1-C1 ).createMesh(_mesh=new mesh_type,_name="R1-C1" );
// auto mesh=R1mesh;
// auto mesh=C1mesh;
auto mesh=R1mC1mesh;
auto Vh = Pch<1>( mesh );
auto u = Vh->element();
auto v = Vh->element();
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=id(v));
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
_expr=gradt(u)*trans(grad(v)) );
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=constant(0.) );
a.solve(_rhs=l,_solution=u);
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().

Results

The program is named feelpp_doc_laplacian_with_holes.


Generated on Sun Oct 20 2013 08:25:06 for Feel++ by doxygen 1.8.4