0 K11 K12 0

ˆ

K2 = ⎢ ⎥ (14.49)

⎣ ¦

2 2

0 K21 K22 0

0 0 0 0

⎡ ¤

0 0 0 0

⎢ ⎥

⎢ ⎥

0 0 0 0

ˆ

K3 = ⎢ ⎥, (14.50)

⎣ ¦

3 3

0 0 K11 K12

3 3

0 0 K21 K22

such that

ˆ ˆ ˆ

K = K1 + K2 + K3, (14.51)

leading to

⎡ ¤

1 1

K11 K12 0 0

⎢ ⎥

+ K11

1 1 2 2

⎢ ⎥

K21 K22 K12 0

K=⎢ ⎥. (14.52)

K22 + K11

⎣ ¦

3 3

2 2

0 K21 K12

3 3

0 0 K21 K22

ˆ

In computer codes, however, the individual matrices K i are never formed explic-

ˆ

itly. The non-trivial components of K i are supplied directly to the appropriate

position within K. This process of assembling the global matrix K based on the

contributions of the individual element matrices K e is called the assembly process.

For the right-hand side the same procedure is followed. For element 2 it can be

written analogously:

2 2

f1 f1

wT f 2 = =

w2 w2 w2 w3

∼2 2 2

1 2 f2 f2

∼

⎡ ¤

0

⎢ f2 ⎥

⎢1 ⎥ ˆ

= ⎥ = wT ∼2 .

⎢2 f (14.53)

w1 w2 w3 w4

⎣ f2 ¦ ∼

0

fˆ

∼2

245 14.5 Galerkin approximation

Following the same procedure for the other elements and adding up the contribu-

tion for the internal source term for all (three) elements gives

ˆ ˆ ˆ

f = f 1 + f 2 + f 3. (14.54)

∼ ∼ ∼ ∼

This leads to

⎡ ¤

1

f1

⎢ ⎥

f2 + f1

1 2

⎢ ⎥

w f int =

T

⎢ ⎥. (14.55)

w1 w2 w3 w4

f2 + f1

⎣ ¦

3

2

∼ ∼

3

f2

What remains is the term B in Eq. (14.25). The effect of this boundary term B may

be included via

B = ’w( a) pa + w( b) pb = wT f ext , (14.56)

∼ ∼

where f ext contains pa and pb at the appropriate positions, according to

∼

⎡ ¤

’pa

⎢ ⎥

⎢0⎥

B = w1 w2 w3 w4 ⎢ ⎥. (14.57)

⎣0¦

pb

This ¬nally leads to an equation of the form:

wT Ku = wT ( f int + f ext ) . (14.58)

∼

∼ ∼ ∼ ∼

Using the fact that Eq. (14.58) must hold for ˜all™ w, this results in the so-called

∼

discrete set of equations:

¤⎡ ¤ ⎡ ¤

⎡

f1 ’ pa

1 1 1

K12 0 0 u1

K11

⎥⎢ ⎥⎢ ⎥

⎢ K1 K1 + K2 f2 + f1

2 1 2

⎥⎢ ⎥⎢ ⎥

⎢ 21 K12 0 u2

⎥=⎢

⎥⎢ ⎥ . (14.59)

⎢ 22 11

K22 + K11 K12 f2 + f1

¦⎣ ¦⎣ ¦

⎣0 3 3 3

2 2 2

K21 u3

f2 + pb

3 3 3

0 0 K21 K22 u4

Assume that in the example of Fig. 14.4 at x = x1 , an essential boundary condi-

tion is prescribed u( x1 ) = U. This means that pa = pu is unknown beforehand.

At x = x4 a natural boundary condition is prescribed, so pb = P is known

beforehand.

Then

¤⎡ ¤⎡ ¤

⎡

f1 ’ pu

1 1 1

K12 0 0 U

K11

0 ⎥ ⎢ u2 ⎥ ⎢ f2 + f1 ⎥

⎥⎢ ⎥ ⎢1

⎢ K1 K1 + K2 2 2

⎥

⎢ 21 K12

=⎢ 2

3 ⎥⎢ u ⎥ 3 ⎥ . (14.60)

⎢ 22 11

2 + K3

⎣ f2 + f1 ¦

11 K12 ¦ ⎣ 3 ¦

⎣0 2

K21 K22

f2 + P

3 3 3

0 0 K21 K22 u4

Numerical solution of one-dimensional diffusion equation

246

It is clear that in Eq. (14.60) the unknowns are u2 , u3 , u4 on the left-hand side

of the equation, and pu on the right-hand side of the equation. So both columns

can be divided in a known part and an unknown part, depending on the essential

and natural boundary conditions that have been prescribed. The next section will

outline how this equation is partitioned to facilitate the solution process.

14.6 Solution of the discrete set of equations

Let Eq. (14.60) be written as

Ku=f (14.61)

∼ ∼

where f = f int + f ext .

∼ ∼ ∼

The unknowns can be partitioned into two groups. First, some of the compo-

nents of the column u will be prescribed. This subset of u is labeled up . The

∼ ∼ ∼

remaining components of u are the actual unknowns, labeled uu . In a similar

∼ ∼

manner K and f can be partitioned. Consequently, Eq. (14.60) can be rewritten as

∼

f

K uu K up uu

=

∼ ∼u . (14.62)

fp

up

K pu K pp ∼ ∼

It is emphasized that the right-hand side partition f u associated with the unknowns

∼

uu will be known, and reversely, f p will be unknown as up is known. Since up is

∼ ∼ ∼

∼

known, the actual unknowns uu can be solved from

∼

K uu uu = f u ’ K up up . (14.63)

∼ ∼

∼

Notice that at the part of the boundary where u is prescribed, the associated

external load fp is unknown.

The components of f p can be obtained by simple multiplication, after having

∼

solved uu from Eq. (14.63):

∼

f = K pu uu + K pp up . (14.64)

∼ ∼

∼p

14.7 Isoparametric elements and numerical integration

In Section 14.4 the concept of shape functions has been introduced. Within each

element uh has been written as

uh | = N T ( x) ue . (14.65)

∼

∼

e

247 14.7 Isoparametric elements and numerical integration

The shape functions are simple polynomial expressions in terms of the coor-

dinate x. For instance for a linear interpolation the shape functions are linear

polynomials, according to

x ’ x1 x ’ x1

N1 = 1 ’ N2 =

, , (14.66)

x2 ’ x1 x2 ’ x1

where x1 and x2 denote the position of the nodes of the element. In this case the

shape functions are linear functions of the global coordinate x. It is appropriate

in the context of a generalization to more-dimensional problems to introduce a

local coordinate ’1 ¤ ξ ¤ 1 within each element such that ξ = ’1 and ξ =

1 correspond to the edges of the element. With respect to this local coordinate

system, the shape functions may be written as

1 1

N1 = ’ ( ξ ’ 1) , N2 = ( ξ + 1) . (14.67)

2 2

This is visualized in Fig. 14.5.

Computation of components of the element coef¬cient matrix and the element

load array requires the evaluation of integrals of the form

f ( x) dx . (14.68)

e