In part (a) of
exercise 3 we will
construct the element mass matrix
M. To do this we can
simply take the code of exercise 1a, add three functions and make a
minor change to the integration function. These are: the familliar
dmatrix,
a function called
basis to collect the basis functions
and a function called
assem which calls the different
fuctions to construct
M. the new variables (as compared
with 1a) are:
P for the order of the polynomial
expansion,
Mdim for the size of the matrix (which here is
equal to
P+1),
phi1 and
phi2
for the basis functions and of course
M for the mass
matrix. Also the function
func is not necessary. The
change to the integration function to make it capable of integrating
the product of two functions is rather easy to implement. To collect
the basis functions as given in equation (2.38) in section 2.3.3.3 of
the
book the routine
jacobfd is called from the function
basis.
jacobfd is documented
here.
Make sure that
the cases
p=0 and
p=P are taken care of
properly.
Here is my
solution. Because of
the many 'near-zero-numbers' it is probably difficult to see the
structure as shown in figure 2.6a in the book so it can be advantageous
to build a function which puts near-zero's to zero.
When using Lagrangian interpolants, as seen in equation (2.41) of the
book, the code looks like
this.
The only
change is in the basis function. To collect the right polynomials use
the polylib function
hgll which is the '
alpha=beta=0
version' of
hglj which is documented
here.
It can be seen
that when running the code with
np=P+1 we get a diagonal
M
(again neglecting the almost-zeros).
We now find ourselves at a crucial point, exercise 3c will be the first
exercise where a matrix will be inverted. This cost me quite a lot of
trouble but with my explanation it will not be difficult. First
however, the right hand side vector
f must be
constructed. It is not difficult to do this, it only takes an extra
command outside the inner loop in the
M assembly function
(which I called
assem). At this point for me the trouble
started. That's why I will now take some time to discuss several
features of the Lapack library. If you are already familliar with
Lapack it might be faster to
skip
this
introduction.
Lapack introduction
The
Lapack
library is a
library that provides a lot of matrix/vector operation routines. It is
assumed that Lapack is installed on your workstation and you have
access to it (which, as I am told, is quite common on a Linux
computer). Lapack is written in Fortran77 so you must remember to put
an underscore _ behind the function in the call. For the Lapack
functions there exists an entire
naming
system
which can be found in the Lapack
user's
guide.
To invert the matrix in exercise 3c I will, unlike the suggestion in
the exercise, use the LU-decomposition and solve routines
dgetrf_
and
dgetrs_. In these functions the first letter, thus
d,
means that the routines are of double precision. The next two letters,
ge,
mean that the routine can be used for general matrices (i.e.
unsymmetric). Finally the letters
trf mean factorisation
and
trs means use the factorisation and solve. To use a
routine from Lapack this routine must be declared in a special way, I
did this like this:
extern "C" {extern void dgetrf_(int *, int *,
double (*), int *, int [], int *);}
Thus without a semicolon in the end! Now the function
dgetrf
can be used just like a normal function (only with the extra
underscore). Because in Fortran functionarguments must be adresses we
must take care that in our code this is also the case. Arguments must
thus be pointers or normal variables preceded by the adress operator
&.
To find an explanation of the Lapack functions I just used
Google and typed (for
example)
"dgetrf
arguments", now the first (with dgetrs the third) reference is a page
with a good description. To use dgetrf we can enter in our code:
dgetrf_(&Mdim, &Mdim, M[0],
&Mdim, ipiv, &INFO);
The first, second and fourth argument are the adresses of the dimension
of the matrix. The third argument is the matrix itself.
ipiv
is an integer pointer the lenght of the matrix which will store the
pivot indices and
INFO is an integer which will get the
value
0 if the factorisation is completed succesfully. To
solve I used:
dgetrs_(&TRANS, &Mdim, &NRHS,
M[0], &Mdim, ipiv, f, &Mdim, &INFO);
TRANS is an unsigned character precision variable that is
set to the capital letter
T, thus:
TRANS = 'T';
which basically states that Lapack must use the transpose of A because
in Fortran matrices are stored in a transpose way as compared to C and
C++.
NRHS is an integer that gives the number of columns
of the right hand side vector (in our case 1) and
f is
the right hand side vector which after the solve is filled with our
solution vector
u. When we now compile the code we must
not forget to link lapack and g2c (g2c translates fortran to c), the
compiling commands now look like
this:
g++ -c polylib.c
g++ -c exercise3c.cpp
g++ -c exercise3c exercise3c.o polylib.o -llapack
-lg2c
If everything goes as planned we should now have the solution vector
u!
This now ends the
Lapack introduction
section and we will continue with exercise 3c again.
We now have the weigths of the different
polynomials in the vector
f, so to construct the solution
u_delta we only have to multiply them by their different
basis functions and sum. For this I made a function called
sol.
For the allocation of the integervector
ipiv I've made a
function
ivector. The output of
my
code is organised in
such a way that when copy-pasting it to Matlab and running it a plot
will be generated of the solution.
Imposing Dirichlet boundary conditions (bc's) is the aim of exercise
3d. To do this we decompose the solution into an unknown homogeneous
solution
u_H and a known solution which satisfies the
Dirichlet bc's
u_D. The weights for the solution
u_H
can be found by solving the matrix inversion where the basis functions
for the mass matrix
M run from
p=1 to
p=P-1
(thus the size of the matrix will become
Mdim=P-1).
Construction of
u_D will be done by a new function called
rhbc. Don't forget to change the right-hand-side vector to
the right form and the change in the summation in the construction of
u_H.
The two-element-vector
bc contains the bc's,
bc[0]
is the left bc and
bc[1] the right bc.
Here is my solution.
The only change that has to be made for part (e) is the inclusion of
the mapping function
chi as defined in
exercise
1b and the insertion
of the Jacobian
Jac at the right places. don't forget to
enter the
right bc's!
Oh,
my code.