sci.math.num-analysis SUNA03, Convex Quadrilaterals ============================= Let's consider next the simplest-but-one element shape in 2-D: a quadrilateral. Function behaviour may be approximated inside such a 3 quadrilateral by a _bilinear_ interpolation between / \ the function values at the vertices, or nodal points. / \ Let T be such a function, and x,y coordinates, / \ then: 1 0,0 4 T = A + B.x + C.y + D.x.y <--there--\ \ / | \ / As an example, consider the quadrilateral depicted --/ \ / with nodes given by the 2'nd and 3'th row of the matrix 2 that originates when specifying T for the nodes: | T1 | | 1 -1 0 0 | | A | The last column of our matrix is zero. | T2 | = | 1 0 -1 0 | | B | Thus A,B,C and D cannot be found here! | T3 | | 1 0 +1 0 | | C | This means that a method as has been | T4 | | 1 +1 0 0 | | D | used for the linear triangle can not 1 x y x.y be employed for higher order elements. There is, however, another way out. Let's assume that the same expression holds for the function T as well as for the coordinates x and y : an ISOPARAMETRIC transformation again. The next step is to find a good _parent_ element. Well, the above one surely did'nt work out. Let's try the following: 3---------4 T = A + B.h + C.k + D.h.k gives: | | | | | T1 | | 1 0 0 0 | | A | | | | T2 | = | 1 1 0 0 | | B | F.E. <- F.D. | | | T3 | | 1 0 1 0 | | C | transition 1---------2 | T4 | | 1 1 1 1 | | D | h k Solving for A,B,C,D | A | | +1 0 0 0 | | T1 | isn't too | B | = | -1 +1 0 0 | | T2 | F.D. <- F.E. difficult: | C | | -1 0 +1 0 | | T3 | transition | D | | +1 -1 -1 +1 | | T4 | So we have two equivalent representations: T = T1 + (T2 - T1).h + (T3 - T1).k + (T1 - T2 - T3 + T4).h.k : F.D.-like T = (1 - h).(1 - k).T1 + h.(1 - k).T2 + (1 - h).k.T3 + h.k.T4 : F.E.-like Partial derivatives are: dT/dh = T2 - T1 dT/dk = T3 - T1 and: d^2(T)/dh.dk = T1 - T2 - T3 + T4 Shape functions: (1 - h).(1 - k) ; h.(1 - k) ; (1 - h).k ; h.k at nodal point: 1 2 3 4 Global coordinates are related to local coordinates inside a general curved quadrilateral by the isoparametric transformation: x = (1 - h).(1 - k).x1 + h.(1 - k).x2 + (1 - h).k.x3 + h.k.x4 y = (1 - h).(1 - k).y1 + h.(1 - k).y2 + (1 - h).k.y3 + h.k.y4 The differential form of an isoparametric transformation is often even more important (here d*/d* are partial derivatives): dx = dx/dh.dh + dx/dk.dk | dx | | dx/dh dx/dk | | dh | dy = dy/dh.dh + dy/dk.dk or: | | = | | | | | dy | | dy/dh dy/dk | | dh | Inverse differential transform: | dh | | dy/dk -dx/dk | | dx | | | = | | / J | | | dk | | -dy/dh dx/dh | | dx | The inverse transform exists if and only if the Jacobian determinant is never zero. Because J will be a function of (h,k) this also means that it should not switch its sign. We shall demand that it is positive: J = dx/dh.dy/dk - dx/dk.dy/dh > 0 Here: dx/dh = (1 - k).(x2 - x1) + k.(x4 - x3) dy/dh = (1 - k).(y2 - y1) + k.(y4 - y3) dx/dk = (1 - h).(x3 - x1) + h.(x4 - x2) dy/dk = (1 - h).(y3 - y1) + h.(y4 - y2) Substitution gives: J = (1 - h).(1 - k).J1 + h.(1 - k).J2 + (1 - h).k.J3 + h.k.J4 Where: J1 = (x2 - x1).(y3 - y1) - (x3 - x1).(y2 - y1) J2 = (x2 - x1).(y4 - y2) - (x4 - x2).(y2 - y1) J3 = (x4 - x3).(y3 - y1) - (x3 - x1).(y4 - y3) J4 = (x4 - x3).(y4 - y2) - (x4 - x2).(y4 - y3) Therefore the Jacobian is expressed by the quadrilateral shape functions into its values Jk at the nodal points, as it should be. Well, not really. Look at the Finite Difference representation: J = J1 + (J2 - J1).h + (J3 - J1).k + (J1 - J2 - J3 + J4).h.k The last term is zero. This is easiest seen by considering the fact that the area of the quadrilateral is equal to the sum of the areas of the triangles (1) and (4), but also to the sum of the areas of the triangles (2) and (3). Therefore J1 + J4 = J2 + J3 which proves the above claim. We may conclude | that J is not a bi-linear but only a _linear_ expression in the Jk 's. Now it was demanded that the Jacobian must be positive over the whole element. This can only be the case if all its values at the vertices are postive. Let's take vertex (1) as an example. Here: J1 = |_12_|.|_13_|.sin(phi) . _1x_ is the vector which points from 1 to (x), (phi) is the angle between the vectors _12_ and _13_. This means that sin(phi) > 0 or 0 < phi < Pi . Thus obtuse angles are forbidden. Quadrilateral elements must be convex. Disclaimer: the above is still well known. To be continued ... - * 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