 Welcome to the 16th lecture. As we have discussed in the previous lecture, we will be developing a FEM code in Sylab. So, Sylab is a freeware software. So, you can download it from any place, you can install it and you can use it and it is similar and the syntax is similar to MATLAB. So, the same code you can use in MATLAB as well. So, whatever may be the codes that we are demonstrating here that can be used in MATLAB as well ok. Now, going further into the code. So, as I told in the previous lecture, the steps involved in the code will start with discretization of the geometry. So, first we will discuss how we will be discretizing a one dimensional geometry ok. So, here you can see the first two commands are CLC and clear off. What does it mean is? So, you will be using the scientific computing software very oftenly some previous code data may be saved in the console. So, that has to be erased. Otherwise, if the previous code variables and this code variables are common, they may interfere each other and that will give some error results ok. That is why it is better at the starting of every code you can add these two statements CLC and clear off that will make you safe and then discretization. So, one dimensional geometry means just a straight line, is it not? One dimensional geometry is just a straight line. A straight line can be divided into segments ok, segment. How a segment is defined? Segment is defined by its starting node and ending node, is it not? So, discretization of a one dimensional domain involves nodes and segments. Are you getting? Is it clear? Now, you can and one more thing I want to tell here is this is not the only way to code. The code that I am demonstrating here is not the only way to code. You can do it in n number of phase ok. So, this is one way to code ok see that I want to make clear. So, starting with discretization here you can see that the n is equal to 11 is the first step. So, this I am defining as number of nodes. So, discretization of one dimensional domain you can start with either defining segments or defining nodes, but that algorithm will be different this algorithm will be different, but here I am choosing the way to start with number of nodes choosing the number of nodes. So, first I am choosing n equal to 11 ok. So, that n equal to 11 will discretize my one dimensional geometry into 10 segments ok. In one dimensional geometry if it is a in one dimensional geometry which is a straight line if you have n number of nodes that n number of nodes will divide the geometry into n minus 1 segments. Are you getting in the previous lecture also you can see there are 4 nodes and number of elements are 3. This is a simple logic. So, once you define the number of nodes your number of segments are also defined. So, the number of segments will be n minus 1 ok. So, here this is how the discretization will be and then I am defining the global coefficient matrices. So, in the matrix equation what we have seen in the previous lecture here there are 3 matrices is it not? a, b, b. So, now our job is to fill these matrices. So, these 3 matrices I am defining first. So, a is equal to 0 of n comma n. So, 0's of n comma n this command will define the matrix size and it will make it as 0's completely. Now, those that matrix we will update with the element coefficient matrices ok. So, as we have seen in the previous lecture size of element coefficient matrices will be in case of a and d it is n by n that is why we have used n comma n in case of source matrix it is n comma 1 ok that this you can correlate directly with the previous lectures matrices. Next thing is what are all the coordinates of the discretization coordinates in the discretization ok. So, there are points what are all the positions of those points ok. See, here I am choosing the coordinates to be equi-space is it not? Why because this will simplify my discretization you can you can choose more number of nodes here more number of nodes here less number of nodes in between that can that is also possible, but that will complicate your code that involves choosing of coordinates. So, that is why here for simplicity I am choosing nodes which are equi-space ok. Let us say in case of for n equal to 11 this this command x equal to 0 to 1 by n minus 1 is a column 1 what does this mean is it will divide it will give the points in between 0 and 1 with a step size of 1 by n minus 1 this is the step size the middle one this you can see now the middle one 1 by n minus 1 is the step size. So, this command will result in x equal to a row matrix whose elements are 0 starting from 0 with a step size of 0.1. So, you got the positions of 11 nodes ok. Now, once we form the details of nodes where the positions are we need to know how these nodes are connected it is obvious in case of one dimensional domain they are connected by a straight line, but in case of 2D it will be difficult. So, that is why we that is why here also we are framing a T matrix to make you familiar about what is T matrix ok. So, T matrix is called as connectivity matrix. So, this matrix contains the information how each node is connected to other nodes. So, you can see this is my disk size domain the numbers which are encircled are element number ok. So, for the first element number the starting node is 1, ending node is 2 in the second element starting node is 2 and ending node is 3 this is 3 and for the third element starting node is 3 ending node is 4. So, what is the common thing here the element number is same as its starting number ok. This will simplify my coding steps ok. So, for each and every element I am assigning its starting number and ending number and ending number is plus one of its element number. So, that is what we have seen here ok. So, here the T matrix. So, here for the T matrix for each element the starting node number is nothing but its element number here I denotes the element number ok. Starting node number is element number and ending node number is 1 plus its element number. So, with this for loop for every I, I means element number for every I we will get starting node number in the first row and first row ok. I hope you are familiar with the matrices indication in case of Sylab and Matlab which is same in both the cases ok. So, the first index is row number, second index is column number. So, I is going to the column number that means element number is nothing but column number and first row will have starting row number and second row will have ending node number. Ending node number is 1 plus element number. Is it clear up to here? Now, how this matrix is a connectivity matrix where the connectivity information is involved in this matrix. So, you can see in this matrix. So, if you run this for loop, if you run this for loop we will get this T matrix as this form where the column numbers are element numbers and first row is starting node number and second row is ending node number. So, the entries in this T matrix are 1, 2, 2, 3, 3, 4, 4, 5. Is it not? Now, by looking at this we can know the connectivity. Is it not? So, the ending node of first element is nothing but the starting node of second element. So, through second node, first element and second element are connected. Similarly, through third node, second element and third element are connected. Similarly, through fourth node, third element and so on. So, such information is involved here and this T matrix is very important that in formation of element coefficient matrix and in converting element coefficient matrix to global coefficient matrix we will be using this very often. So, we need to know each and every part of this connectivity matrix. And one more thing I want to highlight is, so here the first row indicates the local node number 1 and second row indicates the local node number 2. That means for element 1, local node number 1 is 1 and local node number 2 is 2. Similarly, for element 2, local node number 1 is 2 and local node number 2 is 3. Are you getting? So, first row of T matrix indicates the elements, first local node number of each element. Second row of T matrix indicates the second local node of all elements. It is clear now up to here. So, we have completed the discretization part which is the first step. Now, the second step is formation of local coefficient matrix. So, in the previous lecture also we have seen local coefficient of matrix of every element is calculated individually. And since there are only 3 matrices, we wrote it side by side. Let us say we have 100 element, then we will get 100 local coefficient matrix. So, that means what in case of code you need to define 100 variables to store those 100 matrices which is very difficult and you cannot make the code generalized also. Is it not? If suddenly I want to use 1000 matrices, I have to change my complete code for 1000 variables which is very difficult task. Is it not? So, that is why here there is you know matrices with rows and columns, you know matrices with only rows, you know matrices with only columns. But in both Syllab and Matlab, we can define three-dimensional matrices. Are you getting? We can define three-dimensional matrices. That means n number of 2 by 2 matrices, n number of rows by column matrices, we can stack them one after the other in a single variable. That is how we can store all the element coefficient matrices in a single variable even though if you change the number of elements. And this logic is same in case of 1D, 2D as well. So, in the two-dimensional, in the code of 2D also, we can see the same three-dimensional matrix and the three-dimensional matrix will look like this. This is the three-dimensional matrix. So, element 1, the local coefficient matrix of element 1 will be stored like this. And after that element 1 matrix, element 2 matrix will be stored. After that, elementary matrix will be like that, all the matrices will be stacked like this. Now, going further, so now we will see how these element coefficient matrices are formed and how they are stored. So, we have seen in the previous lecture the formula for element coefficient matrix is this. We have seen in the previous lecture only. So, it is 1 by L 1 minus 1 minus 1 1. So, by looking at this, what are all the common things we have? That we have to identify first. That identification of such things will make your coding simpler. So, here what are the common things? Diagonal elements are just 1 by L and out-diagonal elements are just minus 1 by L. Is it not? So, here this complete thing is constant. Only thing you need to calculate is L for each segment. Since here we have chosen the equispacial dissipation, this L is also constant. But if you change your dissipation for non-equispacial thing, then this L will be a variable. So, to be more general, I have used this as variable only. This will be calculated automatically for every element. Now, we need to get the geometrical details of the corresponding element. So, here geometrical details means nothing but coordinates only. Is it not? Geometrical details means nothing but coordinates only. So, I have to get x coordinates of two nodes of each element. So, for the first element, for element equal to 1, local node number 1 is 1 and local node number 2 is 2. That we have seen in the T matrix explanation. This you have to be very much clear. See here, local node number 1 is 1. First row indicates the local node number. Second row indicates the local node number 2. So, local node number is 1. Local node number 2 is 2 in case of element 1. So, here T of 1, element will give 1 here. This T of 1, element will be 1 and T of 2, element will be 2 for in case of first element. So, x of 1 is nothing but 0, x of 2 will be 0.1. So, these two lines of code will give you like this for different element. We are getting the coordinates of each and every node inside the element using the connectivity matrix information. So, connectivity matrix information will give you global node number of each local node. Are you getting? That will give me the coordinate and then I am calculating length using ABS xn of 2 minus xn of 1. See this ABS is absolute value. See here, my second node is greater than my first node. This will give me positive value. If someone start naming the nodes in element from higher node to lower node, like if I choose my starting node as 2 and ending node as 1, this will go to negative. Is it not? To avoid that problem, I have used ABS. So, this will give me always a positive value because length should be positive. Now, length we got, only thing left out is to form this matrix. So, this one basic thing I want to highlight here. If you want to form a matrix whose dimensions are which is having rows and columns, you need to use two for loops. One for loop is for rows, second for loop is for column. If you want to form a row vector or column vector, you need only one for loop. So, since this is a two dimensional matrix, I am using two for loops. Are you getting? So, this i indicates rows and j indicates columns. If i equal to j, then a of element, i, j is 1 by L. If else, it is minus 1 by L. So, here I want to highlight the first difference between Sylab and MATLAB. In case of MATLAB, if you are familiar with MATLAB, this then is not required. In Sylab, you need to add this then. This is one of the differences I want to highlight. In this complete code, this is the only difference. When you use this code in MATLAB, you have to delete this then. That is it. If you use in Sylab, you can directly use this code. Like this, I have formed the element coefficient matrices in single variable, which is a three dimensional matrix, which will be stacked like this. The same two by two matrix which is having constant entries can be entered directly in a matrix form of this type. So, here if you enter the numbers in a row, one after in the bracket, one after the other, then the entries will go inside a row. And to change the row number, you just have to keep a semicolon that will go to the second row. Like that, if the matrix is simple with constant entries, you can use this way. If the matrix, so like, so in entire FEM procedure, either it is 2D or 1D, this D matrix will be simply constant number. We will have simply constant number. That is why in both 2D and 3D, sorry 2D and 1D FEM formulations, we will enter the element level D matrices like this directly. L by 6, bracket 2, 1, semicolon 1, 2. Now, we have completed the formation of element level D matrix. Is it clear? Now, we have to form the element level source matrix. So, we have seen in the previous lecture, the element level B matrix expressions are these two. This corresponds to node number 1, local node number 1 and this corresponds to local node number 2. Here, you can see these expressions are little complicated and they are not common to both the nodes. They are different for different nodes. Is it not? So, that is why we are entering those entries separately one after the other. For the first node, I am entering separately. For the second node, I am entering separately. So, here I want to highlight if your functions are common for all the nodes, if your functions are common, then you can use for loops to form the matrix. Otherwise, you have to enter one after the other like this. So, this is the third way of entering, but it will be used if you have a complicated function in x and y. So, now, as we have discussed earlier, the element level B matrix is a 2 by 1 or a 1 by 2 matrix. Is it not? We have seen in the previous lecture. It is like this, 2 by 1 or 1 by 2. So, to store all these element level matrices in a single variable, we can use simply a rectangular matrix. So, here I am using a rectangular matrix of this kind. So, here in each row, I am entering the element level B matrices. Are you getting? So, here this first row corresponds to the first element. So, here you can see the row index corresponds to element. When element is when element number equal to 1, that will correspond to the first row, local node number 1 entry. When element number 2, when element number is equal to 2, that will correspond to the local node number 1 entry for element 2. That will go in the second row. So, like this, we will completely enter the B matrix information in a single matrix. This is simply a rectangular matrix. Is this clear up to here? So, here you can see the first row corresponds to the first row corresponds to first element. Here this is the element index. Second row corresponds to second element. Third row corresponds to third element. And this first column corresponds to first local node and second row corresponds to second local node. Like this, we will enter the B matrix. Up to here, we have completed the formation of element level coefficient matrices. So, we have completed the second step. I hope this is up to here it is clear. Now, we will now we will go to the third step which is formation of global coefficient matrix using element coefficient matrices. So, actually he till now we have n minus 1 element level coefficient matrices. We have to combine all the element coefficient matrices into a single global coefficient matrix. Is it not? So, that we will see now how we will do. So, here we are forming the global global level coefficient matrix using the connectivity matrix ok. So, here we will use the connectivity matrix information to do that ok. So, I will just recap what is the connectivity matrix. Here the first row of the connectivity matrix is corresponds to the local node number 1 and second row of the connectivity matrix corresponds to the local node number 2. It is clear enough? Now, for each and every element as I have told earlier every operation in FEM we will do element wise. So, we will start with for element 1 to n minus 1 ok for element 1 to n minus 1. We are taking the global node numbers of each element global node numbers in each element ok. So, here the size of nodes is in this case it will be 2 by 1 first row corresponds as I have told earlier first row corresponds to local node number 1 second row corresponds to local node number 2. So, the nodes matrix the entries of the nodes matrix in case of element 1 will be 1 2 in case of element 2 it will be simply 2 3. What does it mean is the local node number 1 of element 2 is 1 local node number 2 of element 2 is 2 for element 1 it may be confusing. So, I will repeat the same thing for element 2. The local node number 1 of element 2 is 2 local node number 2 of element 2 is 3 ok this you have to understand. Next what is the logic that we are using here is the information already we have defined the global coefficient matrices is it not A D B if you see in the first slide very first slide we have already defined the 3 matrices. Now, we have to update these 3 matrices with the entries of element coefficient matrices. Element coefficient matrices also we have formed ok. So, we know the corresponding global node numbers of each local node. So, we will take the entries corresponding to each elemental matrix and we will update the update it in the global matrix corresponding to the position in the global matrix ok. Since our element coefficient matrices are 2 by 2 matrices and 2 dimensional matrices we are using 2 for loops ok. And this is a 2 by 2 matrix which has 2 rows and 2 columns here we are running from i equal to 1 to 2 and j equal to 1 to 2 ok for i equal to 1 to 2 and for j equal to 1 to 2. If it is a 2 dimensional FEM analysis which has triangular element whose nodes are 3 then it will go from 1 to 2 and 1 to 3 ok that we will see later. Now, for this A matrix what does it written here A of nodes i comma nodes j is getting updated with see here in the right hand side in the left hand side and the right hand side if we have the same thing and we are adding it with some other number that means we are updating it are you getting here we are just writing x equal to x plus some x 1 ok. Here this is these 2 are same is it not we are just adding this number with another number and we are keeping that same number in x that means we are updating. So, here this statement means the this statement means we are updating the position of nodes i comma nodes j with this number ok are you getting. So, what does this do that I will explain here for element 1 it will be very a bit confusing. So, I will explain it for element 2 for element 2 let us say this all code started from here. So, in the first loop let us say we have taken for element 2 ok for element 2 what are all the nodes nodes in element 2 they are nothing but 2 comma 3. So, local node number 1 is 2 and local node number 2 is 3. So, nodes are nodes is equal to 2 entry for i equal to 1 and j equal to 1 this is the first step that will be running for element 2 for i equal to 1 and j equal to 1 A of 2 comma 2 that means nodes of 1 is what nodes of i is here 1 nodes of 1 is 2 nodes of j is also 2. So, 2 comma 2 position of A is getting updated with second element 1 comma 1 position is it not second element 1 comma 1 position is this that 3 dimensional matrix that I have told. So, second element 1 comma 1 position is getting added to 2 comma 2 position of global matrix here I have written not equal to 0 because this got updated in element 1 ok. Now, we will go to element 1 this is clear now what we have done here ok. So, we are just so we know the global coordinates we know the local coordinates we know the we have the local entries we are updating the global positions positions in the global coefficient matrix with the corresponding local matrix entries ok. Now, for i equal to 1 and j equal to 2 ok nodes of i is what 2 nodes of j is 3 for i equal to 1 and j equal to 2 the second step are you getting. So, nodes of i is 2 and nodes of j is 3. So, 2 comma 3 position of A matrix is getting updated with 1 comma 2 position of second element ok. Similarly, for i equal to 2 and j equal to 3. So, like that we will be updating the A matrix similar procedure will be followed for D matrix as well any doubts here because this will be there in every code this part of the code will be there in every code to convert local coefficient matrices into global coefficient matrix this is the code we will use in everything only thing is here the numbers will change ok. So, here we have converted A and D matrices from local coefficient matrix to global coefficient matrix. Now, we will see it for B matrix ok. So, element level B matrix will be like this for element 2 similarly, we will see it for element 2 for element 2 nodes will be 2 comma 3 ok for i equal to 1 are you getting for i equal to 1 we are at node 1 element level sorry local node 1. So, local node 1 of B matrix B of 2 comma 1 see here 2 represents element. So, 2 whether the position of element 2 it will go to the position of element 2 and chooses the entry corresponds to node 1 and it will add it and it will update it to global entry 2 because global node number is 2 is it clear. So, here the complexity has reduced it to 1 because this is a 1 dimensional matrix this is a 2 dimensional matrix that is why here we have only i here we have i and j similarly for i equal to 2 for i equal to 2 global node number is 3. So, third position in the global matrix is getting updated with second with the entry of second node in element 2 ok are you getting. So, global position of 3. So, here local node number 2 is 3 ok global global node number is 3. So, third position of global matrix B is getting updated with second entry of element 2 element 2 data is here. So, second entry means this. So, this will get updated for this ok this is how we will convert the source matrix also. So, this these three steps inside this for loops will convert our element level coefficient matrices into global coefficient matrices ok. So, once we get the global coefficient matrices the next step is simplifying the matrices and solving the matrix equation. So, now as we have seen in the previous lecture we will reduce we got A U minus D U is equal to B and we are taking U taking out U as common that will reduce it to A minus D into U is equal to B and K I am assuming K is equal to A minus D. So, I am calculating K equal to A minus D ok and then now I will be solving K U equal to B after imposing boundary conditions. So, in the previous lecture we have seen imposing boundary conditions means we have reduced the size of matrix from 4 by 4 to 2 by 2 is it not. Also those boundary nodes are placed at first node and the last node the which we know, but in case of complicated geometries and complicated systems we do not even know the node numbers for which the boundary conditions will be applied. So, it will be very difficult to do such reduction in matrix size. So, in that case we what we will do is here the further row that corresponds to the boundary node we will make all the entries as 0 and we will make the diagonal entry as 1 that means we are imposing the boundary conditions. See here what is the step that we have followed for the first node and the nth node we are applying boundary condition. So, for the first row we are making all the entries as 0 is it not here and then we are keeping its diagonal entry as 1 and we are making its diagonal entry as 1 and as boundary condition is 0 we are we are updating that B matrix entry corresponding to first node as 0. Similarly for nth row because nth row corresponds to nth node. So, we are making nth row as completely 0 and we are making its diagonal entry as 1 and we are imposing the boundary condition in the B matrix. So, this can be visualized with this following example. This example we have seen in the previous lecture is it not. So, this 4 by 4 matrix you have converted into 2 by 2 matrix this is the 2 by 2 matrix this part is the 2 by 2 matrix. So, instead what I am doing here is here as according to this step first I made this first row completely as 0 and I kept 1 at the diagonal entry. Similarly, for the nth row I am making complete row as 0 and I am making this diagonal entry as 1. If you solve this equation simultaneously that means you are solving this equation the these 2 by 2 matrix equation how. So, here there is a typo this is u 4 this is this should be b 4 this should be u 4 ok. So, here if you see if you convert if you multiply this matrix with this this will give u 1 is equal to 0 and u 4 is equal to 0 that means in this system you have you have imposed the required boundary condition. Here if you want to visualize if you want to correlate this method with the previous method of reducing this you can utilize like this. So, if you have done some matrix thing Gaussian reduction those inverse taking methods. So, I can add this row to this row. So, I mean R 2 row 2 is equal to row 2 plus row 1 if I do that then what I will get 1 0 0 0 this minus 1 plus 1 will be 0 and 2 minus 1 0 ok. And I am making row 3 is equal to row 3 plus row 4 I am updating row 3 also. So, that will be 0 minus 1 2 0 ok and this is nothing but 0 0 0 1 that means these two these set these two equations got independent and this is similar to the matrix that we have seen in the previous lecture ok. This is this is how you can visualize this way of imposing boundary condition as similar to the wave what we have seen in the previous lecture also I want to make a note that making the row corresponding to that node as 0 and assigning diagonal entry as 1 this procedure we will follow in the 2D code as well ok. So, since here only we have only two nodes we have done separately when we have many number of nodes we will use a for loop for this as well ok. Now, we have we have imposed the boundary conditions also now we are ready to take the inverse ok. Now, after taking the inverse so, this inverse is the inverse function is there in all scientific computing software like Matlab and Seiler. So, you can use directly this k backstroke v is equal to u then u plot x versus u you will get the variation as shown in the red color plot this red color curve which is shown corresponds to 1D FEM ok. Since we have chosen less number of elements only 10 elements you can see that there is an error in between the elements. So, in between the nodes there is an error why because here we are minimizing the energy at nodes in case of FEM not over the whole domain ok. So, that is why there is a minimum error at the nodes ok. Now, I have compared the same thing with exact solution this blue one is exact solution and this black one is the approximations approximated solution with second order approximation which we have seen in the previous lecture. As you know this solution is for the is for the one dimensional PDE phi double dash plus phi plus x equal to 0 with boundary conditions phi equal to 0 at x equal to 0 comma x equal to 1. So, you can see here at 0 and 1 the boundary conditions are imposed perfectly ok. And this is the exact solution this we have seen in the previous lecture and this is the approximate solution with second order approximation this also we have seen in the previous lecture. So, in the previous lecture we have seen to improve the accuracy we have we have increased the order of approximation from second order to third order that made us to redo the complete process again ok. So, which is which is more complicated calculations compared to second order calculations ok. The involved mathematics is more compared to the second order is it not. So, as we as we keep on increasing the order the complexity in calculations will increase and also it is not a general procedure also whereas, the code that we have developed here I have shown it for 10 elements. Now, in the code I will change it to 100 elements by changing only the first the very first command this one here I will change this to 21 and then 31 and I will not touch any of the steps inside the code ok. So, what I am what I am trying to say is FEM procedure this kind of coding will be will be general only thing you need to change is discretization ok that I will show now ok. This is the CYLAB interface and this is the code in CYLAB in CYLAB interface. So, it is ok the same code what we have seen in the presentation is here I will be changing only this n number previously we have n equal to 11 I am running it for n equal to 11. So, you can see that solution in red color which is very much approximate and more error ok are you getting. Now, I am changing it to 101 21 first you can see both solutions got overlapped with little error. Now, I am improving it to 31 it got almost overlapped ok. So, the complexity in computation is also got reduced is it not one we have wrote we wrote code for 11 elements and then we got solution for whatever it may be I can make this for 1000 also it got overlapped you cannot see that red line blue line is overlapping that red line you can see red line here now if I very much if I zoom in too much ok. So, like this if you code you can make it make things more general and it will help you in improving your accuracy in computations. Thank you.