In this segment we'll move on to looking at the assemble system function for the homework four template. So let's look at that code now. You'll notice first off we have two new objects. We have this fe_values and fe_face_values. These objects are DO2 objects that hold information about the basis functions, the basis function gradients, information about the quadrature points, the Jacobians, and all this. But we no longer have to calculate the Jacobians. We don't have to write out the basis functions themselves. We don't have to write out the quantiger points It does all that for us, all right? So if you'll look, the first two are just fe and quadrature formula so that it knows what basis function order we're using and it knows what quadrature rule we're using. The third input value is actually a series of input values, of flags that tell fe values what information we're going to be using and that's the information that it'll be updating. So it makes a little bit faster or saves on memory that we don't need to update information that we won't actually be using. All right, so for fe values, which is what we use for volume integrals. We're going to be updating the values, which are the values of the basis function. Update gradients, that's again the basis function gradients. And then, this JxW, that's supposed to be J times W. J stands for the determinate of the Jacobian, and W stands for the quadrature weights, so it's all three quadrature weights multiply together in 3D. All right, so it takes the care of all of that 4s. The gradients we use in Klocal, the values we would use in Flocal if we had a body force. Okay, and course J times W we would use in any volume integral. If we move on to fe face values, you'll see I'm updating values. Quadrature points in J times W values. J times W values is as before. Only now you'll note that since it's a surface integral, it'll be the determinant of the Jacobian mapping from the 2D by unit domain to the 2D surface domain, right? update_values is, of course, still the basis function value on that surface. update_quadrature_points, this actually will give you the position vector of each quadrature point in the real domain, all right? So we need that in this problem because the traction for a Neumann boundary condition depends on the X1 component. Okay, so it varies as X1 varies on that Neumann surface Z equals 1. Okay, so that's why we're updating those values. So, let's scroll down into our element loop. Now, you can see here the first thing is that we do fe_values.reinit(elen), so fr_values is reinitializing for this particular element. So it's getting all the correct values for the quadrature, for the Jacobians and so on. The first step that we're going to look at is defining Klocal. Now before you look any any further at the code, let's look at the board at what that general form is. Okay, so in class we, or in the lectures we looked at Klocal and we had four indices that we were dealing with. Klocal AB ij. And the idea was that AB run over the nodes in your element. Okay, so I will go from 0 up to just less than the number of nodes, In the element. Okay, I and J are nodal degrees of freedom. And so that we run it from 0, 1, and 2 in this case for the 3D elasticity. All right, and the way we pictured that is that we had our Klocal matrix and then inside, we had these little sub-matrices. All right. And so, again there are eight nodes in a hexahedral element and so, these went from 0, 1, 2, 3, up to 7, all right, the same on the side, 0, 1 to 7. Now within the submatrices, we had indices 0, 1, 2, 0,1, 2. Now those correspond to the degrees of freedom at that particular node, okay? Now of course Klocal is not a matrix or matrices in d0l 2. It's just a matrix that's 24 by 24, all right? So we can look at this instead. This makes it, as 0, 1 2, 3, 4, 5, and so on up to 21, 22, 23. Okay, so let me label these. These would be our element degrees of freedom here. Going from 0 to 23, then we have our element nodes, Here. Going from zero to seven. And then we have our nodal degrees of freedom. Here. Going from zero, one, and two, all right? So you will notice when you look at the code again that our loops are actually looping over the element nodes and the nodal degrees of freedom. However, the indices in Klocal will be in terms of the element degrees of freedom. Okay, so how do we do that conversion? Because the indices start at zero it actually makes it simpler in this case. So we would do it like this. Klocal ABij corresponds to [3A + i, 3B + j]. I guess I should do square brackets on both of those, right? And you can see that that's true. Let's do it for 22 here. So that would be A is equal to 7, so 3 times 7 is 21. The nodal degree of freedom would be 1 at that point. So 21 plus 1 gives us an element degree freedom of 22, okay? And we'll work the same way for Flocal. For F local, it'll be F local A sub i, would be the same as F local of 3A, 3 times 8 plus 9, we'll look at that again in a second. Now let's look at how we actually calculate Klocal ABij. Let me write out the formula here. So we have Klocal ABij = the integral over the domain of the element, basis function A, the derivative to with respect to X of j, Our elasticity tensors Cijjl times our basis function B, ldV. All right and notice, we have repeated indices here. So there's an implied summation over j and l. For j and l going from 0. One and two, okay. All right, so let's look at what that will be in our node, in our code over here. All right, so first off, in our loops, we have a loop over quadrature points. Notice it's a single loop over quadrature points. d0l2 has combined all three loops since this is a 3D case into a single loop, all right? And it's keeping track of that for us. Now we have our loop over A and i, which is the loop over nodes A and then loop over nodal degrees of freedom. B and k again looping over element nodes and element or nodal degrees of freedom. Now, actually, I need to come back here and make a small change that I noticed. That's actually not the way I've written it here. That's not for Klocal ABij. It's actually for Klocal ABik, and the reason you should be able to see that is because j has already been removed through the summation. ink are three variables here, and so that's what should show up in the indexes of Klocal, all right? So that's why I have A and i grouped together, then B andk. j and l,, we have a loop there because of the implied summation over j and l. Okay, and now in here you'll define Klocal, using of course this integral that we're written out. Now one thing that's important for you to know is that when we use fe_values to get the gradient, which we do using fe_values.shape_grad() It's actually given us the gradient with respect to x and the real domain. You'll notice in the previous assignments when we wrote the basis gradient functions, the basis function gradients we wrote the gradient with respect to c and the bi unit domain and then we had to find the Jacobian. Do the inverse and so on, right? That's all taken care of here. Okay, so we don't have to deal with the Jacobian in order to make that a gradient with respect to the real domain, that's already done for us. Okay, so in order to access that gradient, you do fe_values.shape_grad, and then you'd input the element degree of freedom. Okay now, notice that's a little bit different than what we've written on the board here. In the lectures we always use the element node number to designate what basis function we're using. Here in d0l2, it's the element degree of freedom, okay? So it will be the same as whatever index you're using in Klocal. So this wouldn't be, for example, this wouldn't be A, this would be 3A plus i. And this wouldn't be B, it would be 3B plus j, for d0l2, okay? Sorry not j, plus k, all right. Again, you'll need to use the elasticity tensor function that I created before. And in order to use det to get your dertiminate of j times the quadrature weights we'll use fe_values.Jxw(q)*/ for this particular quadrature point. Notice also that we aren't inputting the value of c at the quadrature point, we're just inputting this index q which tells you what quadrature point we're at in the loop. Okay, also, to clarify, fe_values.shape_grad gives you a position vector. It's actually a d0l2.object, but you can think of it as a position vector. Or sorry, it's actually a first order tensor, which is very similar to a point and d0l team. Okay, so if you want to use a particular component, which of course you will, you can use just these square brackets. i or j or whatever it may be, okay? So, that should cover it for creating Klocal. Let's move down to Flocal. Now we don't have a body force in this problem or any forcing function like that. But we do have Neumann boundary condition. And so that will involve an integral, Over the Neumann boundary, okay? So let's write that out here. So for Flocal A sub i, which again is the same as in our code, it will be 3A or times A since is equal to 3, + i = this integral over a surface and so I'm going to designate that with this partial omega. And I'll use T to specify its attraction. It's the Neumann condition, and it's for a particular element, okay? And we'll be integrating hiI, so h is the traction, i is the component of the traction, the traction is the first order tense of our vector, times NA, and again it's a surface integral. All right? Okay, so let's go back to the code. You'll notice first, I'm doing a loop over the faces of the element, the current element. Of course there are six faces. So we loop over those and we're going to update fe_face_values for this, not only this current element, but also for the current face. That will put this quadrature points on the face itself that we're on. Now we're going to check to see if that face is at our Neumann boundary, okay. So this lm arrow face arrow center, that gives us a position vector at the center of the current face. And since I'm interested in the z component, I use square brackets too. Okay so here in this if statement, I am checking to see the point at the center of this face at the boundary z = 1, which is our Neumann boundary, okay? If it is, then I'll perform this surface integral, okay? If not then I will just move on. Either to another face, and if none of those faces were Neumann faces I'd move onto another element, and Flocal would be zero for this element, okay? But once we are, once we have found that Neumann face and we are performing that surface integral, we'll move inside and loop over our face quadrature points. Okay, notice that, the number of face quadrature points. Here I've extracted for you the value of X, the x-coordinate of the current quadrature point, okay, and I've done that using this fe_face_values.quadrature_point(1)[0]. If I wanted the z component then I would just change that zero to a two. If I wanted a y, I would change it to a one, okay. But then what you need to take that value of x and specify what the value of the traction vector h is at that quadrature point, okay. Now, once you've done that we'll move inside our loop over the nodes and the nodal degrees of freedom, okay. And, once inside there you will again use this integral that we've written out on the board to define Flocal. Now you'll notice I am looping over all eight nodes of the element. Even though we're only doing a surface integral. How are the quadrature points themselves, D02, is placed on the appropriate, correct, surface the correct face. So, we are only integrating over the surface itself. All right, notice again, we'll be using fe_face_values to find the basis function value. Okay, and you are still passing in the element degree of freedom. So that's 3 times A plus i. Not just A, okay. And again when you're finding JxW, the determinant of J times the quadrature weights. Again use fe_face_values, all right. Everything involved with the service integral will be using fe_face_values, not fe_values at all, all right? So once you've filled that in, you'll have Flocal created as well, and once you have Flocal and Klocal. A symbol system will be very straightforward. It will be exactly the same as previous assignments, all right. Applied boundary conditions will be applied in the same way. Again, you've already specified your knowing your boundary condition. Fixing all degrees of freedom, X equals 0. And I'll scroll down a little bit more just to finish up quickly. Solve is exactly the same and output results again is the same. Outputting our displacement results as a .vtk file which you can then open up using para view or visit to look at the displacements. All right so that concludes our this segment and it concludes our discussion of the template for homework four
