Re-posted from: http://www.stochasticlifestyle.com/julia-ifem-1-mesh-generation-and-julias-type-system/
My first project on the quest for a Julia finite element method is a simple homework problem. Just some background, this is for UC Irvine’s graduate Computational PDEs 226B course where in the first quarter we did all sorts of finite difference methods and now is our first foray into finite element methods. The purpose of the project is to grasp the data structure enough to use simple tools (i.e. mesh creation and plotting) to create a finite element solver for Poisson’s equation in 2D and check the performance differences. For those who haven’t programmed much this is a great learning exercise, but being a pretty standard exercise the MATLAB code took no time to writeup and so I started thinking: how does Julia compare to MATLAB for solving simple PDEs with the finite element method?
Testing this would take a few steps. The code base is all in MATLAB, so there is a step of (partial) porting MATLAB to Julia. Then there is a fiasco as how to plot the data structures given that our display commands are different (i.e. is there a way to plot in Julia that’s almost exactly like MATLAB?). Then we have to optimize the Julia code and see how it fares against the MATLAB code. This post will focus on how I did the first two steps (and the problems I encountered), the last will be its own post later.
Semi-Porting to Julia: MATLAB.jl and number casting
The moment I looked at Julia code I was struck by how similar it looked to MATLAB code. The first part of the problem is to make a uniform square grid, which is easy enough to be just a few MATLAB commands, so my first response to porting it into Julia was to simply copy-paste it, find out what caused errors, and beat through the errors until it worked. So I downloaded the Julia+Juno IDE bundle and started hacking away at it. The first thing I noticed was the Julia did not have a command for meshgrid in its base library. However, not surprisingly enough people want this command that someone already implemented it and ndgrid. So I took that file and it played as a template for the syntax as well.
With meshgrid in hand, things went really smooth. Array dereferencing in Julia uses square brackets, so I had to go through and change X(:,1) to X[:,1] and the like, but that was easy enough. Another issue is that in MATLAB you can drop parts of the array using
t2nidxMap(topNode) = [];
whereas in Julia I used the deleteat function:
t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode));
You may be asking, what is that peculiar “collect” function? Well in Julia arrays defined by ranges (i.e. A=1:100) return a range object and not an array. In most cases it seems that they can be used just like an array (with better performance and less memory usage), but this seems like a case where they could not, and so the collect function builds the array from the range object. This gave me the code to make the square mesh as:
function squaremesh(square,h) #square = [x0 x1 y0 y1], h is stepsize x0 = square[1]; x1= square[2] y0 = square[3]; y1= square[4] x,y = meshgrid(x0:h:x1,y0:h:y1) node = [x[:] y[:]]; ni = size(x,1); # number of rows N = size(node,1); t2nidxMap = 1:N-ni; topNode = ni:ni:N-ni; t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode)) k = t2nidxMap; elem = [k+ni k+ni+1 k ; k+1 k k+ni+1] return node,elem end
whereas the MATLAB code was
%% Generate nodes x0 = square(1); x1 = square(2); y0 = square(3); y1 = square(4); [x,y] = meshgrid(x0:h:x1,y0:h:y1); node = [x(:),y(:)]; %% Generate elements ni = size(x,1); % number of rows N = size(node,1); t2nidxMap = 1:N-ni; topNode = ni:ni:N-ni; t2nidxMap(topNode) = []; k = (t2nidxMap)'; elem = [k+ni k+ni+1 k; k+1 k k+ni+1];
As you can see, porting that function over took just a few minutes and very few things changed. It would take less than a minute if I knew about the range and deleteat issues. In fact, for next many things I did in Julia, those (and square brackets) were really the only differences, other than the plots.
Plotting via matplotlib
The next part of the project is to plot the triangulation. The MATLAB code for the “showmesh” throws this directly into MATLAB’s trisurf function. Thus I Google’d Julia trisurf and the first hit was Julia’s PyPlot package. It didn’t tell me how to do it, but the package simply calls matplotlib, so I went over to the matplotlib page and and found their function conveniently named plot_trisurf. All I had to do was pick a better color scheme than the default and I was good to go. So while the MATLAB code called
trisurf(elem(:,1:3),node(:,1),node(:,2),zeros(size(node,1),1));
to get the code on the project webpage, the Julia code called
plot_trisurf(node[:,1],node[:,2],zeros(size(node,1)),cmap=get_cmap("ocean"))
to give me this nifty 3D plot which shows the triangalization:
Simple enough.
Getting lazy: MATLAB time?
Okay, things are pretty easy. Now we have to make a circle mesh. So I go into the code only to find that it calls a much larger MATLAB package for mesh generation. This is the part where most people give up on the new language: this part is done in MATLAB and looks like it would take more than the next 15 minutes to port… could I just call MATLAB? If you’ve ever dealt with MEX or similar things that “glue together languages”, you they are a mess. But I found the MATLAB.jl package so I decided to give it a try. The directions are dead simple: just stick the folder from the GitHub repository into the spot they tell you and you’re done. So I fire it up and… it didn’t work. But recall that in MATLAB numbers always promote to the highest state, i.e.
int8(5)+4.4
returns 9. We will return to this later when talking about efficiency, but the problem I was having boiled down this problem. All I had to do was pass in the same array explicitly telling it the numbers were floats by putting a . after it (i.e. 1 vs 1.) and all of the MATLAB codes worked. I was able to mix and matching generating meshes in Julia/MATLAB and plotting them in Julia/MATLAB. An example looks like this. Here I use the iFEM package functions and some of their Julia translations I described before, and mix and match calling them in Julia and MATLAB.
# Testing Square Mesh Generation node,elem = squaremesh([0 1 0 1],.1) showmesh(node,elem); # Testing MATLAB call node2,elem2 = mxcall(:squaremesh,2,[0. 1. 0. 1.],.1) #Array has periods to turn to Float showmesh(node2,elem2) # Testing MATLAB showmesh mxcall(:showmesh,0,node,elem) #surpress return of function handle, cannot be saved STDIO for Julia return mxcall(:showmesh,0,node2,elem2) # Testing Quadralateral showmesh C = .5*zeros(length(node[:,1]),length(node[:,2])) D = zeros(length(node[:,1])) p = pcolormesh(node[:,1],node[:,2],C,edgecolor=D,cmap="ocean") # Doing circle and refine via MATLAB, pyplot plotting node,elem = mxcall(:circlemesh,2,0.,0.,1.,0.2) showmesh(node,elem) node,elem = mxcall(:uniformrefine,2,node,elem) showmesh(node,elem) ## 3D Mesh Generation and Plotting node3,elem3 = mxcall(:cubemesh,2,[0. 1. 0. 1. 0. 1.],.25) showmesh(node3,elem3)
You notice that to do MATLAB calls you simply use mxcall, tell it the function (with a colon in front, in Julia this means you converted it to the symbol type), the number of outputs, and give it the inputs. Be careful about integer vs floating point problems and everything works.
Conclusion: Moving to Julia is dead simple, and you don’t have to move too fast
This really appealed to me in two ways. First of all, Julia’s syntax was very natural to pick up and within minutes I was fine. In fact, the only errors that were hard to diagnose were the ones due to MATLAB! That leads me to the second point, Julia’s language interfacing is so good that I didn’t even have to stray too far. If I wanted to use old MATLAB code, it took 2 minutes to get it setup and I was calling it from Julia. The same seems to be true for calling Python and calling R as well. Not only that, but using C and Fortran code is built right into the core of Julia and has the same syntax structure as mxcall (instead you use ccall), meaning it’s easier to use than MATLAB’s MEX interface. Julia melds together different languages so nicely that you can act like you have the packages of all of them! All of this ease of use comes down Julia’s type system, which I will start looking into in my next post.
So was it worth it? As of right now I can’t say that, but I can say it didn’t cost me more than 30 minutes to get up and running in Julia. But will it be better? We will look at Julia’s type system and performance next time.
The post Julia iFEM1: Porting Mesh Generation appeared first on Stochastic Lifestyle.