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