sci.math.numanalysis
SUNA46, Linear tetrahedron
==========================
Let's consider the simplest nontrivial finite element shape in 3 dimensions,
the linear tetrahedron. Function behaviour is approximated inside a tetrahedron
by a linear interpolation between the function values at the vertices, or nodal
points. Let T be such a function, and x,y,z the coordinates, then:
T = A.x + B.y + C.z + D
Hence: T1 = A.x1 + B.y1 + C.z1 + D
Substract: T  T1 = A.(x  x1) + B.(y  y1) + C.(z  z1)
Then the constants A , B and C are determined by solving the equations:
 T2  T1   x2  x1 y2  y1 z2  z1   A 
 T3  T1  =  x3  x1 y3  y1 z3  z1   B 
 T4  T1   x4  x1 y4  y1 z4  z1   C 

matrix J
The determinant of the matrix called J represents 6 times the volume of the
tetrahedron (1,2,3,4). Write the substraction as:
T  T1 = h.(T2  T1) + k.(T3  T1) + l.(T4  T1)
Or: T  T1 =  T2  T1 T3  T1 T4  T1   h 
 k 
 l 
Then it may be a theorem that h , k , l too are a solution of the following 3
equations, an isoparametric transformation, namely:
 x  x1   x2  x1 x3  x1 x4  x1   h 
 y  y1  =  y2  y1 y3  y1 y4  y1   k 
 z  z1   z2  z1 z3  x1 z4  z1   l 

matrix J' where ' = transposed: see J above.
It's consistent: T  T1 =  T2  T1 T3  T1 T4  T1   h 
 k  =
 l 
=  A B C  . J' * (J')^1 .  x  x1  =  A B C   x  x1 
 y  y1   y  y1 
 z  z1   z  z1  Q.E.D.
Thus again the _same_ expression holds for the function T as well as for the
coordinates x , y , z , which an isoparametric transformation is meant to be.
h , k , l can be interpreted as volume coordinates (was: area coordinates).
Instead of volumecoordinates, we prefer to talk about the _local coordinates_
h , k , l of an element, in contrast to the _global coordinates_ x , y . z .
It is possible that local coordinates coincide with the global coordinates.
A tetrahedron for which this is the case is called a _parent element_.
Let's reconsider for a moment the expression:
T  T1 = h.(T2  T1) + k.(T3  T1) + l.(T4  T1)
Partial differentiation to h , k , l gives:
dT/dh = T2  T1 dT/dk = T3  T1 dT/l = T4  T1
So: T = T1 + h.dT/dh + k.dT/dk + l.dT/dl
This is part of a Taylor series expansion around node (1). Such Taylor series
expansions are very common in Finite Difference analysis.
Now rewrite as follows:
T = (1  h  k  l).T1 + h.T2 + k.T3 + l.T4
Here the functions (1  h  k) , h , k , l are called the SHAPE FUNCTIONS of
the Finite Element. Shape functions Nk have the property that they are unity
in one of the nodes (k), and zero in all other nodes. In our case:
N1 = 1  h  k ; N2 = h ; N3 = k ; N4 = l
So we have two representations, which are allmost trivially equivalent:
T = T1 + h.(T2  T1) + k.(T3  T1) + l.(T4  T1) : Finite Difference
T = (1  h  k).T1 + h.T2 + k.T3 + l.T4 : Finite Element
Partial derivatives dT/d(x,y,z) will be considered. In order to take care of
their discretization, transformation to isoparametric local coordinates (h,k,l)
at a Finite Element, being a general F.E. technique, seems to be advantageous.
Derivatives to the global coordinates must then be expressed in derivatives to
the local coordinates. According to the chain rules, written in matrix form:
 dT/dh   dx/dh dy/dh dz/dh   dT/dx 
 dT/dk  =  dx/dk dy/dk dz/dk   dT/dy 
 dT/dl   dx/dl dy/dl dz/dl   dT/dz 
According to Cramer's rule, the solution of the above equations is:
 dT/dh dy/dh dz/dh 
dT/dx = Det(  dT/dk dy/dk dz/dk  ) / J
 dT/dl dy/dl dz/dl 
 dx/dh dT/dh dz/dh 
dT/dy = Det(  dx/dk dT/dk dz/dk  ) / J
 dx/dl dT/dl dz/dl 
 dx/dh dy/dh dT/dh 
dT/dz = Det(  dx/dk dy/dk dT/dk  ) / J
 dx/dl dy/dl dT/dl 
 dx/dh dy/dh dz/dh 
Where: J = Det(  dx/dk dy/dk dz/dk  ) is the Jacobian determinant.
 dx/dl dy/dl dz/dl 
Thus solving for dT/d(x,y,z) we find general formulas which are valid for any
3D element. ^^^^^^^
It is remarked that the first three determinants can be derived from the last
one by substituting x := T , y := T , z := T respectively. More easy: due to
isoparametrics, only the nodal values of these quantities have to be replaced.
This means that the function behaviour of partial derivatives, being a quotient
of two alike determinants, can be studied by considering the Jacobian matrix or
determinant alone. This will be heavily used in the sequel. For the tetrahedron
(1,2,3,4) we trivially find:
T = T1 + h.(T2  T1) + k.(T3  T1) + l.(T4  T1)
dT/dh dT/dk dT/dl
Substitutution of T := x , y , z gives the determinant of the matrix:
 x2  x1 y2  y1 z2  z1 
 x3  x1 y3  y1 z3  z1 
 x4  x1 y4  y1 z4  z1 

* Han de Bruijn; Applications&Graphics  "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354  would be NO idleness in * Oil
* 2600 AJ Delft; The Netherlands.  Mathematics" (HdB). * for
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood