sci.math.num-analysis SUNA47, Isoparametric brick =========================== Let's consider next the simplest-but-one element shape in 3-D: a brick. It is assumed that the parent element is a UNIT CUBE. Take the origin 0 of the local coordinate system in the middle of the cube: (h,k,l) = (1/2,1/2,1/2) . The cube's vertices are numbered as summarized in the following BASIC program: 5 DIM COORD(8,3) 10 FOR I=1 TO 8 : PRINT I , : FOR N=1 TO 3 20 COORD(I,N)=.5*( 2*( ((I-1)\2^(N-1)) MOD 2 ) -1) 30 PRINT COORD(I,N), : NEXT N : PRINT : NEXT I Output: 1 -.5 -.5 -.5 2 .5 -.5 -.5 3 -.5 .5 -.5 4 .5 .5 -.5 5 -.5 -.5 .5 6 .5 -.5 .5 7 -.5 .5 .5 8 .5 .5 .5 nr x y z Function behaviour may be approximated inside such a cube by a Taylor expansion around the origin. Since there are eight vertices, 8 terms must be retained: T = T(0) + dT(0)/dh.h + dT(0)/dk.k + dT(0)/dl.l 2 2 2 + d T(0)/dh.dk .h.k + d T(0)/dh.dl .h.l + d T(0)/dk.dl .k.l 3 + d T(0)/dh.dk.dl .h.k.l The interpolation between the function values at the nodes is thus tri-linear. Let's assume from the beginning that the same expression holds for the function T as well as for the coordinates x and y : an isoparametric transformation. Choose abbreviations for the names for the partial derivatives: T = A1 + A2.h + A3.k + A4.l + A5.h.k + A6.h.l + A7.k.l + A8.h.k.l Next specify the above expression for each of the vertices. This results in: | T1 | | 1 -1/2 -1/2 -1/2 +1/4 +1/4 +1/4 -1/8 | | A1 | | T2 | | 1 +1/2 -1/2 -1/2 -1/4 -1/4 +1/4 +1/8 | | A2 | | T3 | | 1 -1/2 +1/2 -1/2 -1/4 +1/4 -1/4 +1/8 | | A3 | | T4 | = | 1 +1/2 +1/2 -1/2 +1/4 -1/4 -1/4 -1/8 | | A4 | F.E. <- F.D. | T5 | | 1 -1/2 -1/2 +1/2 +1/4 -1/4 -1/4 +1/8 | | A5 | transition | T6 | | 1 +1/2 -1/2 +1/2 -1/4 +1/4 -1/4 -1/8 | | A6 | | T7 | | 1 -1/2 +1/2 +1/2 -1/4 -1/4 +1/4 -1/8 | | A7 | | T8 | | 1 +1/2 +1/2 +1/2 +1/4 +1/4 +1/4 +1/8 | | A8 | ------------------------------------------- 1 h k l h h k h k l l k l Before making any attempt to solve these equations, it should be remarked that all of the matrix columns are _mutually orthogonal_. This is easy to check out. At first it is seen that all columns sum up to zero: as many + as - signs. The inner product of the first and the other columns are just the other columns themselves. The inner products (h.k) , (h.l) , (k.l) are also other columns of the matrix. There remains to investigate cases like h.h.k , a positive column h.h times as many + as - signs in k . Same argument for things like h.h.k.l and h.h.k.k.l . Q.E.D. The inverse of an orthogonal matrix is the same as it's transpose. There is one little problem, however: the lengths of the matrix columns are not equal to 1 . But the latter is a rather simple matter to deal with: divide each row of the transposed matrix by its (length)^2 . The result is: | A1 | | +1/8 +1/8 +1/8 +1/8 +1/8 +1/8 +1/8 +1/8 | | T1 | | A2 | | -1/4 +1/4 -1/4 +1/4 -1/4 +1/4 -1/4 +1/4 | | T2 | | A3 | | -1/4 -1/4 +1/4 +1/4 -1/4 -1/4 +1/4 +1/4 | | T3 | | A4 | = | -1/4 -1/4 -1/4 -1/4 +1/4 +1/4 +1/4 +1/4 | | T4 | F.E. <- F.D. | A5 | | +1/2 -1/2 -1/2 +1/2 +1/2 -1/2 -1/2 +1/2 | | T5 | transition | A6 | | +1/2 -1/2 +1/2 -1/2 -1/2 +1/2 -1/2 +1/2 | | T6 | | A7 | | +1/2 +1/2 -1/2 -1/2 -1/2 -1/2 +1/2 +1/2 | | T7 | | A8 | | -1 +1 +1 -1 +1 -1 -1 +1 | | T8 | ------------------------------------------------ matrix C So we find explicit discretizations for the partial derivatives of the function T at the the origin 0 . The zero'th derivative is, for example: T(0) = 1/8.(+ T1 + T2 + T3 + T4 + T5 + T6 + T7 + T8) Let's further abbreviate by introducing the functions L(i): L1 L2 L3 L4 L5 L6 L7 L8 1 h k l h.k h.l k.l h.k.l And adopt the Einstein convention, which is summation over equal indices. Then: T = L(i).A(i) = L(i).C(i,j).T(j) Shape functions are defined by: T = N(j).T(j) Consequently: N(j) = L(i).C(i,j) . We can work it out. For example: N1 = + 1/8 - 1/4.(h + k + l) + 1/2.(h.k + h.l + k.l) - h.k.l Guess what: N1 = (1/2 - h).(1/2 - k).(1/2 - l) ; N2 = (1/2 + h).(1/2 - k).(1/2 - l) N3 = (1/2 - h).(1/2 + k).(1/2 - l) ; N4 = (1/2 + h).(1/2 + k).(1/2 - l) N5 = (1/2 - h).(1/2 - k).(1/2 + l) ; N6 = (1/2 + h).(1/2 - k).(1/2 + l) N7 = (1/2 - h).(1/2 + k).(1/2 + l) ; N8 = (1/2 + h).(1/2 + k).(1/2 + l) Global coordinates are related to local coordinates with general curved bricks by the isoparametric transformation: x = N(i).x(i) ; y = N(i).y(i) ; z = N(i).z(i) The differential form of an isoparametric transformation is even more important (here d*/d* are partial derivatives): | dx | | dx/dh dx/dk dx/dl | dh | | dy | = | dy/dh dy/dk dy/dl | dk | | dz | | dz/dh dz/dk dz/dl | dl | The inverse differential transform can be found by hand (or by using MAPLE). Now here comes a piece of theory that cannot be generalized from 2-D to 3-D. The Jacobian determinant, along with partial derivatives dT/d(x,y,z) to the global coordinates, come with shape-functions that are far more _complicated_ than the ones they are derived from. But, we do not consider them as relevant. The only interesting case occurs when things are specified for the vertices of the brick. We take vertex (1) as an example. Using MAPLE now: ^^^^^^^^ > with(linalg); > N1 := (1/2 - h)*(1/2 - k)*(1/2 - l); > N2 := (1/2 + h)*(1/2 - k)*(1/2 - l); > N3 := (1/2 - h)*(1/2 + k)*(1/2 - l); > N4 := (1/2 + h)*(1/2 + k)*(1/2 - l); > N5 := (1/2 - h)*(1/2 - k)*(1/2 + l); > N6 := (1/2 + h)*(1/2 - k)*(1/2 + l); > N7 := (1/2 - h)*(1/2 + k)*(1/2 + l); > N8 := (1/2 + h)*(1/2 + k)*(1/2 + l); > x := x1*N1 + x2*N2 + x3*N3 + x4*N4 + x5*N5 + x6*N6 + x7*N7 + x8*N8; > y := y1*N1 + y2*N2 + y3*N3 + y4*N4 + y5*N5 + y6*N6 + y7*N7 + y8*N8; > z := z1*N1 + z2*N2 + z3*N3 + z4*N4 + z5*N5 + z6*N6 + z7*N7 + z8*N8; > xh := diff(x,h) ; xk := diff(x,k) ; xl := diff(x,l) ; > yh := diff(y,h) ; yk := diff(y,k) ; yl := diff(y,l) ; > zh := diff(z,h) ; zk := diff(z,k) ; zl := diff(z,l) ; > xh := subs(h=-1/2, k=-1/2, l=-1/2 , xh); > yh := subs(h=-1/2, k=-1/2, l=-1/2 , yh); > zh := subs(h=-1/2, k=-1/2, l=-1/2 , zh); > xk := subs(h=-1/2, k=-1/2, l=-1/2 , xk); > yk := subs(h=-1/2, k=-1/2, l=-1/2 , yk); > zk := subs(h=-1/2, k=-1/2, l=-1/2 , zk); > xl := subs(h=-1/2, k=-1/2, l=-1/2 , xl); > yl := subs(h=-1/2, k=-1/2, l=-1/2 , yl); > zl := subs(h=-1/2, k=-1/2, l=-1/2 , zl); > A := array([[ xh , xk , xl ], > [ yh , yk , yl ], > [ zh , zk , zl ]]); [ - x1 + x2 - x1 + x3 - x1 + x5 ] [ ] A := [ - y1 + y2 - y1 + y3 - y1 + y5 ] [ ] [ - z1 + z2 - z1 + z3 - z1 + z5 ] > quit; It is seen that the above matrix is the same as would have been obtained by assuming a linear _tetrahedron_ with vertices (1,2,3,5) and origin (1) : T = T1 + (T2 - T1).h + (T3 - T1).k + (T5 - T1).l dT/dh dT/dk dT/dl ; where T := x , y , z . Consequently, a brick is just a composite of 8 linear tetrahedra: (1,2,3,5) , (2,4,1,6) , (3,1,4,7) , (4,3,2,8) (5,7,6,1) , (6,5,8,2) , (7,5,8,3) , (8,6,7,4) with their origins at the vertices (nodes) of the brick. - * 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