This course is an introduction to the finite element method as applicable to a range of problems in physics and engineering sciences. The treatment is mathematical, but only for the purpose of clarifying the formulation. The emphasis is on coding up the formulations in a modern, open-source environment that can be expanded to other applications, subsequently. The course includes about 45 hours of lectures covering the material I normally teach in an introductory graduate class at University of Michigan. The treatment is mathematical, which is natural for a topic whose roots lie deep in functional analysis and variational calculus. It is not formal, however, because the main goal of these lectures is to turn the viewer into a competent developer of finite element code. We do spend time in rudimentary functional analysis, and variational calculus, but this is only to highlight the mathematical basis for the methods, which in turn explains why they work so well. Much of the success of the Finite Element Method as a computational framework lies in the rigor of its mathematical foundation, and this needs to be appreciated, even if only in the elementary manner presented here. A background in PDEs and, more importantly, linear algebra, is assumed, although the viewer will find that we develop all the relevant ideas that are needed. The development itself focuses on the classical forms of partial differential equations (PDEs): elliptic, parabolic and hyperbolic. At each stage, however, we make numerous connections to the physical phenomena represented by the PDEs. For clarity we begin with elliptic PDEs in one dimension (linearized elasticity, steady state heat conduction and mass diffusion). We then move on to three dimensional elliptic PDEs in scalar unknowns (heat conduction and mass diffusion), before ending the treatment of elliptic PDEs with three dimensional problems in vector unknowns (linearized elasticity). Parabolic PDEs in three dimensions come next (unsteady heat conduction and mass diffusion), and the lectures end with hyperbolic PDEs in three dimensions (linear elastodynamics). Interspersed among the lectures are responses to questions that arose from a small group of graduate students and post-doctoral scholars who followed the lectures live. At suitable points in the lectures, we interrupt the mathematical development to lay out the code framework, which is entirely open source, and C++ based. Books: There are many books on finite element methods. This class does not have a required textbook. However, we do recommend the following books for more detailed and broader treatments than can be provided in any form of class: The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, T.J.R. Hughes, Dover Publications, 2000. The Finite Element Method: Its Basis and Fundamentals, O.C. Zienkiewicz, R.L. Taylor and J.Z. Zhu, Butterworth-Heinemann, 2005. A First Course in Finite Elements, J. Fish and T. Belytschko, Wiley, 2007. Resources: You can download the deal.ii library at dealii.org. The lectures include coding tutorials where we list other resources that you can use if you are unable to install deal.ii on your own computer. You will need cmake to run deal.ii. It is available at cmake.org.
03.08ct. Coding assignment 1 (main1.cc, overview of C++ class in FEM1.h)
Loading...
From the course by University of Michigan
The Finite Element Method for Problems in Physics
215 ratings
From the lesson
3
In this unit, you will write the finite-dimensional weak form in a matrix-vector form. You also will be introduced to coding in the deal.ii framework.
Meet the Instructors
Krishna Garikipati, Ph.D.Professor of Mechanical Engineering, College of Engineering - Professor of Mathematics, College of Literature, Science and the Arts
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.
Explore our Catalog
Join for free and get personalized recommendations, updates and offers.
Coursera provides universal access to the world’s best education,
partnering with top universities and organizations to offer courses online.
© 2018 Coursera Inc. All rights reserved.
