 14 lecture we saw finite element method procedure through the discretization method and there we saw how do we calculate the energy associated with the entire domain. We had taken three element example with 4 nodes and then we calculated energy and then we also saw properties of shape functions of the various nodes. Now if you really look at the procedure that we followed for calculating the energy, it was really not amenable for coding particularly for 2D or 3D. Why? Because we calculated you know energy for the entire domain at once is it not? We were integrating from 0 to 3L, you remember that? We are integrating from 0 to 3L and then depending upon whether shape function is 0, not 0 and accordingly we were, is it not? We were either taking some terms, we are not taking some terms. But that is not good from the point of view of coding. From coding point of view you know it will be easy if we go element wise, we calculate energy element wise and then form first element coefficient matrix and then combined element coefficient matrices. It is one and the same thing we are doing the same. Finally we are integrating the energy but we are first calculating the element coefficient matrices because we are calculating the energy of each element with contributions of shape functions and the corresponding potentials associated with only that element. And then we you know combine element coefficient matrices to form global coefficient matrix. That is what we will do now. So here that is why we are taking only first element say phi 1, phi 2 with node numbers 1 and 2 and this is element number 1. So now we will calculate only energy associated with only this first element. That is why it is here only integral 0 to L and then it is we know the what is the expression for energy it is phi dash integral phi dash square minus integral phi hdx is it not. So phi dash here is now we have seen in the previous lecture that your you know phi for this first element will be simply n 1 phi 1 plus n 2 phi 2 is it not. For the first element potential will be n 1 phi 1 plus n 2 phi 2 because n 3 and phi 3, n 3 and n 4 are 0 for element 1. So then simply for you know this will have only 2 term n 1 of 1 dash phi 1 plus n 2 1 dash, dash means it is derivative. Last time we discussed that derivative will be valid only for the shape function because they are the functions of x. phi 1 and phi 2 they are not functions of x because this is a variational procedure we are varying the values of phi at every x. So this will be the you know phi dash square and this is n hdx sorry phi hdx this is phi phi hdx and h stands for the source term. So now we will form the element level matrices. So now you will get a 1 1 and a 2 2 by using the same procedure this is not going to be the same. Instead we will what we have done is we have just one element level. So then you will get a 1 1 and a 2 2 of element 1 as 1 over l. So here in the previous case we had got you know when it was this a 2 2 total was contribution when we considered a node 2 the contribution was 2 by l is it not because that 2 was common to both element 1 and element 2. But here since we are taking only first element a 2 2 will be only 1 by l. When we take the second element again a 2 that for the second element this you know it will be the this node number will be corresponding to node number 1 of that second element is it not. So there again it will come 1 by l and that 1 by l and this 1 by l will get added to form to get 2 by l and that will be the total contribution of node 2 that will give the corresponding energy of because of potential at node 2 is it not is this point clear. So here that is why a 2 2 is only 1 by l because the other contribution of potential at node 2 in the second element is that to be accounted for that will get accounted when we consider the second element is it not. But this makes it very simple now you can see imagine if you go element wise in coding when you develop your own code it becomes very easy you just consider you know an element you set up a do loop and go element wise calculate element coefficient matrix and then later on you form the global coefficient matrix. So how do we do that that we will see and then this a 1 2 and a 2 1 will remain same as earlier minus 1 by l. So these are the four terms and b 1 and b 2 will again remain h l by 2 so there is no change. Similarly for now element 2 which will be you know so if this is 2 then this will be phi 2 and phi 3 with these this will be node 2 and this will be node 3 and then you will get n 2 dash phi 2 plus n 3 dash phi 3 here is it not and similarly here and then this will be a 2 2 2 and a 3 3 2 will be 1 by l and then you know you will have this a 2 3 and a 3 2 as minus 1 by l and then correspondingly you know these terms will remain same. So now if we go further and then similarly for element 3 it will be a 3 3 and a 4 4 now those will be you know again 1 by l so all these element coefficient matrices are identical in this case that is the advantage right. So in this case it becomes very simple. So the three element coefficient matrices will be 1 by l and this matrix similarly the b matrices which are representing the source term will be again h l into 0.5 0.5 right. So now what we have to do is we have to combine the element coefficient matrices to form the global coefficient matrix which will actually basically add the energies of all the elements together and then the contributions of the nodal potentials will automatically get added is it not. So then here you have this these are the three element coefficient matrices 3 a matrices are these three except you know 1 by l is not considered here 1 by l of course is there but that is common so that is not considered. So these three matrices now these are what are written here these are the node numbers these are sort of global node numbers because element number 1 is between 1 and 2 element number 2 is between 2 and 3 and element number 3 is between 3 and 4 right. So then you know you have in element number 1 this is a 22 for element number 3 this is a 22 is it not. So when we combine this a 22 and this a 22 has to get added because remember the energy overall energy is just one scalar number is it not it is scalar quantity. So all these things will get added to Q only one number energy so all these terms will get added so that means you have to you know collect terms with you know same coefficient so that is why this and this gets added that is why the you know the contribution of node 2 becomes 2 in the global coefficient matrix whereas contribution of node number 1 remains 1 because a 11 does not appear anywhere here except this a 11 similarly a 33 appears here as well as here so that gets added you get 2 a 4 4 is only appearing 1 that is why you get 1 here right and then you know you are off diagonal elements will remain as here so this is like a 12 this is a 21 so a 12 is minus 1 a 21 is minus 1 so that is corresponding to these two elements. So these how we have formed we have got the same global coefficient matrix as earlier but it was now so easy if we had gone if we if we go as per this element wise it becomes very easy so you have to calculate you have to set up a do-loop in when you code and you know you have to calculate element wise this coefficient matrix you know you have to just set up a do-loop and do it and then combine later on when we actually see the code will see how easy it is so just wait till that time okay so in fact in the next after this is over we are going to show you one code return in Sylab which will make things more clear similarly for the B matrix on the right hand side which is the source again you know it is 3 B matrices again HL is not considered here because it is just common so 0.5 0.5 so now again same thing this is B2 here also it is B2 B3 B3 so B2 and B3 0.5 plus 0.5 becomes 1 the contribution becomes 1 for node number 2 and 3 again why because those nodes node numbers 2 and 3 are common right and that is how you get this B matrix now these two matrices are same as earlier right okay so now we will see one more thing that is what we have seen till now this procedure was for Poisson's equation because we considered Poisson's equation now you know we will consider non-homogeneous wave equation where you have phi term here right and then this x comes on this side and that represent the source earlier we were taking H as the source so H was constant that means everywhere in the domain there is same source now here source is there and is function of x is it not so source is function of x it is not constant right and so and we have earlier seen this the form will be when we actually saw the global approximation where in you know whole domain approximation we saw there it was AC minus DC is equal to B remember that AC minus DC is equal to B where C where the unknown coefficients in that whole polynomial expression so that now what we have done we have gone from that approach whole domain approach to elemental level approach so now those C coefficients are replaced by the potentials at nodes so that is why that final equation becomes A phi minus D phi equal to B so now A is already evaluated in the previous lecture because this del square phi is there whether it is in wave equation or in the Poisson's equation or even in the Laplace equation when you follow the FEM procedure that will result into the same global coefficient matrix A if geometry and material properties are the same so what actually it needs to be you know what we need to form matrix corresponding to is phi because we have already worked out for del square phi that is A matrix in the previous you know 2, 3 slides we saw that is it not and the previous lecture so now we have to worry about this phi and then also about x because earlier the source was constant it was not varying with x now source is function of x right so now how do we you know form the D matrix so in the functional the extra term that corresponds to D matrix is half phi square DX remember the functional you take phi on this side and then you know multiply by phi that is the phi square integral phi square DX is the and multiplied by half is the term that appears in the functional remember the functional for this when it was you know we discussed this in one of the previous slide when you study previous lecture you will find that this becomes when this comes on this side it will be minus phi minus phi becomes phi minus phi square in the functional half minus phi square with a minus sign but minus sign will consider later when you form the this equation so basically we want to find the phi square half integral phi square DX so now this phi square for element number 1 element number 1 again it is n 1 phi is replaced by n 1 phi 1 plus n 2 phi 2 and square of that so we have to just evaluate this integral so it becomes simply this is it not if you expand this square you will get this so this integral expression can be further elegantly written in this form using matrices so what are these matrices and what are these particularly these elements we will see now so same expression is written rewritten here right so now that original x integral this integral is now rewritten here with one change instead of twice phi 1 phi 2 we have just split that into 2 terms phi 1 phi 2 into these 2 terms and again phi 1 phi 2 into these 2 terms right so now this d 1 2 in this matrix is given by integral 0 to L x 2 minus x upon L into x minus x 1 by L DX right so it is obvious because here if you expand this you know matrix equation matrix expression you can easily correlate that off diagonal term d 1 2 will be this right now if you evaluate this integral you will get the answer as L by 6 now these all these off diagonal terms they basically represent d 1 2 d 2 1 for element 1 d 2 3 and d 3 2 for element 2 and d 3 4 and d 4 3 for element 3 right the diagonal terms are given by for example d 1 1 will be given by this d 1 1 will be integral of this term so that is integral 0 to L x 2 minus x by L square d x which comes equal to L by 3 and in the corresponding you know if you evaluate it for various elements then you will get d 1 1 and d 2 2 for element 1 d 2 2 d 2 d 3 3 for element 2 and d 3 3 and d 4 4 for element 3 remember that this is a 3 element example with 4 nodes 1 2 3 and 4 that is the element level d matrix then will be L by 6 2 1 1 2 because this is d 1 1 d 1 2 d 2 1 d 2 2 and their values are L by 6 for diagonal off diagonal elements and L by 3 for diagonal elements now let us again go back to you know the expression that we saw in the previous slide so now this whole expression involving product of 3 matrices can be expanded to give you this now these 4 terms are just you know split here right now remember this is this is part of functional you know expression is it not this is part of functional expression so this is the functional is it not and this is corresponding to element number 1 so this is representing only energy of element 1 and the corresponding you know functional energy functional for the element 1 is written in this form same thing we are writing it here so this is the energy for element 1 that is split into 4 terms now in the process of FEM method what we are doing we are doing energy minimization so then we will have to differentiate with respect to first phi 1 then if you differentiate this with respect to phi 1 you will get this right and noting that d 2 1 is equal to d d 1 2 you will simply get this as d 1 1 phi 1 plus d 1 2 phi 2 right similarly differentiating with respect to phi 2 you will get you know d 2 1 phi 1 plus d 2 2 phi 2 so combining this then we will get the overall after minimization this whole term involving 3 matrices we will just reduce to d times phi and d matrix is d 1 1 d 1 2 d 2 1 d 2 2 into phi 1 phi 2 right so this is what after minimization you get only d phi now let us see the source matrix now just to summarize we started with this F the overall functional right and the overall functional was this and after minimization we would get a phi minus d phi equal to b a we have already seen which is a standard global coefficient matrix d we just saw how to derive 2 by 2 coefficient matrix d for element 1 and similarly it can be done for other elements now we will see how to evaluate these b matrix elements b matrix terms which represent a source so now the source term now again this is for you know first element that is it is 0 to L we are showing it only for the first element so the source term will be integral x phi dx again just going back here the source term is x here when you take it on this right hand side it will become minus x in the functional expression it will become x phi dx right so integral x phi dx so that is what we are doing integral x phi dx so now this again x into phi expression for phi is n 1 phi 1 plus n 2 phi 2 where n 1 this is the shape function for node 1 shape function for node 2 right and then again this can be written as elegantly in this form product of 2 matrices 1 row vector 1 column vector so now this b 1 small b 1 is given by integral from x 1 to x 2 now instead of 0 to L I have made here x 1 to x 2 so we can generalize although we are doing this for the first element then this can be you know if you are doing for the second or third element corresponding first nodal coordinate and second the coordinate corresponding to the second node can be you know put here at the same expression can be used so and that is why I am instead of 0 to L because 0 to L will be valid for the first element so I am generalizing it so x 1 to x 2 first node coordinate and second node coordinate of the corresponding element so b 1 will be integral from x 1 to x 2 and then this term x into x 2 minus x by L dx and this phi 1 anyway we are taking out in such matrix form so b 1 is given by this after you know evaluating this integral we will get b 1 expression as this similarly b 2 will be integral x 1 to x 2 x minus x 1 upon L into x this x dx and after the evaluation you will get this so now you can write this phi 1 phi 2 matrix into b 1 b 2 column matrix you can write it like this right so these expressions for b 1 and b 2 are written here so this is b 1 and this is b 2 after minimization that means when we differentiate it first with respect to phi 1 and then with respect to phi 2 you will simply get b 1 and we will simply get b 1 and b 2 here right so this is the b matrix so thus we have now understood how do we calculate all these matrices a d and b and these matrices particularly this a and b we have already seen in case of Poisson's equation but this d matrix would be required when you are solving either wave equation or diffusion kind of equation so now what we will do is we will now start with the explanation of code and here my PAD student Mr. B Sairam he has developed that code so he will explain that code to you as we have discussed in the previous slide we will be developing a one dimensional FEM code to solve a partial differential equation one dimensional partial differential equation that we have been seeing in the previous lectures so that one dimensional partial differential equation is as shown in the slide so which is nothing but del square phi plus phi plus x equal to 0 with same boundary condition so the same equation that we have solved using the whole domain approach and different approaches in the previous lecture manually we will be now solving it using a code so after applying FEM procedure this partial differential equation will convert into a linear system of equation or a matrix equation of this form we will convert into a matrix equation of this form a phi minus d phi is equal to b ok so in the previous lecture and then this will again simplified into this form so I am taking out this phi as common and then that will result into m minus d and that m minus d will I can I am calling it as k ok and then phi I am assuming it as u are getting so what I am doing is we are as we have seen in the previous lecture the partial differential equation will be converted into a matrix form and then that matrix form is then simplified into a x equal to b form here a is nothing but k and x unknown matrix is nothing but phi and source matrix is b ok so this is clear now now in the previous lecture we have seen an example of three element example is it not in the three element example we have three elements and four nodes ok so we have only three element coefficient matrices using which we formed a simple 4 by 4 matrix are you getting using three 2 by 2 element level matrices we have formed a 1 4 by 4 matrix and that 1 4 by 4 matrix we have reduced it to a 2 by 2 matrix after applying boundary condition. Since the number of elements are less it is obvious and it is simple to do manually but if we increase the number of elements then it will be very tough job to do it to form those matrices manually even though if you form matrices if it is a if it is simply let us say if I have chosen some five elements let us say with six nodes five elements and six nodes then my global coefficient matrix will be of size 6 by 6 and after applying boundary condition my size reduces to 4 by 4 then how can you take inverse of that manually doing that is very tough job so such kind of operations we can easily do by using scientific computing software called like Sylab, MATLAB, Python anything here in this course we will be seeing codes we will be we will be developing code based on Sylab only why we have chosen Sylab is like the syntax most of the syntax of Sylab and MATLAB are similar. So, even in the in the code that I am going to demonstrate you can see you can use most of the parts of that code in MATLAB directly only a few changes one or two changes only will be there the same code you can run in Sylab the same code you can run in MATLAB as clear. So, that being said we will now see what are all the steps involved in developing a code from the first lecture of FEM module we have seen FEM starts with choosing a geometry discretizing the geometry and then formulating element coefficient matrices and then joining global coefficient matrix applying boundary condition and solve it. You choose any partial differential equation that you want to solve using FEM these are the only steps any partial differential equation these are the only only in case of transient we will be solving it every time, but the procedure will be same. So, this code also I will be demonstrating in the same steps. So, in the next lecture we will see the complete coding procedure in Sylab. Thank you. So, we will stop here and continue next time. Thank you.