Finite-Element Methods for PDEs
One of the common classes of equations that is encountered in several branches of science is partial differential equations. So in this article, I look at a software package called FreeFem++ that is designed to help you calculate these partial differential equations.
One popular method of solving these types of equations, and the one FreeFem++ uses, is the finite-element method. The basic idea with this method is to take the entire problem domain and subdivide it into a mesh of smaller regions. You then apply a simplified version of the partial differential equations that still is locally valid. This makes the problem a tractable one that actually can be solved in a reasonable amount of time.
FreeFem++ is available in several different flavors. The earlier versions were named freefem, freefem3D and freefem+. The latest version is called FreeFem++ and is written in C++. It can be compiled and runs on all three major operating systems, Windows, Mac OS X and Linux. Since this is Linux Journal, I focus on Linux here. On Debian-based distributions, installation is as easy as:
sudo apt-get install freefem++
The other versions also are available as the freefem and freefem3d packages.
Although I am specifically discussing solving electromagnetic fields in this article, FreeFem++ is a general finite-element method solver. This means it should be able to deal with most partial differential equations.
Once it is installed, you will want to start using it. The first thing to be aware of is that FreeFem++ is designed to be used for "production-level" work—meaning active, high-level research work. As such, it does not have a pretty front end to help new users through their first attempts at using it. It is best to think of FreeFem++ as a programming language that you use to write a program to solve your problem.
The first step is to define the geometry of your problem. This step is usually called meshing, or building a mesh. FreeFem++ does not include any CAD functionality to build geometries directly. You can, however, generate meshes automatically from a set of border descriptions. FreeFem++ can take a set of equations describing your geometry and generate a mesh based on the Delaunay-Voronoi algorithm. As a simple example, let's say you wanted to build a square object. You could create the related mesh with these commands:
int sides = 8;
mesh Th = square(sides, sides);
From this first example, you should notice that the language FreeFem++ uses is very similar to C/C++. Variables are typed and can store values only of that type. The language is also polymorphic, so commands and operations will do different things based on the types of the parameters or operands.
You can define more complex shapes with border functions. An example looks like this:
border aa(t=0, 2*pi) {x=cos(t); y=sin(t);}
You then can hand these types of objects in to the
buildmesh()
function
to generate the same type of mesh you would have received from the
higher-level square()
function.
If you want to see what these equations actually
look like (to verify that you have it correct), you can use the
plot()
function to visualize them.
You can hand in the border functions you may
have created or even the mesh object you get from
buildmesh()
.
Once the
mesh is created, you need to create a set of two-dimensional spaces from this
mesh in order to solve your problem. This is done with the
fespace
function:
fespace Vh(Th, P1);
Vh u,v;
The P1
parameter tells
fespace
what type of finite element you want,
whether it is continuous or discontinuous, smooth, linear or with a bubble.
(A rather large set of possibilities is available that will be left as
an exercise for the dear reader.)
Next, you need to define the finite-element functions within this newly generated space—in this case
u
and
v
.
Now that you have all of the background scaffolding built, you need to
define your problem and solve it within your problem geometry. You
can use the problem
type to define a more complex problem. As an example,
you might want to look at the cooling of a hot plate. The following would
set up the problem for you:
mesh Th=square(30,5,[6*x,y]);
fespace Vh(Th, P1);
Vh u=u0, v, uold;
problem thermic(u,v)=int2d(Th)(u*v/dt + k*(dx(u)
↪* dx(v) + dy(u) * dy(v)))
+ int1d(Th,1,3)(alpha*u*v)
- int1d(Th,1,3)(alpha*ue*v)
- int2d(Th)(uold*v/dt) + on(2,4,u=u0);
The variable name thermic
is now a function call. When you issue the
command thermic
, FreeFem++ will go ahead and solve this problem that you
defined. The purpose for this method is to be able to define your problem
and make alterations and adjustments before actually solving it.
If the
problem is simpler to define, you can use the solve
command to define
your problem and do the solving step immediately. For example, if you
wanted to model motion on a membrane, you could use something like this:
solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w)
↪+ dy(phi)*dy(w))
- int2d(Th)(f*w) + on(Gamma1, phi=z);
where the appropriate finite-element space and functions have been defined.
Once FreeFem++ has solved your problem, you can use
plot
with the finite-element functions to visualize the actual results of this numerical
solution.
Although being able to visualize the results of your work immediately is
important, you need to have a way of saving this work so you don't
need to repeat any calculations unnecessarily. You can save meshes
that are generated with the savemesh
function. You simply need to hand in
the mesh to save and a filename to use:
savemesh(Th, "my_mesh.msh");
You can reload this mesh at a later time with the
readmesh
command, for
example:
mesh Sh = readmesh("my_mesh.msh");
Outputting results is a bit more of a hassle. You have access to the
standard C++ input/output streams, specifically cin and cout, so you can
dump out the numerical results that way. You also can create a new output
stream with ofstream
to write things out to a specific file, rather than
to what standard output is redirected to. In this way, you have full
control over what data gets saved, and what format this file and its
data uses.
Now that you've read this introduction to FreeFem++, you should take a look at other tutorials and documentation on the Web. Several good examples are available that should give you at least a starting point to solve the specific problem you are investigating. If the problem you are trying to solve is especially large, FreeFem++ also has MPI support available. In this way, you can spread your calculations over potentially hundreds of CPUs and hopefully get even more work done.