File: Create a C Sharp file called "fem1.cs", which will contain your 1D finite element program.
First and last message: To start, the program should print:
(date and time)
FEM1
C Sharp version
One dimensional finite element code.
and at the end of execution the message
FEM1
Normal end of execution.
(date and time)
Problem parameters: The first part of the program should define some parameters and print their values. To begin with, define the following scalar parameters:
EN, integer, number of elements.
EO, integer, order of each element
XL, real, the left endpoint;
XR, real, the right endpoint;
NN, the number of nodes.
and to begin with, set these to have the values
EN = 8
EO = 1
XL = 0.0
XR = 1.0
NN = EN * EO + 1
Include statements that print the name and value of each parameter.
Node coordinates: Create a real array called XN. Give it size NN. Set the entries of XN so that the first value is XL, the last value is XR, and the intermediate values are evenly spaced. After you set XN, print out the index and value of each element of XN:
XN array:
0 0.0
1 0.125
2 0.250 ...
Element connectivity: Create an integer two-dimensional array called EC. The array EC should have EN rows and EO+1 columns. The first row of of EC should contain the indices of the EO+1 nodes that are "connected" to element 0, that is, nodes 0 and 1. The second row should contain the nodes connected to element 1, that is, nodes 1 and 2, and so on. After you set EC, print out the element index and the nodes:
EC array:
0 0 1
1 1 2
2 2 3...
Program run: Run your program and save the output, which should now consist of
Starting message
parameter values
XN values
EC values
Ending message
Bring a PRINTED copy of your program.
Bring a PRINTED copy of the output of your program.
For our programs, a set of basis functions will be defined. The basis functions we will use are called nodal basis functions. A nodal basis function is a function whose value is 1 at one node and zero at all the other nodes. The nodal basis function will be identified by the index of the node where it is nonzero. If we use the generic name Phi for any basis function, then a specific basis function would have a name like Phii(x). This function is 1 at node XN(i) and zero at all other nodes. Since there are exactly NN nodes, there will be exactly NN nodal basis functions.
Piecewise linear basis functions: One of the most simplest ways to define a set of nodal basis functions is to use the piecewise linear definition. In that case, the function Phii(x) is defined as follows:
Formulas: What is a mathematical formula for the value V of basis function Phii(x)?
Piecewise Linear Basis Functions: We need a procedure, which we can call pwl(), to evaluate finite element basis functions of piecewise linear type. The information that this procedure will need is:
X, the evaluation point.
I, the index of the basis function, which must be between 0 and NN-1.
NN, the number of nodes.
XN, the node coordinates.
The function should return V, the value of the basis function.
Write a function of the form
double pwl ( double x, int i, int nn, double xn[] )
which evaluates the I-th piecewise linear basis function at X.
Piecewise Linear Basis Function Plot: We want to include pwl() in the fem1 file, but before we do, let's test it, and practice using gnuplot. Create a calling program that defines NN = 5 nodes at XN = [ 0.0, 1.0, 3.0, 3.5, 10.0 ]. Then evaluate the basis function associated with node 3 (that is, the one at 3.5), at evenly spaced points between 0 and 10, and create a plot using gnuplot.
Bring a PRINTED copy of your pwl() function. Bring a PRINTED copy of the plot.
When we write out the finite element equations, it turns out that we will need to evaluate the derivatives of the basis functions as well. Since we plan to start with piecewise linear basis functions, the derivatives are simple. The derivative function dPhii(x)/dx is "going up" (derivative is a positive constant) in the interval [x(i-1),x(i)], "going down" (derivative is a negative constant) in the interval [x(i),x(i-1)], and zero elsewhere. Technically, we don't exactly know how to define a value for the derivative at x(i-1), x(i), or x(i+1), but that's OK because piecewise linear functions don't have a derivative at breakpoints. If we need to, we can set the derivative to zero at these points.
Formulas: We need a formula for the derivative. But that's easy, because we just need to differentiate the formulas we got for Phii(x), and work out the nonzero constant value in the two intervals. What are these values?
A procedure: Make a new procedure to evaluate the derivative, called pwld(). Simply make a copy of pwl(), because we will need the same input quantities. The formulas used in pwld() are simply the derivatives of the formulas used in pwl(). Now we should have a new function:
double pwld ( double x, int i, int nn, double xn[] )
which evaluates the derivative of the I-th piecewise linear basis
function at X.
Plot basis function and derivative: Create a calling program called "pwld_test" that defines NN = 5 nodes at XN = [ 0.0, 1.0, 3.0, 3.5, 10.0 ]. Have plwl_test call pwl() and pwld() to evaluate the basis function and derivative associated with node 3 (that is, the one at 3.5), at evenly spaced points between 0 and 10, and plot the function and its derivative together, on one plot, using gnuplot.
Bring a PRINTED copy of your pwld() function. Bring a PRINTED copy of the plot.
The mathematical form of the finite element equation involves many integrals. For the Poisson equation,
-u''(x) = f(x) for a < x < b
the finite element method asks us to evaluate integrals such as
Integral ( a < x < b ) phi'(i,x) * phi'(j,x) dx
and
Integral ( a < x < b ) phi(i,x) * f(x) dx
Although we might have some idea how to handle the first integral, the
second involves a function f(x) that could be fairly complicated. We
prefer to handle all the integrals using a numerical method known as
quadrature.
A quadrature rule typically describes how to estimate the integral of some function f(xi) over some standard interval such as [-1,+1]. The quadrature rule tells us to evaluate the function at several points, multiply the values by certain weights, and take the sum as the integral estimate:
Integral ( -1 < xi < +1 ) f(xi) d xi approx Sum ( 1 <= iq <= nq ) w(iq) * f(xi(iq))
Our integration intervals don't run from -1 to +1, so we need to make an adjustment to the interval [a,b]. First, we note that as the variable xi runs from -1 to +1, we can make the variable x run from a to b using the formula:
x(xi) = ( a + b + ( b - a ) * xi ) / 2
We also need to adjust for the fact that the quadrature rule is defined
on an interval of width 2, but our subinterval only has width
( b - a ).
Our quadrature rule adjusted for the interval [a,b] is
Integral ( a < x < b ) f(x) dx approx (b-a)/2 * Sum ( 1 <= iq <= nq ) w(iq) * f((a+b+(b-a)*xi(iq))/2)
Three quadrature rules for [-1,+1]:
The integral of f(x) = 1/x over [1,4] is ln(4). Write a program which estimates this integral using quadrature rules 1, 2, and 3. For each rule, print the estimate and the error.
Actually, in finite elements, our finite element mesh has divided [a,b] into subintervals of the form [xn(k),xn(k+1)]. Thus,
Integral ( a < x < b ) f(x) dx
= Sum ( 1 <= k < NN-1 ) Integral ( xn(k) < x < xn(k+1) ) f(x) dx )
So it will be convenient to be able to evaluate integrals over
[xn(k),xn(k+1)]. To do this, we just replace a by xn(k), and b
by xn(k+1):
x(iq) = (xn(k)+xn(k+1)+(xn(k+1)-xn(k))*xi(iq))/2
Integral ( xn(k) < x < xn(k+1) ) f(x) dx approx (xn(k+1)-xn(k))/2 * Sum ( 1 <= iq <= nq ) w(iq) * f(x(iq))
Now the formula seems very complicated, but the idea is still quite simple. To compute an estimate q for the integral of f(x) over [a,b], using n intervals, we could do something like this:
q = 0.0
for k = 1 to n
q = q + integral estimate for integral of f(x) over interval k
Write a second program to estimate the integral of 1/x over [1,4], using n equal subintervals. Break the integral up into n subintegrals, and use the quadrature rules on each subinterval. For the values n = 1, 2, 4, 8, 16. and rules 1, 2 and 3, print the integral estimates and errors.
Bring a PRINTED copy of your two programs and output.
Here is a standard example problem for the finite element method:
Over the interval:
a <= x <= b,
solve the partial differential equation:
- u'' = x
with boundary conditions:
u(a) = ua, u(b) = ub.
Suppose we break the interval [a,b] up, using 5 nodes, and define piecewise linear basis functions b1(x) through b5(x). We consider multiplying the PDE by b2, b3 or b4, integrating, and applying integration by parts. We assume that u is a linear combination of the basis functions, with unknown coefficients c. Then we arrive at the following linear system:
Row 1: 1 0 0 0 0
Row 2: Int b2' b1' Int b2' b2' Int b2' b3' 0 0
A = Row 3: 0 Int b3' b2' Int b3' b3' Int b3' b4' 0
Row 4: 0 0 Int b4' b3' Int b4' b4' Int b4' b5'
Row 5: 0 0 0 0 1
with unknown vector c and right hand side f defined by
c1 ua
c2 Int b2 * x
c = c3 f = Int b3 * x
c4 Int b4 * x
c5 ub
We can approximate all those integrals, because we can use quadrature.
Now we need to solve the linear system A*c=f.
In many cases, we can use an iterative method. If we use a direct method, we can take advantage of special features of the matrix, such as sparseness, bandedness, symmetry, or positive definiteness. Instead, we will use the basic method of Gauss elimination.
Gauss elimination should be familiar to you already. There are several variations. In our version, we have an NxN matrix A, an N vector F, and we seek the N vector solution C. We do this in two steps, forward elimination, and backward substitution.
In the following description, we assume that, if F is an N vector, its first element is F[0] and its last element is F[N-1]. Similarly, if A is an MxN matrix, its first element is A[0,0] and last element is A[M-1,N-1].
Forward elimination:
Start with column J = 0, and go to J = N - 2.
Start with row I = J, and go to I = N - 1.
Remember the row K where the largest element |A[I,J]| occurs.
Switch rows I and K of the matrix A.
Switch entries I and K of the right hand side F.
Start with row I = J + 1, and go to I = N - 1.
Let T = - A[I,J] / A[J,J]
Add T*row J to row I. This zeros out A[I,J].
Add T*F[J] to F[I]
At the end of forward elimination, A should be upper triangular,
and F should have changed.
Backward substitution:
Start with row I = N - 1, and go backwards to I = 0.
C[I] = F[I]
Start with column J = I + 1, and go to J = N - 1.
C[I] = C[I] - A[I,J] * C[J]
C[I] = C[I] / A[I,I]
C should now contain the solution of the linear system.
A Test Linear System: Here is a simple linear system, whose solution is known, and which can be used as a first check that your linear solver is correct:
------A------ * --C-- = --F--
1 0 0 0 0 1 1
-1 2 -1 0 0 2 0
0 -1 2 -1 0 * 3 = 0
0 0 -1 2 -1 4 0
0 0 0 0 1 5 5
Give your linear solver the matrix A and right hand side F, and see if
it can compute the correct value of C.
Bring a PRINTED copy of your programs and its output.
Consider the following problem:
Over the interval:
0 <= x <= 1,
solve the partial differential equation:
- u'' = -16
with boundary conditions:
u(0.0) = 3, u(1.0) = 1.
with exact solution
u(x) = 8x^2 - 10x + 3
For our program, we have XL = 0.0, XR = 1.0, UL = 3.0, UR = 1.0.
Divide the interval [XL,XR] into EN = 4 elements. Let the element order EO = 1 (piecewise linears.) Create an EC array that lists the pairs of nodes associated with each element.
Create a mesh of points XN, using NN=EN*EO+1=5 nodes, over the interval [0,1].
Create an NNxNN = 5x5 matrix A, whose entries are:
1 0 0 0 0
Int P2' P1' Int P2' P2' Int P2' P3' 0 0
0 Int P3' P2' Int P3' P3' Int P3' P4' 0
0 0 Int P4' P3' Int P4' P4' Int P4' P5'
0 0 0 0 1
where the expression "Int P2' P1'" represents the integral,
from XL to XR, of the derivative of basis function 2 times
the derivative of basis function 1. (Here I am labeling the
basis functions starting at 1, rather than 0.)
Create an NN = 5 vector F, whose entries are:
UL
Int P2 * 16
Int P3 * 16
Int P4 * 16
UR
Each matrix entry is an integral, which must be approximated by quadrature. Although each integral is defined over the entire interval, most of the integrands will be zero except over one or two subintervals. In our example, there will be two patterns, the diagonal and off-diagonal patterns.
Each diagonal entry is an integral of a basis function derivative times itself, for example: Int P3' P3'. The P3 basis function and its derivative are nonzero over two subintervals: [X2,X3] and [X3,X4]. It makes sense to apply a quadrature rule to each subinterval separately, and then add the results to get our approximate entry.
Each subdiagonal entry is an integral of a basis function derivative times a neighboring basis function derivative, for example: Int P2' P1'. There is only one subinterval in which the P1 and P2 basis functions are nonzero, namely [X1,X2]. Thus, to approximate an off-diagonal entry, we need to identify the single subinterval of interest and apply a quadrature rule.
The obvious way to fill in the matrix A is to do it one matrix entry at a time. Entry A(3,2) is the integral Int P3' P2', for instance, and so just from the index we know what to do.
An alternate way to fill in the matrix is to look at one element at a time, set up the quadrature rule there, and compute all the partial integrals we can. To see what this means, let's imagine that the two boundary condition rows of A have been set up, but the rest of the matrix is zeroed out. We begin by processing subinterval 1. The following are the only integrals that can be nonzero over subinterval [X1,X2]:
Int(X1,X2) P1' P1' <-- don't need, because A(1,1) is used for BC.
Int(X1,X2) P1' P2' <-- This should be added to A(1,2)
Int(X1,X2) P2' P1' <-- This should be added to A(2,1).
Since it's equal to Int(X1,X2) P1' P2', we
have already computed this value.
Int(X1,X2) P1 f(x) <-- Add to right hand side F(1)
Int(X1,X2) P2 f(x) <-- Add to right hand side F(2)
Now we move to subinterval [X2,X3]:
Int(X2,X3) P2' P1' <-- add to A(1,2) and to A(2,1)
Int(X2,X3) P2' P2' <-- This should be added to A(2,2)
Int(X2,X3) P2' P3' <-- add to A(2,3) and to A(3,2)
Int(X2,X3) P2 f(x) <-- Add to right hand side F(2)
Int(X2,X3) P3 f(x) <-- Add to right hand side F(3)
By doing the integrals this way, we only have to set up the quadrature
rule once in each interval, and we can see all the cases where
computing the integral of P2' P3' (which is added to A(2,3)) means we
have also computed the integral of P3' P2' (which is added to A(3,2)).
Creating the linear system is called "assembling" the system. Because the finite element method was originally used in structural mechanics, the matrix itself is often called the "stiffness" matrix.
Solve the linear system A * C = F.
Create a plot of the (approximate) solution Uh(X) from XL to XR, using the fact that Uh(X) = C(1) * P1(X) + C(2) * P2(X) + C(3) * P3(X) + C(4) * P4(X) + C(5) * P5(X).
Create a second plot of the exact solution, U(X) = 8*X^2-10X+3.
How close do the plots look?
Here is a standard example problem for the finite element method:
Over the interval:
0 <= x <= 1,
solve the partial differential equation:
- u'' + u = x
with boundary conditions:
u(0.0) = 0, u(1.0) = 0.
The exact solution is
u(x) = x - sinh(x) / sinh(1.0)
Notice that the left hand side of the differential equation is more complicated than in our last example, because we have the additional term "u". How does this affect our calculation? Really, all it does is add a term to the stiffness matrix. For example, the matrix entry A(2,3) will change from
Int P2' P3'
to
Int P2' P3' + P2 P3
and the right hand side F won't change at all.
So we can set up and solve this problem in a way that is similar to the previous problem. But now it is time to ask about accuracy. How good is our solution? Looking at a plot gives us a "feeling" that our answer is good or bad, but to be objective, we need a number that evaluates our approximation.
Since we have called the exact solution u(x), we will need to call our approximate solution from the finite element method something else, and the traditional name is uh(x), where the "h" is meant to suggest that this solution comes from elements of size h. The two items u(x) and uh(x) are both functions, defined at every point in [0,1]. Error refers to the difference between these two functions, so we can even talk about the error function, e(x) = u(x) - uh(x). But we'd much rather measure a single number that represents the error, and we'd like it to have reasonable properties, such as being zero if we had a perfect approximation.
A numeric measurement of error is known as an error norm, represented by ||e(x)||. Here are three error norms commonly used in finite element calculations:
Note that the H1 seminorm is called a seminorm because it is not quite a norm. If the norm of e(x) is 0, this means e(x) itself must be zero. But if the H1 seminorm of e(x) is 0, this only means that e'(x) is zero, which means e(x) is a constant...and not necessarily the 0 constant! Nonetheless, the H1 seminorm is of great interest because it measures how well the finite element method is recovering the derivative information, something more than just the value of the solution itself.
So our question about accuracy can be answered, once we have computed the approximate solution uh(x), by computing the three norms of e(x).
sqrt ( Integral ( u(x) - uh(x) )^2 dx )
and we can do that by considering each element in turn,
setting up a quadrature rule in that element, and at each
quadrature point xq evaluating the exact function u(xq)
and the finite element function uh(xq), subtracting, squaring,
multiplying by the quadrature weight.
uh'(x) = c0 * phi'(0,x) + c1 * phi'(1,x) + ...
or
uh'(x) = C1 * P1'(x) + C2 * P2'(x) + ... if you use 1-based indexing
Using 5 nodes and 4 elements, apply the finite element method to the test problem, and compute the max norm, L2 norm, and H1 seminorm of the error.
Now we know how to make measurements of the approximation error between a known exact solution, and a given finite element solution that uses a given number of elements of typical size H. It is natural to assume that, the harder we try, the better our approximation will become. That is, as H goes to zero, we expect the error norm to go to zero as well. If this is true, we can say that for the corresponding sequence of finite element approximations converges to the exact solution of the problem.
In order to carry out the tasks described in this section, we need to have a version of the finite element program which behaves logically like this:
NE --> Program --> MXerr, L2err, H1err.
That is, the program should have as input NE, the number of elements.
It should set up and solve the model problem from the previous section,
and return the maximum, L2 and H1 error norms. The purpose of this
section is to observe what happens to the errors as we change NE.
To begin with, we expect the error to go down as H decreases. Does it? Run the program for NE = 1, 2, 4, 8, 16 elements, and make a table of the three error measures. (I will presume your error generally decreases for this exercise.)
The next question is, what is the best way to plot the behavior of errors? Do we plot against NE, the number of elements, or H = 1/NE, the size of the elements over the interval [0.0,1.0]? Let's concentrate on just the L2 error, and make plots of L2err as a function of H, and as a function of NE. (I presume that the plots will not be very interesting to look at, because the initial error will be large, and the other errors won't show up.)
Plots using H or NE as the "x" variable really display the same information, but somehow, the H plots are psychologically confusing. This is because we read from left to right, and we expect actions to move from left to right. In a plot of error versus H, as we proceed, we are decreasing H, and hence our graph is moving from right to left. Moreover, if we add more H data, it will be crowded into a smaller and smaller space near the origin. Of course, using NE as our "x" variable fixes the psychological problem, and our most recent data doesn't get crowded into a corner, but we still can't see what's going on since the latest error data tends to be too small to see.
Now we need to ask, how fast does the error go down?. If we can answer this question, then we can estimate what to do to H if we want a new finite element approximation whose error is half what we got with the current one. And this will turn out to help us to fix our plotting problem as well.
For a moment, forget that we are approximating the function u(x) by solving a finite element problem. Assume we are simply trying to approximate u(x) by a function ui(x) using interpolation with piecewise linear functions. Then we would guess that the approximation error at a point could be described by
| u(x) - ui(x) | <= h^2 * u''(x*) <= h^2 ||u''|| <= H^2 ||u''||
where h is the size of the subinterval containing x, x* is some point in that
subinterval, and ||u''|| can be taken as the infinity norm of u'' over the
entire interval, and H is the maximum length of any subinterval (so H = h
in the common uniform case.)
It turns out that the finite element approximation uh(x) computed using piecewise linear functions is almost as good as the interpolant, that is
| u(x) - uh(x) | <= C * H^2 * ||u''||
where C is some number bigger than 1. This means, in turn, that
the L2 error will behave similarly:
||e(x)|| = ||u(x)-uh(x)|| <= L * C * H^2 * ||u''||
where L is the total length of the interval. A similar bound holds
for the maximum error. Meanwhile, it turns out that the H1 seminorm
error will lose one power in H:
||e'(x)|| = ||u'(x)-uh'(x)|| <= L * C * H * ||u''||
What's important that, for piecewise linear finite elements:
Max error is O(H^2) = O(1/NE^2)
L2 error is O(H^2) = O(1/NE^2)
H1 error is O(H) = O(1/NE)
Here, O(H^2) is to be read as "of the order of H^2".
The error relationships tells us what we can expect as we decrease H. It says that, for a given problem, doubling the number of elements should decrease the max and L2 errors by a factor of 4. (Strictly speaking, it's the bound on the error that decreases, but it is common to speak as though the error itself behaves this way, and it often does.)
Now we need to think again about our error plots. If there is a simple power relationship expected between H or NE and the error norm, then a logarithmic plot is the best way to display this. The exponent of the relationship can be used to draw a line of the correct slope. We would then expect the error curve to have approximately this slope. In particular, plotting log(NE) versus log(L2 error) we would hope to see behavior that tends to a straight line with slope -2, because L2 error behaves like NE^(-2).
Create a log/log plot of NE versus L2 error for the sequence NE = 1, 2, 4, 8, 16, 32, 64. If you can force the axes to use the same scale, then you may see that the (approximately) straight line has slope -2.
Now if we know the rate at which the error decreases with H, then we can do convergence plots or tables. This is essentially computing the slope of each of the individual line segments that make up the plot of NE versus L2 error. To do this, compute the Max, L2 and H1 errors for NE = 1, 2, 4, 8, 16, 32 and 64. Then fill in the following table:
MAX L2 H1
error( 1)/error( 2) | | | |
error( 2)/error( 4) | | | |
error( 4)/error( 8) | | | |
error( 8)/error(16) | | | |
error(16)/error(32) | | | |
error(32)/error(64) | | | |
We expect to see numbers approach 2, 2, and 1 respectively in each column.
So far, we've been using piecewise linear basis functions. There are good reasons for this: the functions are simple to program, the derivatives are easy to compute. However, it is important to realize that there is a wide variety of basis functions available, especially in 2D and 3D. In order to get an idea of what changes and what stays the same, we will look at the use of piecewise quadratic basis functions for our 1D problem.
The following diagram suggests how an interval [XL,XR] might be divided into 4 elements that use piecewise quadratic basis functions. Note that each element now consists of two (shared) boundary nodes and one interior node.
X: XL---------------------------------------------XR
Elements: [----E1----][----E3----][----E3----][----E4----]
Nodes: N1----N2----N3----N4----N5----N6----N7----N8----N9
Now it will become clear that, instead of defining our problem in terms of the number of nodes NN, we should think in terms of the number of elements NE. One very simple reason for this is that we can define a mesh with any number of elements, but we cannot specify a mesh by picking a number of nodes that is even. Moreover, we will see that, in a mesh of piecewise quadratic elements, nodes on the boundary between two elements will have different properties than nodes in the interior.
In particular, boundary nodes will allow jumps in the first derivative, they may be spaced irregularly, and their basis functions extend over two elements and hence interact with basis functions at 5 nodes. On the other hand, at interior nodes, all derivatives will be continuous, these nodes should be evenly spaced (that is, halfway between the two boundary nodes), and their basis functions will only interact with basis functions at 3 nodes, namely, themselves and their two boundary nodes.
Another thing to note is that, since we are using piecewise quadratic basis functions, we may need to use a more accurate quadrature rule to approximate our integrals.
Piecewise quadratic basis functions: Let's use the symbol Psii(x) for our piecewise quadratic basis functions. We want to work out the formula for evaluating Psii(x) for a given point x.
Psi = ( x - x(i-1) ) / ( x(i) - x(i-1) )
* ( x - x(i+1) ) / ( x(i) - x(i+1) )
Psi = ( x - x(i-2) ) / ( x(i) - x(i-2) )
* ( x - x(i-1) ) / ( x(i) - x(i-1) )
Psi = ( x - x(i+1) ) / ( x(i) - x(i+1) )
* ( x - x(i+2) ) / ( x(i) - x(i+2) )
Piecewise Quadratic Basis Functions: We need a procedure, which we can call pql(), to evaluate finite element basis functions of piecewise quadratic type. The information that this procedure will need is:
X, the evaluation point.
I, the index of the basis function, which must be between 0 and NN-1.
NN, the number of nodes (which must be odd).
XN, the node coordinates.
The function should return V, the value of the basis function.
Write a function of the form
double pql ( double x, int i, int nn, double xn[] )
which evaluates the I-th piecewise quadratic basis function at X.
Piecewise Quadratic Basis Function Plot: Create a calling program that defines NN = 5 nodes at XN = [ 0.0, 1.0, 3.0, 3.5, 10.0 ]. Plot, over the interval [0,10], the basis function associated with the node at X = 1.0, and the basis function associated with the node at X = 3.0.
Bring a PRINTED copy of your pql() function. Bring a PRINTED copy of the plots of the basis functions.
Today we will try to quickly wrap up the remaining pieces needed in order to create a 1D FEM code that uses piecewise quadratric basis functions. The issues include:
Derivatives of Piecewise quadratic basis functions: As before, let Psii(x) stand for the piecewise quadratic basis function associated with node i. The formula for the derivative of Psi is easy:
Psi' = 1 / ( x(i) - x(i-1) )
* ( x - x(i+1) ) / ( x(i) - x(i+1) )
+ ( x - x(i-1) ) / ( x(i) - x(i-1) )
* 1 / ( x(i) - x(i+1) )
Psi' = 1 / ( x(i) - x(i-2) )
* ( x - x(i-1) ) / ( x(i) - x(i-1) )
+ ( x - x(i-2) ) / ( x(i) - x(i-2) )
* 1 / ( x(i) - x(i-1) )
Psi' = 1 / ( x(i) - x(i+1) )
* ( x - x(i+2) ) / ( x(i) - x(i+2) )
+ ( x - x(i+1) ) / ( x(i) - x(i+1) )
* 1 / ( x(i) - x(i+2) )
Improved Quadrature: if we are using piecewise quadratic basis functions, then the integrals we need to approximate will be of higher polynomial degree than when we were using piecewise linear basis functions. A reasonable precaution is to increase the precision of the quadrature rule accordingly. The following 3-point rule will precisely integrate any polynomial up to degree 5, so that means that our stiffness and mass matrices would be exact, at least:
nq = 3
w = [ 0.5555555, 0.8888888, 0.5555555 ]
xi = [ -0.7745966, 0.0000000, 0.7745966 ].
Matrix assembly: as in the piecewise linear case, it is recommended that the matrix assembly be done by completely processing each element, creating an "element matrix", and then adding each element's contribution to the global matrix.
In a simple case of 3 elements, the matrix would be 7x7, and after processing the first element, the contributions from element 1 can be shown as:
1 2 3 4 5 6 7
1: 1 1 1 . . . .
2: 1 1 1 . . . .
3: 1 1 1 . . . .
4: . . . . . . .
5: . . . . . . .
6: . . . . . . .
7: . . . . . . .
Moving to element 2, we can "see" three nodes: 3, 4, and 5. Assuming we are simply working on a Poisson problem, the local element matrix is therefore:
3 4 5
3: Int psi'(3) psi'(3) Int psi'(3) psi'(4) Int psi'(3) psi'(5)
4: Int psi'(4) psi'(3) Int psi'(4) psi'(4) Int psi'(4) psi'(5)
5: Int psi'(5) psi'(3) Int psi'(5) psi'(4) Int psi'(5) psi'(5)
When we add the contributions from element 2 to the global matrix, we have:
1 2 3 4 5 6 7
1: 1 1 1 . . . .
2: 1 1 1 . . . .
3: 1 1 1+2 2 2 . .
4: . . 2 2 2 . .
5: . . 2 2 2 . .
6: . . . . . . .
7: . . . . . . .
and the symbol "1+2" means that matrix entry (3,3) contains contributions
from elements 1 and 2; but that's because the stiffness matrix entry for
basis function 3 needs to extend over both elements.
Completing the computation with element 3, we will have:
1 2 3 4 5 6 7
1: 1 1 1 . . . .
2: 1 1 1 . . . .
3: 1 1 1+2 2 2 . .
4: . . 2 2 2 . .
5: . . 2 2 2+3 3 3
6: . . . . 3 3 3
7: . . . . 3 3 3
and thus you can see that, when piecewise quadratic basis functions
are used, the matrix structure changes from what we saw for the
piecewise linear case. Now it is clear that the contributions from
each element show up as a block, and that consecutive blocks overlap
at a single matrix entry.
Try to do what's necessary to get a new version of your program set up, now using piecewise quadratic basis functions.
Apply your program to the same standard problem as before:
Over the interval:
0 <= x <= 1,
solve the partial differential equation:
- u'' + u = x
with boundary conditions:
u(0.0) = 0, u(1.0) = 0.
The exact solution is
u(x) = x - sinh(x) / sinh(1.0)
As a simple comparison, run your piecewise linear code with 8 elements (=9nodes) and your piecewise quadratic code using 4 elements (=9 nodes). Print a table that lists the exact solution and the two approximate solutions at the 9 nodes.
You can return to the FEP_2015 web page.