$
\newcommand{\dq}[2]{\displaystyle \frac{\partial #1}{\partial #2}}
\newcommand{\oq}[2]{\partial #1 / \partial #2}
\def \J {\Delta}
$
Let's consider the simplest non-trivial finite element shape in 3-D, which is
a linear tetrahedron. Function behaviour is approximated inside such a
tetrahedron by a linear interpolation between the function values at the
vertices, also called nodal points. Let $T$ be such a function, and $x,y,z$
coordinates, then:
$$
T = A.x + B.y + C.z + D
$$
Where the constants A, B, C, D are yet to be determined.
Substitute $ x=x_k $ , $ y=y_k $ , $ z=z_k $ with $ k=0,1,2,3 $.
Start with:
$$
T_0 = A.x_0 + B.y_0 + C.z_0 + D
$$
Clearly, the first of these equations can already be used to eliminate the
constant $D$, once and forever:
$$
T - T_0 = A.(x - x_0) + B.(y - y_0) + C.(z - z_0)
$$
Then the constants $A$ , $B$ , $C$ are determined by:
$$ \begin{array}{ll}
T_1 - T_0 = A.(x_1 - x_0) + B.(y_1 - y_0) + C.(z_1 - z_0) \\
T_2 - T_0 = A.(x_2 - x_0) + B.(y_2 - y_0) + C.(z_2 - z_0) \\
T_3 - T_0 = A.(x_3 - x_0) + B.(y_3 - y_0) + C.(z_3 - z_0)
\end{array} $$
Three equations with three unknowns. A solution can be found:
$$
\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]
$$
It is concluded that $A,B,C$ and hence $(T-T_0)$ must be a linear expression
in the $(T_k-T_0)$:
$$
T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)
$$ $$ =
\left[ \begin{array}{ccc} \xi & \eta & \zeta \end{array} \right]
\left[ \begin{array}{c} T_1-T_0 \\ T_2-T_0 \\ T_3-T_0 \end{array} \right]
$$
See above:
$$ =
\left[ \begin{array}{ccc} \xi & \eta & \zeta \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]
\left[ \begin{array}{c} A \\ B \\ C \end{array} \right]
$$
See above:
$$ = T - T_0 =
\left[ \begin{array}{ccc} x-x_0 & y-y_0 & z-z_0 \end{array} \right]
\left[ \begin{array}{c} A \\ B \\ C \end{array} \right]
$$
Hence:
$$ \begin{array}{ll}
x - x_0 = \xi .(x_1 - x_0) + \eta.(x_2 - x_0) + \zeta.(x_3 - x_0) \\
y - y_0 = \xi .(y_1 - y_0) + \eta.(y_2 - y_0) + \zeta.(y_3 - y_0) \\
z - z_0 = \xi .(z_1 - z_0) + \eta.(z_2 - z_0) + \zeta.(z_3 - z_0)
\end{array} $$
But also:
$$
T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)
$$
Therefore the same expression holds for the function $T$ as well as for
the coordinates $x,y,z$. This is called an isoparametric transformation.
It is remarked without proof that the local coordinates $\xi,\eta,\zeta$
within a tetrahedron can be interpreted as sub-volumes, spanned by the vectors
$ \vec{r}_k-\vec{r}_0 $ and $\vec{r}-\vec{r}_0 $ where $\vec{r}=(x,y,z) $ and
$k=1,2,3$.
Reconsider the expression:
$$
T - T_0 = \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0)
$$
Partial differentiation to $ \xi $ , $ \eta $ , $ \zeta $ gives:
$$
\oq{T}{\xi} = T_1 - T_0 \quad ; \quad \oq{T}{\eta} = T_2 - T_0
\quad ; \quad \oq{T}{\zeta} = T_3 - T_0
$$
Therefore:
$$
T = T(0) + \xi \dq{T}{\xi} + \eta \dq{T}{\eta} + \zeta \dq{T}{\zeta}
$$
This is part of a Taylor series expansion around node (0). Such Taylor series
expansions are very common in Finite Difference analysis.
Now rewrite as follows:
$$
T = (1 - \xi - \eta - \zeta).T_0 + \xi.T_1 + \eta.T_2 + \zeta.T_3
$$
Here the functions $ (1-\xi-\eta-\zeta),\xi,\eta,\zeta $ are called the
shape functions of a Finite Element. Shape functions $N_k$ have the property
that they are unity in one of the nodes (k), and zero in all other nodes.
In our case:
$$
N_0 = 1-\xi-\eta-\zeta \quad ; \quad N_1 = \xi \quad ; \quad N_2 = \eta
\quad ; \quad N_3 = \zeta
$$
So we have two representations, which are allmost trivially equivalent:
$$ \begin{array}{ll}
T = T_0 + \xi.(T_1 - T_0) + \eta.(T_2 - T_0) + \zeta.(T_3 - T_0) \quad
& \mbox{: Finite Difference} \\
T = (1 - \xi - \eta - \zeta).T_0 + \xi.T_1 + \eta.T_2 + \zeta.T_3 \quad
& \mbox{: Finite Element}
\end{array} $$
%
What kind of terms can be discretized at the domain of a linear tetrahedron?
In the first place, the function $T(x,y,z)$ itself, of course. But one may also
try on the first order partial derivatives $\oq{T}{(x,y,z)}$. We find:
$$
\oq{T}{x} = A \quad ; \quad \oq{T}{y} = B \quad ; \quad \oq{T}{z} = C
$$
Using the expressions which were found for $A,B,C$:
$$
\left[ \begin{array}{c} \oq{T}{x} \\ \oq{T}{y} \\ \oq{T}{z} \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]
$$
It is seen from this formula that one must determine the inverse of the above
matrix first. This can be done with Cramer's rule, as follows:
subroutine dinv3(mat,det)
* =====
* Direct inversion of 3 x 3 matrix
* --------------------------------
real mat(3,3),sub(3,3)
*
sub(1,1)=+mat(2,2)*mat(3,3)-mat(2,3)*mat(3,2)
sub(2,1)=-mat(1,2)*mat(3,3)+mat(1,3)*mat(3,2)
sub(3,1)=+mat(1,2)*mat(2,3)-mat(1,3)*mat(2,2)
*
sub(1,2)=-mat(2,1)*mat(3,3)+mat(2,3)*mat(3,1)
sub(2,2)=+mat(1,1)*mat(3,3)-mat(1,3)*mat(3,1)
sub(3,2)=-mat(1,1)*mat(2,3)+mat(1,3)*mat(2,1)
*
sub(1,3)=+mat(2,1)*mat(3,2)-mat(2,2)*mat(3,1)
sub(2,3)=-mat(1,1)*mat(3,2)+mat(1,2)*mat(3,1)
sub(3,3)=+mat(1,1)*mat(2,2)-mat(1,2)*mat(2,1)
*
det=mat(1,1)*mat(2,2)*mat(3,3)-mat(1,1)*mat(2,3)*mat(3,2)
+ -mat(2,1)*mat(1,2)*mat(3,3)+mat(2,1)*mat(1,3)*mat(3,2)
+ +mat(3,1)*mat(1,2)*mat(2,3)-mat(3,1)*mat(1,3)*mat(2,2)
*
if(det.eq.0.) stop 'dinv3: det=0'
*
do 10 i=1,3
do 10 j=1,3
10 mat(i,j)=sub(j,i)/det
*
return
end
While carrying out this algorithm, terms can be rewritten as in:
$$
(x_1 - x_0)(y_2 - y_0) - (x_2 - x_0)(y_1 - y_0) =
$$ $$
= x_1 y_2 + x_2 y_0 + x_0 y_1 - x_2 y_1 - x_0 y_2 - x_1 y_0 =
\left| \begin{array}{ccc} 1 & x_0 & y_0 \\
1 & x_1 & y_1 \\
1 & x_2 & y_2 \end{array} \right|
$$
The following shorthand notation will be used for all this:
$$
(x,y;0,1,2) =
\left| \begin{array}{ccc} 1 & x_0 & y_0 \\
1 & x_1 & y_1 \\
1 & x_2 & y_2 \end{array} \right|
$$
Let's do the real thing now. Inverting the above matrix then results in:
$$
\left[ \begin{array}{c} \oq{T}{x} \\ \oq{T}{y} \\ \oq{T}{z} \end{array}\right]
= \left[ \begin{array}{ccc} (y,z;2,3,0) & (y,z;3,0,1) & (y,z;0,1,2) \\
(z,x;2,3,0) & (z,x;3,0,1) & (z,x;0,1,2) \\
(x,y;2,3,0) & (x,y;3,0,1) & (x,y;0,1,2)
\end{array}\right] / \J
\left[ \begin{array}{c} T_1-T_0 \\ T_2-T_0 \\ T_3-T_0 \end{array}\right]
$$
Here $\J$ is the determinant of the original matrix as a whole.
In order to find the coefficients belonging to $T_0$, add up the columns of
the inverted matrix and provide the sum with a minus sign:
$$
- \left[ \begin{array}{c} (y,z;2,3,0) + (y,z;3,0,1) + (y,z;0,1,2) \\
(z,x;2,3,0) + (z,x;3,0,1) + (z,x;0,1,2) \\
(x,y;2,3,0) + (x,y;3,0,1) + (x,y;0,1,2)
\end{array}\right] T_0 =
\left[ \begin{array}{c} (y,z;1,2,3) \\ (z,x;1,2,3) \\ (x,y;1,2,3)
\end{array}\right] T_0
$$
The righthand side being the result of a little exercise in elementary algebra.
We have found a $3 \times 4$ Differentiation Matrix, which represents
the gradient operator $\oq{}{(x,y,z)}$ for the function values $T_{0,1,2,3}$
at a linear tetrahedron:
$$
\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
$$
It is seen that each column of the differentiation matrix corresponds with a
component of the normal (cross product) belonging to the area of the triangle
opposite to the vertex where the accompanying temperature is defined, thereby
everything being divided by the volume of the tetrahedron as a whole.
This for example means that the gradient of the temperature field is described
in a sensible way as a flux through all triangle boundaries of the tetrahedron.
The above procedure of finding the differentiation matrix entries for $T_0$
thus means that the normal on triangle $(1,2,3)$ is equal to minus the sum of
the normals on the other triangles. Which in turn means that the (vector) sum
of all normals is equal to zero, which indeed should be the case for any
closed (tetrahedal) surface.