Next: Neuman and Cauchy Conditions Up: The Numerical Inclusion of Previous: Payne-Irons Method

# Calculation of Direction Cosines

The main constituent of both the Neumann and Cauchy boundary conditions is the rate of change of the scalar field in the direction of the outward normal to the boundary. This has been denoted
 (4.13)
In many problems the value of this will be known in some form; in fact in many cases it will be a constant. However there will be occasions where the explicit value of is not known but is given as the rates of change of with respect to some set of global axes.

For example consider a two-dimensional problem formulated in a Cartesian coordinate system . The boundary values of may well be given with respect to this system of axes. Thus on a boundary

 (4.14)
These are adequate provided the boundary is parallel to one of the coordinate directions.

In the example region shown in Figure 4.1, clearly along AD ()
 (4.15)
and along DE ()
 (4.16)
However, problems present themselves whenever the boundary is not parallel to the coordinate axes, for example along EF. The normal derivative must now be constructed from the known parts. Clearly this will be
 (4.17)
where and are the directions cosines of the outward normal at some point along EF, given by
 (4.18)
i and j being the unit vectors associated with the coordinate directions and . (4.15) and (4.16) can be easily derived from (4.17) and (4.18) and the whole analysis extended to three-dimensional problems.

A straight boundary not parallel with the coordinate axes does not present real difficulty since and will be constants. However in the most general case, where the boundary could be any smooth curve, and will be functions of position. Within the context of finite element analysis the direction cosines will be required on sides of elements.

In general the element in global coordinates is generated from a parent element by means of some isoparametric transformation. The most used transformation is

 (4.19) (4.20) (4.21)

where the are a set of shape functions (see Chapter 2) and are the coordinates of the element nodes in global coordinates. The are defined in terms of a coordinate system natural to the parent element (most often called local coordinates).

As an example, consider the triangular elements in Figure 4.2. Figure 4.2(a) shows the elements in the global coordinate system and Figure 4.2(b) shows the parent element. The boundary is approximated by a series of element sides; for each side the outward normal is required.

The construction of these direction cosines is based on two results (Rutherford, 1962). Firstly any element side can be presented as part of a set of equipotential surfaces in the coordinate system. These are characterised by

 (4.22)
where is some continuous function possessing continuous partial derivatives.

The second is a result from coordinate geometry:

If at any point P on the surface , then is perpendicular to the surface, passing through P,
where
 (4.23)
It is clear from this that the direction cosines of the outward normal are proportional to the components of :
 (4.24)
and just as this applies to the elements in the global system it will also apply to the parent element in the coordinate system.

In the local coordinate system the direction cosines of the outward normal to any of the sides are simple to calculate. As in the global case these are proportional to the components of , where is the equation of one of the element sides; hence

 (4.25)
The chain rule governing partial differentiation reads:
 (4.26)
and using (4.20) the direction cosines of the outward normal to the parent element can be transformed into the outward normal to an element in the global coordinate system.

Only the direction of the normal is given so finally the vector is normalised. The components of this normalised vector will be the direction cosines and . Within the context of a finite element program (Segment 3.2 is the best example) the code to calculate the direction cosines is contained within the boundary integration loop since the normal direction will be required at the Gauss points on the boundary:


CALL MATMUL(LDER, ILDER, JLDER, GEOM, IGEOM, JGEOM, JAC,
*     IJAC, JJAC, DIMEN, NODEL, DIMEN, ITEST)
CALL MATINV(JAC, IJAC, JJAC, JACIN, IJACIN, JJACIN, DIMEN,
*     DET, ITEST)
CALL DCSTRI(JACIN, IJACIN, JJACIN, SIDNUM, COSIN, ICOSIN,
*     ITEST)

The routine DCSTRI uses the current transformation Jacobian and the side number ( SIDNUM) on which the direction cosines are required to make the calculation. The direction cosines are return in COSIN. The code above forms the Jacobian ( JAC) using (2.25) and then inverts it using MATINV. The routines DCSQUA and DCSBRK perform similar operations for rectangular and brick elements.

Next: Neuman and Cauchy Conditions Up: The Numerical Inclusion of Previous: Payne-Irons Method
Chris Greenough (c.greenough@rl.ac.uk): September 2001