Differential equations  

Navier–Stokes differential equations used to simulate airflow around an obstruction.  
Classification  
Types


Relation to processes 

Solution  
General topics 

Solution methods 

The finite element method (FEM) is a numerical method for solving problems of engineering and mathematical physics. Typical problem areas of interest include structural analysis, heat transfer, fluid flow, mass transport, and electromagnetic potential. The analytical solution of these problems generally require the solution to boundary value problems for partial differential equations. The finite element method formulation of the problem results in a system of algebraic equations. The method approximates the unknown function over the domain.^{[1]} To solve the problem, it subdivides a large system into smaller, simpler parts that are called finite elements. The simple equations that model these finite elements are then assembled into a larger system of equations that models the entire problem. FEM then uses variational methods from the calculus of variations to approximate a solution by minimizing an associated error function.
Studying or analyzing a phenomenon with FEM is often referred to as finite element analysis (FEA).
YouTube Encyclopedic

1/5Views:325 65857 70927 3484 46651 458

✪ Mod01 Lec03 Introduction to Finite Element Method

✪ Finite Element Method  Differential Equations in Action

✪ Finite element analysis ( FEA) formulation  One dimensional heat transfer

✪ Isoparametric Elements in Finite Element Method

✪ 8.3.1PDEs: Introduction to Finite Element Method
Transcription
In the last class we saw a number of applications of finite element analysis. We saw really the breadth of applications. We saw a number of case studies from different fields and I am sure that you are all convinced that this is a very good technique which can be applied to lot of practical problems. From now on, we will study some of the basis or theoretical basis of finite element analysis. You may have to have background in various fields in order to understand exactly how this technique works. But as far as possible, in this course, we will see to it that all these things are covered and the course is completely self contained. We stopped with the definition of stiffness. We said that or we said that it is easy to understand what stiffness is as far as the spring is concerned, but how do you define stiffness for other components? We considered that there is force P that is acting on the spring and if the displacement at this point is u, then we said that it is possible to define a term called stiffness K such that P is equal to Ku. Now the question is how am I going to extend this concept to other components or other bodies which are not necessarily nice springs like this? As a first case let us consider a simple bar and let us say that I apply a force P. Can I talk something as stiffness for this particular bar? Is it possible that we can define it? Yes; how do we define it? May be all of you are familiar with it. We know what is the relationship between stress and strain? What is In a very elementary level you would have defined stress as P by A where A is the area of this particular bar. So, P by A is equal to E into change in length. What is change in length in this case, say for example? No; as far as the strain is concerned what is change in length by original length; that is how you would have defined. Let me call this change in length as u, in a similar fashion as we say that the displacement of the spring is u divided by L. It is very simple to now say that P is equal to AE divided by L into u; AE divided by L into u. Now compare this and this. It is quite clear that I can call AE by L as the stiffness of the bar. Stiffness straight away is not a concept which is restricted only to a spring but can be extended to say for example a bar. We are going to extend this for all other types of components as well. Now let us look at another problem. Let us see whether we can extend this kind of concept to a problem, where a simple case a very simple problem, where I have a P, force that is acting on a bar whose, for example, cross sections are say A1 and A2. This poses slightly more difficulty when compared to this problem or it is slightly more difficult. You can, looking at it immediately say that “Sir it is made up of two bars”? It is made up of two bars. So, I can say that why not you consider it as two bars? Once I start considering it as two bars you come into the realm of discretizing this particular piece; you have discretized it into two bars. Now on we take over or we come into the realm of finite elements. Yes, that is nice answer, but there is a small problem. What is the problem? Here again, like this case, we had stiffness to be only one number AE by L. Because one end was fixed it was nice to me to say that the displacement is considered only at the bottom of the bar and we arrived at this point. On other hand look at this; it is not like that. I have now two bars. If you look at these two bars which are joined at this point, this point starts moving, this point also starts moving. One end is not fixed; easy to handle, but slightly more complex than this. You can look at another problem where there are more number of bars. There is a technique to handle this. In order to use that technique first let us define what is called as a bar element? What is a bar element? Let us take an element. I now started calling this as an element. Let me take that element and put it out. Immediately you see that there are two guys sitting on either side which are called as nodes. So, let me call those nodes as say i and j; a very generic term say i and j. Now let me define or determine the stiffness matrix. Note this, we are now going to define stiffness matrix of this particular bar which are defined by these nodes i and j. What is also important is the area. Let me call this as A, area of this bar as well as the length, this length which I would call as L and the Young’s modulus which I will call as E. Let us start with what we know. In this bar, let me apply a displacement which I call as say ui and fix j. Now let me call the force that is required to cause this displacement as Fi. From my previous result I know what Fi is. What is Fi? AE by L into ui. When I fix this j, this j is going to take a reaction or in another words there is also a force that is acting at j. Let me call this as Fj. From elementary strength of materials I know that Fj is a reaction and that it is, what is that? Opposite to Fi and of equal magnitude. That’s right; so, Fj is equal to minus Fi; that’s the first step. Now let me again repeat this with j free or applying a displacement of uj and fixing i. What is that I am going to do? I am now going to release this ….. and say that let j move from the original position by an amount uj and let me fix this position. In this situation what happens to Fj? That is straight forward and that happens to be AE by L into uj and now what happens to my Fi? That is minus of Fj or in other words that is minus AE by L into uj. Is that clear? Now coming back to my problem here, this particular bar, this particular bar now has two nodes and is in the position where both the nodes would start displacing or ui and uj exist for this element. Since I have derived in a very simple fashion what is Fi and Fj when one or the other are fixed, it is possible for me to straight away get to a situation where both ui and uj exists. It is nothing but a superposition of these two results. So, when both ui and uj exist, I can say that Fi is equal to AE by L into ui minus uj. Is that clear? What happens to Fj? Fj is equal to AE by L into minus ui plus uj. Is that clear? This can be expressed in a matrix form. What do I mean by that? I can express this as AE by L into 1, minus 1, minus 1 and 1 multiplied by ui uj is equal to Fi Fj. Look at this particular equation quite closely. You can immediately recognize what is that you can recognize? You can recognize that there is a force on one side and there is a displacement on the other side. In other words this is a relationship between force and displacement and that relationship has been brought about by this particular matrix. Is that clear? Comparing that equation with my first equation here where it is P is equal to Ku for a simple spring where I called K to be stiffness and similarly there is a relationship between force and displacement, I can call that matrix which I derived now to be what is called as the stiffness matrix. That is the stiffness matrix. Let us understand what this stiffness matrix means. Can I extend the concept of stiffness, the physical behavior or the physical understanding of stiffness what I have for a spring, to the stiffness matrix which we derived? If I had asked you to define stiffness and not write down an equation, you would have said that stiffness is equal to force per unit displacement. There was no problem because there is only one place where force is applied and one place where the displacement happens to be, there. What about in this situation? How do you now define here or how do you physically get a picture of what the stiffness matrix is? That is the question. How do you define it? Can someone tell that or I would call this as a K11 K12 K21 and K22, along with AE by L. The same matrix can be written for us to interpret as K11 K12 K21 and K22 that multiplied by say u1 and u2 where i and j’s I have replaced by u1 and u2. If you want you can keep this as i and j itself. If you are getting confused you can keep it as i and j itself. My question is what does K11 indicate? What is the physical meaning of K11? What does K12 indicate and so on because here I said K indicates force per unit displacement. In order to understand that, in order to give you a clue let me say that I apply a unit displacement. Here it is force per unit displacement. Keeping that spirit let me say that ui is equal to 1 and let me fix the other side uj is equal to zero and then look at the complete equation now. What is that in the other side? Fi and Fj. What does the first equation indicate? First equation; K11 is equal to Fi and what does K21 indicate? K21 is what? Fj. Is that clear? Can someone give me an interpretation of what is K11 and what is K21? One of you; yes, so it is nothing but, yes, in other words just to repeat that what it means is K11 is the force that happens to be present at i or the force to be applied at i in order to have, note that, unit displacement at i and zero displacement at j and K21 is the force that is required to be applied at j in order to have a unit displacement at 1 and zero displacement at j. The same thing, same logic can be repeated for j as well; you can put this as zero and this as 1. You can again interpret this Kij. I can have a physical interpretation of Kij’s which are entrées here. i is equal to 1, j is equal to 1; it becomes K11. i is equal to 2, j is equal to 1, it becomes K21 and so on. So, I can have physical interpretation for Kij. What is the physical interpretation? What does Kij indicate? It is the force that has to be applied at j in order to keep ….. Now look at this and tell me. See, K21 means, what is the displacement I have applied for i? It is unit displacement. So, if I have to apply a unit displacement at 1, and a zero displacement at j then K21 is obtained from Fj. If there is a unit displacement at j, then Kij indicates the force that has to sit at i. It is the force that is required to make or to give a unit displacement at j and the force being measured at i, the force that is required at i for a unit displacement of j. So, that is a simple physical interpretation. Yes, zero displacement at i; force that is required to be present at i for a unit displacement at j, zero displacement at i. Now let us look at this matrix closely. What are the other things that you see in this matrix? What are the other things that you see in this matrix? Symmetric, very good; so, this matrix is a symmetric matrix. It is not necessary that every matrix that we are going to come across is symmetric matrix; but this particular matrix because we are just dealing with a linear elastic problem, this matrix happens to be a symmetric matrix. Very rarely do we come across matrices which are not symmetric and in this course most of the problems except when we go to contact we will only be dealing with symmetric matrices. Some of the large deformation problems have nonsymmetric matrix, let us not worry about it. But as you correctly said, the next thing that you observe is that this matrix is a symmetric matrix. It is very simple now to understand that I can have say two bars which are here or two sections in a bar to have different A’s, different E’s and different L’s. The only thing is that A, E and L would change. It will become A1, E1, L1; A2, E2, L2 and so on. Now let us see how we can determine a composite or an assembled matrix of a bar which looks like this. This particular bar which we considered as i and j have what are called as how many degrees of freedom? Two degrees of freedom or i and j whose displacements are ui and uj give the freedom for this bar to be displaced in the X direction. So, there are two degrees of freedom or if you consider each node there is one degree of freedom. What do I mean by degree of freedom? It means that this particular node i can be displaced only in the X direction and j can be displaced in the X direction. As an element this has two degrees of freedom, one degree of freedom per node. Now let us see how this concept can be applied to a composite bar and how assembly can be done so that I get a composite stiffness of this bar. For that the first thing I have to do is to determine the stiffness of individual bar. But before that I have already done one job. What is it that I have done? I have already discretized this particular component into two elements. I did it by looking at it and seeing that there is a change in A, change in E and so on, if it happens to be there, and so I have used that knowledge in order to discretize it. But there are lots of other considerations when I discretize it, but that we will see as we go on. But right now it is important to understand that discretization is very important and can depend upon a number of factors, one factor being a nice physical picture or an engineering sense and that is what has helped me to discretize it. Let me now call this as element 1 and let me call this as element 2 and let me right now worry about fixing this particular node. Let me call that node as 1, sorry 1, 2 and 3. Is that clear? Let me call A1, This is say L1; let me give some lengths to it. Let me say that it is L1 and this length to be say, L2. Just to generalize this problem, let me call the Young’s modulus to be E1 of this bar and this to be E2. Since AE by L comes up as a unit, let me define K1 to be A1E1 by L1 and K2 to be A2E2 by L2. 2640K1 is equal to A1E1 by L1 and K2 be A2E2 by L2. First what I can do is to write down the stiffnesses of individual elements and look at this kind of assembly for individual elements or how does this equation look like? How does it look like for the first element? I can change that here itself. How does it look like? It will be, for the first element, K1, minus K1; K2 does not come here, we are only looking at the first element minus K1 and again K1. Now what happens to i and j? Now what happens to i and j? u1 and u2. What happens to Fi and Fj and that happens to be F1 and F2. I am not worried about whatever the values of F1, how it looks like and so on; just I am repeating it for this particular element. Any questions? i and j was generic indices, now I have replaced it with the actual index. Let us now look at the second element. The second element can also be written in the same fashion. I will do it in the same matrix or maybe we will write it separately so that things become clear. What do I get? Instead of K1, K2; so, K2, minus K2, minus K2 and K2. Instead of now 1 and 2 what are our i’s and j’s? u2 and u3; what happens to this guy here? F2 and F3; F2 and F3. Now I have to assemble them. Please bear in mind what are the steps that we are doing? We discretized it first step. I will write it down at the end of this exercise but nevertheless as we go along please bear in mind what we are trying to do. We discretized it in terms of what we called as elements. In this case, the behavior is bar. We will come to what is behavior in a minute, but in this case we call this as bar element, so, we discretized it. The next step is we wrote down the stiffnesses. That is the next step; we wrote down the stiffnesses of individual element. The third step is I have to assemble them. How am I going to assemble or what does this term assemble mean? In order to assemble them let us see or let us imagine that we are building this bar. Let us not worry about the actual manufacturing process; it may be turning …… let us not worry about that. We will assemble them as if first I will put this and later put that. In fact, it may well turn out to be a welding process because E1 and E2 are different; like we may imagine that they are welded together. So, first I take this piece, the first element, put it there and identify these guys, these nodes there. Let us imagine that that is what happens. Since I have already identified a node at 3, it is imperative that it also has a displacement u3. Let us say that it has a displacement u3, though more importantly this element does not support any force there. If I apply a force here, the force happens to be in the air and this particular element does not support a force. In other words there is no stiffness. See, when can you apply a force? Only when there is a resistance to it or in other words there is stiffness to it; stiffness to the body on which you are applying. Since it is in vacuum, let us imagine that; let us not worry about anything else. So, I cannot put a force at that place or even if I put a force, there is no resistance due to this bar. If I have to now expand u1 as u1 u2 and u3, how will my stiffness matrix look like from our physical understanding of K? Let us see how it would like? So, I will have K1, no doubt; I will have minus K1. These two guys are for u1 and u2 and what will happen to my third term there? Very good, because you may imagine or remember our Kij definition; so, here I will have minus K1 and then followed by K1 0 0 0 0; that is very good. Let me go and put the other element here, weld them together or in other words weld them together. As soon as I do that, this force what we have mentioned here will be supported. Let us not worry about that right now, but the process is quite clear and hence some of the zeroes are now going to be replaced by the stiffness we are giving, by introducing the next element which means that I will have K2 coming to the aid of the first element; so, K1 plus K2, minus K2, minus K2. Please note that this is the last term in the matrix that is in the second row and 2, 3 and K2. One of the things you had not noticed in the last stiffness; when I said that you just said symmetric. I will look at it but you have not noticed one very important thing in this stiffness matrix which I will come to in a minute but what is important here is that we have assembled the matrix. First thing is that look at the matrix size. Now there are three nodes and there is a 3 by 3 stiffness matrix and there are u1 u2 u3 that are not known and what is the right hand side? What is the right hand side? Correct. It is F1 F2 and F3; remember that they are nothing but the forces that can be applied. It is not necessary that they should exist but that can be applied to 1, 2 and 3. So, this is the complete equation which I am interested in after assembling. What is my next step? I have to apply what is called as boundary condition. I have to apply what is called as boundary condition. But before we go to boundary condition let us look at this matrix closely and observe one more point. Again physically there is lot of mathematical implications to it. At this stage let us not worry about the mathematical implications but physically let us look at it and see whether we know why it happens to be in this fashion. Look at that matrix. You will see that there are some positive terms and negative terms. K1 is positive, K2 is positive; so, there are some positive terms and there are some negative terms. See that there is no negative term in the diagonals. All the diagonal terms happens to be positives. Apart from that you also see one more thing. What is that you see? Zeroes; so, there are positive terms, negative terms as well as some zeroes and it is very clear that the diagonal terms are positive, off diagonal terms, there are negative terms as well as zero terms. In this particular problem there are negative terms as well as there are zero terms. The question is, is it that these terms that are sitting in the diagonals are positive only for this problem it happens to be positive or is it by design or it is by chance? Is it that it is positive only for this problem or for any other problem? Is the question clear? What do you think will be the answer? One minute; let me give you a clue. In order to answer that question, let us follow the same logic by applying say unit displacement at say for example 1 and then zero displacements and so on. Now you can see or you can say whether this is positive by design and not by chance. Let us look at my first equation, say K1 say u1 ……. not a unit displacement at least u1 is equal to F1. So, u1 is equal to F1 by K1. u1 is equal to F1 by K1. What does this tell? If I apply a force F1, then u1 would follow the force in the same direction as long as K1 is positive. So, if I apply a force I know that the displacement has to be in the same direction; it cannot be in the opposite direction. So that particular condition is satisfied as long as K1 is positive. On the other hand if K1 happens to be negative, what happens? Say for example I take this bar and then I start applying a force here. It will move in this direction if K1 is negative. That is not possible. The bar has to move in the same direction. What does it indicate? It indicates that K1 has to be positive and the next thing is that, is it clear? That is very simple thing. On the other hand can K1 be zero; then u1 becomes infinity. So, necessarily these diagonal terms have to be positive. We will come to a much larger condition of positive definiteness of stiffness matrices later. Let us not worry about it because we have lot to learn from the physics of the problem before we go into the mathematics of the problem. I will not define all those things right now. Ultimately we will come to that. What is my next step? What is my next step, what do you think is the next step? Boundary condition, very good; so, what is the boundary condition of this problem? This problem states that one is fixed and that a force is applied at this point and that is the secondary thing; we will come to that later but first thing is that u1 is fixed. Please note one thing carefully that when I say boundary condition and say u1 is fixed, I just state that u1 is known. It does not mean many students make a mistake, that boundary condition means u1 is equal to zero; a particular degree of freedom, u1 or whatever be the degree of freedom that happens to be zero; not necessary. It just means that the displacement for that degree of freedom is stated. In this problem, it is stated to be zero. Suppose in another problem I say that this displacement is 0.5 mm; that is also a boundary condition. Boundary condition simply means in this particular problem displacement happens to be zero. What are the types of boundary conditions and what are the ways in which to express them, all these things we will see as we go along. But, for this particular problem it is important to understand that the displacement happens to be a boundary condition and that, that displacement happens to be zero. Why I am stating this point is because a same configuration, as far as geometry is concerned, the same configuration can be looked at as a beam as well. The only thing is if I now apply a force in this direction, it is no more a bar, a beam; it becomes a beam which means that it can bend. A bar can take only axial loads; you know from strength of materials. On the other hand a beam can take a bending load. Is that clear? Again, I have different types of displacements, a displacement, which looks something like this and I have different types of boundary condition. What are the boundary conditions? Look at the question, what are the boundary conditions? The previous case for bar we defined what is the boundary condition? There was only one; now we just modify the problem and say that what are the boundary conditions? Displacement and rotation that is right. So, displacement as well as rotation, both is specified at that place. In other words the same node or the same position, same point is specified in a slightly different fashion; the geometry happens to be the same. But, why is this difference? That is because of the governing differential equation for these two problems; they happen to be different. Before we go further, before we go to that point let us now worry only about this problem and get a complete physical picture of finite element analysis. So, we will continue with this. As you correctly said u1 is specified and u1 happens to be zero. So, how many unknowns are there? Now two unknowns; u2 and u3 and u1 is already known and happens to be zero. Since u1 is equal to zero, I eliminate these two lines or in other words I eliminate that row and the column and restrict my attention to this particular piece. Now I know something else also in this equation. What is it that I know? What is it that I know in this equation? F3; correct. So, I know that F3 happens to be the applied force P and what happens to F2? Do I know it or I do not know it? Is it that we do not know? Note that carefully; this is the mistake many students make. No, no; please note very carefully. Look at the problem again, look at this problem again. Look at this problem, look at this. Now, where is that node? Node is sitting here. It is nice that you have made that answer because this is usually the confusion. Please note that what we are doing here is to put down the external forces that are applied. We are not worried about the corresponding internal forces that are developed; we are looking at the external forces. These are the external forces. Please do not get confused between that and the internal forces which are developed due to stresses. This is the usual confusion for students and hence this is an external force, P is an external force. Correspondingly is there an external force at node 2? No; so, this happens to be zero. Now I have two equations and two unknowns. Just to summarize, what did we do? We discretized, we divide or in other words we divided this into a number of elements. Then we wrote the stiffness matrix for each of these elements. We assembled them, that is what we did here and then we wrote down boundary conditions carefully and also put down the forces that happen. One other thing before we go ahead and solve this problem; it is very important to note from this equation. What is it? What is that you see; another thing that you see? You please note that at the place where we give a boundary condition we are not specifying the force. It is not possible to specify both a boundary condition as well as a force. What does it mean? What is it mean? It means that when I fix a particular node, the force that acts on the node is decided by the system. So, if I say P, it is very clear from this problem before even solving it that the reaction for P has to be taken at the place where I have fixed. So, the force that is going to happen or going to be present at this place is going to be decided by the force that you give here. I cannot specify both a force as well as a boundary condition at the same point. When I specify a force at this point, the displacement at this point is decided by this system, by the stiffness matrices and when I specify a displacement at this point the force that is required to sustain that effect is specified again by the system, by equilibrium and so on. We will have more interpretation of this equation and how to solve them or how to solve this equation in the next class.
Contents
 1 Basic concepts
 2 History
 3 Technical discussion
 4 Discretization
 5 Various types of finite element methods
 5.1 AEM
 5.2 Generalized finite element method
 5.3 Mixed finite element method
 5.4 hpFEM
 5.5 hpkFEM
 5.6 XFEM
 5.7 Scaled boundary finite element method (SBFEM)
 5.8 SFEM
 5.9 Spectral element method
 5.10 Meshfree methods
 5.11 Discontinuous Galerkin methods
 5.12 Finite element limit analysis
 5.13 Stretched grid method
 5.14 Loubignac iteration
 6 Link with the gradient discretisation method
 7 Comparison to the finite difference method
 8 Application
 9 See also
 10 References
 11 Further reading
 12 External links
Basic concepts
The subdivision of a whole domain into simpler parts has several advantages:^{[2]}
 Accurate representation of complex geometry
 Inclusion of dissimilar material properties
 Easy representation of the total solution
 Capture of local effects.
A typical work out of the method involves (1) dividing the domain of the problem into a collection of subdomains, with each subdomain represented by a set of element equations to the original problem, followed by (2) systematically recombining all sets of element equations into a global system of equations for the final calculation. The global system of equations has known solution techniques, and can be calculated from the initial values of the original problem to obtain a numerical answer.
In the first step above, the element equations are simple equations that locally approximate the original complex equations to be studied, where the original equations are often partial differential equations (PDE). To explain the approximation in this process, FEM is commonly introduced as a special case of Galerkin method. The process, in mathematical language, is to construct an integral of the inner product of the residual and the weight functions and set the integral to zero. In simple terms, it is a procedure that minimizes the error of approximation by fitting trial functions into the PDE. The residual is the error caused by the trial functions, and the weight functions are polynomial approximation functions that project the residual. The process eliminates all the spatial derivatives from the PDE, thus approximating the PDE locally with
 a set of algebraic equations for steady state problems,
 a set of ordinary differential equations for transient problems.
These equation sets are the element equations. They are linear if the underlying PDE is linear, and vice versa. Algebraic equation sets that arise in the steady state problems are solved using numerical linear algebra methods, while ordinary differential equation sets that arise in the transient problems are solved by numerical integration using standard techniques such as Euler's method or the RungeKutta method.
In step (2) above, a global system of equations is generated from the element equations through a transformation of coordinates from the subdomains' local nodes to the domain's global nodes. This spatial transformation includes appropriate orientation adjustments as applied in relation to the reference coordinate system. The process is often carried out by FEM software using coordinate data generated from the subdomains.
FEM is best understood from its practical application, known as finite element analysis (FEA). FEA as applied in engineering is a computational tool for performing engineering analysis. It includes the use of mesh generation techniques for dividing a complex problem into small elements, as well as the use of software program coded with FEM algorithm. In applying FEA, the complex problem is usually a physical system with the underlying physics such as the EulerBernoulli beam equation, the heat equation, or the NavierStokes equations expressed in either PDE or integral equations, while the divided small elements of the complex problem represent different areas in the physical system.
FEA is a good choice for analyzing problems over complicated domains (like cars and oil pipelines), when the domain changes (as during a solid state reaction with a moving boundary), when the desired precision varies over the entire domain, or when the solution lacks smoothness. FEA simulations provide a valuable resource as they remove multiple instances of creation and testing of hard prototypes for various high fidelity situations.^{[3]} For instance, in a frontal crash simulation it is possible to increase prediction accuracy in "important" areas like the front of the car and reduce it in its rear (thus reducing cost of the simulation). Another example would be in numerical weather prediction, where it is more important to have accurate predictions over developing highly nonlinear phenomena (such as tropical cyclones in the atmosphere, or eddies in the ocean) rather than relatively calm areas.
History
While it is difficult to quote a date of the invention of the finite element method, the method originated from the need to solve complex elasticity and structural analysis problems in civil and aeronautical engineering. Its development can be traced back to the work by A. Hrennikoff^{[4]} and R. Courant^{[5]} in the early 1940s. Another pioneer was Ioannis Argyris. In the USSR, the introduction of the practical application of the method is usually connected with name of Leonard Oganesyan.^{[6]} In China, in the later 1950s and early 1960s, based on the computations of dam constructions, K. Feng proposed a systematic numerical method for solving partial differential equations. The method was called the finite difference method based on variation principle, which was another independent invention of the finite element method. Although the approaches used by these pioneers are different, they share one essential characteristic: mesh discretization of a continuous domain into a set of discrete subdomains, usually called elements.
Hrennikoff's work discretizes the domain by using a lattice analogy, while Courant's approach divides the domain into finite triangular subregions to solve second order elliptic partial differential equations (PDEs) that arise from the problem of torsion of a cylinder. Courant's contribution was evolutionary, drawing on a large body of earlier results for PDEs developed by Rayleigh, Ritz, and Galerkin.
The finite element method obtained its real impetus in the 1960s and 1970s by the developments of J. H. Argyris with coworkers at the University of Stuttgart, R. W. Clough with coworkers at UC Berkeley, O. C. Zienkiewicz with coworkers Ernest Hinton, Bruce Irons^{[7]} and others at the Swansea University, Philippe G. Ciarlet at the University of Paris 6 and Richard Gallagher with coworkers at Cornell University. Further impetus was provided in these years by available open source finite element software programs. NASA sponsored the original version of NASTRAN, and UC Berkeley made the finite element program SAP IV^{[8]} widely available. In Norway the ship classification society Det Norske Veritas (now DNV GL) developed Sesam in 1969 for use in analysis of ships.^{[9]} A rigorous mathematical basis to the finite element method was provided in 1973 with the publication by Strang and Fix.^{[10]} The method has since been generalized for the numerical modeling of physical systems in a wide variety of engineering disciplines, e.g., electromagnetism, heat transfer, and fluid dynamics.^{[11]}^{[12]}
Technical discussion
The structure of finite element methods
A finite element method is characterized by a variational formulation, a discretization strategy, one or more solution algorithms and postprocessing procedures.
Examples of variational formulation are the Galerkin method, the discontinuous Galerkin method, mixed methods, etc.
A discretization strategy is understood to mean a clearly defined set of procedures that cover (a) the creation of finite element meshes, (b) the definition of basis function on reference elements (also called shape functions) and (c) the mapping of reference elements onto the elements of the mesh. Examples of discretization strategies are the hversion, pversion, hpversion, xFEM, isogeometric analysis, etc. Each discretization strategy has certain advantages and disadvantages. A reasonable criterion in selecting a discretization strategy is to realize nearly optimal performance for the broadest set of mathematical models in a particular model class.
There are various numerical solution algorithms that can be classified into two broad categories; direct and iterative solvers. These algorithms are designed to exploit the sparsity of matrices that depend on the choices of variational formulation and discretization strategy.
Postprocessing procedures are designed for the extraction of the data of interest from a finite element solution. In order to meet the requirements of solution verification, postprocessors need to provide for a posteriori error estimation in terms of the quantities of interest. When the errors of approximation are larger than what is considered acceptable then the discretization has to be changed either by an automated adaptive process or by action of the analyst. There are some very efficient postprocessors that provide for the realization of superconvergence.
Illustrative problems P1 and P2
We will demonstrate the finite element method using two sample problems from which the general method can be extrapolated. It is assumed that the reader is familiar with calculus and linear algebra.
P1 is a onedimensional problem
where is given, is an unknown function of , and is the second derivative of with respect to .
P2 is a twodimensional problem (Dirichlet problem)
where is a connected open region in the plane whose boundary is nice (e.g., a smooth manifold or a polygon), and and denote the second derivatives with respect to and , respectively.
The problem P1 can be solved directly by computing antiderivatives. However, this method of solving the boundary value problem (BVP) works only when there is one spatial dimension and does not generalize to higherdimensional problems or to problems like . For this reason, we will develop the finite element method for P1 and outline its generalization to P2.
Our explanation will proceed in two steps, which mirror two essential steps one must take to solve a boundary value problem (BVP) using the FEM.
 In the first step, one rephrases the original BVP in its weak form. Little to no computation is usually required for this step. The transformation is done by hand on paper.
 The second step is the discretization, where the weak form is discretized in a finitedimensional space.
After this second step, we have concrete formulae for a large but finitedimensional linear problem whose solution will approximately solve the original BVP. This finitedimensional problem is then implemented on a computer.^{[13]}
Weak formulation
The first step is to convert P1 and P2 into their equivalent weak formulations.
The weak form of P1
If solves P1, then for any smooth function that satisfies the displacement boundary conditions, i.e. at and , we have
(1)
Conversely, if with satisfies (1) for every smooth function then one may show that this will solve P1. The proof is easier for twice continuously differentiable (mean value theorem), but may be proved in a distributional sense as well.
We define a new operator or map by using integration by parts on the righthandside of (1):
(2)
where we have used the assumption that .
The weak form of P2
If we integrate by parts using a form of Green's identities, we see that if solves P2, then we may define for any by
where denotes the gradient and denotes the dot product in the twodimensional plane. Once more can be turned into an inner product on a suitable space of once differentiable functions of that are zero on . We have also assumed that (see Sobolev spaces). Existence and uniqueness of the solution can also be shown.
A proof outline of existence and uniqueness of the solution
We can loosely think of to be the absolutely continuous functions of that are at and (see Sobolev spaces). Such functions are (weakly) once differentiable and it turns out that the symmetric bilinear map then defines an inner product which turns into a Hilbert space (a detailed proof is nontrivial). On the other hand, the lefthandside is also an inner product, this time on the Lp space . An application of the Riesz representation theorem for Hilbert spaces shows that there is a unique solving (2) and therefore P1. This solution is apriori only a member of , but using elliptic regularity, will be smooth if is.
Discretization
P1 and P2 are ready to be discretized which leads to a common subproblem (3). The basic idea is to replace the infinitedimensional linear problem:
 Find such that
with a finitedimensional version:
 (3) Find such that
where is a finitedimensional subspace of . There are many possible choices for (one possibility leads to the spectral method). However, for the finite element method we take to be a space of piecewise polynomial functions.
For problem P1
We take the interval , choose values of with and we define by:
where we define and . Observe that functions in are not differentiable according to the elementary definition of calculus. Indeed, if then the derivative is typically not defined at any , . However, the derivative exists at every other value of and one can use this derivative for the purpose of integration by parts.
For problem P2
We need to be a set of functions of . In the figure on the right, we have illustrated a triangulation of a 15 sided polygonal region in the plane (below), and a piecewise linear function (above, in color) of this polygon which is linear on each triangle of the triangulation; the space would consist of functions that are linear on each triangle of the chosen triangulation.
One hopes that as the underlying triangular mesh becomes finer and finer, the solution of the discrete problem (3) will in some sense converge to the solution of the original boundary value problem P2. To measure this mesh fineness, the triangulation is indexed by a real valued parameter which one takes to be very small. This parameter will be related to the size of the largest or average triangle in the triangulation. As we refine the triangulation, the space of piecewise linear functions must also change with . For this reason, one often reads instead of in the literature. Since we do not perform such an analysis, we will not use this notation.
Choosing a basis
To complete the discretization, we must select a basis of . In the onedimensional case, for each control point we will choose the piecewise linear function in whose value is at and zero at every , i.e.,
for ; this basis is a shifted and scaled tent function. For the twodimensional case, we choose again one basis function per vertex of the triangulation of the planar region . The function is the unique function of whose value is at and zero at every .
Depending on the author, the word "element" in "finite element method" refers either to the triangles in the domain, the piecewise linear basis function, or both. So for instance, an author interested in curved domains might replace the triangles with curved primitives, and so might describe the elements as being curvilinear. On the other hand, some authors replace "piecewise linear" by "piecewise quadratic" or even "piecewise polynomial". The author might then say "higher order element" instead of "higher degree polynomial". Finite element method is not restricted to triangles (or tetrahedra in 3d, or higher order simplexes in multidimensional spaces), but can be defined on quadrilateral subdomains (hexahedra, prisms, or pyramids in 3d, and so on). Higher order shapes (curvilinear elements) can be defined with polynomial and even nonpolynomial shapes (e.g. ellipse or circle).
Examples of methods that use higher degree piecewise polynomial basis functions are the hpFEM and spectral FEM.
More advanced implementations (adaptive finite element methods) utilize a method to assess the quality of the results (based on error estimation theory) and modify the mesh during the solution aiming to achieve approximate solution within some bounds from the exact solution of the continuum problem. Mesh adaptivity may utilize various techniques, the most popular are:
 moving nodes (radaptivity)
 refining (and unrefining) elements (hadaptivity)
 changing order of base functions (padaptivity)
 combinations of the above (hpadaptivity).
Small support of the basis
The primary advantage of this choice of basis is that the inner products
and
will be zero for almost all . (The matrix containing in the location is known as the Gramian matrix.) In the one dimensional case, the support of is the interval . Hence, the integrands of and are identically zero whenever .
Similarly, in the planar case, if and do not share an edge of the triangulation, then the integrals
and
are both zero.
Matrix form of the problem
If we write and then problem (3), taking for , becomes
 for (4)
If we denote by and the column vectors and , and if we let
and
be matrices whose entries are
and
then we may rephrase (4) as
 (5)
It is not necessary to assume . For a general function , problem (3) with for becomes actually simpler, since no matrix is used,
 , (6)
where and for .
As we have discussed before, most of the entries of and are zero because the basis functions have small support. So we now have to solve a linear system in the unknown where most of the entries of the matrix , which we need to invert, are zero.
Such matrices are known as sparse matrices, and there are efficient solvers for such problems (much more efficient than actually inverting the matrix.) In addition, is symmetric and positive definite, so a technique such as the conjugate gradient method is favored. For problems that are not too large, sparse LU decompositions and Cholesky decompositions still work well. For instance, MATLAB's backslash operator (which uses sparse LU, sparse Cholesky, and other factorization methods) can be sufficient for meshes with a hundred thousand vertices.
The matrix is usually referred to as the stiffness matrix, while the matrix is dubbed the mass matrix.
General form of the finite element method
In general, the finite element method is characterized by the following process.
 One chooses a grid for . In the preceding treatment, the grid consisted of triangles, but one can also use squares or curvilinear polygons.
 Then, one chooses basis functions. In our discussion, we used piecewise linear basis functions, but it is also common to use piecewise polynomial basis functions.
A separate consideration is the smoothness of the basis functions. For second order elliptic boundary value problems, piecewise polynomial basis function that are merely continuous suffice (i.e., the derivatives are discontinuous.) For higher order partial differential equations, one must use smoother basis functions. For instance, for a fourth order problem such as , one may use piecewise quadratic basis functions that are .
Another consideration is the relation of the finitedimensional space to its infinitedimensional counterpart, in the examples above . A conforming element method is one in which the space is a subspace of the element space for the continuous problem. The example above is such a method. If this condition is not satisfied, we obtain a nonconforming element method, an example of which is the space of piecewise linear functions over the mesh which are continuous at each edge midpoint. Since these functions are in general discontinuous along the edges, this finitedimensional space is not a subspace of the original .
Typically, one has an algorithm for taking a given mesh and subdividing it. If the main method for increasing precision is to subdivide the mesh, one has an hmethod (h is customarily the diameter of the largest element in the mesh.) In this manner, if one shows that the error with a grid is bounded above by , for some and , then one has an order p method. Under certain hypotheses (for instance, if the domain is convex), a piecewise polynomial of order method will have an error of order .
If instead of making h smaller, one increases the degree of the polynomials used in the basis function, one has a pmethod. If one combines these two refinement types, one obtains an hpmethod (hpFEM). In the hpFEM, the polynomial degrees can vary from element to element. High order methods with large uniform p are called spectral finite element methods (SFEM). These are not to be confused with spectral methods.
For vector partial differential equations, the basis functions may take values in .
Various types of finite element methods
AEM
The Applied Element Method, or AEM combines features of both FEM and Discrete element method, or (DEM).
Generalized finite element method
The generalized finite element method (GFEM) uses local spaces consisting of functions, not necessarily polynomials, that reflect the available information on the unknown solution and thus ensure good local approximation. Then a partition of unity is used to “bond” these spaces together to form the approximating subspace. The effectiveness of GFEM has been shown when applied to problems with domains having complicated boundaries, problems with microscales, and problems with boundary layers.^{[14]}
Mixed finite element method
The mixed finite element method is a type of finite element method in which extra independent variables are introduced as nodal variables during the discretization of a partial differential equation problem.
hpFEM
The hpFEM combines adaptively, elements with variable size h and polynomial degree p in order to achieve exceptionally fast, exponential convergence rates.^{[15]}
hpkFEM
The hpkFEM combines adaptively, elements with variable size h, polynomial degree of the local approximations p and global differentiability of the local approximations (k1) in order to achieve best convergence rates.
XFEM
The extended finite element method (XFEM) is a numerical technique based on the generalized finite element method (GFEM) and the partition of unity method (PUM). It extends the classical finite element method by enriching the solution space for solutions to differential equations with discontinuous functions. Extended finite element methods enrich the approximation space so that it is able to naturally reproduce the challenging feature associated with the problem of interest: the discontinuity, singularity, boundary layer, etc. It was shown that for some problems, such an embedding of the problem's feature into the approximation space can significantly improve convergence rates and accuracy. Moreover, treating problems with discontinuities with XFEMs suppresses the need to mesh and remesh the discontinuity surfaces, thus alleviating the computational costs and projection errors associated with conventional finite element methods, at the cost of restricting the discontinuities to mesh edges.
Several research codes implement this technique to various degrees: 1. GetFEM++ 2. xfem++ 3. openxfem++
XFEM has also been implemented in codes like Altair Radioss, ASTER, Morfeo and Abaqus. It is increasingly being adopted by other commercial finite element software, with a few plugins and actual core implementations available (ANSYS, SAMCEF, OOFELIE, etc.).
Scaled boundary finite element method (SBFEM)
The introduction of the scaled boundary finite element method (SBFEM) came from Song and Wolf (1997).^{[16]} The SBFEM has been one of the most profitable contributions in the area of numerical analysis of fracture mechanics problems. It is a semianalytical fundamentalsolutionless method which combines the advantages of both the finite element formulations and procedures, and the boundary element discretization. However, unlike the boundary element method, no fundamental differential solution is required.
SFEM
The SFEM, Smoothed Finite Element Methods, are a particular class of numerical simulation algorithms for the simulation of physical phenomena. It was developed by combining meshfree methods with the finite element method.
Spectral element method
Spectral element methods combine the geometric flexibility of finite elements and the acute accuracy of spectral methods. Spectral methods are the approximate solution of weak form partial equations that are based on highorder Lagragian interpolants and used only with certain quadrature rules.^{[17]}
Meshfree methods
Discontinuous Galerkin methods
Finite element limit analysis
Stretched grid method
Loubignac iteration
Loubignac iteration is an iterative method in finite element methods.
Link with the gradient discretisation method
Some types of finite element methods (conforming, nonconforming, mixed finite element methods) are particular cases of the gradient discretisation method (GDM). Hence the convergence properties of the GDM, which are established for a series of problems (linear and non linear elliptic problems, linear, nonlinear and degenerate parabolic problems), hold as well for these particular finite element methods.
Comparison to the finite difference method
The finite difference method (FDM) is an alternative way of approximating solutions of PDEs. The differences between FEM and FDM are:
 The most attractive feature of the FEM is its ability to handle complicated geometries (and boundaries) with relative ease. While FDM in its basic form is restricted to handle rectangular shapes and simple alterations thereof, the handling of geometries in FEM is theoretically straightforward.
 FDM is not usually used for irregular CAD geometries but more often rectangular or block shaped models.^{[18]}
 The most attractive feature of finite differences is that it is very easy to implement.
 There are several ways one could consider the FDM a special case of the FEM approach. E.g., first order FEM is identical to FDM for Poisson's equation, if the problem is discretized by a regular rectangular mesh with each rectangle divided into two triangles.
 There are reasons to consider the mathematical foundation of the finite element approximation more sound, for instance, because the quality of the approximation between grid points is poor in FDM.
 The quality of a FEM approximation is often higher than in the corresponding FDM approach, but this is extremely problemdependent and several examples to the contrary can be provided.
Generally, FEM is the method of choice in all types of analysis in structural mechanics (i.e. solving for deformation and stresses in solid bodies or dynamics of structures) while computational fluid dynamics (CFD) tends to use FDM or other methods like finite volume method (FVM). CFD problems usually require discretization of the problem into a large number of cells/gridpoints (millions and more), therefore cost of the solution favors simpler, lower order approximation within each cell. This is especially true for 'external flow' problems, like air flow around the car or airplane, or weather simulation.
Application
A variety of specializations under the umbrella of the mechanical engineering discipline (such as aeronautical, biomechanical, and automotive industries) commonly use integrated FEM in design and development of their products. Several modern FEM packages include specific components such as thermal, electromagnetic, fluid, and structural working environments. In a structural simulation, FEM helps tremendously in producing stiffness and strength visualizations and also in minimizing weight, materials, and costs.^{[19]}
FEM allows detailed visualization of where structures bend or twist, and indicates the distribution of stresses and displacements. FEM software provides a wide range of simulation options for controlling the complexity of both modeling and analysis of a system. Similarly, the desired level of accuracy required and associated computational time requirements can be managed simultaneously to address most engineering applications. FEM allows entire designs to be constructed, refined, and optimized before the design is manufactured. The mesh is an integral part of the model and it must be controlled carefully to give the best results. Generally the higher the number of elements in a mesh, the more accurate the solution of the discretised problem. However, there is a value at which the results converge and further mesh refinement does not increase accuracy.^{[20]}
This powerful design tool has significantly improved both the standard of engineering designs and the methodology of the design process in many industrial applications.^{[22]} The introduction of FEM has substantially decreased the time to take products from concept to the production line.^{[22]} It is primarily through improved initial prototype designs using FEM that testing and development have been accelerated.^{[23]} In summary, benefits of FEM include increased accuracy, enhanced design and better insight into critical design parameters, virtual prototyping, fewer hardware prototypes, a faster and less expensive design cycle, increased productivity, and increased revenue.^{[22]}
In the 1990s FEA was proposed for use in stochastic modelling for numerically solving probability models^{[24]} and later for reliability assessment.^{[25]} The stochastic finite element method has since been applied to many branches of engineering,^{[26]} often being applied to characterise variability in material properties.^{[27]}
See also
 Applied element method
 Boundary element method
 Céa's lemma
 Computer experiment
 Direct stiffness method
 Discontinuity layout optimization
 Discrete element method
 Finite difference method
 Finite element machine
 Finite element method in structural mechanics
 Finite volume method
 Finite volume method for unsteady flow
 Infinite element method
 Interval finite element
 Isogeometric analysis
 Lattice Boltzmann methods
 List of finite element software packages
 Meshfree methods
 Movable cellular automaton
 Multidisciplinary design optimization
 Multiphysics
 Patch test
 Rayleigh–Ritz method
 Space mapping
 Tessellation (computer graphics)
 Weakened weak form
References
 ^ Daryl L. Logan (2011). A first course in the finite element method. Cengage Learning. ISBN 9780495668251.
 ^ Reddy, J. N. (2006). An Introduction to the Finite Element Method (Third ed.). McGrawHill. ISBN 9780071267618.
 ^ "Finite Elements Analysis (FEA)". www.manortool.com. Retrieved 20170728.
 ^ Hrennikoff, Alexander (1941). "Solution of problems of elasticity by the framework method". Journal of Applied Mechanics. 8.4: 169–175.
 ^ Courant, R. (1943). "Variational methods for the solution of problems of equilibrium and vibrations". Bulletin of the American Mathematical Society. 49: 1–23. doi:10.1090/s000299041943078184.
 ^ "СПб ЭМИ РАН". emi.nw.ru. Retrieved 17 March 2018.
 ^ Hinton, Ernest; Irons, Bruce (July 1968). "Least squares smoothing of experimental data using finite elements". Strain. 4 (3): 24–27. doi:10.1111/j.14751305.1968.tb01368.x.
 ^ "SAPIV Software and Manuals". NISEE eLibrary, The Earthquake Engineering Online Archive.
 ^ Gard Paulsen; Håkon With Andersen; John Petter Collett; Iver Tangen Stensrud (2014). Building Trust, The history of DNV 18642014. Lysaker, Norway: Dinamo Forlag A/S. pp. 121, 436. ISBN 9788280712561.
 ^ Strang, Gilbert; Fix, George (1973). An Analysis of The Finite Element Method. Prentice Hall. ISBN 9780130329462.
 ^ Olek C Zienkiewicz; Robert L Taylor; J.Z. Zhu (31 August 2013). The Finite Element Method: Its Basis and Fundamentals. ButterworthHeinemann. ISBN 9780080951355.
 ^ Bathe, K.J. (2006). Finite Element Procedures. Cambridge, MA: KlausJürgen Bathe. ISBN 9780979004902.
 ^ Smith, I.M.; Griffiths, D.V.; Margetts, L. (2014). Programming the Finite Element Method (Fifth ed.). Wiley. ISBN 9781119973348.
 ^ Babuška, Ivo; Banerjee, Uday; Osborn, John E. (June 2004). "Generalized Finite Element Methods: Main Ideas, Results, and Perspective". International Journal of Computational Methods. 1 (1): 67–103. doi:10.1142/S0219876204000083.
 ^ P. Solin, K. Segeth, I. Dolezel: HigherOrder Finite Element Methods, Chapman & Hall/CRC Press, 2003
 ^ Song, Chongmin; Wolf, John P. (5 August 1997). "The scaled boundary finiteelement method  alias consistent infinitesimal finiteelement cell method  for elastodynamics". Computer Methods in Applied Mechanics and Engineering. 147 (3–4): 329–355. Bibcode:1997CMAME.147..329S. doi:10.1016/S00457825(97)000212.
 ^ "Spectral Element Methods". State Key Laboratory of Scientific and Engineering Computing. Retrieved 20170728.
 ^ "What's The Difference Between FEM, FDM, and FVM?". Machine Design. 20160418. Retrieved 20170728.
 ^ Kiritsis, D.; Eemmanouilidis, Ch.; Koronios, A.; Mathew, J. (2009). "Engineering Asset Management". Proceedings of the 4th World Congress on Engineering Asset Management (WCEAM): 591–592.
 ^ "Finite Element Analysis: How to create a great model". Coventive Composites. 20190318. Retrieved 20190405.
 ^ Naghibi Beidokhti, Hamid; Janssen, Dennis; Khoshgoftar, Mehdi; Sprengers, Andre; Perdahcioglu, Emin Semih; Boogaard, Ton Van den; Verdonschot, Nico (2016). "A comparison between dynamic implicit and explicit finite element simulations of the native knee joint" (PDF). Medical Engineering & Physics. 38 (10): 1123–1130. doi:10.1016/j.medengphy.2016.06.001. PMID 27349493.
 ^ ^{a} ^{b} ^{c} Hastings, J. K., Juds, M. A., Brauer, J. R., Accuracy and Economy of Finite Element Magnetic Analysis, 33rd Annual National Relay Conference, April 1985.
 ^ McLarenMercedes (2006). "McLaren Mercedes: Feature  Stress to impress". Archived from the original on 20061030. Retrieved 20061003.
 ^ Peng Long; Wang Jinliang; Zhu Qiding (19 May 1995). "Methods with high accuracy for finite element probability computing". Journal of Computational and Applied Mathematics. 59 (2): 181–189. doi:10.1016/03770427(94)00027X.
 ^ Haldar, Achintya; Mahadevan, Sankaran (2000). Reliability Assessment Using Stochastic Finite Element Analysis. John Wiley & Sons. ISBN 9780471369615.
 ^ Arregui Mena, J.D.; Margetts, L.; et al. (2014). "Practical Application of the Stochastic Finite Element Method". Archives of Computational Methods in Engineering. 23 (1): 171–190. doi:10.1007/s1183101491393.
 ^ Arregui Mena, J.D.; et al. (2018). "Characterisation of the spatial variability of material properties of Gilsocarbon and NBG18 using random fields". Journal of Nuclear Materials. 511: 91–108. doi:10.1016/j.jnucmat.2018.09.008.
Further reading
 G. Allaire and A. Craig : Numerical Analysis and Optimization: An Introduction to Mathematical Modelling and Numerical Simulation
 K. J. Bathe : Numerical methods in finite element analysis, PrenticeHall (1976).
 Thomas J.R. Hughes : The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, PrenticeHall (1987).
 J. Chaskalovic : Finite Elements Methods for Engineering Sciences, Springer Verlag, (2008).
 Endre Süli  Finite Element Methods for Partial Differential Equations
 O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu : The Finite Element Method: Its Basis and Fundamentals, ButterworthHeinemann, (2005).
External links
Wikimedia Commons has media related to Finite element modelling. 
 IFER – Internet Finite Element Resources – describes and provides access to finite element analysis software via the Internet
 NAFEMS – International Association Engineering Modelling
 Mathematics of the Finite Element Method