Tip:
Highlight text to annotate it
X
In this lecture, we are going to continue our 1D programming assignments. We had stopped
at the assembly of the element equations their application of the boundary conditions and
the point load for the one dimensional problem. We had raised a question of how to solve the
resulting system that we had, which is of the form, K U is equal to F and K was of dimension
nndofs by nndofs and it was symmetric.
As I had mentioned in passing at the end of the last lecture, there are 2 types of solvers
which are basically popular in engineering analysis, especially when we use FEM. One
class of solvers is called the Direct Solvers. The other class is called the Iterative Solvers.
The direct solvers work by direct elimination of the unknowns to get the solution. The iterative
solver, it works by first taking an initial guess of this U and then trying to correct
the solution then add correction to this initial guess to get to the solution of the problem.
Each one has its own advantage but for our case, for the system that we are interested
in, let us look at only the direct solvers. Within this class of direct solvers, we will
look at an important category of LU decomposition based solvers. What do we mean by LU decomposition?
Let us see this through an example.
Let us take this matrix problem which is a 3 by 3 problem given by this kind of a system.
Deliberately, I have made the matrix A or the K here a non symmetric, this is a general
problem . When we have to do Gauss elimination then we simply go and subtract; I want to
eliminate first the x1 variable from the remaining two equations not from the first one. I simply
subtract twice the first equation from the second equation. To do this job, I do this
and I get this one . I subtract once the first equation from the third equation. By doing
that, I have reached this form, where I have eliminated x1 from this equation and this
equation that is, from the second and third equation. First variable is eliminated from
all the remaining equations. This is how, if we remember, the Gauss elimination is done.
Once we have eliminated x1 from the remaining two equations then we want to eliminate x2
from the third equation so we are going downwards. From the second equation, we take this diagonal
and we would like to eliminate everything under this diagonal. Like in the first equation,
we have eliminated everything under the first diagonal. These diagonal elements are called
the pivots and we are going to simply eliminate this one from here, by simply taking the second
equation and subtracting it from the third one. If I subtract the second equation from
the third one, I will get to this form.
1, 2, 1, 0, minus 2, minus 1, 0, 0, 3, x is 1, 0, 2. If we see this one, this has 0 entries
in the part of the matrix below the diagonal. This type of a matrix where everything below
the diagonal is 0 is called an upper triangular matrix. We will call this by a name U. We
want to know how through a series of operations can we take A and get U. We have done those
operations. How can we write these operations in terms of matrix manipulations?
If we go back, here the first operation that we did was we took twice the first equation
and subtracted it from the second one.
This operation I would like to write through a matrix that I leave the first equation unchanged.
I subtract twice the first equation from the second one. I do not do anything third one
here and I am leaving the third one fixed. This if I multiply it with A does the same
job as subtracting twice the first equation from the second one of the matrix A. This
one I am going to call with a name matrix D1. Similarly, if I go and look at the second
part, second part from is here. To eliminate x1 from the third equation, I take once the
first equation that is the first equation itself and subtract it from the third equation.
This part will be given by a matrix; I have the first equation which remains unchanged.
I am doing nothing to the current second equation. It remains unchanged and I subtract once the
first equation from the third one. This is the subtraction process this I am going to
call as the matrix D2.
Similarly, the third part was, after we have done this, we have a resultant modified matrix
A. This is matrix A . I go and take the second equation and subtract it from the third equation.
How do I write that? This subtraction process can be written in the form. I leave the first
equation unchanged. I leave the second equation unchanged and in third equation, I am going
to subtract the second one from it. So I take the third equation from it to subtract the
second one. This is the process of doing this job. This is given by the matrix D3. In terms
of these three matrices, now can I write U? Yes, I can. If I take the matrix A, first
I operate on it with the matrix D1 that is the first subtraction of twice the first row
from the second row. Then I do the second operation D2 that is, subtraction of the first
row from the third row. Then I do the third operation D3, subtraction of the second row
from the third row. This job should give me the upper triangular form U. If I want to
now write A in terms of U is equal to, I will take the inverse of this set. Inverse of this
will be D1 inverse into D2 inverse into D3 inverse into U, finding the inverses of these
matrices is very simple. Let us see how simple. This quantity I am going to call by a name.
I am subtracting this quantity from the essentially, tells me by a how much amount, I multiplied
the first equation to subtract from the second one. This one I am going to call l21. Similarly,
this quantity I am going to call as l31.
This quantity I am going to call as l32, it tells me from the third equation. This is
the amount by which I multiply the second equation and subtract. It turns out that if
I look at D1 inverse, it is quite simple. D1 inverse will remain 1, 0, 0. We had subtracted
twice the first equation from the second. I simply have l21, 1, 0, 0, 0, 1.
Here I have taken D1 and simply taken this part which corresponds to subtraction and
taken the negative of that, put it there, back into the matrix that becomes the inverse
which as simple that.
Similarly D2 inverse will be 1, 0, 0. Here I will have l31, 0, 1. The third one will
be 1, 0, 0, 0, 1, 0, 0, l32, 1. So finding the inverse of these matrices D1, D2, D3 is
very easy. Once I know what the operation through which I can to this form.
These matrices have a nice property that these are only concentrated in the lower triangular
part of the matrix that is, everything above the diagonal is 0, everything below the diagonal
is non zero. The product of these three things will be matrix L which will have once on the
diagonal l21 in the half diagonal here.
l31 in the half diagonal here. l32 in the half diagonal here, 1 here. Now I have obtained
a lower triangular matrix and the diagonals are all 1. This is the beautiful property
of this matrix L . I can write the matrix A as a lower triangular matrix L into an upper
triangular matrix U. Let me get back and look at what was U.
If we see U, the diagonal entries of U are not 0. So we can make by a suitable modification,
the diagonal entries of U also 0. How can we do that?
We can now say matrix U is equal to U11, 0, 0, 0, U22, 0, 0, 0, U33 into 1, U12 by U11,
U13 by U11 here I will have 0, here 1, U23 by U22, 0, 0, 1. This I am going to call a
diagonal matrix D. I have simply taken the diagonal entries and written it in a diagonal
form here and if I now multiply I should get back the U. This matrix, I am going to call
it as U bar.
I can write A as equal to L D U bar. So this is called the L D U decomposition of A. This
is something that we can also do. Let us get back. What is the advantage of this? Even
U bar has one on the diagonals. Before go ahead now we see how this decomposition is
going to be advantageous. If I have A x is equal to b; instead of A, I write L U. I will
call y is equal to U into x, implies we will have L y is equal to b. Here we have the lower
triangular matrix into y is equal to b, solving for y is very easy, because here it is explicit.
We do not have to do any elimination explicitly. We will first start with a first term, find
y1 put that into the second equation, find y2 and so on; solving this is very easy.
Once I have obtained y then I will obtain x as U x is equal to y. Instead of operating with
A, now we can explicitly operate with the L and U. This makes life very simple because
if I have a stiffness matrix, I want to change the loads applied on the structure keeping
the boundary conditions of the material fixed. Then the stiffness matrix is going to remain
the same. As the load is changing, the load vector which will be this side is going to
change. When the load vector changes then I do not have to eliminate A, do a Gauss elimination
on A all the time. I can have the L and U in store and simply take the new load vector
operate on it with the L and U and I get a solution. So that is something we would like
to have.
Let me now get back to the L D U decomposition. We had said A is equal to L D U bar. This
turns out to be unique. There is only such decomposition which is present that is, if
I have another decomposition and that decomposition equals A then each of these terms have to
be equal. We see something very interesting; when A
is symmetric then A is equal to A transpose due to symmetry. This is equal to U bar transpose
D L transpose, but we know that there is a unique decomposition of A . We know that for
a symmetric matrix U bar transpose is equal to L and L bar transpose is equal to U bar
is just consistent. So then for these cases, for symmetric matrices, the decomposition
becomes very simple. It becomes L D L transpose. This we can write as is equal to L D to the
power of half into L D to the power half transpose. I take this part. Take the transpose of that
multiply by 2. I get this decomposition for A is symmetry. This decomposition is called
Cholesky decomposition. To take this information about what is L U D decomposition base solver
is going to do and take a commercial version of this solver which is available in various
places.
Again we can go to netlib.org, get the symmetric form because the matrix is symmetric that
is L D L decomposition of the source code for this decomposition and in whichever form
we want, Fortran, C or C++ or we can go to the book. This is a fantastic book in numerical
analysis, numerical recipes. This is available in various versions. It is for Fortran, C
and C++. On all these versions, this book is available. In these numerical recipes,
source codes are given; the programs are given. We can take the program from there for the
L U D decomposition; I think it comes under the name ludcmp. I would recommend use of
this book to find the source codes. Once we have now obtained the solver which is the
L U decomposition based certainly for us, we know our stiffness matrix is going to be
banded in nature. There is going to be only a certain number of elements on either side
of the main diagonal which are going to be non zero and everything else is 0. In that
case, I can further reduce my cost of doing the elimination or finding the solution. We
can find in these books, the banded solvers. Essentially what happens, when we have this
kind of a property that the matrices sparse. If there are a lot of 0s in the matrix then
I can simply store only the part which is non zero. There are these solvers which take
advantage of that and the elimination also I will do using that information. So the cost
of computation then comes down. One can also use the banded solvers.
We obtained the solution that is we have obtained this degrees of freedom U which implies that
I should be able to now get the finite element solution, anywhere in the domain. Specifically,
in an element k, I can get the solutions in terms of the master element coordinates. uFE
at any point x in the element k is equal to sigma i is equal to 1 to P plus 1 which is
ndof for the element uig Ni hat where psi corresponds to the x that we taken. Now what
is ig? ig if we remember is nothing but ieldofs (k,
i), for the kth element the ith row. It is very easy to construct the solution in each
element. We can do piece by piece and also duFE by dx in the element k is equal to nothing
but sigma i is equal to 1 to P plus 1 uig dNi hat d psi, whole thing into d psi dx which
is 2 by hk. This will give me the derivative of the finite element solution at any point
in the element. These are the degrees of freedom which we have already solved for . This we
will obtain from the shape function routine. This I will get from the solver. This will
come from the shape Ni hat, give me any point.
As far as plotting is concerned, take the master element divide it into some number
of sub divisions. I say psi1 is equal to minus 1. This is the psi1, psi2, psi3 so on. This
becomes psiN plus 1 is equal to plus 1. I can find the size. At this size, I will find
the solution, by calling the shape function routine, telling it to give me the shape functions
in the derivatives and using that to construct uFE at the point psii. These I can output,
I can loop element by element that is, I go over the various elements in the domain - elements
1, 2, 3 and 4, for each element put this uniformly distributed set of points, output the value
of the displacement and the derivatives at this point. Let us say I get something like
this as the out. I can connect these by a line to get a plot of what the finite element
solution will look over the whole domain. As far as the plotting is concerned, I can
output the values of the finite element solution at these points connect them by a curve and
I am done. Similarly, I can do the same thing for the
derivative. We see that the uFE is very good at the nodes. What are the nodes? Nodes here
we mean as the extremities of the elements. These are our nodes. It is very good. What
do you mean by very good? It is very close to the exact value, if I had an access to
that at these points. There is a rate of conversion for these essentially for such problems that
we have taken, it goes as order h to the power of 2P that is the error between the finite
element solution and the exact solution at this point.
Further, if I have this then I can see that essentially, if I plot the difference between
the exact solutions, if I had access to it and the finite element solution in an element,
it will essentially look like this ; this is u minus uFE. Imagine that the size of the
element tends to 0 that is I am shrinking the size of the element. Then I can show that
this is essentially equal to some constant aP plus 2 NP plus 2 hat that is it is essentially
given by the next higher bubble functions. So I have taken all the bubble functions of
order P. Now this one corresponds to the bubble function of order
P plus 1. If we see that somehow from here I am not
doing it explicitly, essentially it can be written in terms of what we had given as the Legendre polynomials
of order P. If I take the derivative, this is all I am doing in the master element. If
I take the derivative of this finite element solution in the master element, the error
in the finite element solution, it will essentially be dominated by this term which is the term
corresponding to the P th order of Legendre polynomial. If this term has to be 0, what
do we mean by this term being 0? This is 0, d d psi of u minus uFE as hk tends to 0.
It tends to 0 at the points which are the roots of the Legendre polynomial of degree
P. It is as simple as this. Let us take an example. If I have taken linear approximation,
this derivative, if I have to find d dx of u minus uFE, I simply multiply it by 2 by
hk. That is not a big problem. Wherever, this is 0, the actual derivatives are also going
to be 0. If I take the linear approximation that is P is equal to 1, then the derivative
in the finite element solution should tend to be exact at the root of the first Legendre
polynomial. First Legendre polynomial is some constant into psi. So it should vanish at
psi is equal to 0. In fact that is true that the error in the derivative is 0, essentially
at the centre of the element. Similarly, if I go to P equal to 2, I will
see that the derivative vanishes at the 2 roots of P2. What are these roots of P1 and
P2? These are nothing but the Gauss points corresponding to the 1 point rule, the integration
points. This is at the points corresponding to the 2 point rule. I can essentially if
I have obtained the solution using a P order approximation then I go to the P order or
P point rule and find those points corresponding to the P point rule. At those points, I know
that the finite element solution, the derivative of it is going to be very good. These points
are called Super Convergence points. If I have obtained super convergence points,
I should, at the nodes that are at the extremities of the element, the derivative obtained from
the finite element solution is the worst that I can get. In fact, I should try not to use
that information . How can I use this information about these good values of the derivatives
of the finite element solution, to obtain a better approximation of the derivative by
some post processing? What we are actually doing? With the previous transparency also
we have done that thing. Here also we had said the plotting. We are actually in the
regime of the post processor. Once the solution is obtained, we are talking of how best to
interpret our data, to represent our data and how best to extract something better from
it. I would like here to get a better value of du star dx from the duFE dx. I know that
I can use the values of the finite element solution, the derivative of it, at these points
which correspond to the roots of the Legendre polynomial of order P which are nothing but
the P point integration rule.
If I have given these points, let me call these points as super convergence points,
super convergence of derivative of duFE dx at roots of the Legendre polynomial of order
P where P is the order of approximation which corresponds to P point integration rule. I
already have this information stored with me. Where are these points? Simply take these
points and I will them as points psii super convergence, i going from 0 to P. Given these
points, if I go to the elements by mapping back from the master to the physical element, I will have some points in each of
the elements where the derivative of the finite element solution is going to be good. If I
plot these derivatives in the elements, we will see they will do like this and there
will be a jump in the value here at the interface; derivative will do this.
Now this is something bad, I do not want this, so this is d dx of u minus uFE where u is
the exact solution. Here I am talking of the derivative of the error and we want this to
be essentially a straight line 0. I do not want the derivative of the error to be big
anywhere. Further at the interfaces, I get a jump in the derivative of the finite element
solution is the something which I do not want physically, because it leads to a jump in
the axial force that we obtained. How best can we now use this information about this
point to create better derivative information? The simplest possible way is the following.
Let us say this is the element k . This is its neighbor k minus 1 and k plus 1. Let us
say, here I am doing the one point rule but we can have the various things. This is corresponding
to x1 k minus 1 super convergence here. This is x1 k plus 1 super convergence here . These
are the mapped physical points corresponding to the super convergence points in these elements.
Let us say, I have taken linear approximation for the time being. At these points, I know
the derivative is going to be good. Let us say the derivative here is this, derivative
here is this, derivative here is this . To obtain d u star dx in element k, because I
have used linear approximation, I am going to now represent du star dx by a linear in
this element. So I am going to represent as a0 plus a1x in this element and I also extended
by the same thing in the neighbors that is the k minus 1 element and k plus 1.
I would like to get these coefficients such that I hope that it is going to be a better
derivative information in this element then we have out of d uFE dx. I would like to obtain
essentially using this information, the straight line fit and I hope that this straight line
fit is essentially equal to du dx. How do I construct this straight line fit? I take
for the element k, remember that I am doing it for the element k. I have to do it for
each element separately. Find these coefficients a0 a1 from there construct the value of the
variation of du star dx in the element k.
To do that again I take this neighborhood; it is k, k minus 1, k plus 1. I am going to
define this functional
as 1/2 sum over all the i(s), i is equal to 1 to 3 into (3 actually should be 3P of) duFE
dx evaluated at the point x for P is equal to 1, i s minus du star dx evaluated at the
point x1 i s whole squared. I am taking this term is given by a0 plus a1x. I am taking
this value of x from super convergence point for the neighbors, putting it here I am doing
discrete sum. This is the functional J and I am looking for the coefficients a0 and a1
which minimize. So find a0, a1 such that J is minimized . J, I would say is defined for
the element k. This definition, I have to do for any generic element k. I have to go
from the first element to the second one, to third one. I will find this coefficients
a0, a1 corresponding to the element. From there I can construct du star divided by dx
at the element level. What happens at the various places? If I have an element at the
boundary, that is the first element. In this case, I do not have the 0th element; I only
take this neighborhood. For the last element, again I will take the n minus 1th element
and the nth element. So this is nelem and this is nelem minus one . This is how I am
going to construct these so called recovered coefficients.
When I piece them together, I will see that essentially the du star dx it will look like
a piecewise linear. If I am at a material interface, if this is the element which has
this material interface, I am sitting at the material interface then I do not take the
neighbor. If this is the element k, I do not take the k plus 1 neighbor because it is lying
in a different material. I will only take the neighbors from the same material neighborhood.
This is how we can construct the recovered du star dx in each element k. This is called
the Patch Recovery method and this indeed gives us the values of the derivatives which
are very good. If we do this kind of a post-processing, we
can construct very good recovered derivatives without doing any extra computation. The computation
here is minimal and it is done at the element level. I will do for each element separately.
Out of what I did here, I will get the system of equations in terms of a0, a1 up to ap.
I will solve that system of equation for the element. I have to solve it for all the elements
separately one by one and I have constructed the du star dx information for all the elements.
This one said to be a Super Convergence Recovery method.
There is another version of it where I can actually modify J to have EA at the point
x1, i comma s. Here I have taken 1 but in our sum we should have all these integration
points in the three elements which form the neighborhood of the super convergence points.
Here also I will have EA at x1, i comma s. If I add this, this is another way, an alternative
definition for getting these super convergent recoveries. Especially when the EA is non
uniform, but not transitional one it is simply if it is a tapered beam for example, then
this will do the job.
We have done these recoveries -- certain things; once we have done this, our post-processing
has also been completed. After we have done all this, we have written a program, how do
we check our program? We check it using so-called Patch Tests. Patch tests or benchmark tests,
the basic idea in all these are that we take problems for which we can get the exact solution.
Then we go and try to solve this problem and we see whether our finite element solution
is close to the exact solution. It should not be too far away. If it is too far away,
there is some problem. In fact, in the patch tests, we can construct
at least for the case k0 equal to 0, EA is equal to constant. If we take our f(x) is
equal to 1 or x or x square, we know our exact solution in this case will be of the type
x squared by 2 is a polynomial . Here it will be x cube by 3 and here it will be x to the
power of 4 by 4 plus some constants. In these cases, when k0 equal to 0, is EA is equal
to C, taking either this as the load or this as the load or this as the load , correspondingly,
I will get the quadratic, the cubic and the fourth order solution. We can simply pay for
this solution. Let us say I have taken f(x) is equal to 1,
I will take P is equal to 2, the quadratic approximation and with this quadratic approximation
I should get the exact solution to the problem back. Similarly, if I take P is equal to 3
and take f(x) is equal to x, I should get uFE is equal to u exact and for P is equal
to 4, I do this. If somehow the code does not give the finite element solution equal
to the exact solution then there has to be a bug in whatever we have done. We have to
go back and debug our program and set it right so that till at least the machine precision
that is specified. Generally, we do double precision calculations. The numbers will match
with the exact solution. So with this, we have come to the end of the
one dimensional programming assignment. With all the tools given, we can develop the one
dimension code. I tried to give a feel of what will be coming or needed in the 2 and
3 dimensions. Essentially, the data structure that we have constructed here, the kind of
procedure we followed here, will carry over almost unaltered into the 2 and 3 dimensional
cases. It is for you to write the program, test it for this patch test or for benchmark
test for which you construct the exact solutions and then use it for other problems for which
you would like to get a solution for which getting a exact solution is difficult.
Thank you.