Fourier analysis, yielding smoothing and convergence factors which give an idea of the ... A discussion of storage andwork requirements is also included.
/,u-/y22'_"
NASA
Technical
{N_SA-TM-87761) S_RUCTURAZ HC A04/MF
Memorandum
_ULTIGRID
MECHANZCS A01
[NASA)
METBODS 60
87761
IN
N86-28465
p CSCI.
20K
G3/39
MULTIGRID METHODS IN STRUCTURAL MECHANICS
I. S. Raju, C. A. Bigelow, S. Ta'asan and M. Y. Hussaini
July 1986
National Aeronautics and Space Administration Langley Research Center Hampton, Virginia 23665
Unclas 43401
_7
SIJMMAI_Y
Although elasticity have
and
been
are
algorithm
on
the
are
finest
grid
level,
algorithm
yielded
and
better
sequential
of
prolongation an
convergence
rates. the
to
ordering
prolongation required
balance With
multigrid
equivalent
red-black
energy
yielded work
on
standard
of
the
of
applications
multigrid deflections
of
the
of
a
multigrid
the
than
fine
with
cycle.
prolongation the
inter-grid
machine
equivalent
to
about
of
accuracy 4
iterations
1
t
in
was
based
on required
to
improved
fine
with
grid.
a The
energy-conserving
a single the
led
2
operator.
accuracy
the
or
on
than
sequential
machine
linear
200
greater
also
on
to
transfers
and
iterations either
n
restriction
transfers
results
with
of
the
prolongation
Gauss-Seidel
solutions
of
used.
linear
operators inter-grid
the
algorithm
Also
V(n,1)
elements
were
equivalent
with
32
test
work
the
yielded
relaxation
with
cycle
transpose
To
multigrid
V(1,n)
be
in the
the
grid.
during
the
beam.
the
and
with
approximations
ordering, accuracy
levels,
the
initial
energy-conservlng
45
such
the
aspects
half
energy
algorithm
few
of
coarseness)
random
restriction
operator
equations
work,
analysis
model
ordering,
Conserving
Maintaining
which
the
to
machine
rates
present
various
(or
sequential
of
the
suggested,
the
and
to
in detail.
used
iterations
principles.
ordering,
and
been
element
beam,
severely,
convergence
Derivation
the
finite
were
methods
In
explained
were
Gauss-Seidel
prolongation
work
more
which
has
grid-fineness
prolongation
results
standard
that
and six
linear
the
Bernoulli-Euler
studied
multigrid
literature.
to
study,
yielded
energy
in the
applied
of
mechanics
this
multigrid With
structural
supported
In
application
reported
techniques simply
the
V(1,1)
finest
cycle,
grid
level.
INTRODUCTION Multigrid
methods are known to be efficient
certain classes of partial fLeld of fluid application
differential
solution
techniques for
equations, and their advantages in the
dynamics have been clearly
to the equations of elasticity
established
[1,2].
and structural
Although their
mechanics has been
suggested, there have been few attempts to implement multigrid
methods in
these cases. Certain basic ideas underlying multigrid of Southwell [3].
However, the work of Fedorenko [4] should be regarded as
the forerunner of the multigrid difference
methods can be found in the work
methods. Fedorenko solved the elliptic
equations, using concepts of error smoothing by Jacobi relaxation
and calculation
of corrections
on coarser grids.
He also obtained an
asymptotic estimate of the computational complexity for the Poisson equation on a rectangle.
Fedorenko's technique was generalized by Bakhvalov [5] to
include the advection diffusion
equation.
[6] proposed a two-grid method for elliptic methods remained virtually work of Brandt [7], 1970's.
Specifically,
of multigrid
At about the sametime, Wachpress systems.
However, multigrid
unused for almost a decade until
the pioneering
Nicolaides [8] and Hackbush[9] revived them in the midBrandt's early work established
methods, and demonstrated the capability
the actual efficiency of these methods to
treat nonlinear equations, general domains, local grid refinements, and solution
adaptivity.
He demonstrated that the maximumnumber of computer
operations required to solve a discrete system of n equations in a Poisson problem was on the order of n computer operations.
Brandt also showedhow a
local Fourier analysis could be used for the theoretical smoothing rates and optimization
of procedures.
investigation
of
The work of Nicolaides [8] is
the first
systematic study of multigrid
element discretization
of elliptic
background of multigrid
procedures relating
equations.
to the finite
A complete historical
methods may be found in Stuben and Trottenburg [10],
and Hackbush [11]. Initially,
multigrid
methods were successfully
equations to which they are particularly been extended to parabolic Developments in multigrid
applied to elliptic
well suited.
[12] and hyperbolic
Recently, their
use has
[2] equations as well.
methods since the mid-1970's can be found in
proceedings edited by Hackbush and Trottenburg [13],
Paddonand Holstein
[14],
and McCormickand Trottenburg [15]. Manyproblems in elasticity elliptic fields
equations; hence, it remains largely
and structural
mechanics are governed by
is quite surprising
unexplored.
use in these
The present work is an attempt to
introduce these methods to elasticians mechanicists who have had little
that their
and computational structural
or no exposure to multigrid
techniques.
The
focus of the present work is on a simple, one-dimensional example of a Bernoulli-Euler Fourier analysis,
beamfor the following yielding
idea of the efficiency description
reasons:
it
is amenable to local
smoothing and convergence factors which give an
to strive
of various multigrid
for;
it permits the derivation
and
elements without the complications of higher
dimensions; and it serves as a building
block for implementing the algorithm
in higher dimensions. Specifically, the relevant iterative equations.
this paper describes the finite
element discrctizdtLon
of
beamequations and discusses the performance of the Gauss-Seidel
technique for the solution A particular
convergence (usually
of the resulting
geometric multigrid
called the correction
linear
system of
schemefor accelerating scheme) is presented.
the
A general
procedure for defining as restriction) prolongation)
the fine-to-coarse
and coarse-to-fine
transfer
transfer
of residuals
of corrections
is described from a physical view-point.
(referred
(referred
to
to as
A discussion of
storage and work requirements is also included. The convergence characteristics floating
of the solution
and the number of
point operations (work) necessary to achieve convergence are
presented.
Results are given for two different
and for two different of it_rdtions
prolongation
at each grid level
orderings of the iterations
schemes. The effect
of varying the number
is also studied. NOMENCLATURE
Fourier
series
coefficients,
bk
element
length
on
level
B
semi-bandwidth
of
stiffness
EI
flexural
rigidity
value
of
direct
error
in
iteratlve
A
.
J
d
.
J
accuracy
e m
ew.
, i
error
eo.
of
in w
(k =
= I to
N)
I to M)
matrix
at
node
j,
(j =
F
force
vector
8 solutions
at
node
i,
(i
by
maximum
displacement
th m
component
h
element
K
structural
[k]
element
k
entry
of
length
on
F finest
stiffness
stiffness in m
half-length
column, of
level
matrix
matrix
th
mj
÷ I)
z I to
N + e
normalized
m
e
machine
and
error
f
I to N
approximation
i
e
L
k,
solution
the
(j
beam
.th 3 row
of
matrix
K
I)
M
total
M0
concentrated
M,
nodal
number
of
grid
levels
applied
moment
moments
3 N N
number
of
degrees
number
of
elements
of
freedom
e
Ni
shape
nI
number leg of
of iterations the V-cycle
performed
on
the
downward
n2
number leg of
of iterations the V-cycle
performed
on
the
upward
P
prolongation
matrix
PO
concentrated
applied
qo
maximum
value
q(x)
loading
function
R
restriction
R. 3
nodal .th l
r i r W.
3
,
functions,
to
beam,
(0
_ x
_ L)
Q residuals
at
r_sidual node
j,
vector (j =
I to N
+ I) e
3
solution
V.
applied
the
U
V
q(x)
of
L2-norm
1
force
matrix
component
llrli2
U.
of
4)
forces
w and
r@.
(i : I to
of
the
residuals
vector
.th l component
of
U
approximation
to
solution
ith
of
V
component
vector
U
1
WM
work
done
by
concentrated
W
work
done
by
nodal
Wp
work
done
by
concentrated
w
displacement
(%
arbitrary
moment
forces
n
constants
5
force
in displacement
function
or' first
or
second
rotation strain
ir
Superscripts f
energy
:
c fine
and
coarse
grid
iteration
cycle
transpose
of
functions,
respectively
S
T a matrix
FORMULATION
OF
THE
One-Dimensional Consider distributed
a simply loading
symmetric
about
considered.
supported
q(x),
the
as
center
Moment
of
the
equilibrium
d2w
of
PROBLEM
in Figure beam,
Problem
length
only
requires
I
vector
Beam
beam
shown
or
2L
I.
Since
one
half
where the
El
is
loading
the
qo Lx
flexural
function
The
boundary
The
exact
at
rigidity
for
qoxL 3
To
solve
idealized displacement
by
the
a number w and
of
0
of
center for the
the
of the
rotation
the
beam
is needs
to
be
-< x SL
(I)
qo
is
the
maximum
value
of
beam.
simply
loading
and
supported
shown
beam
in Figure
are
w(O)
= w(2L)
= O.
I is
2
using beam
beam
the
4
problem
problem
6L 2
conditions
solution
of
the
x2
( 2
the
the
to
that
) EI --dx 2 =
subjected
the
finite
elements.
dw/dx
(=
element For
9) are
the used
method, two-noded as
the
the
structure
beam
nodal
element,
parameters.
is the The
relevant
boundary conditions are w(0) = 0 and _(L)
matrix [k],
O.
The element stiffness
obtained from the shape functions given in Appendix A, for a beam
element of length b is given below.
12/b 3 [k]
6/b 2
6/b 2
= EI
-12/b 3
element
stiffness
matrix
K and,
for
solved
can
expressed
be :
KU
any
-6/b 2
2/b
-6/b 2
12/b 3
-6/b 2
2/b
-6/b 2
4/b
matrices
are
assembled
loading
and
boundary
given
in matrix
vector,
and
is
the
Conventional
equation
such (4).
vector
finite as
thus,
iterative
technique.
an
of
unknown
element
widely
iteration.
Gauss-Seidel
Iteration with
Gauss-Seidel component
the
stiffness system
alternative
iterative
U 1 as
iteration of
to
be
equation:
matrix,
F is
displacements
and
typically or
use
Cholesky
a direct
time
can
to
solve
become efficient
Jacobi
iteration
Techniques
The
Gauss-Seidel
method
is used
the
initial
approximation
following
7
solution
is an
are
U:
rotations.
solution
techniques
the
load
decomposition,
solution
uses
the
direct
solution
to
Solution
th m
following
discretization,
attractive
used
nodal
elimination
increasing
Gauss-Seidel
U,
a structural
conditions,
square
analyses
Iterative
Starting
the
symmetric,
Gaussian
With
critical;
Two
as
into
(4)
is a positive-definite,
techniques,
form
(3)
F
wher'e K
U
6/b 2
4/b
6/b 2
The
-12/b 3
to
the
algorithm
to
in
the
vector
present
of
calculate
and work.
unknowns s+l um , the
u
s+l
1 = -_ [ k mm
m
where
u
s and m
iteration
the
solution
To shown
The
the
m
N _
-
s .u.] mJ j
k
j=m+l
components
iteration
is
less
is
or than
the
assuming as
convergence
I, with
The
calculating
(m
= 1,2 ....
of
U and
continued
until
the
some
given
F,
and
until
change
s
,N)
the in
indicates
(5)
required
the
the number
current
of
estimate
of
tolerance.
one
of
the
Gauss-Seidel
method,
the
half
the
beam
with
elements,
convergence L2-norm
all
of
of the
elements
are
the
modeled
Gauss-Seidel
residuals, of
equal
iteration
denoted length
32
as b,
beam
was
llr112.
the
was
monitored In
L2-norm
problem
by
discrete
can
be
follows:
IIri12 --
the
u s+l mj j
Study
considered.
where
k
j=l
completed
vector
Figure
expressed
m
the
been
study
in
form,
are
m
has
Convergence
m-I }.
-
th f
number.
iterations
f
residual
r
> r2i
N i:I
(6)
is defined
i
as
follows:
N
r1 : fl-
In be
the
thought
will
have
element
L k..u.
j=1
present of
the
length.
as
finite
the
same
(i = I to
N)
(7)
zj j
element
unbalanced
dimensions,
A non-dimensional
context,
forces the
and
moment L2-norm
the
residuals
moments. residuals can
then
So
that
were be
in each all
equation the
normalized
computed
as
can
residuals by
b,
follows:
the
l lrl te
where by
r w and
the
_
_b_e
r 8 are
factor
2 wi(X)÷ rr°!x)/b 2iJ /(q°L) I
-
the
w
and
q0 L
(the
total
the
Gauss-Seidel
Figure
2 shows
O residuals. load
on
the
The beam)
L2-norm to
has
produce
(8)
also
been
divided
a dimensionless
quantity. To was
start
chosen.
iterations
for
performed of
the
other
two
These
procedure,
all
the
are
designated
then
dashed
black
line
first
50
for
iterations,
results
are
so-called
is
this
I (at
the
shown
by
are of
both
of
for still
llne
was
iterations
were
33
the
center
2.
In the
(at
in Figure
used.
In
this
and
all
the
odd-numbered
performed
on
all
the
red
the
the
red-black
the
ordering
L2-norm
results
than
the
However,
even
after
quite
in error.
high,
B
first,
shown
although
by
the
in the the
sequential 1600
indicating
Appendix
nodes
rapidly
slows,
is
inherent
the the
these
high
low
frequency
frequency
levels high
in the
where
frequency
Gauss-Seidel errors
errors, they errors
method,
but since
become can
not
the
they higher
then
be
are
U
ordering,
nodes
are
drops
convergence
better
is still
grossly
the
for
of
red
iterations.
both
number
to node
solid
ordering
approximation
the
procedure,
support)
methods,
slightly
versus
the
designated
initial
that
further
the
examines
the
problem.
on coarser
levels,
first
thereafter,
numbers
L2-norm
However,
For
produces
removes
approximated
2.
plotted
the
results
cycles;
large
behavior
efficiently
are
The
in Figure
of
This
nodes
a random
red-black
nodes.
solution
convergence
coarser
node
Iterations
the
iterative
In
black.
ordering
ordering
errors.
the
iteration
red-black
L2-norm
from
even-numbered
the
the
procedures.
sequentially
beam).
iteration,
since low
Gauss-Scidel
frequency
smooth,
frequency removed
can
be
errors. more
On
efficiently.
Even if
iterations
the errors are still
smooth on coarser levels,
on the coarser levels require much less work due to the decrease in
the number of unknowns. These ideas are exploited ir_ a recursive
in the multigrid
procedure
manner to improve the convergence rate. MULTIGRIDPROCEDURES
In this section, terminology,
the multigrid
relaxation
method is described.
is used to meaniterations
In multigrid
of some smoothing
technique, such as the Gauss-Seidel described above; thus, to perform relaxations
on a set of equations meansto iterate
Often, the terms relaxation multigrid
literature.
and iteration
are used interchangeably
In the remainder of this
used to meanGauss-Seidel iteration,
using a smoothing method. in
paper, the term relaxation
and n relaxations
meansn iterations
is on
the system of equations. The multigrid
techniques are iterative
grids and a simple relaxation stated, such relaxation errors on a given grid.
schemesuch as Gauss-Seidel.
is very efficient
As previously
for removing the high frequency
The remaining low frequency errors becomehigher
frequency errors on coarser grids, The key elements of the multigrid the coarse grid correction. generally
methods that use a sequence of
where they can be effectively procedure are the relaxation
The multigrid
called a geometric multigrid
smoothed. technique and
algorithm described here is
method.
Coarse Grid Correction Consider the interplay
between a coarse grid and a fine grid,
coarse grid has half as many node points as the fine grid. equation be the discrete a
fine
finite
element representation
grid.
10
where the
Let the following
of the given problem on
Kfu f _ If
Vf
is the
sweeps,
Ff
(9)
approximation
then
the
error
obtained
e f _ U f.-
on
the
fine
V f satisfies
grid
the
after
a
following
few
relaxation
equation:
Kfe f _ r f where
rf
(10)
the
residual
on
If e f is
a smooth
function,
function
e
c
which
the
satisfies
fine
grid
is defined
e f can
the
be
as
r f _ F f - Kfv f
approximated
following
by
a coarse
grid
equation:
KCe c _ r c where
rc Here
grid
has
solve
R
= R is
half
equation
some
method,
fine
grid
the
P
the
fine
right
many
(F f - Kfv f)
to
coarse
grid
grid
points
as
than
(10).
approximate
solution
is the
_ R
fine
(11)
V f ÷ Vf Here
(r f)
the as
(11)
using
the
and
ec
fine
operator.
grid,
it
Since
is more
solving
equation
(11)
is used
to accelerate
the
coarse
economical
to
approximately convergence
by of
the
the
following: (12)
+ p(e c)
prolongation
grid
the
After
error
restriction
operator
arrow
÷ means
that
interpolates
that
the
left
from
side
the
coarse
is replaced
grid
by
to
the
side. c The
results
process to
the
equation
(9)
multiple
grids.
of fine
calculating grid
efficiently, This
is
r
c , solving
called
the
the
above
recursive
for
coarse
procedure
procedure
e
, and
grid is
interpolating
correction.
applied
is described
To
recursively in
the
the solve using
following
section. Multigrid Figure supported
3 shows beam
where
the the
sequence element
of
Cycle
grid
size
bk
11
levels, satisfies
k
_ I to the
M,
for
relation
a simply b k = 2k-lh
and
h is the element size on the finest k is written
grid level
(k =I).
The equation on level
as follows:
Kkuk = Fk
(13)
For k = I (the finest stiffness
matrix
algorithm
must
on
the
the
finest
and
progressively
solutions
in
the it
technique
is
the
Starting
on
(13),
If e I denotes through
the
k
the
= I.
coarser
in practice,
equation
right-hand
determine grid,
level),
form is
the
error
following
levels
the
out
means
by
problem.
is
construct This of
, respectively
algebraic idea
to
= F
original
fundamental
(13).
The
system
to
use
other
procedure
relaxation.
basic
solutions
finite
all
subsequent
residual
from
level, the
k
= I,
let
V 1 be
on
the
finest
residual
U 1 - V I,
rk and,
= R
therefore,
called
smoothing;
A commonly
used
the
error
and
the
an
approximation
level
as
residual
to
r I = F I - KIv I .
are
related
(k
(r k-1
> I),
level
the
residual
r
k - I defined
by
k
is
the
the
restriction
following
of
- Kk-1 e k-1 )
equation
to
be
cycle
for If
obtaining
an
(15) solved
is
now
the
error
equation:
k = I, an
(16)
approximate
improving
the
equation:
Kke k = r k Given
U I in
(14)
previous
the
from
equation:
levels
the
(13)
element
KIeI = r I For
given
in equation
the
is
the
method.
finest
define
the
of
equation
Gauss-Seidel
and
of
solution The
carried
the
side
grid
of
K1 _ K and F!
solution
this
Vk
approximation
perform
n I relaxation
approximate
solution
to
the
can sweeps
V I.
equation
be
executed
on
Calculate
12
(13),
recursively
equation the
the
(13),
residual
multigrid as
follows.
thereby r I _ F I - KIv I .
Down-Sweep: For these
I < k
steps
< M,
repeat
defined
(I)
On
as the
steps
(I)
and
(2);
for
k
_ M,
go
to
step
(3),
wit_
follows: current
level
k,
start
with
an
initial
approximation
k e
= O,
and
obtaining (2)
a new
Restrict
(3)
perform
approximation
residual
sweeps
the
next
r k÷1
_ R( r k - Kke k ), and
set
k
If k
= M,
(16)
solve
equation
on
equation
(16),
e
to
the
the
n I relaxation
coarser
level,
k +
I, where
= k + I. exactly,
set
k
_ M
-
I, and
stdrb
up-sweep.
Up-Sweep: For where
M
> k
these
> I,
steps
repeat
are
steps
defined
(4)
as
and
(5);
when
k
= I,
go
to
step
(6),
follows: k
(4)
Update by
the
adding
approximation the
previous
prolongated
coarser
ek + ek (5)
Perform
to
the
error
correction
e
on
of
the
the
current
error
from
level the
level:
+ P (ek÷1 )
n 2 relaxations
(17) on
equation
(16),
starting
with
the
k updated
error
e
from
equation
(17).
This
will
result
in _n
k improved (6)
For VI
The
indicated
cycle
by
relaxation V-cycle
calculate
+ V I + p(e2), an
described
above
the
n2
k = I,
obtaining
sweeps and
approximation
notation
indicates
the
an
and
improved
at
each
number
e
.
Set
updated
perform
k
= k
- I.
approximation
n 2 relaxations
approximation
in steps
V(nl,n2)
performed
to
(I)
to
through
on
(6)
indicates
the
level
on
the
or
13
solution
equation
is called
nI
relaxation
the
(14),
V I.
, where
of
to
first sweeps
number downward
performed
a V-cycle,
of leg at
of each
th_
leve_
on
drawing shown
the of
second
or
a V-cycle
in Figure
upward
for
leg
six
grid
that
the
the
operators.
to
the
These
From are
foregoing
keys
as
and
the P
ef
multigrid
(16),
in
the
and
f
as
the
the
procedure
are the
it can grid
one
schematic
V-cycle
is
to
residual
Matrices
the
each
scheme,
restriction other
relations
on
as the
it
is apparent
and
prolongation
explained coarse
below.
and
fine
grid
be
represented can
(]8) on
be related
a coarser
level.
using
prolongation
the
In
that
operator
= Pe c
(19) r
c and
r
can
be
related
through
the
restriction
operator
R
follows: c
_ Rr
f (2O) c
equations
(18),
the
errors,
e
f and
e
, can
be
thought
of
as
representing c
the be
case,
equation:
residuals
r In
for
correction
are
related
errors
f and
A flowchart
a
Kfe f = r f
fine
following e
4 shows
follows:
ks smooth,
coarse
Figure
Prolongation of
KCe c = r c If
cycle.
levels.
description
operators
equations
written
the
5. Restriction
From
of
displacements thought
of
of
as
the of
the
strain
energy
the
erro_
equations
c
Fro_ fine
grids
respective
forces the by
required
coarse
the
and
following
grids,
be
(18),
produce
fine
grids
rewritten
as
and
the
residuals,
these can
be
r
displacements. expressed
in
f and
r
, can
Thus, terms
of
expressions:
_
19)
while
to
I _ (eC)Tr c
equations may
the
f
(20),
follows:
14
I = _
(e f)
T f r
the
strain
(21)
energies
of
the
coarse
and
c
I . f _ _ _eC)TRr (22)
f
I
eC)TpT f
If the strain must equal _f .
energy is to remain constant during intergrid
e
transfers,
values of ec and r f, R must equal pT to
Thus, for arbitrary
maintain an energy equivalence between grid levels. One-Dimenslonal BeamExample Before applying the multigrid the restriction Restriction
cycle to the solution
matrix R and the prolongation
matrix P must be defined.
Matrix
The restriction
matrix R is defined to transfer
from a fine grid to a coarse grid. relations
of the beamproblem,
To form this matrix,
must be formulated for transferring
equations (A2) and (A3), it
the loads and moments the consistent
load
the loads and the moments. With
is possible to determine the nodal forces in the
beamelement equivalent
to a concentrated force and momentapplied at the
center of the element.
The details
Concentrated
PO
applied
work
done
at by
the the
force.
The
center
of
an
concentrated
of this formulation
nodal
forces
element force
to
equivalent
of
length
the
work
are given below. to
b are done
a concentrated
found
by
the
by
equating
unknown
force
the
nodal
forces.
Wp, times
the
the
done
displacement
We The
work
= P0
substitution
by
the
of
that
concentrated
force,
is
equal
to
the
force
force:
× w(x_b/2) of
x
_ b/2
PO
(23) into
equations
(A2)
and
(A3)
yields
the
following
result:
Wp
= P0X{Wl/2
+
61(b/8)
+ w2/2
15
- e2(b/8)}
(54)
Similarly,
the work done by the nodal forces acting at the two end nodes
can be written
as follows:
Wn_ RlWl + R2w2÷ MI@ I + M2@2
(25)
where RI, R2 are the nodal forces and MI, M2 are the nodal moments. Wpshould equal Wn for any arbitrary the consistent
values of Wl, w2, 01, and 02.
Thus,
forces are
RI _ R2 = PO/2 M1 =Pob/8
Concentrated concentrated
moment
equating
work
unknown
the nodal
W M, the
the
moment
To respect
M 0 and = M0
determine to
×,
@
moment. M0
As
before,
applied
done
(26)
at
the
the
nodal
center
of
forces the
by
the
concentrated
moment
by
the
concentrated
moment,
equivalent
beam
to
the
element work
to a are
done
found
by
by
the
forces.
work
WM
M2 = -M I
done the
rotation
dw dx
by
that
to
the
product
8
in an the
Wl[-
6x b2
of
moment:
× @(x=b/2)
yielding
....
caused
is equal
(27)
element,
equdtion
following
+
6x 2 ] b3
(A2)
must
be
differentiated
with
equation:
+ el[l
-
4x 3x 2 ] __ + b 7
(28) +
Substituting
x = b/2
concentrated
moment
WM
,6___x W2[b2
into M can
= Mo×[(-3/2b)w
6x 2] b3
equation be
+
(28)
written
I -
@I/4
@2[_
as
+
and
2x --_
3x 2] + b2
simplifying,
the
_ x
work
_ b
due
to
the
follows:
(3/2b)w 2 -
16
0
@2/4]
(29)
As written
where
for
the
as
be
+ MI01
+ M202
R2
are
forces
and
WM
the
nodal
should
equivalent as
RI
equal nodal
Wn
for
forces
= - 3M0/(2b
)
residual
in equations in
the
To
form
consider
two
node
to
the
nodal
forces
is
points
one
element.
and
b,
the
arbitrary the
case
(26)
the
and
restriction
grid
levels
two
The
(31)
of
The
and
moments
from
forces
and
moments
is
M2
= MI
are
nodal
values of
the
moments.
of
w I, w 2,
concentrated
el,
and
moment
to
shown
using
can
the
the
fine
the
in Figure
force
6.
coarse
elements
grid
determine
forces
the
and
moments
residual
weights
transfer.
in Figure
while
nodal
in
Here grid
the R
to
coarse
6 and
moment
the has
fine
matrix the
and
fine two
and
is defined grid.
defined
weights, grid
node
coarse to
has
points
and
grids
restrict
The
three
are
transfer
of
below.
where f
T = {R i
Mi
Rj
Mj
Rk
Mk }
17
c r
T -
{Ri,
Mi,
Rk,
Mk,}
b/2
the
r c = R[r f]
r
O2 .
(31)
consistent
used
restriction
the
shown
The
matrix
as
lengths
forces
= -R I
inter-grid
elements,
respectively.
R2
transfer.
energy-conserving
and
(30)
follows:
Consistent
use
due
M I, M 2 are
any for
M I = - MO/4
fo_
work
= R1w I + R2w 2
expressed
given
the
Wn
Again, the
force,
follows:
RI,
Thus,
concentrated
the
R
This
_
0 0I
I/2 I/2 b/4
0
0
-b/4
restriction
coarser
matrix
grid
The
prolongation
displacements
the
two
grid
on
the
same
as
shown
coarse
grid
at
nodes
as
defined
nodes
method,
only
0 0 -] Ol
-I/4
0
I
when
the
at
that
node
the
i'
and
defines
coarse
in
at
i,
Figure
and
k
the
element
size
of
j
on
(32)
is
on
the
Here
doubled
the
these
on
the
fine
grid
can
i'
and
mesh
node
j
mesh.
k'
i and
on
k
the
are
found
are
the
This
fine
at
each
w
be
found
the
coarse
two
average
amounts
to
given
consider
and
% are
displacements,
of
in
grid
again
displacements
From
nodes
the
matrix,
k'.
fine
coarse
this
and
at
at
obtain
on
i'
nodes
the
displacements
To
6.
displacements
displacements _'
the
grid.
nodes
j
displacements
displacements
assumes
the
levels
at
the
matrix on
displacements
either
valid
0I 0
Matrix
the
known
is
3/(4b) -3/(4b) -I/4
level.
Prolongation
In
0 0I
fine
grid.
ways.
of
the
linear
in
the
two
grid
ways.
are
the
The
One
method
displacements interpolation
below:
c w
= PLw f ]
w f
=
where
p
{ w i
0i
w 3
0j
Wk
8k}T
I
0
0
0
0
I
0
0
I/2
0
I/2
0
0
I/2
0
I/2
0
0
I
0
0
0
0
]
wc
=
{ wl,
0i,
Wk,
8k,}
T
(33)
18
in the and
(A3)
the
element
are
equations mesh.
other used
method,
the
for
interpolation.
in the (A2)
The
the
coarse_mesh
and
(A3).
prolongation
e
matrix
of
use
grid
levels.
size
is doubled
I/2
b/4
I/2
-b/4
I/4
3/(4b)
-I/4
0
I
0
0
I
-_0
matrix
obtained
of
restriction
the
at
this each
restriction
problems.
The direct degrees matrix
K.
solution also
of
be
node
center
of
into
j on
the
fine
shape
functions
in equation
energy
(32).
equivalence
is valid
only
when
is Thus,
between the
element
level.
prolongation
procedure
used
operators here
can
derived be
here
applied
to
are two-
Requirements
storage
solution
the
matrix
the
(34)
element
an
at
(A2)
below:
given
maintains
grid and
the
matrix
prolongation
coarser
three-dimensional Work
matrix
using
in equations
x = b/2 at
given
0
problem,
and
is
0
a one-dimensional
Storage
obtained
I
prolongation
the
displacements
0
Again,
Although
substituting
0
transpose this
thus
the
by
given
displacements
0
prolongation
the
The
found
are
functions
0
- 0
the
shape
I
0
exactly
are
These
-3/(4b)
The
element
requirements techniques
freedom Since
technique,
in the
necessary are
the
problem
multigrid the
easily
storage
and
method
to
solve
calculated the
addressed.
19
equation
from
semi-bandwidth
is being
requirements
the
proposed of
the
the
KU=F
total of
as
the an
multigrid
using
number
of
stiffness
alternate method
should
for or
Obviously,
the
solution
techniques
original
stiffness
additional
grid
elements, total
the
number
multigrid if all matrix
levels
coarser of
method
entries must
must
be
be
levels
elements
requires the
stored
as
M
Ne/2,
levels
storage
bandwidth
the
accommodated.
have
for
within
more
are
finest
If
the
Ne/4 , ...,
space
grid finest
etc.,
than
stored.
The
level grid
direct
and has
N
elements.
e
Thus,
the
is
i__)
The
total
the
present
are
(N e+1)
N e + Ne/2
+ Ne/4
number
degrees
of
+ Ne/8
of
one-dimensional nodes,
each
required,
assuming
following
expression:
...
freedom
case,
with
+
2
for
degrees
a constant
N
the
+
This which
may
must
be
some
zero
to
result
in
in equation be
_torcd.
1} 2M
assumes the
may
of
freedom.
the
B,
can
be
fashion.
finest
Hence,
the
found
(35)
< 2Ne
grid, total
from
In
there
storage
the
solution
become
How_:ver,
is always storage
of
of
several
÷
(36)
that
storlng
less
than
2B(2Ne÷M).
all
entries
within
entries.
These
zero
techniques
nonzeroes.
shows
= 2B(M
) < 2B(2Ne÷M)
the
storage
(5)
a similar
on
+ ..
required
in direct
entries
N e elements
•
storage
computation
stored
iteration n_d
maximum
in
- 2M
N
2Ne{1_
Thus,
_ 2Ne(1
is computed
seml-bandwidth
N 2B(
+ Ne/(2M-I)
since
during
Examination
of
the the
only
the
nonzero
terms
only
the
nonzero
entrles
2O
the
bandwidth,
zero
entries
solution
phase
Gauss-Seidel of
the
matrices
requlrcs
the:
additional
storage of pointer arrays for each row of the stiffness
Many finite-element
stiffness
matrix.
matrices are sparsely populated and, in such
cases, the storage of only the nonzero terms, with the pointer arrays,
can
considerably reduce the storage requirements comparedwith banded or profile storage schemes. With the definition
of the coarser grids used here, the total
required by the stiffness
matrices of all
that needed by the finest
level.
grid must be stored,
the restriction
levels should be less than twice
Additionally,
the solution
as well as the load vectors,
Depending upon the structure
storage
vectors for each
or right-hand sides.
of the problem, it may be advantageous tostore
and prolongation
The work (number of floating
matrices as well. point operations)
done in a single V(nl,n 2)
cycle is equivalent to 2(ni+ n2) times the work of a single relaxation finest
grid level.
n1+ n2 relaxations relaxations
This is determined as follows. are done on the finest
done on the lower levels
on the
In the V(nl,n 2) cycle,
level while the total
sum of the
is always less than or equal to n_+ n2.
Because each subsequently coarser grid level has half as many degrees of freedom as the one before, the total levels (36)).
number of degrees of freedom in all
is always less than twice the number in the finest The work done iterating
on each level
is directly
number of degrees of freedom that must be relaxed.
grid (see equation related
to the
Thus, the total
work done
in a single V(nl,n 2) cycle is less than the work of 2(n1+n2) relaxations the finest
on
grid level. RESULTS ANDDISCUSSION
In this section the performance of the multigrid
algorithm
simply supported beamproblem of Figure I is evaluated.
21
in solving the
Six grid levels are
used, where the finest
mesh (level
the next coarser grid (level on until
I) has 32 elements modeling half
the beam,
2) has 16 elements modeling half the beam, and_so
the coarsest grid (level
6), which has one element modeling half the
beam. The two prolongation operators mentioned earlier simple linear
interpolation
and a prolongation
are considered:
a
operator based on a constant
energy in each grid level. Results are presented for a V(nl,n 2) cycle. solution
and the work necessary to achieve convergence are presented.
Convergenceof the solution residuals,
was monitored with the L2-norm of the
defined in equation (8).
exact solution definition, further
The convergence of the
within
While the displacements converge to the
machine accuracy, the L2-norm, because of its
converges with less accuracy than the displacements.
Appendix C
explains this behavior.
Results are given first
for sequential ordering with linear
prolongation,
for varying values of nI and n2, and then with the energy-conservlng prolongation.
Finally,
the effects
of red-black ordering are presented.
Sequential Orderlng Linear Prolongation Figure 7 presents the L2-norm plotted against an increasing number of Vcycles for both random (solid tions.
These results
solved exactly.
line)
and zero (dashed line)
initial
approxima-
are for a V(2,1) cycle where the coarsest level was
As seen in the figure,
at the end of 20 V-cycles,
is on the order of 10-3 for the random initial
the L2-norm
approximation and on the order
--.6
of
10
solution The
r_ces
for
the
zero
is close of
to
initial the
convergence
approximation,
direct
solution
(indicated
by
indicating obtainable
the
22
slope
on of
the
that
the
multigrid
the
finest
grid
curves)
are
level.
almost
identical
for both the zero and random initial
approximations.
algorithm more severely, only the random initial remainder of this work.
To test the
approximation is used in the
The work necessary to achieve the accuracy shown in
Figure 7 is 2(ni+ n2) times 20 V-cycles, or equivalent to the work of 120 relaxations
on the finest
grid level.
Varying nI and n2 To study the effect
of varying the number of relaxations
the V-cycle, two cases with linear case, the number of relaxations
prolongation were considered.
leg of
In the first
on the downward leg of the cycle was held
constant at nI _ I, while the number of relaxations cycle was increased from I to 5 (n2 _ I to 5). a function of increasing n2.
on either
on the upward leg of the
Figure 8 shows the L2-norm as
As expected, the rate of convergence improves -7
with
increasing
after of
about
216
second
9.
As
increasing
For
the
18 V-cycles.
relaxations
The Figure
n 2.
on
The
the
n2 _
the
previous
numbers
of
cycle,
work
finest
case, in
V(1,5)
of
grid
I while
18
the V(1,5)
relaxations.
is on
cycles
is
increased
the
rate
For
the
of
from
order
of
of
the
cycle
after
previous that
about
case.
to
of
Figures
8 and
9 shows
V(n,1)
cycle
when
following
of
to
10
that
convergence
V(5,1)
cycle,
5,
is
shown
improved
w_th
the
L2-norm
in
is on
216
n
18 V-cycles, As
before,
relaxations that
is
the
the on
V(1,n)
greater
compared
than
the
work
finest
cycle 2.
of
to
10
18
V(1,5)
grid
yields
This
can
for
cycles
level.
better be
the
V(1,5) is
A comparison
convergence
explained
in the
frequency
errors
than
manner.
A prolongation level
of
-7
10
equivalent
the
order
equivalent
I to
-6 the
the
level.
n I was
case,
L2-norm
k-1.
interpolation.
This
from is
Level
level
k to
illustrated k represents
k-1
introduces
schematically a coarse
23
high
in Figure mesh
with
10
for
at
linear
2 elements,
AC
and
CE.
Level
k-1
linear
wB =
represents
prolongation,
(w A
+ Wc)/2
rotations than
at
at
and
E,
as
(osc_llatory)
a
k-1
with
hence
the
better
Figure
a
in
shows the
For
solution a work
the
variation
V(1,n I)
cycle The
n I equal within
For
n I equal
to
for
a work
expenditure
the
finest
grid.
to
the
equivalent
for
plus
total
can
cycle
the
at
high
errors
with
and
Using
found
for
a
as
the
B and
D are
error
can
larger
be
frequency be
does
V(n,1)
DE.
D are
expressions
The
V(1,n)
CD,
B and
displacements
error
The
BC,
at
similar
frequency
compared
k-1
10(b).
(smooth)
relaxations.
prolongation).
converges.
Figure
high
level
with
in the
the
These
on
AB,
easily
eliminated
precisely
this,
on
and
cycle.
Prolongation
for
conserving
4 elements,
+ WE)/2
frequency
few
with
errors
convergence
11
V-cycles
required
The
shown
low
Energy-Conserving
with
= (w C
D.
as
mesh
displacements
wD
errors.
level
exact
fine
the
B and
C and
represented
of
a
2,
3,
with
4 or
accuracy
s'lightiy
convergence,
less
the
of
the
work
as
of
and
5 and the
60
n I equal is needed
of
P
the
_ RT
faster
results
converge
were
on
required
iterations, to
I or
than
the
solution to
the
8 V-cycles finest
for
grid.
convergence
respectively,
2,
when
number
(energy-
the
in approximately
relaxations
V-cycles
48
function
n I,
multigrid
64
a
...,
machine
or
more
for
value
the
80
to
although
L2-norm
n I = I, 2,
5,
96,
equivalent
Thus,
the
larger
to about
I and
of
more
V-cycles
n I is equal
to
on are 3,
4 or
5. The cycle.
energy-conserving Two
cases
are
virtually
the
energy-conserving
_rrors
from
are
identical
the
coarse
prolongation
compared to
in Figure
those
for
prolongation level
(P
to
the
= R T) was
12. V(1,n)
introduces the
fine
24
level.
The
also
results
cycle. little
used for
This or
no
with the
is
high
the
V(n,1)
V(n,1)
cycle
expected
since
frequency
Figure 13 compares the exact solution
to the multigrid
energy-conserving prolongation with the V(2,1) cycle. ee, the errors in w and 8, respectively,
Figure 13 shows ew and
exact
w
0 max
w
and
max
shows
e w and
After
each
e
max
e 0 at
seventh
cycle,
the
maximum
the
end
of
(within
machine
the
the
max
are
V-cycle,
the
maximum
multigrld
values
study
multigrid the
error
both
the
within
for
has
this
fifth,
and
seventh
significantly,
converged
to
Figure
and
the
V-cycles.
after
exact
13
the
solution
and
the
red-black nodes
designated
followed
in
by
The
excellent problem
were
the
At each
black
V-cycle. with
scheme
designated
Thus,
a work
of
is
and
the can
in as
the
level,
solution
with
the to
V(nl,n 2)
and all
all the
cycle.
red For
converged
red-black only
red-black be
the
red
in a V(1,1)
equivalent
performance expected
grid
nodes,
prolongations,
a single
is obtained
ordering
black.
energy-conserving
accuracy
grid.
of
4 relaxations
ordering
explained
by
scheme a
cyclic
method. the
(red
node)
nodes
i-I
be
were
one-dimensional
reduction
solution.
Ordering
even-numbered
first,
convergence
fine
In
the
nodes
relaxed
machine
the
performance all
linear
ordering, on
the
cycle,
were
exact
third,
drops
solution
the
accuracy).
odd-numbered
nodes
of
first,
Red-Black To
as
e(x) - e(x)
exact
Here
for
where ew and ee are calculated
w(x) - w(x) ew
solution
cyclic can
and
eliminated
be
reduction expressed
i+I
(black
from
the
method, in
nodes), system.
terms
the of
displacement the
rotation
displacements
and
the
node
This
can
be
25
and
and
i displacement
done
for
all
the
at
node
rotations and
rotation
red
nodes,
i
at can thus
producing a muchsmaller set of equations. repeatedly until
The sameprocedure can be employed
only one or two nodes remain.
system is then found by reversing
the process.
The solution
for the original
This type of process, which
can be applied only to one-dimensional problems, always produces an exact solution
in one cycle.
For this
ordering is exactly equivalent
problem, the current process with red-black to this cyclic
reduction scheme, and thus
converges in one cycle. CONCLUDING REMARKS Multigrid
procedures were developed for use in one-dimenslonal finite
element analysis.
Several aspects of the multigrid
procedure were studied and
explained in detail. The key elements of the multigrid
method are the restriction
prolongation
operators.
principles.
A general procedure for obtaining
outlined.
These operators were derived based on energy the restriction
This procedure considered the residuals
as the loads on the fine grid and transferred conserving manner, onto the coarse grid. used with any multigrid
and the
finite
from the iterative
solution
these loads, in an energy-
The procedure outlined here can be
element procedure.
Various aspects of the V-cycle multigrid analysis of the deflections
matrix was
algorithm were studied in the
of a one-dimensional, simply supported Bernoulli-
Euler beam. From this study, for six grid levels with 32 elements in the finest
grld, I,
the following
With linear
conclusions were drawn.
prolongation and sequential ordering,
algorithm yielded results equivalent finest
the multigrid
which were of machine accuracy with work
to about 200 standard Gauss-Seidel relaxations
grid.
26
on the
2.
With linear
prolongation
and sequential ordering,
the V(1,n) cycle
with n greater than 2 yielded better convergence rates than the V(n,1) cycle. 3.
Maintaining an energy balance during the inter-grid
transfers
required that the prolongation operator be the transpose of the restriction 4.
operator.
The multigrid
algorithm showed improved convergence rates when energy
was conserved in the inter-grid
transfers.
With the energy-
conserving prolongation and sequential ordering, algorithm yielded results equivalent finest 5.
the multigrid
which were of machine accuracy with a work
to about 50 standard Gauss-Seidel relaxations
on the
grid.
The red-black ordering of the relaxation energy-conservlng prolongation accuracy in a single V(I,1) about 4 relaxations
with either
yielded solutions
cycle.
on the finest
27
the linear
with machine
This is equivalent grid level.
or
to performing
REFERENCES Brandt, A.: Multigrid Techniques: 1984 Guide with Applications to Fluid Dynamics. Monograph. Available as GMD-StudieNo. 85. from GMD-FIT, Postfach 1240, D-5205, St. Augustin I, W. Germany. Jameson, A.: Solution of the Euler equations for two dimensional transonic flow by a multlgrid method. Applied Math. and Comp., v. 13, pp. 327-355, 1983. Southwell, R. V.: Stress calculation systematic relaxation of constraints. Ser. A, 151, pp. 56-95, 1935.
in frameworks by the method of I, II, Proc. Roy. Soo. London
Fedorenko, R. P.: The speed of convergence of an iteratlve process. U.S.S.R. Computational Math. and Math. Phys., v. 4, no. 3, PP. 227235, 1964. Bakhvalov, N. S.: A relaxation method for solving elliptic equations. U.S.S.R. Computational Math. and Math. Phys., v. I, no. 5, pp. 1092-96, 1962. Wachpress, C. B.: Iterative Hall Inc., Englewood Cliffs,
Solution of Elliptic N. J., 1966.
Systems. Prentice-
7.
Brandt, A.: Multi-level adaptive solutions to boundary value problems. Math. Comp.31, ICASEReport 76-27, pp. 333-390, 1977.
8.
Nicolaldes, R. A.: On the _2 convergence of an algorithm for solving finite element equations. Math. Comp., 31, pp. 892-906, 1977.
.
Hackbush, W.: On convergence of a multi-grid iteration applied finite element equations. Report 77-8, Institut f_r Angewandte Mathematik, Universit_t K61n, 197'7.
to
10.
Stuben, K. and Trottenburg, algorithms, model problem Methods. W. Hackbush and Math. 960, Springer-Verlag,
U.: "Multigrid methods: fundamental analysis and applications." Multlgrid U. Trottenburg (eds.), Lecture Notes in pp. 1-176, 1982.
I_.
Hackbush, W.: "Multigrid convergence Hackbush and U. Trottenburg (eds.), Springer-Verlag, pp. 177-219, 1982.
12.
Hackbush, W.: "Parabolic multi-grid methods," Computing Methods in Applied Sciences and Engineering VI, R. Glowinski and J.-L. Lions (eds.), North-Holland, pp. 189-198, 1984.
13.
Hackbush, W. and Trottenburg, U. (eds.): Multigrid Notes in Math. 960, Springer-Verlag, 1982.
la.
Holstein, H. and and differential
Paddon, D. equations.
theory," Multigrld Methods. Lecture Notes in Math. 960,
(eds.): Multigrid Clarendon Press,
28
Methods,
methods 1985.
for
Wo
Lecture
integral
15. McCormick, S. and Trottenburg, U. (eds.): Applied Math. and Comp., Special Issue: Multigrid Methods, v. 13, no. 3 and 4, 1983.
29
APPENDIXA - BEAMELEMENT SHAPEFUNCTIONS For an element with two nodes (four unknowns), the deflected be assumedto be given by a cubic function, w = _I Equation
(AI)
cor_responding
can to
be
+ _2 x written
w. and J
w(x)
÷ a3x
0. J
in (j
2
= 1,2)
= w I N I + 61
N2
i.e.,
+ _4x3
terms
of as
+ w2
shape w can
0 shape
functions
_ x
_ b
N i (i
(AI) z I to
4)
follows:
N3
+ _2
N4
(A2)
where
3x 2 NI = I
2x 3
b2 +
2x 2
b3
N2
x3
I x - -_-- + b-_ 0
N3 =
Here
b
is
The given
in
the
element
element the
text
3x 2
2x 3
b2
b3
x2 N4
