15102017 20 Finite Element Method in MATLAB, 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! You may also be interested in an article on FEM PIC. 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 | 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. | 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 02 0 02 a 08 of 1 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. o 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 © |t then solves Poisson’s equation using the Matlab command U = KF hips ihwmepartlencel zon/20Y2imatab- fend 1015/1072017 20 Finite Element Method in MATLAB © 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, —Uzz — Uy, = 4, where U,z is the second x derivative and U,, is the second y derivative. This has known solution U =1~ 2? — y? 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 stil 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. ARAL LLLP VPS SOTA I NAVAN AVAVAVAVAN Cee SY AiO Figure 2. Triangular mesh of a unit circle from the distmesh tool We created this mesh using the following distmesh commands in MATLAB. Fa=8(p) sart(sum(p."2,2)) 13 [pst] = distmesh2d(fd,@hunsform,2.2,(-1,-151,2],L])5 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], tps siwurpartcteincel.com/2012mattab-fem! 218151012017 20 Finite Element Method in MATLAB, ¢ = boundedges(p,t); unique(e) 5 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. % Fencode2.0 % [p.t,b] from distmesh tool % make sure your matlab path includes the directory where distmesh is installed. Fa=@(p) sqrt(sum(p.*2,2))-25 [pst ]=distmeshad(fd,@huniforn,e.2,[-1,-151,1),(1)5 beunique(boundedges(p,t)); % [KF] = assenble(p,t) % K and F for any mesh of triangles: Linear phi's Nesize(p,1);Tesize(t,1); % number of nodes, number of triangles % p Lists x,y coordinates of N nodes, t lists triangles by 3 node numbers Kesparse(N,N); % zero matrix in sparse format: zeros(N) would be “dense” Fezeros(N,1); % load vector F to hold integrals of phi's tines load (x,y) for e=1:T % integration over one triangular element at a tine nodest(e,: 5 % row of t = node numbers of the 3 corners of triangle Pe=[ones(3,1),p(nodes, :)]; % 3 by 3 matrix with rows=[1 xcorner ycorner] Areacabs(det(Pe))/2; % area of triangle ¢ = half of parallelogram area Ceinv(Pe); % columns of C are coeffs in atbxicy to give phi=1,0,0 at nodes % now compute 3 by 3 Ke and 3 by 1 Fe for elenent © grad=c(2:3,:);Ke=Area*grad' grad; % elenent matrix from slopes b,c in grad Festrea/3*4; % integral of phi over triangle is volume of pyranid: #(x,y)=4 % multiply Fe by 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 assenbled into K and F % [Kb,Fb] = dirichLet(k,F,) % assembled K was singular! K*ones(N1)=@ X Implenent Dirichlet boundary conditions U(b)=@ at nodes in list b K(b, 2): K(b,b)=speye(Length(b),length(b)); % put I into boundary submatrix of K Kbek; FbaF; % Stiffness matrix Kb (sparse format) and load vector Fb 5 K(24b)=0; F(b)=0; % put zeros in boundary rows/colunns of K and F X Solving for the vector U will produce U(b)=9 at boundary nodes UskbFD; X The FEM approximation is Ut phil... + ULN phi_N tps siwuracpartcleincel.com/2012!matlabter 3815102017 20 Finite Element Method in MATLAB, % Plot the FEM approximation U(x,y) with values U1 to ULN at the nodes ‘trisurf(&,p(:,2),p(1,2),0*P(2,1),U, 'edgecolor','k’, 'facecolor’ , "interp"); view(2),axis([-1 1-1 1]),axis equal, colorbar 08 02 25 ° 05 1 Figure 3. Solution of the Poisson's equation on a unit circle We can compare this result to the known solution u — 1—2?—y? 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 retuned by distmesh2d command w= a= plsyt).82 = p(52).82 trisure(t,p(:,2),P(:.2),0*P(:,1).u, 'edgecolor" view(2),axis([+1 1 +1 1]),axis equal, colorbar "“Facecolor’ , “interp'); tps siwuraparcteincel.com/2012mattab-fem! 4815102017 20 Finite Element Method in MATLAB, os 02 mn 05 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- 2) gives L1 norm equal to 0.10568, norm(U-u,2) gives L2 norm equal to 9.215644 nnorm(U-u," if") gives infinity norm equal to 0.073393 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. IMesh Resolution|Node Count|L1 Norm |L2 Norm |L infinity Norm 0.20 ee 0.10568 [0.015644 |0.0073393 0. 162 0.10614 |a.012909 ]o0.0032610 0.10 362 0.083466|0.0073064|0.0016123 Solving 2D Laplace on Unit Circle with nonzero boundary conditions MATLAB Next we will solve Laplaces equation with nonzero dirichlet boundary conditions in 2D using the Finite Element Method. We will solve U;,. + Uy = 0 on region bounded by unit circle with sin(3@) as the boundary value at radius 1. Just like in the previous example, the solution is known, u(r, 0) = r§ sin(30) tps siwuraparcteincel.com/2012mattab-fem! 58ssn0n017 20 Finite Element Metiod in MATLAB 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(30). The following code changes are required The line FeeArea/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 FesArea/3*0; % integral of phi over triangle is volume of pyramid: #(x,y)=@, load is zero in Laplace Also the line for setting the Dirichlet boundary conditions to zero K(b,:)=05 K(:4b)=05 F(b)=95 % put zeros in boundary rows/columns of K and F is changed to K(b,:)=05 [theta,rho] = cart2pol(p(b,:))3 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, 0) = r° sin(3@) 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. tps siwurpartcteincel.com/2012mattab-fem! 68151012017 20 Finite Element Method in MATLAB, [theta,Rho] = cart2pol(p(:,1),p(,2))3 % convert from cartesian to polar coordinates u = Rho.*5.*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 os ° 05 Figure 6. Analytical solution ‘Comparing the Finite Element solution to Laplace equation with known solution produced the following results in MATLAB. norm(U-u,2) gives L1 norm equal to 0.20679 norm(U-u,2) gives L2 norm equal to 2.035458 norm(Usu,"inf') gives infinity norm equal to 6.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 tum 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. Related Articles: Get results faster with Java multithreading tps siwuraparcteincel.com/2012mattab-fem! 781sn012017 20 Finite Element Method in MATLAB, Interpolation using an arbitrary quadrilateral Potential Solver for Composite Dielectrics 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. tps siwuraparcteincel.com/2012mattab-fem! ae