### Archive

Archive for the ‘Numerical Analysis’ Category

## FreeFem++ Tutorial – Part 2

Here are a few tricks, once you know the basics of FreeFem. If your plan is straightforward: define the domain, build the mesh, define the problem, solve the problem, plot the result… then things are rather easy. If you want to do stuff in a loop, like an optimization problem, things may get more complicated, but FreeFem still has lots of tricks up its sleeves. I’ll go through some of them.

1. Defining the geometry of the domain using border might be tricky at the beginning. Keep in mind that the domain should be on the left side of the curves definining the boundaries. This could be acheived in the parametrization chosen or you could reverse the parametrization of a particular part of the domain in the following way. Suppose you have the pieces of boundary C1,C2,C3,C4. You wish to build a mesh with the command mesh Th = buildmesh(C1(100)+C2(100)+C3(100)+C4(100));but you get an error concerning a bad sens on one of the boundaries. If you identify, for example, that you need to change the orientation of C3, you can acheive this by changing the sign of the integer defining the number of points on that part of the boundary: mesh Th = buildmesh(C1(100)+C2(100)+C3(-100)+C4(100));You could make sure that you have the right orientations and good connectivities for all the boundaries by running a plot command before trying mesh: something like

plot(C1(100)+C2(100)+C3(-100)+C4(100));

should produce a graph of all your boundaries with arrows showing the orientations. This is good as a debug tool when you don’t know where the error comes from when you define the domain.

2. Keep in mind that there are also other ways to build a mesh, like square, which meshes a rectangular domain with quite a few options and trunc which truncates or modifies a mesh following some criteria. Take a look in the documentation to see all the options.
3. Let’s say that you need to build a complex boundary, but with many components with similar properties. Two examples come to mind: polygons and domains with many circular holes. In order to do this keep in mind that it is possible to define a some kind of “vectorial boundary”. Let’s say that you want to mesh a polygon and you have the coordinates stored in the arrays xs,ys. Furthermore, you have another array of integers ind which point to the index of the next vertex. Then the boundary of the polygon could be defined with the following syntax: border poly(t=0,1; i){ x=(1-t)*xx[i]+t*xx[ind(i)]; y=(1-t)*yy[i]+t*yy[ind(i)]; label=i; } Now the mesh could be constructed using a vector of integers NC containing the number of desired points on each of the sides of the polygon:
 mesh Th = buildmesh (poly(NC)); 
4. You can change a 2D mesh using the command adaptmesh. There are various options and you’ll need to search in the documentation for a complete list. I use it in order to improve a mesh build with buildmesh. An example of command which gives good results in some of my codes is:
 Th = adaptmesh(Th,0.02,IsMetric=1,nbvx=30000); 
The parameters are related to the size of the triangles, the geometric properties of the griangles and the maximal number of vertices you want in your mesh. Experimenting a bit might give you a better idea of how this command works in practice.
5. There are two ways of defining the problems in FreeFem. One is with solve and the other one is with problem. If you use solve then FreeFem solves the problem where it is defined. If you use problem then FreeFem remembers the problem as a variable and will solve it whenever you call this variable. This is useful when solving the same problem multiple times. Note that you can modify the coefficients of the PDE and FreeFem will build the new problem with the updated coefficients.
6. In a following post I’ll talk about simplifying your code using macros and functions. What is good to keep in mind is that macros are verbatim code replacements, which are quite useful when dealing with complex formulas in your problem definition or elsewhere. Functions allow you to run a part of the code with various parameters (just like functions in other languages like Matlab).

I’ll finish with a code which computes the eigenvalue of the Laplace operator on a square domain with multiple holes. Try and figure out what the commands and parameters do.

int N=5; //number of holes
int k=1; //number of eigenvalue

int nb = 100;  // parameter for mesh size
verbosity = 10; // parameter for the infos FreeFem gives back

real delta = 0.1;
int nbd = floor(1.0*nb/N*delta*2*pi);

int bsquare = 0;
int bdisk   = 1;

// vertices of the squares
real[int] xs(N^2),ys(N^2);

real[int] initx = 0:(N^2-1);
for(int i=0;i&lt;N^2;i++){
xs[i] = floor(i/N)+0.5;
ys[i] = (i%N)+0.5;
}

cout &lt;&lt; xs &lt;&lt; endl;
cout &lt;&lt; ys &lt;&lt; endl;

int[int] Nd(N^2);

Nd = 1;
Nd = -nbd*Nd;

cout &lt;&lt; Nd &lt;&lt; endl;

// sides of the square
border Cd(t=0,N){x = t;y=0;label=bsquare;}
border Cr(t=0,N){x = N;y=t;label=bsquare;}
border Cu(t=0,N){x = N-t;y=N;label=bsquare;}
border Cl(t=0,N){x = 0;y=N-t;label=bsquare;}

border disks(t=0,2*pi; i){
x=xs[i]+delta*cos(t);
y=ys[i]+delta*sin(t);
label=bdisk;
}

plot(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

mesh Th = buildmesh(Cd(nb)+Cr(nb)+Cu(nb)+Cl(nb)+disks(Nd));

plot(Th);

int[int] bc = [0,1]; // Dirichlet boundary conditions

fespace Vh(Th,P2);
// variables on the mesh
Vh u1,u2;

// Define the problem in weak form
varf a(u1,u2) = int2d(Th)
(dx(u1)*dx(u2) + dy(u1)*dy(u2))+on(1,u1=0)//on(C1,C2,C3,C4,u1=1)
+on(bc,u1=0);
varf b([u1],[u2]) = int2d(Th)(  u1*u2 ) ;
// define matrices for the eigenvalue problem
matrix A= a(Vh,Vh,solver=Crout,factorize=1);
matrix B= b(Vh,Vh,solver=CG,eps=1e-20);

// we are interested only in the first eigenvalue
int eigCount = k;
real[int] ev(eigCount); // Holds eigenvalues
Vh[int] eV(eigCount);   // holds eigenfunctions
// Solve Ax=lBx
int numEigs = EigenValue(A,B,sym=true,sigma=0,value=ev,vector=eV);

cout &lt;&lt; ev[k-1] &lt;<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;"></span>&lt; endl;

plot(eV[k-1],fill=1,nbiso=50,value=1);


Here is the mesh and the solution given by the above program:

## Mesh a hollow sphere

This is an interesting experiment I’ve done today and I’d like to share it. There are nice software out there which allow you to mesh regular surfaces, like the sphere, torus or other examples.

I needed something different: I wanted to mesh a sphere with a hole and an inner surface. This kind of surfaces may be useful for those who want to design objects for 3D printing. When using 3D printers, the cost is usually proportional to the volume of material used, since this also consumes impression time. If you manage to make a hollow object, the cost goes down immediately (as the resistance of the object…). However, to make a hole, you need to consider an inner parallel surface. I’ll not handle this here, since it’s more complicated. I’ll just mesh a hollow sphere.

I did this in Matlab using the nice library Distmesh, which is free and really easy to use. The whole difficulty was making a function whose zero level set is the hollow sphere. Recall that it is possible to find the union of two shapes using the minimum of two level set function and the intersection using the maximum. This is all that you need to understand what I did below. You can see the details in the Matlab code below.

function [p,t] = Hole_Spher_Mesh(mh)

R = 0.5;
r = 0.4;
rh = 0.1;

fh=@(p) sqrt(min(p(:,1).^2+p(:,2).^2,ones(size(p(:,1)))));
fh=@(p) sqrt((p(:,1).^2+p(:,2).^2));
[p,t]=distmeshsurface(@(x)fd(x,R,r,rh),@huniform,mh,[-1.1,-1.1,-1.1;1.1,1.1,1.1]);
points = p';
npt = prod(size(points))/3

clf
patch('Faces',t,'Vertices',p,'FaceColor',[0,0.7,0.9],'EdgeColor','none','FaceAlpha',0.5);

function res = fd(p,R,r,rh)

res1 = sqrt(sum(p.^2,2))-R;
res2 = (rh-sqrt(p(:,1).^2+p(:,2).^2));
res3 = p(:,3);
res2 = max(-res2,res3);
res4 = sqrt(sum(p.^2,2))-r;
res2 = min(res2,res4);
res = max(res1,-res2);


And here is the resulting surface for the parameter mh = 0.01 (if you put a higher precision computations will take longer).

## Project Euler Problem 285

Another quite nice problem from Project Euler is number 285. The result of the problem depends on the computation of a certain probability, which in turn is related to the computation of a certain area. Below is an illustration of the computation for k equal to 10.

To save you some time, here’s a picture of the case k=1 which I ignored and spent quite a lot of time debugging… Plus, it only affects the last three digits or so after the decimal point…

Here’s a Matlab code which can construct the pictures above and can compute the result for low cases. To solve the problem, you should compute explicitly all these areas.


function problem285(k)

N = 100000;

a = rand(1,N);
b = rand(1,N);

ind = find(abs(sqrt((k*a+1).^2+(k*b+1).^2)-k)<0.5);

plot(a(ind),b(ind),'.');
axis equal

M = k;
pl = 1;

for k=1:M
if mod(k,100)==0
k
end
r1 = (k+0.5)/k;
r2 = (k-0.5)/k;

f1 = @(x) (x<=(-1/k+r1)).*(x>=(-1/k-r1)).*(sqrt(r1^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f1 = @(x) f1(x).*(f1(x)>=0);
f2 = @(x) (x<=(-1/k+r2)).*(x>=(-1/k-r2)).*(sqrt(r2^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f2 = @(x) f2(x).*(f2(x)>=0);

if k == pl
thetas = linspace(0,pi/2,200);
hold on
plot(-1/k+r1*cos(thetas),-1/k+r1*sin(thetas),'r','LineWidth',2);
plot(-1/k+r2*cos(thetas),-1/k+r2*sin(thetas),'r','LineWidth',2);
plot([0 1 1 0 0],[0 0 1 1 0],'k','LineWidth',3);
hold off
axis off
end

A(k) = integral(@(x) f1(x)-f2(x),0,1);

end

xs = xlim;
ys = ylim;

w = 0.01;
axis([xs(1)-w xs(2)+w ys(1)-w ys(2)+w]);

sum((1:k).*A)



## FreeFem to Matlab – fast mesh import

October 14, 2016 1 comment

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.

## FreeFem++ Tutorial – Part 1

First of all, FreeFem is a numerical computing software which allows a fast and automatized treatment of a variety of problems related to partial differential equations. Its name, FreeFem, speaks for itself: it is free and it uses the finite element method. Here are a few reasons for which you may choose to use FreeFem for a certain task:

1. It allows the user to easily define 2D (and 3D) geometries and it does all the work regarding the construction of meshes on these domains.
2. The problems you want to solve can be easily written in the program once we know their weak forms.
3. Once we have variables defined on meshes or solutions to some PDE, we can easily compute all sorts of quantities like integral energies, etc.

Before showing a first example, you need to install FreeFem. If you are not familiar with command line work or you just want to get to work, like me, you can install the visual version of FreeFem which is available here. Of course, you can find example programs in the FreeFem manual or by making a brief search on the internet.

I’ll present some basic stuff, which will allow us in the end to solve the Laplace equation in a circular domain. Once we have the structure of the program, it is possible to change the shape of the domain in no time.

## Even small errors can be fatal

Machines and computers represent numbers as sequences of zeros and ones called bits. The reason for doing this is the simplicity of constructing circuits dealing with two states. This fact coupled with limited memory capacities means that from the start we cannot represent all numbers in machine code. It is true that we can get as close as we want to any number with great memory cost when precision is important, but in fact we always use a fixed precision. In applications this precision is fixed to a number of bits (16, 32, 64) which correspond to significant digits in computations. Doing math operations to numbers represented on bits may lead to loss of information. Consider the following addition using ${15}$ significant digits:

$\displaystyle 5.00000000000002+6.00000000000003 = 11.0000000000000,$

notice that using only the first ${15}$ significant digits we have made an error of ${5\times 10^{-14}}$. This may seem small, but if we do not pay attention and we let errors of this kind accumulate we may have a final error which is unacceptable. This is what happened in the case of the Patriot missile case which I’ll discuss briefly below.

## Identifying edges and boundary points – 2D Mesh – Matlab

April 21, 2015 1 comment

A triangulation algorithm often gives as output a list of points, and a list of triangle. Each triangle is represented by the indexes of the points which form it. Sometimes we need extra information on the structure of the triangulation, like the list of edges, or the list of boundary points. I will present below two fast algorithms for doing this.

Finding the list of edges is not hard. The idea is to go through each triangle, and extract all edges. The algorithm proposed below creates the adjacency matrix in the following way: for each triangle ${[i,j,k]}$ we set the elements ${a_{ij},a_{jk},a_{ik}}$ (and their symmetric correspondents) to be equal to ${1}$.

In order to find the points on the boundary (in two dimensions), it is enough to look for the edges which are sides to only one triangle. We can do this using the adjacency matrix. Note that if ${A}$ is the adjacency matrix, then ${A^2=(b_{ik})}$ stores the number of paths of length ${2}$ (two sides) between two points of the triangulation. Note that any edge which is not on the boundary will contain the starting and ending point of two paths of length ${2}$. If ${[i,j,k]}$ is a triangle such that points ${i,j}$ are on the boundary, then ${b_{i,j}=1}$ (there is one path of length ${2}$ going through ${i,k,j}$. We also have ${a_{i,j} = 1}$. Conversely, if ${a_{i,j} = 1}$ and ${b_{i,j} = 1}$ then ${i,j}$ are connected, and there is a unique path of length ${2}$ going from ${i}$ to ${j}$. Thus, ${i,j}$ is an edge on the boundary. Therefore, we just need to identify the indexes ${i}$ such that there exists ${j}$ with ${a_{i,j} b_{i,j} = 1}$.

Below are two short Matlab codes doing these two algorithms. I guess they are close to being optimal, since only sparse and vectorized operations are used.


%p is the list of points
%T is the list of triangles, ap is the number of points