Posts Tagged ‘programming’

FreeFem to Matlab – fast mesh import

October 14, 2016 Leave a 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.

Read more…

ICPC 2015 World Final Problem B

May 21, 2015 Leave a comment

This the the solution to Problem B from the International Collegiate Programming Contest. The list of problems can be found here. The idea is to find the maximal area of the intersection of two moving polygons. The inputs give the initial positions of two convex polygons and the vectors giving their speed in the plane.

One idea of a solution is as follows: take a maximal time {T} and a discretization of the interval {[0,T]} by dividing it into {N} parts. Iterate over the translations at times given by this discretization, compute at each step the area of the intersection of the two polygons (if there is any intersection at all), and in the end find the time {t_0} for which this area is maximized. Now, even if the discretization step {T/N} is large (greater than the demanded precision of {10^{-3}}), we can conclude that {t_0} is an approximation of the final time with an error smaller than {T/N}. This is due to the fact that the function representing the area of the intersection has two monotonicity intervals, as a consequence of the fact that the polygons are convex. On the first interval, the intersection area is increasing, on the second one it is decreasing. Thus, once we have a discrete maximum, we are close to the real maximum.

Now, all we are left to do in order to achieve any desired precision is to refine the search near this discrete maximum, find a new, better approximation, and refine again, until we have enough precision. Of course, one first problem with this algorithm is the initial search using a maximal time {T}. It is possible that if {T} or {N} are not large enough, then we do not detect any intersection of the two polygons. Thus, an initial guess, based on a mathematical argument is needed in order to reduce the search of the optimal time to an interval which is small enough to have enough initial precision to detect the discrete maximum.

The algorithm presented below, uses Matlab’s predefined function polybool which can compute the intersection of two polygons. Of course, this is an overkill, since dealing with intersections of convex polygons is not that complicated (but still, I didn’t have enough time to play with the problem, in order to provide a more optimized version). I do not treat the search for the initial time interval. As I think about it, I guess a n argument based on finding some line intersections should give us a narrow enough time interval (with some care for the case when the two speed directions are collinear). The algorithm presented below solves the sample cases, but could fail in more general situations.

function Prob1_2015(p1,p2)

% choose initial polygons; see the text of the problem
p1 = [6 3 2 2 4 3 6 6 6 7 4 6 2 2 2];
p2 = [4 18 5 22 9 26 5 22 1 -2 1];

%p1 = [4 0 0 0 2 2 2 2 0 1 1];
%p2 = [4 10 0 10 2 12 2 12 0 -1 1];

np1 = p1(1);
np2 = p2(1);
cp1 = p1(2:1+2*np1);
cp2 = p2(2:1+2*np2);

vp1 = p1(end-1:end);
vp2 = p2(end-1:end);

cp1 = cp1(:);
cp1 = reshape(cp1,2,np1);

xp1 = cp1(1,:);
yp1 = cp1(2,:);
cp2 = cp2(:);
cp2 = reshape(cp2,2,np2);

xp2 = cp2(1,:);
yp2 = cp2(2,:);

%set precision
prec = 1;

%here there should be a clever initial choice of 
%the starting time and stopping time
%this choice works well in this case
start = 0;
stop  = 5;

while prec>1e-6
n = 100;

timex = linspace(start,stop,n);
areas = zeros(size(timex));

for i = 1:n
   t = timex(i);
   xt1 = xp1+t*vp1(1);
   yt1 = yp1+t*vp1(2);
   xt2 = xp2+t*vp2(1);
   yt2 = yp2+t*vp2(2);
   [x,y] = polybool('intersection',xt1,yt1,xt2,yt2);
   areas(i) = polyarea(x,y);
[m,I] = max(areas);
if m<1e-6
   % if no polygonal intersection is detected
   % there is nothing further to do

tapp = timex(I);
fprintf('Precision %f | Time %f\n',prec,tapp);
start = tapp-prec;
stop  = tapp+prec;
prec  = prec/10;
%this draws the final position of the polygon
%as well as their intersection
axis equal
hold on
axis equal
hold off

Here’s what you get for the first sample test provided in the questions:polygons3

 >> Prob1_2015
    res = 4.193548

ICPC 2015 World Final Problem A

May 21, 2015 Leave a comment

This is the solution to Problem A from the International Collegiate Programming Contest. The list of problems can be found here. This first problem consists simply of reading the parameters of the function defined below, and computing its values on the set {\{1,2,...,n\}}. Then, you need to find the maximum decrease in the function values.

The inputs are parameters {p,a,b,c,d,n}, the function is

\displaystyle \text{price}(k) =p(\sin(ak+b)+\cos(ck+d)+2).

and the values to be considered are {p(1),p(2),...,p(n)}. Below is a Matlab code which works, at least for the given Sample Inputs. This should be optimized by removing the for loop. It is instant for n=1000 fg=000000 but it should be modified to work until n=10^6 fg=000000.

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;

vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);

mat = zeros(n,1);
for i = 1:n
mat(i) = vals(i)-min(vals(i:end));

res = max(mat)

The sample outputs are:

 >> Prob1_2015(42,1,23,4,8,10)  
 >> Prob1_2015(100,7,615,998,801,3)
 >> Prob1_2015(100,432,406,867,60,1000)

New version which works for large n. It suffices to eliminate the computation of the min at each of the phases of the iteration. This solves the problem in 0.2 seconds for n= 10^6

function res = Prob1_2015(p,a,b,c,d,n)
dis = 1:n;
vals = p*(sin(a*dis+b)+cos(c*dis+d)+2);
mat       = zeros(n,1);
lastmax   = vals(1);
for i = 1:n
  lastmax   = max(lastmax,vals(i));
  mat(i)    = lastmax-vals(i);
res = max(mat)
Categories: matlab, Olympiad Tags: , ,

Solving Poisson’s equation on a general shape using finite differences

October 14, 2014 3 comments

One of the questions I received in the comments to my old post on solving Laplace equation (in fact this is Poisson’s equation) using finite differences was how to apply this procedure on arbitrary domains. It is not possible to use this method directly on general domains. The main problem is the fact that, unlike in the case of a square of rectangular domain, when we have a general shape, the boudary can have any orientation, not only the orientation of the coordinate axes. One way to avoid approach this problem would be using the Finite Element Method. Briefly, you discretize the boundary, you consider a triangulation of the domain with respect to this discretization, then you consider functions which are polynomial and have support in a few number of triangles. Thus the problem is reduced to a finite dimensional one, which can be written as a matrix problem. The implementation is not straightforward, since you need to conceive algorithms for doing the discretization and triangulation of your domain.

One other approach is to consider a rectangle {D} which contains the shape {\Omega} and add a penalization on the exterior of your domain {\Omega}. The problem to solve becomes something like:

\displaystyle (-\Delta +\mu I) u = 1 \text{ on }D

where {\mu} is defined by

\displaystyle \mu(x) = \begin{cases} 0 & x \in \Omega \\ + \infty & x \notin \Omega\end{cases}. \ \ \ \ \ (1)

Note that doing this we do not need to impose the boundary condition {u=0} on {\partial \Omega}. This is already imposed by {\mu}, and the fact that {u} is forced to be zero outside {\Omega}.

Read more…

%d bloggers like this: