sci.math.num-analysis

SUNA02, Triangle Isoparametrics
===============================

Isoparametric transformation is the standard approach where the Finite Element
Method relies on when it has to deal with curved geometries: it's strong point.
But let's consider first the simplest non-trivial finite element shape in two
dimensions, the linear triangle, and see where isoparametric transformations
                       do come from. Function behaviour is approximated inside
          3            such a triangle by a _linear_ interpolation between the
        /   \          function values at the vertices, or nodal points.
      /       \        Let  T  be such a function, and  x,y  coordinates,
    /           \      then:
  1 ------------- 2               T = A.x + B.y + C

Where the constants  A, B, C  are determined by:

      T1 = A.x1 + B.y1 + C
      T2 = ...

But wait. The first of these equations can already be used to eliminate the
constant  C , once and forever:

      T - T1 = A.(x - x1) + B.(y - y1)

Then the constants  A  and  B  are determined by:

      T2 - T1 = A.(x2 - x1) + B.(y2 - y1)
      T3 - T1 = A.(x3 - x1) + B.(y3 - y1)

Two equations with two unknowns. The solution is found by Cramer's rule:

      A = [ (y3 - y1).(T2 - T1) - (y2 - y1).(T3 - T1) ]/J
      B = [ (x2 - x1).(T3 - T1) - (x3 - x1).(T2 - T1) ]/J

Where:  J = (x2 - x1).(y3 - y1) - (x3 - x1).(y2 - y1) = 2 * triangle area .
So:
         T - T1 = h.(T2 - T1) + k.(T3 - T1)
Where:
              h = [ (y3 - y1).(x - x1) - (x3 - x1).(y - y1) ]/J
              k = [ (x2 - x1).(y - y1) - (y2 - y1).(x - x1) ]/J
Or:
          | h |   | + (y3 - y1)    - (x3 - x1) |   | x - x1 |
          |   | = |                            |/J |        |
          | k |   | - (y2 -y1)     + (x2 - x1) |   | y - y1 |

The inverse of the following problem is recognized herein:

     | x - x1 |   | x2 - x1     x3 - x1 | | h |
     |        | = |                     | |   |
     | y - y1 |   | y2 - y1     y3 - y1 | | k |
Or:
        x - x1 = h.(x2 - x1) + k.(x3 - x1)
        y - y1 = h.(y2 - y1) + k.(y3 - y1)
But:
        T - T1 = h.(T2 - T1) + k.(T3 - T1)

Typical! The _same_ expression holds for the function  T  as well as for the
coordinates x and y . Which is precisely what people mean by an ISOPARAMETRIC
("same parameters") transformation.

Now recall the formulas which express h and k into x and y:

    h = [ (y3 - y1).(x - x1) - (x3 - x1).(y - y1) ]/J
    k = [ (x2 - x1).(y - y1) - (y2 - y1).(x - x1) ]/J

Thus h can be interpreted as: area of the sub-triangle spanned by the vectors
(x - x1 , y - y1) and (x3 - x1 , y3 - y1) divided by the whole triangle area.
And k can be interpreted as: area of the sub-triangle spanned by the vectors
(x - x1 , y - y1) and (x2 - x1 , y2 - y1) divided by the whole triangle area.
This is the reason why  h  and  k  are sometimes called _area-coordinates_.

There are even _three_ of these coordinates in litterature (: Zienkiewicz).
But one of them can be allways be safely discarded. I think it's a bad habit
to retain a "coordinate" like (1 - h - k), being dependent on the other two.

Instead of area-coordinates, we prefer to talk about _local coordinates_
h and k  of an element, in contrast to the _global coordinates_  x and y .

 3              It is possible that local coordinates coincide with the
 |\             global coordinates. A triangle for which this is the case
 |  \           is called a _parent element_. The portrait of the parent
 |    \         triangle is displayed here: rectangular, two sides equal.
 |______\
 1       2

Let's reconsider for a moment the expression:

     T - T1 = h.(T2 - T1) + k.(T3 - T1)

Partial differentiation to  h and k  gives:

        dT/dh = T2 - T1
        dT/dk = T3 - T1
So:
        T = T1 + h.dT/dh + k.dT/dk

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).T1 + h.T2 + k.T3

Here the functions  (1 - h - k) , h , k  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

So we have two representations, which are allmost trivially equivalent:

     T = T1 + h.(T2 - T1) + k.(T3 - T1)      : Finite Difference like
     T = (1 - h - k).T1 + h.T2 + k.T3        : Finite Element like

The two sets of functions (1,h,k) and (1-h-k,h,k) are mutually related as
follows:
             F.D. <- F.E.                             F.E <- F.D.
    | 1 |   | 1   1   1 | | 1-h-k |      | 1-h-k |   | 1  -1  -1 | | 1 |
    | h | = | 0   1   0 | | h     |      | h     | = | 0   1   0 | | h |
    | k |   | 0   0   1 | | k     |      | k     |   | 0   0   1 | | k |

The two sets of functions can be conceived as two coordinate systems in an
abstract three dimensional space. Then the above matrices are _transition_
matrices, belonging to F.D. <-> F.E. base transitions. Yes, the above is an
instance of the Unification Principle (Re: SUNA, Manifesto).

Not too difficult, eh? Well, not yet. To be continued ...
-
* 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