Finite-Element Methods for PDEs

 in

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.

______________________

Joey Bernard has a background in both physics and computer science. This serves him well in his day job as a computational research consultant at the University of New Brunswick. He also teaches computational physics and parallel programming.