FreeFem to Matlab – fast mesh import

I recently wrote a brief introduction to FreeFem++ in this post. FreeFem is a software designed for the numerical study of partial differential equations. It has the advantage of being able to easily define the geometry of the domain, construct and modify meshes, finite element spaces and solve problems on these meshes.

I use Matlab very often for numerical computations. Most of the numerical stuff I’ve done (take a look here if you want) was based on finite differences methods, fundamental solutions and other classical techniques different from finite elements methods. Once I started using finite elements I quickly realized that Matlab is not that easy to work with if you want some automated quality meshing. PDEtool is good, but defining the geometry is not easy. There is also a simple tool: distmesh which performs a simple mesh construction for simple to state geometries. Nevertheless, once you have to define finite element spaces and solve problems things are not easy…

This brings us to the topic of this post: is it possible to interface Matlab and FreeFem? First, why would someone like to do this? Matlab is easier to code and use than FreeFem (for one who’s not a C expert…), but FreeFem deals better with meshes and solving PDE with finite elements. Since FreeFem can be called using system commands, it is possible to call a static program from Matlab. FreeFem can save meshes and variables to files. Let’s see how can we recover them in Matlab.

There is a tool called “FreeFem to Matlab” developed by Julien Dambrine (link on Mathworks). There’s also a longer explanation in this blog post. I recently tried to use the tool and I quickly found that it is not appropriate for large meshes. It probably scans the mesh file line by line which makes the loading process lengthy for high quality meshes. Fortunately there’s a way to speed up things and I present it below. I will not cover the import of the data (other than meshes) since the function importdata from the FreeFem to Matlab tool is fast enough for this.

First, let’s see what’s the structure of the .msh files obtained using the savemesh command in FreeFem. The first line of a .msh file contains three integers representing the number of points, the number of triangles and the boundary edges in the mesh. The following lines contain the coordinates of the points one point by line, the indexes of the triangles one triangle by line and the indexes of the boundary edges again one by line. A useful info is that the boundary segments are already in the good order forming cycles on the mesh.

Once we know the structure of the mesh file we can devise a simple strategy to extract all the information. First we look at the first line which gives information about the whole file (number of point lines, triangle lines and boundary edge lines). To avoid looking at a file line by line which can be time consuming, as a file associated to a large mesh can easily get to lines we can just import all the data in the file using the command importdata in Matlab. This command returns, in this case, a matrix containing all the info. A quick inspection shows that the first line of the matrix contains the integers from the first line of the .msh file. So far so good… Looking further we notice that some of the lines of the matrix contains NaN values. One can quickly see that these lines are in fact separators so we simply remove them in a first step.

After removing the lines containing NaNs we are left with a matrix with exactly the same structure as the .msh file (counts/coordinates/triangles/boundary). From here on the extraction of the mesh characteristics is straightforward, using basic array manipulation in Matlab. The Matlab code is presented below. The time consumption is hugely improved: what took more than a minute before becomes almost instant. This is again a confirmation of the fact that long loops should be avoided in Matlab. Replacing them with operations on arrays increases the speed in non negligible ways.

Matlab code below:

function [p,t,b] = import_mesh(filename)
rawData = importdata(filename); % import the whole file (quick)
row1 = rawData(1,:); % the first row
np   = row1(1); % number of points
nt   = row1(2); % number of triangles

% this removes the lines containing NaNs
rawData(any(isnan(rawData),2),:) = [];

% now we can just extract the remaining infos from the
% rawData array
p = rawData(2:np+1,1:2);  % points
t = rawData(np+2:np+1+nt,:); % triangle (3 indices per row)
b = rawData(np+2+nt:end,1:2); % boundary edges (2 indices per row, cycles)