Programming The Finite Element Method With Matlab.pdf-资料库

weissb-10689171-4744302543419133179.pdf-第1页.png

第1页 / 共45页

weissb-10689171-4744302543419133179.pdf-第2页.png

第2页 / 共45页

weissb-10689171-4744302543419133179.pdf-第3页.png

第3页 / 共45页

weissb-10689171-4744302543419133179.pdf-第4页.png

第4页 / 共45页

weissb-10689171-4744302543419133179.pdf-第5页.png

第5页 / 共45页

weissb-10689171-4744302543419133179.pdf-第6页.png

第6页 / 共45页

weissb-10689171-4744302543419133179.pdf-第7页.png

第7页 / 共45页

weissb-10689171-4744302543419133179.pdf-第8页.png

第8页 / 共45页

Programing the Finite Element Method with Matlab Jack Chessa ∗ 3rd October 2002 1 Introduction The goal of this document is to give a very brief overview and direction in the writing of ﬁnite element code using Matlab. It is assumed that the reader has a basic familiarity with the theory of the ﬁnite element method, and our attention will be mostly on the implementation. An example ﬁnite element code for analyzing static linear elastic problems written in Matlab is presented to illustrate how to program the ﬁnite element method. The example program and supporting ﬁles are available at http://www.tam.northwestern.edu/jfc795/Matlab/ 1.1 Notation For clarity we adopt the following notation in this paper; the bold italics font v denotes a vector quantity of dimension equal to the spacial dimension of the problem i.e. the displacement or velocity at a point, the bold non-italicized font d denotes a vector or matrix which is of dimension of the number of unknowns in the discrete system i.e. a system matrix like the stiﬀness matrix, an uppercase subscript denotes a node number whereas a lowercase subscript in general denotes a vector component along a Cartesian unit vector. So, if d is the system vector of nodal unknowns, uI is a displacement vector of node I and uIi is the component of the displacement at node I in the i direction, or uI · ei. Often Matlab syntax will be intermixed with mathematical notation ∗Graduate Research Assistant, Northwestern University (j-chessa@northwestern.edu) 1

which hopefully adds clarity to the explanation. The typewriter font, font, is used to indicate that Matlab syntax is being employed. 2 A Few Words on Writing Matlab Programs The Matlab programming language is useful in illustrating how to program the ﬁnite element method due to the fact it allows one to very quickly code numerical methods and has a vast predeﬁned mathematical library. This is also due to the fact that matrix (sparse and dense), vector and many linear algebra tools are already deﬁned and the developer can focus entirely on the implementation of the algorithm not deﬁning these data structures. The extensive mathematics and graphics functions further free the developer from the drudgery of developing these functions themselves or ﬁnding equivalent pre-existing libraries. A simple two dimensional ﬁnite element program in Matlab need only be a few hundred lines of code whereas in Fortran or C++ one might need a few thousand. Although the Matlab programming language is very complete with re- spect to it’s mathematical functions there are a few ﬁnite element speciﬁc tasks that are helpful to develop as separate functions. These have been programed and are available at the previously mentioned web site. As usual there is a trade oﬀ to this ease of development. Since Matlab is an interpretive language; each line of code is interpreted by the Matlab command line interpreter and executed sequentially at run time, the run times can be much greater than that of compiled programming languages like Fortran or C++. It should be noted that the built-in Matlab functions are already compiled and are extremely eﬃcient and should be used as much as possible. Keeping this slow down due to the interpretive nature of Matlab in mind, one programming construct that should be avoided at all costs is the for loop, especially nested for loops since these can make a Matlab programs run time orders of magnitude longer than may be needed. Often for loops can be eliminated using Matlab’s vectorized addressing. For example, the following Matlab code which sets the row and column of a matrix A to zero and puts one on the diagonal for i=1:size(A,2) A(n,i)=0; end for i=1:size(A,1) A(i,n)=0; end 2

A(n,n)=1; should never be used since the following code A(:,n)=0; A(:,n)=0; A(n,n)=0; does that same in three interpreted lines as opposed to nr +nc+1 interpreted lines, where A is a nr×nc dimensional matrix. One can easily see that this can quickly add signiﬁcant overhead when dealing with large systems (as is often the case with ﬁnite element codes). Sometimes for loops are unavoidable, but it is surprising how few times this is the case. It is suggested that after developing a Matlab program, one go back and see how/if they can eliminate any of the for loops. With practice this will become second nature. 3 Sections of a Typical Finite Element Pro- gram A typical ﬁnite element program consists of the following sections 1. Preprocessing section 2. Processing section 3. Post-processing section In the preprocessing section the data and structures that deﬁne the particular problem statement are deﬁned. These include the ﬁnite element discretiza- tion, material properties, solution parameters etc. . The processing section is where the ﬁnite element objects i.e. stiﬀness matrices, force vectors etc. are computed, boundary conditions are enforced and the system is solved. The post-processing section is where the results from the processing section are analyzed. Here stresses may be calculated and data might be visualized. In this document we will be primarily concerned with the processing section. Many pre and post-processing operations are already programmed in Matlab and are included in the online reference; if interested one can either look di- rectly at the Matlab script ﬁles or type help ’function name’ at the Matlab command line to get further information on how to use these functions. 3

4 Finite Element Data Structures in Matlab Here we discuss the data structures used in the ﬁnite element method and speciﬁcally those that are implemented in the example code. These are some- what arbitrary in that one can imagine numerous ways to store the data for a ﬁnite element program, but we attempt to use structures that are the most ﬂexible and conducive to Matlab. The design of these data structures may be depend on the programming language used, but usually are not signiﬁcantly diﬀerent than those outlined here. 4.1 Nodal Coordinate Matrix Since we are programming the ﬁnite element method it is not unexpected that we need some way of representing the element discretization of the domain. To do so we deﬁne a set of nodes and a set of elements that connect these nodes in some way. The node coordinates are stored in the nodal coordinate matrix. This is simply a matrix of the nodal coordinates (imagine that). The dimension of this matrix is nn × sdim where nn is the number of nodes and sdim is the number of spacial dimensions of the problem. So, if we consider a nodal coordinate matrix nodes the y-coordinate of the nth node is nodes(n,2). Figure 1 shows a simple ﬁnite element discretization. For this simple mesh the nodal coordinate matrix would be as follows:   . 0.0 0.0 2.0 0.0 0.0 3.0 2.0 3.0 0.0 6.0 2.0 6.0 nodes = (1) 4.2 Element Connectivity Matrix The element deﬁnitions are stored in the element connectivity matrix. This is a matrix of node numbers where each row of the matrix contains the con- nectivity of an element. So if we consider the connectivity matrix elements that describes a mesh of 4-node quadrilaterals the 36th element is deﬁned by the connectivity vector elements(36,:) which for example may be [ 36 42 13 14] or that the elements connects nodes 36 → 42 → 13 → 14. So for 4

the simple mesh in Figure 1 the element connectivity matrix is   . 1 2 3 2 4 3 4 5 2 6 5 4 elements = (2) Note that the elements connectivities are all ordered in a counter-clockwise fashion; if this is not done so some Jacobian’s will be negative and thus can cause the stiﬀnesses matrix to be singular (and obviously wrong!!!). 4.3 Deﬁnition of Boundaries In the ﬁnite element method boundary conditions are used to either form force vectors (natural or Neumann boundary conditions) or to specify the value of the unknown ﬁeld on a boundary (essential or Dirichlet boundary conditions). In either case a deﬁnition of the boundary is needed. The most versatile way of accomplishing this is to keep a ﬁnite element discretization of the necessary boundaries. The dimension of this mesh will be one order less that the spacial dimension of the problem (i.e. a 2D boundary mesh for a 3D problem, 1D boundary mesh for a 2D problem etc. ). Once again let’s consider the simple mesh in Figure 1. Suppose we wish to apply a boundary condition on the right edge of the mesh then the boundary mesh would be the deﬁned by the following element connectivity matrix of 2-node line elements right Edge = 2 4 4 6 . (3) Note that the numbers in the boundary connectivity matrix refer to the same node coordinate matrix as do the numbers in the connectivity matrix of the interior elements. If we wish to apply an essential boundary conditions on this edge we need a list of the node numbers on the edge. This can be easily done in Matlab with the unique function. nodesOnBoundary = unique(rightEdge); This will set the vector nodesOnBoundary equal to [2 4 6]. If we wish to from a force vector from a natural boundary condition on this edge we simply loop over the elements and integrate the force on the edge just as we would integrate any ﬁnite element operators on the domain interior i.e. the stiﬀness matrix K. 5

4.4 Dof Mapping Ultimately for all ﬁnite element programs we solve a linear algebraic system of the form Kd = f (4) for the vector d. The vector d contains the nodal unknowns for that deﬁne the ﬁnite element approximation uh(x) = NI(x)dI (5) nn I=1 where NI(x) are the ﬁnite element shape functions, dI are the nodal un- knowns for the node I which may be scalar or vector quantities (if uh(x) is a scalar or vector) and nn is the number of nodes in the discretization. For scalar ﬁelds the location of the nodal unknowns in d is most obviously as follows dI = d(I), (6) but for vector ﬁelds the location of the nodal unknown dIi, where I refers to the node number and i refers to the component of the vector nodal unknown dI, there is some ambiguity. We need to deﬁne a mapping from the node number and vector component to the index of the nodal unknown vector d. This mapping can be written as f : {I, i} → n (7) where f is the mapping, I is the node number, i is the component and n is the index in d. So the location of unknown uIi in d is as follows uIi = df (I,i). (8) There are two common mappings used. The ﬁrst is to alternate between each spacial component in the nodal unknown vector d. With this arrange- ment the nodal unknown vector d is of the form d = (9)   u1x u1y ... u2x u2y ... unn x unn y ... 6

where nn is again the number of nodes in the discretization. This mapping is n = sdim(I − 1) + i. (10) With this mapping the i component of the displacement at node I is located as follows in d dIi = d( sdim*(I-1) + i ). (11) The other option is to group all the like components of the nodal un- knowns in a contiguous portion of d or as follows  u1x u2x ... unx u1y u2y ...  d = The mapping in this case is n = (i − 1)nn + I (12) (13) So for this structure the i component of the displacement at node I is located at in d dIi = d( (i-1)*nn + I ) (14) For reasons that will be appreciated when we discuss the scattering of element operators into system operators we will adopt the latter dof mapping. It is important to be comfortable with these mappings since it is an operation that is performed regularly in any ﬁnite element code. Of course which ever mapping is chosen the stiﬀness matrix and force vectors should have the same structure. 5 Computation of Finite Element Operators At the heart of the ﬁnite element program is the computation of ﬁnite element operators. For example in a linear static code they would be the stiﬀness matrix K = Ω BT C B dΩ (15) 7

and the external force vector f ext = Γt Nt dΓ. (16) The global operators are evaluated by looping over the elements in the dis- cretization, integrating the operator over the element and then to scatter the local element operator into the global operator. This procedure is written mathematically with the Assembly operator A K = Ae Ωe BeT C Be dΩ (17) 5.1 Quadrature The integration of an element operator is performed with an appropriate quadrature rule which depends on the element and the function being inte- grated. In general a quadrature rule is as follows ξ=1 ξ=−1 f (ξ)dξ = q f (ξq)Wq (18) where f (ξ) is the function to be integrated, ξq are the quadrature points and Wq the quadrature weights. The function quadrature generates a vector of quadrature points and a vector of quadrature weights for a quadrature rule. The syntax of this function is as follows [quadWeights,quadPoints] = quadrature(integrationOrder, elementType,dimensionOfQuadrature); so an example quadrature loop to integrate the function f = x3 on a trian- gular element would be as follows [qPt,qWt]=quadrature(3,’TRIANGULAR’,2); for q=1:length(qWt) xi = qPt(q); % quadrature point % get the global coordinte x at the quadrature point xi % and the Jacobian at the quadrature point, jac ... f_int = f_int + x^3 * jac*qWt(q); end 8

资料库

Programming The Finite Element Method With Matlab.pdf

相关推荐

课程资源

热门标签

最新资料