sci.math.numanalysis
SUNA56, 3D Potential Flow
==========================
Potential flow in three dimensions is described analytically by the equation of
Laplace:
[ (d/dx)^2 + (d/dy)^2 + (d/dz)^2 ] f = 0
Here x , y , z = space coordinates ; f = unknown function ; d = partial .
In order to arrive at a Finite Element formulation, the equation is weighted
weighted with a test function F (Galerkin method). Partial integration over
a space volume then results in:
///
  [ dF/dx.df/dx + dF/dy.df/dy + dF/dz.df/dz ].dx.dy.dz + : bulk
///
//
+  F.[ df/dx.nx + df/dy.ny + df/dz.nz ].dA = 0 : boundary
// _
n = (nx, ny, xz) = surface normal ; dA = surface area element
//
The boundary integral is also written as:  F.df/dn.dA = 0
//
We could refer to Zienkiewicz "The Finite Element Method" (3th edition) chapter
3.3 . The derivation there is a bit clumsy, however (especially formula 3.22).
But the following conclusion is right:
a) "natural" boundary conditions with df/dn = 0 are automatically satisfied.
Many ingredients for our own finite element formulation have been developed in
previous SUNA posters. Especially the theorem that overlapping tetrahedra in a
brick are equivalent with numerical integration at the vertices of that brick.
We have seen that the latter is an integration method which makes F.E.M. fully
equivalent with a Finite Volume Method (: SUNA50).
Let the unit cube be the parent element of the brick, and numbering according
to our well known convention:
number  1 = vertex coordinates (z,y,x)
selected from {0,1} and conceived as bits.
Then the tetrahedra under consideration are enumerated as follows:
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,8,5,3 , 8,6,7,4
                               
z 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0
y 0 0 1 0 0 1 0 0 1 0 1 1 1 1 0 1 0 1 0 0 0 0 1 0 1 1 0 1 1 0 1 1
x 0 1 0 0 1 1 0 1 0 0 1 0 1 0 1 1 0 0 1 0 1 0 1 1 0 1 0 0 1 1 0 1
Given the number of the brick's corner (1,...8), the other vertexnumbers of a
tetrahedron can be found by making them _one_ bit different from the first one.
For reasons of comparison, also the classical partitions of an eight cornered
brick into tetrahedra were taken into account (see Zienkiewicz' book par 6.3).
There exist two alternatives for a subdivision into five tetrahedra; and there
exist four alternatives for a subdivision into six. Care must be taken that in
the classical F.E. case everything should be according to ideology: the volume
of tetrahedra is 1/6 times a determinant, and there exist no diagonals crossed
at the brick faces. The latter means that the two types of subdivisions into 5
tetrahedra should be used in a mixed fashion: one type for uneven and the other
type for even numbered brick elements. This works for systematical numberings,
and if the total number of bricks in every mesh direction is uneven. For the 6
tetrahedra subdivisions, a mere enumeration of all vertex quadruples would have
been sufficient. Instead, this enumeration is calculated by a routine "coding",
which could do the same job for more than 3 dimensions (iff such ever occurs).
Last, but not least, also the classical trilinear interpolation with numerical
integration at (Gauss & other) points was implemented. The SUNA case above is
found back here for integration points positioned at the brick's corner nodes.
The SUNA case, together with the above classical tetrahedron partitions, are
implemented in the two programs, which are associated with this article, by a
variable called "MODE". As follows:
MODE Kind of brick element
 
0 known analytical solution
1 tetrahedra at the corners (SUNA)
2 composite of 5 tetrahedra (2 types)
3,4,5,6 different composites of 6 tetrahedra
7,8,9 numerical integration at Gauss points,
at vertices (SUNA) & close to midpoint.
The composite elements are implemented as "elem002.for", bricks with numerical
integration are implemented as "elem004.for". The Fortran statement function
kube(i,n) surely works also for hypercubes. The theory behind the elements is
handled by several SUNA posters, but is repeated here shortly for convenience.
For the case of a tetrahedron, the basic formula's are (with "d" partial):
df/dh = dx/dh.df/dx + dy/dh.df/dy + dz/dh.df/dz
df/dk = dx/dk.df/dx + dy/dk.df/dy + dz/dk.df/dz
df/dl = dx/dl.df/dx + dy/dl.df/dy + dz/dl.df/dz
here h,k,l = local element coordinates, x,y,z = global coordinates.
Discretizations for d(x,y,z,f)/d(h,k,l) can be found by elementary means.
In order to arrive at expressions for df/d(x,y,z) , the system of equations
must be inverted, resulting in well known formulas for d(h,k,l)/d(x,y,z) .
And:
df/dx = dh/dx.df/dh + dk/dx.df/dk + dl/dx.df/dl
df/dy = dh/dy.df/dh + dk/dy.df/dk + dl/dy.df/dl
df/dz = dh/dz.df/dh + dk/dz.df/dk + dl/dz.df/dl
Here df/d(h,k,l) = f2  f1 , f3  f1 , f4  f1 . Thus:
 df/dx   dh/dxdk/dxdl/dx dh/dx dk/dx dl/dx   f1 
 df/dy  =  dh/dydk/dydl/dy dh/dy dk/dy dl/dy   f2 
 df/dz   dh/dzdk/dzdl/dz dh/dz dk/dz dl/dz   f3 
  f4 
Approximation matrix
Multiplying the transpose of the Approximation matrix with itself results in
the 4 x 4 element matrix for a tetrahedron, describing a potential problem.
For numerically integrated brick elements, the problem is splitted up into 8
subproblems, which in fact means that an elementmatrix is constructed for
each of the integration points. Expressions d(x,y,z)/d(h,k,l) are specified
at an integration point before the Jacobian matrix is inverted. Derivatives
df/d(h,k,l) fall apart into coefficients for (f1, f2, f3, f4), resulting in
a slightly different procedure for calculating the 4 x 3 Approximation matrix
of trilinear bricks for the 3D potential problem.
The 9 elements for potential flow were tested for two different cases, where
also the analytical solution is well known (the first one thanks to Usenet):
a. Potential flow around a (unit) sphere
b. Potential flow around a (unit) cylinder
The analytical potential function for case (a) is given by: (r + 1/2r^2).z/r
The analytical potential function for case (b) is given by: (r + 1/r).x/r
Case (a) is implementented as the Fortran program "suna56a.for".
Case (b) is implementented as the Fortran program "suna56b.for".
The layout of the finite element mesh is calculated in place. Attention is paid
to a "fancy statement function", which actually relates Finite Difference grids
to Finite Element meshes. It is defined in "suna56b" as:
nr(i,j,k)=NDZ*NDT*(i1)+NDZ*(j1)+k
Meshing in case of the sphere, instead of the cylinder, is much more tricky,
especially at the North and South Poles (:). Also the accompanying statement
function is somewhat more complicated.
Basically, spherical coordinates are used for problem (a):
x = r.cos(theta).cos(phi)
y = r.cos(theta).sin(phi) with discrete values for (r, phi, theta).
z = r.sin(theta)
And, of course, cylinder (or polar) coordinates for (b):
x = r.cos(psi)
y = r.sin(psi) with discrete values for (r, psi, z).
z = z
Degenerate bricks are used where appropriate: at the pole region of the sphere
problem. Take a look at the "suna56a.inp" AVS UCD output to see what I mean.
Explicit boundary conditions are:
 Potential = 0 for one of the boundaries. Due to this B.C. the numerical
solution would be very boring if we used more than just _one_ element in
the direction parallel to that boundary.
 Velocity = 1 for one of it's components, at "infinity" (outer boundary).
Due to this B.C. the potential along the outer boundary can be calculated
as a function of one of the coordinates.
Implicit boundary conditions are given by the fact that the normal derivative
of the potential is zero, if nothing else is done.
A comparison is made between the numerical solutions and the exact one. Well,
don't shoot at me for that! Look at the outcomes yourself, at "suna56*.out".
Especially in case of the sphere potential problem, for certain elements, the
result is ultimately dissapointing. I think degeneration of bricks in certain
regions of the mesh is responsible for it. As expected, the results for MODE=
1 are identical to those for MODE=8 . The SUNA element is among most reliable
elements. It seems to behave well even in the degenerate case (: where some of
it's tetrahedra are simply eliminated for "det=0". OK, so that's the trick ;).
A piece of postprocessing is done also, in order to complete the AVS UCD file
(which is meant for a course I'm preparing) with "cell data". Three equations
are set up, according to the theory of "Stream Function & Flow Potential" (14),
when extended to three dimensions. All three are of the form:
P2  P1 = (x2  x1).u + (y2  y1).v + (z2  z1).w
In order to have a reasonable condition of the 3 equations, opposite midpoints
of the brick's faces are used for evaluating the potential function values.
This kind of "numerical differentiation" seems to work very well.

* 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