 # 11.09ct. 1. Coding Assignment 4 – I

In this segment, we’ll start looking at
the coding template for homework five. So let’s turn and look at the code now. Again, homework five is the transient
heat conduction problem. And you’ll notice some differences
here in the main.cc file. One is that I’ve created this
constant alpha, which is an input for the constructor for
our FEM problem object. Alpha is to specify what
Euler method we’re using. Again, alpha can range
anywhere from zero to one. For the homework, you’ll be doing it for
alpha equals zero, one half, and one. But, other than that,
the mesh generation is very similar, and setup_system and assemble_system. However, now we also
have to solve functions. We have solve_steady,
which solves our steady state problem, and then solve_trans which, obviously,
a function that calculates the L2 norm of the difference between
the steady state solution and the transient solution at a given time. And so
I’m also outputting to the screen here, some of those results of the L2 norms,
there will be four L2 norms, It will calculate the L2 norm at time
equal 0, 1,000, 2,000, and 3,000. Those are some of the results that
we’ll be using in the grading. Let’s shift over to the header file now. Just as in homework four,
I’ve defined order and quadRule as global
variables at the top here. Again, when you turn in your assignment, use order equals one and
quadRule equals two. But for your own edification, you can [LAUGH] change that order and
quadRule if you’d like. The L2, again, will take care of
any any order our quadrature rule. We’ll scroll down and look at
the declaration of functions and objects. The first four solution steps are the
same, but then we have this solve_steady and we have apply_initial_conditions,
since it is a transient problem or includes a transient problem,
we have to include initial conditions. So you’ll be defining those
within this function. That function is called within
the solve_trans function itself. I now have two output results functions,
one for the steady state result and
one for the transient results. Notice for the output_trans_results,
I have as an input an index. So that’s just to distinguish
between the output results for, say time equal zero and time 1,000, it
just puts an integer onto the file name, so that it doesn’t overwrite
the solution file. And again now we have another
L2 norm function again, so if you’ve missed that from homework two,
it’s back. You’ll notice here we only have
the quadrature formula for volume integral; since we don’t
have any Neumann conditions, we aren’t doing any surface
integrals on this problem. You’ll notice I have more matrices and
vectors here than before. We have capitol M,
which is our global mass matrix. There is again K and
then we have a generic system matrix. That’s because when you’re doing the
transient problem, you have to do a couple of matrix inversions to solve,
for example, the VM plus one. And so, system matrix is a generic
matrix that you’ll be defining later on, in the solve_trans function. We also have D_steady, which is
the steady state temperature profile. D_trans, which would be the transient
temperature profile at Tn or Tn plus one. V_trans would also again be similar. It’s the transient time derivative
of the temperature at each node. Either for time Tn or Tn + 1, depending
on what part of solve_trans you’re in. And then again our vector F and RHS. RHS is a generic right hand side vector,
similar to system matrix that you’ll be populating within solve_trans to solve for
your Vn + 1 vector. Okay? We also have this standard vector
of doubles, called L2norm_results. This is going to store the L2norm
at various time steps. And then alpha. Alpha’s just the input,
from our constructor, ok? Specifying again what
Euler method we’re using. So if we scroll down to the constructor, again you can see we’re inputting alpha
and storing that within our class and very well alpha And
if we scroll down to generate_mesh, this again is exactly the same as
previous assignments where again, you just need to define
the limits of your domain. Define boundary conditions, would be
very similar to previous assignments. Very similar to homework three, which was
your previous heat conduction, your study state heat conduction problem, except now
you’ll notice that I have two maps here. Boundary_values of D and
boundary_values of V. I actually skipped over those higher up. They are class variables, class objects. And since we have these two
fields that we’re dealing with, we’re dealing with
the temperature distribution and we’re dealing with the field of
the time derivative of the temperature. And so we have boundary
condition on both of those. The boundary conditions on D are 300
at one end and 310 at the other. They are fixed temperatures, so what does that tell you about
the boundary conditions on V? Which is the time derivative
of the temperature. Tells you, of course, that the derivative
of D would be zero at those boundaries, where temperature is fixed, okay? So you’ll be doing a similar loop to check where your node is to see if it
is at an Dirichlet boundary. If it is at a boundary, you’ll put
that information into these two maps, using the globalNodeIndex, which is the same as degree of
freedom index for this scalar problem. And you’ll passing the prescribed
temperature value and the temperature time derivative,
which again we said would be zero. When you have,
at fixed temperature, alright? So that will take care of
your boundary conditions. Again set up system is very similar to
previous problems, there’s nothing for you to change here. Let’s look at a symbol system very
quickly, before we end this segment. This is going to be sort of a mixture of
the previous two homework assignments, your heat conduction problem and your, and using FEvalues in the linear
elasticity problem. Again, we have FEvalues here,
we’re going to update the values, update the gradients, and the Jacobean determinative of the Jacobean
times the quad digital weight values. I set up this variable row and
that’s the specific heat per unit volume. It’s a constant that we’ll be giving
you in your homework assignment, that you’ll use to define M local, and
then to get your global mass matrix. All right so
let’s move inside the element loop. And first we’re going to
be populating M local. Let me write out that integral for
you here on the board. MlocalAB is equal to the integral over the element’s domain, of row, which again
is your specific heat pre-unit volume. Now, potentially that could
vary across the domain. For this homework problem we’re
keeping it as a constant, so I can pull that out of the integral. But then, it’s simply the basis
function BA times the basis function B. Integrated. Okay, so this is fairly straight forward,
and from your experience doing Klocal and Flocal this should be pretty
straight forward for you. Pretty simple. Okay, so
you’ll edit that there to find that here. We scroll down and
we have to create Klocal. This will be the same Klocal as
you’ve created in homework three. The only difference is you don’t have to
with respect to the real domain. So that simplifies our
quadrature loops a little bit. All right so we have looping over AB. The quadrature points, which again all three quadrature loops are grouped
into a single quadrature loop with. Then we do i and
j to perform the multiplication with 0ur conductivity tensor, kappa. Alright? We don’t have a foreseen function in
this problem, so Flocal will be zero, and so we don’t even bother in this
assembly to assemble F, our global vector. But you will be assembling K and
M very similar to previous assignments. Right. The solve_steady function
acts very similar as the solve functions in
previous assignments. The small difference that I’ve done, is