An Introduction To Physically Based Modeling

An Introduction to Physically Based Modeling: Differential Equation Basics
Andrew Witkin and David Baraff Robotics Institute Carnegie Mellon University
Please note: This document is 1997 by Andrew Witkin and David Baraff. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Differential Equation Basics

Andrew Witkin and David Baraff School of Computer Science Carnegie Mellon University
Initial Value Problems
Differential equations describe the relation between an unknown function and its derivatives. To solve a differential equation is to nd a function that satises the relation, typically while satisfying some additional conditions as well. In this course we will be concerned primarily with a particular class of problems, called initial value problems. In a canonical initial value problem, the behavior of the system is described by an ordinary differential equation (ODE) of the form x = f (x, t), where f is a known function (i.e. something we can evaluate given x and t,) x is the state of the system, and x is xs time derivative. Typically, x and x are vectors. As the name suggests, in an initial value problem we are given x(t0 ) = x0 at some starting time t0 , and wish to follow x over time thereafter. The generic initial value problem is easy to visualize. In 2D, x(t) sweeps out a curve that describes the motion of a point p in the plane. At any point x the function f can be evaluated to provide a 2-vector, so f denes a vector eld on the plane (see gure 1.) The vector at x is the velocity that the moving point p must have if it ever moves through x (which it may or may not.) Think of f as driving p from point to point, like an ocean current. Wherever we initially deposit p, the current at that point will seize it. Where p is carried depends on where we initially drop it, but once dropped, all future motion is determined by f . The trajectory swept out by p through f forms an integral curve of the vector eld. See gure 2. We wrote f as a function of both x and t, but the derivative function may or may not depend directly on time. If it does, then not only the point p but the the vector eld itself moves, so that ps velocity depends not only on where it is, but on when it arrives there. In that case, the derivative x depends on time in two ways: rst, the derivative vectors themselves wiggle, and second, the point p, because it moves on a trajectory x(t), sees different derivative vectors at different times. This dual time dependence shouldnt lead to confusion if you maintain the picture of a particle oating through an undulating vector eld.
Numerical Solutions
Standard introductory differential equation courses focus on symbolic solutions, in which the functional form for the unknown function is to be guessed. For example, the differential equation x = kx, where x denotes the time derivative of x, is satised by x = ekt .
B1
The derivative function
x = f (x,t) forms a vector field.
Vector Field
Figure 1: The derivative function f (x, t). denes a vector eld.
Start Here
Follow the vectors
Initial Value Problem

Figure 2: An initial value problem. Starting from a point x0 , move with the velocity specied by the vector eld.
An Introduction to Physically-Based Modeling
B2
Witkin/Baraff/Kass
Simplest numerical solution method Discrete time steps Bigger steps, bigger errors.
x(t + t) = x(t) + t f(x,t)

Euler's Method
Figure 3: Eulers method: instead of the true integral curve, the approximate solution follows a polygonal path, obtained by evaluating the derivative at the beginning of each leg. Here we show how the accuracy of the solution degrades as the size of the time step increases. In contrast, we will be concerned exclusively with numerical solutions, in which we take discrete time steps starting with the initial value x(t0 ). To take a step, we use the derivative function f to calculate an approximate change in x, x, over a time interval t, then increment x by x to obtain the new value. In calculating a numerical solution, the derivative function f is regarded as a black box: we provide numerical values for x and t, receiving in return a numerical value for x. Numerical methods operate by performing one or more of these derivative evaluations at each time step.
2.1
Eulers Method
The simplest numerical method is called Eulers method. Let our initial value for x be denoted by x0 = x(t0 ) and our estimate of x at a later time t0 + h by x(t0 + h), where h is a stepsize parameter. Eulers method simply computes x(t0 + h) by taking a step in the derivative direction, x(t0 + h) = x0 + h x(t0 ). You can use the mental picture of a 2D vector eld to visualize Eulers method. Instead of the real integral curve, p follows a polygonal path, each leg of which is determined by evaluating the vector f at the beginning, and scaling by h. See gure 3. Though simple, Eulers method is not accurate. Consider the case of a 2D function f whose integral curves are concentric circles. A point p governed by f is supposed to orbit forever on whichever circle it started on. Instead, with each Euler step, p will move on a straight line to a circle of larger radius, so that its path will follow an outward spiral. Shrinking the stepsize will slow the rate of this outward drift, but never eliminate it.
B3
Witkin/Baraff/Kass
Inaccuracy: Error turns x(t) from a circle into the spiral of your choice.
Instability: off to Neptune!
Two Problems
Figure 4: Above: the real integral curves form concentric circles, but Eulers method always spirals outward, because each step on the current circles tangent leads to a circle of larger radius. Shrinking the stepsize doesnt cure the problem, but only reduces the rate at which the error accumulates. Below: too large a stepsize can make Eulers method diverge. Moreover, Eulers method can be unstable. Consider a 1D function f = kx, which should make the point p decay exponentially to zero. For sufciently small step sizes we get reasonable behavior, but when h > 1/k, we have | x| > |x|, so the solution oscillates about zero. Beyond h = 2/k, the oscillation diverges, and the system blows up. See gure 4. Finally, Eulers method isnt even efcient. Most numerical solution methods spend nearly all their time performing derivative evaluations, so the computational cost per step is determined by the number of evaluations per step. Though Eulers method only requires one evaluation per step, the real efciency of a method depends on the size of the steps it lets you takewhile preserving accuracy and stabilityas well as on the cost per step. More sophisticated methods, even some requiring as many as four or ve evaluations per step, can greatly outperform Eulers method because their higher cost per step is more than offset by the larger stepsizes they allow. To understand how we go about improving on Eulers method, we need to look more closely at the error that the method produces. The key to understanding whats going on is the Taylor series: Assuming x(t) is smooth, we can express its value at the end of the step as an innite sum involving the the value and derivatives at the beginning: x(t0 + h) = x(t0 ) + h x(t0 ) + h2 h3 hn n x x(t0 ) + (t0 ) + . . . + x + ... 2! 3! n! t n
As you can see, we get the Euler update formula by truncating the series, discarding all but the rst two terms on the right hand side. This means that Eulers method would be correct only if all derivatives beyond the rst were zero, i.e. if x(t) were linear. The error term, the difference
B4
Witkin/Baraff/Kass
between the Euler step and the full, untruncated Taylor series, is dominated by the leading term, (h 2 /2) (t0 ). Consequently, we can describe the error as O(h 2 ) (read Order h squared.) Suppose x that we chop our stepsize in half; that is, we take steps of size h . Although this produces only about 2 one fourth the error we got with a stepsize of h, we have to take twice as many steps over any given interval. That means that the error we accumulate over an interval t0 to t1 depends linearly upon h. Theoretically, using Eulers method we can numerically compute x over an interval t0 to t1 with as little error as we want, by choosing a suitably small h. In practice, a great many timesteps might be required, depending on the error and the function f .
2.2
The Midpoint Method
If we were able to evaluate x as well as x, we could acheive O(h 3 ) accuracy instead of O(h 2 ) simply by retaining one additional term in the truncated Taylor series: x(t0 + h) = x(t0 ) + h x(t0 ) + h2 x(t0 ) + O(h 3 ). 2 (1)
Recall that the time derivative x is given by a function f (x(t), t). For simplicity in what follows, we will assume that the derivative function f does depends on time only indirectly through x, so that x = f (x(t)). The chain rule then gives x= f x = f f. x
To avoid having to evaluate f ,which would often be complicated and expensive, we can approximate the second-order term just in terms of f , and substitute the approximation into equation 1, leaving us with O(h 3 ) error. To do this, we perform another Taylor expansion, this time of the function of f , f (x0 + x) = f (x0 ) + x f (x0 ) + O( x2 ). (2) We rst introduce x into this expression by choosing x= so that f (x0 + h h h f (x0 )) = f (x0 ) + f (x0 ) f (x0 ) + O(h 2 ) = f (x0 ) + x(t0 ) + O(h 2 ), 2 2 2 h f (x0 ) 2
where x0 = x(t0 ). We can now multiply both sides by h (turning the O(h 2 ) term into O(h 3 )) and rearrange, yielding h2 h x + O(h 3 ) = h( f (x0 + f (x0 )) f (x0 ). 2 2 Substituting the right hand side into equation 1 gives the update formula x(t0 + h) = x(t0 ) + h( f (x0 + h f (x0 )). 2
This formula rst evaluates an Euler step, then performs a second derivative evaluation at the midpoint of the step, using the midpoint evaluation to update x. Hence the name midpoint method. The
B5
Witkin/Baraff/Kass
a b
a. Compute an Euler step x = t f(x,t) b. Evaluate f at the midpoint f mid = f x + x , t + t 2 2
c. Take a step using the midpoint value
x(t + t) = x(t) + t f mid
The Midpoint Method

Figure 5: The midpoint method is a 2nd-order solution method. a) an euler step is computed, b) the derivative is evaluated again at the steps midpoint, and the second evaluation is used to calculate the step. The integral curvethe actual solutionis shown as c. midpoint method is correct to within O(h 3 ), but requires two evaluations of f . See gure 5 for a pictorial view of the method. We dont have to stop with an error of O(h 3 ). By evaluating f a few more times, we can eliminate higher and higher orders of derivatives. The most popular procedure for doing this is a method called Runge-Kutta of order 4 and has an error per step of O(h 5 ). (The Midpoint method could be called Runge-Kutta of order 2.) We wont derive the fourth order Runge-Kutta method, but the formula for computing x(t0 + h) is listed below: k1 = h f (x0 , t0 ) k1 h k2 = h f (x0 + , t0 + ) 2 2 k2 h k3 = h f (x0 + , t0 + ) 2 2 k4 = h f (x0 + k3 , t0 + h) 1 1 1 1 x(t0 + h) = x0 + k1 + k2 + k3 + k4 . 6 3 3 6
Adaptive Stepsizes
Whatever the underlying method, a major problem lies in determing a good stepsize. Ideally, we want to choose h as large as possiblebut not so large as to give us an unreasonable amount of error, or worse still, to induce instability. If we choose a xed stepsize, we can only proceed as fast as the worst sections of x(t) will allow. What we would like to do is to vary h as we march
B6
Witkin/Baraff/Kass
forward in time. Whenever we can make h large without incurring too much error, we should do so. When h has to be reduced to avoid excessive error, we want to do that also. This is the idea of adaptive stepsizing: varying h over the course of solving the ODE. Here well be present adaptive stepsizing for Eulers method. The basic idea is as follows. Lets assume we have a given stepsize h, and we want to know how much we can consider changing it. Suppose we compute two estimates for x(t0 + h). We compute an estimate xa , by taking an Euler step of size h from t0 to t0 + h. We also compute an estimate xb by taking two Euler steps of size h/2, from t0 to t0 + h. Both xa and xb differ from the true value of x(t0 + h) by O(h 2 ). That means that xa and xb differ from each other by O(h 2 ). As a result, we can write that a measure of the current error e is e = |xa xb | This gives us a convenient estimate to the error in taking an Euler step of size h. Suppose that we are willing to have an error of as much as 104 per step, and that the current error is only 108 . Since the error goes up as h 2 , we can increase the stepsize to 104 108
1 2
h = 100h.
Conversely, if we currently had an error of 103 , and could only tolerate an error of 104 , we would have to decrease the stepsize to 104 103
1 2
h .316h.
Adaptive stepsizing is a highly recommended technique.
Implementation
The ODEs we will want to solve may represent many thingsfor instance, a collection of masses and springs, some rigid bodies, or a deformable object. We want to implement ODE solvers and the models on which they operate in a way that isolates each from the internal details of the other. This will make it possible to change solvers easily, and also make the solver code reusable. Fortunately, this kind of modularity is not difcult to acheive, since all solvers can be expressed in terms of a small, stereotyped set of operations. Presumably, the system of ODE-governed objects will be embodied in a structure of some kind. The approach is to write type-specic code that operates on this structure to perform the standard operations, then to implement solvers in terms of these generic operations. From the solvers viewpoint, the system on which it operates is a black-box function f (x, t). The solver needs to be able to evaluate f , as required, at any values of x and t, and then to install the updated x and t when a time step is taken. To support these operations, the object that represents the ODE being solved must be able to handle these requests from the solver: Return dim (x). Since x and x may be vectors, the solver must know their length, to allocate storage, perform vector arithmetic ops, etc. Get/set x and t. The solver must be able to install new values at the end of a step. In addition, a multi-step method must set x and t to intermediate values in the course of performing derivative evaulations.
B7
Witkin/Baraff/Kass
Evaluate f at the current x and t. In an object-oriented language, these operations would naturally be implemented as generic functions that are handled in a type-specic way. In a non-object-oriented language generic functions would be faked by installing pointers to type-specic functions in structure slots, or simply by passing the function pointers as arguments to the solver. Later on we will consider in detail how these operations are to be implemented for specic models such as particle-and-spring systems.
References
[1] W.H. Press, B.P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, England, 1988.
B8
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: Particle System Dynamics

Andrew Witkin Robotics Institute Carnegie Mellon University
Please note: This document is 1997 by Andrew Witkin. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Particle System Dynamics

Andrew Witkin School of Computer Science Carnegie Mellon University
Introduction
Particles are objects that have mass, position, and velocity, and respond to forces, but that have no spatial extent. Because they are simple, particles are by far the easiest objects to simulate. Despite their simplicity, particles can be made to exhibit a wide range of interesting behavior. For example, a wide variety of nonrigid structures can be built by connecting particles with simple damped springs. In this portion of the course we cover the basics of particle dynamics, with an emphasis on the requirements of interactive simulation.
Phase Space
The motion of a Newtonian particle is governed by the familiar f = ma, or, as we will write it here, x = f/m. This equation differs from the canonical ODE developed in the last chapter because it involves a second time derivative, making it a second order equation. To handle a second order ODE, we convert it to a rst-order one by introducing extra variables. Here we create a variable v to represent velocity, giving us a pair of coupled rst-order ODEs v = f/m, x = v. The position and velocity x and v can be concatenated to form a 6-vector. This position/velocity product space is called phase space. In components, the phase space equation of motion is [x1 , x2 , x3 , v1 , v2 , v3 ] = [v1 , v2 , v3 , f 1 /m, f 2 /m, f 3 /m], which, assuming force is a function of x and t, matches our canon ical form x = f(x, t). A system of n particles is described by n copies of the equation, concatenated to form a 6n-long vector. Conceptually, the whole system may be regarded as a point moving through 6n-space. We can still visualize the phase-space ODE in terms of a planar vector eld, though only for a 1D particle, by letting one axis represent the particles position and the other, its velocity. If each point in the phase plane represents a pair [x, v], then the derivative vector is [v, f /m]. All the ideas of integral curves, polygonal approximations, and so forth, carry over intact to phase space. Only the interpretation of the trajectory is changed.
Basic Particle Systems
In implementing particle systems, we want to maintain two views of our model: from outside, especially from the point of view of the ODE solver, the model should look like a monolitha point in a high-dimensional space, whose time derivative may be evaluated at will. From within, the model should be a structureda collection of distinct interacting objects. This duality will be recurring theme in the course.
C1
A particle simulation involves two main partsthe particles themselves, and the entities that apply forces to particles. In this section we consider just the former, deferring until the next section the specics of force calculation. Our goal here is to describe structures that could represent a particle and a system of particles, and to show in a concrete way how to implement the generic operations required by ODE solvers. Particles have mass, position, and velocity, and are subjected to forces, leading to an obvious structure denition, which in C might look like: typedef struct{ float m; float *x; float *v; float *f; } *Particle; /* /* /* /* mass */ position vector */ velocity vector */ force accumulator */
In practice, there would probably be extra slots describing appearance and other properties. A system of particles could be represented in an equally obvious way, as typedef struct{ Particle *p; /* array of pointers to particles */ int n; /* number of particles */ float t; /* simulation clock */ } *ParticleSystem; Assume that we have a function CalculateForces() that, called on a particle system, adds the appropriate forces into each particles f slot. Dont worry for know about what that function actually does. Then the operations that comprise the ODE solver interface could be written as follows: /* length of state derivative, and force vectors */ int ParticleDims(ParticleSystem p){ return(6 * p->n); }; /* gather state from the particles into dst */ int ParticleGetState(ParticleSystem p, float *dst){ int i; for(i=0; i < p->n; i++){ *(dst++) = p->p[i]->x[0]; *(dst++) = p->p[i]->x[1]; *(dst++) = p->p[i]->x[2]; *(dst++) = p->p[i]->v[0]; *(dst++) = p->p[i]->v[1]; *(dst++) = p->p[i]->v[2]; } }
C2
Witkin/Baraff/Kass
/* scatter state from src into the particles */ int ParticleSetState(ParticleSystem p, float *src){ int i; for(i=0; i < p->n; i++){ p->p[i]->x[0] = *(src++); p->p[i]->x[1] = *(src++); p->p[i]->x[2] = *(src++); p->p[i]->v[0] = *(src++); p->p[i]->v[1] = *(src++); p->p[i]->v[2] = *(src++); } } /* calculate derivative, place in dst */ int ParticleDerivative(ParticleSystem p, float *dst){ int i; Clear_Forces(p); /* zero the force accumulators */ Compute_Forces(p); /* magic force function */ for(i=0; i < p->n; i++){ *(dst++) = p->p[i]->v[0]; /* xdot = v */ *(dst++) = p->p[i]->v[1]; *(dst++) = p->p[i]->v[2]; *(dst++) = p->p[i]->f[0]/m; /* vdot = f/m */ *(dst++) = p->p[i]->f[1]/m; *(dst++) = p->p[i]->f[2]/m; } } Having dened these operations, and assuming some utility routines and temporary vectors, an Euler solver be written as void EulerStep(ParticleSystem p, float DeltaT){ ParticleDeriv(p,temp1); /* get deriv */ ScaleVector(temp1,DeltaT) /* scale it */ ParticleGetState(p,temp2); /* get state */ AddVectors(temp1,temp2,temp2); /* add -> temp2 */ ParticleSetState(p,temp2); /* update state */ p->t += DeltaT; /* update time */ } The structures representing a particle and a particle system are shown visually in gures 1 and 2. The interface between a particle system and a differential equation solver is illustrated in gure 3.
Forces
All particles are essentially alike. In contrast, the objects that give rise to forces are heterogeneous. As a matter of implementation, we would like to make it easy to extend the set of force-producing
C3
Witkin/Baraff/Kass
x v f
m
Position Velocity
Position in Phase Space
Force Accumulator mass
Particle Structure
Figure 1: A particle may be represented by a structure containing its position, velocity, force, and mass. The six-vector formed by concatenating the position and velocity comprises the points position in phase space.
particles n time particles n time
x x x x x v v v v v f f f f f
m m m m m
x v f
m
Particle Systems
Figure 2: A bare particle system is essentially just a list of particles.
C4
Witkin/Baraff/Kass
Solver Interface
Particle System
particles n time particles n time
Dim(State) Get/Set State
Diffeq Solver
6n x1 v1 x2 v2 f f v 1 m1 v 2 m2 1 2
xn vn f v n mn n
Deriv Eval
Figure 3: The relation between a particle system and a differential equation solver. objects without modifying the basic particle system model. We accomplish this by having the particle system maintain a list of force objects, each of which has access to any or all particles, and each of which knows how to apply its own forces. The Calculateforces function, used above, simply traverses the list of force structures, calling each of their ApplyForce functions, with the particle system itself as sole argument. This leaves the real work of force calculation to the individual objects. See gures 4 and 5 Forces can be grouped into three broad categories: Unary forces, such as gravity and drag, that act independently on each particle, either exerting a constant force, or one that depends on one or more of particle position, particle velocity, and time. n-ary forces, such as springs, that apply forces to a xed set of particles. Forces of spatial interaction, such as attraction and repulsion, that may act on any or all pairs of particles, depending on their positions. Each of these raises somewhat different implementation issues. We will now consider each in turn.
4.1
Unary forces
Gravity. Global earth gravity (as opposed to particle-particle attraction) is trivial to implement. The gravitational force on each particle is f = mg, where g is a constant vector (presumably pointing down) whose magnitude is the gravitational constant. If all particles are to feel the same gravity, which they need not in a simulation, then gravitational force is applied simply by traversing the
C5
Witkin/Baraff/Kass
particles
time
forces nforces
x x v v f f
m m
x v f
m
F F F
A list of force objects to invoke
Particle Systems, with forces

Figure 4: A particle system augmented to contain a list of force objects. Each force object points at the particles that it inuences, and contains a function that knows how to compute the force on each affected particle.
x x v v ff
x x v v ff
x 1 x v v ff
m m
m m m m
Clear Force Accumulators
F F F F F F
Invoke apply_force functions

x x v v ff
m m
x x v v ff
x x v v ff
3
Return [v, f/m,] to solver.
m m m m
Deriv Eval Loop

Figure 5: The derivative evaluation loop for a particle system with force objects.
C6
Witkin/Baraff/Kass
Force Law:
Particle system k
f drag = -kdragv
F
sys
x v f
m
apply_fun
p->f -= F->k * p->v
A Force Object: Viscous Drag

Figure 6: Schematic view of a force object implementing viscous drag. The object points at the particle to which drag is being applied, and also points to a function that implements the force law for drag. systems particle list, and adding the appropriate force into each particles force accumulator. Gravity is basic enough that it could reasonably be wired it into the particle system, rather than maintaining a distinct gravity object. Viscous Drag. Ideal viscous drag has the form f = kd v, where the constant kd is called the coefcient of drag. The effect of drag is to resist motion, making a particle gradually come to rest in the absence of other inuences. It is highly reccommended that at least a small amount of drag be applied to each particle, if only to enhance numerical stability. Excessive drag, however, makes it appear that the particles are oating in molasses. Like gravity, drag can be implemented as a wired-in special case. A force object implementing viscous drag is shown in gure 6.
4.2
n-ary forces
Our canonical example of a binary force is a Hooks law spring. In a basic mass-and-spring simulation, the springs are the structural elements that hold everything together. The spring forces between a pair of particles at positions a and b are l l l , fb = fa , (1) |l| |l| where fa and fb are the forces on a and b, respectively, l = a b, r is the rest length, ks is a spring constant, and kd is a damping constant. the time derivative of l, is just va vb , the difference l, between the two particles velocities. In equation 1, the spring force magnitude is proportional to the difference between the actual length and the rest length, while the damping force magnitude is proportional to a and bs speed of fa = ks (|l| r ) + kd
C7
Witkin/Baraff/Kass
Force Law:
vx x f 1 = - ks( x - r ) + kd x x f 2 = -f 1
x v f
m
Particle system
p1
x v f
m
kd ks
sys apply_fun
p2
Damped Spring
Figure 7: A schematic view of a force object implementing a damped spring that attaches particles p1 and p2 . approach. Equal and opposite forces act on each particle, along the line that joins them. The spring damping differs from global drag in that it acts symmetrically on the two particles, having no effect on the motion of their common center of mass. Later, we will learn a general procedure for deriving this kind of force expression. A damped spring can be implemented straightforwardly as a structure that points to the pair of particles it connects. The code that applies the forces according to equation 1 fetches the positions and velocities from the two particle structures, performs its calculations, and sums the results into the particles force accumulators. In an object-oriented environment, this operation would be implemented as a generic function. In bare C, the force object would contain a pointer to an ordinary C function. A force object for a damped spring is shown in gure 7
4.3
Spatial Interaction Forces
A spring applies forces to a xed pair of particles. In contrast, spatial interaction forces may act on any pair (or n-tuple) of particles. For local interaction forces, particles begin to interact when they come close, and stop when they move apart. Spatially interacting particles have been used as approximate models for uid behavior, and large-scale particle simulations are widely used in physics [1]. A complication in large-scale spatial interaction simulations is that the force calculation is O(n 2 ) in the number of particles. If the interactions are local, efciency may be improved through the use of spatial buckets.
C8
Witkin/Baraff/Kass
User Interaction
An interactive mass-and-spring simulation is an ideal rst implementation project in physically based modeling, because such simulations are relatively easy to implement, and because interactive performance can be acheived even on low-end computers. The main ingredients of a basic massand-spring simulation are model construction and model manipulation. Model construction can be a simple matter of mouse-clicking to create particles and connect them with springs. Interactive manipulation requires no more than the ability to grab and drag mass points. Although there is barely any difference mathematically between 2D and 3D simulations, supporting 3D user interaction is more challenging. Most of the implementation issues are standard, and will not be dealt with here. However, we give a few useful tips: Controlled particles. Particles whose motion is not governed by forces provide a number of interesting possibilities. Fixed particles serve as anchors and pivots. Particles whose motion is procedurally controlled (e.g. moving on a circle) can provide dynamic elements such as motors. All that need be done to implement controlled particles is to prevent the ODE solver from updating their positions. One subtle point, though, is that the velocities as well as positions of controlled particles must be maintained at their correct values. Otherwise, velocity-dependent forces such as damped spring forces will behave incorrectly. Structures. A variety of interesting non-rigid structuresbeams, blocks, etc.can be built out of masses and springs. By allowing several springs to meet at a single particle, these pieces can be connected with a variety of joints. With some experimentation and ingenuity it is possible to construct entire mechanisms, complete with motors, out of masses and springs. The topic of regular mass-and-spring lattices as an approximation to continuum models will be discussed later.[2] Mouse springs. The simplest way to manipulate mass-and-spring models is to use the mouse directly to control the positions of grabbed particles. However, this method is not recommended because very rapid mouse motions can lead to stability problems. These problems can be avoided by coupling the grabbed particle to the mouse position using a spring.
Energy Functions
Generically, the position-, velocity-, and time-dependent formulae that we use to calculate forces are known as force laws. Forces laws are not laws of physics. Rather, they form part of our description of the system we are modeling. Some of the standard ones, like linear springs and dashpots, represent time-honored idealizations of the behavior of real materials and structures. However, if we wanted to accurately model the behavior of a pair of particles connected by, say, a strand of gooey taffy, the resulting equations would probably be even messier than the taffy. Often, we can regard force laws as things we design to hold things in a desired conguration for instance a spring with nonzero rest length makes the points it connects want to be a xed distance apart. In many cases it is possible to specify the desired conguration by giving a function that reaches zero exactly when things are happy. We can call this kind of function a behavior function. For example, a behavior function that says that two particles a and b should be in the same place is just C(a, b) = a b (which is a vector expression each of whose components is supposed to vanish.) If instead we want a and b to be distance r apart, then a suitable behavior function is C(a, b) = |a b| r (which is a scalar expression.) Later on, when we study constrained dynamics, we will use this kind of function as a way to
C9
Witkin/Baraff/Kass
specify constraints, and we will consider in detail the problem of maintaining such constraints accurately. For now, we will be content to impose forces that pull the system toward the desired state, but that compete with other forces. These energy-based forces can be used to impose approximate, sloppy constraints. However, attempting to make them accurate by increasing the spring constant leads to numerical instability.[3] Converting a behavior function C(x1 . . . , xn ) into a force law is a pure cookbook procedure. We rst dene a scalar potential energy function E= ks C C, 2
where ks is a generalized stiffness constant. Since the force due to a scalar potential is minus the energy gradient, the force on particle xi due to C is fi = E C = ks C xi xi .
In general C is a vector, and this expression denotes its product with the transpose of the Jacobian matrix C/ xi . We will look much more closely at this kind of expression when we study constraint methods, and in particular Lagrange multipliers. For now, it is sufcent to think of the forces f i as generalized spring forces that attract the system to states that satisfy C = 0. When a behavior function depends on a number of particles positions, we get a different force expression for each by using Cs derivative with respect to that particle. The force we just dened isnt quite the one we want: in the absence of any damping, this conservative force will cause the system to oscillate about C = 0. To add damping, we modify the force expression to be C , (2) fi = (ks C kd C) xi where kd is a generalized damping constant, and C is the time derivative of C. Note that when you derive expressions for C, you will be using the fact that xi = vi . So, in a trivial case, if C = x1 x2 , = v1 v 2 . it follows that C As an extremely simple example, we take C = x1 x2 , which wants the points to coincide. We have C C = I, = I, x1 x2 where I is the identity matrix. The time derivative is C = v1 v2 . So, substituting into equation 2, we have f1 = ks (x1 x2 ) kd (v1 v2 ), f2 = ks (x1 x2 ) + kd (v1 v2 ),
which is just the force law for a damped zero-rest-length spring. As another example, we use the behavior function C = |l| r , where l = x1 x2 , which says the two points should be distance r apart. Its derivative w.r.t. l is C l = , l |l| An Introduction to Physically-Based Modeling C10 Witkin/Baraff/Kass
a unit vector in the direction of l. Then, since l = x1 x2 , C C = , x1 l The time derivative of is C C = . x2 l
l l l (v1 v2 ) C= = . |l| |l|
These expressions are then substituted into the general expression of equation 2 to get the forces. You should verify that this produces the damped spring force of equation 1.
Particle/Plane Collisions and Contact
The general collision and contact problem is difcult, to say the least. Later in the course we will examine rigid body collision and contact. Here we only consider, in bare bones form, the simplest case of particles colliding with a plane (e.g. the ground or a wall.) Even these simple collision models can add signicant interest to an interactive simulation.
7.1
Detection
There are two parts to the collision problem: detecting collisions, and responding to them. Although general collision detection is hard, particle/plane collision detection is trivial. If P is a point on the plane, and N is a normal, pointing inside (i.e. on the legal side of the barrier,) then we need only test the sign of (X P) N to detect a collision of point X with the plane. A value greater than zero means its inside, less than zero means its outside (where it isnt allowed to be) and zero means its in contact. If after an ODE step a collision is detected, the right thing to do is to solve (perhaps by linear interpolation between the old and new positions) for the instant of contact, and roll back the whole system to that time. A less accurate but easier alternative is just to displace the point that has collided.
7.2
Response
To describe collision response, we need to partition velocity and force vectors into two orthogonal components, one normal to the collision surface, and the other parallel to it. If N is the normal to the collision plane, then the normal component of a vector x is xn = (N x)x, and the tangential component is xt = x xn . The simplest collision to consider is an elastic collision without friction. Here, the normal component of the particles velocity is negated, whereafter the particle goes its merry way. In an inelastic collision, the normal velocity component is instead multiplied by r , where r is a constant between zero and one, called the coefcient of restitution. At r = 0, the particle doesnt bounce at all, and r = .9 is a superball.
7.3
Contact
If a particle is on the collision surface, with zero normal velocity, then it is in contact. If a particle is pushed into the contact plane (N f < 0) a contact force fc = (N f)f is exerted, exactly canceling
C11
Witkin/Baraff/Kass
the normal component of f . However, if the applied force points away from the contact plane, no contact force is exerted (unless the surface is sticky,) the particle begins to accelerate away from the surface, and contact is broken. In the very simplest linear friction model, the frictional force is k f (f N)vt , a drag force that acts in the tangential direction, with magnitude proportional to the normal force. To model a perfectly non-slippery surface, vt is simply zeroed.
References
[1] R.W Hocknew and J.W. Eastwood. Computer Simulation Using Particles. Adam Hilger, New York, 1988. [2] Gavin S. P. Miller. The motion dynamics of snakes and worms. Computer Graphics, 22:169 178, 1988. [3] Andrew Witkin, Kurt Fleischer, and Alan Barr. Energy constraints on parameterized models. Computer Graphics, 21(4):225232, July 1987.
C12
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: Energy Functions and Stiffness

Michael Kass Pixar
Please note: This document is 1997 by Michael Kass. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Energy Functions and Stiffness

Michael Kass Pixar
In earlier sections of these notes, we remarked that there are two main types of problems that arise in practical numerical solution of differential equations: accuracy and stability. Accuracy refers to the small errors introduced when trying to make a numerical approximation to the solution of a continuous differential equation. Stability refers to the problem of whether or not the approximations result in a method which converges at all. In some engineering contexts, high accuracy is very important, but for physically-based computer animation, accuracy tends to be much less important than stability. Small errors in simulations go largely un-noticed in computer graphics, but instabilities usually make the results of a simulation completely unusable. In addition, stability is very often the limiting factor in the computational speed of simulation because we can usually increase the time step of the differential equation solver until the limit of stability to make a simulation run as fast as possible. Here we will examine some of the causes of instability in differential equations resulting from physically-based modeling and discuss ways of dealing with it. For illustrative purposes, instead of considering correct Newtonian mechanics governed by the law F = ma, we will analyze a simpler physical model guided by the law F = mv. The simpler model yields physics in which nothing moves unless a force is acting on it, in accordance with the theories of Aristotle in early Greece. The advantage of the simpler model is that we can construct interesting behavior directly from energy functions without the need to introduce frictional and damping forces. As a result, the model is particularly easy to analyze, yet it contains the important elements required to understand the source of instability. The simplest physical system which exhibits the key problems that cause instability is a particle of unit mass in two dimensions affected by a potential energy E(x, y). The force on the particle is given by the vector (E/x, E/y) so using Aristotelian physics, the equations of motion of the particle are simply: X/t = (E/x), and Y/t = (E/y). The behavior of the system depends, of course, on the precise nature of the energy function. Suppose we consider the behavior of the above system near an energy minimum. If we consider a small enough region, the energy function will be reasonably well approximated by the rst few terms in its Taylor expansion. The constant term in the Taylor expansion does not affect the behavior of the particle since only relative changes in potential energy are signicant. The gradient vanishes at the minimum, so the rst derivative term has no effect. Thus the second derivative term is the rst term which importantly affects the behavior near a minimum. Since the second derivative term gives rise to a quadratic energy, it is particularly interesting to look at how the system behaves when the energy function is a quadratic. Let us begin by considering the quadratic energy function E = (x2 + y2 )/2. The energy has a minimum at zero where its value is zero. The resulting force on the particle is simply F = E = (x, y). Hence the force on the particle will be a vector from the position of the particle back to the origin. Under Aristotelian dynamics, the particle will take a straight line path from wherever it starts to the origin. The particle will slow down as it approaches because the gradient gets smaller
D1
Figure 1: Symmetric energy function E = (x2 + y2 )/2
D2
Witkin/Baraff/Kass
Figure 2: Stretched energy function E = ( x2 + y2 )/2 and smaller. If the particle is placed at the origin, its velocity will be zero since the gradient vanishes. Suppose we try to compute the motion of the particle using Eulers method. If P0 represents the starting position of the particle, its position at the end of the rst time step is given by P1 = P0 + h E = (1 h)P0 where h is the time step. Obviously, if h is very small, it will take a large number of iterations before the particle moves close to the minimum at (0, 0). As h gets larger, the particle moves further in one iteration, so the simulation proceeds more rapidly. In the range 0 < h < 1, the iteration yields exponential decay with time constant proportional to 1/ h. If h = 1, the particle moves directly to the minimum in one step, which is the fastest possible simulation. If h is increased beyond this point, the simulation oscillates, ipping from one side of the origin to the other on each iteration. Nonetheless, as long as h < 2, the position of the simulated particle converges towards the energy minimum at (0, 0). On each iteration, the sign of the vector from the origin to the particle changes, but the magnitude decrease. If h > 2, however, the iteration diverges. On each iteration, the position of the particle reects around the origin and gets further away. Figure 1 shows some level curves for the energy function and some paths to the minimum. If the energy function is a symmetric quadratic as above and the step size is set appropriately, Eulers method works quite well. How could we complain about being able to jump to the minimum in a single step? Even if we do not know the correct step size to use, it should be possible to get to the minimum in a relatively small number of step. It would be even better if the iteration were stable for all step sizes, but this can nonetheless be regarded as a success for Eulers method. Something very different occurs if we take the original symmetric energy function and stretch it out in one direction. Suppose we have E = ( x2 + y2 )/2. Then the resulting force on the particle is F = E = ( x, y) and the force no longer points to the minimum unless x = 0, y = 0, or = 1. Now consider what happens if we use Eulers method. We begin at some point P0 = (x0 , y0 ). The next point is P1 = ((1 h )x0 , (1 h)y0 ). Note that the y component of the Euler iterates in this case will diverge if h > 2. Thus, even at the largest step size that will still maintain stability, the value of x decreases each iteration by only a factor of (1 2 ), which describes exponential decay with a time constant proportional to 1/ . Note the profound change that resulted from just a simple scaling of the energy function. With the symmetric energy, we can get to the minimum in a single iteration if we pick the right step size. Otherwise, if we pick a stable step size, we approach the minimum with a time constant proportional the the inverse of the step size. With the stretched energy, however, even if we pick the largest pos-
D3
Witkin/Baraff/Kass
sible stable step size, we approach the minimum with time constant proportional to 1/ . How bad is this? Arbitrarily bad. In real problems, 1/ can be on the order of 106 or 109 or worse. Clearly this is a serious difculty. In order to understand the nature of the difculty, consider what the stretched energy looks like. As gets very small, the energy function stops looking like a symmetric ball, and looks increasingly like a steep ravine. At the bottom of the ravine, there is a very small slope towards the minimum, but anywhere else, the dominant slope is towards the bottom of the ravine. As a result, a simulation using Eulers method will primarily send the particle back and forth, up and down the sides of the ravine, while very, very slowly drifting towards the minimum. This is shown in g. 2. The problem is that Eulers method is based on a local planar approximation of the energy function. When the energy function begins to look like a highly curved ravine, that approximation becomes very poor unless the neighborhood and hence the step size is very small. As a result, Eulers method becomes extremely slow in these cases. Note that when we stretch the energy function, we change the time constants on the true solutions to the differential equations. In the symmetric case, the solution is x = x0 et , y = y0 et . In the stretched case, the solution is x = x0 e t , y = y0 et . The time constants for the symmetric case are the same, but in the stretched case, the time constants have a ratio of e. This is the usual root of stiffness in differential equations. As we saw, the maximum step size is set by the process with the fastest time constant, and this can make simulation of the longer time constant unbearably slow. In order to deal with the difculties of stiffness, we have two main choices. We can either use a method based on a more accurate local model of the energy landscape, or we can try to formulate a different but related set of differential equations capable of solving our original problem. Reformulating the problem to eliminate the stiffness is preferred when possible, but very often there is no choice besides simulating the stiff equation. In many uses of physically-based modeling, people have been tempted to use springs of various sorts to enforce constraints. The problem is that keeping constraints satised with springs requires making the spring constants very large, which gives rise to stiff systems of equations. In many cases, Lagrange multipliers or changes of variables can make it possible to achieve the same results with differential equations that are much less stiff. If it is not possible to reformulate the differential equations to remove the stiffness, the only practical way to solve them is to use a technique which combats stiffness using a better local model of the forces than Eulers method. In general, these are known as implicit methods, one of which is treated in detail in Continuum Dynamics for Computer Graphics elsewhere in these notes. The details of implicit differential equation solvers are treated in a wide variety of texts. Here we attempt to give a little intuition about how they work. Remember that the gradient of the symmetric energy points towards the minimum, but not the gradient of the stretched energy. This is a large part of the difculty. If we keep moving in the gradient direction, we will get to the minimum very slowly. In order to take a very large step, therefore, a differential equation solver should move directly towards the minimum, even when faced with the stretched energy. How can a differential equation solver know which way to go to get to the minimum? In the general case, this can be arbitrarily difcult, but for quadratic energies, it is particularly simple. Quadratic energies have linear gradients. The minimum occurs when the gradient vanishes, so nding the point where this occurs can be done by solving a set of linear equations. For the case of the stretched energy, the gradient is ( x, y) which clearly vanishes only at the point (0, 0). Thus, by solving an appropriate linear equation, a differential equation solver can move towards the minimum of a quadratic energy function. In general, the best linear estimate for the gradient of an energy function near a point x0 is given
D4
Witkin/Baraff/Kass
by the expression E(x) E(x0 ) + H (x x0 ) where H is the matrix of second derivatives or Hessian matrix given by H = 2 E/xi x j . Since H is symmetric, by a suitable orthogonal change of basis, it is always possible to turn H into a diagonal matrix. For the symmetric energy function, H is already diagonal: it is simply a two by two identity matrix. For the stretched energy, H is also already diagonal with entries ( , 1). In the basis that diagonalizes H, the condition number of H is dened to be the ratio between the magnitudes of the largest and smallest diagonal entries of H. If the condition number is small, then we have a nearly symmetric energy function and the equations will not be stiff. If the condition number is large, it means that if we consider the two dimensions which gave rise to the smallest and largest diagonal entries of H, we have the situation of the stretched energy function with 1/ corresponding to the condition number. Clearly, if the condition number is large, we are potentially in a troublesome situation. The symmetric and stretched energy functions, point out a canonical difculty in solving differential equations as well as a canonical solution. If we look at a suitable neighborhood of any smooth energy function, we can approximate its behavior as a quadratic. If we pick the right coordinates, the Hessian matrix will be diagonal, so we can characterize the local behavior in those coordinates by the diagonal entries of the Hessian. If the entries are all of approximately the same magnitude, then the energy function is reasonably symmetric at least locally, and our numerical problems are relatively minor. If the entries are of wildly different magnitudes, then if we pick the two coordinates corresponding the the smallest and largest magnitude entries, we have an energy function that looks locally like a steep ravine. In order to take large steps in an energy function of this shape, we need to use an implicit method which will examine the matrix H to nd the axis of the ravine, so it can travel down the ravine directly towards the minimum. We will see an example of one such method in An Introduction to Continuum Dynamics for Computer Graphics.
D5
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: An Introduction to Continuum Dynamics for Computer Graphics
Michael Kass Pixar
Please note: This document is 1997 by Michael Kass. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
An Introduction to Continuum Dynamics for Computer Graphics

Michael Kass Pixar
1 Introduction
Mass-spring systems have been widely used in computer graphics because they provide a simple means of generating physically realistic motion for a wide range of situations of interest. Even though the actual mass of a real physical body is distributed through a volume, it is often possible to simulate the motion of the body by lumping the mass into a collection of points. While the exact coupling between the motion of different points on a body may be extremely complex, it can frequently be approximated by a set of springs. As a result, mass-spring systems provide a very versatile simulation technique. In many cases, however, there are much better simulation techniques than directly approximating a physical system with a set of mass points and springs. The motion of rigid bodies, for example, can be approximated by a small set of masses connected by very stiff springs. Unfortunately, the very stiff springs wreak havoc on the numerical solution methods, so it is much better to use a technique based on the rigid-body equations of motion. Another important special case arises with our subject here: elastic bodies and uids. In many cases, they can be approximated by regular lattices of mass points and springs. While these systems can be simulated with simple general methods, the performance is often very disappointing. The reason for the poor performance is relatively easy to understand because the regular lattices are particularly amenable to analysis. Understanding the difculties leads to formulations of solution methods which exploit the regular structure and thereby achieve far more efcient behavior.
2 Spring with Mass

In idealized mass-spring systems, the springs are considered to have zero mass. This is primarily because such systems are easy to simulate, not because real springs are massless. Here we will use idealized springs to simulate a more realistic spring with mass distributed along its length. With this change, the behavior of the simulated spring becomes a great deal more complex than an idealized spring and includes full wave behavior along its length, much like a Slinky. Figure 1 shows a large zero-rest length spring modeled as a collection of n smaller springs in series. If the mass of the entire spring is m, then the masses between the smaller springs are mi = m/(n + 1). For simplicity and ease of exposition, we will initially consider the one-dimensional problem in which all forces and motion are limited to the x-axis. Let xi denote the position of the ith mass point. If the spring constant of the ith spring is ki , the energy of the collection of small springs
E1
m0
m1
m2
mn1
mn
...
Figure 1: Single spring (top) is simulated by a collection of smaller springs connected in series (bottom).
E2
Witkin/Baraff/Kass
can be written E= 1 2
ki (xi xi1 )2
i=1
(1)
When the small springs are equally stretched out, their energy should sum to the energy of the entire spring. This will happen when ki = nk where k is the spring constant of the entire spring. Substituting into the energy expression yields E= nk 2
n
(xi xi1 )2
i=1
(2)
3 Continuous limit
Additional insight into this energy function can be achieved by examining the limit as the number of springs n goes to innity. Then x becomes a continuous function of a parameter u instead of being indexed by i . In the limit, the energy function can be written E = lim nk n 2
n i=1
ki (xi xi1 )2 = lim
k n 2n
n i=1
x(i/n) x((i 1)/n) 1/n
(3)
The term being summed on the right-hand side of converges to the square of the derivative of x. The summation converges to an integral. As a consequence, the energy function can be written E= k 2
1 0
x(u) u
(4)
in the limit as n goes to innity. It is also possible to arrive at this energy expression directly from arguments about the innitesimal properties of springs. The continuous energy expression in equation 4 is very similar to a continuous energy expression of great importance in computer graphics: the energy k E= 2
1 0
2 x(u) u2
(5)
of a natural spline. The only difference is that the spring energy of equation 4 depends on rst derivatives, while the spline energy of equation 5 depends on second derivatives. When these energy functions are extended to surfaces, the spring energy corresponds to the behavior of a membrane, while the spline energy corresponds to the energy of a thin plate.
4 Discretization
Suppose we are given a continuous energy function such as equation 4 or 5 and we wish to simulate its behavior with a computer. Before we can simulate the behavior, we must decide on a representation for the continuous function x(u). There are two basic choices: we can represent x(u) by its values at a set of sample points, or we can represent it by a collection of coefcients of local basis functions (e.g. the control points of a spline). The rst is known as the nite-difference method and the second as the nite-element method. Here we will examine the nite-difference method.
E3
Witkin/Baraff/Kass
If the function x(u) is represented by n + 1 samples evenly distributed on the interval from zero to one with xi = x(i/n), 0 i n, the derivatives can be approximated by the expressions x u and 2 x u2
x=ih
xi+1 xi h
(6)
x=ih
xi+1 2xi + xi1 h2
(7)
where h = 1/n is the spacing between sample values of u. The approximation of equation 7 is derived by using the approximation of equation 6 twice, once for each differentiation. If the approximation of equation 6 is combined with the energy function of equation 4 and the integral in the energy function is converted into a sum over the sample values, then the result is exactly the energy expression of 2. We therefore have two different interpretations of the spring energy in equation 2. One interpretation is the approximation of a large massive spring by a collection of smaller springs connected to point masses. The other interpretation is that the energy in equation 2 is a discrete approximation to the continuous energy function of 4. If the approximation of equation 7 is combined with the continuous energy function of 5, the result is the following discrete energy function: n2 k E= 2
n
(xi+1 2xi + xi1 )2

i=1
(8)
Like the energy function in equation 2, this energy has an interpretation in terms of springs, but it is not nearly so simple. The energy corresponds to springs between xi and the midpoint between xi1 and xi+1 .
5 Eulers Method
In order to simulate the motion of physical system from a discrete energy like equation 2, we have to compute the force from the energy expression and then use Newtons second law F = ma to calculate the accelerations. We then integrate the accelerations through time to compute the velocities and integrate the velocities to compute the positions of the mass points. The force on the ith mass point is given by the partial derivative of the energy with respect to xi . For the energy of 2, we have if i = 0; nk(x1 x0 ) E if i = n; Fi = = nk(xn xn1 ) (9) xi nk(xi+1 + 2xi xi1 otherwise. The reason for the special cases is that there are two springs pulling on each mass point except for the rst and the last. Note that if the spring is stretched out evenly (xi xj = c(i j)) then the force on all the internal mass points (0 < i < n) is zero because the force from the spring on the left precisely balances the force from the spring on the right. With the above force expressions, we must still choose an integration method for the differential equation F = ma. Let us rst consider the behavior of the simplest integration method, Eulers method. In order to apply Eulers method, we have to transform the differential equation into the
E4
Witkin/Baraff/Kass
canonical form dY/dt = f (Y ). Then if we are given the value of Y at time t, we can compute the value of Y at time t + t from the Euler update formula Y (t + t) = Y (t) + ( t) f (Y ). (10)
The equation F = ma is not in canonical form because the acceleration a is a second derivative and the canonical form requires relationships involving rst derivatives of the state variables. This difculty is easily remedied by introducing new variables to represent the velocities. Then we can write two sets of equations involving rst derivatives instead of one involving second derivatives. Let vi denote the velocity of the ith mass point. Then F = ma in canonical form can be written as follows: dvi dt dxi dt = Fi mi (11) (12) Using this canonical form, the Euler update formula can be written vi (t + t) = vi (t) + ( t)Fi (t)/mi xi (t + t) = xi (t) + ( t)vi (t).
= vi .
(13)
If the time step t is small enough, the above equations provide a means for simulating the behavior of the spring with mass shown in . Unfortunately, the required time step is very small, particularly if n, the number of small springs, is very large. In fact, as we will see shortly, the time step t must be of order 1/n to ensure that the iteration in equation 13 does not diverge. Suppose that the spring is initially at rest with xi = i/n and appropriate forces are applied to both ends in order to keep it at rest. That is, an external force nk(x1 x0 ) is applied to the left end, and an external force nk(xn xn1 ) is applied to the right end. The net force on each mass point is zero, so the system is in equilibrium. Now suppose that a large force P is suddenly applied to the left-most mass point m0. The sudden application of the force will cause a compression wave to move from the left edge towards the right. Physical intuition tells us that the speed of the wave depends on the relationship between the stiffness of the little springs and the mass of the points in between. If the springs are very stiff and the masses are light, we expect the wave to travel very rapidly. Conversely, if the springs are oppy and the masses are heavy, we expect the wave to travel slowly. We will make this notion precise in the next section. For now, we will look at how Eulers method simulates the effect. Suppose the large force P is applied to the left-most mass point at time zero. In the rst Euler step, only the left-most mass point m0 will accelerate. All the other mass points have zero net force, so they remain at rest. At the end of the rst Euler step m0 will have a non-zero velocity, but it still will not have moved. Only after the second Euler step will m0 have moved from its equilibrium position. All the other mass points remain in their original positions. Looking at the force expression in equation 12, we see that the force is still zero except on mass point m1 . The reason is that the force on each mass point depends only on the positions of its immediate neighbors. After two more Euler steps, m1 will have moved and m2 will feel a non-zero force. Only after 2n Euler steps will the rightmost mass point feel any effect from the initial disturbance on the left-most mass point. Let be the time that it should take a disturbance to propagate from one end of the spring to the other in the real physical situation. We know that it has to take at least 2n Euler steps for this to happen in the simulation, so the time step t must be less than /2n in order for the simulation to work.
E5
Witkin/Baraff/Kass
If we try to simulate the behavior with a larger time step, the Euler iteration will simply diverge. Now the problem with Eulers method should be clear. As we increase the number of samples along the spring, we have to keep reducing the time step because information cannot travel faster than one sample per Euler step. Since each Euler step requires a xed amount of work per sample, the amount of computational work needed to simulate the disturbance propagating down the length of the spring is n2 because it requires n iterations each with work proportional to n. If we want to simulate elastic surfaces with thousands of samples, we clearly need a better simulation technique.
6 The Wave Equation

We can gain further insight into the nature of the problem by examining the continuous equations that we are approximating with the collection of springs. Away from the ends of the whole spring, the acceleration Fi /mi of a mass point is given by Fi n(n + 1) = k(xi+1 + 2xi xi1 ) mi m (14)
The right hand side of equation 14 has the familiar look of the second derivative approximation we saw in equation 7. In the limit as n goes to innity, the second derivative returns Fi k 2 x = n m i m u2 lim and we are left with the wave equation: 2 x k 2 x = t2 m u2 (16) (15)
It is also possible to derive equation 16 directly from the energy expression of equation 4 using a mathematical technique known as the calculus of variations. The wave equation has been studied extensively for a long time because of its importance in a number of areas of physics. Its solution is well known. x(u, t) = T1 u t k/m + T2 u + t k/m (17)
The function x(u, t) is given by the sum of a waveform T1 translating in the positive x direction and another waveform T2 translating in the negative x direction. The translation velocity is given by the square root of the constant of proportionality between the second spatial derivative and the second temporal derivative, in this case k/m. If this solution is substituted back into equation 16, it is easy to verify that it satises the differential equation. Now we have another perspective on the difculties with Eulers method. If the spring constant k is very large compared to the mass, then the effect of a force applied at one point in the spring will rapidly spread up and down the length of a spring. A very small time step will be required to combat the problem that Eulers method limits propagation speed to no more than one sample per iteration.
7 Implicit Integration
Not all the difculty in simulating the massive spring is due to Eulers method. Part of it is inherent in the problem statement. Consider for a moment the energy landscape given by equation 4. If the
E6
Witkin/Baraff/Kass
springs are all equally stretched out, the energy is at a local minimum. If any one mass point is moved from this position by itself, the energy increases very rapidly (approximately as the square of the displacement). On the other hand, if all the mass points are moved together by the same amount, the energy does not change at all. While it is somewhat difcult to visualize because it exists in an n dimensional space, this is the description of an energy landscape that acts like a ravine. If you move in most directions, the energy increases rapidly. If you pick a direction near the axis of the ravine, the energy varies very slowly. Note that the sides of the ravine can be made arbitrarily steep by increasing the spring constant k. The correct solution to the differential equation will tend to follow along the ravine in the energy function unless large external forces are pulling it away. The problem with Eulers method is that it tries to follow along the ravine without knowing that a ravine is present. Its local knowledge of the energy function is limited to computing the force at the current state. Since the force is the gradient of the energy, Eulers method computes its next step using a local planar model of the energy function. Unfortunately, for Eulers method, a plane is a very poor model of a ravine, so the method ends up doing a very bad job. The only way to make the planar model a good approximation of a ravine is to consider a very small region. As a result, Eulers method will only work on this problem when very small step sizes are used. To improve the performance, we need to incorporate a better model of the energy function into the differential equation solver. Since rst derivatives of the energy are insufcient, the obvious improvement is to look at second derivatives. In this case, second derivatives of the energy correspond to rst derivatives of the forces.
7.1 General Implicit Formulation

Consider once again the canonical differential equation form dY/dt = f (Y ). Instead of assuming (as Eulers method does) that the derivative dY/dt throughout the time interval is simply f (Y0 ), let us assume that it is some weighted average of the derivative f (Y0 ) at the beginning of the interval and f (Y (t + t)) at the end of the interval. Then we can write the update Y (t + t) = Y (t) + ( t) (1 ) f (Y (t)) + f (Y (t + t)) . (18)
where is a constant between zero and one. Note that when = 0, this reverts to the ordinary Euler update of equation 10. Any update of this form where f (Y (t + t)) appears on the right is known as an implicit update formula. When = 1, this equation is known as the backwards Euler or implicit Euler update. Since we do not know Y (t + t) at the beginning of the step, it may not be clear that equation 18 is helpful in calculating Y (t + t). This is where we make use of the additional derivative information. Around the current state Y0 , we have the approximation dY f (Y0 ) + (Y Y0 )( f )|Y=Y0 dt Substituting this approximation into equation 18 we have Y = ( t) f (Y0 ) + Y ( f )|Y=Y0 Solving for Y, we obtain the update formula Y 1 t I ( f )|Y=Y0 = f (Y0 ) (21) (20) (19)
E7
Witkin/Baraff/Kass
where I denotes the identity matrix. Note that computing the update for Y over a time step now requires solving a linear system. This is the price of using the extra derivative information. For the greatest accuracy, we should use = 1/2 since the best estimate of the derivative in the middle of the interval is halfway between its value at the beginning and the end. There is another criterion that affects our choice, however, and that is the issue of stability. To understand the stability of the update, let us consider the the single variable linear case dY/dt = cY where c is a positive constant and where the initial condition is Y (0) = . Then we have Y (t + t) Y (t) + t[(1 )Y (t) + Y (t + = Y (t) c t[(1 )Y (t) + Y (t + (1 c t(1 ))Y (t) = 1 + c t (1 c t(1 )) t/ t . = 1 + c t t)] t)] (22) (23) (24) (25)
Clearly, the update will diverge if the fraction in equation 25 is greater in magnitude than one. Since c is positive, this can occur only when t > 2/[c(1 2)]. When = 0, we have the ordinary Euler update, and the calculated Y (t) diverges if t > 2/c. In other words, the maximum time step is inversely proportional to the constant of proportionality between Y and Y . As increases, the scheme gets more and more stable. If 1/2, the iteration cannot diverge for this linear differential equation. When applied to non-linear differential equations, however, there is still a possibility of divergence. If is increased from 1/2 to one, the iteration gets even more stable. For severely nonlinear problems, the additional stability, even though it comes at the expense of some reduction in accuracy, can be very useful.
7.2 Application To Springs

When we apply equation 21 to the spring system, there are some simplications because the force depends only on the positions of the mass points and not on their velocities (unless we add damping forces). As a result, we end up with the same equation for the velocities as in the Euler step, but a very different update formula for the positions: ( x) A = v where A= where H= 1 t I H (26)
(27)
2 E xi x j
(28)
is the Hessian matrix of the spring energy. For some problems, this solution method could take even more work than Eulers method because solving a general linear system requires computational work proportional to n3 . In this case, however, the special structure of the problem makes it possible to solve the linear system with computational work proportional to n. As a result, the method is far faster than Eulers method for cases of practical interest.
E8
Witkin/Baraff/Kass
Remember that the force on the ith mass point depends only on the positions of mass points i 1 through i + 1. As a result, the derivative of the force on the ith mass point is zero, except with respect to mass points i 1, i, and i + 1. This means that every row of the matrix A has at most three nonzero elements: the element on the diagonal and one entry on either side. a0 b0 b0 a1 b1 0 .. . b1 a2 .. .. .. A= (29) . . . .. . an3 bn3 0 bn3 an2 bn2 bn2 an1 Matrices with this special form are known as tridiagonal matrices, and linear systems involving them are particularly easy to solve. Numerical Recipes [1], for example, includes a seventeen line program which solves tridiagonal matrices in time proportional to n. The algorithm works by factoring the matrix into the product A = LU of a lower triangular matrix L and an upper-triangular matrix U where r0 s0 c0 d0 c1 r1 s1 0 0 .. d1 c2 . r2 .. .. .. .. . . L= U = . . .. . cn3 rn3 sn3 0 dn3 cn2 0 r s
n2 n2
rn1 (30) It then solves the equation A X = B by rst solving LY = B for Y and then solving Y = U X for X. Consider solving LY = B for Y. The rst component is easy: Y0 = B0 /c0 . After that, we can solve for the remaining components of Y using the recurrence Yi = (Bi Yi1 di1 )/ci . To get an idea of what this recurrence does, let us examine the special case where di = 1/a and ci = (a 1)/a . This is representative of the behavior away from the boundaries. In this case, we can rewrite the recurrence as Yi = Bi + (1 )Yi1 . (31) Some will recognize equation 31 as a simple recursive lter. The next output value Yi is a blend of the last output value Yi1 and the current input value Bi . The effect of any particular input decays exponentially as it gets repeatedly blended in with the subsequent inputs.The result is a low-pass lter with an effective smoothing width determined by . The output of the lter when given an impulse as input is shown in b. When the subsequent linear system Y = U X is solved for X, the result is to do essentially the same thing but in the reverse direction. Figure 2c shows the effect of the two combined ltering operations is shown. The initial spike is blurred out symmetrically in both directions. Note that the behavior will be somewhat different near the boundaries because di and ci will not be constant. Now it should be clear how the implicit solution method of avoids the problems with Eulers method. When a force is applied to a single mass point, the low-pass ltering applies the effect of the force to all the neighboring mass points. If k is large, then the coefcients of the recursive lter
dn2
cn1
E9
Witkin/Baraff/Kass
Figure 2: Effect of recursive lter. Input (a) is ltered forward (b) and and then backward (c). become such that the width of the lter is correspondingly large. In the limit as k grows very large, the lter will essentially move all the mass points together as one. With this solution method, the number of iterations required for a disturbance to propagate from one end of the spring to the other is independent of n. The coefcients of the lters adjust to take account of the different number of samples. Since the computational work per iteration is proportional to n, the entire simulation takes computational work proportional to n instead of the n2 work required by Eulers method. This can be the difference between having a practical simulation and having an impractical one.
8 Fluids
With a very small modication, the spring with mass can be made to simulate the behavior of shallow water. While the derivation of the equations is fairly different, the end result is much the same. If we divide a volume of water into a collection of vertical columns as in gure 3 and make use of a set of approximations known as the shallow water equations, we can describe the motion of the surface by a wave equation in which the wave velocity is proportional to the square root of the depth: 2 h 2 h = gd 2 . t2 x (32)
where h is the height of the column of water, d = h b is its depth, and g is the gravitational constant. In [2], the discretization 2 h t2 = g di1 + di 2( x)2 (hi hi1 ) + g di + di+1 2( x)2 (hi+1 hi ) (33)
is presented. The main difference between these uid equations and the spring equations we have been examining is that the spring constant is no longer uniform. As a result, the differential equation becomes nonlinear and it becomes useful to have = 1.
E10
Witkin/Baraff/Kass
h0
h1
h2
h3
...
hn3
hn2
hn1
b0 u0
b1 u1
b2
b3
...
bn3 un3
bn2
bn1 un2
Figure 3: Discrete two-dimensional height-eld representation of the water surface h, the ground bottom b, and the horizontal water velocity u. The iteration required to solve equation 33 with the implicit method of equation 21 consists of solving a tridiagonal linear system. Precise formulas for the entries in the linear system and many further details can be found in [2].
9 Surfaces
So far, we have been considering the one-dimensional problem of a spring with mass. If we use a rectangular lattice of springs, we can simulate an elastic membrane surface. For the most part, the one-dimensional analysis still applies. We still have difculties with forces propagating over large numbers of samples. On a surface, however, the problem is more severe because the matrix equations we get from the implicit solution technique are no longer tridiagonal. The reason is that mass points are tied to their neighbors along both rows and columns. No matter how the mass points are indexed in the matrix, it will come out with a more complicated structure than the tridiagonal situation in one dimension. One effective solution to this difculty is to consider rows and columns of the spring lattice alternately. Then the one-dimensional technique can be used on each row or column separately. This works very well, for example, in the case of shallow uids.
10 Conclusion
Regular lattices of masses and springs have some special properties that are important to be aware of for the purposes of simulation. If the number of nodes in the lattice is large, then the choice of solution technique is very important. A naive use of Eulers method can result in computational cost
E11
Witkin/Baraff/Kass
proportional to the square of the number of nodes. This is because the effects of forces can never propagate faster than one sample per iteration. Using implicit techniques can avoid these problems and produce a computational cost which is linear in the number of nodes.
References
[1] W. Press, B. Flanner, S. Teukolsky and W. Vetterling,. Numerical Recipes: The Art of Scientic Computing Cambridge University Press, Cambridge 1986. [2] Michael Kass and Gavin Miller. Rapid, Stable Fluid Mechanics for Computer Graphics. In Proceedings of SIGGRAPH 90, pages 4957, August 1990.
E12
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: Constrained Dynamics

Andrew Witkin Robotics Institute Carnegie Mellon University
Please note: This document is 1997 by Andrew Witkin. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Constrained Dynamics
Andrew Witkin School of Computer Science Carnegie Mellon University
Beyond penalty methods
The idea of constrained particle dynamics is that our description of the system includes not only particles and forces, but restrictions on the way the particles are permitted to move. For example, we might constrain a particle to move along a specied curve, or require two particles to remain a specied distance apart. The problem of constrained dynamics is to make the particles obey Newtons laws, and at the same time obey the geometric constraints. As we learned earlier, energy functions provide a sloppy, approximate constraint mechanism. A spring with rest length r makes the particles it connects want to be distance r apart. However, the spring force competes with all other forces acting on the particlesgravity, other springs, forces applied by the user, etc. The constraint can only win this tug-of-war if its spring constant is large enough to overpower all competing inuences, so that very small displacements induce large restoring forces. As we saw in the last section, this is really no solution because it gives rise to stiff differential equations which are all but numerically intractible. The use of extra energy terms to impose constraints is known as the penalty method. If we want both accurate constraints and numerical tractibility, then penalty methods will not ll the bill. Penalty constraints work, to the extent they do, because the restoring forces cancel applied forces that would otherwise break the constraints. The fundamental difculty with penalty constraints is that the applied forces and restoring forces communicate only indirectly, through displacements. In effect, the displacements produced by applied forces act as signals that tell the constraint what restoring force is required. This is not a good communication mechanism because it is impossible to acheive accuracy without stiffness. The basic approach to avoiding this problem is to directly calculate the forces required to maintain the constraints, rather than relying on displacements and restoring forces to do the job. The job of these constraint forces is to cancel just those parts of the applied forces that act against the constraints. Since forces inuence acceleration, this means specically that the constraint forces must convert the particles accelerations into legal accelerations that are consistent with the constraints.
A bead on a wire
We introduce the approach using the simple example of a 2D particle constrained to move on the unit circle. We can express the constraint by writing a scalar behavior function, as we did in Chapter C to create energy functions, 1 (1) C(x) = (x x 1), 2
F1
C=0
C=0
Point-on-circle constraint:
C=
1 2
( xx - 1 )
C=0
C = 0 legal position C = 0 legal velocity C = 0 legal acceleration Add in a constraint force that ensures legal acceleration.
Maintaining Constraints Differentially

Figure 1: If the initial position and velocity are consistent with the constraints, then the constraint can be maintained by ensuring that the acceleration is always legal thereafter. so that the legal positions of x are all those that satisfy C(x) = 0. The function C is an implicit function for the constraint. If x is a legal position, then the legal velocities are all those that satisfy C = x x = 0. In turn, the legal accelerations are all those that satisfy C = x x + x x = 0. (3) (2)
If we start out with a legal position and velocity, then to maintain the constraint, in principal, we need only ensure that equation 3 is satised at every instant thereafter. See gure 1 The particles acceleration is f+ f x= , (4) m where f is the given applied force, and is the as yet unknown constraint force. Substituting for x in f equation 3 gives f+ f C= x + x x = 0, (5) m or x = f x m x x. f (6)
2.1
The principal of virtual work
We have only one equation and two unknownsthe two components of fso we cannot solve for the constraint force without an additional condition. We get that condition by requiring that
F2
Witkin/Baraff/Kass
the constraint force never add energy to nor remove energy from the system, i.e. the constraint is passive and lossless. The kinetic energy is T = and its time derivative is m x x, 2
T = m x x = mf x + m x, f
which is the work done by f and Requiring that the constraint not change the energy means that f. the last term must be zero, i.e. that the constraint force does no work. A subtle point is that we are enjoining the constraint force from ever doing work, rather than saying that it happens not to be at the moment. We therefore require thet x vanish for every legal x, i.e. f x = 0, | x x = 0. f x This condition simply states that x must point in the direction of x, so we can rewrite the constraint force as = x, f where is an unknown scalar. Substituting for in equation 6, and solving for gives f = f x m x x . xx (7)
Having solved for , we calculate = x, then x = (f + f)/m, and proceed with the simulation in f the usual way. In its general form, the condition we imposed on is known as the principal of virtual work.[2] f See gure 2 for an illustration.
2.2
Feedback
If we were solving the differential equation exactly, this procedure would keep the particle exactly on the unit circle, provided we began with valid initial conditions. In practice, we know that numerical solutions to ODEs drift. We must add an extra feedback term to prevent this numerical drift from accumulating, turning the circle into an outward spiral. The feedback term can be just a damped spring force, pulling the particle back onto a unit circle. The feedback force needs to be added in after the constraint force calculation, or else the constraint force will dutifully cancel it! We will discuss feedback in more detail in the next section.
2.3
Geometric Interpretation
= f x x, f xx
When the system is at rest, the constraint force given in equation 7 reduces to (8)
which has clear geometric interpretation as the vector that orthogonally projects f onto the circles tangent. This interpretation makes intuitive sense because the force component that is removed by this projection is the component that points out of the legal motion direction. The orthogonality of the projection also makes sense, because it ensures that the particle will not experience gratuitous accelerations in the allowed direction of motion.
F3
Witkin/Baraff/Kass
fc =
C x
Restrict constraint force to the normal direction. Orthogonal to all legal displacements. No work, no energy gain or loss. One DOF:
Point-on-circle
Constraint Forces
Figure 2: In the case of a point-on-circle constraint, the principle of virtual work simply requires the constraint force to lie in a direction normal to the circle. When x is nonzero, we unfortunately lose this simple geometric picture, but we can still interpret the addition of the constraint force as a projection. Rather than a force projection, it is the projection of the acceleration onto the set of legal accelerations. The velocity-dependent term in equation 7 ensures that the curvature of the particles trajectory matches that of the circle. some effort, we could try to regain the geometric picture by visualizing a projection in phase space, but this is hardly worth the trouble.
The general case
In the last section we derived the constraint force expression for a single particle subject to a single scalar constraint. Our goal in this section is to extend this special case to the general one of a whole system of particles, collectively subjected to a number of constraints. The derivation follows the more detailed one presented in [5]. The key to making this a managable task is to adopt a uniform, monolithic view, much as we do in solving ODEs. Rather than considering each particle separately, we lump their positions into a single state vector, which we will call q. Unlike the phase-space vector that we hand to the solver, this one contains positions only, not velocities,so, in 3D, it has length 3n. To collapse all the particles equations of motion into one global equation, we next dene a mass matrix, M, whose diagonal elements are the particles masses, and whose off-diagonal elements are zero. The diagonal mass matrix for n 3D points is a 3n 3n matrix whose diagonal elements are [m 1 , m 1 , m 1 , m 2 , m 2 , m 2 , . . . , m n , m n , m n ]. implementation, a diagonal matrix may be represented as a vector. Multiplication of a vector by the matrix is then just element-by-element multiplication. The inverse of a diagonal matrix is just the element-by-element reciprocal.
F4
Witkin/Baraff/Kass
Finally, we concatenate the forces on all the particles, just as we do the positions, to create a global force vector, which we denote by Q. Now we can write the global equation governing the particle system as q = WQ, where W is the inverse of M. We will also use global notation for the constraints, concatenating all the scalar constraint functions to form a single vector function C(q). If we have n 3D particles, subject to m scalar constraints, then the output of this global constraint function is an m-vector, and its input is a 3n-vector. At this point, you may well be wondering how all this global notation will ever actually apply to a real network of particles and constraints. This is an example of just the same kind of duality that we encountered in applying ODE solvers to mass-and-spring models. On the one hand, we want to build particle-and-constraint models as networks of distinct interacting objects. On the other, we want to allow the code that calculates constraint forces to act as if the system on which it operates really were a structureless monolith, just as an ODE solver does. Soon, we will show very concretely how this dual view can be maintained. As in the point-on-circle example, we assume that the conguration q and the velocity q are both initially legal, i.e. that C = C = 0. Then our problem is to solve for a constraint force Q that, added to the applied force Q, guarantees that C = 0. To do this, we will need to take derivatives of C. In the previous section, we had a specic algebraic expression for the constraint function, so we were able to derive expressions for their derivatives as well. Now, since we are regarding C as an anonymous function of state, we will be writing derivatives generically, using expressions such as C . To keep things down to earth, you q should think of expressions such as these, as well as C itself, as things we are able to evaluate numerically by invoking functions. 1 By the chain rule, C q. C= q The matrix C/q is called the Jacobian of C. We will denote it henceforward by J. Differentiating again w.r.t. time gives C = Jq + Jq. The quantity J, the time derivative of the Jacobian, might be a little puzzling. By the chain rule we could write it as J = J/qq. However, taking the derivative of a matrix w.r.t. a vector yields a rank 3 tensor (essentially, a 3D array). We can avoid introducing this new kind of object by writing, equivalently J = C/q. Assuming we have an expression for C, this entails only differentiating a vector expression w.r.t. a vector, which is less menacing. Next, we use the systems equations of motion to replace q by a force expression, giving C = Jq + JW(Q + Q). Setting C to zero and re-arranging gives JWQ Jq JWQ,
1
(9)
Speaking of derivatives, it is important to understand what a quantity such as C/q is. Since both C and q are vectors, the derivative of one with respect to the other is a matrix, obtained by taking the scalar derivative of each component of C w.r.t. each component of q.
F5
Witkin/Baraff/Kass
which is the counterpart, in general form, to equation 7. As in the point-on-circle example, we have more unknowns than equations, and once again we introduce the principle of virtual work. The legal velocities (i.e. the ones that dont change C) are all the ones that satisfy J = 0. To ensure x that the constraint force does no work, we therefore require that Q x = 0, | J = 0. x x
All and only vectors Q that satisfy this requirement can be expressed in the form Q = JT , where is a vector with the dimension of C. To understand what this expression means, it helps to regard the matrix J as a collection of vectors, each of which is the gradient of one of the scalar constraint functions comprising C. Since our fundamental requirement is that C = 0, these gradients are normals to the constraint hypersurfaces, representing the state-space directions in which the system is not permitted to move. The vectors that have the form JT are the linear combinations of these gradient vectors, and hence span exactly the set of prohibited directions. Restricting the constraint force to this set ensures that its dot product with any legal displacement of the system will be zero, which is exactly what the principle of virtual work demands. In matrix parlance, the set of vectors JT is known as the null space complement of J. The null space of J is the set of vectors v that satisfy Jv = 0. The null space vectors are the legal displacements, while the null space complement vectors are the prohibited ones. The components of are known as Lagrange multipliers. These quantities, which determine how much of each constraint gradient is mixed into the constraint force, are the unknowns for which we must solve. To do so, we replace Q by JT in equation 9, which gives JWJT = Jq JWQ, (10)
which is a matrix equation in which all but are known. The matrix JWJT is a square matrix with the dimensions of C. Once is obtained, it is multiplied by JT to obtain Q, which is added to the applied force before calculating acceleration. We already noted the need for a feedback term to prevent the accumulation of numerical drift. This term can be incorporated directly into the constraint force calculation. Instead of solving for C = 0, as we did above, we solve for C = ks C kd C, where ks and kd are spring and damping constants. By adding this term, we make the constraint force perform the extra function of returning, with damping, to a valid state should drift occur. The nal constraint force equation, with feedback, is JWJT = Jq JWQ ks C kd C. (11)
The values assigned to ks and kd are not critical, since this term only plays the relatively undemanding role of absorbing drift. See [1, 3, 5] for further discussion.
F6
Witkin/Baraff/Kass
Tinkertoys: Implementing Constrained Particle Dynamics
The general formula of equation 11 is just a skeleton. To actually simulate anything, a specic constraint function C(q) must be provided, in a form that lets us evaluate the function itself and its various derivatives. One way to esh out the skeleton would be to write down an expression for such a function, symbolically take the required derivatives, substitute these expressions into equation 11 and, after simplifying and massaging the resulting mess, write code that performs the numerical evaluations required for a simulation. This is essentially what we did in section 2, but only for a trivially simple example. As an exercise, you might try working through a somewhat more complicated example, say a double pendulum. If you actually do try it, you will probably be able to carry it through to a working implementation, but it will be readily apparent that this derive-and-implement methodology does not scale up! Instead of hand-deriving and hand-coding models, we want to build models interactively by snapping the pieces together, drawing freely from a set of useful pre-dened constraints, such as distance and point-on-curve constraints. The main problem we must solve to acheive this goal is to evaluate the global matrices and vectors that comprise equation 11. The evaluations must be quick, and also dynamic in the sense that we can freely change the structure of the model on the y. Naturally, we want constraints to be modular just as forces are. In this section we describe an architecture for constrained particle dynamics simulation that meets these objectives. The approach is to represent individual constraints using objects similar to those that represent forces. Each constraint object is responsible for evaluating the constraint function it represents, and also that functions derivatives. These evaluations produce fragments of the vectors and matrices comprising equation 11. The fragments are then combined dynamically by the constrained particle system.
4.1
The constrained simulation loop
The machinery required to support constraints ts neatly into our basic particle system architecture. From the standpoint of the ODE solver, the main job of the particle system, in both the unconstrained and constrained case, is to perform derivative evaluations. The sequence of steps that the particle system must perform to evaluate the derivative is nearly the same, with one important addition: calculating the constraint force. This is how the extra step ts in: 1. Clear forces: zero each particles force accumulator. 2. Calculate forces: loop over all force objects, allowing each to add forces to the particles it inuences. 3. Calculate constraint forces: On completion of the previous step, each particles force accumulator contains the total force on that particle. In this step, the global equation 11 is set up and solved, yielding a constraint force on each particle, which is added into the applied force. 4. Calculate the derivative: Divide force by mass to get acceleration, and gather the derivatives into a global vector for the solver. In this section, we are concerned exclusively with the third step in this sequence.
F7
Witkin/Baraff/Kass
C xi xj
Each constraint contributes one or more blocks to the matrix. Sparsity: many empty blocks. Modularity: let each constraint compute its own blocks.
Matrix Block Structure

Figure 3: The block-sparse Jacobian matrix. The constraints are shown above, and the particles along the right side. Each constraint/particle pair corresponds to a block in the matrix. A block is non-zero only if the constraint depends on the particles. Here, we show a binary constraint, such as a distance constraint, that connects two particles. The two shaded blocks are that constraints only contributions to the sparse Jacobian.
4.2
Block-structured matrices
Each individual constraint contributes a slice to the global constraint vector C, just as each particle contributes a slice to the global state vector q. In addition to a list of particles, and a list of forces, a constrained particle system must also maintain a list of constraints. Evaluating the global vectors such as C and C is straightforward, assuming that each constraint points to a function that performs its own portion of these evaluations. We simply loop over the constraints, invoking the functions, and placing the results in the global vectors with appropriate offsets. This is essentially the same gather operation that we use for communication with the ODE solver. The main new ingredient is that we must evaluate global matrices as well as vectors. Whereas constraints and particles each occupy a slice of their respective global vectors, each constraintparticle pair occupies a block of the global derivative matrix. The vector slice that is owned by a constraint can be described by an offset and a length, say i and ilength. Similarly, a particles global station can be described by j and jlength. While jlength is always the dimension of the space that the particles live in, ilength may vary from constraint to constraint. The derivative of a constraint with respect to a particle occupies an ilength jlength block of the Jacobian matrix J, with originn at position (i, j). See gure 3. A typical constraint inuences at most just a few particles. The value of the constraint function depends on only these particles; its derivative with respect to all other particles is zero. This means that the matrices J and J are typically very sparse. Of the n blocks per constraint in an n-particle system, a unary constraint contributes only one non-zero block, a binary one two non-zero blocks,
F8
Witkin/Baraff/Kass
etc. Given this structure, a natural way to represent the sparse matrices is by lists of the non-zero blocks. In this scheme, each matrix block is represented by a structure that species the blocks origin, (i, j), and dimensions, (ilength, jlength), and that contains an ilength jlength oat array holding the blocks data, e.g. struct{ int i; int j; int ilength; int jlength; float *data; }; To support constrained particle dynamics using block-sparse matrices, as we will soon see, we must implement only two operations: matrix times vector, and matrix-transpose times vector. Both are simple operations, looping over the matrix blocks, and performing an ordinary matrix multiplication for each block, using the blocks i and j as offsets into the destination and source vectors.
4.3
Building the matrices
In addition to holding lists of particles and forces, the constrained particle system will hold a list of constraints, block-sparse matrices to represent J and J, and vectors to hold C, C, etc. The structures that represent the constraints may be similar in many respects to the structures we use to represent simple forces, i.e. they point to the particles on which they depend and they point to functions that perform their type-specic operations. When a constraint is instantiated, matrix blocks must be created for each particle on which the constraint depends, and the blocks must be added to the global matrices. Since the number and shape of blocks involved varies with the constraint type, this initialization may be handled by the constraint in a type-specic way. Thereafter, the constraint must be able to evaluate its portions of the vectors C and C, and of the matrices J and J. The results of the matrix evaluations are placed in the matrix blocks that were created by the constraint on initialization. All the required global quantities can then be computed simply by looping over the constraints, and invoking the functions that perform these evaluations.
4.4
Solving the linear system
The solution of sparse linear systems is a eld unto itself. Of the many available options, we give one that is simple and readily available. A matrix equation of the form Mx = b may be solved iteratively by nding a vector x that minimizes (Mx b) (Mx b). A conjugate gradient algorithm that solves this problem is given in Numerical Recipes [4], Chapter 2. The conjugate gradient algorithm offers the advantage that it gives a least-squares solution for over-determined systems, and tolerates redundant constraints. The solver takes as arguments two routines which contitute its only access the matrix: vector-times-matrix, and vector-transpose-times matrix. Sparsity is exploited by implementing these routines efciently. The routine requires O(n) iterations to solve an n n matrix, and the cost of each iteration is O(m), where m is the number of non-zero entries in the matrix.
F9
Witkin/Baraff/Kass
The matrix of equation 11 is JWJT , where J is block-sparse and W is diagonal. We need never actually calculate the matrix. Instead we need only calculate JWJT x, given a vector x. We do this by calculating JT x, using the block-sparse matrix-transpose multiply routine described above, then performing an element-by-element multiplication of the result by the vector representing the diagonal W. Finally, the resulting vector is multiplied by J. Since the compound matrix is symmetric, we do not need a separate function for multiplication by the transpose. Evaluating the right hand side vector of equation 11 is a straightforward application of the block-sparse matrix routines, and standard vector operations. Finally, once the linear system has been solved, the vector is multiplied by JT to produce the global constraint force vector Q, which is then scattered into the particles force accumulator.
4.5
Summary
To introduce constraints into a particle system simulation, we add the additional step of constraint force calculation to the derivative evaluation operation. After the ordinary applied forces have been calculated, but before computing accelerations, we perform the following steps: Loop over the constraints, letting each evaluate its own portion of C, C, J and J. Each constraint points to one or more matrix blocks that receive its contributions to the global matrices. Form the right-hand-side vector of equation 11. Invoke the conjugate gradient solver to obtain the Lagrange multiplier vector, . Multiply by JT to obtain the global constraint force vector, and scatter this force to the particles.
Lagrangian Dynamics: modeling objects other than particles
In the previous sections we have seen how to constrain the behavior of a particle system through the use of constraint forces. Our starting point for the derivation was the use of implicit functions functions of state that are supposed to be zeroto represent the constraints. Each scalar implicit function denes a hypersurface in state space, and the legal states of the system are those that lie on the intersection of all the hypersurfaces. Suppose instead that we represented the constraints using a parametric functiona function q(u), with dim u < dim q, so that q(u) species all and only the legal states. In the case of a unit circle, the parametric function would of course be x = [cos , sin ], leaving as the single degree of freedom. In order to use parametric functions to represent constraints, we need to express the constrained systems equations of motion in terms of the new, constrained degrees of freedom, u, rather than the unconstrained q. These new equations, which we will derive in this section, are known as Lagranges equations of motion.[2]. A clear advantage of the parametric constraint representation is that the extra degrees of freedom are actually removed from the system, rather than being neutralized through the use of constraint forces. This, as one would expect, can lead to better performance. However, Lagrangian dynamics has a very serious drawback: it is often difcult or impossible to nd a parametric function that
F10
Witkin/Baraff/Kass
captures the desired constraints. Moreover, in contrast to the implicit form, there is no automatic way to combine multiple constraints. Lagrangian dynamics is therefore unsuitable as a vehicle for interactive model building. Its important role is as an off-line tool for dening new primitive objects that are more complex than particles. As before, we begin with a collection of particles whose positions are described by a global state vector q, a diagonal mass matrix M, and a global applied force vector Q. We also retain the idea of a constraint force vector Q that satises the principle of virtual work. Now, however, the qs are not independent variables, but are given by a function q(u). Our goal is to solve for u, accounting for the constraint forces. In developing the constraint force formulation we made extensive use of the Jacobian of the implicit constraint function. The Jacobian of the parametric function, q , u has a different meaning but is equally important. By the chain rule, the legal particle velocities are given by q = Ju. J= The principle of virtual work therefore requires that QT Ju = 0, u,
which simply means that JT Q = 0. As before, we can write the unconstrained equations of motion as Mq = Q + Q. Now, however, instead of solving for the constraint force, we can simply make it go away, by multiplying both sides of the equation by JT , giving JT Mq JT Q = 0. (12)
Since q is a function of u, we can now remove q from the expression, leaving u as the unknown. Once again invoking the trusty chain rule, q = Ju + Ju. Substituting this expression into equation 12 gives JT MJu + JT MJu JT Q = 0. (13)
which is a matrix equation to be solved for u. Although you will usually see it expressed in a supercially quite different form, equation 13 is equivalent to the classical Lagrangian equation of motion. As weve expressed it here, its close relation to the constraint force formulation should be strikingly clear.
5.1
Hybrid models
The Lagrangian dynamics formulation is well-suited to creating compound objects off-line, while constraint force methods are well-suited to creating constrained models on the y. In [5] we describe an architecture that combines both methods, allowing constraints to be applied dynamically to complex objects that had been pre-dened using Lagrangian dynamics. In [6], Lagrangian dynamics is used to create simplied non-rigid bodies.
F11
Witkin/Baraff/Kass
References
[1] Ronen Barzel and Alan H. Barr. A modeling system based on dynamic constaints. Computer Graphics, 22:179188, 1988. [2] Herbert Goldstein. Classical Mechanics. Addision Wesley, Reading, MA, 1950. [3] John Platt and Alan Barr. Constraint methods for exible models. Computer Graphics, 22:279 288, 1988. [4] W.H. Press, B.P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, England, 1988. [5] Andrew Witkin, Michael Gleicher, and William Welch. Interactive dynamics. Computer Graphics, 24, 1990. Proc. 1990 Symposium on 3-D Interactive Graphics. [6] Andrew Witkin and William Welch. Fast animation and control of non-rigid structures. Computer Graphics, 24(4):243252, July 1990. Proc. Siggraph 90.
F12
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: Rigid Body Simulation IUnconstrained Rigid Body Dynamics
David Baraff Robotics Institute Carnegie Mellon University
Please note: This document is 1997 by David Baraff. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Rigid Body Simulation

Introduction
This portion of the course notes deals with the problem of rigid body dynamics. To help get you started simulating rigid body motion, weve provided code fragments that implement most of the concepts discussed in these notes. This segment of the course notes is divided into two parts. The rst part covers the motion of rigid bodies that are completely unconstrained in their allowable motion; that is, simulations that arent concerned about collisions between rigid bodies. Given any external forces acting on a rigid body, well show how to simulate the motion of the body in response to these forces. The mathematical derivations in these notes are meant to be fairly informal and intuitive. The second part of the notes tackles the problem of constrained motion that arises when we regard bodies as solid, and need to disallow inter-penetration. We enforce these non-penetration constraints by computing appropriate contact forces between contacting bodies. Given values for these contact forces, simulation proceeds exactly as in the unconstrained case: we simply apply all the forces to the bodies and let the simulation unfold as though the motions of bodies are completely unconstrained. If we have computed the contact forces correctly, the resulting motion of the bodies will be free from inter-penetration. The computation of these contact forces is the most demanding component of the entire simulation process.1
Collision detection (i.e. determining the points of contact between bodies) runs a close second though!
G1
Part I. Unconstrained Rigid Body Dynamics
1 Simulation Basics
This portion of the course notes is geared towards a full implementation of rigid body motion. In this section, well show the basic structure for simulating the motion of a rigid body. In section 2, well dene the terms, concepts, and equations we need to implement a rigid body simulator. Following this, well give some code to actually implement the equations we need. Derivations for some of the concepts and equations we will be using will be left to appendix A. The only thing you need to be familiar with at this point are the basic concepts (but not the numerical details) of solving ordinary differential equations. If youre not familiar with this topic, youre in luck: just turn back to the beginning of these course notes, and read the section on Differential Equation Basics. You also might want to read the next section on Particle Dynamics as well, although were about to repeat some of that material here anyway. Simulating the motion of a rigid body is almost the same as simulating the motion of a particle, so lets start with particle simulation. The way we simulate a particle is as follows. We let a function x(t) denote the particles location in world space (the space all particles or bodies occupy during d simulation) at time t. The function v(t) = x(t) = dt x(t) gives the velocity of the particle at time t. The state of a particle at time t is the particles position and velocity. We generalize this concept by dening a state vector Y(t) for a system: for a single particle, Y(t) = x(t) v(t) . (11)
When were talking about an actual implementation, we have to atten out Y(t) into an array. For a single particle, Y(t) can be described as an array of six numbers: typically, wed let the rst three elements of the array represent x(t), and the last three elements represent v(t). Later, when we talk about state vectors Y(t) that contain matrices as well as vectors, the same sort of operation is done to atten Y(t) into an array. Of course, well also have to reverse this process and turn an array of numbers back into a state vector Y(t). This all comes down to pretty simple bookkeeping though, so henceforth, well assume that we know how to convert any sort of state vector Y(t) to an array (of the appropriate length) and vice versa. (For a simple example involving particles, look through the Particle System Dynamics section of these notes.) For a system with n particles, we enlarge Y(t) to be x1 (t) v1 (t) . Y(t) = . (12) . xn (t) vn (t)
G2
Witkin/Baraff/Kass
where xi (t) and vi (t) are the position and velocity of the ith particle. Working with n particles is no harder than working with one particle, so well let Y(t) be the state vector for a single particle for now (and when we get to it later, a single rigid body). To actually simulate the motion of our particle, we need to know one more thingthe force acting on the particle at time t. Well dene F(t) as the force acting on our particle at time t. The function F(t) is the sum of all the forces acting on the particle: gravity, wind, spring forces, etc. If the particle has mass m, then the change of Y over time is given by d Y(t) = dt x(t) v(t) = v(t) F(t)/m . (13)
Given any value of Y(t), equation (13) describes how Y(t) is instantaneously changing at time t. A simulation starts with some initial conditions for Y(0), (i.e. values for x(0) and v(0)) and then uses a numerical equation solver to track the change or ow of Y over time, for as long as were interested in. If all we want to know is the particles location one second from now, we ask the solver to compute Y(1), assuming that time units are in seconds. If were going to animate the motion of 1 2 the particle, wed want to compute Y( 30 ), Y( 30 ) and so on. The numerical method used by the solver is relatively unimportant with respect to our actual implementation. Lets look at how wed actually interact with a numerical solver, in a C++-like language. Assume we have access to a numerical solver, which well generically write as a function named ode. Typically, ode has the following specication:
typedef void *dydt_funcdouble t, double y void , double ydot ;
odedouble y0 , double yend , int len, double t0, double t1, dydt_func dydt;
We pass an initial state vector to ode as an array y0. The solver ode knows nothing about the inherent structure of y0. Since solvers can handle problems of arbitrary dimension, we also have to pass the length len of y0. (For a system of n particles, wed obviously have len = 6n.) We also pass the solver the starting and ending times of the simulation, t0 and t1. The solvers goal is to compute the state vector at time t1 and return it in the array yend. We also pass a function dydt to ode. Given an array y that encodes a state vector Y(t) and a d time t, dydt must compute and return dt Y(t) in the array ydot. (The reason we must pass t to dydt is that we may have time-varying forces acting in our system. In that case, dydt would have to know what time it is to determine the value of those forces.) In tracing the ow of Y(t) from t0 to t1, the solver ode is allowed to call dydt as often as it likes. Given that we have such a routine ode, the only work we need to do is to code up the routine dydt which well give as a parameter to ode. Simulating rigid bodies follows exactly the same mold as simulating particles. The only differd ence is that the state vector Y(t) for a rigid body holds more information, and the derivative dt Y(t) is a little more complicated. However, well use exactly the same paradigm of tracking the movement of a rigid body using a solver ode, which well supply with a function dydt.
2 Rigid Body Concepts

The goal of this section is to develop an analogue to equation (13), for rigid bodies. The nal differential equation we develop is given in section 2.11. In order to do this though, we need to dene
G3
Witkin/Baraff/Kass
a lot of concepts rst and relations rst. Some of the longer derivations are found in appendix A. In the next section, well show how to write the function dydt needed by the numerical solver ode to d compute the derivative dt Y(t) developed in this section.
2.1 Position and Orientation

The location of a particle in space at time t can be described as a vector x(t), which describes the translation of the particle from the origin. Rigid bodies are more complicated, in that in addition to translating them, we can also rotate them. To locate a rigid body in world space, well use a vector x(t), which describes the translation of the body. We must also describe the rotation of the body, which well do (for now) in terms of a 3 3 rotation matrix R(t). We will call x(t) and R(t) the spatial variables of a rigid body. A rigid body, unlike a particle, occupies a volume of space and has a particular shape. Because a rigid body can undergo only rotation and translation, we dene the shape of a rigid body in terms of a xed and unchanging space called body space. Given a geometric description of the body in body space, we use x(t) and R(t) to transform the body-space description into world space (gure 1). In order to simplify some equations well be using, well require that our description of the rigid body in body space be such that the center of mass of the body lies at the origin, (0, 0, 0). Well dene the center of mass more precisely later, but for now, the center of mass can be thought of as a point in the rigid body that lies at the geometric center of the body. In describing the bodys shape, we require that this geometric center lie at (0, 0, 0) in body space. If we agree that R(t) species a rotation of the body about the center of mass, then a xed vector r in body space will be rotated to the worldspace vector R(t)r at time t. Likewise, if p0 is an arbitrary point on the rigid body, in body space, then the world-space location p(t) of p0 is the result of rst rotating p0 about the origin and then translating it: p(t) = R(t) p0 + x(t). (21)
Since the center of mass of the body lies at the origin, the world-space location of the center of mass is always given directly by x(t). This lets us attach a very physical meaning to x(t) by saying that x(t) is the location of the center of mass in world space at time t. We can also attach a physical meaning to R(t). Consider the x axis in body space i.e. the vector (1, 0, 0). At time t, this vector has direction 1 0 R(t) 0 in world space. If we write out the components of r xx r xy R(t) = r xz then R(t) as r yx rzx r yy rzy , r yz rzz (22)
r xx 1 R(t) 0 = r xy 0 r xz
(23)
G4
Witkin/Baraff/Kass
body space y
world space
y'
x(t)
z'
p0 0 x
y z x x'
p(t)
Figure 1: The center of mass is transformed to the point x(t) in world space, at time t. The xed x, y, and z axes of the body in body space transform to the vectors x = R(t)x, y = R(t)y and z = R(t)z. The xed point p0 in body space is transformed to the point p(t) = R(t) p0 + x(t). which is the rst column of R(t). The physical meaning of R(t) is that R(t)s rst column gives the direction that the rigid bodys x axis points in, when transformed to world space at time t. Similarly, the second and third columns of R(t), r yx rzx r yy rzy and r yz rzz are the directions of the y and z axes of the rigid body in world space at time t (gure 2).
2.2 Linear Velocity

For simplicity, well call x(t) and R(t) the position and orientation of the body at time t. The next thing we need to do is dene how the position and orientation change over time. This means we need expressions for x(t) and R(t). Since x(t) is the position of the center of mass in world space, x(t) is the velocity of the center of mass in world space. Well dene the linear velocity v(t) as this velocity: v(t) = x(t). (24)
If we imagine that the orientation of the body is xed, then the only movement the body can undergo is a pure translation. The quantity v(t) gives the velocity of this translation.
2.3 Angular Velocity

In addition to translating, a rigid body can also spin. Imagine however that we freeze the position of the center of mass in space. Any movement of the points of the body must therefore be due to the
G5
Witkin/Baraff/Kass
world space
y'
R(t) = [ x' y' z' ] z'
y z x x'
Figure 2: Physical interpretation of the orientation matrix R(t). At time t, the columns of R(t) are the world-space directions that the body-space x, y, and z axes transform to. body spinning about some axis that passes through the center of mass. (Otherwise the center of mass would itself be moving). We can describe that spin as a vector (t). The direction of (t) gives the direction of the axis about which the body is spinning (gure 3). The magnitude of (t), |(t)|, tells how fast the body is spinning. |(t)| has the dimensions of revolutions/time; thus, |(t)| relates the angle through which the body will rotate over a given period of time, if the angular velocity remains constant. The quantity (t) is called the angular velocity. d For linear velocity, x(t) and v(t) are related by v(t) = dt x(t). How are R(t) and (t) related? (Clearly, R(t) cannot be (t), since R(t) is a matrix, and (t) is a vector.) To answer this question, lets remind ourselves of the physical meaning of R(t). We know that the columns of R(t) tell us the directions of the transformed x, y and z body axes at time t. That means that the columns of R(t) must describe the velocity with which the x, y, and z axes are being transformed. To discover the relationship between (t) and R(t), lets examine how the change of an arbitrary vector in a rigid body is related to the angular velocity (t). Figure 4 shows a rigid body with angular velocity (t). Consider a vector r(t) at time t specied in world space. Suppose that we consider this vector xed to the body; that is, r(t) moves along with the rigid body through world space. Since r(t) is a direction, it is independent of any translational effects; in particular, r (t) is independent of v(t). To study r (t), we decompose r(t) into vectors a and b, where a is parallel to (t) and b is perpendicular to (t). Suppose the rigid body were to maintain a constant angular velocity, so that the tip of r(t) traces out a circle centered on the (t) axis (gure 4). The radius of this circle is |b|. Since the tip of the vector r(t) is instantaneously moving along this circle, the instantaneous change of r(t) is perpendicular to both b and (t). Since the tip of r(t) is moving in a circle of radius b, the instantaneous velocity of r(t) has magnitude |b||(t)|. Since b and (t) are perpendicular, their cross product has magnitude |(t) b| = |(t)| |b|. (25)
G6
Witkin/Baraff/Kass
(t)
v(t) x(t) y z x
Figure 3: Linear velocity v(t) and angular velocity (t) of a rigid body. Putting this together, we can write r (t) = (t) (b). However, since r(t) = a + b and a is parallel to (t), we have (t) a = 0 and thus r (t) = (t) b = (t) b + (t) a = (t) (b + a). Thus, we can simply express the rate of change of a vector as r (t) = (t) r(t). (27) (26)
Lets put all this together now. At time t, we know that the direction of the x axis of the rigid body in world space is the rst column of R(t), which is r xx r xy . r xz At time t, the derivative of the rst column of R(t) is just the rate of change of this vector: using the cross product rule we just discovered, this change is r xx (t) r xy . r xz The same obviously holds for the other two columns of R(t). This means that we can write r xx r yx rzx R = (t) r xy (t) r yy (t) rzy . r xz r yz rzz
(28)
G7
Witkin/Baraff/Kass
(t)
(t) b
b a r(t) = a + b
Figure 4: The rate of change of a rotating vector. As the tip of r(t) spins about the (t) axis, it traces out a circle of diameter |b|. The speed of the tip of r(t) is |(t)||b|. This is too cumbersome an expression to tote around though. To simplify things, well use the following trick. If a and b are 3-vectors, then a b is the vector a y bz b y az a x bz + b x az . ax b y bx a y Given the vector a, let us dene a to be the matrix 0 az a y az 0 a x . a y a x 0 Then2 a y bz b y az 0 az a y bx a b = a z 0 a x b y = a x bz + b x az = a b. a y a x 0 bz ax b y bx a y Using the notation, we can rewrite R(t) more simply as r xx r yx rzx R(t) = (t) r xy (t) r yy (t) rzy . r xz r yz rzz
2
(29)
(210)
This looks a little too magical at rst. Did someone discover this identity accidentally? Is it a relation that just happens to work? This construct can be derived by considering whats known as innitesimal rotations. The interested reader might wish to read chapter 4.8 of Goldstein[10] for a more complete derivation of the a matrix.
G8
Witkin/Baraff/Kass
By the rules of matrix multiplication, we can factor this into r xx r yx R(t) = (t) r xy r yy r xz r yz
rzx rzy rzz
(211)
which is a matrix-matrix multiplication. But since the matrix on the right is R(t) itself, we get simply that R(t) = (t) R(t). (212)
This, at last, gives us the relation we wanted between R(t) and (t). Note the correspondence = (t) R(t) for the rotation matrix. between r (t) = (t) r(t) for a vector, and R(t)
2.4 Mass of a Body

In order to work out some derivations, well need to (conceptually) perform some integrations over the volume of our rigid body. To make these derivations simpler, were going to temporarily imagine that a rigid body is made up of a large number of small particles. The particles are indexed from 1 to N. The mass of the ith particle is mi , and each particle has a (constant) location r0i in body space. The location of the ith particle in world space at time t, denoted ri (t), is therefore given by the formula ri (t) = R(t)r0i + x(t). The total mass of the body, M, is the sum
N
(213)
M=
i=1
mi .
(214)
(Henceforth, summations are assumed to be summed from 1 to N with index variable i.)
2.5 Velocity of a Particle

The velocity r i (t) of the ith particle is obtained by differentiating equation (213): using the relation R(t) = R(t), we obtain r i (t) = R(t)r0 i + v(t). We can rewrite this as r i (t) = (t) R(t)r0i + v(t) = (t) (R(t)r0i + x(t) x(t)) + v(t) = (t) (ri (t) x(t)) + v(t) using the denition of ri (t) from equation (213). Recall from the denition of the operator that (t) a = (t) a for any vector a. Using this, we can simply write r i (t) = (t) (ri (t) x(t)) + v(t). (217)
(215)
(216)
Note that this separates the velocity of a point on a rigid body into two components (gure 5): a linear component v(t), and an angular component (ri (t) x(t)).
G9
Witkin/Baraff/Kass
(t)
x(t) y z x v(t) + (t) ( ri(t) x(t)) v(t) v(t) ri(t)
(t) ( ri(t) x(t))
Figure 5: The velocity of the ith point of a rigid body in world space. The velocity of ri (t) can be decomposed into a linear term v(t) and an angular term (t) (ri (t) x(t)).
2.6 Center of Mass

Our denition of the center of mass is going to enable us to likewise separate the dynamics of bodies into linear and angular components. The center of mass of a body in world space is dened to be mi ri (t) M (218)
where M is the mass of the body (i.e. the sum of the individual masses mi ). When we say that we are using a center of mass coordinate system, we mean that in body space, 0 mir0i 0 . (219) =0= M 0 Note that this implies that mir0i = 0 as well. We have spoken of x(t) as being the location of the center of mass at time t. Is this true? Yes: since the ith particle has position ri (t) = R(t)r0i + x(t) at time t, the center of mass at time t is miri (t) = M mi (R(t)r0i + x(t)) R(t) = M mir0i + M mi x(t) = x(t) mi = x(t). M
Additionally, the relation mi (ri (t) x(t)) = is also very useful. mi (R(t)r0i + x(t) x(t)) = R(t) mir0i = 0 (220)
G10
Witkin/Baraff/Kass
i(t) = ( ri(t) x(t)) Fi(t)
x(t) y z x
ri(t) Fi(t)
Figure 6: The torque i (t) due to a force Fi (t) acting at ri (t) on a rigid body.
2.7 Force and Torque

When we imagine a force acting on a rigid body due to some external inuence (e.g. gravity, wind, contact forces), we imagine that the force acts on a particular particle of the body. (Remember that our particle model is conceptual only. We can have a force act at any geometrical location on or inside the body, because we can always imagine that there happens to be a particle at that exact location.) The location of the particle the force acts on denes the location at which the force acts. We will let Fi (t) denote the total force from external forces acting on the ith particle at time t. Also, we dene the external torque i (t) acting on the ith particle as i (t) = (ri (t) x(t)) Fi (t). (221)
Torque differs from force in that the torque on a particle depends on the location ri (t) of the particle, relative to the center of mass x(t). We can intuitively think of the direction of i (t) as being the axis the body would spin about due to Fi (t), if the center of mass were held rmly in place (gure 6). The total external force F(t) acting on the body is the sum of the Fi (t): F(t) = while the total external torque is dened similarly as (t) = i (t) = (ri (t) x(t)) Fi (t). (223) Fi (t) (222)
Note that F(t) conveys no information about where the various forces acted on the body; however, (t) does tell us something about the distribution of the forces Fi (t) over the body.
G11
Witkin/Baraff/Kass
2.8 Linear Momentum

The linear momentum p of a particle with mass m and velocity v is dened as p = mv. (224)
The total linear momentum P(t) of a rigid body is the sum of the products of the mass and velocity of each particle: P(t) = mir i (t). (225)
From equation (217), the velocity r i (t) of the ith particle is r i (t) = v(t) + (t) (ri (t) x(t)). Thus, the total linear momentum of the body is P(t) = = = mir i (t) mi v(t) + (t) (ri (t) x(t)) mi v(t) + (t) mi (ri (t) x(t)) . (226)
Because we are using a center of mass coordinate system, we can apply equation (220) and obtain P(t) = mi v(t) = mi v(t) = Mv(t). (227)
This gives us the nice result that the total linear momentum of our rigid body is the same as if the body was simply a particle with mass M and velocity v(t). Because of this, we have a simple transformation between P(t) and v(t): P(t) = Mv(t) and v(t) = P(t)/M. Since M is a constant, v(t) = P(t) . M (228)
The concept of linear momentum lets us express the effect of the total force F(t) on a rigid body quite simply. Appendix A derives the relation P(t) = F(t) (229)
which says that the change in linear momentum is equivalent to the total force acting on a body. Note that P(t) tells us nothing about the rotational velocity of a body, which is good, because F(t) also conveys nothing about the change of rotational velocity of a body! Since the relationship between P(t) and v(t) is simple, we will be using P(t) as a state variable for our rigid body, instead of v(t). We could of course let v(t) be a state variable, and use the relation v(t) = F(t) . M (230)
However, using P(t) instead of v(t) as a state variable will be more consistent with the way we will be dealing with angular velocity and acceleration.
G12
Witkin/Baraff/Kass
2.9 Angular Momentum

While the concept of linear momentum is pretty intuitive (P(t) = Mv(t)), the concept of angular momentum (for a rigid body) is not. The only reason that one even bothers with the angular momentum of a rigid body is that it lets you write simpler equations than you would get if you stuck with angular velocity. With that in mind, its probably best not to worry about attaching an intuitive physical explanation to angular momentumall in all, its a most unintuitive concept. Angular momentum ends up simplifying equations because it is conserved in nature, while angular velocity is not: if you have a body oating through space with no torque acting on it, the bodys angular momentum is constant. This is not true for a bodys angular velocity though: even if the angular momentum of a body is constant, the bodys angular velocity may not be! Consequently, a bodys angular velocity can vary even when no force acts on the body. Because of this, it ends up being simpler to choose angular momentum as a state variable over angular velocity. For linear momentum, we have the relation P(t) = Mv(t). Similarly, we dene the total angular momentum L(t) of a rigid body by the equation L(t) = I(t)(t), where I(t) is a 3 3 matrix (technically a rank-two tensor) called the inertia tensor, which we will describe momentarily. The inertia tensor I(t) describes how the mass in a body is distributed relative to the bodys center of mass. The tensor I(t) depends on the orientation of a body, but does not depend on the bodys translation. Note that for both the angular and the linear case, momentum is a linear function of velocityits just that in the angular case the scaling factor is a matrix, while its simply a scalar in the linear case. Note also that L(t) is independent of any translational effects, while P(t) is independent of any rotational effects. The relationship between L(t) and the total torque (t) is very simple: appendix A derives L(t) = (t), analogous to the relation P(t) = F(t). (231)
2.10 The Inertia Tensor

The inertia tensor I(t) is the scaling factor between angular momentum L(t) and angular velocity (t). At a given time t, let ri be the displacement of the ith particle from x(t) by dening ri = ri (t) x(t). The tensor I(t) is expressed in terms of ri as the symmetric matrix I(t) = mi (ri 2 + ri 2 ) y z miri y ri x miri z ri x mi ri x ri y mi (ri 2 + ri 2 ) x z miriz ri y miri x ri z mi ri y ri z mi (ri 2 + ri 2 ) x y (232)
For an actual implementation, we replace the nite sums with integrals over a bodys volume in world space. The mass terms mi are replaced by a density function. At rst glance, it seems that we would need to evaluate these integrals to nd I(t) whenever the orientation R(t) changes. This would be prohibitively expensive to do during a simulation unless the bodys shape was so simple (for example, a sphere or cube) that that the integrals could be evaluated symbolically. Fortunately, by using body-space coordinates we can cheaply compute the inertia tensor for any orientation R(t) in terms of a precomputed integral in body-space coordinates. (This integral is typically computed before the simulation begins and should be regarded as one of the input parameters
G13
Witkin/Baraff/Kass
describing a physical property of the body.) Using the fact that ri T ri = ri 2 + ri 2 + ri 2 , we can rewrite x y z I(t) as the difference mi ri 2 mi ri x ri y miri x riz 1 0 0 x T m i ri 2 miri y riz (233) I(t) = m i ri ri 0 1 0 m i ri y ri x y 0 0 1 m i ri z ri x m i ri z ri y m i ri 2 z Taking the outer product multiplication of ri with itself, that is ri x T ri ri = ri y ri z ri x ri y ri z = ri y ri x riz ri x ri 2 x ri x ri y ri x ri z ri 2 y riz ri y (234) ri x ri z 2 ri z
and letting 1 denote the 3 3 identity matrix, we can express I(t) simply as I(t) = mi ((ri ri )1 ri ri )
T T
(235)
How does this help? Since ri (t) = R(t)r0i + x(t) where r0i is a constant, ri = R(t)r0i . Then, since R(t)R(t) T = 1, I(t) = = = = mi ((ri ri )1 ri ri ) mi ((R(t)r0 i ) T (R(t)r0 i )1 (R(t)r0 i )(R(t)r0 i ) T ) mi (r0 iT R(t) T R(t)r0i 1 R(t)r0i r0 iT R(t) T ) mi ((r0 iT r0i )1 R(t)r0i r0 iT R(t) T ). (236)
T T
Since r0 iT r0 i is a scalar, we can rearrange things by writing I(t) = = mi ((r0 iT r0i )1 R(t)r0i r0 iT R(t) T ) mi (R(t)(r0 iT r0i )R(t) T 1 R(t)r0i r0 iT R(t) T ) mi ((r0 iT r0 i )1 r0 i r0 iT ) R(t) T . (237)
= R(t) If we dene Ibody as the matrix Ibody =
mi ((r0 iT r0i )1 r0ir0 iT )
(238)
then from the previous equation we have I(t) = R(t)Ibody R(t) T . (239)
Since Ibody is specied in body-space, it is constant over the simulation. Thus, by precomputing Ibody for a body before the simulation begins, we can easily compute I(t) from Ibody and the orientation matrix R(t). Section 5.1 derives the body-space inertia tensor for a rectangular object in terms of an integral over the bodys volume in body space.
G14
Witkin/Baraff/Kass
Also, the inverse of I(t) is given by the formula I 1 (t) = R(t)Ibody R(t) T = R(t) T =
1
1 1 Ibody R(t)1 1 R(t)Ibody R(t) T T
(240)
since, for rotation matrices, R(t) T = R(t)1 and R(t) T during the simulation.
1 = R(t). Clearly, Ibody is also a constant
2.11 Rigid Body Equations of Motion

Finally, we have covered all the concepts we need to dene the state vector Y(t)! For a rigid body, we will dene Y(t) as x(t) R(t) (241) Y(t) = P(t) . L(t) Thus, the state of a rigid body is its position and orientation (describing spatial information), and its linear and angular momentum (describing velocity information). The mass M of the body and bodyspace inertia tensor Ibody are constants, which we assume we know when the simulation begins. At any given time, the auxiliary quantities I(t), (t) and v(t) are computed by v(t) = The derivative
d dt Y(t)
P(t) , M is
I(t) = R(t)Ibody R(t) T
and
(t) = I(t)L(t).
(242)
v(t) x(t) R(t) (t) R(t) d d = Y(t) = F(t) dt dt P(t) L(t) (t)
. (243)
d The next section gives an implementation for the function dydt that computes dt Y(t). One nal note: rather than represent the orientation of the body as a matrix R(t) in Y(t), it is better to use quaternions. Section 4 discusses using quaternions in place of rotation matrices. Briey, a quaternion is a type of four element vector that can be used to represent a rotation. If we replace R(t) in Y(t) with a quaternion q(t), we can treat R(t) as an auxiliary variable that is computed directly from q(t), just as (t) is computed from L(t). Section 4 derives a formula analogous to R(t) = (t) R(t), that expresses q(t) in terms of q(t) and (t).
3 Computing
d dt Y(t)
Lets consider an implementation of the function dydt for rigid bodies. The code is written in C++, and well assume that we have datatypes (classes) called matrix and triple which implement, respectively, 3 3 matrices and points in 3-space. Using these datatypes, well represent a rigid
G15
Witkin/Baraff/Kass
body by the structure

struct RigidBody * Constant quantities * double mass; * mass M * matrix Ibody, * Ibody * 1 Ibodyinv; * Ibody inverse of Ibody * * State variables * triple x; matrix R; triple P, L;
* * * *
x(t) * R(t) * P(t) * L(t) *
* Derived quantities auxiliary variables * matrix Iinv; * I 1 (t) * triple v, * v(t) * omega; * (t) * * Computed quantities * triple force, * F(t) * torque; * (t) * ;
and assume a global array of bodies

RigidBody Bodies NBODIES ;
The constant quantities mass, Ibody and Ibodyinv are assumed to have been calculated for each member of the array Bodies, before simulation begins. Also, the initial conditions for each rigid body are specied by assigning values to the state variables x, R, P and L of each member of Bodies. The implementation in this section represents orientation with a rotation matrix; section 4 describes the changes necessary to represent orientation by a quaternion. We communicate with the differential equation solver ode by passing arrays of real numbers. Several bookkeeping routines are required:
* Copy the state information into an array * void State_to_ArrayRigidBody *rb, double *y *y++ = rb- x 0 ; *y++ = rb- x 1 ; *y++ = rb- x 2 ; * x component of position * * etc. *
forint i = 0; i 3; i++ * copy rotation matrix * forint j = 0; j 3; j++ *y++ = rb- R i,j ;
G16
Witkin/Baraff/Kass
*y++ = rb- P 0 ; *y++ = rb- P 1 ; *y++ = rb- P 2 ; *y++ = rb- L 0 ; *y++ = rb- L 1 ; *y++ = rb- L 2 ;
and
* Copy information from an array into the state variables * void Array_to_StateRigidBody *rb, double *y rb- x 0 rb- x 1 rb- x 2 = *y++; = *y++; = *y++;
forint i = 0; i 3; i++ forint j = 0; j 3; j++ rb- R i,j = *y++; rb- P 0 rb- P 1 rb- P 2 rb- L 0 rb- L 1 rb- L 2 = *y++; = *y++; = *y++; = *y++; = *y++; = *y++;
* Compute auxiliary variables... * * v(t) = M * rb- v = rb- P

P(t)
mass;
1 * I 1 (t) = R(t)Ibody R(t) T * rb- Iinv = R * Ibodyinv * TransposeR;
* (t) = I 1 (t)L(t) * rb- omega = rb- Iinv * rb- L;
Note that Array_to_State is responsible for computing values for the auxiliary variables Iinv, v and omega. Well assume that the appropriate arithmetic operations have been dened between real numbers, triples and matrixs, and that Transpose returns the transpose of a matrix.
G17
Witkin/Baraff/Kass
Examining these routines, we see that each rigid bodys state is represented by 3 + 9 + 3 + 3 = 18 numbers. Transfers between all the members of Bodies and an array y of size 18 NBODIES are implemented as
define STATE_SIZE 18
void Array_to_Bodiesdouble y
forint i = 0; i NBODIES; i++ Array_to_State&Bodies i , &y i * STATE_SIZE ;
and
void Bodies_to_Arraydouble y
forint i = 0; i NBODIES; i++ State_to_Array&Bodies i , &y i * STATE_SIZE ;
Now we can implement dydt. Lets assume that the routine

void Compute_Force_and_Torquedouble t, RigidBody *rb;
computes the force F(t) and torque (t) acting on the rigid body *rb at time t, and stores F(t) and (t) in rb- force and rb- torque respectively. Compute_Force_and_Torque takes into account all forces and torques: gravity, wind, interaction with other bodies etc. Using this routine, well dene dydt as
void dydtdouble t, double y , double ydot *
* put data in y into Bodies Array_to_Bodiesy; forint i = 0; i NBODIES; i++
Compute_Force_and_Torquet, &Bodies i ; ddt_State_to_Array&Bodies i , &ydot i * STATE_SIZE ;
The numerical solver ode calls calls dydt and is responsible for allocating enough space for the arrays y, and ydot (STATE_SIZE NBODIES worth for each). The function which does the real work d of computing dt Y(t) and storing it in the array ydot is ddt_State_to_Array:
G18
Witkin/Baraff/Kass
void ddt_State_to_ArrayRigidBody *rb, double *ydot * copy *ydot++ *ydot++ *ydot++

d dt x(t)
= v(t) into ydot *
= rb- v 0 ; = rb- v 1 ; = rb- v 2 ;
* Compute R(t) = (t) R(t) * matrix Rdot = Starrb- omega * rb- R; * copy R(t) into array * forint i = 0; i 3; i++ forint j = 0; j 3; j++ *ydot++ = Rdot i,j ;
*ydot++ = rb- force 0 ; *ydot++ = rb- force 1 ; *ydot++ = rb- force 2 ; *ydot++ = rb- torque 0 ; *ydot++ = rb- torque 1 ; *ydot++ = rb- torque 2 ;
d * dt P(t) = F(t) *
d * dt L(t) = (t) *
The routine Star, used to calculate R(t) is dened as

matrix Startriple a;
and returns the matrix

a 2 a 1
0 a 2 0
a 0
a 1 a 0
. 0
Given all of the above, actually performing a simulation is simple. Assume that the state variables of all NBODIES rigid bodies are initialized by a routine InitStates. Well have our simulation 1 run for 10 seconds, calling a routine DisplayBodies every 30 th of a second to display the bodies:
void RunSimulation double y0 STATE_SIZE * NBODIES , yfinal STATE_SIZE * NBODIES ;
InitStates; Bodies_to_Arrayyfinal;
G19
Witkin/Baraff/Kass
fordouble t = 0; t
10.0; t += 1. 30.
* copy yfinal back to y0 * forint i = 0; i STATE_SIZE * NBODIES; i++ y0 i = yfinal i ; odey0, yfinal, STATE_SIZE * NBODIES, t, t+1. 30., dydt;
d 1 * copy dt Y(t + 30 ) into state variables *
Array_to_Bodiesyfinal; DisplayBodies;
4 Quaternions vs. Rotation Matrices

There is a better way to represent the orientation of a rigid body than using a 3 3 rotation matrix. For a number of reasons, unit quaternions, a type of four element vector normalized to unit length, are a better choice than for rotation matrices[16]. For rigid body simulation, the most important reason to avoid using rotation matrices is because of numerical drift. Suppose that we keep track of the orientation of a rigid body according to the formula R(t) = (t) R(t). As we update R(t) using this formula (that is, as we integrate this equation), we will inevitably encounter drift. Numerical error will build up in the coefcients of R(t) so that R(t) will no longer be precisely a rotation matrix. Graphically, the effect would be that applying R(t) to a body would cause a skewing effect. This problem can be alleviated by representing rotations with unit quaternions. Since quaternions have only four parameters, there is only one extra variable being used to describe the three freedoms of the rotation. In contrast, a rotation matrix uses nine parameters to describe three degrees of freedom; therefore, the degree of redundancy is noticeably lower for quaternions than rotation matrices. As a result, quaternions experience far less drift than rotation matrices. If it does become necessary to account for drift in a quaternion, it is because the quaternion has lost its unit magnitude.3 This is easily correctable by renormalizing the quaternion to unit length. Because of these two properties, it is desirable to represent the orientation of a body directly as a unit quaternion q(t). We will still express angular velocity as a vector (t). The orientation matrix R(t), which is needed to compute I 1 (t), will be computed as an auxiliary variable from q(t). We will write a quaternion s + v x i + v y j + vz k as the pair [s, v].
Any quaternion of unit length corresponds to a rotation, so quaternions deviate from representing rotations only if they lose their unit length. These notes will deal with that problem in a very simplistic way.
3
G20
Witkin/Baraff/Kass
Using this notation, quaternion multiplication is [s2 , v2 ] = [s1 s2 v1 v2 , s1 v2 + s2 v1 + v1 v2 ]. A rotation of radians about a unit axis u is represented by the unit quaternion [cos(/2), sin(/2)u]. In using quaternions to represent rotations, if q1 and q2 indicate rotations, then q2 q1 represents the composite rotation of q1 followed by q2 .4 In a moment, well show how to change the routines of section 3 to handle the quaternion representation for orientation. Before we can make these changes though, well need a formula for q(t). Appendix B derives the formula q(t) = 1 (t)q(t). 2 (42) (41)
where the multiplication (t)q(t) is a shorthand for multiplication between the quaternions [0, (t)] and q(t). Note the similarity between equation (42) and R(t) = (t) R(t). To actually use a quaternion representation, well need to redene the type RigidBody:
struct RigidBody * Constant quantities * double mass; * mass M * matrix Ibody, * Ibody * 1 Ibodyinv; * Ibody inverse of Ibody * * State variables * triple x; quaternion q; triple P, L;
* * * *
x(t) * q(t) * P(t) * L(t) *
* Derived quantities auxiliary variables * matrix Iinv, * I 1 (t) * R; * R(t) * triple v, * v(t) * omega; * (t) * * Computed quantities * triple force, * F(t) * torque; * (t) * ;
4 This is according to the convention that the rotation of a point p by a quaternion q is qpq1 . Be warned! This is opposite the convention for rotation in the original paper Shoemake[16], but it is in correspondence with some more recent versions of Shoemakes article. Writing a composite rotation as q2 q1 parallels our matrix notation for composition of rotations.
G21
Witkin/Baraff/Kass
Next, in the routine State_to_Array, well replace the double loop
forint i = 0; i 3; i++ * copy rotation matrix * forint j = 0; j 3; j++ *y++ = rb- R i,j ;
with
* Assume that a quaternion is represented in terms of elements `r' for the real part, and `i', `j', and `k' for the vector part. * *y++ *y++ *y++ *y++ = = = = rbrbrbrbq.r; q.i; q.j; q.k;
A similar change is made in Array_to_State. Also, since Array_to_State is responsible for computing the auxiliary variable I 1 (t), which depends on R(t), Array_to_State must also compute R(t) as an auxiliary variable: in the section
* Compute auxiliary variables... * * v(t) = P(t) * M rb- v = rb- P
mass;
1 * I 1 (t) = R(t)Ibody R(t) T * rb- Iinv = R * Ibodyinv * TransposeR;
* (t) = I 1 (t)L(t) * rb- omega = rb- Iinv * rb- L;
we add the line

rb- R = quaternion_to_matrixnormalizerb- q;
prior to computing rb- Iinv. The routine normalize returns q divided by its length; this unit length quaternion returned by normalize is then passed to quaternion_to_matrix which returns a 3 3 rotation matrix. Given a quaternion q = [s, v], quaternion_to_matrix returns the matrix 1 2v2 2v2 2v x v y 2svz 2v x vz + 2sv y y z 2v x v y + 2svz 1 2v2 2v2 2v y vz 2sv x . x z 2v x vz 2sv y 2v y vz + 2sv x 1 2v2 2v2 x y In case you need to convert from a rotation matrix to a quaternion,
G22
Witkin/Baraff/Kass
quaternion matrix_to_quaternionconst matrix &m quaternion double tr = m 0,0 iftr = 0 s = q.r s = q.i q.j q.k else int i = 0; ifm 1,1 m 0,0 i = 1; ifm 2,2 m i,i i = 2; switch i case 0: s = sqrtm 0,0 - m 1,1 + m 2,2 + 1; q.i = 0.5 * s; s = 0.5 s; q.j = m 0,1 + m 1,0 * s; q.k = m 2,0 + m 0,2 * s; q.r = m 2,1 - m 1,2 * s; break; case 1: s = sqrtm 1,1 - m 2,2 + m 0,0 + 1; q.j = 0.5 * s; s = 0.5 s; q.k = m 1,2 + m 2,1 * s; q.i = m 0,1 + m 1,0 * s; q.r = m 0,2 - m 2,0 * s; break; case 2: s = sqrtm 2,2 - m 0,0 + m 1,1 + 1; sqrttr + 1; = 0.5 * s; 0.5 s; = m 2,1 - m 1,2 * s; = m 0,2 - m 2,0 * s; = m 1,0 - m 0,1 * s; q; tr, s; + m 1,1 + m 2,2 ;
G23
Witkin/Baraff/Kass
q.k s = q.i q.j q.r
= 0.5 * s; 0.5 s; = m 2,0 + m 0,2 * s; = m 1,2 + m 2,1 * s; = m 1,0 - m 0,1 * s;
return q;
The matrix m is structured so that m 0,0 , m 0,1 and m 0,2 form the rst row (not column) of m. The routines Array_to_Bodies and Bodies_to_Array dont need any changes at all, but note that the constant STATE_SIZE changes from 18 to 13, since a quaternion requires ve less elements than a rotation matrix. The only other change we need is in ddt_State_to_Array. Instead of
matrix Rdot = Starrb- omega * rb- R;
* copy R(t) into array * forint i = 0; i 3; i++ forint j = 0; j 3; j++ *ydot++ = Rdot i,j ; well use
quaternion qdot = .5 * rb- omega * rb- q; *ydot++ = qdot.r; *ydot++ = qdot.i; *ydot++ = qdot.j; *ydot++ = qdot.k;
Were assuming here that the multiplication between the triple rb- omega and the quaternion . rb- q is dened to return the quaternion product [0, rb- omega]q.
5 Examples
5.1 Inertia Tensor of a Block
Let us calculate the inertia tensor Ibody of the rectangular block in gure 7. The block has dimensions x0 y0 z0 . As required, the center of mass of the block is at the origin. Thus, the extent of the block is from x0 to x0 along the x axis, and similarly for the y and z axes. To calculate the inertia 2 2 tensor, we must treat the sums in equation (232) as integrals over the volume of the block. Let us assume that the block has constant unit density. This means that the density function (x, y, z) is always one. Since the block has volume x0 y0 z0 , the mass M of the block is M = x0 y0 z0 . Then, in
G24
Witkin/Baraff/Kass
x0 , y0 , z0 2 2 2
x0 , y0 , z0 2 2 2
Figure 7: A rectangular block of constant unit density, with center of mass at (0,0,0). body space, Ixx = = = = =
x0 2 x0 2 x0 2 x0 2 x0 2 x0 2 x0 2 x0 2 x0 2 x0 2 y0 2 y0 2 y0 2 y0 2 y0 2 y0 2 z0 2 z0 2
(x, y, z)(y + z ) dx dy dz =
2 2 z= z=
z0 2 z0 2
x0 2 x0 2
y0 2 y0 2
z0 2 z0 2
y2 + z2 dx dy dz
z3 y z+ 3
2
dx dy
y2 z0 +
z3 0 dx dy 12
y= y=
y0 2 y0 2
y3 z3 0 z0 + 0 y0 3 12
(51) dx
y3 z0 z3 y0 0 + 0 dx 12 12
x= x=
x0 2 x0 2
y3 z0 z3 y0 = 0 + 0 12 12 Similarly, I yy = are Ixy =

M 2 12 (x0
x0 2 x0 2 y0 2 y0 2
y3 z0 x0 z3 y0 x0 M x0 y0 z0 2 0 + 0 = (y0 + z2 ) = (y2 + z2 ). 0 0 12 12 12 12 0
M 2 12 (x0
+ z2 ) and Izz = 0
z0 2 z0 2
+ y2 ). Now, the off-diagonal terms, such as Ixy , 0

x0 2 x0 2 y0 2 y0 2 z0 2 z0 2
(x, y, z)(xy) dx dy dz =
xy dx dy dz = 0
(52)
G25
Witkin/Baraff/Kass
y z
x(t) x (3,0,2) (3,0,2)
Figure 8: A block acted on by two equal forces F at two different points. (and similarly for the others) because the integrals are all symmetric. Thus, the inertia tensor of the block is 2 y0 + z2 0 0 0 M . (53) Ibody = 0 x2 + z2 0 0 0 12 2 + y2 0 0 x0 0
5.2 A Uniform Force Field

Suppose a uniform force acts on each particle of a body. For example, we typically describe a gravitational eld as exerting a force mi g on each particle of a rigid body, where g is a vector pointing downwards. The net force Fg acting due to gravity on the body then is Fg = mi g = Mg (54)
which yields an acceleration of Mg = g of the center of mass, as expected. What is the torque due g to the gravitational eld? The net torque is the sum (ri (t) x(t)) mi g = mi (ri (t) x(t)) g = 0 (55)
by equation (220). We see from this that a uniform gravitational eld can have no effect on the angular momentum of a body. Furthermore, the gravitational eld can be treated as a single force Mg acting on the body at its center of mass.
5.3 Rotation Free Movement of a Body

Now, let us consider some forces acting on the block of gure 8. Suppose that an external force F = (0, 0, f ) acts on the body at points x(t) + (3, 0, 2) and x(t) + (3, 0, 2). We would expect
G26
Witkin/Baraff/Kass
that this would cause the body to accelerate linearly, without accelerating angularly. The net force acting on the body is (0, 0, 2 f ), so the acceleration of the center of mass is 2f M along the z axis. The torque due to the force acting at x(t) + (3, 0, 2) is 3 3 ((x(t) + 0 ) x(t)) F = 0 F 2 2 while the torque due to the force acting at x(t) + (3, 0, 2) is 3 3 ((x(t) + 0 ) x(t)) F = 0 F. 2 2 The total torque is therefore 0 3 3 3 3 = 0 F + 0 F = ( 0 + 0 ) F = 0 F. 2 2 2 2 2 But this gives 0 0 = 0 0 = 0. f 2 As expected then, the forces acting on the block impart no angular acceleration to the block.
5.4 Translation Free Movement of a Body

Suppose now that an external force F1 = (0, 0, f ) acts on the body at point x(t) + (3, 0, 2) and an external force F2 = (0, 0, f ) acts on the body at point x(t) + (3, 0, 2) (gure 9). Since F1 = F2 , the net force acting on the block is F1 + F2 = 0, so there is no acceleration of the center of mass. On the other hand, the net torque is
G27
Witkin/Baraff/Kass
y z
x(t) x (3,0,2) F (3,0,2)
Figure 9: A block acted on by two opposite forces F1 and F2 = F1 , at two different points.
3 ((x(t) + 0 ) x(t)) F1 + 2 3 3 ((x(t) + 0 ) x(t)) F2 = 0 2 2 0 = 3f + 0
0 3 0 0 + 0 0 f 2 f 0 0 3f = 6f . 0 0
(56)
Thus, the net torque is (0, 6 f, 0), which is parallel to the y axis. The nal result is that the forces acting on the block cause it to angularly accelerate about the y axis.
5.5 Force vs. Torque Puzzle

In considering the effect of a force acting at a point on a body, it sometimes seems that the force is being considered twice. That is, if a force F acts on a body at a point r + x(t) in space, then we rst consider F as accelerating the center of mass, and then consider F as imparting a spin to the body. This gives rise to what at rst seems a paradox: Consider the long horizontal block of gure 10 which is initially at rest. Suppose that a force F acts on the block at the center of mass for some period of time, say, ten seconds. Since the force acts at the center of mass, no torque is exerted on the body. After ten seconds, the body will have acquired some linear velocity v. The body will not have acquired any angular velocity; thus the kinetic energy of the block will be 1 M|v|2 . 2 Now suppose that the same force F is applied off-center to the body as shown in gure 11. Since the force acting on the body is the same, the acceleration of the center of mass is the same. Thus,
G28
Witkin/Baraff/Kass
F (ten seconds later)
F v
Energy: 0
Energy:
1 Mv T v 2
Figure 10: A rectangular block acted on by a force through its center of mass.
(ten seconds later) v F F Energy: 0 Energy:
1 1 Mv T v + T I 2 2
Figure 11: A block acted on by a force, off-center of the center of mass.
G29
Witkin/Baraff/Kass
after ten seconds, the body will again have linear velocity v. However, after ten seconds, the body will have picked up some angular velocity , since the force F, acting off center, now exerts a torque on the body. Since the kinetic energy is (see appendix C)
2 1 2 M|v|
+ 1 T I 2
the kinetic energy of the block is higher than when the force acted through the center of mass. But if identical forces pushed the block in both cases, how can the energy of the block be different? Hint: Energy, or work, is the integral of force over distance.
G30
Witkin/Baraff/Kass

(ten seconds later) v F F Energy: 0 Energy:
1 1 Mv T v + T I 2 2
Figure 12: The path the force acts over is longer than in gure 10. As a result, the force does more work, imparting a larger kinetic energy to the block. Figure 12 shows why the force acting off center results in a higher kinetic energy. The kinetic energy of the block is equivalent to the work done by the force. The work done by the force is the integral of the force over the path traveled in applying that force. In gure 11, where the force acts off the center of mass, consider the path traced out by the point where the force is applied. This path is clearly longer than the path taken by the center of mass in gure 10. Thus, when the force is applied off center, more work is done because the point p at which the force is applied traces out a longer path then when the force is applied at the center of mass.
G31
Witkin/Baraff/Kass
An Introduction to Physically Based Modeling: Rigid Body Simulation IINonpenetration Constraints

Please note: This document is 1997 by David Baraff. This chapter may be freely duplicated and distributed so long as no consideration is received in return, and this copyright notice remains intact.
Part II. Nonpenetration Constraints
6 Problems of Nonpenetration Constraints

Now that we know how to write and implement the equations of motion for a rigid body, lets consider the problem of preventing bodies from inter-penetrating as they move about an environment. For simplicity, suppose we simulate dropping a point mass (i.e. a single particle) onto a xed oor. There are several issues involved here. Because we are dealing with rigid bodies, that are totally non-exible, we dont want to allow any inter-penetration at all when the particle strikes the oor. (If we considered our oor to be exible, we might allow the particle to inter-penetrate some small distance, and view that as the oor actually deforming near where the particle impacted. But we dont consider the oor to be exible, so we dont want any inter-penetration at all.) This means that at the instant that the particle actually comes into contact with the oor, what we would like is to abruptly change the velocity of the particle. This is quite different from the approach taken for exible bodies. For a exible body, say a rubber ball, we might consider the collision as occurring gradually. That is, over some fairly small, but non-zero span of time, a force would act between the ball and the oor and change the balls velocity. During this time span, the ball would deform, due to the force. The more rigid we made the ball, the less the ball would deform, and the faster this collision would occur. In the limiting case, the ball is innitely rigid, and cant deform at all. Unless the balls downward velocity is halted instantaneously, the ball will inter-penetrate the oor somewhat. In rigid body dynamics then, we consider collisions as occurring instantaneously. This means we have two types of contact we need to deal with. When two bodies are in contact at some point p, and they have a velocity towards each other (as in the particle striking the oor), we call this colliding contact. Colliding contact requires an instantaneous change in velocity. Whenever a collision occurs, the state of a body, which describes both position, and velocity, undergoes a discontinuity in the velocity. The numerical routines that solve ODEs do so under the assumption that the state Y (t) always varies smoothly. Clearly, requiring Y (t) to change discontinuously when a collision occurs violates that assumption. We get around this problem as follows. If a collision occurs at time tc , we tell the ODE solver to stop. We then take the state at this time, Y (tc ), and compute how the velocities of bodies involved in the collision must change. Well call the state reecting these new velocities Y (tc )+ . Note that Y (tc ) and Y (tc )+ agree for all spatial variables (position and orientation), but will be different for the velocity variables of bodies involved in the collision at time tc . We then restart the numerical solver, with the new state Y (tc ), and instruct it to simulate forward from time tc . Whenever bodies are resting on one another at some point p (e.g. imagine the particle in contact with the oor with zero velocity), we say that the bodies are in resting contact. In this case, we compute a force that prevents the particle from accelerating downwards; essentially, this force is the weight of the particle due to gravity (or whatever other external forces push on the particle). We call the force between the particle and the oor a contact force. Resting contact clearly doesnt require us to stop and restart the ODE solve at every instant; from the ODE solvers point of view, contact forces are just a part of the force returned by Compute_Force_and_Torque.
G32
Witkin/Baraff/Kass
t0
tc
t0 + t (inter-penetration detected)
Figure 13: At time t0 + t, the particle is found to lie below the oor. Thus, the actual time of collision tc lies between the time of the last known legal position, t0 , and t0 + t. So far then, we have two problems well need to deal with: computing velocity changes for colliding contact, and computing the contact forces that prevent inter-penetration. But before we can tackle these problems we have to deal with the geometric issue of actually detecting contact between bodies. Lets go back to dropping the particle to the oor. As we run our simulation, we compute the position of the particle as it drops towards the oor at specic time values (gure 13). Suppose we consider the particle at times t0 , t0 + t, t0 + 2 t etc.5 and suppose the time of collision, tc , at which the particle actually strikes the oor, lies between t0 and t0 + t. Then at time t0 , we nd that the particle lies above the oor, but at the next time step, t0 + t, we nd the particle is beneath the oor, which means that inter-penetration has occurred. If were going to stop and restart the simulator at time tc , well need to compute tc . All we know so far is that tc lies between t0 and t0 + t. In general, solving for tc exactly is difcult, so we solve for tc numerically, to within a certain tolerance. A simple way of determining tc is to use a numerical method called bisection[14]. If at time t0 + t we detect inter-penetration, we inform the ODE solver that we wish to restart back at time t0 , and simulate forward to time t0 + t/2. If the simulator reaches t0 + t/2 without encountering inter-penetration, we know the collision time tc lies between t0 + t/2 and t0 + t. Otherwise, tc is less than t0 + t/2, and we try to simulate from t0 to t0 + t/4. Eventually, the time of collision tc is computed to within some suitable numerical tolerance. The accuracy with which tc is found depends on the collision detection routines. The collision detection routines have some parameter . We decide that our computation of tc is good enough when the particle inter-penetrates the oor by no more than , and is less than above the oor. At this point we declare that the particle is in contact with the oor (gure 14). The method of bisection is a little slow, but its easy to implement and quite robust. A faster method involves actually predicting the time tc of the collision, based on examining Y(t0 ) and Y(t0 + t). Baraff[1, 2] describes how to make such predictions. How to actually implement all of
5
The ODE solver doesnt have to proceed with equal size time steps though.
G33
Witkin/Baraff/Kass
t0 tc found (within tolerance)

t0 + t (inter-penetration detected)
Figure 14: When the particle is found to be within some tolerance of contacting the oor, then tc is considered to have been computed to within sufcient accuracy. this depends on how you interact with your ODE routines. One might use exception handling code to signal the ODE of various events (collisions, inter-penetration), or pass some sort of messages to the ODE solver. Well just assume that you have some way of getting your ODE solver to progress just up to the point tc . Once you actually reach the time of a collision, or whenever youre in a state Y (t) where no inter-penetration has occurred, a geometric determination has to be made to nd all the points of contact. (Just because you may be looking for the time of collision between two bodies A and B doesnt mean you get to neglect resting contact forces between other bodies C and D. Whenever youre trying to move the simulation forward, youll need to compute the point of contact between bodies and the contact forces at those points.) There is a vast amount of literature dealing with the collision detection problem. For instance, some recent SIGGRAPH papers dealing with the subject are Von Herzen, Barr and Zatz[17] and Moore and Wilhelms[12]; in robotics, a number of papers of interest are Canny[4], Gilbert and Hong[6], Meyer[11] and Cundall[5]. Preparata and Shamos[13] describes many approaches in computational geometry to the problem. In the next section, well briey describe a collision detection philosophy that leads to very efcient algorithms, for the sorts of simulation these course notes are concerned with. Actual code for the algorithms is fairly easy to write, but a little too lengthy to t in these notes. Following this, well move on to consider colliding and resting contact.
7 Collision Detection
The collision detection algorithm begins with a preprocessing step, in which a bounding box for each rigid body is computed (a box with sides parallel to the coordinate axes). Given n such bounding boxes, we will want to quickly determine all pairs of bounding boxes that overlap. Any pair of rigid bodies whose bounding boxes do not overlap need not be considered any further. Pairs of rigid
G34
Witkin/Baraff/Kass
bodies whose bounding boxes do overlap require further consideration. Well rst describe how to efciently check for inter-penetration or contact points between rigid bodies dened as convex polyhedra. Then well show how to perform the bounding box check efciently. As described in section 1, the simulation process consists of the repeated computation of the d derivative of the state vector, dt Y(t), at various times t. The numerical ODE solver is responsible for choosing the values of t at which the state derivative is to be computed. For any reasonably complicated simulation, the values of t chosen are such that the state Y does not change greatly between successive values of t. As a result, there is almost always great geometric coherence between successive time steps. At a time step t0 + t, the idea is to take advantage of the collision detection results computed at the previous time step t0 .
7.1 Convex Polyhedra

Our primary mechanism for exploiting coherence will be through the use of witnesses. In our context, given two convex polyhedra A and B, a witness is some piece of information that can be used to quickly answer the yes/no question are A and B disjoint? We will utilize coherence by caching witnesses from one time step to the next; hopefully a witness from the previous time step will be a witness during the current time step. Since we are considering convex polyhedra, two polyhedra do not inter-penetrate if and only if a separating plane between them exists. A separating plane between two polyhedra is a plane such that each polyhedron lies on a different side of the plane. A given plane can be veried to be a separating plane by testing to make sure that all of the vertices of A and B lie on opposite sides of the plane. Thus, a separating plane is a witness to the fact that two convex polyhedra do not inter-penetrate. If a separating plane does not exist, then the polyhedra must be inter-penetrating. The cost of initially nding a witness (for the very rst time step of the simulation, or the rst time two bodies become close enough to require more than a bounding box test) is unavoidable. A simple way to nd a separating plane initially is as follows. If a pair of convex polyhedra are disjoint or contacting (but not inter-penetrating), then a separating plane exists with the following property: either the plane contains a face of one of the polyhedra, or the plane contains an edge from one of the polyhedra and is parallel to an edge of the other polyhedra. (That is, the separating planes normal is the cross product of the two edge directions, and the plane itself contains one of the edges.) We will call the face or edges in question the dening face or edges. Initially, we simply check all possible combinations of faces and edges to see if one such combination forms a separating plane (gure 15). Although this is inefcient, its done so infrequently that the inefciency is unimportant. For subsequent time steps, all we need to do is form a separating plane from the dening face or edges found during the previous time step, and then verify the plane to see that it is still valid (gure 16). On those (rare) occasions when the cached face or two edges fails to form a valid separating plane (gure 17), faces or edges adjacent to the previously cached face or edges can be examined to see if they form a separating plane; however, this happens infrequently enough that it may be simpler to start from scratch and compute a new separating plane without using any prior knowledge. Once the separating place has been found, the contact region between the two polyhedra is determined, assuming the polyhedra are not disjoint. Contact points between the two polyhedra can only occur on the separating plane. Given the separating plane, the contact points can be quickly and efciently determined by comparing only those faces, edges, and vertices of the polyhedra that are coincident with the separating plane. However, if no separating plane can be found, then the two polyhedra must be inter-penetrating. When two polyhedra inter-penetrate, it is almost always the case that either a vertex of one poly-
G35
Witkin/Baraff/Kass
separating plane
Figure 15: Exhaustive search for a separating plane. Only one face of the two polygons forms a separating plane.
(a)
(b)
defining face
Figure 16: (a) At this time step, the separating plane is dened by a face of one of the polygons. (b) At the next time step, the polygons have moved, but the same face still denes a separating plane.
G36
Witkin/Baraff/Kass
(a)
(b)
Figure 17: The face that has been dening a separating plane no longer does so, and a new separating plane must be found. hedron is inside the other, or an edge of one polyhedron has intersected a face of the other.6 In this case, the inter-penetrating vertex, or intersecting edge and face are cached as a witness to the inter-penetration. Since this indicates a collision at some earlier time, the simulator will back up and d attempt to compute dt Y(t) at some earlier time. Until the collision time is determined, the rst action taken by the collision/contact determination step will be to check the cached vertex or edge and face to see if they indicate inter-penetration. Thus, until the collision time is found, states in which the inter-penetration still exists are identied as such with a minimum of computational overhead.
7.2 Bounding Boxes

To reduce the number of pairwise collision/contact determinations necessary, a bounding box hierarchy is imposed on the bodies in the simulation environment. If two bounding boxes are found not to overlap, no further comparisons involving the contents of the boxes are needed. Given a collection of n rectangular bounding boxes, aligned with the coordinate axes, we would like to efciently determine all pairs of boxes that overlap. A naive pairwise comparison of all pairs requires O(n2 ) work and is too inefcient, unless the number of bodies is small. Computational geometry algorithms exist that can solve this problem in time O(nlogn + k) where k is the number of pairwise overlaps; a general result is that the problem can be solved in time O(nlogd2 n + k) for d-dimensional bounding boxes[13]. Using coherence, we can achieve substantially better performance.
6 An exception is the following. Stack two cubes of equal size atop one another so that their contacting faces exactly coincide. Lower the top one. This produces an inter-penetration such that no vertex is inside either cube, and no edge penetrates through any face.
G37
Witkin/Baraff/Kass
(a)
(b)
I6 I3 I1
I5 I4 I2
b3 b6 b1
e6
e1 b5 b2 e3
b4 e5
e4
e2
Figure 18: The sweep/sort algorithm. (a) When b1 is encountered, the active list contains intervals 3 and 6; interval 1 is reported to overlap with these two intervals. Interval 1 is added to the active list and the algorithm continues. (b) When e3 is encountered, the active list contains intervals 2, 3 and 5. Interval 3 is removed from the active list. 7.2.1 The one-dimensional case
Consider the problem of detecting overlap between one-dimensional bounding boxes, aligned with the coordinate system. Such a bounding box can be described simply as an interval [b, e] where b and e are real numbers. Let us consider a list of n such intervals, with the ith interval being [bi , ei ]. The problem is then dened to be the determination of all pairs i and j such that the intervals [bi , ei ] and [b j , e j ] intersect. The problem can be solved initially by a sort and sweep algorithm. A sorted list of all the bi and ei values is created, from lowest to highest. The list is then swept, and a list of active intervals, initially empty, is maintained. Whenever some value bi is encountered, all intervals on the active list are output as overlapping with interval i, and interval i is then added to the list (gure 18a). Whenever some value ei is encountered , interval i is removed from the active list (gure 18b). The cost of this process is O(n log n) to create the sorted list, O(n) to sweep through the list, and O(k) to output each overlap. This gives a total cost of O(n log n + k), and is an optimal algorithm for initially solving the problem. Subsequent comparisons can be improved as follows. First, there is no need to use an O(n log n) algorithm to form the sorted list of bi and ei values. It is considerably more efcient to start with the order found for bi and ei values from the previous time step; if coherence is high, this ordering will be nearly correct for the current time step. A sorting method called an insertion sort[15] is used to permute the nearly sorted list into a sorted list. The insertion sort algorithm works by moving items towards the beginning of the list, until a smaller item is encountered. Thus, the second item is interchanged with the rst if necessary, then the third item is moved towards the beginning of the list until its proper place is found, and so on; each movement of an item indicates a change in the ordering of two values. After the last item on the list has been processed, the list is in order. Such
G38
Witkin/Baraff/Kass
I6 I3 I1
I5 I4 I2
b3 b6 b1
e6
e1 b5 b2 e3
e5
b4
e4
e2
Figure 19: A coherence-based method of detecting overlaps. The order produced in gure 18 is nearly correct for this arrangement of intervals. Only b4 and e5 need to be exchanged. When the exchange occurs, the change in overlap status between interval 4 and 5 is detected. a sort takes time O(n + c) where c is the number of exchanges necessary. For example, the only difference between gures 19 and 18 is that interval 4 has moved to the right. Starting from the ordered list of bi and ei values of gure 18, only a single exchange is necessary to sort the list for gure 19. The insertion sort is not recommendeded as a sorting procedure in general, since it may require O(n2 ) exchanges; however, it is a good algorithm for sorting a nearly sorted list, which is what occurs in our highly coherent environment. To complete the algorithm, note that if two intervals i and j overlap at the previous time step, but not at the current time step, one or more exchanges involving either a bi or ei value and a b j or e j value must occur. The converse is true as well when intervals i and j change from not overlapping at the previous time step to overlapping at the current time step. Thus, if we maintain a table of overlapping intervals at each time step, the table can be updated at each time step with a total cost of O(n + c). Assuming coherence, the number of exchanges c necessary will be close to the actual number k of changes in overlap status, and the extra O(c k) work will be negligible. Thus, for the one-dimensional bounding box problem, the coherence view yields an efcient algorithm of extreme (if not maximal) simplicity that approaches optimality as coherence increases. 7.2.2 The three-dimensional case
Efcient computational geometry algorithms for solving the bounding box intersection problem in IR3 are much more complicated than the sort and sweep method for the one-dimensional case. However, these algorithms all have in common a step that is essentially a sort along a coordinate axis, as in the one-dimensional case. Each bounding box is described as three independent intervals [bi(x) , ei(x)], [bi(y) , ei(y) ], and [bi(z) , ei(z)] which represent the intervals spanned on the three coordinate axes by the ith bounding box. Thus, our rst thought towards improving the efciency of a computational geometry algorithm for coherent situations would be to sort a list containing the bi(x) and ei(x) values, and similarly for the y and z axes. Again, such a step will involve O(n + c) work, where c is now the
G39
Witkin/Baraff/Kass
total number of exchanges involved in sorting all three lists. However, if we observe that checking two bounding boxes for overlap is a constant time operation, it follows that if we simply check bounding boxes i and j for overlap whenever an exchange is made between values indexed by i and j (on any coordinate axis), we will detect all changes in overlap status in O(n + c) time. Again, we can maintain a table of overlapping bounding boxes, and update it at each time step in O(n + c) time. The extra work involved is again O(c k). For the three-dimensional case, extra work can occur if the extents of two bounding boxes change on one coordinate axis without an actual change of their overlap status. In practice, the extra work done has been found to be completely negligible, and the algorithm runs essentially in time O(n + k).
8 Colliding Contact
For the remainder of these notes, were going to be concerned with examining the bodies in our simulator at a particular instant of time t0 . At this time t0 , we assume that no bodies are interpenetrating, and that the simulator has already determined which bodies contact, and at which points. To simplify matters, well imagine that all bodies are polyhedra, and that every contact point between bodies has been detected. Well consider contacts between polyhedra as either vertex/face contacts or edge/edge contacts. A vertex/face contact occurs when a vertex on one polyhedra is in contact with a face on another polyhedra. An edge/edge contact occurs when a pair of edges contact; it is assumed in this case that the two edges are not collinear. (Vertex/vertex and vertex/edge contacts are degenerate, and are not considered in these notes.) As examples, a cube resting on a plane would be described as four vertex/face contacts, one contact at each corner of the cube. A cube resting on a table, but with its bottom face hanging over the edge of the table would still be described as four contacts; two vertex/face contacts for the vertices on the table, and two edge/edge contacts, one on each edge of the cube that crosses over an edge of the table. Each contact is represented by a structure
struct Contact RigidBody triple
*a, *b; p, n, ea, eb; vf;
* * * * * *
body containing vertex * body containing face * world-space vertex location * outwards pointing normal of face * edge direction for A * edge direction for B *
boolean ; int Ncontacts; Contact *Contacts;
* TRUE if vertex face contact *
If the contact is a vertex/face contact, then the variable a points to the rigid body that the contact vertex is attached to, while b points to the body the face is attached to. Well call these two bodies A and B respectively. For vertex/face contacts, the variable n is set to the outwards pointing unit normal of the contact face of body B, and the variables ea and eb are unused. For edge/edge contacts, ea is a triple of unit length, that points in the direction of the contacting
G40
Witkin/Baraff/Kass
A pa ( t) pa ( t0 ) pb ( t) pb ( t0 )
Figure 20: (a) The points pa (t) and pb (t) for a vertex/face contact. (b) At time t0 , the bodies come into contact at pa (t0 ) = pb (t0 ). edge of body A (pointed to by a). Similarly, eb is a unit vector giving the direction that the contact edge on body B points. For edge/edge contacts, n denotes a unit vector in the ea eb direction. Well adopt the convention that the two contacting bodies are labeled A and B such that the normal direction ea eb points outwards from B, towards A, as it does for vertex/face contacts. For both types of contact, the position of the contact in world space (which is either the contact vertex, or the point where the two edges intersect) is given by p. The collision detection routines are responsible for discovering all the contact points, setting Ncontacts to the number of contact points, and allocating space for and initializing an array of Contact structures. The rst thing well need to do is examine the data in each Contact structure to see if colliding contact is taking place. For a given contact point, the two bodies A and B contact at the point p. Let pa (t) denote the particular the point on body A that satises pa (t0 ) = p. (For vertex/face contacts, this point will be the vertex itself. For edge/edge contacts, it is some particular point on the contact edge of A.) Similarly, let pb (t) denote the particular point on body B that coincides with pa (t0 ) = p at time t0 (gure 20). Although pa (t) and pb (t) are coincident at time t0 , the velocity of the two points at time t0 may be quite different. We will examine this velocity to see if the bodies are colliding or not. From section 2.5, we can calculate the velocity of the vertex point, pa (t0 ) by the formula pa (t0 ) = va (t0 ) + a (t0 ) ( pa (t0 ) xa (t0 )) (81)
where va (t) and a (t) are the velocities for body A. Similarly, the velocity of the contact point on the face of B is pb (t0 ) = vb (t0 ) + b (t0 ) ( pb (t0 ) xb (t0 )). (82)
G41
Witkin/Baraff/Kass
a ( t0 ) pb ( t 0 ) p A n
Figure 21: The vector pa (t0 ) pb (t0 ) points in the same direction as n(t0 ); the bodies are separating. Lets examine the quantity vrel = n(t0 ) ( pa (t0 ) pb (t0 )), (83)
which is a scalar. In this equation, n(t0 ) is the unit surface normal, described by the variable n, for each contact point. The quantity vrel gives the component of the relative velocity pa (t0 ) pb (t0 ) in the n(t0 ) direction. Clearly, if vrel is positive, then the relative velocity pa (t0 ) pb (t0 ) at the contact point is in the positive n(t0 ) direction. This means that the bodies are moving apart, and that this contact point will disappear immediately after time t0 (gure 21). We dont need to worry about this case. If vrel is zero, then the bodies are neither approaching nor receding at p (gure 22). This is exactly what we mean by resting contact, and well deal with it in the next section. In this section, were interested in the last possibility, which is vrel < 0. This means that the relative velocity at p is opposite n(t0 ), and we have colliding contact. If the velocities of the bodies dont immediately undergo a change, inter-penetration will result (gure 23). How do we compute the change in velocity? Any force we might imagine acting at p, no matter how strong, would require at least a small amount of time to completely halt the relative motion between the bodies. (No matter how strong your car brakes are, you still need to apply them before you hit the brick wall. If you wait until youve contacted the wall, its too late...) Since we want bodies to change their velocity instantly though, we postulate a new quantity J called an impulse. An impulse is a vector quantity, just like a force, but it has the units of momentum. Applying an impulse produces an instantaneous change in the velocity of a body. To determine the effects of a given impulse J, we imagine a large force F that acts for a small time interval t. If we let F go to innity and t go to zero in such a way that F t= J (84)
then we can derive the effect of J on a bodys velocity by considering how the velocity would change if we let the force F act on it for t time.
G42
Witkin/Baraff/Kass
A n
a ( t 0 ) pb ( t0 ) p
B contact force
Figure 22: The vector pa (t0 ) pb (t0 ) is perpendicular to n(t0 ); the bodies are in resting contact. A contact force may be necessary to prevent bodies from accelerating towards each other.
A n
a ( t 0 ) pb ( t0 ) p
Figure 23: Colliding contact. The relative velocity pa (t0 ) pb (t0 ) is directed inwards, opposite n(t0 ). Unless the relative velocity is abruptly changed, inter-penetration will occur immediately after time t0 .
G43
Witkin/Baraff/Kass
For example, if we apply an impulse J to a rigid body with mass M, then the change in linear velocity v is simply v= J . M (85)
Equivalently, the change in linear momentum P is simply P = J. If the impulse acts at the point p, then just as a force produces a torque, J produces an impulsive torque of impulse = ( p x(t)) J. (86)
As one would imagine, the impulsive torque impulse also gives rise to a change in angular momentum L of L = impulse. The change in angular velocity is simply I 1 (t0 )impulse, assuming the impulse was applied at time t0 . When two bodies collide, we will apply an impulse between them to change their velocity. For frictionless bodies, the direction of the impulse will be in the normal direction, n(t0 ). Thus, we can write the impulse J as J = jn(t0 ) (87)
where j is an (as yet) undetermined scalar that gives the magnitude of the impulse. Well adopt the convention that the impulse J acts positively on body A, that is, A is subject to an impulse of + jn(t0 ), while body B is subject to an equal but opposite impulse jn(t0 ) (gure 24). We compute j by using an empirical law for collisions. Lets let p (t0 ) denote the velocity of the contact vertex of A prior a to the impulse being applied, and let p+ (t0 ) denote the velocity after we apply the impulse J. Let a p (t0 ) and p+ (t0 ) be dened similarly. Using this notation, the initial relative velocity in the normal b b direction is v = n(t0 ) ( p (t0 ) p (t0 )); a b rel after the application of the impulse, a b v+ = n(t0 ) ( p+ (t0 ) p+ (t0 )). rel The empirical law for frictionless collisions says simply that v+ = v . rel rel (810) (89) (88)
The quantity is called the coefcient of restitution and must satisfy 0 1. If = 1, then v+ = v , and the collision is perfectly bouncy; in particular, no kinetic energy is lost. At the rel rel other end of the spectrum, = 0 results in v+ = 0, and a maximum of kinetic energy is lost. After rel this sort of collision, the two bodies will be in resting contact at the contact point p (gure 25). Calculating the magnitude j of the impulse J = jn(t0 ) is fairly simple, although the equations are a bit tedious to work through. Lets dene the displacements ra and rb as p xa (t0 ), and p xb (t0 ). If we let v (t0 ) and (t0 ) be the pre-impulse velocities of body A, and v+ (t0 ) and + (t0 ) be the a a a a post-impulse velocities, we can write p+ (t0 ) = v+ (t0 ) + + (t0 ) ra a a a (811)
G44
Witkin/Baraff/Kass
j ( t0 ) n
A n( t0 ) p
n j ( t0 )
Figure 24: The impulse between two bodies at a contact point. An impulse of jn(t0 ) acts on A, while an impulse of jn(t0 ) acts on B.
(a) n p a ( t0 ) pb ( t0 ) p
(b) n p
a+( t0 ) pb+ ( t0 ) p
(c) n p
a+( t0 ) pb+ ( t0 ) p
(d) n p
a+( t 0 ) pb+ ( t0 ) p
Figure 25: (a) The relative velocity before application of the impulse. (b) The component of the relative velocity in the n(t0 ) direction is reversed for an = 1 collision. The relative velocity perpendicular to n(t0 ) remains the same. (c) A collision with 0 < < 1. The bodies bounce away in the n(t0 ) direction with less speed than they approached. (d) A collision with = 0. The bodies do not bounce away from each other, but the relative velocity perpendicular to n(t0 ) is unaffected by the collision.
G45
Witkin/Baraff/Kass
along with v+ (t0 ) = v (t0 ) + a a jn(t0 ) Ma and

1 + (t0 ) = (t0 ) + Ia (t0 ) ra jn(t0 ) a a
(812)
where Ma is the mass of body A, and Ia (t0 ) is its inertia tensor. Combining the two previous equations yields p+ (t0 ) = v (t0 ) + a a jn(t0 ) Ma
1 + (t0 ) + Ia (t0 ) ra jn(t0 ) a
ra ra (813)
= v (t0 ) + (t0 ) ra + a a = p + j a
jn(t0 ) Ma
1 + Ia (t0 ) ra jn(t0 )
n(t0 ) 1 + Ia (t0 ) ra n(t0 ) Ma
ra .
It is important to note the form of p+ (t0 ): it is a simple linear function of j. For body B, an opposite a impulse jn(t0 ) acts, yielding b p+ (t0 ) = p j b This yields b a b p+ (t0 ) p+ = p+ (t0 ) p+ + j a
1 Ia (t0 )
n(t0 ) 1 + Ib (t0 ) rb n(t0 ) Mb n(t0 ) n(t0 ) + + Ma Mb ra +

1 Ib (t0 )
rb .
(814)
(815) rb .
ra n(t0 )
rb n(t0 )
To calculate v+ , we dot this expression with n(t0 ). Since n(t0 ) is of unit length, n(t0 ) n(t0 ) = 1, rel and we obtain v+ = n(t0 ) ( p+ (t0 ) p+ ) a b rel a b = n(t0 ) ( p (t0 ) p ) + j
1 n(t0 ) Ia (t0 ) ra n(t0 )
1 1 + + Ma Mb
1 ra + n(t0 ) Ib (t0 ) rb n(t0 )
rb
(816)
= v + j rel
1 1 + + Ma Mb
1 ra + n(t0 ) Ib (t0 ) rb n(t0 )
1 n(t0 ) Ia (t0 ) ra n(t0 )
rb .
By expressing v+ in terms of j and v , we can compute j according to equation (810). If we rel rel substitute equation (816) into equation (810), we get v + j rel 1 1 1 + + n(t0 ) Ia (t0 ) ra n(t0 ) Ma Mb n(t0 )
1 Ib (t0 )
ra + (817) rb = v . rel
rb n(t0 )
G46
Witkin/Baraff/Kass
Finally, solving for j, j=

1 Ma
1 Mb
+ n(t0 )
(1 + )v rel 1 1 Ia (t0 ) ra n(t0 ) ra + n(t0 ) Ib (t0 ) rb n(t0 )
rb
. (818)
Lets consider some actual code (written for clarity, not speed). First, we determine if two bodies are in colliding contact.
* Operators: if `x' and `y' are triples, assume that `x y' is their cross product, and `x * y' is their dot product. * * Return the velocity of a point on a rigid body * triple pt_velocityBody *body, triple p return body- v + body- omega p - body- x;
* Return TRUE if bodies are in colliding contact. The parameter `THRESHOLD' is a small numerical tolerance used for deciding if bodies are colliding. * boolean collidingContact *c triple double ifvrel padot = pt_velocityc- a, p, pbdot = pt_velocityc- b, p; vrel = c- n * padot - pbdot; * moving away * * resting contact * * vrel -THRESHOLD * * p (t0 ) * a * p (t0 ) * b * v * rel
THRESHOLD return FALSE; ifvrel -THRESHOLD return FALSE; else return TRUE;
Next, well loop through all the contact points until all the collisions are resolved, and actually compute and apply an impulse.
void collisionContact *c, double epsilon triple padot = pt_velocityc- a, p, pbdot = pt_velocityc- b, p, n = c- n, * p (t0 ) * a * p (t0 ) * b * n(t0 ) *
G47
Witkin/Baraff/Kass
double
ra = p - c- a- x, * ra * rb = p - c- b- x; * rb * vrel = n * padot - pbdot, * v * rel numerator = -1 + epsilon * vrel;
* We'll calculate the denominator in four parts * double term1 = 1 c- a- mass, term2 = 1 c- b- mass, term3 = n * c- a- Iinv * ra n ra, term4 = n * c- b- Iinv * rb n rb; * Compute the impulse magnitude * double j = numerator triple force = j * n; * ccccterm1 + term2 + term3 + term4;
Apply the impulse to the bodies * a- P += force; b- P -= force; a- L += ra force; b- L -= rb force;
* recompute auxiliary variables * c- a- v = c- a- P c- a- mass; c- b- v = c- b- P c- b- mass; c- a- omega = c- a- Iinv * c- a- L; c- b- omega = c- b- Iinv * c- b- L;
void find_all_collisionsContact contacts boolean had_collision; double epsilon = .5; do had_collision = FALSE;
, int ncontacts
forint i = 0; i ncontacts; i++ ifcolliding&contacts i collision&contacts i , epsilon; had_collision = TRUE; * Tell the solver we had a collision * ode_discontinuous;
G48
Witkin/Baraff/Kass
whilehad_collision == TRUE;
Note several things. First, = .5 was chosen arbitrarily. In a real implementation, wed allow the user to use different values of depending on which two bodies were colliding. Also, every time we nd a collision, we have to rescan the list of contacts, since bodies that were at rest may no longer be so, and new collisions may develop. If there are initially several collisions to be resolved (such as a cube dropped at onto a plane, with all four vertices colliding at once), the order of the contact list may have an effect on the simulation. There is a way to compute impulses at more than one contact point at a time, but it more complicated, and is based on the concepts used for resting contact in the next section. For further information, see Baraff[1]. Incidentally, if you want to have certain bodies that are xed, and cannot be moved (such as 1 oors, or walls), you can use the following trick: for such bodies, let mass be zero; also let the inverse inertia tensor also be the 3 3 zero matrix. You can either special-case the code to check if a body is supposed to be xed, or you can recode the denition of RigidBody to have the variable invmass instead of mass. For ordinary bodies, invmass is the inverse of the mass, while for xed bodies, invmass is zero. The same goes for the inertia tensor. (Note that nowhere in any of the dynamics computations (including the next section) is the mass or inertia tensor ever used; only their inverses are used, so you wont have to worry about dividing by zero.) The same trick can be used in the next section on resting contact to simulate bodies that can support any amount of weight without moving.
9 Resting Contact
The case of resting contact, when bodies are neither colliding nor separating at a contact point, is the last (and hardest) dynamics problem well tackle in these notes. To implement whats in this section, youll have to obtain a fairly sophisticated piece of numerical software, which well describe below. At this point, lets assume we have a conguration with n contact points. At each contact point, bodies are in resting contact, that is, the relative velocity vrel , from section 8, is zero (to within the numerical tolerance THRESHOLD). We can say that this is so, because colliding contact is eliminated by the routine find_all_collisions, and any contact points with vrel larger than THRESHOLD can be safely ignored, since the bodies are separating there. As was the case for colliding contact, at each contact point, we have a contact force that acts normal to the contact surface. For the case of colliding contact, we had an impulse jn(t0 ) where j was an unknown scalar. For resting contact, at each contact point there is some force f i ni (t0 ), where f i is an unknown scalar, and ni (t0 ) is the normal at the ith contact point (gure 26). Our goal is to determine what each f i is. In computing the f i s, they must all be determined at the same time, since the force at the ith contact point may inuence one or both of the bodies of the j contact point. In section 8, we wrote how the velocity of the contact points pa (t0 ) and pb (t0 ) changed with respect to j. Well do the same thing here, but now well have to describe how the acceleration of pa (t0 ) and pb (t0 ) depends on each f i . For colliding contact, we had an empirical law which related the impulse strength j to the relative velocity and a coefcient of restitution. For resting contact, we compute the f i s subject to not one, but three conditions. First, the contact forces must prevent inter-penetration; that is, the contact
G49
Witkin/Baraff/Kass
A n3 ( t0 ) B n1 (t 0 ) p3 n2 ( t 0 ) n4 ( t 0 ) n5 ( t0 ) C
p4
D p1 p2 p5
Figure 26: Resting contact. This conguration has ve contact points; a contact force acts between pairs of bodies at each contact point. forces must be strong enough to prevent two bodies in contact from being pushed towards one another. Second, we want our contact forces to be repulsive; contact forces can push bodies apart, but can never act like glue and hold bodies together. Last, we require that the force at a contact point become zero if the bodies begin to separate. For example, if a block is resting on a table, some force may act at each of the contact points to prevent the block from accelerating downwards in response to the pull of gravity. However, if a very strong wind were to blow the brick upwards, the contact forces on the brick would have to become zero at the instant that the wind accelerated the brick off the table. Lets deal with the rst condition: preventing inter-penetration. For each contact point i, we construct an expression di (t) which describes the separation distance between the two bodies near the contact point at time t. Positive distance indicates the bodies have broken contact, and have separated at the ith contact point, while negative distance indicates inter-penetration. Since the bodies are in contact at the present time t0 , we will have di (t0 ) = 0 (within numerical tolerances) for each contact point. Our goal is to make sure that the contact forces maintain di (t) 0 for each contact point at future times t > t0 . For vertex/face contacts, we can immediately construct a very simple function for di (t). If pa (t) and pb (t) are the contact points of the ith contact, between bodies A and B, than the distance between the vertex and the face at future times t t0 is given by di (t) = ni (t) ( pa (t) pb (t)). (91)
At time t, the function d(t) measures the separation between A and B near pa (t). If di (t) is zero, then the bodies are in contact at the ith contact point. If di (t) > 0, then the bodies have lost contact at the ith contact point. However, if di (t) < 0, then the bodies have inter-penetrated, which is what we need to avoid (gure 27). The same function can also be used for edge/edge contacts; since ni (t)
G50
Witkin/Baraff/Kass
(a)
(b)
(c)
A pa ( t) n( t) pb ( t) B B n( t) A A pa ( t) pb ( t) pa ( t) B n( t) pb ( t)
Figure 27: (a) The displacement pa (t) pb (t), indicated by an arrow, points in the same direction as n(t). Thus, the distance function d(t) would be positive. (b) The distance function d(t) is zero. (c) The displacement pa (t) pb (t) points in the opposite direction as n(t). The distance function d(t) is negative, indicating inter-penetration. points outwards from B towards A (by convention), ni (t) ( pa (t) pb (t)) will be positive if the two contacting edges move so as to separate the bodies. Since di (t0 ) = 0, we have to keep di (t0 ) from decreasing at time t0 ; that is, we have to have i (t0 ) 0. What is di (t0 )? Differentiating, d di (t) = ni (t) ( pa (t) pb (t)) + ni (t) ( pa (t) pb (t)). (92)
Since di (t) describes the separation distance, di (t) will describe the separation velocity at time t. However, at time t0 , pa (t0 ) = pb (t0 ), which means that di (t0 ) = ni (t0 ) ( pa (t0 ) pb (t0 )). This should look familiar: its vrel from the previous section! The function di (t0 ) is a measure of how the bodies are separating, and for resting contact, we know that d(t0 ) is zero, because the bodies are neither moving towards nor away from each other at a contact point. At this point then, we have di (t0 ) = di (t0 ) = 0. Now well look at di (t0 ). If we differentiate equation (92), we get di (t) = ni (t) ( pa (t) pb (t)) + ni (t) ( pa (t) pb (t)) + ni (t) ( pa (t) pb (t)) + ni (t) ( pa (t) pb (t)) = ni (t) ( pa (t) pb (t)) + 2ni (t) ( pa (t) pb (t)) + ni (t) ( pa (t) pb (t)). Since pa (t0 ) = pb (t0 ), we can write di (t0 ) as d(t0 ) = ni (t0 ) ( pa (t0 ) pb (t0 )) + 2ni (t0 ) ( pa (t0 ) pb (t0 )). (94) (93)
G51
Witkin/Baraff/Kass
The quantity di (t0 ) measures how the two bodies are accelerating towards each other at the contact i (t0 ) > 0, the the bodies have an acceleration away from each other, and contact will point p. If d break immediately after t0 . If di (t0 ) = 0, then contact remains. The case di (t0 ) < 0 must be avoided, for this indicates the bodies are accelerating towards each other. Note that if ni (t0 ) is a constant (if i (t0 ) is zero, leading to further simplications. body B is xed), then n Thus, we satisfy our rst condition for contact forces by writing the constraint di (t0 ) 0 (95)
for each contact point. Since the acceleration di (t0 ) depends on the contact forces, this is really a constraint on the contact forces. Lets turn our attention to the second and third constraints. Since contact forces must always be repulsive, each contact force must act outward. This means that each f i must be positive, since a force of f i ni (t0 ) acts on body A, and ni (t0 ) is the outwards pointing normal of B. Thus, we need fi 0 (96)
for each contact point. The third constraint is expressed simply in terms of f i and di (t0 ). Since the contact force f i ni (t0 ) must become zero if contact is breaking at the ith contact, this says that f i must be zero if contact is breaking. We can express this constraint by writing f i di (t0 ) = 0; (97)
if contact is breaking, di (t0 ) > 0 and equation (97) is satised by requiring f i = 0. If contact is not breaking, then di (t0 ) = 0, and equation (97) is satised regardless of f i . In order to actually nd f i s which satisfy equations (95), (96), and (97), we need to express each di (t0 ) as a function of the unknown f i s. It will turn out that we will be able to write each di (t0 ) in the form di (t0 ) = ai1 f1 + ai2 f2 + + ain f n + bi . In matrix parlance, this means we will be able to write f1 d1 (t0 ) . . . = A . + . . fn dn (t0 ) bi . . . bn (98)
(99)
where A is the n n matrix of the aij coefcients of equation (98). Although the code needed to calculate the aij s and the bi s is not too complicated, working out the derivations on which the code is based is somewhat tedious. The derivations are worked out in appendix D, along with code to compute the matrix of aij s and bi s. Appendix D gives an implementation of the routines
void compute_a_matrixContact contacts bigmatrix &a; , int ncontacts,
G52
Witkin/Baraff/Kass
void
compute_b_vectorContact contacts vector &b;
, int ncontacts,
where the types bigmatrix and vector represent matrices and vectors of arbitrary size. The rst routine computes the aij s, while the second routine computes the bi s. Once weve computed all this, we can think about solving equations (95), (96), and (97). This system of equations forms what is called a quadratic program (QP); that is, f i s that satisfy these three equations are found by an algorithm called quadratic programming. Not all quadratic programs can be solved efciently, but because our contact forces are all normal to the contact surfaces (that is, they do not involve friction), it turns out that our QP can always be solved efciently. One interesting thing to note is that QP codes can easily handle the case di (t0 ) = 0 instead of di (t0 ) 0. We use i (t0 ) = 0 (and also drop the constraint f i 0) if we wish to constrain two bodies to never separate d at a contact point. This enables us to implement hinges, and pin-joints, as well as non-penetration constraints during simulation. Quadratic programming codes arent terribly common though; certainly, they are not nearly as common as linear equation codes, and are much harder to implement. The quadratic programming routines used by the author were obtained from the Department of Operations Research at Stanford University. See Gill et al.[7, 8, 9] for further details. More recently, we have been using code described by Baraff[3] to solve the quadratic programs. If you are determined to really implement this, we suggest a thorough study of this paper (excepting for the section on contact with friction). At any rate, lets assume that youve got a working QP solver at your disposal. Well assume that you pass the matrix A, and the vector of bi s to the QP solver, and you get back the vector of f i s. Lets pretend the interface is
void qp_solvebigmatrix &a, vector &b, vector &f;
Lets see how to compute all the resting contact forces. The following routine is presumably called from Compute_Force_and_Torque, after find_collisions has been called.
void compute_contact_forcesContact contacts , int ncontacts, double t
* We assume that every element of contacts represents a contact in resting contact. Also, we'll assume that for each element of Bodies , the `force' and `torque' fields have been set to the net external force and torque acting on the body, due to gravity, wind, etc., perhaps by a call to Compute_External_Force_and_Torque_for_all_Bodiest; * * Allocate n n matrix `amat' and n-vectors `fvec', and `bvec'. *
G53
Witkin/Baraff/Kass
bigmatrix amat = new bigmatrixncontacts, ncontacts; vector bvec = new vectorncontacts, fvec = new vectorncontacts; * Compute aij and bi coefficients * compute_a_matrixcontacts, ncontacts, amat; compute_b_vectorcontacts, ncontacts, bvec; * Solve for f j 's * qp_solveamat, bmat, fvec; * Now add the resting contact forces we just computed into the `force' and `torque' field of each rigid body. * forint i = 0; i double triple RigidBody ncontacts; i++ f = fvec i ; n = contacts i - n; *A = contacts i - a, *B = contacts i - b; * * * *
fi * ni (t0 ) * body A * body B *
* apply the force `f n' positively to A... * A- force += f * n; A- torque += contacts i .p - A- x * f*n; * and negatively to B * B- force -= f * n; B- torque -= contacts i .p - B- x * f*n;
Thats pretty much it! Now that the resting forces have been computed and combined with the external forces, we return control to the ODE solver, and each body goes merrily along its way, in a physically correct manner, without inter-penetration.
G54
Witkin/Baraff/Kass
Appendix A Motion Equation Derivations

In this appendix, well ll in some of the missing details from section 2, with regards to the equa tions P(t) = F(t), L(t) = (t), and L(t) = I(t)(t). The derivation method used here is somewhat nonstandard, and was proposed by Andy Witkin. The derivation in this appendix is (we feel) much shorter and considerably more elegant than the one found in traditional sources such as Goldstein[10]. Weve described the external force acting on a rigid body in terms of forces Fi (t), where Fi (t) is the external force acting on the ith particle. However, for a rigid body to maintain its shape, there must be some internal constraint forces that act between particles in the same body. We will make the assumption that these constraint forces act passively on the system and do not perform any net work. Let Fci (t) denote the net internal constraint force acting on the ith particle. The work performed by Fci on the ith particle from time t0 to t1 is
t1 t0
Fci (t) r i (t) dt
where r i (t) is the velocity of the ith particle. The net work over all the particles is the sum
t1 i t0
Fci (t) r i (t) dt =
t1 t0 i
Fci (t) r i (t) dt,
which must be zero for any interval t0 to t1 . This means that the integrand Fci (t) r i (t) (A1)
is itself always zero for any time t. (Henceforth well just write these expressions as Fci r i = 0.) We can use this fact to eliminate any mention of the constraint forces Fci from our derivations. First, some quick notes about the operator dened in section 2.3: since a b = a b, and a b = b a, we get a b = b a = b a. Since a is an antisymmetric matrix, (a ) T = a . Finally, since the operator is a linear operator, d (a) = (a ) = (a ) dt and for a set of vectors ai a = i ai
(A2)
(A3)
(A4)
(A5)
G55
Witkin/Baraff/Kass
Recall that we can write the velocity r i as r i = v + (ri x) where ri is the particles location, x is the position of the center of mass, and v and are linear and angular velocity. Letting ri = ri x and using the notation, r i = v + ri = v ri . Substituting this into Fci r i , which is always zero, yields Fci (v ri ) = 0.

(A6)
(A7)
Note that this equation must hold for arbitrary values of v and . Since v and are completely independent, if we choose to be zero, then Fci v = 0 for any choice of v, from which we conclude that in fact Fci = 0 is always true. This means that the constraint forces produce no net force. Similarly, choosing v to be zero we see that Fci (ri ) = 0 for any . Rewriting Fci (ri ) as Fci T (ri ) we get that Fci T ri = for any , so
Fci T ri
=0
(A8)
Fci T ri = 0T . Transposing, we have (ri ) T Fci =
(ri ) Fci =
ri Fci = 0
(A9)
which means that the internal forces produce no net torque. We can use the above to derive the rigid body equations of motion. The net force on each particle is the sum of the internal constraint force Fci and the external force Fi . The acceleration ri of the ith particle is ri = d d r i = (v ri ) = v r i ri . dt dt (A10)
Since each individual particle must obey Newtons law f = ma, or equivalently ma f = 0, we have mi ri Fi Fci = mi (v r i ri ) Fi Fci = 0 for each particle. To derive P = F = Fi , we sum equation (A11) over all the particles. We obtain mi (v r i ri ) Fi Fci = 0. Breaking the large sum into smaller ones, mi (v r i ri ) Fi Fci = mi v mi v mi v d dt mir i m i ri mir i

(A11)
(A12)
m i ri m i ri m i ri

Fi Fi Fi
Fci = Fci = Fci = 0. (A13)
G56
Witkin/Baraff/Kass
Since we are in a center-of-mass coordinate system, equation (220) from section 2.6 tells us that d miri = 0, which also means that dt mi ri = 0. Removing terms with miri , and the term Fci from the above equation yields mi v Fi = 0 (A14)
or simply M v = P = Fi as advertised. = = ri Fi , we again start with equation (A11). Multiplying both sides by To obtain L ri yields ri mi (v r i ri ) ri Fi ri Fci = ri 0 = 0. Summing over all the particles, we obtain ri m i v Since

(A15)
ri m i r i
ri m i ri
ri Fi
ri Fci = 0.
(A16)
ri Fci = 0, we can rearrange this to obtain m i ri
m i ri r i
m i ri ri
ri Fi = 0.
(A17)
Using
mi ri = 0, we are left with m i ri r i ri Fi =

m i ri ri
ri Fi = 0
(A18)
or, recognizing that
ri Fi = , m i ri r i

m i ri ri
= .
(A19)
Were almost done now: if we refer back to the matrix dened by the notation, one can easily verify the relation that the matrix a a is equivalent to the matrix (a T a)1 aa T where 1 is the 3 3 identity matrix. (This relation is equivalent to the vector rule a (b c) = ba T c ca T b.) Thus mi ri ri = Substituting into equation (A19), this yields miri r i

mi ((ri ri )1 ri ri ) = I(t).
(A20)
+ I(t) = .
(A21)
The above expression is almost acceptable, as it gives an expression for in terms of , ex cept that it requires us to evaluate the matrix miri r , which is as expensive as computing the i inertia tensor from scratch. Well use one last trick here to clean things up. Since ri = ri and ri = ri , we can write m i r i ri =

mi ( ri ) ( ri ) =
mi ( ri ) ( ri ) = 0.
(A22)
G57
Witkin/Baraff/Kass
Thus, we can add
mir i ri = 0 to equation (A21) to obtain miri r i mi r i ri

+ I(t) = .
(A23)
Finally, since d I (t) = dt we have d I (t) + I(t) = (I(t)) = . dt Since L(t) = I(t)(t), this leaves us with the nal result that L (t) = . (A26) (A25) miri ri =

miri r i mir i ri
(A24)
Appendix B Quaternion Derivations

A formula for q(t) is derived as follows. Recall that the angular velocity (t) indicates that the body is instantaneously rotating about the (t) axis with magnitude |(t)|. Suppose that a body were to rotate with a constant angular velocity (t). Then the rotation of the body after a period of time t is represented by the quaternion [cos |(t)| t |(t)| t (t) , sin ]. 2 2 |(t)|
Let us compute q(t) at some particular instant of time t0 . At times t0 + t (for small t), the orientation of the body is (to within rst order) the result of rst rotating by q(t0 ) and then further rotating with velocity (t0 ) for t time. Combining the two rotations, we get q(t0 + t) = [cos |(t0 )| t |(t0 )| t (t0 ) , sin ] q(t0 ). 2 2 |(t0 )| (B1)
Making the substitution t = t0 + t, we can express this as q(t) = [cos |(t0 )|(t t0 ) |(t0 )|(t t0 ) (t0 ) , sin ]q(t0 ). 2 2 |(t0 )| (B2)
Let us differentiate q(t) at time t0 . First, since q(t0 ) is a constant, let us differentiate [cos At time t = t0 , d |(t0 )|(t t0 ) |(t0 )| |(t0 )|(t t0 ) cos = sin dt 2 2 2 |(t0 )| = sin 0 = 0. 2 |(t0 )|(t t0 ) |(t0 )|(t t0 ) (t0 ) , sin ]. 2 2 |(t0 )|
(B3)
G58
Witkin/Baraff/Kass
Similarly, d |(t0 )|(t t0 ) |(t0 )| |(t0 )|(t t0 ) sin = cos dt 2 2 2 |(t0 )| |(t0 )| = cos 0 = . 2 2 Thus, at time t = t0 , q(t) = |(t0 )|(t t0 ) |(t0 )|(t t0 ) (t0 ) , sin ] q(t0 ) 2 2 |(t0 )| |(t0 )|(t t0 ) |(t0 )|(t t0 ) (t0 ) [cos , sin ] q(t0 ) 2 2 |(t0 )| |(t0 )| (t0 ) = [0, ] q(t0 ) 2 |(t0 )| d dt d = dt [cos = [0, 1 (t0 )] q(t0 ) = 1 [0, (t0 )] q(t0 ). 2 2 The product [0, (t0 )] q(t0 ) is abbreviated to the form (t0 )q(t0 ); thus, the general expression for q(t) is q(t) = 1 (t)q(t). 2 (B6)
(B4)
(B5)
Appendix C Some Miscellaneous Formulas

C.1 Kinetic Energy
The kinetic energy T of a rigid body is dened as T=
1 T 2 mi r i r i .
(C1)
Letting ri = ri x, we have r i = v(t) + ri . Thus T= = =

1 2 1 T 2 mir i r i T 1 2 m i (v + ri ) (v + ri )
mi v v +
T
v m i ri +
T
1 2
mi (ri ) (ri )
T
(C2)
T
= 1 vT 2 Using
T
mi v + v T
m i ri
+ 1 T 2
mi (ri ) ri
mi ri = 0 and (ri ) = ri , we have T = 1 v T Mv + 1 T 2 2 miri ri

= 1 (v T Mv + T I) 2
(C3)
since I = miri ri from appendix A. Thus, the kinetic energy can be decomposed into two terms: a linear term 1 v T Mv, and an angular term 1 T I. 2 2
G59
Witkin/Baraff/Kass
C.2 Angular Acceleration

It is often necessary to compute (t). Since L(t) = I(t)(t), we know (t) = I 1 (t)L(t). Thus, (t) = I 1 (t)L(t) + I 1 (t) L(t). Since we know that L(t) = (t), let us consider I 1 (t). From equation (240),
1 I 1 (t) = R(t)Ibody R(t) T ,
(C4)
so
1 1 I 1 (t) = R(t)Ibody R(t) T + R(t)Ibody R(t) T .
(C5)
Since R(t) = (t) R(t), R(t)T = ((t) R(t)) T = R(t) T ((t) ) T . Since (t) is antisymmetric, (i.e. ((t) ) T = (t)), R(t) T = R(t) T (t) . This yields
1 1 I 1 (t) = R(t)Ibody R(t) T + R(t)Ibody (R(t) T (t) ) 1 = (t) R(t)Ibody R(t) T I 1 (t)(t)
(C6)
(C7)
(C8)
= (t) I 1 (t) I 1 (t)(t) . Then (t) = I 1 (t)L(t) + I 1 (t) L(t) = (t) I 1 (t) I 1 (t)(t) L(t) + I 1 (t) L(t) = (t) I 1 (t)L(t) I 1 (t)(t) L(t) + I 1 (t) L(t). But since I 1 (t)L(t) = (t), the rst term, (t) I 1 (t)L(t) is equivalent to (t) (t), or (t) (t), which is zero. This leaves the nal result of (t) = I 1 (t)(t) L(t) + I 1 (t) L(t) = I 1 (t)(t) L(t) + I 1 (t) L(t) =I
1
(C9)
(t)(L(t) (t)) + I
(t) L(t)
(C10)
= I 1 (t)(L(t) (t) + L (t)).
G60
Witkin/Baraff/Kass
We can see from this that even if no forces act, so that L (t) is zero, (t) can still be non-zero. (In fact, this will happen whenever the angular momentum and angular velocities point in different directions, which in turn occurs when the body has a rotational velocity axis that is not an axis of symmetry for the body.)
C.3 Acceleration of a Point

Given a point of a rigid body with world space coordinate p(t), it is often necessary to compute p(t). Let the body space coordinate that transforms at time t to p(t) be p0 ; then p(t) = R(t) p0 + x(t) If we let r(t) = p(t) x(t), then p(t) = R(t) p0 + x(t) = (t) R(t) p0 + v(t) = (t) (R(t) p0 + x(t) x(t)) + v(t) = (t) ( p(t) x(t)) + v(t) = (t) r(t) + v(t). Then p(t) = (t) r(t) + (t) r (t) + v(t) = (t) r(t) + (t) ((t) r(t)) + v(t). We can interpret this as follows. The rst term, (t) r(t) is the tangential acceleration of the point; that is, (t) r(t) is the acceleration perpendicular to the displacement r(t) as a result of the body being angularly accelerated. The second term, (t) ((t) r(t)) is the centripetal acceleration of the point; this centripetal acceleration arises because the body is rigid, and points on the body must rotate in a circular orbit about the center of mass. The last term, v(t) is the linear acceleration of the point due to the linear acceleration of the center of mass of the body. (C12) (C11)
Appendix D Resting Contact Derivations

If youre determined to implement resting contact in your simulator, youll need the derivations and the code in this appendix. This is probably not a fun appendix to work through; then again, this wasnt a fun appendix to write! The derivations in here are somewhat terse, but the code at the end of the appendix will hopefully make things clearer.
D.1 Derivations
We need to express di (t0 ) in terms of all the unknown f i s. It will turn out that well be able to write i (t0 ) in the form each d di (t0 ) = ai1 f1 + ai2 f2 + + ain f n + bi . (D1)
G61
Witkin/Baraff/Kass
Given i and j, we need to know how di (t0 ) depends on f j , that is, we need to know aij . Also, we need to compute the constant term bi . Lets start by determining aij and ignoring the constant part bi . Well assume the ith contact involves two bodies A and B. From equation (94), we can write di (t0 ) as d(t0 ) = ni (t0 ) ( pa (t0 ) pb (t0 )) + 2ni (t0 ) ( pa (t0 ) pb (t0 )) (D2)
where pa (t0 ) = pi = pb (t0 ) is the contact point for the ith contact at time t0 . The term 2ni (t0 ) ( pa (t0 ) pb (t0 )) is a velocity dependent term (i.e. you can immediately calculate it without knowing the forces involved), and is part of bi , so well ignore this for now. So we only need to know how pa (t0 ) and pb (t0 ) depend on f j , the magnitude of the jth contact force. Consider the jth contact. If body A is not one of the bodies involved in the jth contact, then pa (t0 ) is independent of f j , because the jth contact force does not act on body A. Similarly, if B is also not one of the two bodies involved in the jth contact, then pb (t0 ) is also independent of f j . (For example, in gure 26, the acceleration of the contact points at the rst contact is completely unaf fected by the contact force acting at the fth contact. Thus, d1 (t0 ) would be completely independent 5 (t0 ) is completely independent of f1 .) of f5 . Conversely, d Suppose though that in the jth contact, body A is involved. For deniteness, suppose that in the jth contact, a force of jn j (t0 ) acts on body A, as opposed to jn j (t0 ). Lets derive how pa (t0 ) is affected by the force jn j (t0 ) acting on A. From equation (C12), we can write pa (t) = va (t) + a (t) ra (t) + a (t) (a (t) ra (t)) (D3)
where ra (t) = pa (t) xa (t), and xa (t), va (t), and a (t) are all the variables associated with body A. We know that va (t) is the linear acceleration of body A, and is equal to the total force acting on A divided by the mass. Thus, a force of jn j (t0 ) contributes n j (t0 ) f j n j (t0 ) = fj ma ma to va (t) and thus pa (t). Similarly, consider a (t), from equation (C10):
1 1 a (t) = Ia (t)a (t) + Ia (t)(La (t) a (t))
(D4)
where a (t) is the total torque acting on body A. If the jth contact occurs at the point p j , then the force jn j (t0 ) exerts a torque of ( p j xa (t0 )) f j n j (t0 ). Thus, the angular contribution to pa (t0 ) is
1 f j Ia (t0 ) ( p j xa (t0 )) n j (t0 )
ra .
(D5)
The total dependence of pa (t0 ) on f j is therefore fj n j (t0 ) 1 + Ia (t0 ) ( p j xa (t0 )) n j (t0 ) ma ra .
G62
Witkin/Baraff/Kass
Now, if a force of f j n(t0 ) had acted on A instead, wed get the same dependence, but with a minus sign in front of f j . Clearly, pb (t0 ) depends on f j in the same sort of manner. Once we compute how pa (t0 ) and pb (t0 ) depend on f j , we combine the results together and take the dot product with ni (t0 ), to see how di (t0 ) depends on f j . This gives us aij . Confused? See the code below. We still need to compute bi . We know that di (t0 ) contains the constant term 2ni (t0 ) ( pa (t0 ) pb (t0 )). But we also have to take into account the contributions to pa (t0 ) and pb (t0 ) due to known ex ternal forces such as gravity, as well as the force-independent terms a (t0 ) (a (t0 ) ra ) and 1 (Ia (t0 )(a (t0 ) a (t0 ))) ra . If we let the net external force acting on A be Fa (t0 ), and the net external torque be a (t0 ), then from equations (D4) and (D5), we get that Fa (t0 ) contributes Fa (t0 ) ma and that a (t0 ) contributes
1 Ia (t0 )a (t0 ) ra .
Thus, the part of pa (t0 ) that is independent from all the f j s is Fa (t0 ) 1 1 + Ia (t0 )a (t0 ) ra + a (t0 ) (a (t0 ) ra ) + Ia (t0 )(a (t0 ) a (t0 )) ra ma
and similarly for pb (t0 ). To compute bi , we combine the constant parts of pa (t0 ), pb (t0 ), dot with ni (t0 ), and add the term 2ni (t0 ) ( pa (t0 ) pb (t0 )).
G63
Witkin/Baraff/Kass
D.2 Code
Heres the code to implement the above derivations. Lets start by computing the constant bi terms.
* return the derivative of the normal vector * triple compute_ndotContact *c ifc- vf * vertex face contact *
* The vector `n' is attached to B, so... * return c- b- omega c- n; else * This is a little trickier. The unit normal `n' is n = ea eb . ea eb Differentiating n with respect to time is left as an exercise... but here's some code * triple eadot = c- a- omega ea, ebdot = c- b- omega eb; n1 = ea * eb, z = eadot * eb + ea * ebdot; l = lengthn1; length; l;
* ea * * eb *
double n1 = n1
* normalize *
return z - z * n * n
void compute_b_vectorContact contacts forint i = 0; i
, int ncontacts, vector &b
ncontacts; i++
Contact *c = &contacts i ; Body *A = c- a, *B = c- b; triple n = c- n, ra = c- p - A- x, rb = c- p - B- x;
* ni (t0 ) * * p xa (t0 ) * * p xb (t0 ) *
* Get the external forces and torques * triple f_ext_a = A- force, f_ext_b = B- force, t_ext_a = A- torque,
G64
Witkin/Baraff/Kass
t_ext_b = B- torque; triple a_ext_part, a_vel_part, b_ext_part, b_vel_part;
* Operators: `' is for cross product, `*', is for dot products between two triples, or matrix-vector multiplication between a matrix and a triple. * * Compute the part of pa (t0 ) due to the external force and torque, and similarly for pb (t0 ). * a_ext_part = f_ext_a A- mass + A- Iinv * t_ext_a ra, b_ext_part = f_ext_b B- mass + B- Iinv * t_ext_b rb;
* Compute the part of pa (t0 ) due to velocity, and similarly for pb (t0 ). *
a_vel_part = A- omega A- omega ra + A- Iinv * A- L * A- omega ra; b_vel_part = B- omega B- omega rb + B- Iinv * B- L * B- omega rb;
* Combine the above results, and dot with ni (t0 ) * double k1 = n * a_ext_part + a_vel_part b_ext_part + b_vel_part; triple ndot = compute_ndotc;
* See section 8 for `pt_velocity' definition * double k2 = 2 * ndot * pt_velocityA, c- p pt_velocityB, c- p; b i = k1 + k2;
Computing the aij terms is a little more tricky, because we have to keep track of how the jth contact force affects the ith contact point. The following routine is not the most efcient way to do things, because with a good data structure, you can tell in advance which of the aij s are going to be zero. Still unless youre working with really huge numbers of contacts, not too much extra work will be done.
G65
Witkin/Baraff/Kass
void compute_a_matrixContact contacts
, int ncontacts, bigmatrix &a
forint i = 0; i ncontacts; i++ forint j = 0; j ncontacts; j++ a i,j = compute_aijcontacts i , contacts j ;
double compute_aijContact ci, Contact cj * If the bodies involved in the ith and jth contact are distinct, then aij is zero. * ifci.a != cj.a && ci.b != cj.b && ci.a != cj.b && ci.b != cj.a return 0.0; Body triple *A *B ni nj pi pj ra rb = = = = = = = = ci.a, ci.b; ci.n, cj.n, ci.p, cj.p, pi - A- x, pi - B- x;
* * * *
ni (t0 ) * n j (t0 ) * ith contact point location * jth contact point location *
* What force and torque does contact j exert on body A? * triple force_on_a = 0, torque_on_a = 0; ifcj.a == ci.a * force direction of jth contact force on A * force_on_a = nj; * torque direction * torque_on_a = pj - A- x nj; else ifcj.b == ci.a force_on_a = - nj; torque_on_a = pj - A- x nj;
* What force and torque does contact j exert on body B? * triple force_on_b = 0, torque_on_b = 0;
G66
Witkin/Baraff/Kass
ifcj.a == ci.b * force direction of jth contact force on B * force_on_b = nj; * torque direction * torque_on_b = pj - B- x nj; else ifcj.b == ci.b force_on_b = - nj; torque_on_b = pj - B- x nj;
* Now compute how the jth contact force affects the linear and angular acceleration of the contact point on body A * triple a_linear = force_on_a A- mass, a_angular = A- Iinv * torque_on_a * ra;
* Same for B * triple b_linear = force_on_b B- mass, b_angular = B- Iinv * torque_on_b * rb;
return ni * a_linear + a_angular - b_linear + b_angular;
G67
Witkin/Baraff/Kass
References
[1] D. Baraff. Analytical methods for dynamic simulation of non-penetrating rigid bodies. In Computer Graphics (Proc. SIGGRAPH), volume 23, pages 223232. ACM, July 1989. [2] D. Baraff. Curved surfaces and coherence for non-penetrating rigid body simulation. In Computer Graphics (Proc. SIGGRAPH), volume 24, pages 1928. ACM, August 1990. [3] D. Baraff. Fast contact force computation for nonpenetrating rigid bodies. Computer Graphics (Proc. SIGGRAPH), 28:2334, 1994. [4] J. Canny. Collision detection for moving polyhedra. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8(2), 1986. [5] P.A. Cundall. Formulation of a three-dimensional distinct element modelPart I. A scheme to represent contacts in a system composed of many polyhedral blocks. International Journal of Rock Mechanics, Mineral Science and Geomechanics, 25, 1988. [6] E.G. Gilbert and S.M. Hong. A new algorithm for detecting the collision of moving objects. In International Conference on Robotics and Automation, pages 813. IEEE, 1989. [7] P. Gill, S. Hammarling, W. Murray, M. Saunders, and M. Wright. Users guide for LSSOL: A Fortran package for constrained linear least-squares and convex quadratic programming. Technical Report Sol 86-1, Systems Optimization Laboratory, Department of Operations Research, Stanford University, 1986. [8] P. Gill, W. Murray, M. Saunders, and M. Wright. Users guide for QPSOL: A Fortran package for quadratic programming. Technical Report Sol 84-6, Systems Optimization Laboratory, Department of Operations Research, Stanford University, 1984. [9] P. Gill, W. Murray, M. Saunders, and M. Wright. Users guide for NPSOL: A Fortran package for nonlinear programming. Technical Report Sol 86-2, Systems Optimization Laboratory, Department of Operations Research, Stanford University, 1986. [10] H. Goldstein. Classical Mechanics. Addison-Wesley, Reading, 1983. [11] W. Meyer. Distance between boxes: Applications to collision detection and clipping. In International Conference on Robotics and Automation, pages 597602. IEEE, 1986. [12] P.M. Moore and J. Wilhelms. Collision detection and reponse for computer animation. In Computer Graphics (Proc. SIGGRAPH), volume 22, pages 289298. ACM, August 1988. [13] F.P. Preparata and M.I. Shamos. Computational Geometry. Springer-Verlag, New York, 1985. [14] W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. Numerical Recipes. Cambridge University Press, 1986. [15] R. Sedgewick. Algorithms. Addison-Wesley, 1983. [16] K. Shoemake. Animating rotation with quaternion curves. In Computer Graphics (Proc. SIGGRAPH), volume 19, pages 245254. ACM, July 1985. [17] B. Von Herzen, A. Barr, and H. Zatz. Geometric collisions for time-dependent parametric surfaces. In Computer Graphics (Proc. SIGGRAPH), volume 24, pages 3948. ACM, August 1990.
G68
Witkin/Baraff/Kass

An Introduction To Physically Based Modeling

Uploaded by

Copyright:

Available Formats

An Introduction To Physically Based Modeling

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

An Introduction To Physically Based Modeling

Uploaded by

Copyright:

Available Formats

An Introduction to Physically Based Modeling: Differential Equation Basics

Differential Equation Basics

Initial Value Problems

The derivative function

x = f (x,t) forms a vector field.

Follow the vectors

Initial Value Problem

An Introduction to Physically-Based Modeling

x(t + t) = x(t) + t f(x,t)

An Introduction to Physically-Based Modeling

Instability: off to Neptune!

An Introduction to Physically-Based Modeling

The Midpoint Method

An Introduction to Physically-Based Modeling

a. Compute an Euler step x = t f(x,t) b. Evaluate f at the midpoint f mid = f x + x , t + t 2 2

c. Take a step using the midpoint value

x(t + t) = x(t) + t f mid

The Midpoint Method

An Introduction to Physically-Based Modeling

Adaptive stepsizing is a highly recommended technique.

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

An Introduction to Physically Based Modeling: Particle System Dynamics

Particle System Dynamics

Basic Particle Systems

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

Position in Phase Space

Force Accumulator mass

particles n time particles n time

An Introduction to Physically-Based Modeling

Dim(State) Get/Set State

An Introduction to Physically-Based Modeling

A list of force objects to invoke

Particle Systems, with forces

Clear Force Accumulators

Invoke apply_force functions

Deriv Eval Loop

An Introduction to Physically-Based Modeling

p->f -= F->k * p->v

A Force Object: Viscous Drag

An Introduction to Physically-Based Modeling

Spatial Interaction Forces

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

a unit vector in the direction of l. Then, since l = x1 x2 , C C = , x1 l The time derivative of is C C = . x2 l

l l l (v1 v2 ) C= = . |l| |l|

Particle/Plane Collisions and Contact

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

An Introduction to Physically Based Modeling: Energy Functions and Stiffness

Energy Functions and Stiffness

Figure 1: Symmetric energy function E = (x2 + y2 )/2

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

An Introduction to Physically-Based Modeling

An Introduction to Continuum Dynamics for Computer Graphics

2 Spring with Mass