 So in the previous lecture we saw properties of shape functions and using shape functions how we can define potential at any point in a triangular element using the three shape functions at the node. So we will go further now, we will take, we know already this is the functional for two-dimensional case what we are saying A is again magnetic vector potential and this is the source term. Now here u is taken here and j is on the right hand side del square A equal to 1 upon u del square A is equal to j minus j that is why j gets multiplied by A in the functional expression. So this is the functional. Now we already have this you know various expression for example here A is substituted as summation ni A i because A at any x within the element is n1 A1 plus n2 A2 plus n3 A3 that is what we have substituted is it not. So similarly here this A is substituted here by this expression right is it clear and then we will call this as F1 minus F1 right. So now here we will actually what we will do is this del operator we will actually make it operate only on n1 n2 n3 because they are the functions of x and y right. A i is not function of x and y because A i at ith node in general is varied to minimize the energy. We have seen this in variation calculation is it not. So A i is not function of x and y right. So that is why this A A1 A2 A3 they are got out is it clear okay. So now then this F1 becomes only this. Now if we go further now this is the square the whole thing square. So that can be you know simplified or rewritten using this identity we know A dot A is equal to mod of A square is it not. So that is why then this A square of this whole thing can be written as this expression multiplied with this as dot product. So this A dot A so this is what you know this whole thing square is split into this vector multiplied to itself right. So now this when you actually expand this dot product you will get how many terms 9 terms is it not because these are 3 and 3 so 9 terms. So those 9 terms can be written nicely using two summation i goes from 1 to 3 j goes from 1 to 3 A i del ni dot del nj aj right. And then we are just rearranging this we are just taking A i and aj outside and integral del ni dot del nj as you know inside. So this is there is a typo here these two should not be there. So because we are replacing this square by the dot product right. So now here C ij now what we are going to do this is you know i goes from 1 to 3 j goes from 1 to 3 so this you can write in a nice 3 by 3 matrix Te right. T11 to T33 and then C ij for that elemental matrix is given by this. So remember this what we are doing we are simplified this F1 as this and then it becomes A i then this bracketed term and then aj so this bracketed term involving this integral they are calling it as C ij okay and then it will become 1 upon mu e del ni dot del nj integrated over the element area right. And also remember that we have this mu in general mu is same we can assume mu as same over the element area so mu is constant over the element. So that is why mu does not come inside the integral so mu is coming outside the integral because mu is constant over the element okay. So then you know how do we calculate this C ij just we will see but before that now having replaced C i capital C e as this whole 3 by 3 matrix this F1 can be written as elegantly like this half A e transpose C e capital C e A e where A of that element is A1, A2, A3 column vector right and then remember this is this capital C e is called as element coefficient matrix and it has information about the geometry and the material property how because you know this is del ni dot del nj because ni nj are functions of x and y and it has the information about all the coordinates of the you know all the 3 vertices and that in general it is a function of x and y and x and y being any arbitrary point within that element. So it has information about geometry and material property right. So this geometry and material property in general will determine the energy is it not energy related to that element and then we are actually using this element coefficient matrix will eventually form global coefficient matrix by combining all such element coefficient matrix. Now this C11 is simply if you put i and j equal to 1 and 1 so C11 is small C11 will be simply given by this and del n1 dot del n1 here I have written the expression for n1 as this right. So now actually del n1 will be simply given by this because this is anyway constant so derivative of this with respect to x will be simply y2 minus y3 into the corresponding x hat because del is an vector operator is it not. So and then similarly the second term will be yes again it is just we will get simply this. So now you can further see del n1 dot del n1 will be simply you know further simplified to this right. Now here this is 1 upon 2 delta 1 upon 2 delta so 1 upon 4 delta square that basically gets cancelled with 1 delta because integral ds is equal to area of that element triangle so 1 delta gets cancelled so that is why you have only 4 delta not 4 delta square right and then similarly C12 you can calculate like this now C12 now here there will be two terms because i and j so you will have here del n1 dot del n2 right and then correspondingly if you use the expressions of n2 and n3 like n1 if you have seen earlier you will get this C12 and C13 expression. Then similarly you can find expressions for all other remaining coefficients of this p matrix right. So in general then you can write pij which is this as 1 upon 4 delta ue p i capital p i pj q i plus q i q j where p i p and q i are given by this so you can verify and then the area of triangle will be given by half p2 q3 minus p3 q2. How that you can verify quickly like this we know the area of triangle as half determinant 1 x1 y1 1 x2 y2 1 x3 y3 if you actually expand this determinant you will get this right and I actually if you substitute this delta expression and substitute this p i's and q i's from this you will again get the same expression. So now this you have to remember this actually this expression is very important cij will be using this quite often later on right and also this area of triangle in terms of cij and cij going further now this we are talking of f2 f2 is these f2 and that is related to this j. Also remember when we wrote this expression this is for the whole domain when we went to FEM discretization we went over summation e is it not so that means from whole domain functional we went to element level and that is why f is now summation over all element right and then we actually integrated this element wide right. So we have already gone from whole domain thing to element level right. So now similarly this second term which is representing the so from the whole domain we are going to go to the element level. So again summation over all is and that expression right and in place of a e we are substituting it as summation i goes from 1 2 3 ni ai n1 a1 plus n2 a2 plus n3 a3 that will be a because j a ds so a is given by n1 a1 plus n2 a2 plus n3 a3 is it not. So going further now again this j can be brought outside the integral because we can assume and in most of the problems in magnetic field calculations j over individual element is constant because we have typically a conductor or winding over which you have current density which is imposed is constant in that case when j becomes constant and is independent of x and y. Going further we are just rearranging this you know expression we are interchanging the summation and the integral operator because it does not matter finally this f is what is a scalar variable is it not some scalar number is the energy whether you first integrate and then some or the otherwise it will not matter right. So here it will so we are interchanging this integral and summation operator because the result is not going to be there. So let Bi e be j integral ni ds right we are just you know calling this as Bi e. Now we will use this formula which is you know fairly complicated derivation is there for this but will not go get into derivation but I show proof of this applied to a simple case in the next slide. So now remember this ni's are functions of x and y is it not ni's are functions of x and y. So this will be in general when you have here of course here it is only ni raised to 1 and n2 and n3 are not there but in some other cases you will have in general n1, n2 and n3 raised to something and in that case the formula reduces to this after some complicated derivation. So it is L factorial n factorial upon into n factorial upon L plus n plus n plus 2 factorial into twice delta delta is the area of the element that is in this case our triangle. On the next slide we will see the you know one example for this. Now in this j here it is ni raised to 1. So here L is 1 for i equal to 1 and m and n are 0 because n2 raised to 0 is 1 so that is why they do not they are not adhering there we can say. So L is equal to 1, m is equal to 0 and n is equal to 0 that you substitute in this formula then you get this whole thing as j delta by 3 right after you know substitution and for the simplicity. So then you will get b1e, b2e and b3e all as j delta by 3 right by this and then you know j delta by 3 1, 1, 1 is the column vector. So you get this b e so what is b e? This is the element level contribution by the source what is the source? Source is j. Now j is actually distributed over the entire elemental area right is distributed is it not but that is a portion equally to the 3 nodes. So j which was actually you know distributed over the entire element area is now taken to be equal to j, j delta by 3 at the nodal vertices is it not? So it is a part of distribution. So it is like this j delta by 3, j delta by 3, j delta by 3 will be the corresponding d matrix for the element which represents the source and then f2e the energy related to the element due to the source term can be now written in matrix form as AE transpose Now this f2 has to be added to f1 to calculate the total energy. Now let us see the you know sort of a basic example of that formula. So let us take one triangular element like this I have purposely taken a right angle triangle so that we get a simplified integrated okay. This triangle is between 0 0, 1 0 and 1 1 these are the 3 vertices and nodes are 1, 2 and 3. So area of the triangle is half base into height right simply half right. Now n1 is this we know the expression for n1 right. So now you can actually substitute the values of various coordinates into this n1 expression and then you will get n1 expression as 1-0 right. So now you can also verify that at when x is equal to 0 that is at node 1, x is equal to 0, n1 will be equal to 1 which should be and at the other 2 nodes 2 and 3 x is 1 so 1-x is 0 so at these 2 other nodes n1 goes to 0 that is the basic property of the shape function is it not n1 is 1 at node number 1 0 at node numbers 2 and 3 so that also you know you can easily verify here. Now similarly then using the expressions for n2 and n3 you can very easily find out that n2 is n-y x-y and n3 is just y. Now for node 1 now we are trying to evaluate this what is our final objective here? We want to understand this you know formula and in general integral n i ds dy how do we know get it. So now let us evaluate integral over the element area n1 ds and ds is nothing but ds dy. So now you substitute n1 is 1-x dy now what we are doing is we are integrating from 0 to 1 x because x varies from 0 to 1 here and y varies from 0 to whatever is y on this line. So on that line y is equal to x is it not because the equation of this line is y equal to x. So our integration becomes 0 to 1 then 0 to y equal to x 1-x n1 ds dy. Now actually if you just evaluate this because it is just a function of this 1-x and this is 0 to x so this can be evaluated and you finally get 1 by 3. And then you can also see verification we thought del what was our you know by that formula it was delta by 3 is it not. So delta area of triangle is half so half by 3 again is 1 by 3. So we in fact verify this formula it really works for this right angle triangle if it can work for this it will work for any arbitrary triangle also. Going further similarly for node 2 you can verify it we already got the expression for n2 as x-y and n3 as y x-y and y and then if you do the integration similarly again you will get 1 by 6 and 1 by 7 right the integration result. Okay so now what we have got we already have derived you know element level coefficient matrices 3 by 3 coefficient matrices which how many of them are there 18 because there are 18 elements remember our geometry was this there was a conductor rectangular conductor at the center and this was the boundary. So now here there are 18 elements which are now here see red numbers written in red they are the element numbers right numbers written in blue 1, 2, 3, 1, 2, 3 they are all the local node number 1, 2, 3 and green are the global node number 1, 2, 3, 4, 5, 6, 7, 8, 13, 14, 15, 16 so green are the global node number. So what we have to do now from 18 3 by 3 coefficient matrices which we have derived using this formula is it not CijE so this will give you 3 by 3 matrix for each element right and all these is a function of only the coordinates of vertices and the material property is it not right. So after having calculated those we have to now form one global 16 by 16 coefficient matrix because why 16 by 16 because there are 16 global node numbers. So there are eventually 16 potential variables with respect to which we have to minimize the energy is it not varying which we have eventually minimize the energy out of this 16 you know node numbers some of the potentials at some of these node numbers are known for example on the whole boundary we are going to put A equal to 0. So those potentials will be known so eventually we will impose that boundary condition right. So at this moment there are 16 nodes so our matrix size will be eventually matrix size will be 16 by 16 right. So now what we will do is like we did in 1D example there also we had connected matrix is it not and we in fact saw the 1D code also. So here now the connectivity matrix we can form similarly here we have E1, 2, E18 as the element and 1, 2, 3 the corresponding global node numbers are given. So these are the local node numbers of each of these elements with the global node numbers are for example E1 element number 1 local node numbers will be always like 1, 2, 3 in blue global node numbers will be 1, 6, 5 right node number element number 2 1, 2, 6 element number 2 global node numbers are 1, 2, 6 right similarly you can do for and then global source matrix is B12, B16 and that also we have evaluated is it not this will be individual B1, B2, B3 for each element are evaluated using that we have to form B capital B2, capital B16 into 1. So again we have to form this using those small bi for each element right. Now how do we do that we will see now in the further slide. Now here this is you know by intuition we can do it as we did in 1D. In fact when we saw 1D code we in fact saw how do we you know combine various element coefficient matrices you know entries into global. So same thing is being done here we want to find out C C11 or this capital C is the global coefficient matrix. Capital C11 will be now this is the global coefficient global node number 1. So global node number 1 on that two triangles are incident is it not because this number global node number 1 is vertex of element number 1 as well as element number 2. So capital C11 will be given by small C11 of 1 plus small C11 of 2. Effectively what it means is potential at node number 1 contributes to both elements 1 and 2 in deciding the energy is it not. Similarly you know C22 capital C22 because on two three elements are incident who is common to 234 will have corresponding three contributions right. Global node number 6 as contribution from six entries these are the six and you guys common to six elements and you have to remember here some diagonal entries will have addition of multiple terms and some of diagonal entries will have addition of two terms. For example capital C16 will be C12 of 1 plus C13 of 2. So in case of 1D you know this was not there because 1D there was no like A there was nothing only nodes were common but there was no H being common there in one dimension case only nodes were common between two adjacent elements. But here now not only are nodes common but edges are also common and any edge will be common to only two elements the capital C13 will be 0 because there is no direct connection between nodes 1 and 3. So capital C13 will be equal to 0 right. Similarly now the right hand side matrix those matrix B1 will be B1 B1 E1 plus B1 E2. Now these small Bi's we have already calculated J delta by 3 is it not. There those we have already calculated. Similarly B2 will have three contributions and so on. That means at the end of this step what we have got we have got we would have got one global coefficient matrix of 16 by 16 and one global source matrix which is 16 by 1 right by combining all this you know element coefficient matrices and element source matrices. So now your F the total energy is F1 minus F2 and that is now you can write it as half A transpose C A minus A transpose B where now A and there it was AE it was AE is it not. So it was only 3 by 1 now it this is A transpose will be 1 by 16 right remember always all these matrices are the column vectors. So A transpose means it becomes a row vector. So A transpose will be 1 into 16 right C is 16 by 16 A is 16 by 1 is it not. So the total 1 16 16 and 16 1 so the whole this product will be just 1 by 1 right similarly this A transpose we just saw it is 1 by 16 A transpose is 1 by 16 B is 16 by 1. So again this A transpose B will be 1 by 1 and it should be because this F is simply energy the scalar just one number final is it not. So it should be 1 by 1 okay then you know the fourth step of this whole FEM procedure we go towards energy minimization and then impose the boundary condition. So daba F by daba A equal to 0 now this A is the 16 by 1 column vector that means this amounts to 16 equations is it not with respect to every A i right and that what is written here those 16 equations. So this you can then eventually get when you actually differentiate all those Fs capital Fs by individual A i's you will get the equation of this form CA minus B equal to 0 and then CA is equal to B but here now B has contribution only due to source till now B has contribution only due to the source. Now there will be additional contribution due to boundary conditions we already you know discussed this point that right hand side matrix B will have contribution from two things one is the source or there is a boundary condition right. So now we will impose the boundary conditions and modify this system of equations right. So now boundary conditions are all those you know outside outermost restangles all those potential are equal to 0 these also I mentioned earlier this is just we are taking it as 0 as a you could take instead of 0 you could take something else also but the answer will not change because flux is integral 8.dl and for 2d case it will be simply difference of potentials at two points. So only difference matter if you scale this boundary condition to some other value the difference is going to remain same between any two potentials is it not and that is why the flux will remain same between the between any two points right. So then B matrix can be modified by transferring terms containing known voltages and the system of equation can get reduced to 4 by 4 right. But however you know it is not amenable for coding and in previous code when it was explained Sairam has already explained you both these cases in this case how will we get 4 by 4 and how will we how we can continue with by 16 by 16 matrix and go further so both things are explained with respect to that one decode is it not. So here we will continue to solve 16 by 16 system because that is easy for on the point of a coding rather than converting into 4 by 4 because if you want to convert it to 4 by 4 you have to do some row and row operation is it not you have to you know add or subtract some row from the other and then simplify it so it involves lot of operations extra operations so rather than doing that operate the whole matrix because that 16 by 16 matrix although it looks bigger matrix many of those rows are simply you know statement like something is 1 or something is 0 as you know will be obvious here. Now for example here we want to impose condition that A5 which is on the outermost boundary is 0 because we want to impose the condition A5 is equal to 0 how will we do that the fifth row of this equation CA is equal to B what you do you make all the half diagonal element entries equal to 0 and only this A capital C55 make it 1 and then corresponding fifth entry here in the B matrix also you make it 0 so when you expand this you will get A5 equal to 0 effectively you you have imposed boundary condition A5 equal to 0 right so now effectively what will happen on these 12 there are how many node boundary nodes there are 12 nodes is it not so 4 plus 4 8 plus 4 12 so 12 rows will get modified like this so in 12 rows you will make half diagonal entries as 0 and the corresponding diagonal entry as 1 and on the right hand side corresponding capital B entry you will make it 0 right and then once you do that and then you invert the matrix is it not so finally you have to invert this you have to take C inverse sorry yeah C inverse you have to take and then you will get the A as C inverse B and then we will get this solution now here this solution is shown here with course mesh now you can see here the field lines are not smooth one of the very first lectures I mentioned that field distribution while just looking at the field distribution you can know whether the you are you are meshing is good enough or not so here field contours are not smooth or circular that means there is scope for improvement in the mesh and now you make the mesh very fine you will get you can see now the field contours are quite smooth but still you can you see here since we have got this boundary quite close can you see here this contour sorry actually suddenly you know it becomes flat here although it is smooth ideally if it was a isolated conductor all these these flux contours particularly this last one near the boundary they should approach more or less circle is it not since our boundary is quite close for it to become an isolated conductor because for an isolated conductor the boundaries should be far off so that the boundary conditions do not affect the the field distribution so the fact that we have not put that boundary far off now you can see here suddenly you know here you are sort of forcing the potential to go to zero otherwise it should have been something like this you know by the circular contour so that also you know you can check while looking at the field distribution and in fact later on we will see by using boundary which is closer we are not really getting the correct answer because when we calculate the inductance of this bar we will find in this case the error is high when the boundary is closer whereas when we take the boundary far off the inductance value that we calculate is more correct as compared to the actual value okay so with this end we will end this 18 lecture I hope you have now understood finite element procedure because we have seen 1D as well as now 2D including you know the procedure now in the next lecture what we will do we will like we saw for 1D we will see 2D code but not in detail because we have already seen 1D code okay thank you