Stable Real Time Deformations
Stable Real-Time Deformations
Matthias M¨uller
Julie Dorsey
Leonard McMillan
Robert Jagnow
Barbara Cutler
ETH Z¨urich
Massachusetts Institute of Technology, Laboratory for Computer Science
typically done offline – that is, computers spend minutes or hours
to arrive at a single answer or a simulation of a few seconds.
Real-time simulation of deformable objects is a younger field.
Abstract
The performance of modern computers and graphics hardware has
made physically-based animation possible in real time. But even
The linear strain measures that are commonly used in real-time an-
with today’s best hardware and most sophisticated techniques [De-
imations of deformable objects yield fast and stable simulations.
bunne et al. 2001; Wu et al. 2001; Zhuang 2000], only a few hun-
However, they are not suitable for large deformations. Recently,
dred elements with small deformations have been simulated in real-
more realistic results have been achieved in computer graphics by
time to date. Since simulating the dynamic behavior of deformable
using Green’s non-linear strain tensor, but the non-linearity makes
objects in real time is an important and challenging task, a great
the simulation more costly and introduces numerical problems.
deal of work has been done in the field and a large variety of tech-
In this paper, we present a new simulation technique that is sta-
niques and methods have been proposed in the last two decades
ble and fast like linear models, but without the disturbing artifacts
[Gibson and Mitrich 1997]. Typical applications for real-time de-
that occur with large deformations. As a precomputation step, a
formable objects include virtual surgery [Debunne et al. 2001; Wu
linear stiffness matrix is computed for the system. At every time
et al. 2001], virtual sculpting, games or any application requiring an
step of the simulation, we compute a tensor field that describes the
interactive virtual environment. A real-time simulator could offer
local rotations of all the vertices in the mesh. This field allows us
artists the option to design and test animations interactively before
to compute the elastic forces in a non-rotated reference frame while
rendering their work offline in higher quality.
using the precomputed stiffness matrix. The method can be applied
An interactive simulation system needs to meet two main re-
to both finite element models and mass-spring systems. Our ap-
quirements. It certainly needs to be fast enough to generate 15 to
proach provides robustness, speed, and a realistic appearance in the
20 frames per second. Speed, however, is not the only requirement.
simulation of large deformations.
Ideally, we want to give the user of the system complete freedom
CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional
of action. Thus, stability and robustness are just as important as the
Graphics and Realism—Animation and Virtual Realiy
frame rate.
With the availability of fast computers, there has been a trend in
Keywords: Physically Based Animation, Finite Elements, Large
real-time animation away from simple models such as mass-spring
Deformations, Elasticity, Stiffness Warping
systems toward the more sophisticated Finite Element approach.
FEM is computationally more expensive, but it is physically more
1
Introduction
accurate, and the object’s deformation behavior can be specified us-
ing a few material properties instead of adjusting a large number of
spring constants. However, because of its computational cost, only
Mathematical and physical modeling of deformable objects has a
the simplest variant of FEM has been used so far – namely tetra-
long history in mechanical engineering and materials science. In
hedral elements with linear shape functions. While not suitable for
those disciplines, the main objective is to model the physical world
engineering analysis, such models are sufficient to obtain visually
as accurately as possible. In graphics applications, the primary con-
plausible results.
cern is usually the computational efficiency of generating plausible
behaviors, rather than the accurate prediction of exact results. The
There is an additional option when choosing an FEM model –
most widely used technique to model deformable objects is to view
namely how strain is measured with respect to the deformation of
material as a continuum. In this case, the constitutive laws yield
an object. Linear elasticity only models small deformations accu-
partial differential equations that describe the static and dynamic
rately, but its computational cost is much lower than the cost of
behavior of the material. These equations are usually solved nu-
a non-linear strain measure. One important feature of the linear
merically using the Finite Element Method (FEM) [Bathe 1982] or
approach is that the stiffness matrix of the system is constant and
finite differences [Terzopoulos et al. 1987]. Such simulations are
numerically well-conditioned, yielding a fast and stable simulation.
Under large rotational deformation, however, objects increase un-
¡
e-mail: muellerm@inf.ethz.ch
naturally in volume because the linear model is only a first order
approximation at the undeformed state (see Fig. 7).
Non-linear elasticity, on the other hand, models large rotational
deformations accurately [Picinbono et al. 2000]. With a non-linear
strain measure, the stiffness matrix is no longer constant. For im-
plicit integration it must be reevaluated at every time step as the
Jacobian of the non-linear function that describes the internal elas-
tic forces. This process slows down the simulation and introduces
numerical instabilities when the Jacobian is evaluated far from the
equilibrium state. This is why models have usually only been sub-
jected to small displacements in demonstrations of real-time sys-
tems thus far. Dramatic deformations are not possible without ei-
ther slowing the simulator down or risking numerical divergence.
because they are simple to implement. However, models that treat
In this paper, we propose a new technique that is as fast and sta-
objects as a continuum have several advantages over simple mass-
ble as linear elasticity while avoiding the artifacts associated with
spring networks. The physical material properties can be described
large deformations. We do this by warping the stiffness matrix of
using a few parameters, which can be looked up in textbooks, and
the system according to a tensor field that describes local rotations
the force coupling between mass elements is defined throughout
of the deformed material. In this way, we can use a precomputed
the volume rather than according to the spring network. As a result,
stiffness matrix. The evaluation of the tensor field is much cheaper
continuous models yield more accurate results. The deformation of
than the cost of a single time step. Our technique is easy to un-
an object in such a model is described by a boundary value partial
derstand and implement, making it practical for a wide variety of
differential equation. For realistic objects, this equation cannot be
applications.
solved analytically. A standard technique to solve it numerically is
the Finite Element Method [Bathe 1982]. Using FEM, an object is
1.1
Related Work
subdivided into elements of finite size – typically polyhedra – and
a continuous deformation field within each element is interpolated
Many methods have been proposed to simulate deformable objects
from the deformation vectors at the vertices. Once the interpolation
in real time. We will discuss just a few recent publications and
functions for all the elements are chosen, the deformation vectors at
papers that describe those techniques similar to ours.
all the vertices describe a piecewise continuous deformation field.
To improve the numerical stability of the simulation, Terzopou-
This field incorporated into the partial differential equation yields a
los et al. [Terzopoulos and Witkin 1988] proposed a hybrid model
set of simultaneous algebraic equations for the deformation vectors
that breaks a deformable object into a rigid and a deformable com-
at the vertices.
ponent. The rigid reference body captures the rigid-body motion
Regardless of the choice of element type and shape functions, the
¢
while a discretized displacement function gives the location of
Finite Element Method yields an algebraic function
that relates
mesh nodes relative to their position within the rigid body. As in
the deformed positions of all the nodes in the object to the internal
their approach, we handle the rotational component of the defor-
elastic forces at all the nodes:
mation separately. However, they use one single rotation matrix for
£¥¤§¦©¨
¢ !#"$!&%('0)
the entire model – namely the one associated with the underlying
(1)
rigid body frame – even if regions of the deformable object undergo
£¥¤§¦©¨1
£32
£04
£6
where
)
)555)
'¥7
contains the internal force vec-
large rotations while other regions don’t rotate at all.
2
4
6
!
§!
)3!
)555)3!
'¥7
!&%
In ArtDefo (Accurate Real Time Deformable Objects) [James
tors of all
nodes, and
and
8
2
4
6
!&%
)3!&%
)555)3!&%
'¥7
and Pai 1999], James et al. used linear elasticity in connection with
represent their deformed and original posi-
the Boundary Element Method (BEM) to deform objects in real
tions, respectively.
6FE
This ev6en holds for mass-spring networks. The
AGC
time. Because of the linearity of the model, many system responses
function ¢@9BADC
is, in general, non-linear and encapsu-
can be precomputed and then combined later in real time. How-
lates the material properties as well as the type of mesh and dis-
ever, the linear model is not accurate for large deformations, as we
cretization used.
already mentioned.
Desbrun et al. [Desbrun et al. 1999] split the forces in mass-
2.1
Dynamic Deformation
spring networks into linear and non-linear (rotational) parts. The
!
rotational part is first neglected to compute a rapid approximation
In a dynamic system, the coordinate vector
is a function of time,
!BH'
of the implicit integration. Then they correct the estimate to pre-
. The dynamic equilibrium equation has the following form:
serve momentum.
Y£
IQP
!SRUTWV
!SRX¢ §! "$!&%('
)
To guarantee a real-time frame rate, Debunne et al. [Debunne
ext
(2)
et al. 2001] use an automatic space and time adaptive level of detail
V
P
where ! and ! are the first and second derivatives of ! with respect
technique. The body is partitioned in a multi-resolution hierarchy
I
T
to time,
is the mass matrix and
the damping matrix [Cook
of tetrahedral meshes. The high resolution meshes are only used in
1981]. Eqn. (2) defines a coupled system of
ordinary differential
regions of high stress. This reduces the number of active elements,
`a8
equations for the
position vectors contained in ! . To solve them,
thus increasing the speed of the simulation. We also use this method
8
!BbHc'
the continuous
-dimensional
2
%
function
is approximated by
in our system to further increase the speed of our simulation.
`a8
!
)3!
)5553!
)555
!
!Bbdfe
a series of vectors
, where
approximates
Wu’s approach [Wu et al. 2001] is very similar to Debunne’s
g
H'
. In a first step, (2) is transformed into a system of
technique. They use progressive meshes to adapt the number of
hpiq`a8
equations of first derivatives:
elements according to the internal stresses.
sr
!
V
(3)
1.2
Overview
r
r
£
V
I
"1T
"p¢S§!#"$!&%t'&R
)
ext
r
In the next section, we introduce linear and non-linear models of
where
is an additional vector of
velocities. Although there are
`a8
static and dynamic deformation and discuss their advantages and
mathematically more accurate integration methods (see [Pozrikidis
disadvantages for real-time simulation. This motivates the need for
1998]), Euler’s first order method is known to better handle dis-
our technique called Stiffness Warping, which we describe in Sec-
continuities (caused, for instance, by collisions) than higher order
tion 3. We propose two ways of computing the rotation field of
methods [Desbrun et al. 1999]. The implicit form of Euler’s method
a deformed mesh along which the stiffness matrix is warped. A
approximates (3) as follows:
comparison of our technique with linear and non-linear approaches
shows the advantages of the method. In the last section, we present
uB2
vuw2
g
r
!
!
R
H
a collection of our results.
uw2
uw2
vuw2
vuw2
g
r
r
r
£
I
I
R
Hx"1T
"p¢ §!
"y!&%t'R
'x5
ext
(4)
2
Modeling Deformation
g
bdRU'
H
In order to find the positions and velocities at time
, a
coupled system of algebraic
uw2
equations
vuw2
needs to be solved, because
r
!
There are a variety of ways to model the behavior of deformable
the unknown values
and
appear on both sides of Eu-
objects. Mass-spring networks are popular in real-time simulators
ler’s implicit equation. To compute positions and velocities at time
dynamic (2) simulations, a non-linear algebraic system of equations
has to be solved. This generally involves the computation of the
Jacobian
of ¢ . Since ¢
is
-dimensional,
is a matrix of di-
`a8
mension
. Even though
is usually sparse, its evaluation is
`a8i`a8
computationally expensive. Moreover, the numerical conditioning
of
deteriorates when evaluated far from the equilibrium state.
2.3
Linear Elasticity
Figure 1: Quadratic stress approximates the real stress-deformation
¢
curve better than linear stress. It is not a full second order approxi-
In linear elasticity,
is replaced by a first order approximation:
mation of the real curve though.
4
y
¢S§!#"$!&%('
e§!#"$!&%t'&RX x
!#"$!&%(
'x)
(6)
¢
!
¢
!&%
where
is the Jacobian
of
evaluated at
, usually
called the stiffness matrix of the system. The stiffness matrix is
computed only once before the simulation is run. At every time
step, a linear system (usually well conditioned) has to be solved.
This is why a linear simulation is faster and more stable than a
simulation based on non-linear elasticity. The drawback of this ap-
proach, however, is that large deformations are not rendered cor-
rectly.
More precisely, linear elastic forces are invariant under
Figure 2: With a method to prevent material from over-stretching,
translations but not under rotations. This raises the question of
the linear model can fit a real stress curve appropriately.
whether it is possible to work with a constant linear stiffness matrix
and extract the rotational part of the deformation. The next sec-
tion describes our new technique called Stiffness Warping, which is
g
bdwRU'
H
we use implicit integration because it is stable for much
based on this idea.
larger time steps than explicit integration [Witkin and Baraff 1998],
g
d
H
which only uses quantities at time
. (For a detailed discussion
3
Stiffness Warping
of implicit and explicit methods see [Witkin and Baraff 1997].)
In linear elasticity, the elastic forces for a single tetrahedral element
2.2
Non-Linear Elasticity
in 3D are evaluated as follows:
£¥¤¥¦©¨cy
e§!#"$!&%t'x)
In order for a strain measure to be accurate for large deformations,
(7)
it should not include the rigid body motions of the simulation el-
2§4ae&2¥4
d
£
¤§¦©¨
A
)3!
ements. This can be achieved by defining strain as the change in
where
2¥4
is the element’s stiffness matrix and
!&%
A
length of an infinitesimal material vector going from the original
and
contain the elastic forces, the displaced positions
configuration to the deformed configuration. In 3D, it is more prac-
and the original positions of the four vertices of the tetrahedron. As
!
tical to measure the change of the squared length of a vector, be-
long as the deformed shape
is only stretched and translated with
!%
cause the squared length is merely the dot product of a vector with
respect to the original shape
, the linear approach yields plausible
itself. This is why the Green-Lagrange strain tensor [Bathe 1982]
results. If the transformation from !&% to ! contains a rotation, the
is defined via the expression
artifacts associated with a linear model emerge.
Let us assume now that we know a global rotational component
4
4
fhg
b
t'
"qb
%t'
e
of the rigid body transformation of the element, where fig
)
4
(5)
ADC
C
b
%('
is a 3D
2§4ae&2¥4
(orthogonal) rotation matrix. We can then construct
h
¤j
f
A
, which contains four copies of fhg along its diagonal
where
% and
are corresponding infinitesimal vectors in the
and zeros everywhere else:
undeformed and deformed configuration respectively. The omis-
fig
k
sion of the square root when measuring the length of a vector
fig
yields quadratic strain-displacement and stress-displacement rela-
¤
f
(8)
fig
tionships. This is a nice side effect because for some materials,
k
fig
quadratic stress approximates the real displacement-stress curve
better than linear stress (Fig. 1). Green-Lagrange strain is not a
This matrix rotates quantities of all four nodes of the tetrahedron
full second order approximation of the real stress curve though, be-
fig
by the same matrix
. If we compute the elastic forces as
cause, as with the linear model, there is only one coefficient (i.e.
2
Young’s modulus
) to fit the curve.
£0¤§¦©¨cbv
¤l
f
efDm
! "$!&%('x)
¤
(9)
Desbrun [1999] describes a method to prevent material from
over-stretching. He approximates the real stress curve with a piece-
we get the same forces as if fhg was not present in ! (Fig. 3). We
!
wise linear function (see Fig. 2). When combined with a linear
first rotate the deformed positions
2
back to their original coordi-
¤
stress measure, this method yields realistic results. Thus, the reason
nate frame using the inverse f
m
. 2The forces are then computed in
enf
!p"o!&%t'
¤
m
why one would use a quadratic strain tensor in computer graphics
this coordinate frame as
and then rotated back
¤
f
and real-time simulations is not because a linear deformation-stress
using
.
©q
psct
relationship would not yield plausible results, but because linear
Let
be the
sub-matrix of
containing entries
,
p
`Dir`
dB"
d
"
strain tensors are not invariant under rigid body transformations,
with
and
. Using (9), we get for
`
hvu$wGuY`
`ax
hvu$yzuY`ax
£
d
and therefore are inappropriate for rendering rotational deforma-
the force
at vertex :
6
tions correctly.
2
¢
With a quadratic stain tensor, the function
describing the inter-
£§B
|q
q
q
fig
bf#m
!
"$!&%
'0)
g
(10)
p
nal elastic forces becomes non-linear. Thus, in both static (1) and
q3{B2
located at vertices immediately adjacent to vertex d . Therefore, the
f
rotation matrix
is only used in the local neighborhood of vertex
d
. In this way, the force at vertex d is computed exactly as in linear
FEM, b2ut as if the local neighborhood of vertex d were rotated back
f
by
m
.
q
£
f
We also tried to use the individual
’s to compute
6
2
£w
©q
q
q
f
bf#m
!
"y!&%
'x)
q
(13)
p
q3{B2
but observed that instability may emerge when more than one rota-
£
d
tion matrix is involved in the computation of
and that the stability
f
depends on the method used to compute
.
}~
Computing the elastic forces as in Eqn. (12) yields fast and ro-
bust simulations. However, the forces are not guaranteed to yield
Figure 3: If the rotational part fhg of the deformation ! is known,
2
zero total momentum as elastic forces should. Errors come from the
f
!
the forces can be computed with respect to a deformation
m
g
fact that the same rotation matrix is used in a finite size environment
that only contains translation and stretching. Here, the original el-
and also depend on the way the rotation matrices are computed (see
ement (a) is deformed (b), and then rotated back into the original
next section). Even though the errors in the force vectors at individ-
coordinate frame (c).
ual vertices are tiny and don’t show as long as objects are anchored,
their sum – if non-zero – acts as a ghost force on free floating ob-
jects. In [Desbrun et al. 1999] Desbrun shows how to solve this
problem by performing a simple and computationally cheap cor-
rection step after every time step. We used the same technique in
our simulator.
3.1
Rotation Tensor Field
We now have to answer the question of how to estimate the local
rotations of a deformed mesh. Extracting the rotational part of a
©t
¡¢a£
mapping between two arbitrary sets of vectors is not straightfor-
ward and not unique if the two sets are not related via a pure 3D
rotation. One approach to finding an optimal rotation matrix is to
Figure 4: Instead of using a single rotation matrix fS¤ from an un-
f
minimize an error function using a least squares method. This, how-
derlying rigid body frame (a), we compute local matrices
for
ever, requires the ability to take derivatives with respect to a matrix.
every vertex (b).
Lasenby et al.[1998] describe an elegant alternative that uses ge-
ometric algebra [Hestenes and Sobczyk 1984]. In the geometric
algebra notation, rotations can be represented by multivectors (ro-
q
q
d
0f555l¥¦'
!
!&%
where
and
and
are the displaced and original
tors). Given two sets of vectors, the theory allows for minimizing
positions of vertex . If we use the same approach for the entire
x
with respect to such rotors and for finding optimal rotations.
mesh, we get the following formula for the elastic force at vertex d :
§«
r&§«
For two given sets
e
of vectors
and
with cardinality
©(ª
©
¬
ADC
C
6
a matrix ¢
is formed:
2
£§B
|q
q
q
f
bf#m
!
"$!&%
'0)
¤
¤
(11)
p
q3{B2
©q¯®
q
r
¢
±
e
'0§±
e
'x5
°
°
(14)
ª
{B2
|q
¨
d
0§555
'
where
°
6ne
6
. The
’s are now sub-matrices of
8
p
ADC
C
, the stiffness matrix of the entire mesh. This raises the
2
4
where the vectors ±
)3±
and ±
are orthonormal basis vectors of
C
f
question of what
¤
– the mesh’s rotation – should be in this case.
ADC
¢
. In a second step,
is decomposed by SVD (singular value
If we kept track of a global rigid body frame associated with the de-
³²µ´·¶
¢
7
decomposition [Golub and Loan 1996]), which yields
.
formable body as in Terzopoulos’ model [Terzopoulos and Witkin
The rotation matrix f
is then simply given by the product
f
1988], we could derive
¤
from this rigid body rotation. For stiff
materials with little deformation but arbitrary rigid body motion,
7
¸¶²
f
5
(15)
this model would yield acceptable results.
Large deformations
other than the rigid body modes would still yield the typical arti-
For our simulator, we have also used a simpler and faster tech-
facts of a linear model, such as growth in volume.
nique to compute local rotations. We found that the stability of Eqn.
A natural extension of the rigid body approach is to use indi-
(12) is not sensitive to the choice of the rotation field and that even a
f
d
vidual rotation matrices
for every vertex
in the mesh (Fig. 4).
very simple approach can yield stable and fast simulations. Figure
Hence, instead of rotating the stiffness matrix
, we warp it along
5 illustrates our faster approximation procedure.
£
f
)d
0§555
'
a rotation field described by the matrices
. For
,
8
For corresponding vertices in the undeformed and deformed
we now get:
2
4
§±
)3±
)3±
'
mesh, we compute orthonormal frames of vectors
and
C
±¹
)3±¹
)3±¹
'
¹
2
4
6
based on a selection of outgoing edges
and
, re-
C
º
º
2
2
±
£§B
©q
q
q
f
bf#m
!
"$!&%
'05
spectively. More specifically,
is computed as the normalized
(12)
p
q3{B2
average of three deterministically chosen edges. The second vector
4
2
±
is evaluated as the cross product of ±
and the direction of a
|q
2
4
±
±
±
The only nonzero
in (12) are those for which there is an edge
chosen edge. The last vector
is the cross product of
and
.
C
p
£
2
4
bda)
'
§±
)3±
)3±
'
in the mesh. Thus, the quantities used to compute
are all
These three vectors form a matrix
. The same
C
x
¬
â$ã
ÎÏ
ä
0.02
Linear
ÐGÑ
ÇtÈ
Ò
0.018
ÁÂ
Warped
ÙGÚ
Û
Non-linear
0.016
)
3
ÓGÔ
ËÍÌ
Õ
ßrà
ÜrÝ
ÃcÄ
0.014
á
Þ
me (m
Å(Æ
ÉÊ
ÖG×
Ø
Volu 0.012
¿½À
åæ·çvèéèéèwê
»½¼S¾¾
ëUìYí
îpï
ð·ñvòrï
ó©ò$ï
óòrï
õ÷öùø
ô
0.01
0.008
Figure 5: A fast way of estimating the rotational part of the defor-
1
11
21
31
41
51
mation at a node is to compute the relative rotation between two
Time (50 ms)
orthonormal vector frames that are based on the directions of adja-
cent edges.
Figure 6: The volume of a bar that deformes under gravity simu-
lated using linear, warped and non-linear stress measures.
procedure applied to the deformed mesh yields a matrix
¹
. It is
¬
¹
important that
is computed using the exact same edges, but in
¬
4
Results
¹
their deformed directions
. The rotation matrix we are looking
º
for can now be evaluated as
4.1
The Bars
¹
7
f
5
(16)
¬
¬
To demonstrate the advantages of our approach, we compare it to
In the case of a rigid body, where the two meshes are related
a linear and a non-linear model. In all three cases, we use implicit
via rotations and translations only, this simple approach yields the
Euler integration and lumped inertia and damping matrices. A Con-
correct constant rotation matrix for all the vertices.
jugate Gradients solver [Pozrikidis 1998] is used for Eqn. 4.
¥
¥
We animate a rectangular bar of
vertices containing
i
i
¥
3k
tetrahedral elements. The block is fixed to a wall on one side
3.2
The Algorithm
and deforms under the influence of gravity (Fig. 7). In the linear
Let us summarize the entire simulation algorithm:
case, the stiffness matrix of the object is evaluated once and used
6ne
6
throughout the simulation to compute the internal elastic forces. In
dûýü(þ
{
d
C
C
A
ú
(
)
ÿ
ÿ
¡
the warped stiffness case, we use the same constant stiffness matrix
ücÿ
%
%
and warp it along a rotation field. This field is computed as shown
û
r
û
£¢
!
!&%
d
¥f555
'
ú
;
for all
8
in Fig. 5. In the non-linear case, a new stiffness matrix is computed
at every time step as the Jacobian of the non-linear force function
û
H
k
ú
¢
based on
4
Green’s strain tensor. We use an elastic modulus
of
3k
k5
and a Poisson ratio of
, meaning the volume of the
ú
loop
¬
``
material should not change substantially during the simulation.
– evaluate f
for all d
0§555
'
Fig. 6 depicts the volume of the block versus time. The linear
8
– solve
model shows the typical growth artifact under deformation. As with
6
buw2
buw2
2
r
Yr
r
|q
R
¥¤
"
"Df
bf
§!
R
the non-linear model, our method does not exhibit this problem.
m
q
q3{B2
p
¦¨§©
buw2
buw2
The time to compute one time step is 5 ms in the linear case, 6
g
r
q
£
H
'&"$!&%
'R
¤
q
g
ms for stiffness warping and 12 ms for the non-linear simulation.
buw2
r
for all unknown
0§555
'
, d
The simulated time step is 10 ms. The experiment shows that our
8
uw2
uw2
g
û
r
!
!
R
H
d
0§555
'
approach is nearly as fast as the linear model but as accurate as the
– set
for all
8
non-linear model in terms of volume conservation (try our applet at
û
– H
H·R
graphics.lcs.mit.edu/simulation/warp/).
To demonstrate the stability of stiffness warping, we repeated
ú
end loop
¥
¥
!
the simulation with a longer bar of
vertices and
i
i
"!#
6
6
E
C
C
The function ¢ §!B'X9A
A
describes internal elastic
tetrahedra (Fig. 7). The linear and warped stiffness methods still
!
forces given the deformed coordinates
of all
vertices of a mesh.
yield stable simulations with a time step of 10 ms while the non-
8
This function does not necessarily need to stem from a Finite El-
linear techniques diverges even with bars slightly longer than in the
ement discretization – it can also be defined by a spring network.
previous example.
¢
First, the Jacobian
of
is evaluated. Then the positions and ve-
locities of all the vertices are initialized and the time is set to zero.
4.2
A Simple Tube
In the simulation loop, the
rotation tensor field is evaluated based
!
on the actual coordinates
as described in section
The tube depicted in Fig. 8 is composed of a thousand tetrahedra
uw2
3.1. Then, the
r
linear system for the unknown new velocities
is solved. This
and it is 50 cm x 13 cm in size. For its material we chose a density
%
C
system is derived by substituting Eqn.12 into Eqn.4 for implicit in-
of 1
and a Poisson Ratio of 0.33. The user interacts with
$
©q
tegration. Note that the
are
sub-matrices of
containing
the system by grabbing the tube at one vertex. This vertex is then
p
`vi#`
psct
df"
d
"
entries
with
and
. Note
attached to the mouse via a spring. In the first experiment, only the
`
h#uYw$u½`
`ax
h#uYy
u½`ax
I
also that we lump the mass matrix
in Eqn.4 to the vertices, i.e.
upper part of the tube is included in the simulation while the lower
replace it by its diagonal entries
. The positions of the vertices
part remains fixed to the ground plane. When the entire model is
are then updated using the new velocities before going to the next
animated, the user can pick it up. It bends due to inertial forces.
time step.
The tube shows deformations and vibrations without the artifacts of
a linear model. When dropped from a diagonally oriented position
References
50 cm above the ground, the impact causes deformations that can
&
lead to instabilities in the simulation. The following table shows
BATHE, K. J. 1982. Finite Element Procedures in Engineering
the largest time step we were able to use before the system became
Analysis. Prentice-Hall, New Jersey.
instable. This value depends on the stiffness (Young’s Modulus
)
of the material.
COOK, R. D. 1981. Concepts and Applications of Finite Element
Analysis. John Wiley & Sons, NY.
4
[3k'
]
2.0
1.0
0.5
0.2
0.1
¬
D
Warp
30 ms
20 ms
10 ms
10 ms
10 ms
EBUNNE, G., DESBRUN, M., CANI, M. P., AND BARR, A. H.
2001. Dynamic real-time deformations using space & time adap-
Non-Linear
5 ms
5 ms
2 ms
1 ms
1 ms
tive sampling. In Computer Graphics Proceedings, Annual Con-
As the results show, the simulation using the warped stiffness
ference Series, ACM SIGGRAPH 2001, 31–36.
technique can be further accelerated by choosing larger time steps
than in the non-linear case. Smaller elastic moduli cause larger
DESBRUN, M., SCHR ¨
ODER, P., AND BARR, A. 1999. Interactive
deformations after the impact and smaller time steps need to be
animation of structured deformable objects. Graphics Interface,
taken.
1–8.
GIBSON, S. F., AND MITRICH, B. 1997. A survey of deformable
4.3
The Bunny
models in computer graphics. Technical Report TR-97-19, Mit-
subishi Electric Research Laboratories, Cambridge, MA.
To generate the animation depicted in Fig. 9, we used a volumetric
mesh of 3kkk tetrahedra. The mesh is composed of a bone core and
GOLUB, G. H., AND LOAN, C. F. V. 1996. Matrix Computations,
a layer of skin tetrahedra. Only the bunny’s head, composed of h!"!#
Third Edition. The Johns Hopkins Univ. Pr., Baltimore and Lon-
bone tetrahedra and ! skin tetrahedra is animated. We treat all
don.
bone tetrahedra as one rigid body. This rigid skull can rotate about
a fixed axis and is attached to the mouse via a spring. The skin
HESTENES, D., AND SOBCZYK, G. 1984. Cliffold Algebra to
tetrahedra follow the movement of the skull dynamically.
Geometric Calculus: A unified language for mathematics and
We use the deformation field of the tetrahedral mesh to animate
physics. D. Reidel, Dordrecht.
a triangle surface mesh with higher resolution (3kkk triangles). Ev-
ery vertex in the surface mesh is associated with a tetrahedron in the
JAMES, D., AND PAI, D. K. 1999. Artdefo, accurate real time
volumetric mesh and uses its barycentric coordinate with respect to
deformable objects. In Computer Graphics Proceedings, Annual
that tetrahedron to interpolate its position.
Conference Series, ACM SIGGRAPH 99, 65–72.
LASENBY, J., FITZGERALD, W. J., DORAN, C. J. L., AND
4.4
The Great Dane
LASENBY, A. N. 1998. New geometric methods for computer
vision. Int. J. Comp. Vision 36(3), 191–213.
As our last example, we animate the floppy skin of a Great Dane
(Fig. 10). As in the bunny example, we simulate the bone core as
PICINBONO, G., DELINGETTE, H., AND AYACHE, N.
2000.
a rigid body and let the skin layer follow its movements, but in this
Real-time large displacement elasticity for surgery simulation:
case, the entire model (i.e.
bone and
¥¥
skin tetrahedra) is
Non-linear tensor-mass model. Third International Conference
"
`
h
animated. The
4
elastic modulus
of the skin in the Dane’s face
on Medical Robotics, Imaging And Computer Assisted Surgery:
3k¦C
is
– much lower than in the previous examples – which
MICCAI 2000 (Oct.), 643–652.
¬
makes the surface lag noticeably behind the skull movement. The
3kkk
visible surface mesh is formed with
triangles, the vertices of
POZRIKIDIS, C. 1998. Numerical Computation in Science and
which are interpolated using the underlying tetrahedral mesh.
Engineering. Oxford Univ. Press, NY.
TERZOPOULOS, D., AND WITKIN, A. 1988. Physically based
5
Conclusions
models with rigid and deformable components. IEEE Computer
Graphics & Applications (Nov.), 41–51.
In this paper, we have presented a new technique to animate de-
formable objects in real-time. By warping the constant stiffness
TERZOPOULOS, D., PLATT, J., BARR, A., AND FLEISCHER, K.
matrix of the system used in linear approaches along a rotation field,
1987. Elastically deformable models. In Computer Graphics
we eliminate the visual artifacts while the simulator remains as sta-
Proceedings, Annual Conference Series, ACM SIGGRAPH 87,
ble and fast as a linear one, even for large rotational deformations.
205–214.
In contrast to a non-linear approach, the stiffness matrix needs only
WITKIN, A., AND BARAFF, D. 1997. Physically based modeling:
to be computed once. The same matrix can be used throughout the
Principles and practice. Siggraph Course Notes (Aug.).
entire simulation for implicit integration, making the system fast
and robust. We have also proposed a fast way of estimating a rota-
WITKIN, A., AND BARAFF, D. 1998. Large steps in cloth simu-
tion field along which the stiffness matrix is warped at every time
lation. In Computer Graphics Proceedings, Annual Conference
step.
Series, ACM SIGGRAPH 98, 43–54.
Our examples show that stiffness warping makes possible real-
time animation of detailed models in an interactive environment.
WU, X., DOWNES, M. S., GOKTEKIN, T., AND TENDICK, F.
In the future, we would like to incorporate material fracture into
2001. Adaptive nonlinear finite elements for deformable body
our simulator. Stiffness warping works with a constant system ma-
simulation using dynamic progressive meshes.
Eurographics
trix. This matrix changes when the structure of the underlying mesh
(Sept.), 349–358.
changes. Fortunately, local changes in the mesh only cause local
changes in the coefficients of the global stiffness matrix. Such up-
ZHUANG, Y. 2000. Real-time Simulation of Physically Realistic
dates can be done incrementally and will not slow down the simu-
Global Deformation. Ph. D. thesis of Univ. of California, CA.
lation significantly.
Stable Real-Time Deformations: M. M¨uller, J. Dorsey, L. McMillan, R. Jagnow, B. Cutler
Figure 7: Three bars attached to a wall under the influence of gravity. They are simulated using non-linear (green), warped (blue) and linear
(red) strain measures. Longer bars more noticeably show the artifacts with linear FEM.
(a)
(b)
(c)
(d)
Figure 8: A tube is bent under user-applied forces (a), inertial forces (b) and collision forces with low (c) and high (d) elasticity modulus.
Figure 9: The bone core (white) is animated as a rigid body while the bunny’s skin follows it dynamically.
Figure 10: The great dane’s skin has a low elastic modulus, which makes the surface lag noticeably behind the skull movement.
