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:

• Finite elements in 50 lines of MATLAB
• femcode.m

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.

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.

Did you find this article useful? If so, please share the link with your friends, coworkers, and colleagues. And don't forget to subscribe to the newsletter and follow us on Google+ and Twitter. To find out more about how we can assist you on your project, please see the request for proposal page.