sci.math.numanalysis
SUNA03, Convex Quadrilaterals
=============================
Let's consider next the simplestbutone element shape in 2D: 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:
34 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
12  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 bilinear 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
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood