In exercise 4 the
global matrix
M_G
will be assembled for a ten-element projection problem. The most
involved new concept this exercise will be the (generally) non-square
integer mapping matrix which I will call
map. To allocate
memory for the m by n matrix map I made a function called imatrix which
performs this task when the command
map = imatrix(m, n);
is entered. I started my coding from the one-element code of exercise
3c. First the mapping
chi has to be inserted. Then I made
a function called
mapping which constructs the
map
matrix and a function
elbound which fills the vector
bound
with the x-positions of the element boundaries (in this case
bound
will just be
{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} because
the
Nel=10 elements are equispaced). The variables
M_G
and
f_G are the global mass matrix and right-hand-side
vector and now
Mdim is the size of the global mass matrix
(thus for now
Mdim=Nel*P+1). Now the Lapack routines must
be given the global variables as arguments and also the function which
calculates the (global) solution
u_delta,
sol,
must be changed appropriately. Remember that every two neighbouring
elements have one gridpoint in common when you construct the global
solution. I also made a global x-vector called
x_G, which
is filled in the function
xGlob. The exercise (and
my code) use a sine as forcing
function but
I think you can just use any function you like. From now on all my
output (that is the output of my codes) will be suited for input in
Matlab.
To finalize the projection problem in part (b) we will implement
Dirichlet bc's. Doing this involves changing the first and last element
mass matrices, changing the mapping matrix, constructing the known
solution
u_D which satisfies the bc's and changing the
construction of the homogeneous solution
u_H in the
sol
function. To construct the known (local) solution
u_D I
defined a function called
rhbc. This function only has to
be called during construction of the first and last element. The
assembly of these elements will happen in the function
assembound
which thus will call
rhbc.
assembound will
be called when the global assembly is at the point of inserting element
eln=0 or element
eln=Nel-1 (which of course
are the first and the last element). The
sol function has
to be changed considerably and I also defined a global
u_D
called
u_DG. Don't forget to change the size of the
global matrix to
Mdim=Nel*P-1, the 2-element vector
bc
will again contain the (right!) bc's. Here is
my
solution.