N8 = (1 ’ ξ ) (1 ’ ·2 ) . (17.35)

2

Other examples may be found in Zienkiewicz [18] and Hughes [10].

17.5 Numerical integration

’ I be some function, and assume that the integral:

Let f : R

e

f ( x) dx, (17.36)

e

over the domain e of an element is to be computed. In ¬nite element computa-

tions there is a mapping from the x-space to the ξ -space, such that (see Section

16.7, on isoparametric elements)

1 dx

f ( x) dx = f ( ξ) ( ξ ) dξ . (17.37)

dξ

’1

e

φ(ξ )

This integral can be approximated with a numerical integration rule:

nint

1

g( ξ ) dξ ≈ g( ξi ) Wi , (17.38)

’1 i=1

where ξi denotes the location of an integration point and Wi the associated weight

factor.

In Fig. 17.9 an interpretation is given of the above numerical integration rule.

At a discrete number of points ξi within the interval ξ ∈[ ’1, +1] the function

value g( ξi ) is evaluated. Related to each point ξi a rectangle is de¬ned with height

g( ξi ) and width Wi . Note that it is not necessary that the point ξi is located on the

symmetry line of the rectangle. By adding up the surfaces g( ξi ) Wi of all rectangles

an approximation is obtained of the total surface underneath the function, which

is the integral. It is clear that the weight factor Wi in Eq. (17.38) can be interpreted

as the width of the interval around ξi .

The integration rule that is mostly used is the Gaussian quadrature. In that

case the locations of the integration points and weight factors are chosen so as

Shape functions and numerical integration

306

Table 17.1 Gaussian quadrature up to nint =3.

ξi

nint Wi

ξ1 = 0 W1 = 2

1

’1 1

ξ1 = √ , ξ2 = √ W1 = W2 = 1

2

3 3

5 8

3 3

ξ1 = ’ 5 , ξ2 = 0, ξ3 = W1 = W3 = 9 , W2 = 9

3 5

g(ξ)

Wi

ξi

“1 1

ξ

Figure 17.9

Numerical integration of a function φ( ξ ).

to obtain optimal accuracy for polynomial expressions of g( ξ ). In Table 17.1 the

location of the integration points and the associated weight factors are given up to

nint = 3. For two-dimensional problems the above generalizes to

1 1

= f ( x( ξ , ·) , y( ξ , ·) ) j( ξ , ·) dξ d·

f ( x, y) d

’1 ’1

e

1 1

= g( ξ , ·) dξ d·, (17.39)

’1 ’1

with j( ξ , ·) according to Eq. (16.50) and

nint nint

1 1

g( ξ , ·) dξ d· ≈ g( ξi , ·j ) Wi Wj . (17.40)

’1 ’1 i=1 j=1

The above integration scheme can be elaborated for the 9-node rectangular

Lagrangian element in Fig. 17.10 using Table 17.2.

As was already remarked in Section 17.3 integration over a triangular domain is

not trivial. For a triangular element that is formed by degeneration from a quadri-

lateral element the integration can be performed in the same way, with the same

integration points, as given above.

In the case that triangular coordinates are used, the evaluation of the integrals is

far from trivial. If the triangle coordinates »1 and »2 are maintained by eliminating

307 17.5 Numerical integration

Table 17.2 Integration points in the square and associated weight

factors.

nint point location of the integration points Wi

ξ ·

’0.7745966692 ’0.7745966692

9 1 0.3086420047

’0.7745966692

2 0.7745966692 0.3086420047

3 0.7745966692 0.7745966692 0.3086420047

’0.7745966692

4 0.7745966692 0.3086420047

’0.7745966692

5 0 0.4938271818

6 0.7745966692 0 0.4938271818

7 0 0.7745966692 0.4938271818

’0.7745966692

8 0 0.4938271818

9 0 0 0.7901234686

4 7 3

·

ξ

8 9 6

1 5 2

Figure 17.10

Position of the integration points in the square.

2

1

»2

3 1

0

0 1

»1

Figure 17.11

Mapping of triangle in a »1 , »2 -coordinate system.

Shape functions and numerical integration

308

Table 17.3 Integration points in the triangles.

nint point location of the integration points weight factors

»1 »2 »3 Wi

3 1 0.5 0.5 0 0.16667

2 0 0.5 0.5 0.16667

3 0.5 0 0.5 0.16667

7 1 0.3333333333 0.3333333333 0.3333333333 0.11250

2 0.0597158717 0.4701420641 0.4701420641 0.00662

3 0.4701420641 0.0597158717 0.4701420641 0.00662

4 0.4701420641 0.4701420641 0.0597158717 0.00662

5 0.7974269853 0.1012865073 0.1012865073 0.00630

6 0.1012865073 0.7974269853 0.1012865073 0.00630

7 0.1012865073 0.1012865073 0.7974269853 0.00630

»2

»2

6

2

2 1 4

1

7

»1 »1

3 5

3

Figure 17.12

Position of the integration points in the triangles.

the coordinate »3 these can represented in a rectangular coordinate system as given

in Fig. 17.11. It can easily be seen that the surface integral of a function φ( »1 , »2 )

can be written as

1 1’»1

φ( »1 , »2 ) d»2 d»1 . (17.41)

0 0

Despite a problem with a variable integral limit, it appeared possible to derive

numerical integration rules for triangles. In Table 17.3 the position of integration

points in triangular coordinates as well as weight factors are given for two higher-

order elements.

309 Exercises

Exercises

17.1 A MATLAB script to calculate the shape functions and the derivatives with

respect to the isoparametric coordinates of 4-noded isoparametric bi-linear

elements in some point of the elements may consist of the following code:

% Shape = program to calculate shape functions

xi= ;

eta= ;

x=[ , , , ]™;

y=[ , , , ]™;

N=[0.25*(1-xi)*(1-eta);

0.25*(1+xi)*(1-eta);

0.25*(1+xi)*(1+eta);

0.25*(1-xi)*(1+eta)]

dNdxi=[-0.25*(1-eta);

0.25*(1-eta);

0.25*(1+eta);

-0.25*(1+eta)]

dNdeta=[-0.25*(1-xi);

-0.25*(1+xi);

0.25*(1+xi);

0.25*(1-xi)]

y

1 (10, 6)

(5, 6)

32

2 4

(0, 4)

(8, 5)

4

1

23

(6, 2)

1

(0, 0)

x

To study some of the properties of the shape functions the script will be

completed by considering the elements in the ¬gure above and extended

in the following. In this ¬gure the local node numbering has also been

indicated.

Shape functions and numerical integration

310

(a) Use the program to test whether

n