x

x1 x2

Figure 14.5

Shape functions with respect to the global x-coordinate.

N1 N2

1

ξ

1

“1

Figure 14.6

Shape functions with respect to the local ξ -coordinate.

Numerical solution of one-dimensional diffusion equation

248

These integrals may be transformed into integrals using the local coordinate

system, according to

1 dx

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

dξ

’1

e

This requires the computation of the derivative dx/dξ . For this purpose the concept

of isoparametric elements is introduced. The shape functions Ni introduced for

the interpolation of the unknown function uh are also used for the relation between

the coordinates x and the coordinates ξ within an element:

n

x( ξ ) = Ni ( ξ ) xi = N T xe ,

e

(14.70)

∼∼

i=1

e

where xi are the coordinates of the nodes of the element. As a result the derivative

dx/dξ is obtained easily:

dN T

dx

= ∼ xe . (14.71)

dξ ∼

dξ

The element coef¬cient matrix requires the derivatives of the shape functions with

respect to the coordinate x. For this purpose

dNi dNi dξ

= , (14.72)

dx dξ dx

where

’1

dξ dx

= , (14.73)

dx dξ

is used, which is easily obtained from Eq. (14.71).

The integral on the right-hand side of Eq. (14.69) can be approximated by

means of numerical integration. The numerical integration of an arbitrary function

g( ξ ) over the domain ’1 ¤ ξ ¤ 1 is approximated by

nint

1

g( ξ ) dξ = g( ξi ) Wi , (14.74)

’1 i=1

where ξi denotes the location of the i-th integration point, nint is the total number

of integration points and Wi a weighting factor, i.e. the length of the ξ -domain,

associated with this integration point.

A simple example of a numerical integration is the trapezoidal integration

scheme, as illustrated in Fig. 14.7 (a). Integration of g( ξ ) over the domain

’1 ¤ ξ ¤ 1 using the trapezoidal integration rule yields:

1

g( ξ ) dξ ≈ g( ξ = ’1) + g( ξ = 1) , (14.75)

’1

249 14.7 Isoparametric elements and numerical integration

Table 14.1 Gaussian quadrature points up to nint = 3.

ξi

nint Wi

ξ1 = 0 W1 = 2

1

’1 1

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

2

3 3

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

3 5 5 9 9

g (ξ)

g (ξ)

ξ ξ

“1 +1 +1

“1 1 1

“

√3 √3

(a) (b)

Figure 14.7

(a) Trapezoidal integration (b) 2-point Gauss integration.

which corresponds to the shaded area in Fig. 14.7(a). For trapezoidal integration

the integration point positions ξi are given by

ξ1 = ’1, ξ2 = 1, (14.76)

while the associated weighting factors are

W1 = 1, W2 = 1. (14.77)

The trapezoidal integration rule integrates a linear function exactly. A 2-point

Gaussian integration rule, as depicted in Fig. 14.7(b) may yield a more accurate

result since this integration rule integrates up to a third order function exactly

using two integration points only. In this case the integral is approximated by

’1

1 1

g( ξ ) dξ ≈ g ξ = √ +g ξ=√ . (14.78)

3 3

’1

The location of the Gaussian integration (quadrature) points and the associated

weighting factors are summarized in Table 14.1.

Application to element coef¬cient matrix Use of the local coordinate sys-

tem, with isoparametric formulation and numerical integration to the element

coef¬cient matrix yields

Numerical solution of one-dimensional diffusion equation

250

dN dN T

Ke = ∼

c ∼ dx

dx dx

e

’1

1 dN dξ dN T dξ dξ

= ∼

c∼ dξ

dξ dx dξ dx dx

ξ =’1

nint

dN dN T dξ

≈ ∼

c∼ Wi . (14.79)

dξ dξ dx ξ =ξi

i=1

14.8 Basic structure of a finite element program

The objective of a ¬nite element program is, to compute the coef¬cient matrix

K and the right-hand side array f and eventually to solve the resulting system of

∼

equations taking the boundary conditions into account. To illustrate the typical

data structure and the layout of a ¬nite element program, consider, as an example,

the mesh depicted in Fig. 14.8.

The MATLAB programming language is used for explanation purposes. The

following data input is needed:

• Element topology First of all the domain is divided into a number of elements and

each node is given a unique global number. In this example two elements have been

used, the ¬rst element 1 is a quadratic element connecting nodes 3, 4 and 2 (in that

order) and the second element is a linear element having nodes 1 and 3 (again, in that

order). The node numbers of each element are stored in the topology array top, such

that the i-th row of this array corresponds to the i-th element. In the current example

the topology array would be:

3 4 2

top = .

1 3 0

Besides the node numbers of the element, a number of identi¬ers may be included

for each element, for instance to refer to different material parameters c or different

1 3 4 2

©1

©2

x1 x3 x4 x2

Figure 14.8

Mesh for a one-dimensional problem, consisting of a linear and a quadratic element.

251 14.8 Basic structure of a finite element program

element types, e.g. linear versus quadratic elements. In fact, the MATLAB code pro-

vided to experiment with, has two identi¬ers per element. Please consult the manual of

the code: mlfem_nac.

• Nodal coordinates The nodal coordinates ∼ are stored in the array coord, hence in

x

this example:

⎡ ¤

x1

⎢x ⎥

⎢ 2⎥

x = coord = ⎢ ⎥.

⎣ x3 ¦

∼

x4

The nodal coordinates of the nodes associated with the element e can be retrieved

from coord using top. The node numbers of element e can be extracted from the

array top by:

ii = nonzeros(top(e,:)),

such that the nodal coordinates of the e-th element are obtained via

x = nodcoord = coord(ii,:).

∼e

• Solution array The nodal unknowns, also called degrees of freedom, u are stored in

∼

the array sol:

⎡ ¤

u1

⎢u ⎥

⎢2 ⎥

u = sol = ⎢ ⎥.

⎣ u3 ¦

∼

u4

It is not necessary to store the degrees of freedom in a sequential manner, in fact, any

ordering may be chosen, as long as each array component corresponds to a unique

degree of freedom. To extract the nodal degrees of freedom for element e, ue from

∼

sol, a separate index array is needed: pos. The e-th row of the array pos contains

the location of the nodal degrees of freedom of element e in the array sol. For this

example

⎡ ¤

u3

⎢ ⎥ u1

u 1 = ⎣ u4 ¦ , u2 = ,

∼ ∼

u3

u2

hence the index array pos should contain

3 4 2

pos = .

1 3 0

Using this array, the nodal degrees of freedom of element e can be extracted from

sol via

ii = nonzeros(pos(e,:))

ue = nodu = sol(ii)

∼

Numerical solution of one-dimensional diffusion equation

252

• Shape functions The shape functions and their derivatives are needed at the integration

points N ( ξi ) and dN /dξ . The shape function values are stored in the array n such that

∼ ∼

at the i-th integration point