. 50
( 67 .)


1 N1 N2

x1 x2

Figure 14.5
Shape functions with respect to the global x-coordinate.

N1 N2


Figure 14.6
Shape functions with respect to the local ξ -coordinate.
Numerical solution of one-dimensional diffusion equation

These integrals may be transformed into integrals using the local coordinate
system, according to
1 dx
f ( x) dx = f ( x( ξ ) ) dξ . (14.69)


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:
x( ξ ) = Ni ( ξ ) xi = N T xe ,
where xi are the coordinates of the nodes of the element. As a result the derivative
dx/dξ is obtained easily:
dN T
= ∼ xe . (14.71)
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
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
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:
g( ξ ) dξ ≈ g( ξ = ’1) + g( ξ = 1) , (14.75)
249 14.7 Isoparametric elements and numerical integration

Table 14.1 Gaussian quadrature points up to nint = 3.
nint Wi
ξ1 = 0 W1 = 2
’1 1
ξ1 = √ , ξ2 = √ W1 = W2 = 1
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
g( ξ ) dξ ≈ g ξ = √ +g ξ=√ . (14.78)
3 3
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

dN dN T
Ke = ∼
c ∼ dx
dx dx
1 dN dξ dN T dξ dξ
= ∼
c∼ dξ
dξ dx dξ dx dx
ξ =’1
dN dN T dξ
≈ ∼
c∼ Wi . (14.79)
dξ dξ dx ξ =ξi

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


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
this example:
⎡ ¤
⎢x ⎥
⎢ 2⎥
x = coord = ⎢ ⎥.
⎣ x3 ¦

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,:).

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

the array sol:
⎡ ¤
⎢u ⎥
⎢2 ⎥
u = sol = ⎢ ⎥.
⎣ u3 ¦

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
⎡ ¤
⎢ ⎥ u1
u 1 = ⎣ u4 ¦ , 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

• 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


. 50
( 67 .)