 Hey guys, welcome back skidstone series episode 26 I'm gonna go nice and slow and clearly in this video topic today is the finite element method We already covered this in great detail in an hour and a half long deep dive I'll put a link to that in the description. This video will cover more of the assembly language implementation for that topic and we're gonna focus on 3d frame elements, so We're finally getting into the actual nitty-gritty engineering topics, which is a good thing First off an example of what this even is now This is a ladder and this ladder is actually the same example that I gave in the previous video just in assembly this time not Matlab and so in green you can see the Undeformed ladder geometry that it's you know pressed up against the wall There's some angle here against the wall and then in red you have the deformed Configuration, so there's someone standing on the top wrong of the ladder not very safe in my opinion Applying a load in that direction and of course, there's also boundary conditions applied on the system So this bottom two nodes where the ladder hits the ground they're being fixed in space both Positionally and rotationally and then at the top of the ladder It's basically constrained that it can't move out of that Plane that it's in it's stuck in in this case the z-direction Basically, that's the ladder is up against the wall and it's not going to move through the wall or off the wall Couple other things about this The lines are elements the circles and x's are nodes talk about that what that is in this video and Lastly obviously the deformation is a function of the forces the moments the boundary conditions the materials the geometry all that but also one other thing that there's a function here is the scale factor and so in reality this deformation mode that you see here is very exaggerated I Multiplied the deformation by a factor of 500 just so you can see So with that out of the way, I'll get into more of the kind of the theory now I'm going to go into too much detail here because we already made the other video, but I'll cover the basics in this video so What is the finite element method? Basically, it's a way that engineers analyze structures among other things Idea is you divide a structure into its simple constituent pieces We call them elements with nodes, but you can imagine you know a bridge being decomposed into its individual trust elements or any kind of house framing or railing or fence or anything like that as long as you can and I give you examples of frames but Obviously any kind of system a table Any kind of body an aircraft a car can be modeled in the same way But in this video we're talking about frame elements, which I'll cover what that is in a second Anyway, the idea is that you divide a structure into pieces So for that ladder each rung if you saw was two elements the rails had a bunch of elements along them, etc Then what you do is you describe how each of those pieces by itself theoretically behaves as a function of its material of its geometries and of its Degrees of freedom so two end points if they move in this way What forces are required to make that happen and vice versa? How can the forces affect those degrees of freedom if I twist it if I bend it if I pull it? How do those forces affect degrees of freedom? That's called those equations are summed up into a matrix and we call that an elemental stiffness matrix Then once you have that matrix set up for every element in the body You can then Based off which elements are touching which other elements you can assemble a global Synthesis matrix, so let's say you had ten elements you put them all in one each of those Individual elemental synthesis matrices you assemble into one Global synthesis matrix in the global system and to do that you take the element connectivity So this is element one this is element two and this is node zero node one and node two Both elements one and two contain node one So the Degrees of freedom have to match so if I apply a force at node zero and a force at node two Those disciplines can move around in space any which way. However at node one it has to match Right you have to have those degrees of freedom from element two and element one match at node one because they're touching You can do that. They can't fly off apart from each other And then lastly once you've characterized the entire system of equations for the entire structure You can pull in all the forces that you're applying as well as all the boundary conditions that you're applying Get that entire Linear system set up and then solve it in our case. We are going to use I Think pivoted l u decomposition, but there are many other ways to solve these systems as well Okay, what's their frame element? Also called the beam element It's basically a long skinny object that you can describe the displacements of the two endpoints and kind of capture the whole system movement as of, you know As a result. So basically here you can see I have this long beam in pink here with node a and node b on it at each node I have a According system set up where x is along the axis of the beam and then Y is up and z is to the side obviously that's all you know Kind of random you can pick whatever you want for that and You can see that not only are there systems at each end. There's also displacements. That's the u There's rotations. That's the theta as well as forces f and moments m at each of those points So you could imagine pushing node a in the x direction Pushing it in the z direction pushing in the y direction Twisting it around the x direction, etc. So many different things you can do both at node a and at node b So there's six degrees of freedom per node 12 degrees of freedom for the entire frame element in 3d Now if it was planar Different story completely, but it's a this is a full 3d frame element in this video and then one last thing you might notice that the Moments have kind of random arbitrary signs so a negative some positive That's all picked arbitrarily to make the shearing Work out but as long as you're consistent with the moments for all your elements in the model You're okay to pick whatever sign that you want as long as you have one negative one positive obviously so positive moment y direction at node a negative at node b Positive at node b for x negative a etc for z positive negative Okay, no big deal What's next so The Elemental forces those things in red there Need to be related to the elemental displacements. Those are in green here So those 12 degrees of freedom in green are being related to those 12 Applied forces by this 12 by 12 pink elemental stiffness matrix k And it looks kind of big and i'm not going to get into the details of each individual Entry here, but this is basically populated with various geometric and material quantities So a is area l is length e is a is a stiffness property as well as g um j and i are inertia terms and That's pretty much the gist of it Um, but basically if you want to get details about how those terms come you can watch the other video on this topic that I made Um, and how does this work? So if you think about it think about force in the x direction at node a that is this guy here so obviously a force In that direction a compressive force an axial force at node a That is related to you think about rows times columns It's related related to basically this degree of freedom here and this degree of freedom here By these coefficients here And they're equal and opposite in sign but basically what you can think of is that a force in the x direction at node a is going to Basically slide this element in the x local x direction and that was only going to affect this value and This value of course the reason why the sign is negative for the second half is because A compressive force in this member Is going to actually be this way at node b and that's going to relate in a negative displacement in x That's why the sign on this coefficient is negative while this one is positive But who cares about that watch the video if you're curious But yeah, anything like that. So for example, um, how about this last one here the moment at node b about the z axis Takes these four coefficients Which relate let's see the second one So the y's as well as the data's It's this one and this one Which makes sense, right? Let's look at the picture. So that's the moment in the z direction at node b That's this one here. So moment like this And like this you can see how that's only going to affect the Angles of the two nodes About z as well as potentially it could push up node a and push up node b That's a displacement in the y direction at both of those Would you if we go back That makes sense Because we've just marked those trees of freedom. So there is a coefficient that relates that force to that particular degree of freedom moving how much well look at the value of this geometric slash material quantity here So it's pretty cool. If you have a system that's just one element, this can perfectly describe how the element Reacts to a applied force And then usually of course You know the applied forces right for our ladder example We knew that the force was you know 200 pounds at the top of the ladder You know what's being applied to the system. You don't know the internal forces But you don't care about that what you actually care about usually most of the time is the strains In the system or stresses, which are just basically the strains times something What is strain strain is how much each thing is moving how much each thing is stretching? That's basically this quantity here. How much we're moving in each direction twisting, etc So you have this kind of system where you have the forces, which you kind of know most of the time So they're the applied forces. This is a geometric slash material quantity. So you can get that by measuring Um You know if it's wood these are the values if it's steel these are the values If it's this thickness and this width you can compute the inertia terms, etc Now this you don't usually know and that's why before I mentioned our objective is to Solve the system of equations while you're solving for this unknown U vector here Now you might know some of them in the case that I gave before the ladder We knew that the degrees of freedom were fixed at the bottom So that node Let's call it node a all of these values here u You know theta u theta u theta We're all fixed to be zero for both of those two points at the bottom of the ladder And then also at the top we constrained the z direction. So that was u z and u z or whatever Where is it? Yeah, well for both nodes. So it was u z for both nodes Um, so you can constrain things and how does that work? Well? I won't get into it in this video check the other video, but you can apply constraints as well to the system Um, and you have to otherwise you get a you know an unsolvable system at you know singular Off the charts condition number, etc Okay, so how do you assemble that now you have each element right the idea was you created this 12 by 12 matrix for every element every wrong in that ladder every You know frame every stud in your house you came up with this 12 by 12 matrix What you do then is you have to assemble those into one big matrix because again where those elements touch you have to constrain the degrees of freedom to match so In this case for this bridge you can see here Well, let's go through the you know process so This entire bridge was broken into multiple frame elements And each element has its own 12 by 12 pink matrix um And each of those matrices needs to be transformed into a global system, but then once you do that you can uh assemble one massive system where Multiple elements that contain the same degree of freedom Have their individual stiffness contributions for that degree of freedom summed up And so how does that work? Well look at node 13 for example of this bridge um It has its own six degrees of freedom six ways in which it can move it can move in the x and the y and the z as well as about the x the y and the z directions and Those degrees of freedom are affected by the stiffness is of element six element seven and element eight remember if element Six is made of rubber and seven is made of rubber and eight is made of steel It didn't matter that the two elements six and seven were very soft. The steel element Was what made the system very stiff And so that's why you have to be able to deform the rubber piece And the steel piece both with the same degrees of freedom and the same forces So you're going to be adding up the contributions of each elemental stiffness entry at each degree of freedom Of course one key detail is the Orientation of each element needs to be picked needs to be the same So If you go back to this we've defined the elemental coordinate system to kind of be x along The the element and then y and z are kind of random When we get to the full up frame system that we have like this bridge, we have to pick one common System one common Coordinate frame to you So maybe call this one x call this one y call this one z or whatever you want to pick It has to match for all the elements So you have to basically transform Each of the local systems into the global system while you're assembling them Not too hard, but definitely you have to worry about that as well so How does this work in our implementation? Well, it's pretty simple. I have this data structure here which encodes pretty much all the relevant information for each 3d frame system And here are the entries. So the first one is the number of nodes So I have a bike frame here as an example You can see I have marked each node of this bike frame node zero one two Three four five and six. So there's a total of seven nodes dq seven Then we have 10 elements. Let's count them You know, there's one here one here one here And there's four green ones this these two salmon ones and then a blue one over here That total is 10 and then number of element types So this is a way to be more efficient so I can say well All these pink elements are more or less the same geometry and same materials so they The length may change But they have the same cross-sectional area the same inertia's the same You know material properties and so all the pink elements are of one type call that type zero The salmon ones here are maybe a little bit thinner Maybe different material or whatever Probably not, but that's its own element. So there's one pink element set one salmon element set One green element set and then one of those blue element sets back there. So there's one two three four element types Then what you do is and this is it's going to have to do manually is You set up an array for the nodes so Every single node node zero one two all seven of them you have to describe its position in space x y and z That's how this node array is set up. And so this is a pointer to that array It's somewhere else in memory. You put the address here Then the element array so in this case You can see here each row of this array is a Quad word of Three things so each row is three quad words. So the first node of the element The next node of the element as well as the element type So you can see here, let's take element zero for example that contains node zero Node one and let's say its element type is zero The pink ones are all element type zero. Let's say so then this element array contains basically The first row of it is just zero one zero And the next one next row would be for element one. It would say node one Node two element type zero and the third and then element two So the third entry here would say okay node two node three element type zero How about element four? Well, it would say node three Node four. Oh, sorry. That's element four node element three is up here It would say node zero node four element type one Or whatever you picked And again the order can change so it could have been zero four could have been four zero It doesn't make a difference as long as you are consistent in how you Make things match up Lastly is the element type matrix now This is basically just a list of all the geometric and material properties for each element type So the pink element we call that type zero. That would be the first row. I have to write down the sickness terms so the E the g the area the inertial terms as well as this v thing What is v? I won't give you the too much detail Basically, it's a a vector that kind of is used to orient the local elements frame in the global system so it's kind of a A vector that we can hit the cross product of the Beam itself with to get the z direction. That's kind of what that's it up. So we'll put in A generic Vector that we can use to orient our elements if you're curious to how that works check the other video Then that's pretty much it for manual stuff For the most part then this function that we implemented here For this video is called the assemble frame elements function Basically what that does is it takes this entire frame data structure and it populates a full up Once you allocate the memory this can populate the Data for this whole matrix this full not not 12 by 12 But the entire system it can be you know a million by million This function will populate the entire system Based off the input data. So the material properties geometry as well as the node positions the element array etc And then of course you're trying to solve this f equals k u function f is the applied forces k is the System that we just made with this function and u is the unknowns that we want to solve for So f has to come from something Again, that's the applied forces and I guess you could say they're always going to be zero But that might not be an interesting problem to solve So in our case the forces for that ladder were you know 200 pounds at the top And you know you could think of other forces to apply as well Let's say you know someone was holding the ladder against the wall another force there. This is basically a Encoding of this so the element well The forces at each degree of freedom in the global system. So you have to make that yourself as well So once you've got this set up, what's the overall process? So first was to define the geometry get all the nodes positions defined in that node array Define the elements themselves. So which nodes are in which elements What the element types may be as well as to find the geometry and the materials That's all these items here Number of nodes elements element types As well as node array element array element type array Okay, then what you do is you use that function I talked about to create up the full you know massive to this matrix And then you apply the boundary conditions. So you may impose Zeros here. You may apply a displacement here You may apply a force here Sure, whatever you want to do you can do then you solve a system Whatever way you want. We're going to use l u decomposition, but you can do whatever you'd like And lastly you post process. So if you want to get the strains you want to get the stresses You want to see how it looks you can do that And for more details again, check the other video So now into the code Let's see how this works So there were two examples for this There was a cantilever beam example and the ladder example that you just saw. So in the cantilever beam example Let's run it first So it's very simple. It's a beam with two elements And three nodes and you can see here that it looks it's constrained on the left So zero everything at that position and then we're applying a downward force At the right hand side So very simple. How does this get set up? So let's see In this example, we have to have a heap obviously because we're going to be Rendering to the screen again that requires buffers to render to the screen As well as a few other things require heap as well So the include of importance here is this well, I guess two of them one is the Actual l u solve Function in this case the pivoted version of that which we use to solve the system of equations and then we have the actual Function of importance for this particular video, which was the one that Assembles the entire in this matrix So again, we have a cursor that we're drawing to the screen as before All the rendering is as it was in previous videos. I'll ignore that just because it's kind of repetitive at this point so Let's take a look at the actual um data structure So here is that 3d frame system If you recall there are indeed three nodes on that beam You know left middle right. There were two elements right left and right and there was one element type Why I defined it and then I have a pointer here for all the node arrays, etc And let's see I have those defined here. So here's the node array. So you can see that the left hand point is at 000 The middle point was at 0.500 the right hand point was at 000 and then the element array down here Shows that element 0 contains nodes 0 and 1 an element 1 Contains nodes 1 and 2 and they're both type 0 And what is type 0? Basically, it's an element that has these two material properties An area in the cross section of 1 and then these three initial terms as well as a generic orientation matrix Or should say vector for how we're going to position that local element frame in the global frame And then I have some space in this binary set aside for the status matrix the forcing matrix and the unknown matrix and here you can see in the forcing matrix It's 18 entries long because there's six degrees of freedom times three nodes and in that vector Looks like the 14th element down Is set to negative 50,000. So i'm applying a massive negative load at What's probably the z direction or the y direction? At the right hand degree of freedom And so how did this work? Well, basically you can see here We have that data structure defined in memory And then these two lines here I put that address of that frame element structure into rdi I call this function that will populate the global synthesis matrix for that frame system And then I go through and I apply the boundary conditions. So here you can see looks like I have set for the first node the six degrees of freedom that can Control its position and rotation are also the zeros And so the way that works is basically i'm looking into the details But you set basically identity matrix for those degrees of freedom and then you zero out all the rows and columns in those rows and columns And then I also Yeah, put the put the one on the diagonal there you can see it put 1.0 in the six rows and columns and Then we get into actually solving. So now we have the entire system set up We have the nodes positions the elements positioned All the material defined the matrix traded the forces applied The displacement is applied Now we solve and so here we're using that plu solve. So the pivoting l u decomposition method And so we have to make a pivoting vector to kind of use to Order our rows and columns if our algorithm requires it. So we've allocated memory for that solve the system And then this basically puts the solution in memory Where we've given this space. So we have this unknown matrix 18 quad words long That's attached to our frame element structure here at the bottom offset of like 64 bytes from the top So our answers are now going to be in that slot essentially So that's what this does here and then lastly Looks like I free some memory for some reason. I don't know why I did that and then I now add Well, I create all the rendering stuff. So I create the green undeformed beam and the red deformed beam as a function of those solved degrees of freedom So yeah, that's how this works. Oh, and then I should say at the very end It has the same old rendering loop as usual so perspective Structure set up and in terms of what's being rendered to the screen. We have a couple of different geometries. We have basically the line segments as well as the points Both for the undeformed and the deformed configuration And so that's total of four things being drawn if you look down here And so just to get to show you The undeformed is in green deformed some red and yes, it works as intended So that's that's a good thing I guess um example B so Most of the same thing Again, I run just giving example how it looks So again the latter example a bunch of nodes It's basically the same thing as before just with more elements, etc, you know So how does this work? Well, it's pretty much the exact same thing with one key difference And that is that the it's it's too hard for me to draw basically to Come up with the node positions manually and so I've parametrized it And so you kind of can see down here We have some parameters that control the latter angle for example as well as the cross sectional area of The rungs and the side rails and the lengths of things and so you could change these values Maybe we should let's make this value a 20 Make them longer. I don't know Make ladder, you know, how about some steeper ladder. Let's say 10 degrees. I don't know Now if we run it You can see it's definitely taller Um, and it's More steep of a ladder, but again, it's the same kind of information mode So the difference here is pretty much one thing and that it's not it's that I can't Define the ladders points and geometry Very easily manually and so I have a function does it does it for me It's called generate ladder system and you can see here It basically creates that entire 3d frame structure itself Based off the input. So if you pass in number of rungs number of No side rail elements between rungs The height the everything that you want to do for that ladder you can kind of parametrize it and call this function So yeah, that's pretty cool. Besides that it's pretty much the exact same thing as the previous one So again Uh You can do pretty cool stuff with this And this is just frame elements. Again, there's shell elements. There's solid elements There's many many more element types out there. These are just two elements per What two nodes per element you could have Seven For elements you could have non linear elements. You can do other things that you'd like as well And again, this is linear analysis. You can do non linear analysis You can do geometric non linearity as well as material linearity And so much more Also, you can couple things together. Let's say you want to let's say this ladder at the top It had a plate so you it would connect These one two three four five six seven elements together with a plate Well, then you could add that shell or that plate element to this model With its own degrees of freedom as long as it's, you know, set up properly It can intermingle with this implementation as well. So Anyway, it's very easy to do this. Um, took me some time to go through it, but it's not that hard And if you look at the file size That entire ladder fea that we just did Define the ladder, you know, the entire rendering pipeline that handles text and 3d solids and stuff all that together is 25 kilobytes so Not a lot and by the way, a lot of that kilobytes is just random Garbage for I can even show you what kind of garbage we're talking about. Let's go Let's go to lib engineer fem And here's that function that I talked about. This is the actual black box function of today's video I'll show you like look at this at the bottom of this program I have defined in the binary in the flat memory model like A thousand bytes for the elemental tip of matrix before it's transformed I have a thousand more bytes just here in the in the file in the binary for The you know The transformation, you know and the transpose and all this different stuff like I've defined and I have a working slot here This could all be on the heap. So between these five things I could cut out and save over five kilobytes But I'm not gonna but I you know, I could so There are ways to make this better more efficient, but it works the way it is and I'll talk more about this in the next video, but With that out of the way, uh, thanks for watching. It's pretty cool stuff. Hope you enjoyed it. See you in the next video