Examples
The strength of the TFEL is to provide a DSL close to the mathematical language typically used to express and describe linear partial differential equations and their solutions.
This page go through several examples, going form the mathematical description of continuous problem to it’s discretization and implementation with the help of the TFEL.
We start from the most basic problem such as the Poisson equation with Dirichlet boundary conditions to more advanced problems related to transport equations and fluid flows.
Poisson problem with Dirichlet boundary conditions
In this example we solve the Poisson equation in a square in the 2D plane. The Poisson equation can be thought as a model for the elastic deformation a thin membrane under a surfacic load. The problem is stated as follow. Let \( \Omega \) be the area occupied by the membrane. For simplicity we well consider a square membrane:
Let \( f:\Omega\to\mathbb R \) a function in \( \mathrm L^2(\Omega) \). The values of the function \( f \) model the force per unit surface applied on the membrane. Then, the vertical displacment of the membrane \( u:\Omega\to\mathbb R \) is a solution of the Poisson equation
For the problem to be well-posed, it is well known that we need additional boundary conditions. Here we consider the homogeneous Dirichlet boundary conditions
where \( \partial\Omega \) is the boundary of the open set \( \Omega \).
We now give the weak formulation of the Poisson problem. We are looking for the vertical displacement \( u \in H_0^1(\Omega) \) such that
for all test function \( v \in H_0^1(\Omega) \). The functional space \( H_0^1(\Omega) \) is the subspace of the Sobolev space \( W^{1, 2} \) whose elements vanish on the boundary of \( \Omega \):
Expanding the intregand, the expression \(\ref{eq:weak-form}\) is equivalent to
We now proceed to the discretization of the weak formulation of the Poisson problem \(\ref{eq:weak-form}\). Let \( \mathcal M_h \) be a conformal triangular mesh of the domain \( \Omega \), where \( h \) is the largest element diameter:
and let \( V_h \) be the function space
We obtain a discrete version of the weak problem with the following Galerkin approximation. The discrete problem is a follow. We look for a function \( u_h \in V_h \) such that
for all test function \( v_h \in V_h \). It is well known that this last problem is well-posed and that the function \( u_h \) is an approximation of \( u \) in the sense that
We know proceed to the implementation of discrete weak formulation
with the help of the TFEL. All the following code should go in a single file,
called poisson.cpp
for example.
We start by including the library header. For clarity, the function
f
is only declared here, and will be defined later in the
source. The constant n
hold the number of subdivision for the
discretization of \(\Omega\). We then open the main
code block.
1
2
3
4
5
6
7
#include <tfel/tfel.hpp>
double f(const double* x);
const std::size_t n = 10;
int main() {
The value of the function f
is the force that we apply on
the membrane, as a function of space coordinates given in the argument
x
. Spaces coordinates are always const double*
in the parameter
list of functions.
We know need to declare the type of mesh that we wish to use, the type of finite element and finite element space that will be used for the discretization of the weak form. In the TFEL, all these model parameters are represented by types and combinaison of types through templates parameters. For clarity, it is a good idea to declare aliases to these new types as follow:
1
2
3
4
5
using cell_type = cell::triangle;
using mesh_type = fe_mesh<cell_type>;
using fe_type = cell_type::fe::lagrange_p1;
using fes_type = finite_element_space<fe_type>;
using quad_type = quad::triangle::qf5pT;
Here we declare cell_type
, the type of cells we use to discretize
the domain \(\Omega\) (triangles), mesh_type
the type of mesh
datatype which is parametrized by the cell type, fe_type
the finite
element we picked (continuous, linear and Lagrange with nodes on the
vertices of the triangles), fes_type
the finite element space type
which describe the space where the solution and test functions live,
and quad_type
the quadrature formula that will be used to integrate
the linear and bilinear forms.
We can then declare a mesh and its boundary submesh:
1
2
mesh_type m(gen_square_mesh(1.0, 1.0, n, n));
submesh<cell_type> dm(m.get_boundary_submesh());
We use the utility function gen_square_mesh
to generate the
structured mesh of a rectangle with side \((1,1)\) into \(2n^2 =
20\) triangles. The variable dm
hold the boundary submesh of m
.
We then allocate a finite element space based on the mesh m
with the
finite element fe_type
, and specify the Dirichlet boundary:
1
2
fes_type fes(m);
fes.add_dirichlet_boundary(dm);
The add_dirichlet_boundary
method of the class fes_type
specify a
homogeneous Dirichlet condition by default, but arbitraty function
defined on \(\partial\Omega\) can be specified.
We can now proceed to the allocation and integration of the bilinear and linear forms:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
bilinear_form<fes_type, fes_type> a(fes, fes); {
auto u(a.get_trial_function());
auto v(a.get_test_function());
a += integrate<quad_type>(d<1>(u) * d<1>(v) +
d<2>(u) * d<2>(v)
, m);
}
linear_form<fes_type> b(fes); {
auto v(b.get_test_function());
b += integrate<quad_type>(f * v, m);
}
On the lines 1 and 10 we allocate bilinear and linear forms
parametrized by the finite element space fes_type
. For convenience
and to avoid name conflicts for the test and trial functions in
particular we wrap the definition of each forms in a code block
{...}
.
Lines 2, 3 and 11 declare instances of the test and trial functions of
the problem, and which are used in the expressions in the integrands of the
integrals that define the linear and bilinear forms. We then use a
natural syntax to express the expression to be integrated on the lines
5-7 and 13. The first argument of the integrate<quad_type>
function
call is an expression to be integrated over the domain discretized by
the mesh m
. This expression is a function of the variables u
(for
the bilinear form) and v
, and must be linear in these variables
(although no test is done at compile time or runtime). Factors and
terms can be composed with the usual C++
operators +
, -
and
*
. Derivatives of the variables u
and v
are expressed with the
function d<position>(u)
where position
is a std::size_t
and
represent the argument of u
to derive with respect to. For exemple
d<1>(u)
is correspond to \(\frac{\partial u}{\partial x_1}\),
where \(x_1\) is the first argument of the function \(u\). You can
compare the expression on the lines 5-7 and 13 with the integrands in
equation \(\ref{eq:weak-form-exp}\).
This notation allow for a high readability of the code and help to
prevent most mistakes in the assembly routines.
Once the two forms are built, we can solve the resulting linear system. So far two linear solver are available: a dense matrix LU solver which call LAPACK routines, and a sparse matrix GMRES solver with an ILU preconditioner implemented with PETSc.
1
2
3
4
5
6
7
8
9
10
dictionary param(dictionary()
.set("maxits", 2000u)
.set("restart", 1000u)
.set("rtol", 1.e-8)
.set("atol", 1.e-50)
.set("dtol", 1.e20)
.set("ilufill", 2u));
solver::petsc::gmres_ilu s(param);
const fes_type::element u(a.solve(b, s));
Here we first declare a dictionary
which keys are the parameter of
the solver declare on line 7. Then we call the solve
method of the
bilinear form with the linear form f
and the solver s
as
parameters.
The solution u
is a finite element function, that is an element of
the finite element space fes
. The solution can eventually be
exported for visualization in the Ensight6 file format and return:
1
2
3
4
5
exporter::ensight6("poisson",
to_mesh_vertex_data<fe_type>(u), "u");
return 0;
}
Since most file format only support data defined on the cell or on the
vertices of a mesh, we first need to convert the finite element
function u
to a mesh_data<double, mesh_type>
. It is done
temporarily on line 2.
The last piece missing is the definition of the function f
. For
simplicity, we choose a constant function \( f(x) = 1\ \forall
x\in\Omega\):
1
2
3
double f(const double* x) {
return 1.0;
}
The full source code corresponding to this example fits in just under 55 lines, and is ready to be copy-pasted here:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include <tfel/tfel.hpp>
double f(const double* x);
const std::size_t n = 10;
int main() {
using cell_type = cell::triangle;
using mesh_type = fe_mesh<cell_type>;
using fe_type = cell_type::fe::lagrange_p1;
using fes_type = finite_element_space<fe_type>;
using quad_type = quad::triangle::qf5pT;
mesh_type m(gen_square_mesh(1.0, 1.0, n, n));
submesh<cell_type> dm(m.get_boundary_submesh());
fes_type fes(m);
fes.add_dirichlet_boundary(dm);
bilinear_form<fes_type, fes_type> a(fes, fes); {
auto u(a.get_trial_function());
auto v(a.get_test_function());
a += integrate<quad_type>(d<1>(u) * d<1>(v) +
d<2>(u) * d<2>(v)
, m);
}
linear_form<fes_type> b(fes); {
auto v(b.get_test_function());
b += integrate<quad_type>(f * v, m);
}
dictionary param(dictionary()
.set("maxits", 2000u)
.set("restart", 1000u)
.set("rtol", 1.e-8)
.set("atol", 1.e-50)
.set("dtol", 1.e20)
.set("ilufill", 2u));
solver::petsc::gmres_ilu s(param);
const fes_type::element u(a.solve(b, s));
exporter::ensight6("poisson",
to_mesh_vertex_data<fe_type>(u), "u");
return 0;
}
double f(const double* x) {
return 1.0;
}
The file poisson.cpp
can be compiled with the command
$ clang++ -std=c++14 poisson.cpp -I$HOME/.local/include -I/usr//local/Cellar/lapack/3.7.1/include/ -I$PETSC_DIR/include -L$HOME/.local/lib -L$PETSC_DIR/lib -L/usr/local/Cellar/lapack/3.7.1/lib -llapacke -ltfel -lpetsc -lmpi -o poisson
(can vary depending on your setup) and run:
$ ./poisson
The result is stored in an Ensight6 case file, which can be visualized
with Paraview. Below is a screenshot of the solution where the
color/x3-coordinate is the value of \(u\), the displacment of the
membrane under the load \(f = 1\) on \(\Omega\).
Transport equation with an SUPG numerical stabilisation method
ToDo
Stokes problem with Dirichlet boundary conditions
ToDo
Navier-Stokes problem with Dirichlet boundary conditions
ToDo