Convection part
$
\newcommand{\dq}[2]{\displaystyle \frac{\partial #1}{\partial #2}}
\newcommand{\oq}[2]{\partial #1 / \partial #2}
\def \J {\Delta}
$
There are several paths to understanding with respect to proper discretization
of the convective terms in the CR equation. The first thing to keep in mind is
that Finite Volume methods are to be preferred and that techniques of upwinding
according to this method are to be considered as common practice. There will be
no dispute about the pros and cons of this approach. Now it is well known that,
with Finite Volume methods, rectangular (equidistant) grids are the easiest to
work with. Having adopted the idea that tetrahedral (sub)elements are necessary
and sufficient, the only thing we have to do is to transform the coordinates of
the tetrahedral element into a rectangular coordinate system and then apply the
Finite Volume tradition within these new coordinates.
Let $(0,1,2,3)$ be the numbering of the vertices of an arbitrary tetrahedron.
Let $(u,v,w)$ be the velocity components in the global $(x,y,z)$ coordinates.
Let $(U,V,W)$ the components of the velocity, as seen in the local coordinate
system spanned by the sides $\overline{01}$, $\overline{02}$, $\overline{03}$
of the tetrahedron. Then the following relationship holds:
$$
\left[ \begin{array}{c} u \\ v \\ w \end{array} \right] =
U \left[ \begin{array}{c} x_1-x_0 \\ y_1-y_0 \\ z_1-z_0 \end{array} \right] +
V \left[ \begin{array}{c} x_2-x_0 \\ y_2-y_0 \\ z_2-z_0 \end{array} \right] +
W \left[ \begin{array}{c} x_3-x_0 \\ y_3-y_0 \\ z_3-z_0 \end{array} \right]
$$
The inverse of this transformation can be written as follows:
$$
\left[ \begin{array}{ccc} U & V & W \end{array} \right] =
\left[ \begin{array}{ccc} u & v & w \end{array} \right]
\left[ \begin{array}{ccc} x_1-x_0 & y_1-y_0 & z_1-z_0 \\
x_2-x_0 & y_2-y_0 & z_2-z_0 \\
x_3-x_0 & y_3-y_0 & z_3-z_0 \end{array} \right]^{-1}
$$
The transformed scheme for convection now simply reads:
$$
U (p_1 - p_0) + V (p_2 - p_0) + W (p_3 - p_0)
$$
In the chapter "Linear Tetrahedron" the following formula was derived:
$$
\left[ \begin{array}{c} A \\ B \\ C \end{array} \right] =
\left[ \begin{array}{ccc} x_1-x_0 & y_1-y_0 & z_1-z_0 \\
x_2-x_0 & y_2-y_0 & z_2-z_0 \\
x_3-x_0 & y_3-y_0 & z_3-z_0 \end{array} \right]^{-1}
\left[ \begin{array}{c} T_1-T_0 \\ T_2-T_0 \\ T_3-T_0 \end{array} \right]
$$
Together with:
$$
\oq{T}{x} = A \quad ; \quad \oq{T}{y} = B \quad ; \quad \oq{T}{z} = C
$$
Combining this with the above, it is inferred that:
$$
\left[ \begin{array}{ccc} U & V & W \end{array} \right]
\left[ \begin{array}{c} p_1-p_0 \\ p_2-p_0 \\ p_3-p_0 \end{array} \right] =
$$ $$
\left[ \begin{array}{ccc} u & v & w \end{array} \right]
\left[ \begin{array}{ccc} x_1-x_0 & y_1-y_0 & z_1-z_0 \\
x_2-x_0 & y_2-y_0 & z_2-z_0 \\
x_3-x_0 & y_3-y_0 & z_3-z_0 \end{array} \right]^{-1}
\left[ \begin{array}{c} p_1-p_0 \\ p_2-p_0 \\ p_3-p_0 \end{array} \right] =
$$ $$
\left[ \begin{array}{ccc} u & v & w \end{array} \right]
\left[ \begin{array}{c} \oq{p}{x} \\ \oq{p}{y} \\ \oq{p}{z} \end{array} \right]
= u \dq{p}{x} + v \dq{p}{y} + w \dq{p}{z}
$$
Just meaning that the theory is quite consistent with respect to the proposed
coordinate transformation.
In order to achieve proper weighting of the convection contribution as compared
with the contribution due to diffusion, it must be examined how the scheme fits
into the same Galerkin procedure as was applied to diffusion:
$$
\iiint f \left( u \dq{p}{x} + v \dq{p}{y} + w \dq{p}{z} \right) \, dx.dy.dz
$$
Again, the integration is carried out numerically, with the integration points
located at the vertices of the finite element brick. The result is:
$$
\frac{1}{n} \sum_{k=1}^n \left| \J_k \right| \left( \mbox{convex} \right)_k
$$
Where "convex" are the convective terms, discretized according to the above, $
n=$number of brick vertices, $\J_k=$determinant of tetrahedron at vertex $(k)$.
So far so good. What has been forgotten so far is just one more thing:
upwinding. Let's recall the transformed scheme for convection:
$$
U (p_1 - p_0) + V (p_2 - p_0) + W (p_3 - p_0) =
U p_1 + V p_2 + W p_3 - (U + V + W) p_0
$$
When applying the conventional Finite Volume upwind technique \cite{SV}, only
velocity components $(U,V,W)$ which have a negative sign are allowed to
contribute to the scheme. There is a small additional complication, which can
be explained as follows. It is assumed that the local grid is orthogonal in the
first place. Attention is paid to the velocity component $U$ only, which means
in effect that we have a one-dimensional problem. Brainless application of the
Galerkin method, without upwinding, would have resulted in a scheme like:
$$
\frac{U_2 - U_1}{2} + \frac{U_3 - U_2}{2} = \frac{U_3 - U_1}{2}
\qquad \mbox{instead of} \qquad (U_2 - U_1)
$$
The infamous & wrong zero-diagonal result is on the left. The correct upwind
scheme is on the right. It is seen that a factor $2$ in the denominator would
give insufficient weight to the latter. Hence we must apply $2 \times$ these
terms.