Diffusion part
$
\newcommand{\dq}[2]{\displaystyle \frac{\partial #1}{\partial #2}}
\newcommand{\oq}[2]{\partial #1 / \partial #2}
\def \J {\Delta}
$
Consider the three-dimensional Laplace-like term, which is defined in Cartesian
coordinates by:
$$
\dq{Q_x}{x} + \dq{Q_y}{y} + \dq{Q_z}{z}
$$
where:
$$
\left[ \begin{array}{c} Q_x \\ Q_y \\ Q_z \end{array} \right] =
\left[ \begin{array}{ccc}
\kappa_{xx} & \kappa_{xy} & \kappa_{xz} \\
\kappa_{yx} & \kappa_{yy} & \kappa_{yz} \\
\kappa_{zx} & \kappa_{zy} & \kappa_{zz} \end{array} \right]
\left[ \begin{array}{c} \oq{p}{x} \\ \oq{p}{y} \\ \oq{p}{z} \end{array} \right]
$$
Which is essentially the first term of the CR equation.
Suppose the contribution is valid in a domain $D$ with boundary $S$.
According to the so called Galerkin method, the contribution is multiplied by
an arbitrary function $f$ and then integrated over the Domain of interest.
Since the function $f$ is completely arbitrary (well, continuous at least),
this is supposed to be equivalent to the original problem:
$$
\iiint f(x,y,z)
\left( \dq{Q_x}{x} + \dq{Q_y}{y} + \dq{Q_z}{z} \right) \, dx.dy.dz
$$
The advantage of the Galerkin formulation is that we are able now to convert
second order derivatives into first order derivatives. To see how this works,
let us recall Green's theorem, or partial integration for triple integrals, by
which the following expression can be substituted for the Galerkin integral:
$$
\oint \! \! \oint f.Q_n \, dS \quad - \quad
\iiint \left( Q_x \dq{f}{x} + Q_y \dq{f}{y} + Q_z \dq{f}{z} \right)
\, dx.dy.dz
$$
The first of these terms incorporates boundary conditions at the surface $S$.
The boundary integral is zero in case the normal derivative $ Q_n = 0 $.
For this reason $Q_n = 0$ is called a natural boundary condition: it is
fulfilled automatically if the first term in the above formulation is
simply discarded. With the CR problem, such a natural boundary condition
occurs at the $ z=0 $ plane (at least in the case of interest for us). Due to
symmetry, the gradient of the pressure normal to it must vanish:
$$ \left( \dq{p}{z}=0 \right)_{z=0} $$
It can be shown that the $\stackrel{\leftrightarrow}{\kappa}$ tensor reads as
follows, when expressed in spherical coordinates:
$$ \left[ \begin{array}{ccc}
\kappa_{rr} & 0 & \kappa_{r\phi} \\
0 & \kappa_{\theta\theta} & 0 \\
\kappa_{\phi r} & 0 & \kappa_{\phi\phi} \end{array} \right]
$$
At $z=0$ the vector $\vec{e}_\theta$ coincides with $-\vec{e}_z$ and the normal
$\vec{n}$. Therefore:
$$ \left( Q_n = 0 \right)_{z=0} \quad \Longleftrightarrow \quad
\stackrel{\leftrightarrow}{\kappa} \nabla p \cdot \vec{n} =
\kappa_{\theta\theta} \dq{p}{\theta} = 0
\quad \Longleftrightarrow \quad \left( \dq{p}{z}=0 \right)_{z=0}
$$
Hence the natural boundary condition at $z=0$ is fulfilled if and only if the
symmetry condition is fulfilled. Implementation of the boundary condition would
have been a zillion times more difficult if this were not the case!
The second term is an integral for the bulk material. Substitute pressure
fluxes herein and watch out for the minus sign:
$$ - \iiint
\left[ \begin{array}{ccc} \dq{f}{x} & \dq{f}{y} & \dq{f}{z} \end{array} \right]
\left[ \begin{array}{ccc}
\kappa_{xx} & \kappa_{xy} & \kappa_{xz} \\
\kappa_{yx} & \kappa_{yy} & \kappa_{yz} \\
\kappa_{zx} & \kappa_{zy} & \kappa_{zz} \end{array} \right]
\left[ \begin{array}{c} \oq{p}{x} \\ \oq{p}{y} \\ \oq{p}{z} \end{array} \right]
dx.dy.dz
$$
Note that the second order derivatives have been removed indeed. It is possible
to handle second order problems with first order (linear) finite elements only.
If the integration domain is subdivided into finite elements E, then the above
integral is splitted up as a sum of integrals over these elements. Integration
is always carried out numerically, by using integration points
[ OC ].
A little bit of innovation is involved in recognizing that it doesn't make much
difference if the integration points are deliberately chosen in such a way that
evaluation is as easy as possible: we always select them at the vertices of any
elements involved. It can be shown that elements which are integrated in such a
way are in fact superpositions of linear tetrahedra. Linear tetrahedra
are so to speak the ultimate 3-D elements and there's no need for anything else
most of the time. At such linear tetrahedra, differentiations $\oq{}{(x,y,z)}$
are given as a matrix operation, with the Differentiation Matrices found
in an earlier stage. Here is the symbolic representation for the element-matrix
contribution belonging to the diffusion term, using differentiation matrices
$\oq{}{(x,y,z)}$:
$$ - \sum_{\mbox{tetrahedra}} \frac{\J}{n}
\left[ \begin{array}{ccc} \oq{}{x} & \oq{}{y} & \oq{}{z} \end{array} \right]
\left[ \begin{array}{ccc}
\kappa_{xx} & \kappa_{xy} & \kappa_{xz} \\
\kappa_{yx} & \kappa_{yy} & \kappa_{yz} \\
\kappa_{zx} & \kappa_{zy} & \kappa_{zz} \end{array} \right]
\left[ \begin{array}{c} \oq{}{x} \\ \oq{}{y} \\ \oq{}{z} \end{array} \right]
$$
Here $n=$ number of element-brick vertices, $\J=$ determinant of tetrahedron.
A formula which can be traced back easily in the Fortran listing of the ELEMENT
subroutine. In [ OC ]
it is suggested that a tensor should be represented as
a diagonal matrix, thus making it necessary to find eigenvalues and -vectors in
the first place. The above analysis shows that such a procedure is too much of
the good. (Well, I didn't believe it at first sight, but experience shows that
the coefficients obtained for the element matrix are always identical with both
methods.)