N T ( ξi ) = n(i,:).

∼

Using these shape functions the coordinate x( ξi ) within element e can be computed at

the i-th integration point:

x( ξi ) = N T ( ξi ) ∼e = n(i,:)*nodcoord.

x

∼

Similarly, the value of the solution at the i-th integration point of element e is

obtained via

u( ξi ) = N T ( ξi ) ue = n(i,:)*nodu.

∼

∼

In a similar fashion the shape function derivatives with respect to the local coordinate

ξ are stored in an array, called dndxi.

Structure of the Finite Element code Typically the structure of a ¬nite element

programme is as follows.

(i) Pre-processing: mesh generation, boundary condition speci¬cation and parame-

ter declaration. This should provide the topology array top, the coordinate array

coord and a number of auxiliary arrays containing boundary conditions and

material parameters.

(ii) Based on the mesh and element types used, the index array pos can be computed.

(iii) Assembly of the coef¬cient matrix q = K and the element array rhs = f . Let

∼

qe = K e and rhse = f e , then the assembly process in a MATLAB environment

∼

would look like

% nelem: the number of elements

for ielem = 1:nelem

% compute qe and rhse

[qe,rhse]=<elementfunction>(ielem,coord,top,.....)

% get the location of the degrees of freedom

% in the solution array

ii = nonzeros(pos(ielem,:));

% add the element coefficient matrix

% and the element right-hand side array

253 14.9 Example

% to the total coefficient matrix q

% and the load array rhs

q(ii,ii) = q(ii,ii) + qe;

rhs(ii) = rhs(ii) + rhse;

end

(iv) Solution of the set of equations taking into account the boundary conditions.

(v) Post-processing based on the solution, for instance by computing associated quanti-

ties such as heat ¬‚uxes or stresses.

14.9 Example

As an example consider the diffusion problem with the following parameter set-

ting. We consider the domain : 0 ¤ x ¤ 1, with prescribed essential boundary

conditions at x = 0 and x = 1. These conditions are: u(0) = 0 and u(1) = 0. There

are no natural boundary conditions. The material constant satis¬es: c = 1 and the

source term: f = 1.

The domain is divided into ¬ve elements of equal length. Fig. 14.9 shows

the solution. The left part displays the computed solution u (solid line) as well as

the exact solution (dashed line). Remarkably, in this one-dimensional case with

0.14 0.5

0.4

0.12

0.3

0.1 0.2

0.1

0.08

u 0

p

0.06

“0.1

“0.2

0.04

“0.3

0.02

“0.4

0 “0.5

0.2 0.4 0.6 0.8 1

0 0 0.2 0.4 0.6 0.8 1

x x

Figure 14.9

Five element solution. Left: (solid line) approximate solution uh ( x), (dashed line) exact solution u( x).

Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).

Numerical solution of one-dimensional diffusion equation

254

the current choice of parameters, the nodal solutions are exact. The right part of

the ¬gure shows the computed ¬‚ux, say ¬‚ux p = c du/dx. Again, the solid line

denotes the computed ¬‚ux p, which is clearly discontinuous from one element

to the next, and the dashed line denotes the exact solution. The discontinuity of

the computed ¬‚ux ¬eld is obvious: the ¬eld u is piecewise linear, therefore the

derivative du/dx is piecewise constant. The ¬‚ux p does not necessarily have to be

piecewise constant: if the parameter c is a function of x, the ¬‚ux p will be varying

within an element.

Mesh re¬nement leads to an improved approximate solution. For instance,

using ten rather than ¬ve elements, yields the results depicted in Fig. 14.10. The

impact of a varying c, say c = 1 + x is depicted in Fig. 14.11.

Changing the interpolation order of the shape functions Ni from linear to

quadratic, also has a signi¬cant impact on the results, in particular on the qual-

ity of the ¬‚ux prediction. For the constant c case, the solution becomes exact.

Also for c = 1 + x a signi¬cant improvement can be observed, as depicted in

Fig. 14.13.

0.5

0.14

0.4

0.12

0.3

0.1 0.2

0.1

0.08

p 0

u

0.06

“0.1

“0.2

0.04

“0.3

0.02

“0.4

“0.5

0

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

x x

Figure 14.10

Ten element solution. Left: (solid line) approximate solution uh ( x), (dashed line) exact solution u( x).

Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).

255 14.9 Example

0.6

0.09

0.08

0.4

0.07

0.2

0.06

0

0.05

u p

0.04 “0.2

0.03

“0.4

0.02

“0.6

0.01

0 “0.8

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

x x

Figure 14.11

Five element solution for c = 1 + x. Left: (solid line) approximate solution uh ( x), (dashed line) exact

solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux p( x).

0.5

0.14

0.4

0.12

0.3

0.1 0.2

0.1

0.08

0

p

u

0.06

“0.1

“0.2

0.04

“0.3

0.02

“0.4

“0.5

0

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

x x

Figure 14.12

Solution for c = 1 using ¬ve quadratic elements. Left: (solid line) approximate solution uh ( x),

(dashed line) exact solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux

p( x).

Numerical solution of one-dimensional diffusion equation

256

0.6

0.09

0.08

0.4

0.07

0.2

0.06

0

0.05

p

u

0.04 “0.2

0.03

“0.4

0.02

“0.6

0.01

“0.8

0

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1

x x

Figure 14.13

Solution for c = 1 + x using ¬ve quadratic elements. Left: (solid line) approximate solution uh ( x),

(dashed line) exact solution u( x). Right: (solid line) approximate ¬‚ux ph ( x), (dashed line) exact ¬‚ux

p( x).

Exercises

14.1 The method of weighted residuals can be used to ¬nd approximations of a

given function. Let f ( x) be a function that one would like to approximate

with a polynomial of the order n in a certain domain, say 0 ¤ x ¤ 1. Let

the polynomial be given by

h( x) = a0 + a1 x + · · · + an xn .

Ideally

g( x) = h( x) ’f ( x) = 0 for all x with 0 ¤ x ¤ 1.

In terms of the weighted residuals equation this may also be expressed as

1