Archive

Archive for the ‘matlab’ Category

Mesh a hollow sphere

April 5, 2018 Leave a comment

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).

HollowSphere

Advertisements

Using parfor in Matlab

February 27, 2018 Leave a comment

We all know that loops don’t behave well in Matlab. Whenever it is possible to vectorize the code (i.e. use vectors and matrices to do simultaneous operations, instead of one at a time) significant speed-up is possible. However, there are complex tasks which cannot be vectorized and loops cannot be avoided. I recently needed to compute eigenvalues for some 10 million domains. Since the computations are independent, they could be run in parallel. Fortunately Matlab offers a simple way to do this, using parfor.

There are some basic rules one need to respect to use parfor efficiently:

  1. Don’t use parfor if vectorization is possible. If the task is not vectorizable and computations are independent, then parfor is the way to go.
  2. Variables used in each computation should not overlap between processors. This is for obvious reasons: if two processors try to change the same variable using different values, the computations will be meaningless in the end.
  3. You can use an array or cell to store the results given by each processor, with the restriction that processors should work on disjoint parts of the array, so there is no overlap.

The most restrictive requirement is the fact that one cannot use the same variables in the computations for different processors. In order to do this, the simplest way I found was to use a function for the body of the loop. When using a matlab function, all variables are local, so when running the same function in parallel, the variables won’t overlap, since they are local to each function.

So instead of doing something like

parfor i = 1:N
   commands ...
   array(i) = result
end

you can do the following:

parfor i=1:N
   array(i) = func(i);
end

function res = func(i)
   commands...

This should work very well and no conflict between variables will appear. Make sure to initialize the array before running the parfor, a classical Matlab speedup trick: array = zeros(1,N). Of course, you could have multiple outputs and the output array could be a matrix.

There is another trick to remember if the parpool cannot initialize. It seems that the parallel cluster doesn’t like all the things present in the path sometimes. Before running parfor try the commands

c = parcluster('local');
c.parpool

If you recieve an error, then run

restoredefaultpath
c = parcluster('local');
c.parpool

and add to path just the right folders for your code to work.

Project Euler Problem 285

March 25, 2017 Leave a comment

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.

pb285_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…

pb285_1

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)

Fireworks – Project Euler 317

March 19, 2017 Leave a comment

You can read the text of the problem here. The idea is that we have an explosion at a given height in a uniform gravitational field (no friction/wind). Supposing that all particles go outwards from the explosion point with a constant initial speed, what is the shape of the body formed by these particles? A hint is given in the picture below with a nice Matlab simulation.

euler317

Now what happens if we add a little bit of wind? Here’s the perturbation obtained when adding a 10m/s uniform wind speed vs the initial configuration. Something like this would be a lot more challenging 🙂

pb317_wind

Project Euler Problem 144 – Laser light escaping an ellipse

February 5, 2017 Leave a comment

ellipse

Project Euler – problem 144 – visualization of the solution in Matlab

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.

Read more…

Print your Matlab models in 3D

December 13, 2015 Leave a comment

The emergence of 3D printers opens a whole new level of creation possibilities. Any computer generated model could be materialized as soon as it can be transformed in a language that the 3D printer can use. This is also the case with objects and structures which emerge from various mathematical research topics. Since I’m working on shape optimization problems I have lots of structures that would look nice printed in 3D. Below you can see an example of a 3D model and its physical realization by a 3D printer.

I want to show below how can  you can turn a Matlab coloured patch into a file which can be used by a 3D printer. The first step is to export the Matlab information regarding the position of the points, the face structure and the colours into an obj file format. This is not at all complicated. Vertex information is stored on a line of the form

\displaystyle v\ x\ y\ z\ R\ G\ B

where {v} is exactly the character {v}, {x,y,z} give the coordinates of the points and {R,G,B} give the colour associated to the point in the RGB format. The face information can be entered in a similar fashion:

\displaystyle f\ t_1\ t_2\ t_3

where {t_1, t_2, t_3} are the indices of the points in the corresponding face. Once such an obj file is created, it can be imported in MeshLab (a free mesh editing software). Once you’re in MeshLab you should be able to export the structure into any file format you want, which can be understood by a 3D printer (like STL). Once you have the stl file, you can go on a 3D printing website like Sculpteo and just order your 3D object.

%d bloggers like this: