 So, welcome to the 19th lecture. In the previous lecture, we saw two-dimensional FEM procedure. Now, we will quickly see the coding aspects. Manual machine. We have already seen 1D code. So, this code we are not going into line by line and only part of you know code is will be explained in detail. So, now here the same the geometry is rectangular conductor in that boundary. Now, what we are doing is we are taking 51 by 51 grid effectively. So, number of nodes along x and y axis are 51 and 51 and then the total number of nodes that will be that will be 51 into 51. Then we are actually initializing global coefficient matrix. All the elements of C are made equal to 0 similarly B. So, source matrix is the entries are made equal to 0. Then we are actually this grid that we will eventually get because of this that actually we are getting the coordinates x and y coordinates of those all those grid points which is simply you know now I am not going into the line by line. It is obvious and you know students are generally good in coding. So, some simple things like this I will not explain. So, basically by using this command our geometry is 0 to 0.1 and 0 to 0.1 in both x and y direction that is a 0.1 it is 100 mm on both x and y direction. So, we are getting x and y coordinates of all the grid points by this command. Then you have this you know after you have formed the grid on x axis will have you know these are the nodes that you will get 1, 2, 51, 52, 2, 102 and so on so forth 2551 to 2601. This will be the node number. So, we are basically forming the P matrix. Now there are we have seen earlier there are three matrices P, E and T. P matrix has information of coordinates right. So, x and y coordinates here. So, now here given domain is discretized into squares and E square is divided into two triangular elements as shown here. One lower triangular element and one upper triangular element right and then we what we are doing we are actually filling the P matrix. Now we have already got the coordinates x and y all these coordinates we have got. So, we are basically populating the P matrix right. So, P matrix you will get like this. So, there are total number of nodes will be 2601 that is 51 into 51 remember 51 nodes. So, 50 segments right and then you will get P matrix as 2 by 2601. So, there are 2601 nodes. So, x and y coordinates of each of those nodes. So, for example, node number 1 it is 00, node number 2 it is 0.002 and 0 and so on. So, first you go like this then like this. So, remember we are doing manual machine. In the next lecture we will use gmesh which is a freeware machine software and days you know work of bookkeeping of element numbers and node numbers can get avoided by using that gmesh freeware software right. So, here again why we are going into this manual machine because then this you are coding fundamental we will get we will get consolidated in this. Then we can go to gmesh and then later on we are not get when we start looking at complex problems we are not going to get into coding for each of those problems. We would have by then explained you at least two or three course then you know you can develop course for any other problem also using such logic. Then next is the T matrix. Now, T matrix has information of where that you know what is the sub domain number and sub domain number of each element and the global node numbers of each element. So, for example, this T matrix will be 4 by number of elements is it not 4 will be rows in which first row first entry will be now this column each column represents one element and the first entry of that of any column is the sub domain number and then these are the three global node numbers of that element right. So, here then we will form the T matrix now what we will do is we will take the lower triangular elements first that means we will first take take all the lower triangular element go to the second row again lower triangular elements finish all lower triangular elements and then go to the upper triangular elements while forming the T matrix right. So, these you know again is you are running this you know index from i 1 to m minus 1 i 1 to n minus 1. So, you know you get the total you span the whole geometry and here we are taking same material everywhere right actually speaking you have two materials here is it not one is the copper conductor around the surrounding the copper conductor you have air, but copper and air both are non magnetic which you are equal to 1 is it not that is why we are taking we are talking of same we will give define them as same sub domain number, but actually speaking you could have divide you could have named conductor as said sub domain number 1 and air as sub domain number 2. So, that if you have if you wanted to derive something exclusively for the conductor to the core then you know you can for example calculate the energy associated with only material number 1 or material number 2 then in that case you have to define the two regions with different you know sub domain numbers is it clear. But here since we are doing we already know what are the coordinates of that central conductor and where it is lying. So, we do not have to you know define them separately by two numbers, but you will find in commercial software as you always you will always have different materials although their properties may be same they are always numbered by different separate numbers that we could have done here also, but just to reduce the coding effort we have named it one by one only because the material properties same Mu R is one for both air as well as copper practically because we have seen copper is diamagnetic Mu R is slightly less than 1, but for engineering purposes Mu R is one for copper. So, now when you execute this code so this sub routine is explained here that when for the first one when H is equal to 1 index is 1 these three commands will result into 1, 2 and 53, 1, 2, 53. So, these are the global node numbers for first element. Next time when you again go into this loop with H is equal to 2 you will get 2, 3, 53 as these three numbers. So, by this way you would have formed E matrix entries or you would have got T matrix entries corresponding to all the lower triangular elements in the whole domain. Similarly, we can do it for upper triangular element. So, these T matrix entries we already got it. So, 1, 2, 5, 53 so that is why 1, 2, 53 so and this set of entries are for the lower triangular elements. And now when you execute again the similar code with you know these commands little bit different than the earlier one you will now get when H is equal to 1 you will get 1, 53 and 52 as the output of this sub routine. So, 1, 53, 52 and similarly then when next time when you execute it will be 1, 2, 54, 53 and likewise. So, then at the end of these two you know codes two sub routines being executed then you would have got all the T matrix. Now, we will go further. Now this command a3n underscore elements is equal to size of T what will happen? Now T matrix is 4 into number of I explained you it is 4 into 4 by number of elements is it not? So, that is why it is 4 into number of elements. So, 4 will get assigned to a3 and number of elements which are there that is number of columns will get assigned to this variable n underscore element by this command. Similarly, you know you for the P matrix a1n nodes n underscore nodes is size of size P that will assign a1 as 2 because P matrix is 2 by number of nodes is it not? And number of nodes will get assigned to n underscore nodes and that is why and then when you display those number of elements will be displayed as 5000 and number of nodes as 2601. So, now we will go further and see how we can form element correlation matrices. So, here we run a for loop from 1 to number of elements and then we take in in nodes where in nodes will be now a matrix containing three entries corresponding to the global node numbers of that element. So, for each element then we will get in nodes those three corresponding global node numbers when we execute this statement because the T matrix 2, 3, 4 entries are the global node number and the first entry is the sub domain number is it not? So, by this command we would get global node numbers of each element like this for element number 1 it will be 1, 2 and 53 are the global node numbers. Having got those nodes we need to get x and y coordinates of those node numbers that now x and y coordinates are in P matrix because P matrix has information about the coordinates. So, by executing these two commands x c is equal to P 1 comma nodes and similarly y we will get in x c and y c the corresponding node the x and y coordinates of these three nodes. So, again x c and y c will be 3 by 1 matrix corresponding to x and y coordinates of these three global node numbers are corresponding to that element number under consideration. So, for each element number we will get this information which will be helpful to form element coefficient matrices as we will see later. Now, these equations we have seen earlier number of times now and safe functions n 1, n 2 and n 3 are given by this is it not? These three expressions. So, now we call this y 2 minus y 3 as P 1 x 3 minus x 2 as P 1 and similarly the other you know two sets of you know expressions then we will call this P 1 to Q 3 as you know given here to simplify our you know equations and we will just in place of this we will just replace them by P i's and Q i's right. So, these as I said earlier also these expressions would be required quite often when we do coding. Now, we want to find because see in our geometry we have this rectangular conductor and then this is a air box. So, as far as the material properties are concerned both air and copper they are non-magnetic. So, mu r is 1 but what is the difference between the two there is a current in this conductor. So, that current we have to impose right. So, we have to find out all those elements which lie within this conductor so that we can assign current density and the correspondingly we find the B matrix for those elements is it not? The x coordinate and the y coordinate of the centroid of an element is the average of the x and y coordinates of the three vertices respectively. And then we know this conductor the coordinate here the left most point and the right most point here top most point top and right and bottom and left most point they are 0.03, 0.03 and 0.07, 0.07 remember we had this whole geometry form 0 to 100 mm whereas this was from 32 70 mm is it not this conductor was from 32 70 mm in x direction. So, we have what we are checking we are checking whether the x coordinate lies between 0.03 and 0.07 similarly if y coordinate lies between 0.03 and 0.07 of that centroid then that element is within this you know rectangular area is it not? For example, let us say we have one element like this triangular element like this. Now, we are calculating the centroid of that triangle which is at the center of that triangle. So, that triangular that centroid if it lies within this then that triangle is lying within the conductor is it not? And that is how if that being the case then we assign mu r mu as 4.2 10 raised to 7 because mu r is 1 mu 0 into mu r is it not mu 0 into mu r and then mu r is 1 and b j is 10 raised to 3 this is just we are taking as an example we are considering current density in this conductor as 10 raised to 3 ampere per meter square. So, here just note what I have written current density or it can be ampere turn density because in a winding now this is like a this conductor is like a 1 ton winding you can say is it not? So, what actually you have to feed in is ampere turn density. So, when we will actually see you know transformer or rotating machines example. So, there we will be defining the current carrying area. So, there we have to define ampere turn density. So, ampere turns divided by the corresponding area. So, that is what is fed there and else means that element is outside this conductor area then again mu r is 1 but b j is equal to 0 because there is no current carrying part in this region. So, b j is 0 this we have already seen I think last time area of the triangle is given by this. You can verify this by using the expression that we have got now for p i's and q i's and then also absolute value of is taken since in case the node numbers are not numbered in anticlockwise that in case this also I have explained earlier and then I have also probably in one of the previous slides we in fact verified this by this expression which is this and formula for area of triangle in terms of determinant both come same this we have seen earlier. Now, we form element coefficient matrix. Now, see the element coefficient matrix is just formed by this one statement it is so simple in terms of coding effort is very very straightforward in just one statement we are forming you know element coefficient matrix. Now, this expression we have seen earlier c ij is this this we have seen. Now, same thing we are putting it here right c element ij see remember this is being run in the for loop for each element. So, suppose element number 1 corresponding to that you will have now 3 by 3 matrix is it not. So, Mr. Sairam explained to you about this thing earlier. So, this is like you know various stacks element 1 3 by 3 matrix element 2 3 by 3 matrix element 3 3 by 3 matrix right. So, these are all those you know stacks of 3 by 3 matrices that would get formed when you run the for loop for each element right. So, c say for example, first element i goes from 1 2 3 j goes from 1 2 3. So, it will be 3 by 3 matrix for element 1 and then the expression is this c ij is this as simple as that right. So, then b b e the corresponding source contribution to that element will be b j j delta by 3 this also we have seen last time. So, the current density if it if it exists in that element it basically gets into each of those nodes it gets approximation to each of the nodes of 3 vertices of that element triangular element equally as j delta by 3 right. So, that is what is being done here. So, this is what that j delta by 3 1 1 1 going further now like we saw in case of 1 d now we will see how we form global coefficient matrix by combining all element coefficient matrices right. So, again we run a for loop starting with 1 2 number of elements for each element we take element by element and go on upending the global coefficient matrix right because why we have to do that each element will have 3 nodes some of the elements some of the nodes of each element may be common to different elements is it not. Also some edges can be common to 2 elements each element we get in this 3 by 1 column matrix named as nodes the corresponding global node number using the t matrix right then we run this for i goes from 1 2 3 j goes from 1 2 3. So, those 9 entries will be here when i and j run from 1 2 3 now depending upon the corresponding global node number that each entry each of those 9 entries will get appended to corresponding entries of the global coefficient matrix is it not because here is nodes i and nodes j these will return the global node number because nodes are nodes is a 3 by 1 matrix. So, nodes i and node j will return global node number. So, for example here element 1 and 1 1 1 of element 1 in this case for example element 1 if you see our bring this in front of you. So, your first element is lower triangular 1 2 53 second is 2 3 54 then all 20 500 elements all lower triangular elements will get completed then the next one 2 5 0 1 element will be 153 52 and all such upper triangular elements right. So, the first element has global node numbers of 1 2 and 53 right and the nodes that is why it is 1 2 and 53 now when we run that these two power loops what we are going to get because this one is the global and that one is also the local node number of element 1 is it not. So, that is why small c 1 1 that means 1 is element its local node numbers 1 1. So, effectively it is c small c 1 1 right that basically gets upended to capital C 1 1 because this local node number 1 of element 1 is a in fact global node number 1 is it not. So, that is how you get and here initial value is 0 for this capital C 1 1 because first time we are putting something there right ok. So, then capital C 1 1 will get filled with small c 1 1 of element 1 similarly when i is equal to 1 and j is equal to 2 1 2 of element 1 will get upended to capital C 1 2 because again the global node number of this local number 2 is again 2 only right. So, it is C 1 so C 1 2 will get upended when i is equal to 1 and j equal to 3. So, the local node number of first element this is 3 but the global node number is 53 is it not. So, small c 1 3 of element number 1 1 2 and 3 are local node numbers. So, 1 3 will correspond to 1 and 53 globally is it not. So, capital C 1 53 will get the value of small c 1 3 of element 1. Now, you go to element number 2501 that means all after having gone through all the lower triangular element we come to the first upper triangular element which is 15352. Now, here you know when we run this or that element now c 1 1 is not 0 because this node number was already encountered when we got this element 1 and already you know you have in capital C 1 1 this value small c 1 1 of element 1 to that now what will get upended added small c 1 1 of element 2501 effectively we have seen earlier now this global node number 1 is common to this lower triangular element and this upper triangular element. So, that is why small c 1 1 will be added twice in formation of global coefficient matrix that is capital C 1 1 here for element number 1 and again capital when we again here it will get added to capital C 1 1 when we consider 2501 is it not. So, now whatever this 53 element 53 global node number will be common to how many elements it will be common to six elements is it not because there is one element like this second then third then this fourth fifth you can see that diagram and you will find this 53 will be common to six triangular elements. Since the global node number 53 is common to six elements the diagonal entry C 53 53 will get upended six times. But this whole thing gets done in just one this statement you may feel that is so complicated but it is just happening due to one this execution of this statement is it not. Similarly, for you know B matrix the source matrix here it is you know it is like one dimensional that is element number 1 and it is it is corresponding local node number. So, BE of 1 local node number 1 of element 1 will get upended to the global B 1 and initially it is 0 first time it will get added second time when it is encountered again for this element 2 5 0 1 again that node global node number 1 will be encountered right. So, again you know here the corresponding contribution will of that element will get upended to capital B 1 is it not and this goes on and this is by the end of this we would have got then both global matrices capital C and capital BJ and why we are calling this as BJ and not B because this has contribution of only source current source right. There will be another contribution and the corresponding modification when we impose boundary condition is it not that is why we are calling this as only BJ. Now, we know we apply the boundary condition. So, what we are doing the outermost boundary we are specifying a equal to 0. So, all these global node numbers which are on the boundary we have to specify they are a equal to 0 right magnetic vector potential and although it is magnetic vector here of course it is a scalar potential because I have already mentioned the direction is already fixed it is two dimension problem. So, A is in only z direction. So, direction is already known what needs to be determined is only the magnitude of A at every point and every point means at each of these all nodes right because we are converted this problem into you know a problem wherein we need needed to calculate potential at every point in the domain to a problem wherein we need to initially calculate potentials at only this grid point. Once you calculate grid point potential then potential at any point within the in any element can be found out by using that formula A 1 A plus B X plus C Y that original our approximation is it not okay. So, now we want to make the diagonal this also was explained to you when we saw the 1D code we have to impose the conditions for example this A 1 equal to 0. So, what you do is you make the corresponding diagonal entry 1 of diagonal entry as 0 right easiest way to do is first you make all the entries as 0 and then make the diagonal entry 1 rather than you know searching which is the diagonal entry and you know this makes it simple. So, this is straight forward is it not. So, C of I I means the diagonal entry is made 1 and all the orthogonal entries then already are becoming 0 by this command and then right hand side you have to make it 0. So, for the last time we I think saw one example A 5 equal to 0 in 1D code is it not do you remember that. So, the diagonal entry 5 we made it 1 and on the right hand side we made it 0 all of diagonal entries also are already 0 there. So, same thing we repeat for so this was bottom edge right. So, it was you know we are running from I 1 to 51 because n and m are 51. So, 1 to 51 we are running then the second edge vertical edge we are taking from only 52 to 2500 because 1 is already assigned and then we are adjusting this index. So, that we get this thing run from only 52 to 2500. That is why you know you have to little bit keep this book keeping because then we need to remember that. But when we use you know next lecture we will see use of freeware called that gmesh then it becomes all these things become very easy right. We do not have to worry about this keeping track of what is the you know node number global node number and where we want to impose what and. So, then same thing procedure that we use for this bottom edge is followed for vertical edge then vertical edge second vertical edge and finally the top edge. So, same procedure is followed only here index for this I is adjusted so that we get proper this number of nodes when we execute this form. So, now we are nearing the you know end of this code. So, now this bj although we are calling this bj now this bj got modified with the boundary conditions is it not bj got modified using boundary condition. Now, it has got contributions from both j source and boundary condition. Now, we know we have this matrix of the form ca is equal to b is it not we have seen this earlier and then what we have to do is we have to invert now c. So, a is c inverse b now in Matlab and Sylab both c inverse is given by this backslash c inverse is basically backslash. Now, we have got solution in a by doing this by executing this command we have got magnetic vector potential at all the unknown all those points within the geometry except boundary point because there we had already imposed the condition. Now, in order to plot we have to convert this a which is a column vector into a two dimensional vector because our geometry is two dimensional is it not we do not that it is of no use for us to have a column vector for the plotting purpose because we need to plot it now as a two dimensional in the two dimensional space. So, that is what we are doing now we are converting this column vector into n by n matrix and similarly two rows of p matrix which are basically coordinates. That means what we have to do at each of the coordinates now again p is a matrix of x and y coordinates. So, that we have to convert it as a two dimensional my matrix. So, for each of those x and y coordinates in two dimensional plane now we have to plot the corresponding a value which is also now obtained those that two dimensional values are obtained by converting the corresponding column vectors into n by n matrix. For example, here by this command vn matrix and this command here. So, we can we are converting this column vector into n by n grid. Similarly, you know this xx and yy by executing this command we get x and y coordinates in that n into n grid. So, similar command in MATLAB for example is reshape. So, those who are having MATLAB commercial software they can use the same code with you know changing this to reshape. So, having done that now it is a we need to plot. So, if you want to have a surface plot then you get remember a is equal to 0 at the boundary and a will be higher at the center. So, by using this search command xx, yy that is the now two dimensional grid and vn potential. If you want to get a counter plot then contour xy, vn and 25 means 25 number of contours. So, this is how you get it. Thank you. So, we will end this 19 lecture here and then we will see how the same coding can be done using gmesh meshing software which is a freeware. Thank you.