such that the approximate solution of the boundary value problems, dealt with in

the previous chapter, generated by a ¬nite element analysis, converges to the exact

solution at mesh re¬nement. The answer is:

(i) The shape functions should be smooth within each element e , i.e. shape functions

are not allowed to be discontinuous within an element.

(ii) The shape functions should be continuous across each element boundary. This con-

dition does not always have to be satis¬ed, but this is beyond the scope of the present

book.

(iii) The shape functions should be complete, i.e. at element level the shape functions

should enable the representation of uniform gradients of the ¬eld variable(s) to be

approximated.

Conditions (i) and (ii) allow that the gradients of the shape functions show ¬nite

jumps across the element interface. However, smoothness in the element inte-

rior assures that all integrals in which gradients of the unknown function, say u,

occur can be evaluated. In Fig. 17.1(a) an example is given of an admissible shape

function. In this case the derivative of the shape function is discontinuous over

the element boundary, however the jump is ¬nite. In Fig. 17.1(b) the discontinu-

ous shape function at the element boundary leads to an in¬nite derivative and the

integrals in the weighted residual equations can no longer be evaluated.

Completeness When the ¬nite element mesh is re¬ned further and further, at

the element level the exact solution becomes more and more linear in the coordi-

nates and its derivatives approach constant values in each element. To ensure that

Shape functions and numerical integration

296

N

N

element 1 element 2 element 1 element 2

x

x

dN

dN

dx

dx

element 1 element 2

element 1 element 2

x

x

’∞

(a) Admissible shape functions (b) Non admissible shape functions

Figure 17.1

Continuity of shape functions.

an adequate approximation can be achieved, the shape functions have to contain

all constant and linear functions. Assuming that the exact solution is described by

an arbitrary linear polynomial, the element interpolation has to be able to exactly

describe this ¬eld. In mathematical terms this means the following. Assume u to

be approximated by

n

uh ( x) = Ni ( x) ui , (17.1)

i=1

with Ni ( x) the interpolation functions and ui the nodal values of uh ( x). Con-

sider the case, where the nodal values ui are selected to be related to the nodal

coordinates xi , yi , zi by

ui = c0 + c1 xi + c2 yi + c3 zi , (17.2)

according to a linear ¬eld. Substitution of Eq. (17.2) into (17.1) reveals:

n n n n

u h = c0 Ni ( x) +c1 Ni ( x) xi + c2 Ni ( x) yi + c3 Ni ( x) zi . (17.3)

i=1 i=1 i=1 i=1

297 17.2 Isoparametric, bilinear quadrilateral element

For the element types described in the current chapter the coordinates within the

elements are interpolated in the following way:

n

x= Ni ( x) xi

i=1

n

y= Ni ( x) yi

i=1

n

z= Ni ( x) zi . (17.4)

i=1

Completeness implies that Eq. (17.3) has to lead to

uh ( x, y, z) = c0 + c1 x + c2 y + c3 z, (17.5)

for arbitrary c0 , . . . , c3 . Hence

n

Ni ( x) = 1. (17.6)

i=1

The last equation can be used to check whether the shape functions may have been

speci¬ed in an adequate way.

17.2 Isoparametric, bilinear quadrilateral element

Consider a four-noded, straight-edged, element e ∈ I 2 as depicted in

R

Fig. 17.2(a). The nodal points are assumed to be numbered counterclockwise from

1 to 4. This element can be mapped onto a square, with local coordinates ξ and

·

(“1, 1) (1, 1)

3

4

3

4

ξ

y

1

2

2

1

x (“1, “1) (1, “1)

(a) Quadrilateral element with (b) Quadrilateral element with

respect to global coordinates respect to local coordinates

Figure 17.2

Quadrilateral element with respect to a global and a local coordinate system.

Shape functions and numerical integration

298

·, with ξ , · ∈[ ’1, 1]. A point with coordinates x, y ∈ is related to a point

e

ξ , · ∈[ ’1, 1] by the mapping

4

x(ξ , ·) = Ni (ξ , ·) xi (17.7)

i=1

4

y(ξ , ·) = Ni (ξ , ·) yi . (17.8)

i=1

The shape functions Ni may be determined by assuming the ˜bi-linear™ expansion:

x(ξ , ·) = ±0 + ±1 ξ + ±2 · + ±3 ξ ·

y(ξ , ·) = β0 + β1 ξ + β2 · + β3 ξ ·. (17.9)

The parameters ±j and βj (j = 0, 1, 2, 3) must be calculated such that the relations

x(ξi , ·i ) = xi

y(ξi , ·i ) = yi , (17.10)

are satis¬ed, where ξi and ·i refer to the local coordinates of the nodes. Applying

the restriction (17.10) to Eq. (17.9) leads to

⎛ ⎞⎛ ⎞

⎞⎛

1 ’1 ’1 ±0

1

x1

⎜x ⎟ ⎜1 ⎟

1 ’1 ’1 ⎟ ⎜ ±1 ⎟

⎜ 2⎟ ⎜ ⎟⎜

⎟=⎜

⎜ ⎟.

⎟⎜ (17.11)

1 ⎠ ⎝ ±2 ⎠

⎝ x3 ⎠ ⎝ 1 1 1

1 ’1 1 ’1 ±3

x4

From this set of equations the coef¬cients ±j can be linearly expressed in the

nodal coordinates x1 , x2 , x3 and x4 . By reorganizing, the eventual expressions for

the shape functions become:

1

N1 (ξ , ·) = (1 ’ ξ ) (1 ’ ·)

4

1

N2 (ξ , ·) = (1 + ξ ) (1 ’ ·)

4

1

N3 (ξ , ·) = (1 + ξ ) (1 + ·)

4

1

N4 (ξ , ·) = (1 ’ ξ ) (1 + ·) . (17.12)

4

The element is called isoparametric as both the spatial coordinates as well as the

element interpolation function uh are interpolated with the same shape functions,

that is

299 17.3 Linear triangular element

n

uh (ξ , ·) = Ni (ξ , ·) ui . (17.13)

i=1

17.3 Linear triangular element

There are two ways to arrive at a triangular element. First, it is possible to coalesce

two nodes of the quadrilateral element, for instance node 3 and 4. This is done by

setting x4 = x3 (with xT = [ x3 , y3 ] and xT = [ x4 , y4 ]) and by de¬ning new shape

∼ ∼ ∼3 ∼4

ˆ i according to

functions N

4

x= Ni xi

∼ ∼

i=1

= N1 x1 + N2 x2 + ( N3 + N4 ) x3

∼ ∼ ∼

ˆ ˆ ˆ

N1 N2 N3

3

ˆ

= Ni xi . (17.14)

∼

i=1

Fig. 17.3 illustrates this operation.

The second method is based on using so-called triangle coordinates. A conve-

nient set of coordinates »1 , »2 , »3 for a triangle can be de¬ned by means of the

following equations:

x = »1 x1 + »2 x2 + »3 x3 (17.15)

y = »1 y1 + »2 y2 + »3 y3 (17.16)

1 = »1 + » 2 + » 3 . (17.17)

To every set »1 , »2 , »3 ∈[ 0, 1] corresponds a unique set of Cartesian coordinates

x, y. The triangle coordinates are not independent, but related by Eq. (17.17).

3&4

3

4

y y

1

1

2

2

x x

Figure 17.3

Degeneration from quadrilateral element to triangular element.

Shape functions and numerical integration

300

Solving Eqs. (17.15) to (17.17) for »i leads to

1

»1 = (( x2 y3 ’ x3 y2 ) + ( y2 ’ y3 ) x + ( x3 ’ x2 ) y)

1

»2 = (( x3 y1 ’ x1 y3 ) + ( y3 ’ y1 ) x + ( x1 ’ x3 ) y)

1

»3 = (( x1 y2 ’ x2 y1 ) + ( y1 ’ y2 ) x + ( x2 ’ x1 ) y) , (17.18)

where

= ( x3 ’ x2 ) ( y1 ’ y2 ) ’ ( y2 ’ y3 ) ( x2 ’ x1 ) . (17.19)

Note that | | is 2 times the square of the triangle surface. Eqs. (17.15) to (17.17)

shows that »i can be (linearly) expressed in x and y: »i = »i ( x, y). The relation

satis¬es

»i ( xj , yj ) = δij . (17.20)

Now, the shape functions for the 3-node, 6-node and 7-node triangular elements