# 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

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, , where is the second x derivative and is the second y derivative. This has known solution

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.

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 |

We can compare this result to the known solution 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 |

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 on region bounded by unit circle with as the boundary value at radius 1. Just like in the previous example, the solution is known,

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

We can compare this result to the known solution 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 |

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.

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

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

Thanks

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

The code was based on code sample provided in section 3.6 of the book “Compuational Science and Engineering”. The matlab code is also available on this books website.

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

This online book about the finite volume method contains sample matlab code.

http://www.zums.ac.ir/files/research/site/ebooks/mechanics/introductory-finite-volume-methods-for-pdes.pdf

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

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

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/

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

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.

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

thanks

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

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?

i want to create a triangle using equation in matlab.

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?

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?

The 3D plot is made using MATLAB 3D plotting commands like mesh(x,y,z) or surf(…) for a 3D mesh or surface plot. See matlab documentation for details such as

http://www.mathworks.com/help/matlab/ref/mesh.html

and

http://www.mathworks.com/help/matlab/visualize/representing-a-matrix-as-a-surface.html

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.

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?

There is a link provided in this article to a matlab tool called distmesh for generating triangular mesh.

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

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…

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.

https://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).

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

What equation are you solving?

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

How can we change the rectangular shape to an L-shape in matlab and then solve the steady state heat problem by Laplace equation using FEM with boundary conditions?

Hi Hazem

You can create an L-shape using “distmesh” mentioned in the article. Go to the distmesh website and look at the polygon example for how to create mesh for an L-shape.

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

It gives the following example for a polygon.

pv=[-0.4 -0.5;0.4 -0.2;0.4 -0.7;1.5 -0.4;0.9 0.1;

1.6 0.8;0.5 0.5;0.2 1;0.1 0.4;-0.7 0.7;-0.4 -0.5];

[p,t]=distmesh2d(@dpoly,@huniform,0.1,[-1,-1; 2,1],pv,pv);

Then get the list of boundary points using

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

You can set the boundary points to whatever value you want your boundary points to have. For example I set the boundary values to the values of a sine wave in one of the examples in the article.

F(b)=sin(theta.*3);

hi,john,I want creat a crack in a rectangle plate.how can do it by distmesh,please.

Hi James

You can find out how to do it by looking at the examples and documentation on the distmesh site. Look at example 4 on page 10 and 11 of this document.

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

The sample code he gives is

);

(4) Polygons. It is easy to create dpoly (not shown here) for the distance to a

given polygon, using MATLAB’s inpolygon to determine the sign. We mesh a regular

hexagon and fix its six corners:

>> phi=(0:6)’/6*2*pi;

>> pfix=[cos(phi),sin(phi)];

>> [p,t]=distmesh2d(@dpoly,@huniform,0.1,[-1,-1;1,1],pfix,pfix);

Note that pfix is passed twice, first to specify the fixed points, and next as a parameter to dpoly to specify the polygon. In plot (4), we also removed a smaller rotated

hexagon by using ddiff.

He uses ddiff to cut out a shape inside another shape. So you can probably use ddiff with a rectangle and a polygon shape for the crack.

I want to do finite element analysis of composite plates.But since i am a fresher to matlab i am unable to do. Please can somebody provide me simple problems solved in matlab.

john

I find the post extremely helpful to mesh the geometry that I am working in. But connecting this with psi-omega equations of the Navier Stoke’s equations with 2D incompressible fluid flow is difficult for me. Any help?

how to analyze numerical weather prediction using FEM?

Hi,

Thanks for this, it has been very useful.

Just a note, this routine is a LOT faster if you create the sparse K matrix AFTER the element loop, as described here under Method 3: http://imechanica.org/node/7552.

Catherine

Hi,how do you determine which nodes are boundary nodes?

Looking at the distmesh function reference

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

you can get the boundary point list as the unique points returned by boundedges function for the boundary edges.

e = boundedges(p,t);

b = unique(e);

b contains the list of boundary points.

How to calculate the u(x,y)? where (x,y) is not the node.

Look at the P1 Linear Elements mentioned in this post.

https://www.particleincell.com/2012/finite-element-examples/

and the continuous Galerkin method mentioned in the post.

U(x,y) = U1 * phi_1 + U2 * phi_2 + U3 * phi_3

where U1, U2, and U3 are the values at the node and U(x,y) is inside the triangle. The basis functions for the triangle, (phi_1, phi_2, phi_3) are calculated as explained in the P1 Linear Elements section.

When setting the Dirichlet boundary conditions to cart2pol(p(b,:)), I recieve a MATLAB error that the coordinate y is not specified. Do you happen to know why this happening?

Error in ==> cart2pol at 22

th = atan2(y,x);

No I don’t know why it is happening. But instead of

[theta,rho] = cart2pol(p(b,:));

you might try using

[theta,rho] = cart2pol(p(b,1),p(b,2));

where you explicitly provide the y coordinate as p(b,2) . Also try looking at the cart2pol documentation in matlab.

http://www.mathworks.com/help/matlab/ref/cart2pol.html

Hi, May I know if it is possible to adapt the code to neumann boundary condition equals to 0? Sorry I am new to FEM