Programing the Finite Element Method with Matlab

1 Programing the Finite Element Method withMatlabJack Chessa 3rd October 20021 IntroductionThe goal of this document is to give a very brief overview and directionin the writing of Finite Element code using Matlab . It is assumed that thereader has a basic familiarity with the theory of the Finite Element Method ,and our attention will be mostly on the implementation. An example finiteelement code for analyzing static linear elastic problems written in Matlabis presented to illustrate how to program the Finite Element Method . Theexample program and supporting files are available NotationFor clarity we adopt the following notation in this paper; the bold italics fontvdenotes a vector quantity of dimension equal to the spacial dimension of displacement or velocity at a point, the bold non-italicizedfontddenotes a vector or matrix which is of dimension of the number ofunknowns in the discrete system matrix like the stiffness matrix,an uppercase subscript denotes a node number whereas a lowercase subscriptin general denotes a vector component along a Cartesian unit vector.

2 So, ifdis the system vector of nodal unknowns,uIis a displacement vector of nodeIanduIiis the component of the displacement at nodeIin theidirection, oruI ei. Often Matlab syntax will be intermixed with mathematical notation Graduate Research Assistant, Northwestern University hopefully adds clarity to the explanation. The typewriter font,font,is used to indicate that Matlab syntax is being A Few Words on Writing Matlab ProgramsThe Matlab programming language is useful in illustrating how to programthe Finite Element Method due to the fact it allows one to very quickly codenumerical methods and has a vast predefined mathematical library. This isalso due to the fact that matrix (sparse and dense), vector and many linearalgebra tools are already defined and the developer can focus entirely onthe implementation of the algorithm not defining these data structures.

3 Theextensive mathematics and graphics functions further free the developer fromthe drudgery of developing these functions themselves or finding equivalentpre-existing libraries. A simple two dimensional Finite Element program inMatlab need only be a few hundred lines of code whereas in Fortran or C++one might need a few the Matlab programming language is very complete with re-spect to it s mathematical functions there are a few Finite Element specifictasks that are helpful to develop as separate functions. These have beenprogramed and are available at the previously mentioned web usual there is a trade off to this ease of development. Since Matlabis an interpretive language; each line of code is interpreted by the Matlabcommand line interpreter and executed sequentially at run time, the runtimes can be much greater than that of compiled programming languageslike Fortran or C++.

4 It should be noted that the built-in Matlab functionsare already compiled and are extremely efficient and should be used as muchas possible. Keeping this slow down due to the interpretive nature of Matlabin mind, one programming construct that should be avoided at all costs is thefor loop, especially nested for loops since these can make a Matlab programsrun time orders of magnitude longer than may be needed. Often for loopscan be eliminated using Matlab s vectorized addressing. For example, thefollowing Matlab code which sets the row and column of a matrixAto zeroand puts one on the diagonalfor i=1:size(A,2)A(n,i)=0;endfor i=1:size(A,1)A(i,n)=0;end2A(n,n)=1;shoul d never be used since the following codeA(:,n)=0;A(:,n)=0;A(n,n)=0;does that same in three interpreted lines as opposed tonr+nc+1 interpretedlines, whereAis anr ncdimensional matrix.

5 One can easily see that this canquickly add significant overhead when dealing with large systems (as is oftenthe case with Finite Element codes). Sometimes for loops are unavoidable,but it is surprising how few times this is the case. It is suggested that afterdeveloping a Matlab program, one go back and see how/if they can eliminateany of the for loops. with practice this will become second Sections of a Typical Finite Element Pro-gramA typical Finite Element program consists of the following sections1. Preprocessing section2. Processing section3. Post-processing sectionIn the preprocessing section the data and structures that define the particularproblem statement are defined. These include the Finite Element discretiza-tion, material properties, solution The processing section iswhere the Finite Element matrices, force , boundary conditions are enforced and the system is solved.

6 Thepost-processing section is where the results from the processing section areanalyzed. Here stresses may be calculated and data might be visualized. Inthis document we will be primarily concerned with the processing pre and post-processing operations are already programmed in Matlaband are included in the online reference; if interested one can either look di-rectly at the Matlab script files or typehelp function name at the Matlabcommand line to get further information on how to use these Finite Element Data Structures in MatlabHere we discuss the data structures used in the Finite Element Method andspecifically those that are implemented in the example code. These are some-what arbitrary in that one can imagine numerous ways to store the data fora Finite Element program, but we attempt to use structures that are the mostflexible and conducive to Matlab .

7 The design of these data structures may bedepend on the programming language used, but usually are not significantlydifferent than those outlined Nodal Coordinate MatrixSince we are programming the Finite Element Method it is not unexpected thatwe need some way of representing the Element discretization of the do so we define a set of nodes and a set of elements that connect thesenodes in some way. The node coordinates are stored in the nodal coordinatematrix. This is simply a matrix of the nodal coordinates (imagine that).The dimension of this matrix isnn sdimwherennis the number of nodesandsdimis the number of spacial dimensions of the problem. So, if weconsider a nodal coordinate matrixnodesthe y-coordinate of thenthnode isnodes(n,2). Figure 1 shows a simple Finite Element discretization.

8 For thissimple mesh the nodal coordinate matrix would be as follows:nodes= .(1) Element Connectivity MatrixThe Element definitions are stored in the Element connectivity matrix. Thisis a matrix of node numbers where each row of the matrix contains the con-nectivity of an Element . So if we consider the connectivity matrixelementsthat describes a mesh of 4-node quadrilaterals the 36th Element is definedby the connectivity vectorelements(36,:)which for example may be[ 3642 13 14]or that the elements connects nodes 36 42 13 14. So for4the simple mesh in Figure 1 the Element connectivity matrix iselements= 1 2 32 4 34 5 26 5 4 .(2)Note that the elements connectivities are all ordered in a counter-clockwisefashion; if this is not done so some Jacobian s will be negative and thus cancause the stiffnesses matrix to be singular (and obviously wrong!)

9 !!). Definition of BoundariesIn the Finite Element Method boundary conditions are used to either formforce vectors (natural or Neumann boundary conditions) or to specify thevalue of the unknown field on a boundary (essential or Dirichlet boundaryconditions). In either case a definition of the boundary is needed. The mostversatile way of accomplishing this is to keep a Finite Element discretizationof the necessary boundaries. The dimension of this mesh will be one orderless that the spacial dimension of the problem ( 2D boundary mesh fora 3D problem, 1D boundary mesh for a 2D problemetc.). Once again let sconsider the simple mesh in Figure 1. Suppose we wish to apply a boundarycondition on the right edge of the mesh then the boundary mesh would be thedefined by the following Element connectivity matrix of 2-node line elementsright Edge=[2 44 6].

10 (3)Note that the numbers in the boundary connectivity matrix refer to the samenode coordinate matrix as do the numbers in the connectivity matrix of theinterior elements . If we wish to apply an essential boundary conditions onthis edge we need a list of the node numbers on the edge. This can be easilydone in Matlab with = unique(rightEdge);This will set the vectornodesOnBoundaryequal to [2 4 6]. If we wish tofrom a force vector from a natural boundary condition on this edge we simplyloop over the elements and integrate the force on the edge just as we wouldintegrate any Finite Element operators on the domain Dof MappingUltimately for all Finite Element programs we solve a linear algebraic systemof the formKd=f(4)for the vectord. The vectordcontains the nodal unknowns for that definethe Finite Element approximationuh(x) =nn I=1NI(x)dI(5)whereNI(x) are the Finite Element shape functions,dIare the nodal un-knowns for the nodeIwhich may be scalar or vector quantities (ifuh(x) isa scalar or vector) andnnis the number of nodes in the discretization.