sci.math.num-analysis

SUNA42, Discretization for 2-D diffusion
========================================

Consider the two-dimensional Laplace equation, which describes for example heat
conduction in a metal plate:
                                 [ (d/dx)^2 + (d/dy)^2 ] f = 0

Here the ASCII "d's" denote partial differentiation, as usual.

The variational principle, associated with the Laplace equation, is well known:

     //
    || [ (df/dx)^2 + (df/dy)^2 ] dx.dy = minimum
   //

If the integration domain is subdivided into finite elements E, then the above
integral is plitted up as a sum of integrals over these elements:

         //
   Sum  || [ (dT/dx)^2 + (dT/dy)^2 ] dx.dy = minimum
    E  //

Take linear triangles as the finite elements, then function behaviour inside
each triangle is like:
                          f = f1 + (f2 - f1).h + (f3 - f1).k
Isoparametrics:           x = x1 + (x2 - x1).h + (x3 - x1).k
                          y = y1 + (y2 - y1).h + (y3 - y1).k

Here (h,k) are the independent local (area) coordinates of a triangle.

Let's work out, according to the chain-rule for partial derivatives:

  df/dh = df/dx.dx/dh + df/dy.dy/dh
  df/dk = df/dx.dx/dk + df/dy.dy/dk

  f2 - f1 = (x2 - x1).df/dx + (y2 - y1).df/dy
  f3 - f1 = (x3 - x1).df/dx + (y3 - y1).df/dy

Define: f21 = f2 - f1 ; f31 = f3 - f1 ; likewise for x21, x31, y21, y31. Then:

  df/dx = (f21.y31 - f31.y21) / (x21.y31 - x31.y21)
  df/dy = (f21.x31 - f31.x21) / (y21.x31 - y31.x21)

The determinant is recognized:  Det = x21.y31 - x31.y21 = 2 * triangle area .
Assume all triangles are non-degenerate & oriented counterclockwise:  Det > 0 .

Rewrite as follows:

  df/dx = [ (y21 - y31).f1 + y31.f2 - y21.f3 ] / Det
  df/dy = [ (x31 - x21).f1 - x31.f2 + x21.f3 ] / Det

Or even neater, with help of matrix notation:

  df/dx = [ f1 f2 f3 ] | y23 |            df/dy = - [ f1 f2 f3 ] | x23 |
                       | y31 | / Det                             | x31 | / Det
                       | y12 |                                   | x12 |

Consequently we find for the integrand belonging to one arbitrary triangle:

  (df/dx)^2 + (df/dy)^2 = constant = 1/Det^2 *

  (y23.f1 + y31.f2 + y12.f3)^2 + (x23.f1 + x31.f2 + x12.f3)^2

Last but not least:    dx = (x2 - x1).dh + (x3 - x1).dk
                       dy = (y2 - y1).dh + (y3 - y1).dk

Consequently:   dx.dy = Det | x21  x31 | . dh.dk
                            | y21  y31 |

      //               //
     || dx.dy = Det . || dh.dk = (x21.y31 - x31.y21) . 1/2
    //               //

The minimum of a variational integral is obtained by (partial) differentiation
to each of the quantities  f , as discretized in the nodal points: f1, f2, ...
This results in the standard general finite element assembly procedure, as has
been demonstrated in the 1-D case: SUNA41. See also the book of Zienkiewicz.
A possible implementation is realized by the following BASIC program-fragment:

240 PRINT "Assembly procedure for";LIMT;"triangles with 1 d.o.f. at each node"
250 REM ======================================================================
260 DIM NR(3),E(3,3)
265 DIM S(NN,NB1),B(NN) : REM NN = number of unknowns , NB1 = bandwidth + 1
266 REM
268 REM Clear global matrix and vector:
270 FOR I=1 TO NN : FOR J=1 TO NB1 : S(I,J)=0 : NEXT J : B(I)=0 : NEXT I
280 REM
290 OPEN "R",#1,"TOPOL"+C$+".BIN",12 : REM Connectivity file
300 FIELD #1, 4 AS F1$, 4 AS F2$, 4 AS F3$
310 FOR F=1 TO LIMT : GET #1,F : PRINT F;
320 NR(1)=CVI(F1$) : NR(2)=CVI(F2$) : NR(3)=CVI(F3$)
330 REM
340 GOSUB 600 : REM Form element matrix
345 REM
350 REM Account for bulk elements:
360 FOR I=1 TO 3 : IG=NR(I)
370 FOR J=I TO 3 : JG=NR(J)
380 II=IG : IF JG<IG THEN II=JG : REM Symmetric
390 JJ=ABS(IG-JG)+1             : REM storage mode
400 S(II,JJ)=S(II,JJ)+E(I,J)
410 NEXT J : NEXT I
440 NEXT F : CLOSE #1 : PRINT

The element-matrix, to be calculated by "GOSUB 600", is defined in our case by:

    [ f1 f2 f3 ] | E11 E12 E13 | | f1 |
                 | E21 E22 E23 | | f2 | = [ (df/dx)^2 + (df/dy)^2 ] * Det/2
                 | E31 E32 E33 | | f3 |    variational integral at triangle

Working out:  [ E(i,j) ] = 1/Det *

  | y23 |                       | x23 |
  | y31 | [ y23  y31  y12 ]  +  | x31 | [ x23  x31  x12 ]  =
  | y12 |                       | x12 |

  | y23.y23 + x23.x23   y23.y31 + x23.x31   y23.y12 + x23.x12 |
  |                     y31.y31 + x31.x31   y31.y12 + x31.x12 | / Det
  | symmetrical                             y21.y21 + x21.x21 |

The Determinants are essential and cannot be weighted out (as with L.S.FEM).

A possible interpretation of the matrix elements is:

  E23 = - (x3 - x1).(x2 - x1) - (y3 - y1).(y2 - y1) =
      = - (r31.r21)  apart from  1/Det
           ^   ^
With other words: any matrix element is (+ or -) the inner product of vectors
                  pointing from one vertex to (one or two) other vertices.

It is remarked that, according to one of Patankar's "Five Basic Rules", namely
"the rule of positive coefficients", off-diagonal terms like  E23  must be less
than zero. This is only the case for positive inner products (r31.r21) .

Which implies that any angle of a Diffusion Triangle must be  <= 90  degrees.

Special case: angle between vectors  r21 and r31  =  90 degrees -->  E23 = 0.
This is the finite difference case on a rectangular grid: "no" triangles, but
five point stars "instead".
-
* 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