Tip:
Highlight text to annotate it
X
In the previous lecture, we had started talking about the two dimensional problem, where,
we had obtained the weak formulation for a given partial differential equation which
arises in the two dimensional case.
We had considered the generalized problem of steady state heat conduction and from there
we went ahead and used all the tools which we had developed for the one dimensional problem
to obtain the variation formulation for the two-D problem. We saw what are the major differences
between the one-dimensional problem and the two-dimensional problem.
Let me recap. The major differences are that here we are talking of a two dimensional domain,
so the geometry of the domain now plays a role and how well we resolve the geometry
of the domain is another issue that has to be handled. For example, I mesh my domain
with triangles. This is my domain. I know that these triangles with straight edges will
only approximate this boundary curve. So I have to ensure that this approximation of
the boundary curve is decent. That is, I have a sufficiently refined mesh, that is, there
is sufficiently large number of small triangles at the boundary so that these curves are taken
care of. We will also see other approaches of exactly representing these boundary curves
wherever it is possible. The geometry was the problem and this leads to the discretisation
error as far as the definition of geometry is concerned. Next thing was that in application
of the boundary conditions in the one-D case, we had talked about the boundaries being two
points, the end points of the domain. Here, the boundary could be a line, a curve, or
a contour. So the boundary conditions have to be now imposed on a close contour. So this
is another difference between the one dimensional and two dimensional cases.
Another thing we had discussed is that if I have a domain like this, you see that because
of this corner here in the domain, there is going to be a singularity in the solution
that you obtain. Irrespective of what the loading is, I am going to get, depending on
the boundary condition on these two edges, I am going to get a singularity in the solution
that you obtain in the vicinity of the corner. So these are certain things we have to keep
in mind. Now let us go ahead and build the finite element approximation for the model
problem we have taken. For the sake of simplicity, here I am going to consider the simplified
version of that model problem, that is, del to the power of 2 u by del x to the power
of 2 minus del to the power of 2 u by del y to the power of 2 is equal to r, which is
what we know as the parson equation. I have deliberately chosen, with respect to our earlier
model problem, a11=1, a22=1, and a12=0. Two is just to fix ideas. Everything else we will
follow. We will take a rectangular domain. In the rectangular domain, we had said, we
are going to make a mesh of triangles. This is the simplest possible domain that I can
have. Let us say it is a square domain not even a rectangular domain. It is a square
domain with sides 'a' and 'a' and I have made these triangles.
What we had said is, the first thing that we are going to do is construct the crudest
possible approximation; the crudest possible approximation was an extension of the linear
approximation that we had in the one-dimensional case. So we would now like to create a linear
approximation in the two dimensional case. The linear approximation in the two dimensional
case would be of this form now. So, the linear in the two-D case will be a0+a1x+a2y and in
this linear approximation, what we are going to do is, we are going to construct the basis
functions over these elements that we have, over this domain, such that, I have a piece-wise
linear approximation in each element, like we had in the one-D case. Over each element
here, I would like to construct a piece-wise linear approximation. Piece-wise linear approximation
means that over the element, my basis functions, the truncation of the basis functions are
linear polynomials. Now how many independent functions do we need in an element in order
to completely define a piece wise linear? The answer is that, since this piece wise
linear incorporates three unknowns, we need three independent functions. In the two-D
case, we required two independent functions in an element. Now we need three independent
functions. How do we define these three independent functions?
We had said that we are going to define these three independent functions with respect to
the three corner vertices of the element or the three corner nodes of the element. These
are the nodes of the element of interest. In the two-D case, I have the triangles as
the elements and the corners of the triangles are the nodes. So we are going to define these
piece-wise linears with respect to these elements. Let us now go and see how we do this job.
What we say is that, let the polynomial, let any function 'u' be represented by a linear
in the element of interest.
So my element would be this and as we had said, these are the corner vertices of the
element. So this is my element top and I give it some name. So what we want, as we had done
in the one-D case is, to interpolate the given function u with a linear in this element.
In order to interpolate this function u with a linear in this element, what we do is: u
is given by some linear which is c1+c2+c3y in this element. Then what is a value of this
linear at these nodes of the element? The nodes of the element have coordinates (x1,
y1), (x2, y2) and (x3, y3). So, we want a representation of the function that we have
taken at these nodes. So the value of the function at these nodes is u at (x1, y1) u
at (x2,y2) and u at (x3 ,y3). What we want now, going by what we had done in the one-D
case is, the linear interpolation of this function to match the value of this function
at the three nodes. When we have this linear interpolation, we want c1 plus c2 x1 plus
c3 y1 to be equal to the value of the function at the point x1 y1. Similarly, I want c1+c2
x2+c3y3 to be the value of the function at the point (x2, y2) and at the third point.
These values of the function at these three points I am going to call as mu1, mu2, and
mu3. I assume that I know my function, I know the values mu1, mu2 and mu3, which can be
written in terms of the coordinates of the three points and these unknown constants.
We want to find these unknown constants c1 c2 and c3. So the three values of u are given
as a combination of these unknown constants c1 c2 c3. Now the job is simple: it is a three
by three matrix, so we have to invert this matrix to get c1 c2 c3 in terms of mu1, mu2
and mu3. That is about all we have to do. What we end up getting is: I call this matrix
as matrix A, so A inverse will be one by twice the area of the triangle, we have called it
tau, into some coefficients alpha1, alpha2, alpha3, beta1, beta2, beta3, gamma1, gamma2,
gamma3. We can write the area of the triangle A tau in terms of these alpha1, alpha2, alpha3.
What do we get then? The unknown constants c become as A inverse into this vector mu
and after expanding this out, the constant c1 is 1 by 2 area of the element tau into
alpha1, mu1 (mu1, is a value of the function mu at the first node which has coordinates
(x1,y1)) plus alpha2, mu2, plus alpha3 mu3. Similarly, I can find c2 and c3, where these
alpha 'i's are given like this in terms of the three nodal coordinates, that is, alpha1
will be x2 y3 minus x3 y2, where, i j k are written in a cyclic permutation. Similarly,
alpha2 will be equal to x3 y1-x1 y3 and so on. Similarly, I can find what is beta1 and
what is gammaI in terms of the x1 x2 x3, y1 y2 y3. Once I have these quantities, I can
rewrite: the linear interpolation of u over the element tau in terms of the values of
the function at the three points u1 u2 u3 as one by twice area of the element tau into
alpha1 u1 plus alpha2 u2 plus alpha3 u3 plus the part multiplying x plus the part multiplying
y.
So this you can say is the part multiplying one, the part multiplying x and the part multiplying
y. So what is the next thing? We would like to group all the coefficients or all the monomials
multiplying u1 together, all the monomials multiplying u2 together, all the monomials
multiplying u3 together. Rewrite u in the element tau in terms of the values of the function at the three nodes
into some polynomial Psii. Psii tau is equivalent to what we know as element shape function
Ni tau. It turns out that Psii tau is given as one by area of the triangle into alphai
tau plus betai tau x plus gammai tau y. These Psii taus are the shape functions in the element.
How many do we have? We have three of them. So we have the three linear shape functions
in the element obtained by doing this procedure. You see that each of the shape functions is
a linear showing that this set is a complete set, that is, it can represent any linear
polynomial exactly by taking a linear combination of these and showing that these are linearly
independent is now relatively an easy job. One can do exactly the same thing that we
had done in the two dimensional case. Let us go further; so what will these functions
look like?
Let me just take an example. Let us say this is node 1, this is node 2 of the element and
this is node 3 of the element. Here I would like to plot the function psi or N1. If I
take N1, as you will see, the N1 has a value one at the node one. So, we will see N1 tau
at x1 tau y1 tau is equal to one. It will have a value one at the first node and it
will have a value zero at the second node and the third node. So this function is going
to be like this. It is like a slanting roof. It has a value one at the first node, zero
at the other two nodes. That is, it is zero along this whole line. Similarly, this is
function N1 tau or as I have written in the previous slide as psi1 tau. Similarly, if
I want to draw N2 tau, N2 tau will be a function, which has value one at the second node and
value zero at the first node and the third node. So this is my function N2 tau. So it
is a roof slanting down from the second and similarly I can draw the third function. These
are our element shape functions and you can see that these shape functions for the ith
node will correspond to a tent like structure. That is, by piecing together all the shape
functions, I get the global basis function phii for the ith node of the mesh that we
have and this function vanishes on these edges. It vanishes on this black edges and it is
one at the node with respect to which it is defined and it is piece-wise linear and each
of the elements shares that node. This, you will remember again, that in the two-dimensional
case, a node may be shared by many more than two elements. In this figure that I have drawn
there are four elements sharing a node. This basis function has to be defined with respect
to all the elements which share that particular node. With this, now, we have obtained our
representation of the linear shape functions. Let us now construct the finite element solution,
if I have to do that in an element. I define the finite element solution globally. So,
given this phii, let me again make a mesh. For our own convenience, let us take an example
of a mesh of eight elements. You can take the simplest possible method. So let us say,
this is my node 1, this is my node 2, I am doing a numbering of the nodes in the mesh,
node 3, node 4, node 5, node 6, node 7, node 8, and node 9.
Similarly, I will number the elements, this is element 1, this is element 2, element 3,
element 4 and elements 5, 6, 7, and 8. So, I have numbered the nodes and the elements.
Now how many basis functions will I have formed by piecing together these piece-wise linears?
I will have nine basis functions corresponding to the 9 nodes of the mesh. If I want to make
a finite element solution, which is a linear approximation over this domain, then, uFE
would be equal to sum of i is equal to 1 to 9 ui phii, which is the function of (x, y).
Now, I want to look at the restriction or the part of this finite element solution over
some generic element. Let us say I am looking at it over this third element. So, uFE over
a generic element tau will be given in terms of the value of this coefficients ui at the
three nodes, which are the N vertices of this particular element, multiplied by the restriction
of these basis functions over this element. So, it can be written as, sigma i is equal
to 1 to 3 ui in the element tau Ni tau (x, y). By this we understand that if I am in
the third element, I will have a local numbering for the nodes. If you recollect, in the one-D
case too we had defined the global numbers and the local numbers. So this will be my local node 1, local node 2 and local node
3 for the element. Then I see that the u1 of the element tau is equivalent to the global
u2, u2 of the element tau is equivalent to the global u3, u3 of the element tau, which
is the third element, is equivalent to the global u6.
For every element, I have to obtain the local to global numbering. So remember that same
data structures have to be continued here. Local to global enumeration has to be done,
which tells me which global degree of freedom does my local degree of freedom correspond
to. Similarly, the N1 tau for this element is equivalent to the phi2 in this element
N2 tau is equivalent to the phi3 and in this element N3 tau is equivalent to the phi6.
We now know the one to one correspondence between what is the global solution and what
is that we are doing locally. Once I have the representation of the finite
element solution in the element then I can go and do the next step. That is, I can do
my element calculations.
What do we do for the element calculations? We should remember, if we take the weak form
that we obtained and put it for this particular model problem that we have taken, the weak
form would be del u del x del v del x plus del u del y del v del y dA would be equal
to integral over the domain r v dA plus if you remember we had talked about the Norman
part of the boundary, the part were the flux conditions are specified, you will get g v
d s. g is the specified boundary flux and v is our test function or the weight function
that we have used. This is going to be our weak form, which is
obtained as a specialization of the general weak form we had obtained earlier. Now, what
do we know from this? We can now partition as sum over all the elements: sum 1 to number
of elements integral over the elements tau of del u del x del v del x plus del u del
y del v del y dA. Similarly, we can do the same for the right hand side. Now from the
finite element point of view, we replace u with uFE . uFE on the element is given by
uFE tau. Similarly, what we do as far as v is concerned is, we take v to be the global
basis function itself. So v would be again given as v over tau. Here we are talking of
the global basis functions that are going to be non-zero in this element or the global
basis functions that correspond to the Ni tau, that is, the shape functions of the element.
Essentially, what happens is, as far as the contribution of the element is concerned,
the element is going to only contribute to the rows in the global stiffness matrix or
in the global equations, which correspond to the phiis, which are non-zero in the element.
Let us go back to our previous figure.
If I am in element three, then, corresponding to the uFE that we have constructed, I will
have a 9 by 9 global system, in terms of the nine unknown uis. So, as far as element 3
is concerned, phi2, phi3, and phi6 are non-zero in element three. All other phiis are zero.
So my element 3 is only going to contribute to the global equations corresponding to phi2,
phi3, and phi6, that is, it is going to contribute to the second, third and the sixth global
equation and similarly for the other equations. It is also going to contribute to the second
row to the third row and the sixth row, because these are the only coefficients or the basis
functions, which are active in this element. With that understanding, for the element,
there are only three active basis functions. Even the uFE of the element is in terms of
these three coefficients. The element stiffness matrix will be a 3 by 3 stiffness matrix as
far as the piece-wise linear approximation is concerned. This 3 by 3 stiffness matrix
for the element is what we would like to obtain. Once I obtain this 3 by 3 stiffness matrix,
then we do the usual process of assembly to get my global system. What does this 3 by
3 stiffness matrix correspond to?
Just like we have done earlier, let us talk of the integral over the element del uFE over
the element tau del x (now here instead of the v, we are taking the basis functions which
are non-zero in this element, which are nothing but) del Ni in the element tau over del x
plus del uFE in the element tau del y del Ni in the element tau del y. These are essentially
corresponding to i=1, 2, 3. We will get the three equations from the element. In the three
equations, when I substitute for the uFE tau, I will get the three columns, which are there
in terms of u1 tau, u2 tau and u3 tau. If I do this, I will get the stiffness matrix
entry in the element tau by expanding uFE tau and is given as integral over tau del
Nj tau del x del Ni tau del x plus del Nj tau del y del Ni tau del y. So this is going
to be my ijth term of the element stiffness matrix, where, ij is equal to 1, 2 and 3.
Everything that we do is exactly the way we had done for one-D. Only thing is, the number
of entries increases here and we are dealing with partial derivatives of the quantities
of interest. If I have to do this for the standard model problem that we had taken,
that is, the generalized one, then we will have to do it with the a11, a12 and a22 sitting
there. The stiffness matrix entries for the element can now be computed in terms of the
derivatives of the element shape functions Ni and Nj, where Ni and Nj have been obtained
from our interpolation that we have done earlier. Now, if I take the piece-wise linears, we
see we see that these derivatives are going to be constant. As far as the piece- wise
linears are concerned and the model problem that we have taken, the integrant here is
going to be a constant integrant. This constant we can easily obtain corresponding to each
i and j. The integral of the constant against the area will be equal to the area times this
constant. As far as this part is concerned, finding
the integral is not difficult at all and one can do it very easily and explicitly. However,
let us go and build some of the other features, which are going to be important from the point
of view of a computational implementation. Let us first talk of
the geometry representation. The geometry is going to be represented by x being given
in terms of, that is, in the interior of an element tau, if I want to find any point x,
it will be given in terms of the coordinates of the end points of the element. It will
be given in terms of, as a linear interpolation, x1 tau N1 tau plus x2 tau N2 tau plus x3 tau
N3 tau. Similarly, y is equal to y1 tau N1 tau plus y2 tau N2 tau plus y3 tau N3 tau.
Now, if you remember, from a computational point of view, it was much easier to work.
We could do everything in the physical element, but we decided to map our physical domain
physical element in a one-D case to a master element. Similarly, in the case of this two
dimensional problem, let us map our physical element into a master element. Let us define
the mapping as a linear map.
So what happens? I have my physical element, which has these nodes x1 tau y1 tau, x2 tau
y2 tau, and x3 tau y3 tau. This, we are going to map to a master element given in terms of the coordinate system psi
eta and defined as this or master element is going to be such that the first node has
coordinates (0,0). The nodes of the master element I am numbering them with hat. The
second node has coordinate (1, 0) and the third node has coordinate (0, 1). So what
have I done? I have simply mapped my original physical triangle linearly to this master
triangle. That is, the straight edges go to straight edges and the area maps linearly.
I have done the linear map and now my next job is to define this map. How do I define
this map? In order to define this map, let us take the linear shape functions and map
them to the shape functions that we have, that is, the corresponding entities over the
master element. I will call this master element tau hat. This is my element tau. So how do
I do this mapping? A linear polynomial in the physical element goes to a linear polynomial
in the master element. If a linear shape function was one at a physical node, it will be also
one at the corresponding master node. If it was zero at a physical node, it will be zero
at the corresponding master node and so on. If I have the N1 tau, the N1 tau in the master
element will become N1 hat. This N1 hat will be a function of psi and eta.
Now, what is this N1 as a function of psi and eta? N1 hat as a function of psi and eta
is quite easy to show, it is a function which is one here and zero along this line because
remember that N1 tau was zero along this line, so similarly N1 hat has to be zero along this
line. This line has an equation psi plus eta equal to one. If it is zero along this line,
the function has to be one minus psi minus eta. Check that this will be one at this point
where it is (0, 0) and along these two points - two and three.
This is how I have defined my N1 hat. Similarly, I can define N2 hat as a function of psi and
eta is equal to psi. Again, the function has a value one at the second point and value
zero at the first and the third point. This becomes psi, one can check. Similarly, N3
hat as a function of psi and eta will be equal to eta. By doing this mapping, I have at least
in principle constructed the corresponding map versions of the shape functions. How can
I define the geometry mapping?
Geometry mapping becomes x is equal to x1 of tau into N1 hat plus x2 of tau into N2
hat plus x3 of tau into N3 hat and y becomes y1 of tau into N1 hat plus y2 of tau into
N2 hat plus y3 of tau into N3 hat. This is how I can give the x and y at any point psi
and eta in the master element, the corresponding value of x and y.
Once we have done this mapping, remember that we have to do master calculations. In the
master element, if you remember, we had an integral of this type. We would like to convert
it to an integral over the master elements. So from the physical I come to the master
where I have the integral over the master element tau hat of this integrant. This integrant
converted into an expression in terms of psi and eta into (note that the area in the physical
element has to be mapped to the area in the master element) and that will be in terms
of something called the Jacobian into dA hat. The question is, what is this and how do I
use the derivatives in the master element to obtain the derivatives in the physical
element? For that we will have to obtain the metrics of the transformation. So how do we
obtain the matrix of the transformation?
If I say that I want dx, dx will be equal to del x, because x is now a function of psi
and eta, del x del psi d psi plus del x del eta d eta. Similarly, dy is equal to del y
del psi d psi plus del y del eta d eta. This is because x and y are functions of psi and
eta. So, change in x is again in terms of change in psi and eta. Now we have to obtain
these quantities. Where are these quantities going to come from? These are going to come
from the definitions of x and y that we have given in terms of the three nodal xs and ys
in the physical element and the definition of the shape functions in the master elements.
One can show from the expression we have there that del x del psi is equal to x2 of tau minus
x1 of tau. Let us go back and see. Here, if you remember, N1 hat was one minus psi minus
eta, N2 hat is psi, N3 hat is eta. If I take this definition and go back, dx
d psi will be equal to x2 minus x1. Similarly, del x del eta is equal to x3 tau minus x1
tau and del y del psi is y2 tau minus y1 tau del y del eta is equal to y3 tau minus y1
tau. Let us now rearrange this thing in a matrix form. So I can write it as, dx dy is
equal to del x del psi , del x del eta , del y del psi, del y del eta, d psi d eta. I have
simply rewritten this expression in a matrix form. Once I have this expression in the matrix
form, then I can talk of the Jacobian.
The Jacobian is nothing but the determinant of this expression of this matrix that we
have written. And we would like to obtain expressions for del psi del x, del psi del
y, del eta del x, del eta del y. That is, since we can write x and y in terms of psi
and eta, we can also write psi and eta in terms of x and y. So we want metrics of the
inverse transformation. How do I do it? By exactly following what we have done there,
we can get d psi d eta is equal to del psi del x, del psi del y, del eta del x, del eta
del y dx dy and this matrix is nothing but the inverse of what I will call as the Jacobian
matrix here. We take the inverse of this, bring it on this side and we will get d psi
d eta in terms of dx dy. So, all we have to do is, given this matrix
J, we need to find the inverse of that matrix. Now what is going to be the inverse of this matrix,
say J inverse?
J inverse is equal to one by Jacobian, we interchange the diagonal entries and so I
will get del y del eta, del x del psi and take the negative of the signs of the half
diagonal entries and so minus del x del eta , minus del y del psi. This is my J inverse
and the J inverse is now going to give me my del psi del x and so on. Del psi del x
is equal to one by Jacobian into del y del eta, del psi del y is equal to minus one by
Jacobian into del x del eta. del eta del x is going to be minus one by Jacobian into
del y del psi and similarly, del eta del y is equal to one by Jacobian into del x del
psi. Obtain the matrix of the inverse transformation. Why do I need it? Because, if you remember,
these quantities are equivalent to saying del Ni hat del psi del psi del x plus del
Ni hat del eta del eta del x. Because I can rewrite my Ni tau as a function of psi and
eta, that is, Ni hat as a function of psi and eta, then the x derivative of Ni tau is
given by this chain rule that del Ni del psi del psi del x del Ni del eta del eta del x.
We now need to obtain these quantities. But these quantities are quite easy because we
have already obtained each one of these quantities from the definition of x and y.
Given these quantities, I can go ahead and now define what is del psi del x, del eta
del x. In this case, you notice that for this triangle, for linear mapping, del y del eta
del x del eta del y del psi del x del psi are all constants. Why? Because x and y was
linear in terms of psi and eta. Take the derivative of the linear and you will get constant. Because
they are all constant, the Jacobian is also a constant. This is true only for the linear
maps of triangles. We will see that when we go to quadrilaterals, Jacobian will no longer
be a constant. So, given these entities now we can plug them in and obtain all the quantities.
For the element, I have to input the three nodal coordinates the x1 x2 x3, y1 y2 y3 and
I can construct the metrics of the transformation. Let us now take a very simple example of the
same domain that we had taken. Let us assume that my elements have edges of size h.
Given these elements of size h, I would like to construct, let us say, for the first element
in my domain, the entries of the stiffness matrix using linear approximation.
In the next class, we are going to construct the stiffness matrix of the elements and then
we will talk a little bit about the assembly for this special problem. Another important
thing that we are going to talk about in the next class is impassing, how to apply the
natural boundary conditions, that is, the force boundary conditions. We will also finally
discuss how to apply the displacement boundary conditions. Once we do that, then we are not
going to be happy with linear approximations in the null element. We would like to extend
it to higher order approximations, quadratic, cubic, fourth order and we will see how we
can construct the basis functions or the shape functions in an element, corresponding to
any order of approximation that we desire and then how to put it in the framework of
whatever we have developed here as far as the element calculations are concerned.