sci.math.numanalysis
SUNA47, Isoparametric brick
===========================
Let's consider next the simplestbutone element shape in 3D: a brick.
It is assumed that the parent element is a UNIT CUBE. Take the origin 0 of the
local coordinate system in the middle of the cube: (h,k,l) = (1/2,1/2,1/2) .
The cube's vertices are numbered as summarized in the following BASIC program:
5 DIM COORD(8,3)
10 FOR I=1 TO 8 : PRINT I , : FOR N=1 TO 3
20 COORD(I,N)=.5*( 2*( ((I1)\2^(N1)) MOD 2 ) 1)
30 PRINT COORD(I,N), : NEXT N : PRINT : NEXT I
Output:
1 .5 .5 .5
2 .5 .5 .5
3 .5 .5 .5
4 .5 .5 .5
5 .5 .5 .5
6 .5 .5 .5
7 .5 .5 .5
8 .5 .5 .5
nr x y z
Function behaviour may be approximated inside such a cube by a Taylor expansion
around the origin. Since there are eight vertices, 8 terms must be retained:
T = T(0) + dT(0)/dh.h + dT(0)/dk.k + dT(0)/dl.l
2 2 2
+ d T(0)/dh.dk .h.k + d T(0)/dh.dl .h.l + d T(0)/dk.dl .k.l
3
+ d T(0)/dh.dk.dl .h.k.l
The interpolation between the function values at the nodes is thus trilinear.
Let's assume from the beginning that the same expression holds for the function
T as well as for the coordinates x and y : an isoparametric transformation.
Choose abbreviations for the names for the partial derivatives:
T = A1 + A2.h + A3.k + A4.l + A5.h.k + A6.h.l + A7.k.l + A8.h.k.l
Next specify the above expression for each of the vertices. This results in:
 T1   1 1/2 1/2 1/2 +1/4 +1/4 +1/4 1/8   A1 
 T2   1 +1/2 1/2 1/2 1/4 1/4 +1/4 +1/8   A2 
 T3   1 1/2 +1/2 1/2 1/4 +1/4 1/4 +1/8   A3 
 T4  =  1 +1/2 +1/2 1/2 +1/4 1/4 1/4 1/8   A4  F.E. < F.D.
 T5   1 1/2 1/2 +1/2 +1/4 1/4 1/4 +1/8   A5  transition
 T6   1 +1/2 1/2 +1/2 1/4 +1/4 1/4 1/8   A6 
 T7   1 1/2 +1/2 +1/2 1/4 1/4 +1/4 1/8   A7 
 T8   1 +1/2 +1/2 +1/2 +1/4 +1/4 +1/4 +1/8   A8 

1 h k l h h k h
k l l k
l
Before making any attempt to solve these equations, it should be remarked that
all of the matrix columns are _mutually orthogonal_. This is easy to check out.
At first it is seen that all columns sum up to zero: as many + as  signs.
The inner product of the first and the other columns are just the other columns
themselves. The inner products (h.k) , (h.l) , (k.l) are also other columns of
the matrix. There remains to investigate cases like h.h.k , a positive column
h.h times as many + as  signs in k . Same argument for things like h.h.k.l
and h.h.k.k.l . Q.E.D.
The inverse of an orthogonal matrix is the same as it's transpose. There is one
little problem, however: the lengths of the matrix columns are not equal to 1 .
But the latter is a rather simple matter to deal with: divide each row of the
transposed matrix by its (length)^2 . The result is:
 A1   +1/8 +1/8 +1/8 +1/8 +1/8 +1/8 +1/8 +1/8   T1 
 A2   1/4 +1/4 1/4 +1/4 1/4 +1/4 1/4 +1/4   T2 
 A3   1/4 1/4 +1/4 +1/4 1/4 1/4 +1/4 +1/4   T3 
 A4  =  1/4 1/4 1/4 1/4 +1/4 +1/4 +1/4 +1/4   T4  F.E. < F.D.
 A5   +1/2 1/2 1/2 +1/2 +1/2 1/2 1/2 +1/2   T5  transition
 A6   +1/2 1/2 +1/2 1/2 1/2 +1/2 1/2 +1/2   T6 
 A7   +1/2 +1/2 1/2 1/2 1/2 1/2 +1/2 +1/2   T7 
 A8   1 +1 +1 1 +1 1 1 +1   T8 

matrix C
So we find explicit discretizations for the partial derivatives of the function
T at the the origin 0 . The zero'th derivative is, for example:
T(0) = 1/8.(+ T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8)
Let's further abbreviate by introducing the functions L(i):
L1 L2 L3 L4 L5 L6 L7 L8
1 h k l h.k h.l k.l h.k.l
And adopt the Einstein convention, which is summation over equal indices. Then:
T = L(i).A(i) = L(i).C(i,j).T(j)
Shape functions are defined by: T = N(j).T(j)
Consequently: N(j) = L(i).C(i,j) .
We can work it out. For example:
N1 = + 1/8  1/4.(h + k + l) + 1/2.(h.k + h.l + k.l)  h.k.l
Guess what:
N1 = (1/2  h).(1/2  k).(1/2  l) ; N2 = (1/2 + h).(1/2  k).(1/2  l)
N3 = (1/2  h).(1/2 + k).(1/2  l) ; N4 = (1/2 + h).(1/2 + k).(1/2  l)
N5 = (1/2  h).(1/2  k).(1/2 + l) ; N6 = (1/2 + h).(1/2  k).(1/2 + l)
N7 = (1/2  h).(1/2 + k).(1/2 + l) ; N8 = (1/2 + h).(1/2 + k).(1/2 + l)
Global coordinates are related to local coordinates with general curved bricks
by the isoparametric transformation:
x = N(i).x(i) ; y = N(i).y(i) ; z = N(i).z(i)
The differential form of an isoparametric transformation is even more important
(here d*/d* are partial derivatives):
 dx   dx/dh dx/dk dx/dl  dh 
 dy  =  dy/dh dy/dk dy/dl  dk 
 dz   dz/dh dz/dk dz/dl  dl 
The inverse differential transform can be found by hand (or by using MAPLE).
Now here comes a piece of theory that cannot be generalized from 2D to 3D.
The Jacobian determinant, along with partial derivatives dT/d(x,y,z) to the
global coordinates, come with shapefunctions that are far more _complicated_
than the ones they are derived from. But, we do not consider them as relevant.
The only interesting case occurs when things are specified for the vertices of
the brick. We take vertex (1) as an example. Using MAPLE now: ^^^^^^^^
> with(linalg);
> N1 := (1/2  h)*(1/2  k)*(1/2  l);
> N2 := (1/2 + h)*(1/2  k)*(1/2  l);
> N3 := (1/2  h)*(1/2 + k)*(1/2  l);
> N4 := (1/2 + h)*(1/2 + k)*(1/2  l);
> N5 := (1/2  h)*(1/2  k)*(1/2 + l);
> N6 := (1/2 + h)*(1/2  k)*(1/2 + l);
> N7 := (1/2  h)*(1/2 + k)*(1/2 + l);
> N8 := (1/2 + h)*(1/2 + k)*(1/2 + l);
> x := x1*N1 + x2*N2 + x3*N3 + x4*N4 + x5*N5 + x6*N6 + x7*N7 + x8*N8;
> y := y1*N1 + y2*N2 + y3*N3 + y4*N4 + y5*N5 + y6*N6 + y7*N7 + y8*N8;
> z := z1*N1 + z2*N2 + z3*N3 + z4*N4 + z5*N5 + z6*N6 + z7*N7 + z8*N8;
> xh := diff(x,h) ; xk := diff(x,k) ; xl := diff(x,l) ;
> yh := diff(y,h) ; yk := diff(y,k) ; yl := diff(y,l) ;
> zh := diff(z,h) ; zk := diff(z,k) ; zl := diff(z,l) ;
> xh := subs(h=1/2, k=1/2, l=1/2 , xh);
> yh := subs(h=1/2, k=1/2, l=1/2 , yh);
> zh := subs(h=1/2, k=1/2, l=1/2 , zh);
> xk := subs(h=1/2, k=1/2, l=1/2 , xk);
> yk := subs(h=1/2, k=1/2, l=1/2 , yk);
> zk := subs(h=1/2, k=1/2, l=1/2 , zk);
> xl := subs(h=1/2, k=1/2, l=1/2 , xl);
> yl := subs(h=1/2, k=1/2, l=1/2 , yl);
> zl := subs(h=1/2, k=1/2, l=1/2 , zl);
> A := array([[ xh , xk , xl ],
> [ yh , yk , yl ],
> [ zh , zk , zl ]]);
[  x1 + x2  x1 + x3  x1 + x5 ]
[ ]
A := [  y1 + y2  y1 + y3  y1 + y5 ]
[ ]
[  z1 + z2  z1 + z3  z1 + z5 ]
> quit;
It is seen that the above matrix is the same as would have been obtained by
assuming a linear _tetrahedron_ with vertices (1,2,3,5) and origin (1) :
T = T1 + (T2  T1).h + (T3  T1).k + (T5  T1).l
dT/dh dT/dk dT/dl ; where T := x , y , z .
Consequently, a brick is just a composite of 8 linear tetrahedra:
(1,2,3,5) , (2,4,1,6) , (3,1,4,7) , (4,3,2,8)
(5,7,6,1) , (6,5,8,2) , (7,5,8,3) , (8,6,7,4)
with their origins at the vertices (nodes) of the brick.

* 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