Hi. In this segment, we'll start looking at the function for generating the mesh, and for defining the boundary conditions. All right, so let's look at the code. Generating a mesh, there's actually not a lot for you to do. I have this global variable L, you remember I declared above, it's a class object, a class variable. So you simply have to define the length of the domain. I'm assuming that the value of x at the left of the domain is zero. And then of course, the maximum value of x would be L. From there, I take the x_min and x_max and I put them I create these DL2 points. A point in DL2 is, it's essentially a location vector. A position vector. The dimension or the number of components in the point is defined by dim in this case, and it holds doubles. I have actually, I don't use point a lot in this template, you won't actually need to use it at all yourself. But that is the input that DL2 requires to create the mesh. All right, so I take the limits of the domain, which you define here with L. And remember the input is an unsigned int called numberOfElements. In main.cc, remember, I had put in 15 as an example. But that just gives the number of elements in our mesh. All right. And then DL2 will use that number to create a mesh. It stores that again in triangulation, which again is a class object. All right. So not a lot for you to do there, but an important step. A lot of it is, again, DL2's functions. All right. Now, in define_boundary_conds, there actually isn't anything that you need to do in this function for this assignment. However, for the future assignments, you will be writing the defined_boundary_conds function yourself. It will be very similar to this function, but again, you have to change anything here, but I want to step through with you now because again, you will be using in future assignments. All right, so here we're using that boundary values map that we talked about, and we already talked a little bit about how the map itself works. However because of the way DL2 does their node numbering, instead of simply saying okay, we know how many nodes there are, so I'll just do boundary values at zero, is equal to the fixed boundary condition of zero. And if I have a Dirichlet boundary condition on the right, then I'll just take the last node. But again, since they aren't numbered sequentially, we can't just take the highest node number. Because that's not necessarily the last node in the domain. Right? So again, as an example, let me take Let's create a domain here. And I will, for simplicity, I'll create four elements. And on the top, I'll do the x coordinate. So, this isn't the real domain. We're going zero. And I'll go from zero to one. That would make these 0.25, 0.5, and 0.75. And just to make my point about the node numbering, I'll make this quadratic, quadratic basis functions. And so we have one mid side node. Down here, I'll put the node number. Of course in DL2 node numbering. Okay, so it's 0, 1, 2, 3, 4, 5, 6, 7, 8. Again so, this node, assuming we're doing two Dirichlet conditions, we want to fix this node. And this, let's say we want to apply a displacement of 0.1, that's a pretty big displacement here actually. Since it's linear elasticity, I'll do 0.01. So let's say, g2 is equal to 0.01. All right, so what we want to do is boundary_values of node zero is equal to zero. And boundary values of the last node, which in this case is actually node 7, not node 8, is equal to 0.01. These are semicolons, and so since we don't necessarily, okay, for the one d case, we can keep track of the way DL2's numbering their nodes. However, once we get to 2D and 3D, it gets more complicated from there. So instead what we're doing and what I'm doing in the code, is I just have a for loop. I'm looping over all of the nodes, so you can see here I have unsigned int globalNode=zero, and I'm looping up through nodeLocation.size. Remember, nodeLocation.size is a vector of the locations of the nodes. All right, so let me write out what nodeLocation, what this vector would look like, in our case, the actual data that it's storing. Okay, so we have nodeLocation. And this is what it's holding. Okay, so its hold for node zero, again looking at the indices are these node numbers down here. So for node zero, we're holding the value at zero. For node 1, the x coordinate is 0.25. For node 2, it would be 0.125. All right, and that goes on, until node 7, which has a value of 1. And then node 8, which would have a value of 0.875. Right. Okay. So those are the entries in nodeLocation, okay? And so the length, you can see it has a size of 9, which is of course equal to the number of nodes in this system, okay? So, here in the function define_boundary_conds, we're simply looping over all of the nodes using this nodeLocation vector. And then we have an if statement. If the value stored in the node location at the current node that we're on is equal to zero, in other words, if we're at the left end of our domain, then we set boundary values of whatever the node index is equal to g1. And actually, I can do here an else if, if I wanted to, else if, or just if nodelocation at our current node, wherever we happen to be, is equal to L. In other words, if we're at the right side of our domain. And if problem equals 1. So problem 1 is where we have two Dirichlet boundary conditions. So if we have a Dirichlet boundary condition in our problem, and we are at the right hand side of our domain, then we will store the Dirichlet value of g2 in our domain. So that's this one right here. Okay, so again, there's nothing you have to edit here. But it's good for you to understand what's happening so that you can use it on the next homework assignments. All right. Let's move on, and in this segment, let's still look at set up system. Okay. So in set up system, first we have to define. Again, we were using g1 and g2 in the previous function, then here's where we actually define what those values are. And you'll have to look in your homework assignment to see what those values will be. Of course in problem 1, we use both g1 and g2, in problem 2, g2 won't actually be used, because its a Dirichlet Neumann problem. But you can still use the same value of g2, right. Okay. This next line is just DL2 distributing the degrees of freedom. It's just keeping track of degrees of freedom. Now here, there are several lines going on here, but essentially all I'm doing is creating this vector node location which holds the x coordinates for each global node, again by its global index. Now that those structures are defined, now we will call this function define_boundary_conds. Okay. Notice that since we are still within another class function itself, we don't have to create like, problemObject. We don't have a class name, since we're not dealing with an actual class object yet. So we don't have this problemObject.define_boundary_conds(). We can simply use the function name itself, because it is also an FEM class function, okay? All right. So we call that function, which is the one we were just talking about. From there on, this next section just has to deal with resizing the matrices and vectors. You can see rather than using the function.resize, which is what standard vectors use, the DL2 vectors and matrices use .reinit, short for reinitialization or reinitialize. And we use this dof_handler.n_dofs. That gives us the total number of degrees of freedom in our domain. Since we have one degree of freedom per node, it's the same as the total number of nodes in this problem. Now, here we define the quadrature rule. I set up a quadrature rule of two quadrature points per element. However, that's not going to be enough. Remember that for the quadrature rule, let's say you have a quadrule, whoops. Let's say you have a quadrule of 2. This is actually good for up to, for an integrand of up to third order. So, remember that Good for and integrand of third order or third order polynomial, right? So remember that this Gaussian quadrature is exact for polynomials, Of order equal to 2 times the quadrature rule, minus 1. So in this case if you have a quadrature rule of 2, it's good for a third order Because 2 times 2 is 4 minus 1. So if you have an integrand of order 3, with polynomials of order 3, then the Gaussian quadrature will integrate exactly for you. Okay. But remember, you'll be using not only linear, but also quadratic and cubic basis functions. And you'll be doing an L2 norm. Okay? So you have to, in each of these cases, don't look actually just at the order of your basis function. Look at the whole order of the integrand that you're integrating. And once you know that polynomial order, then you can decide what quadrature rule order you need. Now having said that, you don't have to necessarily change your quadrule for each problem. You can decide what's the maximum quadrature rule that you need. And that actually will come from the L2 norm using cubic basis functions. You can decide what quadrature rule that is. Is it four, is it five, is it six? Decide what that number is, and then you can use that quadrature rule for all of your problems on homework one. Okay, for whatever basis function ordering. All right. So again, you will define that here. So you would change quadRule to whatever value you want. And then quad_points and quad_weights are, again, just vectors. You put in the, so quad_points stores the values of c for each quadrature point. For two point, it's minus squared one third and positive squared one third. And the corresponding weights are one. Okay. If you go to a higher quadrature rule then you will obviously create a longer vector. You can see that here I do quad_points.resize, and I resize it according to the quadrature rule that you define up here. Okay, so it will automatically go to the correct length. You have to define the correct values. And those can be found online, or in a book on finite elements or on numerical integration, all right. Okay, so let's stop this segment here. And in the next segment, we'll move on to the assemble system function, the assembly function, which involves creating f local and k local.
