In this segment, we'll start looking at the main function in the template file and that's assemble system. So, let's move over and look at the code for a second. The assemble system is only called once, of course. So here we zero out K and F. I've defined several constants and objects up here. So I don't think I've used this before. I have const unsigned int. So what that constant does is, it allows me to define this variable only once, so it's actually not even a variable. It can't vary, it's a constant. So this dofs_per_elem is degrees of freedom per element. And of course that changes with your basis function order. Same as number of notes per element in our case. Then I declare a full matrix. Again, this is another DL2 object, and it stores doubles and as I'm initializing it here, Klocal is dofs_per_elem by dofs_per_elem. Okay, it just defines the size of the matrix again. Flocal is also a DL2 vector. With length dofs_per_elem or nodes per element. Now here I'm creating this standard vector of unsigned ints called local_dof_indices. It also has a length of dofs_per_elem. And this vector will actually update for each element. It's different in the values,and will be different for each element. But what it does is it relates our our local node numbering to the global node numbering. So let's look at that for a minute. Again, I'll create another sample domain. I'll do let's, I guess in this case three elements. For quadratic. So we have one mid side node. And let me do the global node numbering on top. And then on bottom I'll do the local or the element node numbering. Again, both of these are in DL2's syntax. All right, so we have 0, 1, 2, 3, 4, 5, and 6. And now I'll break this up into element 1, it's called it e1, e2, And e3. Okay, so those are our three elements. So locally, these nodes would be labeled 0, 1, 2. Same for all of them, 0, 1, 2, and 0, 1, 2. So local_dof_indices will look different for each element. So local_dof_indices. For element 1, this is what it will look like, local node 0 corresponds to global node 0. Local node 1 corresponds to global node 1. Local node 2 corresponds to global node 2. Okay, so they all match up there, and let me specify here that these are the global nodes or local nodes. All right let me change that. So on the top this is the index, is the local node, and the actual value stored is the global. So the index is local. And the value is the global node. So for e2, again we have 0, 1, 2. But globally now we look at element 2. Local node 0 corresponds to global node 1. Local node 1 corresponds to global node 3. Local node 2 corresponds to global node 4. Now, of course, DO2 is handling all this for us. It will populate this local_ dof_indices vector. But I just wanted you to understand exactly what's what's going on here. So, for element 3. Now we look at local node 0. Here corresponds to global node 3. Local node 1 is global node 5. Local node 2 is global node 6. So if we, let's say we're working with element 2. We've initialized local_dof_indices for element 2 and so now using the C++ syntax, if I put in 1 as my index it would give me 3. If I put in a local node of 0, it would return 1, a global node number of 1. So this vector local_dof_indices will relate our local node number to the global node number. The index to the vector is the local node number. The value stored is the global node number. And we'll need that especially when we're going to assembly. So now let's go and let's see I also have these variables stored here they're all doubles, h_e is the element length, x we'll use to find the value of x at quadrant points. Which we'll need to define our forcing function. And then f is the value of that forcing and function. Now we have a loop, a loop over elements. This is where, remember in our review of C++, in our introduction to C++, we talked about iterators and they're similar to pointers in some ways. And that we can use these iterators and four loops. In this loop over elements, we're actually using an iterator. You can see here, it's an active cell iterator. By the way, in the DL2 name, nomenclature, they call elements cells. So when you're seeing a name that DL2s set up, it usually says cell for the element. So, here, we're actually looping over elements, they're called cells. And now all of this I could have put, this line elem = dof_handler.begin_active() I could have put here in this declaration portion of the for loop. But since it's so long I can declare it ahead of time and that's okay to have an empty spot there in the for loop, declaration. So I've declared two objects here, I have elem, which is my iterator and it's And I'm setting it to the beginning element. I also have this endc, which is the last element. Or really it's a marker after the last element, okay. And so my for loop, I'm already starting at, element equals the first iterator. For the iterator is equal to the first element. And my condition is, as long as I'm not at that marker past the last element, then I'll go through the for loop. And again, this ++elem, it just increments to the next element. Again, even though elem is not an integer, it's an iterator. The ++ still iterates forward, okay? So now that we're within a particular element, we have to update the local_dof_indices vector, okay? So we do it this way, elem->get_dof_indices. Remember this arrow is the way we, for either a pointer or in this case an iterator, this is how we access a function. Okay, so since elem is an iterator, instead of doing elem.get_dof_indices, we have to use this arrow essentially, okay. And so that will populate the correct local_dof_indices with the correct, global, node numbers, okay? All right, now, I'm going to find for you the element length. As I told you before, the way we've defined our mesh, each element will have the same element length. So, you could actually have just taken the length of the whole domain and divided by the number of elements to get element length. This is just a little bit more general in case you ever set up your mesh to have elements of different lengths. What we're doing here is I'm saying h_e = nodeLocation at local node i minus the nodeLocation at local node 0. Okay. Remember, local node, or sorry, local node 1. So local node 1 again is always on the right. You can look back here on our schematic here. The local node 1 is always on the right of the node. Local node 0 is always on the left. If we take the x coordinates of those two and subtract them then that gives us the element length for whatever element we happen to be in. Again I use local_dof_indices though Because I may not be in element 1, I may be in element 2. And I really want the, and remember, nodeLocation, the vector nodeLocation, is indexed according to the global index. And so I would want, in this case, if I wanted the x values at local node 1 and 0, I would need the global node numbers 3 and 1 in order to access the correct x values within my nodeLocation vector. So that local_dof_indices will come up again as well as nodeLocation. But now we have h_e which is our element length. Now, first let's create our local force vector, Flocal. Now there are really two parts to this. There is the portion due to the forcing function, or the body force, and then if you have a Neumann, condition, a non-zero Neumann condition like we do on problem two, then you'll have to add in the effect of the traction at that boundary, okay? So first we will do the portion of Flocal due to the body force. Okay, so here I've set Flocal = 0. Again, we zero it out for each new element. And now in looping over A, so each component of Flocal, so I'm at Flocal[A]. And now I'm looping over my quadrature rule. Okay, I'm looping over each quadrature point, because we will need, remember in Flocal and in Klocal we're performing integrals. All right, so let me actually write out this general form for this portion for Flocal, and again this is in the lectures. Let me write this out. So Flocal[A] is equal to the area the cross-sectional area times the element length over two. Remember, we use that element length over two to convert from bi unit domain to the real domain. We're integrating from negative one to one. Basis function a corresponds to our index here evaluated at node xi, times our forcing function or our body force, f. f we're given as a function of x, but since we're integrating over xi, we need to evaluate everything in terms of xi. Okay, and then we integrate it over our d xi domain. Okay. Now since we're using quadrature, let me write this out using quadrature points. Ahe/2 is still on the outside. However, now, we are going to loop over our quadrature points. And so, I'll just do q equals 0 as long as we're less than quad rule. Okay, and now we're going to evaluate the integrand at each quadrature point. And remember you get the quadrature point from that quadPoints vector. I'm going to distinguish it here just by saying it's xi q. So xi q is the quadrature point. Okay? And f is a function of x, but again we evaluate, you can get x, a value of x for that particular quadrature point. Now it's no longer an integral. It's a summation. And so we also multiply it by the weight. So that's short for the quadrature weight at that particular quadrature point. And again, we also have a vector of quadWeights. Now, this term should be easy enough to evaluate because we already have our basic function that has as an input the node and xi. Now here, you're given a forcing function or a body force. How do we do this term here? We have a value of xi. We need a value of x. So this is also something from the lecture, but just as a reminder, we can interpolate. Over the values of x at our local nodes to find the values of x at xi. Let me specify here that this is the value of x at xi q for example. Okay, so what will we do? We'll start at A = 0, and as long as A is less than the number of nodes, the element, we will take x, A times our basis function evaluated at xi q. And what do I mean by x a? This is the value of x at node A. In other words, which we can get from nodeLocations, okay? NodeLocations, but remember, A is a local element, or a local node number, nodeLocations requires a global node number. And so, to convert we use this local_dof_indices of node A. I've actually written all of this up for you in the code, but I wanted to explain it to you in this one case. Because you will be performing other integrals in this homework, and in other homework assignments, so I want you to understand what's going on here, okay? So, going back to the code you can see, again, we're looping over B=0, B < dofs_per_elem, in other words, less than number of nodes in the element. And x +=, so we're doing this summation, nodeLocation(local_dof_indices[B]) * basis_function(B, quad_points(q). Okay, again I'm just using the next B instead of A but it's the same thing right? Okay, so now you have the value of x here, this is where you have to start doing it yourself. So now that you have x, you'll need to evaluate f, evaluate N A, and essentially just create this integrand, well what was an integrand and is now this inner portion of this summation. And you will use that to define Flocal at A, okay? So you'll need to do that. Now that takes care of the body force. However, in problem two, or part two of the this homework, you also have a Neumann condition. Of course that only applies at the right hand node. And so what I've done here, is if were at problem two, then it's checking if your node location, if the node that your at is on the right hand side, then that's where you'll have to add in the traction times the area which is H, okay? So you're really just adding an h to the appropriate component of Flocal. And again you can review your lectures if that process is a little bit fuzzy. So that will take care of your Flocal. So we've taken care of the body force and when applicable we've taken care of the Neumann boundary condition. Now the next part is to define Klocal. And let me write out the general form for Klocal as well. Now Klocal is of course a matrix so let's look at component A, B of Klocal. We have 2EA, E is the Young's modulus, A is the cross-sectional area over h e where h e is the element length again. And here we've already converted to an integral over the bi-unit domain. That's why you have the h e in there. Now we're dealing with gradients. So this is the partial, I guess the derivative of basis function A with respect to the xi, times the derivative of basis function B with respect to the xi, okay? And so again, these two terms you'll get from your basis gradient functions that we've defined previously, all right? And again, here you'll convert this integral into a summation over quadrature points, multiplied by the appropriate quadrature weight. Okay, so I've set up the for loops for you. We have a loop over A and a loop over B, and a loop over your quadrature points. Now you simply have to write out what is Klocal A B equal to? And of course, since it's a summation, you'll be doing a += just like you did with Flocal and like we did when evaluating the value of x at the quadrature point, right? Okay, so once you've filled in those lines, you'll have your Flocal and you'll have your Klocal, okay? Flocal will already have any Neumann boundary conditions that are involved in the problem. Now you need to go to assemble system, you have to actually assemble them into the global matrix in the global force vector. Now again this was done in the lectures. You can go back and review that there. The real trick of it here is to remember to use local_dof_indices and it makes everything very simple. So for example if for F We have the global nodenumber. So essentially you'll, for each element, you'll be looping over each node in each element, okay? And you'll put that value of Flocal into the corresponding value, that should be a +=, because we don't want to delete any values that were already in F, right? We just want to add on to them. And so you add on the Flocal for the local node number. Okay and the local node number will of course be going from zero up to the number of nodes in the element. And then how do you go from the local node number to the global node number?. Again its just local_dof_indices, okay? All right so, it's actually pretty straight forward. You will do that with F and you'll also do that with K. Of course, K has two indices, so we have two loops here. We have a loop over A and a loop over B. Now, the one thing that you need to be aware of, really, is that since K is a sparse matrix, we can't just do K [A global] just as a short hand and [B global] You can't just do += like this. And of course A and B here are the local. All right, you can't do it this way. We instead do K.add. And it's because it's a sparse matrix and we're not dealing with all of the, not all of the indices are actually available to us, okay. So we do this K.add, you have your A global, B global. So these are the indices, these correspond here. And then this term is the actual value that you're adding. So it would be this here. Okay? So those are equivalent as far as the thinking goes behind it. This is, of course, the actual writing that you'll do. Okay, now this assembly process is really just two for loops, just a few lines. It's good for you to understand what is going on here because you'll be doing this exact same process in all of the remaining homework assignments. It's actually the exact same lines of code. Okay. And so it's good for you to understand what's going on here in this first homework assignment. Now the last function that's called in the symbol system is this apply boundary values. Now we're not going to go through all the mechanics that we learned in the lectures of applying the Dirichlet boundary conditions. You remember in the class we had K, D equals F, right? If we happen to know that boundary conditions, if we had Dirichlet boundary conditions, and say these, excuse me. If we had Dirichlet boundary conditions in these two, on these two nodes, then what would we do? We would actually remove the corresponding rows and columns. Okay, we'd also remove the corresponding entries in F, and that would give us a smaller matrix vector system, right? I think we set these as bar, and then these were actually K, D, and F, or something like that, right? The way deal two does it a little bit differently. Instead of actually deleting these rows, it just modifies the entries. And KM modifies the entries in F to give you the correct values of D, or the correct displacement values, okay? So you can actually look into that a little bit more if you'd like to, but the main point is that the deal two will be taking care of the application of Dirichlet boundary conditions, and it actually does it without modifying the size of K, D, or F. Okay, so that is what's happening here. You can see we're inputting boundary values, which is the map that we created earlier that holds the values for the Dirichlet boundary values, K, D, and F. This false has to do with whether or not deal two is modifying entries in the row only, or if it also modifies entries in the column as well. False means it's only modifying entries in the row. Again, that's not something you'd need to modify, but it's helpful to see what's happening there. Okay? All right, let's stop that segment here. And in the next segments, we'll look at solve and output results, which are really short functions, but then, we'll look mostly at calculating the LT norm of the error.
