MS5019 – FEM 1

7.0. Introduction
z After considering the LST element in previous chapters, we can see that
the development of element matrices and equations expressed in terms
of global coordinate system becomes an enormously difficult task.
z Hence, the isoparametric formulation was developped. This formulation
allows elements to be created that are non-rectangular and have curved
sides.
z Numerous commercial computer FE software adapted this formulation
for their various libraries of elements.
z Here, we will illustrate the isoparametric formulation for the simple bar
element and the rectangular plane stress element.
z Also, we will learn a little bit about numerical integration method for
evaluating the element stiffness matrix.
MS5019 – FEM 2

1
 7.1. Isoparametrid Formulation of Bar Element
z The term isoparametric is derived from the use of the same shape (or
interpolation function) [N] to define the element’s geometric shape as
are used to define the displacement within the element.
z Thus, when the shape function is u = a1 + a2s for the displacement, we
use x = a1 + a2s for the description of the nodal coordinate of a point on
the bar element, and hence the physical shape of the element.
z Isoparametric element equations are formulated using a natural (or
intrinsic) coordinate system s that is defined by element geometry and
NOT by the element orientation in global coordinate system.
z There is a relationship (called a transformation mapping) between the
natural coordinate system s and the global coordinate system x for each
element of a specific structure, and this relationship must be used in the
element equations formulations.
MS5019 – FEM 3

z Step 1 Select Element Type

First, the natural coordinate s is attached to the element, with the origin
located at the center of the element, as shown in Fig. 7-1(a). The s axis
need not be parallel to the x axis – this is only for convenience.
We consider the bar element to have two d.o.f – axial displacement u1
and u2 at each node associated with the global x axis.

s=0
s
s = −1 s = +1 x1 x2
x, u
1 L 2 1 L 2

(a) (b)

Figure 7-1 Linear bar element in

(a) natural coordinate system s and in (b) global coordinate system x

MS5019 – FEM 4

2
 For the special case when the s and x are parallel to each other, the
s and x coordinates can be related by
L
x = xc +
s (7.1.1)
2
where xc is the global coordinate of the centroid.
The shape function used to define a position within the bar are
found in a manner similar to that used in Chapter 3 to define
displacement within a bar. We begin by relating the natural coordinate
to the global coordinate by.

x = a1 + a2 s (7.1.2)

where we note that s is such that - 1 ≤ s ≤ 1.

MS5019 – FEM 5

Solving for ai ' s in terms of x1 and x2 , we obtain

x = 12 [(1 − s ) x1 + (1 + s ) x2 ] (7.1.3)

or, in matrix form, we can express Eq. (7.1.3) as

⎧x ⎫
{x} = [ N1 N 2 ] ⎨ 1 ⎬ (7.1.4)
⎩ x2 ⎭
where the shape functions in Eq. (7.1.4) are
1− s 1+ s
N1 = N2 = (7.1.5)
2 2
The shape function Eq. (7.1.5) map the s coordinate of any point in the
element to the x coordinate (for s = 1 we obtain x = x1). These shape
function are shown in Figure 7-2. N1 represents the physical shape of the
coordinate x when plotted over the length of the element for x1=1 and
x2=0, (and N2 for x1=0 and x2=1). Again, we must have N2 + N2 = 1.
MS5019 – FEM 6

3
 N1 N2
1 1

s s
-1 0 1 -1 0 1
Figure 7-2 Shape function variations with natural coordinates

z Step 2 Select Displacement Functions

The displacement function within the bar is now defined by the same
shape functions, Eqs. (7.1.5), as used to define the element shape; that
is
⎧u ⎫
u = [ N1 N2 ] ⎨ 1 ⎬ (7.1.6)
⎩u 2 ⎭

MS5019 – FEM 7

When a particular coordinate s of the point of interest is substituted into

[N], Eq. (7.1.6) yields the displacement of a point on the bar in terms of
the nodal d.o.f u1 and u2. Since u and x are defined by the same shape
functions at the same nodes, the element is called isoparametric.

z Step 3 Define the Strain/Displacement and Stress/Strain

Relationship
We now want to formulate element matrix [B] to evaluate [k]. We use
the isoparametric formulation to illustrate its manipulations. For a
simple bar element, no real advantage is evident; but for higher-order
elements, the advantage will become more clear because relatively
simple computer program formulation will result.
The displacement u is now a function of s as given by Eq. (7.1.6).
Therefore, we must apply the chain rule of differentiation as follows
MS5019 – FEM 8

4
 du du dx
= (7.1.7)
ds dx ds
We can evaluate (du ds ) and (dx ds) using Eqs. (7.1.6) and (7.1.3).
We seek (du dx) = ε x . Therefore, we solve Eq. (7.1.7) for (du dx) as
⎛ du ⎞
⎜ ⎟
du ⎝ ds ⎠
= (7.1.8)
dx ⎛ dx ⎞
⎜ ⎟
⎝ ds ⎠
Using Eqs. (7.1.6) and (7.1.3) in Eq. (7.1.8), we obtain
⎧u ⎫
{ε x } = du = ⎡⎢− 1 1 ⎤⎥ ⎨ 1 ⎬ (7.1.9)
dx ⎣ L L ⎦ ⎩u2 ⎭
where we have used
dx x2 − x1 L
= = (7.1.10)
ds 2 2
in obtaining Eq. (7.1.9).
MS5019 – FEM 9

Since {ε x } = [ B ]{d }, the strain/dis placement matrix [ B ] is then given in

Eq. (7.1.9) as
⎡ 1 1⎤
[ B ] = ⎢− (7.1.11)
⎣ L L ⎥⎦
For higher - order elements, such as quadratic bar with t hree nodes,
[ B ] becomes a function of natural coordinate s.

z Step 4 Derive Element Stiffness Matrix and Equations

The stiffness matrix is
L
[k ] = ∫ [ B ]T [ D ][ B ] A dx (7.1.12)
0

However, in general, we must transform the coordinate x to s because [ B ] is,

in general, a function of s. This general type of transformation is given by

MS5019 – FEM 10

5
 L 1

∫
0
f ( x ) dx = ∫ f ( s) J ds
−1
(7.1.13)

where J is called the Jacobian. In the one - dimensiona l case, we have J = J.

For the simple bar element, from Eq. (7.1.3), we have
dx L
J = = (7.1.14)
ds 2
because x2 − x1 = L. Observe that in Eq. (7.1.14), the Jacobian relates an
element length in the global coordinate system to an element length in the
natural coordinate s system. In general, J is a function s and depends on the
numerical values of the nodal coordinate s.
Using Eqs (7.1.13) and (7.1.14) in Eq. (7.1.12), we obtain the expression of
stiffness matrix in natural coordinate s
1
L
[k ] = ∫
2 −1
[ B ]T E [ B ] A ds (7.1.15)

MS5019 – FEM 11

where, for the 1 - D case, we have used the modulus of elasticity E = [D]
in Eq. (7.1.15). Substituting Eq. (7.1.11) in Eq. (7.1.15) and performing
the simple integration, we obtain
AE ⎡ 1 − 1⎤
[k ] = (7.1.16)
L ⎢⎣− 1 1 ⎥⎦
which the same as Eq. (3.1.14). For higher-order 1-D elements, the
integration in closed form becomes difficult if not impossible. Even the
simple rectangular elements stiffness matrix is difficult to evaluate in
closed form. In this case, the use of numerical integration, illustrates the
distinct advantage of the isoparametric formulation of the equations.

MS5019 – FEM 12

6
7.2. Rectangular Plane Stress Element
z We will now develop the rectangular plane stress element stiffness
matrix. We will later refer to this element in the isoparametric
formulation of a general quadrilateral element.
z Two advantages of the rectangular element over the triangular element
are:
1. ease of data input
2. simpler interpretation of output stresses

z Step 1 Select Element Type

Consider the rectangular element shown in Figure 7-3 (all interior
angles are 90o) with corner nodes 1, 2, 3, and 4 (again labeled CCW)
and base and height dimensions are 2b and 2h, respectively.

MS5019 – FEM 13

v4 y, v v3
u4 b b
u3
4 3
h
x, u
h
u1 1 2
u2
v1 v2

Figure 7-3 Basic four-node rectangular element

The unknown nodal displacements are now given by

⎧{d1}⎫
⎪{d }⎪
{d } = ⎪⎨ 2 ⎪⎬ = [u1 v1 u2 v2 u3 v3 u4 v4 ]T (7.2.1)
⎪{d 3 }⎪
⎪⎩{d 4 }⎪⎭

MS5019 – FEM 14

7
z Step 2 Select Displacement Functions
For a compatible displacement field, the element displacement functions u
and v must be linear along each edge because only two points (the corner
nodes) exist along each edge. We then select the linear displacements
functions as
u ( x, y ) = a1 + a2 x + a3 y + a4 xy
( 7 . 2 .2 )
v ( x, y ) = a5 + a6 x + a7 y + a8 xy
We can proceed in the usual manner to eliminate the ai ' s from Eqs. (7.2.2)
to obtain
1
u ( x, y ) = [(b − x )( h − y )u1 + (b + x )( h − y )u 2
4bh
+ (b + x )( h + y )u3 + (b − x )( h + y )u 4 ]
(7.2.3)
1
v ( x, y ) = [(b − x )( h − y )v1 + (b + x )( h − y )v2
4bh
+ (b + x )( h + y )v3 + (b − x )( h + y )v4 ]
MS5019 – FEM 15

The displacement expressions in Eq. (7.2.3) can be expressed equivalently

in term of the shape functions and unknown nodal displacemets as
{ψ } = [ N ]{d } ( 7 .2 . 4 )

where the shape functions are given by

(b − x)(h − y ) (b + x)(h − y )
N1 = N2 =
4bh 4bh
(7.2.5)
(b + x)(h + y ) (b − x)(h + y )
N3 = N4 =
4bh 4bh
and the N i ' s are again such that N1 = 1 at node 1 and N 2 = 0 at all the other
nodes, with similar requirements for the other shape functions. In expanded
form, Eq. (7.2.4) becomes

MS5019 – FEM 16

8
 ⎧ u1 ⎫
⎪v ⎪
⎪ 1⎪
⎪u 2 ⎪
⎪ ⎪
⎧u ⎫ ⎡ N1 0 N2 0 N3 0 N4 0 ⎤ ⎪ v2 ⎪
⎨ ⎬=⎢ ⎨ ⎬ ( 7 .2 .6 )
⎩v ⎭ ⎣ 0 N1 0 N2 0 N3 0 N 4 ⎥⎦ ⎪u3 ⎪
⎪ v3 ⎪
⎪ ⎪
⎪u 4 ⎪
⎪v ⎪
⎩ 2⎭

z Step 3 Define the Strain/Displacement and Stress/Strain Relationship

Again the element strains for 2-D stress state are given by

⎧ ε x ⎫ ⎧ ∂∂ux ⎫
⎪ ⎪ ⎪ ∂v ⎪
⎨ ε y ⎬ = ⎨ ∂y ⎬ (7.2.7)
⎪γ ⎪ ⎪ ∂u + ∂v ⎪
⎩ xy ⎭ ⎩ ∂y ∂x ⎭
MS5019 – FEM 17

Using Eq. (7.2.6) in Eq. (7.2.7) and taking the derivatives of u and v as
indicated, we can express the strains in terms of the unknown nodal
displacements as
{ε } = [ B]{d } (7.2.8)
where

⎡− (h − y ) 0 (h − y ) 0
1 ⎢
[ B] = 0 − (b − x) 0 − (b + x)
4bh ⎢
⎢⎣ − (b − x) − (h − y ) − (b + x) (h − y )

(h + y ) 0 − (h + y ) 0 ⎤
0 (b + x) 0 (b − x) ⎥⎥ (7.2.9)
(b + x) (h + y ) (b − x) − (h + y )⎥⎦

MS5019 – FEM 18

9
 From Eq. (7.2.8) and (7.2.9), we observe that ε x is a function of y , ε y is
a function of x, and γ xy is a function of both x and y. The stresses are again
given by the formulation in Eq. (5.3.33) where now [ B ] is that of Eq. (7.2.9)
and {d } is that of Eq. (7.2.1).

z Step 4 Derive the Element Stiffness Matrix and Equations

The stiffness matrix is determined by
h b
[k ] = ∫ ∫ [ B] [ D][ B] t dx dy
T
(7.2.10)
− h −b

where [ D ] given by the usual plane stress or plane strain conditions .

Since the [ B ] matrix is a function of x and y , integratio n of Eq. (7.2.10)
must be performed. The [ k ] matrix for rectangular element is now
of order 8 × 8.

MS5019 – FEM 19

The element force vector is determined by

[ f ] = ∫∫∫[ N ]T { X }[ N ]dV + {P} + ∫∫ [ N ]T {T }dS (7.2.11)
V A

where [ N ] is the rectangular matrix in Eq. (7.2.6), and N1 through N 4 are

given by Eqs. (7.2.5). The element equations are then given by
{ f } = [ k ]{d } (7.2.12)

z Step 5, 6, and 7
Steps 5, 6, and 7, involving assembling the global stiffness matrix and
equations, determining the unknown nodal displacements, and
calculating the stresses, are identical to those in Section 5 for CST.
However, the stresses within each element now vary in both x and y
directions.

MS5019 – FEM 20

10
7.3. Isoparametric Formulation of the Plane
Element Stiffness Matrix
z This element resembles that in Figure 7-3, except that here it can be an
arbitrary quadrilateral as shown in Figure 7-4.
z The natural coordinate system s-t is defined by element geometry and
not by the element orientation in the global coordinate x-y.
z Thus, when the shape function is u = a1 + a2s + a3t + a4st for the
displacement, we use x = a1 + a2s + a3t + a4st for the discription of a
coordinate point in plane element.
z This formulation is general enough to be applied to more-complicated
(high-order) elements such as a quadratic plane element with three
nodes along an edge, which can have straight or curved sides.

MS5019 – FEM 21

y, v
t t
s = + 12 3( x3 , y3 )
Edge t = +1
1 1
t = + 12
4 3 4( x4 , y4 )
1 s
s
1 Edge s = −1
Edge s = +1
1 2 1( x1 , y1 )
2( x2 , y2 )
Edge t = −1
(a) x, u
PARENTS ELEMENT
(b)
ISOPARAMETRIC ELEMENT

Figure 7-4
(a) Linear square element in s-t coordinates and
(b) Square element mapped into quadrilateral in x-y coordinates
MS5019 – FEM 22

11
z Step 1 Select Element Type
¾ First, we consider s-t coordinates are attached to the element, with the
origin at the center of the element, as shown in Figure 7-4 (a). The s and t
coordinate axes NEED NOT be orthogonal, neither has to be parallel to
the x or y axes.
¾ We consider the quadrilateral element that have eight d.o.f, u1, v1, ..., u4,
and v4, associated with the global x and y directions. The element then has
straight sides but is otherwise of arbitrary shape, as shown in Figure 7-
4(b).
¾ As we have shown for a rectangular element, the shape functions that
define the displacement within the element is given in Eq. (7.2.5). These
same shape function will now be used to map the square of Figure 7-4(a)
in isoparametric coordinates s-t to the quadrilateral of Figure 7-4(b) in x-y
coordinates.

MS5019 – FEM 23

The global coordinates x and y then are defined, in matrix form, as

⎧ x1 ⎫
⎪y ⎪
⎪ 1⎪
⎪ x2 ⎪
⎪ ⎪
⎧ x ⎫ ⎡ N1 0 N 2 0 N 3 0 N 4 0 ⎤ ⎪ y 2 ⎪
⎨ ⎬ ⎢= ⎥⎨ ⎬ (7.3.1)
⎩ y ⎭ ⎣ 0 N1 0 N 2 0 N 3 0 N 4 ⎦ ⎪ x3 ⎪
⎪ y3 ⎪
⎪ ⎪
⎪ x4 ⎪
⎪y ⎪
⎩ 4⎭
where
(1 − s)(1 − t ) (1 + s)(1 − t )
N1 = N2 =
4 4
(7.3.2)
(1 + s )(1 + t ) (1 − s )(1 + t )
N3 = N4 =
4 4
MS5019 – FEM 24

12
z The shape function of Eq. (7.3.2) are seen map the s and t coordinates
of any point in the square element of Figure 7-4(a) to those x and y
coordinates in the quadrilateral element of Figure 7-4(b). For instance,
for node 1 with s = -1 and t = -1 becomes x = x1 and y = y1. Similarly
for another nodes.
z We can observe that N1 through N4 have the properties that Ni (i = 1, 2,
3, 4) is equal to one at node i and equal to zero at all other nodes. The
physical shapes of Ni as they vary over the element with natural
coordinates are shown in Figure 7-5.
z Other important characteristic of Ni is that the sum of all shape
functions must equal to one or
n

∑N
i =1
i =1 (i = 1, 2, 3, L, n) (7.3.3)

where n = the number of shape functions.

MS5019 – FEM 25

Figure 7-5
Variations of the shape functions over a linear square element
MS5019 – FEM 26

13
z Step 2 Select Displacement Functions
The displacement functions within an element are now similarly defined by
the same shape functions as used to define the element shape; that is.
⎧ u1 ⎫
⎪v ⎪
⎪ 1⎪
⎪u 2 ⎪
⎪ ⎪
⎧u ⎫ ⎡ N1 0 N 2 0 N 3 0 N 4 0 ⎤ ⎪v2 ⎪
⎨ ⎬=⎢ ⎥⎨ ⎬ (7.3.4)
⎩v ⎭ ⎣ 0 N1 0 N 2 0 N 3 0 N 4 ⎦ ⎪u3 ⎪
⎪ v3 ⎪
⎪ ⎪
⎪u 4 ⎪
⎪v ⎪
⎩ 4⎭
where u and v are displacement parallel to global x and y coordinates,
and the shape functions are given by Eq. (7.3.2).
MS5019 – FEM 27

z Step 3 Define the Strain/Displacement and Stress/Strain Relationship

The displacements are now functions of the s and t coordinates as shown in
Eq. (7.3.4), with the shape functions given by Eq. (7.3.2).

Let f be some function of x and y. For the plane element which is discussed
here, f is either u or v. The chain rule yields
∂f ∂f ∂x ∂f ∂y
= +
∂s ∂x ∂s ∂y ∂s
(7.3.5)
∂f ∂f ∂x ∂f ∂y
= +
∂t ∂x ∂t ∂y ∂t

⎧ f,s ⎫ ⎧ f,x ⎫
or ⎨ ⎬ = [ J ]⎨ ⎬ (7.3.6)
⎩ f ,t ⎭ ⎩ f, y ⎭

MS5019 – FEM 28

14
 where [J] is Jacobian matrix, defined as
⎡ x, s y, s ⎤ ⎡ J11 J12 ⎤
[J ] = ⎢ ⎥=⎢ ⎥ (7.3.7)
⎣ x,t x,t ⎦ ⎣ J 21 J 22 ⎦

The inverse of Eq. (7.3.6) is

⎧ f,x ⎫ ⎧ f , s ⎫ ⎡ Γ11 Γ12 ⎤ ⎧ f , s ⎫

⎨ ⎬ = [Γ]⎨ ⎬ = ⎢ ⎥⎨ ⎬ (7.3.8)
⎩ f, y ⎭ ⎩ f ,t ⎭ ⎣Γ21 Γ22 ⎦ ⎩ f ,t ⎭
where [Γ] = [J]-1
For the four - node element, from Eq. (7.3.1) we obtain
J11 = x, s = N1, s x1 + N 2, s x2 + N 3, s x3 + N 4, s x4
The expression for J12 , J 21 , and J 22 can be obtained similarly using Eq. (7.3.1).
If x = s and y = t , then we have [J] = [Γ] = [I]

MS5019 – FEM 29

We now want to express the element strains as

{ε } = [ B]{d } (7.3.9)
where [ B ] must now be expressed as a function of s and t. Firstly we can write
⎧u, x ⎫
⎧ ε x ⎫ ⎧ ∂∂ux ⎫ ⎡1 0 0 0⎤ ⎪ ⎪
u
{ε } = ⎪⎨ ε y ⎪⎬ = ⎪⎨ ∂∂yv ⎪⎬ = ⎢⎢0 0 0 1⎥⎥ ⎪⎨ , y ⎪⎬ (7.3.10)
⎪γ ⎪ ⎪ ∂u + ∂v ⎪ ⎢0 1 1 0⎥ ⎪ v, x ⎪
⎩ xy ⎭ ⎩ ∂y ∂x ⎭ ⎣ ⎦ ⎪v ⎪
⎩ ,y ⎭
where
⎧u, x ⎫ ⎡ Γ11 Γ12 0 0 ⎤ ⎧u, s ⎫
⎪u ⎪ ⎢ ⎪ ⎪
⎪ , y ⎪ ⎢Γ21 Γ22 0 0 ⎥⎥ ⎪u,t ⎪
⎨ ⎬ = ⎨ ⎬ (7.3.11)
v
⎪ ⎪, x ⎢ 0 0 Γ11 Γ12
⎥ v
⎪ ⎪, s
⎪⎩v, y ⎪⎭ ⎢⎣ 0 ⎥
0 Γ21 Γ22 ⎦ ⎪⎩ v,t ⎪⎭
and

MS5019 – FEM 30

15
 ⎧u1 ⎫
⎪v ⎪
⎪ 1⎪
⎧u, s ⎫ ⎡ N1, s 0 N 2, s 0 N 3, s 0 N 4,s 0 ⎤ ⎪u2 ⎪
⎪u ⎪ ⎢ N ⎪ ⎪
⎪ ,t ⎪ ⎢ 1,t 0 N 2 ,t 0 N 3 ,t 0 N 4 ,t 0 ⎥⎥ ⎪v2 ⎪
⎨ ⎬=⎢ ⎨ ⎬ (7.3.12)
⎪v, s ⎪ ⎢ 0 N1, s 0 N 2, s 0 N 3, s 0 N 4, s ⎥ ⎪u3 ⎪
⎪⎩ v,t ⎪⎭ ⎣ 0 ⎥
N1,t 0 N 2 ,t 0 N 3, t 0 N 4,7 ⎦ ⎪ v3 ⎪
⎪ ⎪
⎪u4 ⎪
⎪v ⎪
⎩ 4⎭

The matrix [ B ] is the product of successive rectangular matrices in Eqs. (7.3.10),

(7.3.11), and (7.3.12).

MS5019 – FEM 31

z Step 4 Derive the Eelement Stiffness Matrix and Equations

We now want to express the stiffness matrix in terms of s-t coordi-nates. For
a constant thickness element, we have
[k ] = ∫∫ [ B]T [ D][ B] t dx dy (7.3.13)
A

However, [B] is now a function of s and t, so we must integrate with respect

to s and t. We need to transform teh variables and the region from x and y to s
and t. This general type of tranformation is given by

∫∫ f ( x, y)dxdy = ∫∫ f (s, t ) [ J ] ds dt
A A
(7.3.14)

Using Eq. (7.3.14) in Eq. (7.3.13), we obtain

1 1
[k ] = ∫ ∫ [ B] [ D][ B] t [ J ] ds dt
T
(7.3.15)
−1 −1

The determinant [ J ] is a polynomial in s and t , and is tedous to evaluate even

for the simplest case of the linear plane element.
MS5019 – FEM 32

16
7.4. Introduction to the Gauss Quadrature
We will describe Gauss’s method, one of the many schemes for
numerical evaluation of definite integrals, because it has proven most
powerfull for FE work. The concept of the Gauss quadrature is first
illustrated in 1-D in the context of approximating the integral
1
I= ∫ f ds
−1
(7.4.1)

where f = f ( s ). By the Gauss integration procedure, the integration

represented by Eq. (7.4.1) can be approximated by
1 n
I= ∫ f ds ≅ ∑ Wi f i = W1 f1 + W2 f 2 + L + Wn f n (7.4.2)
−1 i =1

Thus, to approximate I, f = f(s) is evaluated at several point of

locations of fi. The resulting fi are multiplied by an weight factor Wi,
and the values added.
MS5019 – FEM 33

The Gauss method locates the sampling points (or Gauss points) so that,
for a given number of them, greatest accuracy is achieved.
The sampling points and weights for the Gauss quadrature up to order n =
4 are given in Table 7-1. It is clear that the sampling points are located
symmetrically with respect to the centre of the interval. This symmetry is
illustrated for 1-D problem in Figure 7-6. In general, a polynomial of
degree (2n-1) is integrated exactly by an n-point Gauss quadrature.
Table 7-1 Sampling points and weighing factors for
Gauss quadrature numerical integration
Order n Location, si Weigth, Wi
1 s1 = 0.000... 2.0
2 s1,2 = ± 0.577350269189626 1.0
s1,2 = ± 0.774596669241483 0.555555555555556
3
s3 = 0.000... 0.888888888888889
s1,4 = ± 0.861136311594053 0.347854845137454
4
s2,3 = ± 0.339981043584856 0.652145154862456
MS5019 – FEM 34

17
 f f

1 1 2

s s
−1 +1 −1 +1
(a) One sampling point (b) Two sampling points
Figure 7-6 Sampling points of the Gauss quadrature in one dimension

In 2-D problems, where f = f(s,t), the Gauss quadrature is applied by

integrating with respect to s and then with respect to t as follow:
1 1 1
⎡ ⎤
I= ∫ ∫ f ( s, t ) ds dt = ∫ ⎢⎣∑W
−1 −1 −1 i
i f ( si , t )⎥ dt
⎦
⎡ ⎤
= ∑ W j ⎢∑ Wi f ( si , t j )⎥ = ∑∑ WiW j f ( si , t j ) (7.4.3)
j ⎣ i ⎦ i j

MS5019 – FEM 35

In Eq. (7.4.3), we need not use the same number of Gauss points in ecah
direction (that is, i does not have to equal j), but this is usually done. Thus,
for example, a four-point Gauss rule (often describe as a 2 x 2 rule) is
shown in Figure 7-7. Equation 7.4.3 with i = 1, 2 and j = 1,2 yields
I = W1W1 f ( s1 , t1 ) + W1W2 f ( s1 , t 2 ) + W2W1 f ( s2 , t1 ) + W2W2 f ( s2 , t 2 ) (7.4.4)
where the four sampling points are at si , ti = ± 0.5773... = ± 1 3 , and the
weight are all 1.000.
s = −0.5773...(i = 1) t s = 0.5773...(i = 2)

( s1 , t 2 ) ( s2 , t 2 )
t = 0.5773...( j = 2)
3 4
s
( s1 , t1 ) ( s2 , t1 )
t = −0.5773...( j = 1)
1 2

Figure 7-7 Four-point Gaussian quadrature in two dimension

MS5019 – FEM 36

18
In general, in three dimension, we have
1 1 1
I=∫ ∫ ∫ f (s, t , z ) ds dt dz = ∑∑∑W W W
−1 −1 −1 i j k
i j k f ( si , t j , z k ) (7.4.5)
1 n
Now if the limits are ∫ f ( x)dx = ∑ Wi f ( xi ), then the weights Wi , and locations
0 i =1

xi are different than for limits - 1 to + 1 (as listed in Tabel 7 - 2), but the procedures
are the same. In Table 7 - 2, the weights are W1 = W4 and W2 = W3 .
Table 7-2 Table of Gauss points for four-point
Gaussian quadrature for integration from 0 to 1.

Locations, xi Associated Weight, Wi

0.06943185 0.1739274
0.3300095 0.3260725
0.6699905 0.3260725
0.9305682 0.1739274
MS5019 – FEM 37

7.5. Evaluation of the Stiffness Matrix by

Gaussian Quadrature
We have shown in Section 7.3 that [k] for a quadrilateral element can
be evaluated in terms of a local set of coordinates s-t with limits from
minus one to positive one within the element, and in term of global
nodal coordinates as given by Eq. (7.3.15) and repeated here for
convenient as
1 1
[k ] = ∫ ∫ [ B]
T
( s, t )[ D ][ B]( s, t ) [ J ] t ds dt (7.5.1)
−1 −1

In Eq. (7.5.1), each coefficient of the integrand [ B]T [ D ][ B ] [ J ] must be

evaluated by numerical integration in the same manner as f ( s,t ) was
integrated in Eq. (7.4.4)
A flowchart to evaluate [k] of Eq. (7.5.1) using four-point Gaussian
quadrature is given in Fugure 7-8.
MS5019 – FEM 38

19
 Start

Read in four Gauss points and weights

si, ti = ± 0.5773...; Wi, Wj = 1.00, 1.00

Zero [k(e)]

Do i =1, 4

Compute J ( s, t ) , B( s, t ), D

Compute k = BT D B J t

Compute k ( e ) ← k ( e ) + k Wi W j

End

Figure 7- 8 Flowchart to evaluate k(e) by four-point Gaussian quadrature

MS5019 – FEM 39

Reference:
1. Logan, D.L., 1992, A First Course in the Finite Element Method,
PWS-KENT Publishing Co., Boston.
2. Imbert, J.F.,1984, Analyse des Structures par Elements Finis,
2nd Ed., Cepadues.
3. Zienkiewics, O.C., 1977, The Finite Eelement Method, 3rd ed.,
McGraw-Hill, London.
4. Finlayson, B.A., 1972, The Method of Weighted Residuals and
Variational Principles, Academic Press, New York.

MS5019 – FEM 40

20