sci.math.numanalysis
SUNA49, 3D Theorem Preliminaries
=================================
This poster is a generalization of the 2D theory in SUNA11 to three dimensions.
Consider the following partial differential equation:
dQx/dx + dQy/dy +dQz/dz = 0 ("d" denotes partial differentiation)
Here (Qx,Qy,Qz) can be interpreted as the heatflow in a 3dimensional medium,
or as the velocities (u,v,w) in a 3D incompressible flow field.
According to the Finite Domain technique, a kind of Finite Difference Method,
this equation is first integrated over a certain control volume:
///
 (dQx/dx + dQy/dy + dQz/dz) dx.dy.dz = 0
///
Partial integration (Gauss' theorem) results in an equivalent formulation,
which contains a surfaceintegral:
//
OO (Qx.dAx + Qy.dAy + Qz.dAz) where: (dAx,dAy,dAz) = infinitesimal area
// * normal vector
According to the (Galerkin) Finite Element Method, the conservation equation
must first be multiplied by a weighting function F , and then integrated over
a finite element volume, according to:
///
 F.(dQx/dx + dQy/dy + dQ/dz) dx.dy.dz = 0
///
Partial integration results in an expression with surfaceintegrals, which all
become part of the boundary conditions, and a term for the bulk material:
///
  (dF/dx.Qx + dF/dy.Qy + dF/dz) dx.dy.dz (: mind the minus sign)
///
For the sake of generality, both the Finite Volume and the Galerkin formulation
are assumed to be applicable at a curvilinear grid, consisting of brick shaped
cells or elements.
For the sake of simplicity, we wish to restrict ourselves to regular molecules
and elements, such as the classical nine point star, and cubic bricks.
The general curvilinear grid is mapped upon a regular grid by a transformation
of coordinates x(h,k.l) , y(h,k,l) , z(h,k,l) in such a way that the _whole_
problem may be assumed to be regular in the (h,k,l) plane. This is the Finite
Volume point of view, essentially a _global_ transformation.
From a Finite Element point of view, however, such coordinate transformations
are always carried out _locally_/elementwise: by isoparametrics. Restrictions
of orthogonality do not apply in a Finite Element context.
In either case, integrands of the Finite Volume and the Galerkin formulations
must be transformed first into such regular (h,k,l) coordinates.
The vector (Qx,Qy,Qz) is expressed in components at the (h,k,l) coordinate
system as follows:
 Qx   dx/dh   dx/dk   dx/dl   xh xk xl   Qh 
 Qy  = Qh. dy/dh  + Qk. dy/dk  + Ql. dy/dl  =  yh yk yl   Qk 
 Qz   dz/dh   dz/dk   dz/dl   zh zk zl   Ql 
Here we made use of the following substitutions:
xh = dx/dh ; xk = dx/dk ; xl = dx/dl
yh = dy/dh ; yk = dy/dk ; yl = dy/dl
zh = dz/dh ; zk = dz/dk ; zl = dz/dl
Differentials of the global coordinates are transformed as follows:
 dx   dx/dh dx/dk dx/dl   dh  (with proper interpretation of
 dy  =  dy/dh dy/dk dy/dl   dk  partial and nonpartial "d's")
 dz   dz/dh dz/dk dz/dl   dl 
Last but not least, the chain rules for dT(h,k,l) are just the inverse of:
 dT/dx   dx/dh dy/dh dz/dh  1  dT/dh 
 dT/dy  =  dx/dk dy/dk dz/dk   dT/dk 
 dT/dz   dx/dl dy/dl dz/dl   dT/dl 
Parametrize the top and bottom surfaces of the finite volume by h and k .
Then the two local base vectors at these surfaces are given by:
(dx/dh, dy/dh, dz/dh) and (dx/dk, dy/dk, dz/dk)
The normal at the surfaces is given by +/ the outer product of these vectors,
and the components of the infinitesimal area's (dAx,dAy,dAz) are dh.dk times
the normal:
 dAx   dy/dh.dz/dk  dz/dh.dy/dk 
 dAy  =  dz/dh.dx/dk  dz/dk.dx/dh  .dh.dk
 dAz   dx/dh.dy/dk  dx/dk.dy/dh 
The integrand of the Finite Volume formulation can now be written as follows,
for one of the face surfaces:
 Qx Qy Qz   dAx 
 dAy  =
 dAz 
 Qh Qk Ql   dx/dh dy/dh dz/dh   dy/dh.dz/dk  dz/dh.dy/dk 
 dx/dk dy/dk dz/dk   dz/dh.dx/dk  dz/dk.dx/dh  .dh.dk
 dx/dl dy/dl dz/dl   dx/dh.dy/dk  dx/dk.dy/dh 
With MAPLE, we will calculate the inverse of the above 3 x 3 matrix:
> with(linalg);
> A := array([[ xh, yh, zh ],
> [ xk, yk, zk ],
> [ xl, yl, zl ]]);
> inverse(A);
[ yk zl  zk yl yh zl  zh yl yh zk  zh yk ]
[     ]
[ %1 %1 %1 ]
[ ]
[ xk zl  zk xl xh zl  zh xl xh zk  zh xk ]
[      ]
[ %1 %1 %1 ]
[ ]
[ xk yl  yk xl xh yl  yh xl xh yk  yh xk ]
[     ]
[ %1 %1 %1 ]
%1 := xh yk zl  xh zk yl  xk yh zl + xk zh yl + xl yh zk  xl zh yk
# The determinant of the original matrix is recognized in expression %1 .
# It's precisely the Jacobian determinant J , of course.
> quit;
The last column in MAPLE's result gives us the thing we need. Because It proves
that the product of the Jacobian matrix and the outproduct vector must be equal
to the last column of the unity matrix times the Jacobian determinant:
 Qx Qy Qz   dAx   Qh Qk Ql   0 
 dAy  =  0  . dh.dk = Ql.J.dh.dk
 dAz   J 
Which gives a remarkebly simple result:
 //
 OO (Qx.dAx + Qy.dAy + Qz.dAz) =
 //
 // // //
  Qh.J.dk.dl +  Qk.J.dh.dl +  Ql.J.dh.dk +  surface integral
 // // //  over all six faces
 // // //  of the unit cube,
   Qh.J.dk.dl   Qk.J.dh.dl   Ql.J.dh.dk .  parent of the brick
 // // //
Same kind of procedure for the Galerkin integral:
{ dF/dx.Qx + dF/dy.Qy + dF/dz.Qz}.dx.dy.dz =
[ dF/dh dF/dk dF/dl ]  xh xk xl  1  xh xk xl   Qh 
 yh yk yl  .  yh yk yl   Qk  .J.dh.dk.dl
 zh zk zl   zh zk zl   Ql 
Same kind of remarkebly simple result:
 ///
   (dF/dx.Qx + dF/dy.Qy + dF/dz.Qz) dx.dy.dz =
 ///
 ///
   J. (dF/dh.Qh + dF/dk.Qk + dF/dl.Ql) dh.dk.dl
 ///
In order to shorten notation, we shall make everywhere the replacement:
Qh := J.Qh ; Qk := J.Qk ; Ql := J.Ql

* 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