FEM discretization
$
\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}{c} \oq{T}{x} \\ \oq{T}{y} \\ \oq{T}{z} \end{array} \right]
$$
Here $T = $ temperature, $x,y,z = $ coordinates.
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.
The second term is an integral for the bulk material.
Substitute temperature 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}{c} \oq{T}{x} \\ \oq{T}{y} \\ \oq{T}{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.
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
contributions belonging to the diffusion term, using differentiation matrices
$\oq{}{(x,y,z)}$:
$$
\left[ \begin{array}{ccc} \oq{}{x} & \oq{}{y} & \oq{}{z} \end{array} \right]
\left[ \begin{array}{c} \oq{}{x} \\ \oq{}{y} \\ \oq{}{z} \end{array} \right]
\J
$$
Here $\J=$ determinant (volume) of the tetrahedron.
Using the end-result of the preceding chapter:
$$
\left[ \begin{array}{c} \oq{}{x} \\ \oq{}{y} \\ \oq{}{z} \end{array}\right] =
\left[ \begin{array}{cccc}
(y,z;1,2,3) & (y,z;2,3,0) & (y,z;3,0,1) & (y,z;0,1,2) \\
(z,x;1,2,3) & (z,x;2,3,0) & (z,x;3,0,1) & (z,x;0,1,2) \\
(x,y;1,2,3) & (x,y;2,3,0) & (x,y;3,0,1) & (x,y;0,1,2)
\end{array}\right] / \J
$$
Now we are going to use the accompanying interpretation: the columns of the
differentiation matrix are components of normals $\vec{N}$ on triangles which
are at the boundary of each tetrahedron. Finite element matrices for diffusion
can be written down then in the following form, finally:
$$
\left[ \begin{array}{cccc}
\vec{N}_{123} \cdot \vec{N}_{123} &
\vec{N}_{123} \cdot \vec{N}_{230} &
\vec{N}_{123} \cdot \vec{N}_{301} &
\vec{N}_{123} \cdot \vec{N}_{012} \\
\vec{N}_{230} \cdot \vec{N}_{123} &
\vec{N}_{230} \cdot \vec{N}_{230} &
\vec{N}_{230} \cdot \vec{N}_{301} &
\vec{N}_{230} \cdot \vec{N}_{012} \\
\vec{N}_{301} \cdot \vec{N}_{123} &
\vec{N}_{301} \cdot \vec{N}_{230} &
\vec{N}_{301} \cdot \vec{N}_{301} &
\vec{N}_{301} \cdot \vec{N}_{012} \\
\vec{N}_{012} \cdot \vec{N}_{123} &
\vec{N}_{012} \cdot \vec{N}_{230} &
\vec{N}_{012} \cdot \vec{N}_{301} &
\vec{N}_{012} \cdot \vec{N}_{012}
\end{array}\right] / \J
$$
Where it should be mentioned that, in addition:
$$
\vec{N}_{123} + \vec{N}_{230} + \vec{N}_{301} + \vec{N}_{012} = 0
$$
Hence it's easy to show that each row of the matrix sums up to zero.
If the above result is specialized for a Finite Difference mesh, then the three
normals $\vec{N}_{012}$, $\vec{N}_{301}$ and $\vec{N}_{230}$ will be orthogonal
to each other. It's a nice exercise to show that the well known F.D. scheme for
diffusion can be reconstructed in this case.