Transcript Physically Based Animation

Physically Based Animation
Ordinary Differential
Equations
Particle Dynamics
Rigid-Body Dynamics
Collisions
Lecture 22
6.837 Fall 2001
Particle
A single particle in 2-D moving in a flow field

x1 
Position x   
x 2 

v 1 
dx
Velocity v    , v 
dt
v 2 
x2
The flow field function dictates
particle velocity v  g x , t


x t 

g  x ,t 
x1
Lecture 22
Slide 2
6.837 Fall 2001
Vector Field
The flow field g(x,t) is a vector field that defines a vector for
any particle position x at any time t.
x2
g  x ,t 
x1
How would a particle move in this vector field?
Lecture 22
Slide 3
6.837 Fall 2001
Differential Equations
The equation v = g(x, t) is a first order differential equation:
dx
 g  x,t 
dt
The position of the particle is computed by integrating the
differential equation:
t
x t   x t 0    g  x , t  dt
t0
For most interesting cases, this integral cannot be computed
analytically.
Lecture 22
Slide 4
6.837 Fall 2001
Numeric Integration
Instead we compute the particle’s position by numeric
integration: starting at some initial point x(t0) we step along
the vector field to compute the position at each subsequent
time instant. This type of a problem is called an initial value
problem.
x2
x t 1 
x t 2 
x t 0 
x1
Lecture 22
Slide 5
6.837 Fall 2001
Euler’s Method
Euler’s method is the simplest solution to an initial value
problem. Euler’s method starts from the initial value and takes
small time steps along the flow:
x t  t   x t   t g  x , t 
Why does this work?
Let’s look at a Taylor series expansion of function x(t):
d x t 2 d 2 x
x t  t   x t   t


2
dt
2 dt
Disregarding higher-order terms and replacing the first
derivative with the flow field function yields the equation for
the Euler’s method.
Lecture 22
Slide 6
6.837 Fall 2001
Other Methods
Euler’s method is the simplest numerical method. The error is
proportional to t 2. For most cases, the Euler’s method is
inaccurate and unstable requiring very small steps.
x2
Other methods:




Midpoint (2nd order Runge-Kutta)
Higher order Runge-Kutta (4th order, 6th order)
Adams
Adaptive Stepsize
Lecture 22
Slide 7
x1
6.837 Fall 2001
Particle in a Force Field
What is a motion of a particle in a force field?
The particle moves according to Newton’s Law:
d 2x f

 f  ma 
2
m
dt
The mass m of a particle describes the particle’s inertial
properties: heavier particles are easier to move than lighter
particles. In general, the force field f(x, v, t) may depend on
the time t and particle’s position x and velocity v.
Lecture 22
Slide 8
6.837 Fall 2001
Second-Order Differential Equations
Newton’s Law yields an ordinary differential equation of second
order:
d 2 x t  f  x , v , t 

2
dt
m
A clever trick allows us to reuse the same numeric
differentiation solvers for first-order differential equations. If
we define a new phase space vector y, which consists of
particle’s position x and velocity v, then we can construct a
new first-order differential equation whose solution will also
solve the second-order differential equation.
x 
y   ,
v 
Lecture 22
d y d x / dt   v 



dt d v / dt   f / m 
Slide 9
6.837 Fall 2001
Particle Systems
How would we compute the motion of a particle system
consisting of n particles?
Concatenating the phase space positions of all n particles
yields a large first-order ordinary differential equation. For
particles in 3-D, there will be 6n equations (3 equations for
position and 3 equations for velocity of each particle):
x 11
x 21
x 31
v 11
v 21
v 31
Lecture 22
 x 1   v1 
v   f /m 
 1  1

 x2   v2 
d   


v
f
/
m

dt  2   2
  

  

x
v
n
 n 

 v n   fn / m n 
Slide 10
6.837 Fall 2001
Particle Animation
AnimateParticles(n, y0, t0, tf)
{
y = y0
t = t0
DrawParticles(n, y)
while(t != tf) {
f = ComputeForces(y, t)
dydt = AssembleDerivative(y, f)
{y, t } = ODESolverStep(6n, y, dy/dt)
DrawParticles(n, y)
}
}
Lecture 22
Slide 11
6.837 Fall 2001
Rigid-Body Dynamics
We could compute the motion of a rigid-body by computing
the motion of all constituent particles. However, a rigid body
does not deform and position of few of its particles is sufficient
to determine the state of the body in a phase space. We’ll
start with a special particle located at the body’s center of
mass.
x t 
v t 
Lecture 22
 x t  


? 

y t  
 v t  


 ? 
Slide 12
6.837 Fall 2001
Rigid-Body Equation of Motion
The equation of motion describes the dynamics of a rigid body:
the motion of a body subjected to external forces. We already
know how to write a portion of this coupled first-order
differential equation:
 x t    v t  

 

d
d  ?   ? 
y t  

dt
dt M v t    f t  

 

 ?   ? 
Lecture 22
Slide 13
6.837 Fall 2001
Net Force
f1 t 
f2 t 
x t 
v t 
f t    fi t 
i
Lecture 22
f3 t 
Slide 14
6.837 Fall 2001
Rest of the State
 x t  


? 

y t  
 v t  


 ? 
R t 
ω t 
If we knew the orientation of the body, we could compute the
position of any body point. Let’s represent the orientation with
a rotation matrix R(t). A point pb in body coordinates is
transformed to world coordinates with the following equation:
p t   R t  pb  x t 
The last component is the rate of change of the rotation
matrix, which we’ll call angular velocity ω(t).
Lecture 22
Slide 15
6.837 Fall 2001
Angular Velocity
Angular velocity is a vector ω(t) that encodes both the axis of
rotation (with its direction) and the speed of rotation (with its
speed)*.
ω t 
How are the orientation matrix R(t) and the angular velocity
ω(t) related? The obvious relationship is not sensible:
[matrix]
d
R t   ω t 
dt
[vector]
*This may sound like a quaternion but be careful any vector ω with
magnitude 2π would represent identical orientation. That’s a singularity.
Lecture 22
Slide 16
6.837 Fall 2001
Rotation Matrix
The columns of the rotation matrix describe the world-space
directions of the body-space coordinate axes.
ry t 




R t   rx t  ry t  rz t  


rx t 
rz t 
If we can determine the rate of change of an arbitrary body
vector then we can apply the sample equation to compute the
rate of change of each column in the rotation matrix.
Lecture 22
Slide 17
6.837 Fall 2001
Velocity of Body Vector
ω t 
a
b
r t   a  b
d
r t   ω t   b
dt
 ω t    a  b 
 ω t   r t 
The tip of the vector r moves with the speed of |ω(t)||b| and
the direction of the velocity vector is perpendicular to the
radius b and the axis of rotation ω(t). Therefore, we can write
the velocity vector as the cross product between the angular
velocity ω(t) and the radius vector b. Because the vectors
ω(t) and a are collinear, we can rewrite the velocity as a
function of angular velocity ω(t) and the original vector r(t).
Lecture 22
Slide 18
6.837 Fall 2001
Angular Velocity and Orientation
The orientation rate of change and the angular velocity are
related by applying the same formula on each column of the
rotation matrix. We introduce a special operator to make the
expression shorter to write.




d
d
d
d
R t    rx t 
ry t 
rz t  
dt
dt
dt
 dt









 ω t   rx t  ω t   ry t  ω t   rz t  


 ω t  R t 
Lecture 22
Slide 19
6.837 Fall 2001
Rigid-Body Equation of Motion
Need to relate rate of change of angular velocity to applied
forces. This term must depend on the mass distribution.
 x t   

v t 

 

d
d  R t    t  R t  
y t  


dt
dt M v t   
f t 

 

?
 ω t   

Lecture 22
Slide 20
6.837 Fall 2001
Inertia Tensor
I xx

I t   I yx
 I zx
diagonal terms
I xx  M   y  z
2
V
Lecture 22
2
I xy
I yy
I zy
 dV
Slide 21
I xz 

I yz 
I zz 
off-diagonal terms
I xy  M  xy dV
V
6.837 Fall 2001
Rigid-Body Equation of Motion
 x t   

v t 

 

d
d  R t    t  R t  
y t  


f t 
dt
dt  M v t   

 

τ t 
I t  ω t   

M v t   linear momentum
I t  ω t   angular momentum
Lecture 22
Slide 22
6.837 Fall 2001
Net Torque
f2 t 
f1 t 
p1b t 
x t 
pb2 t 
v t 
pb3 t 


τ t    pbi  x t   fi t 
i
f3 t 
Lecture 22
Slide 23
6.837 Fall 2001
Simulations with Collisions
Simulating motions with collisions requires that we detect them
(collision detection) and fix them (collision response).
y t 3 
y t 2 
y t 1 
y t 0 
Lecture 22
Slide 24
6.837 Fall 2001
Collision Detection
Plane-Particle collision can be detected by computing the
signed-distance function between the boundary and the
particle. In this example, collision detection is deceptively
simple. Fast collision detection between many curved surface
is a difficult problem. Detecting collisions between
bn
deformable surfaces is
bp
p1  b p   bn  0
even harder.
p
2
 b p   bn  0
p
Lecture 22
Slide 25
2
 b p   bn  0
6.837 Fall 2001
Collision Response
The mechanics of collisions are complicated and many
mathematical models have been developed. We’ll just look at
a one simple model, which assumes that when the collision
occurs the two bodes exchange collision impulse
instantaneously.
bn

v
bp

v ?
Lecture 22
Slide 26
6.837 Fall 2001
Frictionless Collision Model

bn  v   bn  v

v



bn
bp
 1
1
 
2
 0
Lecture 22
Slide 27
6.837 Fall 2001
Animation
Animate(y0, t0, tf)
{
y = y0
t = t0
Draw(y)
while(t != tf) {
if (Collision(y)) {
t = CollisionTime(y)
y = ApplyImpulse(y)
}
f = ComputeForces(y, t)
dydt = AssembleDerivative(y, f)
{y, t } = ODESolverStep(y, dy/dt)
Draw(y)
}
}
Lecture 22
Slide 28
6.837 Fall 2001
Animation Design
Suppose we’d like to animate a hat flying through the air and
landing on a coatrack.
We can compute the hat’s motion by techniques described in
this lecture, but this won’t be enough. We also have to
discover what initial position and velocity will result in the hat
landing on the coatrack.
Lecture 22
Slide 29
6.837 Fall 2001
Next Time
Controlling Animation
Lecture 22
Slide 30
6.837 Fall 2001