sci.math.num-analysis

SUNA50, 3-D Unification Theorem
===============================

Using a F.D. method, the grid is seen to be built up from *seven* node stars,
where the central node is enclosed by a control volume. The surface-integral
has to be applied to the boundary of this control volume, which is rectangular.
Employing a F.E. viewpoint, on the contrary, the same grid is seen to be built
up from bricks. The Galerkin integral will be evaluated for each of these brick
elements (rectangular prisms or whatever) in the mesh.

Switching again to the F.D. viewpoint, all the elements contain eight pieces of
a finite domain control volume, each piece belonging to one-eighth of a nodal
point. If we restrict ourselves to unit coordinates, then a F.E. brick is given
by its vertices as:

                     (1) = (-1/2,-1/2,-1/2)   ;  (2) = (+1/2,-1/2,-1/2)
And: dh = 1          (3) = (-1/2,+1/2,-1/2)   ;  (4) = (+1/2,+1/2,-1/2)
     dk = 1          (5) = (-1/2,-1/2,+1/2)   ;  (6) = (+1/2,-1/2,+1/2)
     dl = 1          (7) = (-1/2,+1/2,+1/2)   ;  (8) = (+1/2,+1/2,+1/2)

As a consequence, there are also eight pieces wherein the Finite Domain integral
is splitted up:
                    //
                   || (Qh.dk.dl + Qk.dh.dl + Ql.dh.dk) = 0
                  //

After a little algebra, these area-integrals correspond to the following vector
of necessarely incomplete, equations:

       | - (Qh1 + Qh2) - (Qk1 + Qk3) - (Ql1 + Ql5) |
       | + (Qh2 + Qh1) - (Qk2 + Qk4) - (Ql2 + Ql6) |
       | - (Qh3 + Qh4) + (Qk3 + Qk1) - (Ql3 + Ql7) |
 - 1/8 | + (Qh4 + Qh3) + (Qk4 + Qk2) - (Ql4 + Ql8) | = 0
       | - (Qh5 + Qh6) - (Qk5 + Qk7) + (Ql5 + Ql1) |
       | + (Qh6 + Qh5) - (Qk6 + Qk8) + (Ql6 + Ql2) |
       | - (Qh7 + Qh8) + (Qk7 + Qk5) + (Ql7 + Ql3) |
       | + (Qh8 + Qh7) + (Qk8 + Qk6) + (Ql8 + Ql4) |

On the other hand, for the same brick element, there exists a Galerkin integral.
The shape functions of a brick (cube) element are:

 Ni(h,k,l) = (1/2 % h)(1/2 % k)(1/2 % l)   where % is + or - , depending on i .

An arbitrary function  F  at the cube will be interpolated by these polynomials.
The partial derivatives of  F  are calculated and (since we make use of the 2-D
heuristics) they are immediately specified for the vertices of the brick:

                   1,2:     3,4:     5,6:     7,8:
       dF/dh = ( F2 - F1, F4 - F3, F6 - F5, F8 - F7)
                   1,3:     2,4:     5,7:     6,8:
       dF/dk = ( F3 - F1, F4 - F2, F7 - F5, F8 - F6)
                   1,5:     2,6:     3,7:     4,8:
       dF/dl = ( F5 - F1, F6 - F2, F7 - F3, F8 - F4)

Numerical integration will be employed. Eight integration points are defined,
which are precisely at the vertices of the brick, analogous to the 2-D case:

    dF/dh.Qh + dF/dk.Qk + dF/dl.Ql .dh.dk.dl = 1/8 *

    (F2 - F1).(Qh2 + Qh1) + (F3 - F1).(Qk3 + Qk1) + (F5 - F1).(Ql5 + Ql1) +
    (F4 - F3).(Qh4 + Qh3) + (F4 - F2).(Qk4 + Qk2) + (F6 - F2).(Ql6 + Ql2) +
    (F6 - F5).(Qh6 + Qh5) + (F7 - F5).(Qk7 + Qk5) + (F7 - F3).(Ql7 + Ql3) +
    (F8 - F7).(Qh8 + Qh7) + (F8 - F6).(Qk8 + Qk6) + (F8 - F4).(Ql8 + Ql4)

By partial differentiation to the nodal values of  F , the Galerkin integral
reaches its minimum value. A vector of equations is then associated which each
cubic element, finally resulting in what is called the local "stiffness" matrix
and "load" vector of the element considered.

This results in a set of, necessarily incomplete, equations. This system turns
out to be IDENTICAL to the one derived above for the Finite Difference approach.
Again this is a *unified* F.D.-F.E. discretization scheme for the heat fluxes.

It is well known that, in the case of thermal diffusion, the Galerkin integral
turns out to be equivalent to a variational principle:

    ///
   ||| Kx.(dT/dx)^2 + Ky.(dT/dy)^2 + Kz.(dT/dz)^2 .dx.dy.dz = minimum
  ///

At this point we introduce curvilinear global coordinates (x,y,z). The global
derivatives to (x,y,z) are expressed in their local counterparts as follows:

  | 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 |

The elements of the inverse matrix are *all* expressed in  d(x,y,z)/d(h,k,l) .
Now an isoparametric transformation is characterized by the mere fact that all
approximations, including those of the global coordinates x,y,z , are expressed
in the same uniform manner. Hence the local derivatives of  T , x , y , z  will
have the same appearance. These derivatives, however, can also be interpreted
as approximations in _linear tetrahedra_, as we have seen several times before.
These tetrahedra are situated with their origins (1) at the vertex of each F.E.
brick element.

It is concluded that the Finite Difference equivalent is in fact applicable to a
mesh of bricks, which are in turn subdivided into eight overlapping tetrahedra.
Using a local numbering for each brick, these tetrahedra are:

   (1,2,3,5) , (2,1,4,6) , (3,4,1,7) , (4,3,2,8)
   (5,6,7,1) , (6,5,8,2) , (7,8,5,3) , (8,7,6,4)

It is remarked that the tetrahedra do *not* fill up the whole brick element.
This becomes more clear if we split up the above collection of tetrahedra into
the following two parts, each of which constitute a well known subdivision in
"The Finite Element Method" (O.C. Zienkiewicz, paragraph 6.3):

    (1,2,3,5) , (4,3,2,8) , (6,5,8,2) , (7,8,5,3) , missing: (2,3,5,8)
    (2,1,4,6) , (3,4,1,7) , (5,6,7,1) , (8,7,6,4) , missing: (1,4,6,7)

Each of the subdivisions would occupy the whole space, if nothing was "missing".

Let us resume now the result, while omitting the heuristic intermediate steps:
a Finite Element - Finite Difference Unification Theorem, for 3-D Diffusion.

Theorem:
-------
| Consider the diffusion problem (Laplace's equation) in three dimensions.
| Let this problem be formulated for the general curvilinear case.
| Then the (most common) Finite Difference (domain) Method for it turns out to
| be (pretty) equivalent to the (most common) Finite Element (Galerkin) Method,
| provided the latter is carried out in a rather unusual way, namely for two
| mutually *overlapping* meshes of (linear) tetrahedra: the origin vertices of
| these tetrahedra must be coincident with the vertices of the brick elements
| they are to be associated with. As a consequence, (each of) the above two
| meshes do *not* fill up the whole 3-D space.

Using 3-D solid elements, there seems to be a problem with overstiff equations
in F.E. applications (if I am well informed concerning this).
Since there is no such a rumour heard from the Finite Difference world, I have
the guts to make the following

Prediction
----------
| In order to cure overstiffness in 3-D Finite Element applications, one should
| discretize with brick elements, subdivided into *four* tetrahedra, instead of
| the five that would fill up the whole space. This means that the tetrahedron
| in the middle must be deliberately be "missing", or EMPTY, as indicated above.

Disclaimer: I am myself not in the position to carry out such a test! Because
            these posters were created at my own Personal Computing Centre, at
            home in Barendrecht. They have nothing to do with my job in Delft.
-
* 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
* E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood