Module 4 Lecture 2 Finite Element Method

We are at the lecture 12; what we had discussed till now was how to do the element calculations in the master element. We had transformed our equations from the physical domain to the master domain and in the master domain we had said that we will be defining the shape functions for the element and transform the integrals that we had over this domain. The question was how to convert these integrals into computer implementable form, because for us to do the integration of polynomials and other functions by hand is easy, but how will the computer do that? Computer generally what does it do? It integrates some quantity by replacing it with a sum. A sum over some points 1 to n bar let us say; so sum of the function psi, the f at the point psii into the so-called weights wi. These points psii are called integration points and wi is the corresponding weight. The question is fine, how do we choose points psii and the corresponding weights wi? What we use here is the so called Gauss Legendre Quadrature integration rule; what is the basic idea? I have now my master element spanning from psi equal to minus 1 to psi equal to plus 1. The centre of it is the point psi equal to 0. What I am going to do is I am going to choose points which are symmetrically placed on either side of 0, so I can choose this point and the corresponding point here . These are going to be our choices of integration points such that if I call this point as psi1 this point as psi2, this is psi3, psi4 then psi1 is equal to minus psi4, psi2 is equal to minus psi3. Further, for this set of symmetric located points we will have the same weights. So wi will be equal to w4 in the example that we have taken and w2 will be equal to w3. We will say that this, the feature of this method is (1) the integration points are symmetrically located on either side of zero. (2) The weights for a pair of such symmetric points is the same. How do you design the location and the weights of the points? The basic idea again is that, we would like our choice of points and weights to be such that for a polynomial of a particular order the integral comes out to be exact. If my integrant here is a polynomial, if this is a polynomial then for a particular order of the polynomial this integral should be exact; that is the sum should be exact, therefore psi. Let us now go ahead and construct these sets of points, take the first example. First example is let us take n is equal to 1. If I take n is equal to 1 then I notice that the point psi1 can be nowhere else, but at the 0 location because it is a symmetrically placed point and the corresponding weight will be w1. Now w1 is what we have to determine. Let us take first a constant; I integrate the constant from minus 1 to 1; let us take the constant as 1. This integral will be equal to 2 . Now this I would expect one point rule to at least integrate the constants correctly. It is going to give me from the one point rule F at 0 into w1 which is equal to 1 into w1; implies w1 is equal to 1. This is going to be our one point rule. So one point rule for n equal to 1 implies the following pair at the point is 0 and the weight is 2. For which higher polynomials does this one point rule do an exact job? Let us look at the next higher polynomial. If I take fpsi to be linear in psi, it is a0 plus a1 psi d psi. If I do this integral, I am going to call this integral as I1, is equal to this. This integral will turn out to be 2 a0 . The second part will drop off because if you see this a1 psi is an odd function and the integral of an odd function over this interval is going to be 0. My question is, is it equal to what I get out of the one point rule? Out of the one point rule f at the point psi equal to 0 is a0 into the weight at the point psi equal to 1 which is 2; this integral is exact. The one point rule not only integrates the constants exactly it also integrates the linears exactly. Let us become a little bit more greedy, and ask whether this integral can also be obtained exactly using the one point rule. This integral if you see it will be equal to by what we have done 2a0 plus this part will drop off; this part is an even part it will give me a value 2 by 3 a2. If I look at the one point rule, the one point rule gives what? One point rule will give me this integral equal to f at 0 into w1 which is equal to a0 into 2 which is not equal to I2. The one point rule is only exact for all polynomials which are linear, for higher polynomials it is not exact and it will not give you the exact value of the integral. So for this; let us go to the next step. We go to the two point rule. If I take the two point rule what do I know? The points are psi1 and minus psi1 and the corresponding weight for the two pair, because these are symmetric points, so both will have the same weight. This is my point if I want to call it as psi2 is equal to minus and the weight of w1 is equal to w2. What we want is this two point rule to be atleast exact for the quadratic. Because the one point rule could not do it so I expect the two point rule to do it. Let us take I2 which we got as 2a0 plus two third a2. This will be the function, the integrant evaluated at psi1 plus the integral evaluated at minus psi1 into weight 1. This is equal to fx psi1 plus fx minus psi1 will be equal to if you work it out 2 a0 plus 2 a2 psi1 squared. If we compare the two the coefficient of a0 and a2 I get this whole thing into w1. If I compare the coefficients then I will get 2w1 is equal to 2 and if I look at the second part it will be 2 psi1 squared w1 is equal to 2 by 3. From here w1 is equal to 1 and from here by substituting everything I get psi1 equal to plus-minus root of 1 by 3. So choosing these points psi1 equal to plus-minus root of 1 by 3 and w1 is equal to 1, the summation for the two point rule gives us an exact integral for I2. Our two point rule n is equal to 2 will be equivalent to saying I have the point psi1 equal to minus 1 by root 3 and the weight w1 is equal to 1 and the point psi2 equal to 1 by root 3 and the weight w2 is equal to 1. Let us see whether this two point rule can give us an exact integral for a cubic polynomial. What will the cubic be? This integral if I do it will come out to be exactly the same as the integral I2 or it will be given as 2a0 plus 2 by 3 a2 and we will see by taking the two point rule, I will get the exact integral for n is equal to 2. Let us ask if it will do the same job for I4 which is the integral of a fourth order polynomial. It is a0 plus a1 psi plus a2 psi square plus a3 psi cube plus a4 psi to the power 4 d psi. If I do this integral it will turn out to be 2a0 not plus 2 by 3 a2 plus 2 by 5 a4. And if I do the two point integration rule it will give me the two point integration rule it will give me 2 a0 plus two third a2 this one should check and this is not equal to I4. The two point integration rule can give us an exact integral for all polynomials up to a cubic. But for a fourth order it does not do the job. Then we naturally go to the three point rule. If I take the three point rule what will be the points? I will have psi1, psi2 is equal to 0, psi3 is equal to minus psi1, and here I have weight one, weight two, weight three is equal to weight one. If I look at this here is my psi is equal to 0 ; this is my psi2, this is my psi1 and this is my psi3 equal to minus psi1. If I take this and it should give us an exact integral for the fourth order one; I4 was equal to 2a0 plus two third a2 plus two fifth a4 what I want from the summation is equal to f at psi1 into w1 plus f at psi2 into w2 plus f at psi3 into w3. This will be equal to w1 f at psi1, plus f at minus psi1, plus w2 into f at 0. f at psi1 plus f at minus psi1 you see that this is also a symmetric quantity. This will be equal to 2a0 plus 2a2 psi1 squared plus 2a4 psi1 to the power 4. This will be 2a0 . If I match the coefficients of a0 a2 and a4 you see that by matching the coefficient of this won’t be equal to 2a0, f0 will be equal to a0 by what we have done. If I match the coefficient I will get 2 w1 plus w2 is equal to 2, this is my first relation that comes out. Then coefficient a2 I will get two third is equal to 2 psi1 squared w1 coefficient of a4 will give me two fifth is equal to 2 psi1 to the power of 4 w1. If I take the ratio of these two things this will give me psi1 squared, so divide this by this is equal to what will I get? psi1 squared is equal to 3 by 5; implies psi1 is equal to plus-minus root of 3 by 5. See it is very easy to get these integration point locations though as the number of point’s increases the task becomes more and more laborious. If I have that substituting this expression for psi1 squared from our previous equation we will get that 2 into psi1 squared w1 equal to 2 by 3. So psi1 squared into w1 is equal to 2 by 3, implies w1 is equal to 5 by 9 implies by our first equation. Since 2w1 plus w2 is equal to 2 implies w2 is equal to 8 by 9. We have found the locations of the points and the weights. So psi1 is equal to minus root of 3 by 5, w1 is equal to 5 by 9, psi2 is equal to 0, w2 is equal to 8 by 9, psi3 is equal to plus root of 3 by 5, w3 is equal to 5 by 9; this is three point integration rule. It turns out that if I go and check now for the fifth order polynomial, this was derived using the fourth order polynomial, if I find I5 which will be the fifth order polynomial for this also the three point rule is exact. If I write it down, the order of the rule and the polynomial order P bar for which it is an exact; the one point rule was exact for the linear. The two point rule was exact for all polynomials which were at most cubic; the three point rule is exact for all polynomials which are at most fifth order. This way we can continue, can we come up with a generic rule? Yes, we can come up; you see that if I have n, if I have a polynomial which is 2n minus 1 P bar equal to 2n minus 1, then the n point rule does an exact job. If P bar is less than equal to 2n minus 1, then n point rule is exact. Out of this how do I choose the appropriate integration rule? It is very easy now that you say that if I have a polynomial of order P then n should be greater than or equal to order P bar plus 1 by 2. If I use integration rule with number of points greater than or equal to P bar plus 1 by 2 then I know this integration is going to be give us an exact integral for that particular polynomial. Let us say that if P bar is odd then n terms have to be P bar plus 1 by 2. So n is equal to P bar plus 1 by 2 will do the job. I am going to fix my polynomial integration rule to n equal to P bar plus 1 by 2. If P bar is even, that is the order of the integrant is in even, then I would choose n equal to P bar plus 2 by 2. That is it goes to the next higher integer point. So with this choice of integration points I can expect that whatever integrals I have to do, if they are polynomials I will get an exact value of the integral. This is very important from our numerical computational point of view. Question is, can we construct these integration points easily through some other means? It turns out that yes, these points psii are the roots of Legendre polynomials defined in the intervals minus 1 to plus 1. Here we have an example of the Legendre polynomials that can be used to obtain the locations of the integration points. For example, if you look at the first polynomial it has its root nowhere in this interval from minus 1 to 1. Second one has the root at the point x plus 0 which for us is the point psi equal to 0. This is nothing but the point for the one point rule. If you look at the third one, the second polynomial which is a quadratic polynomial you see that the roots of this appear at plus-minus 1 by root 3. These are nothing but again the points - the coordinates of the integration points for the two point rule. The third one here can be written as half into x into 5x square minus 3. The roots are x is equal to 0 and x equal to plus-minus root of 3 by 5. This is exactly what we had derived using the long approach. Once I have these points then I can go back to what we have done there and obtained the weights for the integration rule. I can go on with the Legendre polynomials and I can get the integration rules. The location of the points corresponding to the integration rule of choice and it is given here in the form of a nice table, you see that if I have the one point rule there is nothing big that have to be done the location is psi is equal to 0 and the weight is 2. If I have two point rule here I have two points with plus-minus Xi as my coordinates. This is the symmetric points and the weights are this. Here the efforts have been made; this has been borrowed from a book, to obtain the points with sufficient precision; if you see 1, 2, 3, 4, 5 so up to 15 decimal places this has been obtained. That is very important from the computational point of view because we would like to do our computation in double precision. If you want to do our computation in double precision we should ensure that these points the gauss points, or the integration points that we have and the corresponding rates are also obtained in double precision. Very nicely all these points are given here you see for the three point rule you have these things and so on. Here in this table, I think the points are given up to n is equal to 10 to the ten point rule. 1 2 3 4, up to the 10 point rule and in whatever we will do in this course and most practical applications we do not need more than a 10 point rule. The question is, given this 10 point rule, you should be able to write a program where if I give you what is the n that I need n specified then you should return the corresponding points psii and the corresponding weights. Let us call this program or routine or the sub program whatever you want to call it as intake points where the input is n and the output is the values of the psii(s) and wi(s). Once we have this program where we input the n, I should be able to output what are the coordinates, all the integration points and corresponding weights. Then I should be able use this in the program. Remember we are going to use this in the finite element computations. More importantly, in the element stiffness matrix and the element load vector calculations; we have to first specify what is n. This question has to be first answered in order to obtain the n which will pass to this routine and this routine will return back the coordinates of the point’s integration points and the corresponding weights. How are we going to determine this n? Let us say that the objective of our program is to handle a material parameter which is linear at most. That is it corresponds to tapering bar. I can have this kind of a situation which I would like to handle, a tapering bar. We will say that our distributed spring support or the elastic support is constant . We do not want this to be anything more than a constant. And when we say constant and linear we mean that to be piecewise. That is I could have this kind of a tapered bar or I could have elastic support in only a particular region of my member . As long as it is satisfied piecewise we are happy. So k(x) has to be piecewise constant we are going to fix it which is what our code is capable of handling, the finite element program that we are going to develop. EA is linear and this distributed load f(x) you say it is almost quadratic. Again when you say quadratic it need not be quadratic in the whole domain, it may be quadratic in one piece, 0 in the other piece and linear in the other piece but not more than quadratic piecewise. If we have these things, then what is it that we have to do at the element level? If you remember that in the stiffness calculation at the element level we had these integrals we had 2 by h of the element d psi plus integral psi equal to minus 1 to 1, if I am are putting k to be constant, it will k0 into which a constant value Ni hat Nj hat into h of the element by 2 d psi. This is from what we have already developed. If I want to now fix the order of the integrants here what are the integrants? This is an integrant, this is an integrant. If you look at the order of the first one what is the order? The order of this one is 1 it is linear, the order of this one is P minus 1, the order of this P minus 1, the order of this one is 0. The total order of this integrant is equal to, for this one, I will have the order of the integrant is 2P minus 1; I simply sum up this order. So 1 plus P, 1 is P 2P minus 1. If I look at it here, this one of order 0 this is of order P, this is of order P, this is of 0. Sum of the orders this integrant is of order 2P. This order I will call it as P1 bar this I will call it as P2 bar. Again let us look at from the element k what is my load vector? Its load vector goes from minus 1 to plus 1 f hat function of psi Ni hat into h k by 2 d psi. By what we had done, this is of order at most 2 this is of order p this is of order 0. This one if I look at the order of the integrant here the order is P3 bar, I will call it, is equal to P plus 2. For our integration rule the order of the integration rule to choose that I will have to take that the maximum of all this polynomial orders, because I would like to integrate the highest order polynomial exactly. What I will get is the order q is equal to max of 2P minus 1, 2P and P plus 2. You see that q is equal to for P is equal to 1. If I take this one is 1, this one is 2, this one is 3. So 3 for P equal to 1. If I take P is equal to 2 then this one is 3, this one is 4, this one is also 4. If I take P is equal to 3, then this one is certainly less than this one the second one, second one is 6 and this one is 5. This is equal to 3 for the choice of the variation of EA F and K that we have chosen and the order of the polynomials that we are picking for the approximations. The integration order, the maximum order of polynomial for which I would like to have exact integrals is given by either 3 when P is equal to 1 have or 2P when P is greater than or equal to 2. Then our n is equal to q plus 1 by 2 will do the job; this is true when q is odd and this is equal to q plus 2 by 2 when q is even. You see that in our approximations apriory, even before doing our computations we can fix what is the order of the integration we need because polynomial order of approximation P the material data that is EA and k and the load data are available to us a priory as input data; using those data I should be able to determine n. The next task is if I know n, I know how to construct the integration points and the rates, given the n. I pass this n to my routine integ- points and it will output to me the integration points psii and wi. Once I have output these integration points in order to do smart computation; you remember that finally all this will be done by a computer. And the computer should not be subjected to repeated tasks over and over again. It is a good idea when I am using sufficiently large number of elements in the mesh and then to construct the shape functions at the integration points. How do we do this? We can use another routine which we can call as Shapall which takes as input the P that we have and which outputs the following things . These arrays I am giving it some name P, psitot, dpsitot where pistot can be an array of size the number of integration points into P plus 1 and dpsitot is also size n and P plus 1. What am I am going to do here? For each integration points psii I am going to load the value of all the shape functions of order P obtained at these integration point. Similarly, here I am going to load at each integration point the values of the derivatives of shape functions with respect to psi at this integration point. How do I do it? You start loop over integration points. For each integration point you have already obtained or programmed, the routine which will give you the shape functions at a given point. So here I will put psi equal to psii which is the integration point; call shape function routine to give Nj psi and dNj d psi at the point psi. For j equal to 1 to P plus 1, I am going to do the following: psitot I, j is equal to Nj at this point psi. Similarly, dpsitot i, j is equal to d Nj d psi at this point psi; then I can end the loops. This is a program which are a routine again which will take the integration points, the order of the shape functions that we want. Which order shape functions? It will go to my shape function routine return the value of the shape functions and the derivatives of the shape functions with respect to the master coordinate at each of the integration points and store it in a array where we give it some name. The first array is the array of the shape functions the psitot and the second one is the array of the derivatives dpsitot. Once we have loaded it then you see that all elements when you do the integration then we do not need to call the shape function routine over and over again before we do the integration. We can simply use these two arrays that we have created and these are not big, these are very small in size. We can simply use them over and over again in each of the element calculation. In the next class or in the next lecture we are going to extend this further. We have if you remember; we have now developed capabilities to get the shape functions and the derivatives at any given point psi. Similarly, we can give you the integration points given the order of the rule n and thirdly given the order P we can get all values of shape functions at integration points. We are going to use this information in the next lecture to obtain the suitable routine which can do the element stiffness matrix and the load vector calculations. Once we have that routine then we would have done a significant bulk of our one D computational program finite element program that we want to write in one dimension. This we start off in the next lecture. Thank you.