Remaining parts
$
\def \J {\Delta}
$
The last term to be discretized is the multiplicative factor $a_0$. Due to the
fact that everything is integrated at the element vertices, this represents no
problem whatsoever. Simply add terms $a_0$ times $\J_k/n$ to the diagonal of
the element matrix under construction.
The entire local discretization is contained in one Fortran subroutine,
called ELEMENT. Two building blocks ("bricks") are implemented herein: a prism
with 6 vertices and a hexahedron with 8 vertices. Both are ultimately built up
from (mutually overlapping) tetrahedra, which are located with their "origins"
$0$ at the vertices of the bricks. A step in determining the Differentiation
Matrices is the inversion of a $3 \times 3$ matrix by subroutine DINV3.
Needed for the element matrices are two things: local geometry information and
the coefficients of the CR equation, which also depend upon the coordinates.
Calculating the geometry for a hemisphere is elementary. It is done by a piece
of in-line code in the MAIN program. Calculating the coefficients of the CR
equation is done in the ELEMENT subroutine itself: subroutine COEFFS.
The local element matrices are assembled into a large a-symmetrical (banded)
global matrix by the traditional finite element Assembly procedure. What's
needed for such a procedure is the so called connectivity or topology of the
mesh. We have decided also to keep the connectivity information in-line.
Use has been made of a Fortran statement function $nr(i,j,k)$, which basically
is an association between the finite element $(nr)$ and the finite difference
$(i,j,k)$ numbering schemes, thus combining (again) the best of both worlds.
There are two types of Boundary Conditions associated with CR. The "natural"
type B.C. which is applicable for the plane $z=0$ was discussed already in the
section "Diffusion part". Dirichlet type boundary conditions occur at the inner
and the outer hemi-spheres, which furthermore delimit the integration domain.
Pressure values must be prescribed there. In the demonstration model, values 0
and 1 are adopted, respectively.