sci.math.num-analysis SUNA46, Linear tetrahedron ========================== Let's consider the simplest non-trivial finite element shape in 3 dimensions, the linear tetrahedron. Function behaviour is approximated inside a tetrahedron by a linear interpolation between the function values at the vertices, or nodal points. Let T be such a function, and x,y,z the coordinates, then: T = A.x + B.y + C.z + D Hence: T1 = A.x1 + B.y1 + C.z1 + D Substract: T - T1 = A.(x - x1) + B.(y - y1) + C.(z - z1) Then the constants A , B and C are determined by solving the equations: | T2 - T1 | | x2 - x1 y2 - y1 z2 - z1 | | A | | T3 - T1 | = | x3 - x1 y3 - y1 z3 - z1 | | B | | T4 - T1 | | x4 - x1 y4 - y1 z4 - z1 | | C | --------------------------- matrix J The determinant of the matrix called J represents 6 times the volume of the tetrahedron (1,2,3,4). Write the substraction as: T - T1 = h.(T2 - T1) + k.(T3 - T1) + l.(T4 - T1) Or: T - T1 = | T2 - T1 T3 - T1 T4 - T1 | | h | | k | | l | Then it may be a theorem that h , k , l too are a solution of the following 3 equations, an isoparametric transformation, namely: | x - x1 | | x2 - x1 x3 - x1 x4 - x1 | | h | | y - y1 | = | y2 - y1 y3 - y1 y4 - y1 | | k | | z - z1 | | z2 - z1 z3 - x1 z4 - z1 | | l | --------------------------- matrix J' where ' = transposed: see J above. It's consistent: T - T1 = | T2 - T1 T3 - T1 T4 - T1 | | h | | k | = | l | = | A B C | . J' * (J')^-1 . | x - x1 | = | A B C | | x - x1 | | y - y1 | | y - y1 | | z - z1 | | z - z1 | Q.E.D. Thus again the _same_ expression holds for the function T as well as for the coordinates x , y , z , which an isoparametric transformation is meant to be. h , k , l can be interpreted as volume coordinates (was: area coordinates). Instead of volume-coordinates, we prefer to talk about the _local coordinates_ h , k , l of an element, in contrast to the _global coordinates_ x , y . z . It is possible that local coordinates coincide with the global coordinates. A tetrahedron for which this is the case is called a _parent element_. Let's reconsider for a moment the expression: T - T1 = h.(T2 - T1) + k.(T3 - T1) + l.(T4 - T1) Partial differentiation to h , k , l gives: dT/dh = T2 - T1 dT/dk = T3 - T1 dT/l = T4 - T1 So: T = T1 + h.dT/dh + k.dT/dk + l.dT/dl 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 - l).T1 + h.T2 + k.T3 + l.T4 Here the functions (1 - h - k) , h , k , l 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 ; N4 = l So we have two representations, which are allmost trivially equivalent: T = T1 + h.(T2 - T1) + k.(T3 - T1) + l.(T4 - T1) : Finite Difference T = (1 - h - k).T1 + h.T2 + k.T3 + l.T4 : Finite Element Partial derivatives dT/d(x,y,z) will be considered. In order to take care of their discretization, transformation to isoparametric local coordinates (h,k,l) at a Finite Element, being a general F.E. technique, seems to be advantageous. Derivatives to the global coordinates must then be expressed in derivatives to the local coordinates. According to the chain rules, written in matrix form: | dT/dh | | dx/dh dy/dh dz/dh | | dT/dx | | dT/dk | = | dx/dk dy/dk dz/dk | | dT/dy | | dT/dl | | dx/dl dy/dl dz/dl | | dT/dz | According to Cramer's rule, the solution of the above equations is: | dT/dh dy/dh dz/dh | dT/dx = Det( | dT/dk dy/dk dz/dk | ) / J | dT/dl dy/dl dz/dl | | dx/dh dT/dh dz/dh | dT/dy = Det( | dx/dk dT/dk dz/dk | ) / J | dx/dl dT/dl dz/dl | | dx/dh dy/dh dT/dh | dT/dz = Det( | dx/dk dy/dk dT/dk | ) / J | dx/dl dy/dl dT/dl | | dx/dh dy/dh dz/dh | Where: J = Det( | dx/dk dy/dk dz/dk | ) is the Jacobian determinant. | dx/dl dy/dl dz/dl | Thus solving for dT/d(x,y,z) we find general formulas which are valid for any 3-D element. ^^^^^^^ It is remarked that the first three determinants can be derived from the last one by substituting x := T , y := T , z := T respectively. More easy: due to isoparametrics, only the nodal values of these quantities have to be replaced. This means that the function behaviour of partial derivatives, being a quotient of two alike determinants, can be studied by considering the Jacobian matrix or determinant alone. This will be heavily used in the sequel. For the tetrahedron (1,2,3,4) we trivially find: T = T1 + h.(T2 - T1) + k.(T3 - T1) + l.(T4 - T1) dT/dh dT/dk dT/dl Substitutution of T := x , y , z gives the determinant of the matrix: | x2 - x1 y2 - y1 z2 - z1 | | x3 - x1 y3 - y1 z3 - z1 | | x4 - x1 y4 - y1 z4 - z1 | - * 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