sci.math.numanalysis
SUNA42, Discretization for 2D diffusion
========================================
Consider the twodimensional 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 chainrule 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 nondegenerate & 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 1D case: SUNA41. See also the book of Zienkiewicz.
A possible implementation is realized by the following BASIC programfragment:
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 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
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood