PIC-C, Particle In Cell Consulting, LLC. Plasma physics and rarefied gas simulation codes.

2D Finite Element Method in MATLAB

This guest article was submitted by John Coady (bio below). Would you like to submit an article? If so, please see the submission guidelines. If your article is on scientific computing, plasma modeling, or basic plasma / rarefied gas research, we are interested!

Finite Element Method in Matlab

The Finite Element Method is one of the techniques used for approximating solutions to Laplace or Poisson equations. Searching the web I came across these two implementations of the Finite Element Method written in less than 50 lines of MATLAB code:

The first one of these came with a paper explaining how it worked and the second one was from section 3.6 of the book “Computational Science and Engineering” by Prof. Strang at MIT. I will use the second implementation of the Finite Element Method as a starting point and show how it can be combined with a Mesh Generator to solve Laplace and Poisson equations in 2D on an arbitrary shape. The MATLAB code in femcode.m solves Poisson’s equation on a square shape with a mesh made up of right triangles and a value of zero on the boundary. Running the code in MATLAB produced the following


Figure 1. Solution of the Poisson’s equation on a square mesh using femcode.m

If you look at the Matlab code you will see that it is broken down into the following steps.

  • First it generates a triangular mesh over the region
  • Next it assembles the K matrix and F vector for Poisson’s equation KU=F from each of the triangle elements using a piecewise linear finite element algorithm
  • After that it sets the Dirichlet boundary conditions to zero
  • It then solves Poisson’s equation using the Matlab command U = K\F
  • Finally it plots the results

This particular problem could also have been solved using the Finite Difference Method because of it’s square shape. One of the advantages that the Finite Element Method (and the Finite Volume Method) has over Finite Difference Method is that it can be used to solve Laplace or Poisson over an arbitrary shape including shapes with curved boundaries.

Solving 2D Poisson on Unit Circle with Finite Elements

To show this we will next use the Finite Element Method to solve the following poisson equation over the unit circle, -U_{xx} -U_{yy} =4, where U_{xx} is the second x derivative and U_{yy} is the second y derivative. This has known solution
U=1-x^2 -y^2 which can be verified by plugging this value of U into the equation above giving the result 4 = 4 . Although we know the solution we will still use the Finite Element Method to solve this problem and compare the result to the known solution.The first thing that Finite Elements requires is a mesh for the 2D region bounded by the arbitrary 2D shape. In order to do this we will be using a mesh generation tool implemented in MATLAB called distmesh. You can find out more about the distmesh tool and how to use it here. We will use distmesh to generate the following mesh on the unit circle in MATLAB.


Figure 2. Triangular mesh of a unit circle from the distmesh tool

We created this mesh using the following distmesh commands in MATLAB.

fd=@(p) sqrt(sum(p.^2,2)) -1;
[p,t] = distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);

The values [p,t] returned from the distmesh2d command contain the coordinates of each of the nodes in the mesh and the list of nodes for each triangle. Distmesh also has a command for generating a list of boundary points b from [p,t],

e = boundedges(p,t);
b = unique(e);

The Finite Element method from the first example requires p, t and b as inputs. We will now modify this first example and to use p, t and b generated by distmesh for the region bounded by the unit circle. This can be accomplished by replacing the mesh generation code from the first part of femcode.m with the mesh creation commands from distmesh. The modified code is shown below and produced the result in Figure 3.

% femcode2.m
 
% [p,t,b] from distmesh tool
% make sure your matlab path includes the directory where distmesh is installed.
 
fd=@(p) sqrt(sum(p.^2,2))-1;
[p,t]=distmesh2d(fd,@huniform,0.2,[-1,-1;1,1],[]);
b=unique(boundedges(p,t));
 
 
% [K,F] = assemble(p,t) % K and F for any mesh of triangles: linear phi's
N=size(p,1);T=size(t,1); % number of nodes, number of triangles
% p lists x,y coordinates of N nodes, t lists triangles by 3 node numbers
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense"
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y)
 
for e=1:T  % integration over one triangular element at a time
  nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e
  Pe=[ones(3,1),p(nodes,:)]; % 3 by 3 matrix with rows=[1 xcorner ycorner]
  Area=abs(det(Pe))/2; % area of triangle e = half of parallelogram area
  C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes
  % now compute 3 by 3 Ke and 3 by 1 Fe for element e
  grad=C(2:3,:);Ke=Area*grad'*grad; % element matrix from slopes b,c in grad
  Fe=Area/3*4; % integral of phi over triangle is volume of pyramid: f(x,y)=4
  % multiply Fe by f at centroid for load f(x,y): one-point quadrature!
  % centroid would be mean(p(nodes,:)) = average of 3 node coordinates
  K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K
  F(nodes)=F(nodes)+Fe; % add Fe to 3 components of load vector F
end   % all T element matrices and vectors now assembled into K and F
 
% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb
 
% Solving for the vector U will produce U(b)=0 at boundary nodes
U=Kb\Fb;  % The FEM approximation is U_1 phi_1 + ... + U_N phi_N
 
% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes
trisurf(t,p(:,1),p(:,2),0*p(:,1),U,'edgecolor','k','facecolor','interp');
view(2),axis([-1 1 -1 1]),axis equal,colorbar

Figure 3. Solution of the Poisson’s equation on a unit circle

We can compare this result to the known solution u = 1 - x^2 - y^2 to our poisson equation which is plotted below for comparison. We generated this plot with the following MATLAB commands knowing the list of mesh node points p returned by distmesh2d command.

u = 1 - p(:,1).^2 - p(:,2).^2
trisurf(t,p(:,1),p(:,2),0*p(:,1),u,'edgecolor','k','facecolor','interp');
view(2),axis([-1 1 -1 1]),axis equal,colorbar

Figure 4. The analytical solution

Next we can calculate the difference between the Finite Element approximation and the known solution to the poisson equation on the region bounded by the unit circle. Using MATLAB norm command we can calculate the L1 norm, L2 norm and infinity norm of the difference between approximated and known solution (U – u), where capital U is the Finite Element approximation and lowercase u is the known solution.

norm(U-u,1) gives L1 norm equal to 0.10568
norm(U-u,2) gives L2 norm equal to 0.015644
norm(U-u,'inf') gives infinity norm equal to 0.0073393

When we repeat this experiment using a finer mesh resolution we get the following results with mesh resolution values of 0.2, 0.15 and 0.1 which are passed as a parameter to the distmesh2d command when generating the mesh. As can be seen from the table below, a mesh with a finer resolution results in a Finite Element approximation that is closer to the true solution.

Mesh Resolution Node Count L1 Norm L2 Norm L infinity Norm
0.20 88 0.10568 0.015644 0.0073393
0.15 162 0.10614 0.012909 0.0032610
0.10 362 0.083466 0.0073064 0.0016123

Solving 2D Laplace on Unit Circle with nonzero boundary conditions in MATLAB

Next we will solve Laplaces equation with nonzero dirichlet boundary conditions in 2D using the Finite Element Method. We will solve U_{xx}+U_{yy}=0 on region bounded by unit circle with \sin(3\theta) as the boundary value at radius 1. Just like in the previous example, the solution is known,
u(r,\theta) = r^3 \sin(3\theta)

We will compare this known solution with the approximate solution from Finite Elements. We will be using distmesh to generate the mesh and boundary points from the unit circle. We will modify the MATLAB code to set the load to zero for Laplace’s equation and set the boundary node values to \sin(3\theta). The following code changes are required:

The line

Fe=Area/3*4; % integral of phi over triangle is volume of pyramid: f(x,y)=4, load was 4 in poisson equation

from assembling F is changed to

Fe=Area/3*0; % integral of phi over triangle is volume of pyramid: f(x,y)=0, load is zero in Laplace

Also the line for setting the Dirichlet boundary conditions to zero

K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F

is changed to

K(b,:)=0; [theta,rho] = cart2pol(p(b,:)); F(b)=sin(theta.*3); % set rows in K and F for boundary nodes.

Running the non-zero boundary FEM code produced the result in Figure 5.


Figure 5. FEM solution for non-zero boundary condition

We can compare this result to the known solution u(r,\theta) = r^3 \sin(3\theta) of the Laplace equation with the given boundary conditions which is plotted below for comparison. We generated this plot with the following MATLAB commands given the list of mesh node points p.

[Theta,Rho] = cart2pol(p(:,1),p(:,2)); % convert from cartesian to polar coordinates
u = Rho.^3.*sin(Theta.*3);
trisurf(t,p(:,1),p(:,2),0*p(:,1),u,'edgecolor','k','facecolor','interp');
view(2),axis([-1 1 -1 1]),axis equal,colorbar

Figure 6. Analytical solution

Comparing the Finite Element solution to Laplace equation with known solution produced the following results in MATLAB.

norm(U-u,1) gives L1 norm equal to 0.20679
norm(U-u,2) gives L2 norm equal to 0.035458
norm(U-u,'inf') gives infinity norm equal to 0.011410

Summary

  • The Finite Element Method is a popular technique for computing an approximate solution to a partial differential equation.
  • The MATLAB tool distmesh can be used for generating a mesh of arbitrary shape that in turn can be used as input into the Finite Element Method.
  • The MATLAB implementation of the Finite Element Method in this article used piecewise linear elements that provided a good approximation to the true solution.
  • More accurate approximations can be achieved by using a finer mesh resolution
  • More accurate approximations can also be achieved by using a Finite Element algorithm with quadratic elements or cubic elements instead of piecewise linear elements.
Subscribe to the newsletter and follow us on Twitter. Send us an email if you have any questions.
Bio: John Coady is an electrical engineer with 20 years experience in software development. His interests include scientific computing, computer simulations and web based programming. He enjoys learning new technologies and adopting them for new uses. Examples include scientific computing algorithms applied to the field of computer vision and using the latest web technologies for simulating and visualizing 3D content on the web.

24 comments to “2D Finite Element Method in MATLAB”

  1. chinki
    October 16, 2012 at 1:44 am

    nice article can u please send something on adjoint formulation of spn equations. It would be really helpful…………

    Thanks

  2. Jean Marie
    November 15, 2012 at 1:31 am

    Documents très intéréssents. Je besoin la détail du code MATLAB de la méthode de volumes finis.

  3. daniel leitner
    January 28, 2013 at 8:20 am

    thanks for the code!

    imho for dirichlet boundary conditions (unequal zero) the code should read:

    % Implement Dirichlet boundary conditions U(b)=d at nodes in list b
    K(b,:)=0; F(b)=d; % put zeros in boundary rows of K
    Kb=K+sparse(b,b,ones(size(b)); Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb

  4. sujal padhiyar
    February 13, 2013 at 2:36 am

    I want to create Matlab program of 2 dimentional frame please tell me how to do it i am engineering student

    • john
      February 16, 2013 at 12:20 pm

      I am not sure what you mean by “create Matlab program of 2 dimentional frame”. If you are asking how to create a mesh of different 2 dimentional shapes using matlab then see the distmesh tool for instructions.

      http://persson.berkeley.edu/distmesh/

  5. palekar shailesh
    May 6, 2013 at 12:16 pm

    i am a phd student i have to generate mixed meshing for my problem in which the total plate is meshed by Q4 elements and at the center of plate there is a crack surrounding this crack i have to use triangular elements,how to generate this type of meshing using matlab code,pls tell me how to do this or give me some source to get this type of coding

  6. John
    May 6, 2013 at 1:47 pm

    For the triangular meshing part you can use distmesh tool in matlab. Have a look at the examples and documentation on this site.

    http://persson.berkeley.edu/distmesh/

    Also look at the function reference and published paper.

    http://persson.berkeley.edu/distmesh/funcref.html

    http://persson.berkeley.edu/distmesh/persson04mesh.pdf

    For your problem you can probably use the dpoly function for the crack shape to list the points on the boundary of the crack and calculate the distance from the crack. As you can see in the examples the author puts shapes inside of other shapes to calculate distances for generating a mesh. You can for instance put a polygon shape like a crack inside a rectangle and then generate a triangular shape something similar for what he does in his examples.

    % Example: (Square, with polygon in interior)
    fd=@(p) ddiff(drectangle(p,-1,1,-1,1),dpoly(p,[0.3,0.7; 0.7,0.5]));
    fh=@(p) min(min(0.01+0.3*abs(dcircle(p,0,0,0)), …
    0.025+0.3*abs(dpoly(p,[0.3,0.7; 0.7,0.5]))),0.15);
    [p,t]=distmesh2d(fd,fh,0.01,[0,0;1,1],[0,0;1,0;0,1;1,1]);

    You can also use another dpoly instead of a drectangle for the outer boundary of the 2D distmesh shape.

    For the Q4 elements you can use some other meshing tool to generate Q4 elements.

    http://people.sc.fsu.edu/~jburkardt/m_src/quad_mesh/quad_mesh.html

    Or you can write your own matlab code for this something like the first few lines of the MIT example in femcode.m

    http://math.mit.edu/cse/codes/femcode.m

    You will need to put a shape inside another shape for the region of the Q4 mesh. The inner boundary for the Q4 mesh will need to match the outer boundary of the triangular distmesh mesh so that you can connect these two meshes together. Distmesh has some routines for telling you the boundary points.

    e = boundedges(p,t); % boundary edges
    b = unique(e); % boundary points

    I am not too familiar with meshing software packages but maybe your school has some other software for doing meshing.

  7. johnny
    May 8, 2013 at 2:42 pm

    is there any sample code of FEM by Cosserat model ( micropolar) ?
    thanks

  8. rayehe
    June 24, 2013 at 8:45 am

    I need to mesh a code for meshing a curved plane with triangle,the equation of the plane is specified,in MATLAB
    I’ve surfed the net but found noting
    please help me

  9. ali
    August 17, 2013 at 12:51 am

    I need to solve poisson equation for a rectangle shape that has a sheet charge in the middle of that. how I can change this code for my structure?

  10. jobayer
    October 1, 2013 at 3:55 am

    i want to create a triangle using equation in matlab.

  11. richard
    November 22, 2013 at 4:35 pm

    I need to write a matlab code to generate arbitrary triangular mesh for a rectangle. Does anyone know how to write the code in matlab?

  12. Baljaa
    December 19, 2013 at 3:19 pm

    Hi,
    I have a question about 3D plot of our solution?
    There is 3D plot of solution above, but I cannot find how to make it. Does someone know how to do it?

  13. habib
    July 28, 2014 at 6:17 am

    I am a PhD student in the heat transfer problem I am solving with MATLAB.
    Conductivity of the matrix is equal to the page below. Constant heat source is applied to the page.
    For the boundary conditions given below with the help of finite element software with 20 hexagonal nodal temperature values ​​get resolved.
    • The program is quite flexible so if Pzbr raise a few hundred or a few thousand German German small change can solve the problem (ie just get the number of elements along the x and y to resolve the issue).
    • nodal values ​​of temperature and temperature distribution contour is given as output.
    • Naming variables and explained the program to understand it fully before.

    • asawa
      July 28, 2014 at 6:25 am

      I need to write a matlab code to generate arbitrary triangular mesh for a rectangle. Does anyone know how to write the code in matlab?

  14. John Doe
    August 1, 2014 at 12:42 am

    Dear John,

    I just can’t seem to figure out why this part

    K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K

    appears in your code. Where do these ones come from?
    When I derive the weak form of the given equations, and use Gauss’ Theorem to make it more simple.. I get one integral over the region and another one over the edge.
    So does this identity matrix come from this edge integral? Although I assumed this would be zero…

  15. John
    August 1, 2014 at 8:47 am

    This came from the finite element code in section 3.6 of the book “Computational Science and Engineering” by Prof Strang at MIT.

    http://www-math.mit.edu/cse/

    If you scroll down this page you will see a lot of MATLAB codes for material in the book and the MATLAB code for section 3.6 is the finite element code over a square region.

    I also wrote another blog on Finite Elements using higher order elements which you might find helpful.

    http://www.particleincell.com/blog/2012/finite-element-examples/

    Here is his video lecture for 2D finite element over a circular region.

    http://ocw.mit.edu/courses/mathematics/18-085-computational-science-and-engineering-i-fall-2008/video-lectures/lecture-27-finite-elements-in-2d-part-2/

    He has a number of other lectures on Finite Elements in his course if you want to see how it is all derived.

    That particular line appears in the original code and it look like for boundary nodes in U from KU = F, he previously zeroed out the whole row and column in K for boundary nodes and then he puts back in K the identity for boundary nodes. This is because U contains boundary nodes and 1*U(b) = F(b) is just saying set the boundary nodes to their boundary values F(b).

  16. raju
    August 19, 2014 at 5:52 am

    i need to solve pde b vp by finite element method please help me

    • August 19, 2014 at 8:28 am

      What equation are you solving?

  17. raju
    September 23, 2014 at 8:12 am

    i need matlab code for finding stress temp displacement 2d 3d heat transfer

Leave a Reply

You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> In addition, you can use$latex ...$ to include equations.