3rd Workshop in Virtual Reality Interactions and Physical Simulation "VRIPHYS" (2006)
C. Mendoza, I. Navazo (Editors)

Position Based Dynamics
Matthias Müller Bruno Heidelberger Marcus Hennix John Ratcliff
AGEIA

Abstract
The most popular approaches for the simulation of dynamic systems in computer graphics are force based. Internaland external forces are accumulated from which accelerations are computed based on Newton’s second law of
motion. A time integration method is then used to update the velocities and finally the positions of the object.
A few simulation methods (most rigid body simulators) use impulse based dynamics and directly manipulate
velocities. In this paper we present an approach which omits the velocity layer as well and immediately works
on the positions. The main advantage of a position based approach is its controllability. Overshooting problemsof explicit integration schemes in force based systems can be avoided. In addition, collision constraints can behandled easily and penetrations can be resolved completely by projecting points to valid locations. We have used
the approach to build a real time cloth simulator which is part of a physics software library for games. Thisapplication demonstrates the strengths and benefits of the method.
Categories and Subject Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry
and Object ModelingPhysically Based Modeling; I.3.7 [Computer Graphics]: Three-Dimensional Graphics andRealismAnimation and Virtual Reality

1. Introduction
Research in the field of physically based animation in computer graphics is concerned with finding new methods forthe simulation of physical phenomena such as the dynamics
of rigid bodies, deformable objects or fluid flow. In contrast
to computational sciences where the main focus is on accuracy, the main issues here are stability, robustness and speed
while the results should remain visually plausible. Therefore, existing methods from computational sciences can not
be adopted one to one. In fact, the main justification for
doing research on physically based simulation in computer
graphics is to come up with specialized methods, tailored tothe particular needs in the field. The method we present falls
into this category.
The traditional approach to simulating dynamic objects
has been to work with forces. At the beginning of each time
step, internal and external forces are accumulated. Examples
of internal forces are elastic forces in deformable objects or
viscosity and pressure forces in fluids. Gravity and collision
forces are examples of external forces. Newton’s second lawof motion relates forces to accelerations via the mass. So usc The Eurographics Association 2006.
°

ing the density or lumped masses of vertices, the forces aretransformed into accelerations. Any time integration scheme
can then be used to first compute the velocities from the accelerations and then the positions from the velocities. Someapproaches use impulses instead of forces to control the animation. Because impulses directly change velocities, one
level of integration can be skipped.
In computer graphics and especially in computer games
it is often desirable to have direct control over positions ofobjects or vertices of a mesh. The user might want to attacha vertex to a kinematic object or make sure the vertex always
stays outside a colliding object. The method we propose here
works directly on positions which makes such manipulations easy. In addition, with the position based approach it ispossible to control the integration directly thereby avoiding
overshooting and energy gain problems in connection with
explicit integration. So the main features and advantages of
position based dynamics are• Position based simulation gives control over explicit integration and removes the typical instability problems.

M. Müller et al. / Position Based Dynamics

Figure 1: A known deformation benchmark test, applied here to a cloth character under pressure.

• Positions of vertices and parts of objects can directly be
manipulated during the simulation.
• The formulation we propose allows the handling of general constraints in the position based setting.
• The explicit position based solver is easy to understandand implement.

2. Related WorkThe recent state of the art report [NMK∗ 05] gives a good
overview of the methods used in computer graphics to simulate deformable objects, e.g. mass-spring systems, the finite
element method or finite difference approaches. Apart fromthe citation of [MHTG05], position based dynamics does notappear in this survey. However, parts of the position based
approach have appeared in various papers without naming it
explicitly and without defining a complete framework.
Jakobsen [Jak01] built his Fysix engine on a position
based approach. His central idea was to use a Verlet integrator and manipulate positions directly. Because velocities
are implicitly stored by current and the previous positions,the velocities are implicitly updated by the position manipulation. While he focused mainly on distance constraints,
he only gave vague hints on how more general constraints
could be handled. In this paper we present a fully generalapproach which handles general constraints. We also focus
on the important issue of conservation of linear and angular momenta by position projection. We work with explicit
velocities instead of storing previous positions which makes
damping and friction simulation much easier.
Desbrun [DSB99] and Provot [Pro95] use constraint projection in mass spring systems to prevent springs from overstretching. In contrast to a full position based approach, projection is only used as a polishing process for those springs
that are stretched too much and not as the basic simulation
method.
Bridson et al. use a traditional force based approach for

cloth simulation [BFA02] and combine it with a geometric collision resolving algorithm based on positions to makesure that the collision resolving impulses are kept within stable bounds. The same holds for the kinematical collision correction step proposed by Volino et al. [VCMT95].
A position based approach has been used by Clavet etal. [CBP05] to simulate viscoelastic fluids. Their approachis not fully position based because the time step appears invarious places of their position projections. Thus, the integration is only conditionally stable as regular explicit integration.
Müller et al. [MHTG05] simulate deformable objects
by moving points towards certain goal positions which arefound by matching the rest state to the current state of theobject. Their integration method is the closest to the one we
propose here. They only treat one specialized global constraint and, therefore, do not need a position solver.
Fedor [Fed05] uses Jakobsen’s approach to simulate characters in games. His method is tuned to the particular problem of simulating human characters. He uses several skeletal
representations and keeps them in sync via projections.
Faure [Fau98] uses a Verlet integration scheme by modifying the positions rather than the velocities. New positions
are computed by linearizing the constraints while we workwith the non linear constraint functions directly.
We define general constraints via a constraint function
as [BW98] and [THMG04]. Instead of computing forces as
the derivative of a constraint function energy, we directly
solve for the equilibrium configuration and project positions.
With our method we derive a bending term for cloth which
is similar to the one proposed in [GHDS03] and [BMF03]
but adopted to the point based approach.
In Section 4 we use the position based dynamics approachfor the simulation of cloth. Cloth simulation has been an active research field in computer graphics in recent years. Instead of citing the key papers of the field individually werefer the reader to [NMK∗ 05] for a comprehensive survey.
c The Eurographics Association 2006.
°

M. Müller et al. / Position Based Dynamics

3. Position Based Simulation
In this section we will formulate the general position based
approach. With cloth simulation, we will give a particularapplication of the method in the subsequent and in the results
section. We consider a three dimensional world. However,
the approach works equally well in two dimensions.
3.1. Algorithm Overview
We represent a dynamic object by a set of N vertices and M
constraints. A vertex i ∈ [1, . . . , N] has a mass mi , a position
xi and a velocity vi .
A constraint j ∈ [1, . . . , M] consists of• a cardinality n j ,• a function C j : R3n j → R,
• a set of indices {i1 , . . . in j }, ik ∈ [1, . . . N],• a stiffness parameter k j ∈ [0 . . . 1] and• a type of either equality or inequality.
Constraint j with type equality is satisfied if
C j (xi1 , . . . , xin j ) = 0. If its type is inequality then it issatisfied if C j (xi1 , . . . , xin j ) ≥ 0. The stiffness parameter k j
defines the strength of the constraint in a range from zero toone.
Based on this data and a time step ∆t, the dynamic object
is simulated as follows:in exact correspondence with a Verlet integration step anda modification of the current position [Jak01], because theVerlet method stores the velocity implicitly as the difference
between the current and the last position. However, working
with velocities allows for a more intuitive way of manipulating them.
The velocities are manipulated in line (5), (6) and (16).
Line (5) allows to hook up external forces to the system ifsome of the forces cannot be converted to positional constraints. We only use it to add gravity to the system in whichcase the line becomes vi ← vi + ∆tg, where g is the gravitational acceleration. In line (6), the velocities can be damped
if this is necessary. In Section 3.5 we show how to add global
damping without influencing the rigid body modes of theobject. Finally, in line (16), the velocities of colliding vertices are modified according to friction and restitution coefficients.
The given constraints C1 , . . . ,CM are fixed throughout thesimulation. In addition to these constraints, line (8) generates
the Mcoll collision constraints which change from time stepto time step. The projection step in line (10) considers both,the fixed and the collision constraints.
The scheme is unconditionally stable. This is because theintegration steps (13) and (14) do not extrapolate blindly
into the future as traditional explicit schemes do but movethe vertices to a physically valid configuration pi computedby the constraint solver. The only possible source for instabilities is the solver itself which uses the Newton-Raphson
method to solve for valid positions (see Section 3.3). However, its stability does not depend on the time step size buton the shape of the constraint functions.

(1) forall vertices i(2)initialize xi = x0i , vi = v0i , wi = 1/mi
(3) endfor(4) loop(5)
forall vertices i do vi ← vi + ∆twi fext (xi )(6)
dampVelocities(v1 , . . . , vN )The integration does not fall clearly into the category(7)
forall vertices i do pi ← xi + ∆tviof implicit or explicit schemes. If only one solver iteration(8)
forall vertices i do generateCollisionConstraints(xi → pi ) is performed per time step, it looks more like an explicit(9)
loop solverIterations times
scheme. By increasing the number of iterations, however, a(10)
projectConstraints(C1 , . . . ,CM+Mcoll , p1 , . . . , pN )constrained system can be made arbitrarily stiff and the al(11) endloop
gorithm behaves more like an implicit scheme. Increasing(12) forall vertices ithe number of iterations shifts the bottleneck from collision(13)vi ← (pi − xi )/∆tdetection to the solver.
(14)
xi ← pi(15) endfor
3.2. The Solver(16) velocityUpdate(v1 , . . . , vN )(17) endloop
The input to the solver are the M + M
constraints andcollLines (1)-(3) just initialize the state variables. The coreidea of position based dynamics is shown in lines (7), (9)(11) and (13)-(14). In line (7), estimates pi for new locationsof the vertices are computed using an explicit Euler integration step. The iterative solver (9)-(11) manipulates these
position estimates such that they satisfy the constraints. It
does this by repeatedly project each constraint in a GaussSeidel type fashion (see Section 3.2). In steps (13) and (14),the positions of the vertices are moved to the optimized estimates and the velocities are updated accordingly. This isc The Eurographics Association 2006.
°

the estimates p1 , . . . , pN for the new locations of the points.
The solver tries to modify the estimates such that they satisfy all the constraints. The resulting system of equations
is non-linear. Even a simple distance constraint C(p1 , p2 ) =|p1 − p2 | − d yields a non-linear equation. In addition, theconstraints of type inequality yield inequalities. To solve
such a general set of equations and inequalities, we use a
Gauss-Seidel-type iteration. The original Gauss-Seidel algorithm (GS) can only handle linear system. The part weborrow from GS is the idea of solving each constraint independently one after the other. However, in contrast to GS,

M. Müller et al. / Position Based Dynamics

solving a constraint is a non linear operation. We repeatedly iterate through all the constraints and project the particles to valid locations with respect to the given constraint
alone. In contrast to a Jacobi-type iteration, modifications topoint locations immediately get visible to the process. This
speeds up convergence significantly because pressure waves
can propagate through the material in a single solver step, aneffect which is dependent on the order in which constraints
are solved. In over-constrained situations, the process canlead to oscillations if the order is not kept constant.
3.3. Constraint Projection
Projecting a set of points according to a constraint means
moving the points such that they satisfy the constraint. The
most important issue in connection with moving points directly inside a simulation loop is the conservation of linear
and angular momentum. Let ∆pi be the displacement of vertex i by the projection. Linear momentum is conserved if∑ mi ∆pi = 0,(1)iwhich amounts to conserving the center of mass. Angularmomentum is conserved if∑ ri × mi ∆pi = 0,(2)iwhere the ri are the distances of the pi to an arbitrary common rotation center. If a projection violates one of these constraints, it introduces so called ghost forces which act like external forces dragging and rotation the object. However, onlyinternal constraints need to conserve the momenta. Collision
or attachment constraints are allowed to have global effects
on the object.
The method we propose for constraint projection conserves both momenta for internal constraints. Again, thepoint based approach is more direct in that we can directly
use the constraint function while force based methods derive
forces via an energy term (see [BW98, THMG04]). Let us
look at a constraint with cardinality n on the points p1 , . . . , pn
with constraint function C and stiffness k. We let p be theconcatenation [pT1 , . . . , pTn ]T . For internal constraints, C isindependent of rigid body modes, i.e. translation and rotation. This means that rotating or translating the points does
not change the value of the constraint function. Therefore,
the gradient ∇pC is perpendicular to rigid body modes because it is the direction of maximal change. If the correction∆p is chosen to be along ∇Cp both momenta are automatically conserved if all masses are equal (we handle different
masses later). Given p we want to find a correction ∆p suchthat C(p + ∆p) = 0. This equation can be approximated by
C(p + ∆p) ≈ C(p) + ∇pC(p) · ∆p = 0.

(3)

Restricting ∆p to be in the direction of ∇pC means choosing a scalar λ such that∆p = λ ∇pC(p).

∆p2

d

∆p1

m2
p2

m1
p1
Figure 2: Projection of the constraint C(p1 , p2 ) = |p1 −p2 | − d. The corrections ∆pi are weighted according to theinverse masses wi = 1/mi .

Substituting Eq. (4) into Eq. (3), solving for λ and substituting it back into Eq. (4) yields the final formula for ∆p∆p = −C(p)∇pC(p)|∇pC(p)|2(5)

which is a regular Newton-Raphson step for the iterative solution of the non-linear equation given by a single constraint.
For the correction of an individual point pi we have∆pi = −s ∇pi C(p1 , . . . , pn ),(6)where the scaling factors=C(p1 , . . . , pn )|∇∑ j p j C(p1 , . . . , pn )|2(7)

is the same for all points. If the points have individual
masses, we weight the corrections ∆pi by the inverse masses
wi = 1/mi . In this case a point with infinite mass, i.e. wi = 0,
does not move for example as expected. Now Eq. (4) is replaced by∆pi = λ wi ∇pi C(p) yieldings=C(p1 , . . . , pn )∑ j w j |∇p j C(p1 , . . . , pn )|2(8)

for the scaling factor and for the final correction∆pi = −s wi ∇pi C(p1 , . . . , pn ).

(9)

To give an example, let us consider the distance constraintfunction C(p1 , p2 ) = |p1 − p2 | − d. The derivative with respect to the points are ∇p1 C(p1 , p2 ) = n and ∇p2 C(p1 , p2 ) =−p2−n with n = |pp1 −p| . The scaling factor s is, thus, s =12|p1 −p2 |−dw1 +w2 and the final corrections

w1
p1 − p2(|p1 − p2 | − d)w1 + w2|p1 − p2 |w2
p1 − p2∆p2 = +(|p1 − p2 | − d)w1 + w2|p1 − p2 |∆p1 = −

(10)
(11)

which are the formulas proposed in [Jak01] for the projection of distance constraints (see Figure 2). They pop up as a
special case of the general constraint projection method.

(4)c The Eurographics Association 2006.
°

M. Müller et al. / Position Based Dynamics

We have not considered the type and the stiffness k ofthe constraint so far. Type handling is straight forward. If
the type is equality we always perform a projection. If
the type is inequality, the projection is only performed if
C(p1 , . . . , pn ) < 0. There are several ways of incorporating
the stiffness parameter. The simplest variant is to multiply
the corrections ∆p by k ∈ [0 . . . 1]. However, for multiple
iteration loops of the solver, the effect of k is non-linear.
The remaining error for a single distance constraint after
ns solver iterations is ∆p(1 − k)ns . To get a linear relationship we multiply the corrections not by k directly but by
k0 = 1 − (1 − k)1/ns . With this transformation the error becomes ∆p(1 − k0 )ns = ∆p(1 − k) and, thus, becomes linearly
dependent on k and independent of ns as desired. However,
the resulting material stiffness is still dependent on the time
step of the simulation. Real time environments typically use
fixed time steps in which case this dependency is not problematic.
3.4. Collision Detection and Response
One advantage of the position based approach is how simply
collision response can be realized. In line (8) of the simulation algorithm the Mcoll collision constraints are generated.
While the first M constraints given by the object representation are fixed throughout the simulation, the additional Mcoll
constraints are generated from scratch at each time step. The
number of collision constraints Mcoll varies and depends on
the number of colliding vertices. Both, continuous and static
collisions can be handled. For continuous collision handling,
we test for each vertex i the ray xi → pi . If this ray enters an
object, we compute the entry point qc and the surface normal
nc at this position. An inequality constraint with constraintfunction C(p) = (p − qc ) · nc and stiffness k = 1 is added
to the list of constraints. If the ray xi → pi lies completely
inside an object, continuous collision detection has failed atsome point. In this case we fall back to static collision handling. We compute the surface point qs which is closest topi and the surface normal ns at this position. An inequality
constraint with constraint function C(p) = (p − qs ) · ns andstiffness k = 1 is added to the list of constraints. Collision
constraint generation is done outside of the solver loop. This
makes the simulation much faster. There are certain scenarios, however, where collisions can be missed if the solver
works with a fixed collision constraint set. Fortunately, according to our experience, the artifacts are negligible.
Friction and restitution can be handled by manipulating
the velocities of colliding vertices in step (16) of the algorithm. The velocity of each vertex for which a collision constraint has been generated is dampened perpendicular to thecollision normal and reflected in the direction of the collisionnormal.
The collision handling discussed above is only correct forcollisions with static objects because no impulse is transferred to the collision partners. Correct response for two dyc The Eurographics Association 2006.
°

namic colliding objects can be achieved by simulating bothobjects with our simulator, i.e. the N vertices and M constraints which are the input to our algorithm simply represent
two or more independent objects. Then, if a point q of one
objects moves through a triangle p1 , p2 , p3 of another object,
we insert an inequality constraint with constraint function
C(q, p1 , p2 , p3 ) = ±(q − p1 ) · [(p2 − p1 ) × (p3 − p1 )] whichkeeps the point q on the correct side of the triangle. Since
this constraint function is independent of rigid body modes,
it will correctly conserve linear and angular momentum.
Collision detection gets slightly more involved because thefour vertices are represented by rays xi → pi . Therefore thecollision of a moving point against a moving triangle needs
to be detected (see section about cloth self collision).

3.5. Damping
In line (6) of the simulation algorithm the velocities aredampened before they are used for the prediction of the
new positions. Any form of damping can be used and manymethods for damping have been proposed in the literature(see [NMK∗ 05]). Here we propose a new method with someinteresting properties:
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
(9)

xcm = (∑i xi mi )/(∑i mi )vcm = (∑i vi mi )/(∑i mi )L = ∑i ri × (mi vi )I = ∑i r̃i r̃Ti miω = I−1 L
forall vertices i∆vi = vcm + ω × ri − vivi ← vi + kdamping ∆vi
endforHere ri = xi − xcm , r̃i is the 3 by 3 matrix with the propertyr̃i v = ri × v, and kdamping ∈ [0 . . . 1] is the damping coefficient. In lines (1)-(5) we compute the global linear velocity
xcm and angular velocity ω of the system. Lines (6)-(9) then
only damp the individual deviations ∆vi of the velocities vi
from the global motion vcm + ω × ri . Thus, in the extreme
case kdamping = 1, only the global motion survives and theset of vertices behaves like a rigid body. For arbitrary valuesof kdamping , the velocities are globally dampened but without
influencing the global motion of the vertices.

3.6. Attachments
With the position based approach, attaching vertices to staticor kinematic objects is quite simple. The position of the vertex is simply set to the static target position or updated at every time step to coincide with the position of the kinematic
object. To make sure other constraints containing this vertex
do not move it, its inverse mass wi is set to zero.

M. Müller et al. / Position Based Dynamics
n1

p2

p2

n

n1

n
q

p1,2
ϕ

p4
p3

q

p1

p3
p3

n2

n2

h

p1

p2

p4

p1

p3

Figure 4: For bending resistance, the constraint function
C(p1 , p2 , p3 , p4 ) = arccos(n1 · n2 ) − ϕ0 is used. The actual
dihedral angle ϕ is measure as the angle between the normals of the two triangles.

Figure 5: Constraint function C(q, p1 , p2 , p3 ) = (q − p1 ) ·n − h makes sure that q stays above the triangle p1 , p2 , p3
by the the cloth thickness h.

4.2. Collision with Rigid Bodies
4. Cloth Simulation
We have used the point based dynamics framework to implement a real time cloth simulator for games. In this section
we will discuss cloth specific issues thereby giving concrete
examples of the general concepts introduced in the previous
section.

For collision handling with rigid bodies we proceed as described in Section 3.4. To get two-way interactions, we apply
an impulse mi ∆pi /∆t to the rigid body at the contact point,
each time vertex i is projected due to collision with thatbody. Testing only cloth vertices for collisions is not enoughbecause small rigid bodies can fall through large cloth triangles. Therefore, collisions of the convex corners of rigid
bodies against the cloth triangles are also tested.

4.1. Representation of Cloth
Our cloth simulator accepts as input arbitrary triangle
meshes. The only restriction we impose on the input meshis that it represents a manifold, i.e. each edge is shared by atmost two triangles. Each node of the mesh becomes a simulated vertex. The user provides a density ρ given in mass perarea [kg/m2 ]. The mass of a vertex is set to the sum of one
third of the mass of each adjacent triangle. For each edge,
we generate a stretching constraint with constraint function
Cstretch (p1 , p2 ) = |p1 − p2 | − l0 ,stiffness kstretch and type equality. The scalar ł0 is the initiallength of the edge and kstretch is a global parameter provided
by the user. It defines the stretching stiffness of the cloth. For
each pair of adjacent triangles (p1 , p3 , p2 ) and (p1 , p2 , p4 )we generate a bending constraint with constraint functionµ
acos

Cbend (p1 , p2 , p3 , p4 ) =(p2 − p1 ) × (p3 − p1 ) (p2 − p1 ) × (p4 − p1 )·|(p2 − p1 ) × (p3 − p1 )| |(p2 − p1 ) × (p4 − p1 )|¶− ϕ0 ,

stiffness kbend and type equality. The scalar ϕ0 is the initial dihedral angle between the two triangles and kbend is a
global user parameter defining the bending stiffness of thecloth (see Figure 4). The advantage of this bending term
over adding a distance constraint between points p3 and p4
or over the bending term proposed by [GHDS03] is that it isindependent of stretching. This is because the term is independent of edge lengths. This way, the user can specify clothwith low stretching stiffness but high bending resistance for
instance (see Figure 3).
Eqns. (10) and (11) define the projection for the stretching constraints. In the appendix A we derive the formulas toproject the bending constraints.

4.3. Self Collision
Assuming that the triangles all have about the same size,
we use spatial hashing to find vertex triangle collisions[THM∗ 03]. If a vertex q moves through a triangle p1 , p2 ,p3 , we use the constraint function(p2 − p1 ) × (p3 − p1 )− h,|(p2 − p1 ) × (p3 − p1 )|(12)
where h is the cloth thickness (see Figure 5). If the vertex
enters from below with respect to the triangle normal, theconstraint function has to beC(q, p1 , p2 , p3 ) = (q − p1 ) ·(p3 − p1 ) × (p2 − p1 )−h|(p3 − p1 ) × (p2 − p1 )|(13)
to keep the vertex on the original side. Projecting these constraints conserves linear and angular momentum which is essential for cloth self collision since it is an internal process.
Figure 6 shows a rest state of a piece of cloth with self collisions. Testing continuous collisions is insufficient if clothgets into a tangled state, so methods like the ones proposedby [BWK03] have to be applied.
C(q, p1 , p2 , p3 ) = (q − p1 ) ·4.4. Cloth Balloons
For closed triangle meshes, overpressure inside the mesh caneasily be modeled (see Figure 7). We add an equality constraint concerning all N vertices of the mesh with constraintfunction!
Ãn
trianglesC(p1 , . . . , pN ) =∑ (pt × pt ) · pt

i=1

i
1

i
2

i
3

− kpressureV0(14)

c The Eurographics Association 2006.
°

M. Müller et al. / Position Based Dynamics

Figure 3: With the bending term we propose, bending and stretching are independent parameters. The top row shows
1
1(kstretching , kbending ) = (1, 1), ( 21 , 1) and ( 100, 1). The bottom row shows (kstretching , kbending ) = (1, 0), ( 21 , 0) and ( 100, 0).

Figure 7: Simulation of overpressure inside a character.
Figure 6: This folded configuration demonstrates stable self
collision and response.

and stiffness k = 1 to the set of constraints. Here t1i ,t2i and t3i
are the three indices of the vertices belonging to triangle i.
The sum computes the actual volume of the closed mesh. It
is compared against the original volume V0 times the overpressure factor kpressure . This constraint function yields thegradients∇pi C = ∑ (pt j × pt j ) + ∑ (pt j × pt j ) + ∑ (pt j × pt j )j:t1j =i2

3

j:t2j =i3

1

j:t3j =i

1

2

(15)
These gradients have to be scaled by the scaling factor given
in Eq. (7) and weighted by the masses according to Eq. (9)to get the final projection offsets ∆pi .
5. Results
We have integrated our method into Rocket [Rat04], a gamelike environment for physics simulation. Various experic The Eurographics Association 2006.
°

ments have been carried out to analyze the characteristics
and the performance of the proposed method. All test scenarios presented in this section have been performed on a
PC Pentium 4, 3 GHz.
Independent Bending and Stretching. Our bending term
only depends on the dihedral angle of adjacent triangles, noton edge lengths, so bending and stretching resistances canbe chosen independently. Figure 3 shows a cloth bag withvarious stretching stiffnesses, first with bending resistanceenabled and then disabled. As the top row shows, bending
does not influence stretching resistance.
Attachments with Two Way Interaction. We can simulate both, one way and two way coupled attachment constraints. The cloth stripes in Figure 8 are attached via one
way constraints to the static rigid bodies at the top. In addition, two way interaction is enabled between the stripes andthe bottom rigid bodies. This configuration results in realistically looking swing and twist motions of the stripes. The
scene features 6 rigid bodies and 3 pieces of cloth which are
simulated and rendered with more than 380 fps.

M. Müller et al. / Position Based Dynamics

7. Future Work

Figure 8: Cloth stripes are attached via one way interaction
to static rigid bodies at the top and via two way constraints
to rigid bodies at the bottom.

A topic we have not treated in this paper is rigid body simulation. However, the approach we presented could quite easily be extended to handle rigid objects as well. Instead ofcomputing a set of linear and angular impulses for the resolution of collisions as regular rigid body solvers typicallydo, movements and rotations would be applied to the bodies at the contact points and the linear and angular velocities
would have to be adjusted accordingly after the solver has
completed.

Real Time Self Collision. The piece of cloth shown in
Figure 6 is composed of 1364 vertices and 2562 triangles.
The simulation runs at 30 fps on average including self collision detection, collision handling and rendering. The effectof friction is shown in Figure 9 where the same piece of cloth
is tumbling in a rotating barrel.
Tearing and stability. Figure 10 shows a piece of clothconsisting of 4264 vertices and 8262 triangles that is torn
open by an attached cube and finally ripped apart by athrown ball. This scene is simulated and rendered with 47 fps
on average. Tearing is simulated by a simple process: Whenever the stretching of an edge exceeds a specified thresholdvalue, we select one of the edge’s adjacent vertices. We thenput a split plane through that vertex perpendicular to the edge
direction and split the vertex. All triangles above the split
plane are assigned to the original vertex while all triangles
below are assigned to the duplicate. Our method remains stable even in extreme situations as shown in Figure 1, a scene
inspired by [ITF04]. An inflated character model is squeezed
through rotating gears resulting in multiple constraints, collisions and self collisions acting on single cloth vertices.
Complex Simulation Scenarios. The presented method
is especially suited for complex simulation environments(see Figure 12). Despite the extensive interaction with animated characters and geometrically complex game levels,
simulation and rendering of multiple pieces of cloth can still
be done at interactive speed.
6. Conclusions
We have presented a position based dynamics frameworkthat can handle general constraints formulated via constraintfunctions. With the position based approach it is possible tomanipulate objects directly during the simulation. This significantly simplifies the handling of collisions, attachment
constraints and explicit integration and it makes direct andimmediate control of the animated scene possible.
We have implemented a robust cloth simulator on top of
this framework which provides features like two way interaction of cloth with rigid bodies, cloth self collision and response and attachments of pieces of cloth to dynamic rigid
bodies.
c The Eurographics Association 2006.
°

M. Müller et al. / Position Based Dynamics

Figure 9: Influenced by collision, self collision and friction, a piece of cloth tumbles in a rotating barrel.

Figure 10: A piece of cloth is torn open by an attached cube and ripped apart by a thrown ball.

Figure 11: Three inflated characters experience multiple collisions and self collisions.

Figure 12: Extensive interaction between pieces of cloth and an animated game character (left), a geometrically complex gamelevel (middle) and hundreds of simulated plant leaves (right).

c The Eurographics Association 2006.
°

M. Müller et al. / Position Based Dynamics

References

Appendix A:[BFA02] B RIDSON R., F EDKIW R., A NDERSON J.: Robust treatment of collisions, contact and friction for cloth animation. Proceedings of ACM Siggraph (2002), 594–603.

Gradient of the Normalized Cross Product[BMF03] B RIDSON R., M ARINO S., F EDKIW R.: Simulation ofclothing with folds and wrinkles. In ACM SIGGRAPH Symposium on Computer Animation (2003), pp. 28–36.
[BW98] BARAFF D., W ITKIN A.: Large steps in cloth simulation. Proceedings of ACM Siggraph (1998), 43–54.
[BWK03] BARAFF D., W ITKIN A., K ASS M.: Untangling cloth.
In Proceedings of the ACM SIGGRAPH (2003), pp. 862–870.
[CBP05] C LAVET S., B EAUDOIN P., P OULIN P.: Particle-based
viscoelastic fluid simulation. Proceedings of the ACM SIGGRAPH Symposium on Computer Animation (2005), 219–228.
[DSB99] D ESBRUN M., S CHRÖDER P., BARR A.: Interactive
animation of structured deformable objects. In Proceedings ofGraphics Interface ’99 (1999), pp. 1–8.

Constraint functions often contain normalized cross products. To derive the projection corrections, the gradient of the constraint function is needed. Therefore it is useful to know the gradient of a normalized cross product with respect to both arguments. Given the×p2
normalized cross product n = |pp1 ×p, the derivative with respect to1
2|
the first vector is ∂ nx∂ nx∂ nx 
∂n=
∂ p1 

[Jak01] JAKOBSEN T.: Advanced character physics Ű the fysix
engine. www.gamasutra.com (2001).
[MHTG05] M ÜLLER M., H EIDELBERGER B., T ESCHER M.,
G ROSS M.: Meshless deformations based on shape matching.
Proceedings of ACM Siggraph (2005), 471–478.
[NMK∗ 05] N EALEN A., M ÜLLER M., K EISER R., B OXERMAN
E., C ARLSON M.: Physically based deformable models in computer graphics. Eurographics 2005 state of the art report (2005).
[Pro95] P ROVOT X.: Deformation constraints in a mass-springmodel to describe rigid cloth behavior. Proceedings of Graphics
Interface (1995), 147Ű–154.
[Rat04] R ATCLIFF J.: Rocket - a viewer for real-time phyics simulations. www.physicstools.org (2004).
[THM∗ 03] T ESCHNER M., H EIDELBERGER B., M ÜLLER M.,
P OMERANERTS D., G ROSS M.: Optimized spatial hashing forcollision detection of deformable objects. Proc. Vision, Modeling, Visualization VMV 2003 (2003), 47–54.
[THMG04] T ESCHNER M., H EIDELBERGER B., M ÜLLER M.,
G ROSS M.: A versatile and robust model for geometrically complex deformable solids. Proceedings of Computer Graphics International (CGI) (2004), 312–319.
[VCMT95] VOLINO P., C OURCHESNE M., M AGNENATT HALMANN N.:
Versatile and efficient techniques for
simulating cloth and other deformable objects. Proceedings of
ACM Siggraph (1995), 137–144.

∂ p1z
∂ ny
∂ p1z
∂ nz
∂ p1z





(16)

p2z
0
−p2x



−p2y
T

p2x
+ n(n × p2 )0
(17)

Shorter and for both arguments we have
¡
¢
∂n
1
=
−p˜2 + n(n × p2 )T∂ p1|p1 × p2 |
¡
¢
1
∂n
=−
−p˜1 + n(n × p1 )T∂ p2|p1 × p2 |[Fed05] F EDOR M.: Fast character animation using particle dynamics. Proceedings of International Conference on Graphics,
Vision and Image Processing, GVIP05 (2005).

[ITF04] I RVING G., T ERAN J., F EDKIW R.: Invertible finite elements for robust simulation of large deformation. In Proceedings of the ACM SIGGRAPH Symposium on Computer Animation(2004), pp. 131–140.

∂ p1y
∂ ny
∂ p1y
∂ nz
∂ p1y


0
1
 −p2z=|p1 × p2 |p2y[Fau98] FAURE F.: Interactive solid animation using linearized
displacement constraints. In Eurographics Workshop on Computer Animation and Simulation (EGCAS) (1998), pp. 61–72.

[GHDS03] G RINSPUN E., H IRANI A., D ESBRUN M.,
S CHRODER P.:
Discrete shells.
In Proceedings of the
ACM SIGGRAPH Symposium on Computer Animation (2003).

∂ p1x
∂ ny
∂ p1x
∂ nz
∂ p1x

(18)
(19)
(20)

where p̃ is the matrix with the property p̃x = p × x.

Bending Constraint Projection
The constraint function for bending is C = arccos(d) − ϕ0 , whered = n1 · n2 = nT1 n2 . Without loss of generality we set p1 = 0 and get
×p3
4 . With d arccos(x) =for the normals n1 = |pp2 ×p
and n2 = |pp2 ×pdx
23|
2 ×p4 |
1
we get the following gradients:
−√
1−x2

¶
µ
1
∂ n1 T∇p3 C = − √n2 )(1 − d 2 ∂ p3µ¶
1
∂ n2 T∇p4 C = − √(n1 )1 − d 2 ∂ p4
µ
¶T
µ
¶
1
∂ n1
∂ n2 T∇p2 C = − √
(
n2 +
n1 )
∂ p2
1 − d 2 ∂ p2
∇p1 C = −∇p2 C − ∇p3 C − ∇p4 C

(21)
(22)
(23)
(24)

Using the gradients of normalized cross products, first computep2 × n2 + (n1 × p2 )d|p2 × p3 |p2 × n1 + (n2 × p2 )dq4 =|p2 × p4 |p3 × n2 + (n1 × p3 )d p4 × n1 + (n2 × p4 )dq2 = −−|p2 × p3 ||p2 × p4 |q3 =q1 = −q2 − q3 − q4

(25)
(26)
(27)
(28)

Then the final correction is√
wi 1 − d 2 (arccos(d) − ϕ0 )∆pi = −qi
∑ j w j |q j |2(29)

c The Eurographics Association 2006.
°