Hi. In this segment we will begin looking at the coding template for the first coding assignment which is homework two. All right. So, before we look at the code itself, think to yourself about what are some of the functions and data objects that you'd want to include in your C plus plus program. All right, think about that on your own for a minute. Okay? Now, let's write down some of those. First, we would have our basis functions themselves. Right, so those would be the functions that are interpolated for you to create your finite element solution. And along with the Basis Functions would also be the gradients of the Basis Functions. If you look at your first assignment, you will also see that another function that we need to create is a function that will calculate the L2 Norm of the error between the actual solution and your finite element solution. So that's another function that you will need to create. Now how about the steps of the finite element method? What are the different steps that we go through? Well, we have to generate a mesh. Or import a mesh of some sort. Right? We need to create F local and K local At some point. And once those are created, then we have to assemble F local and K local into the global stiffness matrix and the global forcing vector. All right, and then we also have to define and apply our boundary conditions. Now those can be Dirichlet boundary conditions or Neumann boundary conditions. Okay. And after that's created have K and F, then we have to solve for D. So we have to have some function that solves, for D. So D = K inverse F. And then you'll probably want to output your result somehow into some format. Alright, so these are some of the functions that we would have. Some of the data objects themselves, what we've already mentioned here, would need our stiffness vector, or our stiffness matrix, K. We need F and D as well. Our foreseen vector, and our solution vector. We would need some sort of object that would relate our local degrees of freedom to our global degrees of freedom. Sometimes we call that a connectivity matrix in general. Although it may not actually be a matrix so I'll just call that a connectivity object. We, of course, would need some sort of object that actually stores the information in the mesh, which would be related to this connectivity object. So, I'll include that here in mesh. All right, there may be others, too. We'll look more at that now as we look at the code, okay? So let's move over and look at the code itself. As we told you in the homework file in the homework description, we're giving you two homework files. There's the .cc which is the source code and a .h file which is the header file. The header file is what includes your finite element class. All right. So, that contains the bulk of the programming. But to begin with, let's look at the main .cc which is the main program that's running. All right. So, here you can see I've included some standard C++ header files, but I've also included this FEM homework to template.h. That's your header file that includes the information about the class. And here on line 12, you can see that I said this using main space deal.II. Now the reason for doing that is because we have a lot of deal.II objects and normally, we would have to type in deal.II, double colon, and then the name of the data object. Just like with standard. For example, standard vectors we type in SDD::vector. If we instead wrote using name space standard, then we wouldn't have to do the standard colon colon every time we wrote vector or c out, for example. In this case, we're using namespace deal.II. All right. And then our main program. Okay? This first part, it just has to do with some of the mechanics going on. Our program really starts here. First we define the order of our basis function. That will be an input to the class constructor. Here I have order equals one but again for the first homework assignment, it can be linear, quadratic or cubic. Okay? So either order one, two, or three. Also there are two parts to the homework problem, part one and part two. And so again that's also an input. The only difference between part one and part two is in the boundary conditions. Problem one has two Dirichlet boundary conditions and problem two has a Dirichlet and a Neumann condition. Alright here on line 25 is where we actually declare our object. The name of the class is FEM. It's a templated class and all I'm going to explain about that is that there is a template parameter that goes into the class that the functions and data objects are based off of. In this case it's the dimension. Okay, so I put in a one here between the angle brackets. Just specify that it's a one dimensional problem that we're working on for this homework assignment. I'm naming it problem object but again that's just a name I'm giving it, it's not the class name, it's just the name of this particular instantiation of the class. And again the inputs are order, the lagrounging basis function order and the problem number. All right, so now we get into some of the actual functions. Remember, in order to call a function for a class, we first include the name of the object and then a decimal and then the name of the function. Alright, so the first function that I'm calling here is generate mesh, and as an example here I am inputting 15. That just tells, or this is where we specify how many elements we want in our mesh. So in this case I am saying a 15 element mesh, of course you can change that. In this case, with the way I've set it up, this will make 15 equal elements. Okay, all the elements are of equal length. Then I call this function, called setup system. Of course this doesn't really come in from the finite element method, but it's mostly for initializing our data objects. Resizing the force vector, the solution vector, the global stiffness matrix. That's also where we will be setting up, or defining our quadrature rule. We'll decide how many quadrature points we need per element, and then we define what those points are and what the corresponding weights are. After that is assemble system, and it's within assemble system that we actually create K local, F local and then assemble those into the global stiffness matrix and the global force vector. All right. We've previously defined the Dirichlet boundary conditions within setup system, the Neumann boundary conditions, you'll see we define as we're creating F local. And then we actually apply the Dirichlet boundary conditions within assembly as well. Solve is a very short function that simply solves for D equals cambers F. Here, I am outputting to the screen, the L2 norm of the error. For this, for the given basis functions, and for the given number of elements. And finally, we output results. We output results as a .vtk file. And for the linear problem, there's actually not much to see, since it's just a straight line. But when we move on to the 2D and 3D problems, then I will show you how to open up one of those files, in a visualization tool, called ParaView. Okay? So let's move on to the header file. Now scroll to the top here. [COUGH] Now in the comments you see at the top. I've specified that anywhere you actually have to make a change or modify the template I've included the word edit in capital letters, okay. So, anywhere you see edit you'll have to edit something, otherwise you can let it be. Also, as you're writing your code don't change any of the function names or inputs that I've already set up for you. You are, however, free to go ahead and create new functions, or additional functions, if you find that helpful. So there are a lot of header files here that you can see in the code. And almost all of them are DO2 header files. We actually won't necessarily be using all of them in these first homework assignments, such as quadrature point. We'll be creating our own quadrature points, we'll be creating our own basis functions. However as the class goes on, in the last two homework assignments you will be using more and more of DO2's actual functions. Okay. Now let's scroll down to the declaration of objects. So I've created everything as public, so that means we can access all of these objects and functions from main.cc. All right, now to begin with here's our constructor, and we've already seen this constructor in main.cc. We have as an input the order, the order the Lagrangian basis functions, again, that's one, two or three. For this homework assignment, of course it could go higher if you chose. And then the problem, again that's part one or part two. Here we also have a class destructor. No inputs to it, but you will see that there is something going on in the destructor. Now, here we get to the functions that I've set up for you. The first function is a function called xi_at_node. And I'll explain a little bit more what that means when we get down into the actual definitions of the functions, here we're just declaring them. But it has to do with, the differences in the local node numbering between DL2 and what we've seen in the lectures already. So this function will help you get the value of xi using the DL2 node numbering. Here we have our basis function and basis gradient. And now down to the solution steps. We have generate mesh, which we solved before already. Define boundary conditions. That's actually called from within setup systems, so that's why we didn't see that in main.cc. Then again setup system, assemble system, solve, and output results. And again, after that, we have the l2norm_of_error functions. So those, the rest of those are functions we saw already in main.cc. That concludes the functions. Now, we have some class objects. These first three are objects, are DO2 objects that they've, that they've defined and we're declaring them here. Triangulation basically stores the mesh. FE and DOFhandler have to do with the connectivity, keeping track of the degrees of freedom, as well as general basis function ideas going on. Okay? Now we get into our own data objects. And these next three have to do with the Gaussian quadrature, that we'll be using to perform the integrals, that occur in K local and F local. All right. So the first is an integer called quadro. It's just the number of quadrature points. In your element. Now, you'll have to decide how many quadrature points, and we'll look at that a little bit more later. And then we define what those quadrature points are, again in 1D, and the corresponding quadrature rates. All right. So, it's an integer and then the next two are standard vectors of doubles that we've looked at before, in our review of C++. Okay. More data structures. The first one is a sparsity pattern. And we use a sparsity pattern with our global stiffness matrix, K. Because as you remember from the lectures, most of K are zeroes. And so, when we're dealing with a large stiffness matrix it doesn't make sense to store a lot of zeroes and so we use what's called a sparse matrix and it only, it essentially stores just the non zero components of K global okay. Sparsity pattern is what keeps track of where those locations actually are in the real, in the full matrix itself, okay. We have a vector of doubles, D and F, that's our solution, solution vector and our global force vector. You'll notice this isn't a standard vector. It's a DO2 Vector. There's a capital V. And so the functions related to vector, the DO2 Vector are a little bit different than the standard vector. And mostly we use that in our solve function. When we're solving for the solution vector D. Okay. Now, I have this standard vector of doubles called node location. This is a vector that just gives you the x coordinate, okay, and the real domain of each global node. And then, a map. A standard map called boundary values. Let me give a little more information on the standard map. This is a C++ data object, the standard map. And as you can infer from the input there, you have to define two data types. In my case, in this case we've defined two number types. This first is an unsigned int, and the second is a double. And so this basically just, as you might expect from the name, it maps an unsigned int to a double in our case. One big difference though from a vector, is that a vector, in a vector the indices start at zero and they just go sequentially up from there. So zero, one, two, three, and so on. With a map, what you might call the index, doesn't have to be sequential. So in a standard map, the standard map of boundary values. So we have a standard map of unsigned int, and that should be. Let me erase. Should actually be an angle bracket. Okay. Unsigned int and double. Called boundary_values. Okay. So in boundary_values, we will be relating the Dirichlet boundary node number, the global node number, and the actual value. So, if we have, if this is our domain, and we're breaking it up into elements. Let's say at the left we're applying a Dirichlet boundary node of g1 = 0 and at the right we're playing a boundary node value of, a Dirichlet boundary value of g2 = 1. All right? So all we need to, and I'll let these be our global node numberings so I'll just go sequentially, zero, one, two, three, four, and five. So in boundary values we want to store the node number, with the Dirichlet value. So the node value is an unsigned int. The Dirichlet value is a double. Okay? And so all we need to do here is we do boundary_values. The index, this first term, you can think of as an index. And so at node zero we have a boundary_value of again zero. Alright so there we've defined that, we've set that into this boundary values map. On the right, at node 5, we have an applied boundary condition of one. Okay, and so now, whereas in a vector we would have had to have a vector of length 6, in order to have an index of 0 and an index of 5. With a map, a standard map, this only, if you were to do a boundary_values.size, it would return a size of 2, because we're only storing two values in this map, okay? So let's go back to the code. After the boundary values map, I've defined several doubles. There's the basis function order, problem, or prob. Those again will store the Lagrangian order that we're using and the problem, part one or two. Then I have L which is the length of the domain. G1 and g2, which again are the Dirichlet boundary values. G1 at the left and g2 at the right. Okay, now under here, solution name array, these are just vectors that have to do with the output, so you actually won't be dealing with these at all. Okay. Let's stop this segment here and in the next segment, we'll start looking at some of these functions themselves.
