game

physics

2012

dbalster@gmail.com

B E E R E N G I N E

agenda

• introduction
• laws of physics and "some" math
• simulation
• particles
• rigid bodies
• collision detection

required mathematics

vectors and scalars

• real numbers
• vectors (x,y), (x,y,z), (x,y,z,w), ...
• add (accumulate forces)
• sub (distance)
• scale and length
• dot and cross product
• reflect
• project
• linear interpolation
• lerp(a,b,t) = a + (b-a)*t

projection of vector u on line v

v

u

proj_v(u)

• shadow of u on v
• strength of u in direction of v

vector projection

proj_v(u) = (u dot v / v dot v) * v

​

• "v dot v" is one if v is normalized
• perpendicular projection

​

reflection of vector u on v

u

• u is called incident vector
• v can be seen as plane, n is the normal of the plane
• u is reflected on "plane" v

u'

v

n

vector reflection

n = constructed from v

reflect(u,n) = u – 2 ∗ n dot u ∗ n

​

• n must already be normalized in order to achieve the desired result.

​

vector identities

for vectors u, v and w

+ = add, ∙ = scalar mul, ⊙ = dot, ⊗ = cross

• u + v = v + u
• u ⊙ v = v ⊙ u
• u ⊗ v = - v ⊗ u
• (u + v) ⊙ w = u ⊙ w + v ⊙ w
• (u + v) ⊗ w = u ⊗ w + v ⊗ w
• u ⊙ (v ⊗ w) = v ⊙ (w ⊗ u) = w ⊙ (u ⊗ v)
• u ⊗ (v ⊗ w) = (u ⊙ w)∙v - (u ⊙ v)∙w

​

can be used to speed up operations

matrices

• unit matrix
• rotate, translate, scale
• matrix-vector multiplication
• transform vector v with matrix m
• matrix-matrix multiplication
• concatenate transforms
• matrix inversion
• solve linear equation

​

euler angles and rotations

• coordinate system
• naive usage
• rotate around X-axis by alpha etc.
• sin, cos, tan
• degree (0°...360°) and radians (-π...π)
• problem of gimbal lock
• if two rotational axes rotate around the third axis equally, the rotation is not reversible

quaternions

• express rotation around axis in four values ("quat"ernion)
• do not have problem of gimbal lock
• q = w + ix + jy + kz = (w,x,y,z)�i,j,k imaginary numbers with properties�i² = j² = k² = -1�ij = k = -ji, jk = i = -kj, ki = j = -ik
• can be stored in vector4, but it's something completely different

quaterion operations

• add, sub, scale, length like vector4
• conjugation and inversion
• it's simple to inverse, compare with 3x3 matrix
• quaternion-vector multiply
• rotate vector around
• quaternion-quaternion multiply
• concatenation of rotations
• "rotate q1 by q2"
• fewer storage than matrixes
• translate=3 floats, rotate=4 floats
• fewer ops to concatenate
• 4x4 muls vs. 3x3x3 muls

mathematical stuff needed

• equations
• functions
• partial sums Σ
• delta increment Δt
• integral ∫
• difference quotient
• differential calculus
• ...

​

let's start with physics

we have to remember/learn some important laws of physics

some units first

• space
• m in meters
• distance
• area
• volume
• time
• t in seconds

​

newton's laws of motion

• Bodies continue to move in a state of constant velocity unless acted upon by an external force
• The change of momentum of a body is proportional to the impulse impressed on the body: F = m a
• Every force has an equal opposing force

movement

• move an object from p₁ to p₂
• unit is in world coordinates

p₁

p₂

velocity

• displacement over time
• contains both speed and direction�_v= dx� dt
• vector quantity
• scalar length is called speed
• unit of speed is m/s

p₁

p₂

acceleration

• _a=dv _{= F/m}� dt
• change of velocity over time
• also vector quantity
• meters per second per second
• m / s²
• first derivative of velocity(t)
• second derivative of position(t)

​

mass

• SI unit is kg
• definition: m = F/a
• inertia describes how hard it is to move a mass (better: to change its momentum)
• reason for gravity

(linear) momentum

• momentum is mv (mass times velocity)
• SI unit is kg m/s
• vector quantiy
• ...e=mc² is related and public knowledge
• conservation of momentum
• change of momentum is called impulse

​

German!

• Momentum = "Impuls"
• Impulse = "Kraftstoß" or "Impulsänderung"

​

force

• SI unit is newton (N)
• F = ma

​

still simple

• sum of all forces is zero�=> velocity is constant
• F = m a
• acceleration = F/m or change of velocity
• velocity = change of position

​

let's handle time

simulation

simulating time

• computer graphics is done in frames
• animation is defined by frames per second
• each frame gets a small time step dt

​

​

double dt = 1.0 / 60.0;

​

while (true)

{

simulate( dt );

render();

}

​

​

​

double t = 0.0;

​

void simulate( double dt )

{

t += dt; // accumulate time

​

if ( t >= 1.0 ) // each second

{

t = 0.0;

}

​

// full rotation in one second

rotation = sin( t * 2PI );

}

​

​

​

double t = 0.0;

​

void simulate( double dt )

{

t += dt; // accumulate time

​

if ( t >= 1.0 ) // each second

{

t -= 1.0; // !!

}

​

// full rotation in one second

rotation = sin( t * 2PI );

}

​

​

​

double dt = 1.0 / 60.0;

​

while (true)

{

simulate( dt ); (1)

render(); (2)

}

​

• what if (1)+(2) take longer than dt ?
• we could measure the time

​

​

double last = now();

​

while (true)

{

double time = now();

double dt = time - last;

simulate( dt );

render();

​

last = time;

}

• do not forget to start your simulation at t=0
• problem of too small/large deltas

​

​

double last = now();

​

while (true)

{

double time = now();

double dt = time - last;

dt = clamp(dt, 0, 1/60);

simulate( dt );

render();

​

last = time;

}

​

​

simulation

• how do you handle rates
• create N objects per second
• emit N bullets per second
• move at N km/h
• ...
• beware of numerical errors

double t = 0.0;

double rate = 1.0 / 10.0; // items per second;

​

void simulate( double dt )

{

t += dt; // accumulate time

​

while ( t >= 0 )

{

t -= rate;

spawn_item(...);

}

...

}

​

fixed / dynamic time steps

• dynamic time steps good for varying processing speeds
• dynamic time steps adopt to system speed
• upper and lower limits
• dynamic time steps can be converted back into fixed time steps
• beware: large dt will likely increase dt next time
• here: small dt means fast frame rate
• normally you want to detach logic and render
• render may need to wait for vsync or triplebuffer
• concurrent logic and render

let's simulate physics!

// F=ma , dv/dt = a = F/m , dx/dt = v

t = 0

dt = 1 / 60

velocity = 0

position = 0

force = 10 // exact values doesn't really matter here

mass = 1

​

while ( t <= 10 ) // simulate for 10 seconds

{

position += velocity * dt

velocity += (force/mass) * dt

t += dt

}

​

print( simulated position after 10 seconds )

what did we just do?

• we just solved an ordinary differential equation
• we used the numerical integration method of Euler

​

Problems:

• Euler is very inaccurate (ignores curvature etc)
• requires a lot of small, fixed time steps (makes simulation very slow)

some numerical integrators

• Closed form (almost never known)
• Euler
• Midpoint
• Runge-Kutta 4th order (RK4)
• Verlet

​

We need a method that

• few iterations with good accuracy

​

ordinary differential equations

= F(t,X), t is time X(0) = X₀ is state at time 0

"dt" on both sides

=> dX = F( t,X(t) )

rewrite dX as difference quotient

=> ( X(t+h) - X(t) ) / h = F( t,X(t) )

move stuff to right

=> X(t+h) = X(t) + h F( t,X(t) )

​

this is how we're going to solve it numerically

​

dX

dt

Runge-Kutta 4th order method

t₀, X₀ is initial time and state(-vector)

h is our step size

a = h F(t₀,X₀)

b = h F(t₀+½h,X₀+½a)

c = h F(t₀+½h,X₀+½b)

d = h F(t₀+h, X₀+c)

X₁ = X₀ + ¹/₆(a+2(b+c)+d) , t₁ = t₀ + h

Runge-Kutta 4th order method

x_n+1 = x_n + h/6 (a+2b+2c+d)

​

a = f( t_n , x_n )

b = f( t_n + h/2 , x_n + h/2 * a )

c = f( t_n + h/2 , x_n + h/2 * b )

d = f( t_n + h , x_n + h * c )

​

...and this is easy to implement for your f(t,X)

Verlet integration method

• specialized, derived from motion laws

​

x_i+1 = 2x_i - x_i-1 + a * dt²

x_i is your position, a is acceleration, dt is timestep

​

• acceleration is constant
• velocity is implicit
• needs to store last position

solve this problem

A car stands still and starts to accelerate with a m/s² for n seconds. Where is the car after this time?

euler method

// euler integration

​

float euler( float dt, float acceleration, float seconds )

{

float velocity = 0.0f; // m/s

float x = 0.0f; // m

​

for (float t=0; t<seconds; t+=dt)

{

x += velocity * dt;

velocity += acceleration * dt;

}

return x;

}

​

printf( "%f\n", euler( 0.1, 10.0, 10.0 ) );

​

verlet method

// verlet integration

​

void verlet( float dt, float acceleration, float seconds )

{

float x = 0.0f; // m

​

float previous = x;

for (float t=0; t<seconds; t+=dt)

{

float saved = x;

x = 2*x - previous + acceleration * dt*dt;

previous = saved;

}

return x;

}

​

printf( "%f\n", verlet( 0.1, 10.0, 10.0 ) );

​

runge kutta method (1)

​

void rungekutta(float dt, float acceleration, float seconds)

{

Vector X; // abusing Vector as state vector

X.x = 0.0f; // position in m

X.y = 0.0f; // velocity in m/s

Function f(acceleration); // later explained

​

for (float t=0; t<seconds; t+=dt)

{

X = rungekutta4(t, X, dt, f);

}

printf("p=%f, v=%f\n",X.x,X.y);

}

​

runge kutta method (2)

// looks like math in C++

// "Vector" has all required properties

​

Vector rungekutta4(float t, Vector X, float h, Function& f)

{

Vector a,b,c,d;

float h2 = h/2.0f;

a = f( t ,X );

b = f( t+h2 ,X + a*h2 );

c = f( t+h2 ,X + b*h2 );

d = f( t+h ,X + c*h );

return X + (a+(b+c)*2+d)*(h/6.0f);

}

​

​

runge kutta method (3)

// C++ supports "function objects"

​

struct Function

{

float acceleration;

Function(float _acceleration) : acceleration(_acceleration)

{

}

​

Vector operator()( float dt, const Vector& X )

{

return Vector(

X.y, // derivation of position is velocity

acceleration // derivation of velocity is acceleration

);

}

};

​

using numerical solvers

• we're solving differential equations
• we solve in small numerical steps
• each step gives us the chance to "fix"
• check for collisions
• check constraints

​

for (int i=0; i<numsteps; ++i)

{

apply forces, impulses, ...

check constraints/collisions and take action

integrate / solve ODE

}

Constraints

• Pin position
• do not apply changes over time
• Dependent on other object
• hinge
• spring
• joint
• ...

Constraint: hinge

• hinge between two "points"
• distance between points is not allowed to change
• at each change of position, make sure that contraints "distance to other position" is valid
• ie apply change, than correct position

hooke's law ("springs")

F = -kx

x displacement of spring in equilibrium (m)

k spring rate constant (N/m)

​

spring constraint

• two bodies are connected
• connection is our spring
• spring constant
• distance between bodies in equilibrium
• check if hooke's law is valid
• check distance and compare to equilibrium
• apply spring constant

simple collision handling

• penetrating is not allowed
• "if (x < 0) x = 0"
• energy not preserved
• "if (x < 0) x = -x"
• energy is not damped
• time of contact

​

Recall

what are we doing

current list of tools

• we have Newton's laws of motion
• we don't care about friction, damping, ...
• we don't care about mass and rotation, ...
• we know how to solve ODEs
• euler, verlet, runge-kutta
• we can do naive constraints/collisions

​

Let's do it!

~ live demo ~

sphere collision: reflection

contact

contact surface

opposing contact surface normals

sphere collision: handling velocities with vector projection

original velocities

projected velocities

sphere collision: handling velocities with vector projection

original velocities

projected velocities

new velocities

summary of live demo

• our bodies have no mass or extents
• collisions behave strange
• collisions of spheres are easy
• what about arbitrary shapes?
• looks already nice and natural
• perfect for balls, bridges, ropes, cloth, hair
• problem of resting
• problem of contacts
• problem of energy conservation
• impulse, friction, resistance, ... more laws needed

advanced collision detection

sphere, aabb, obb, capsule, ...

collision detection workflow

• start broad then narrow
• for two objects we are interested in
• do they intersect?
• contact points?
• revert penetration
• contact normals (penetration direction)

collision detection problems

• cheap vs expensive algorithm
• data consumption per primitive
• cpu cost per primitive

​

continuous collision detection

• continuous movement makes it very likely that nearby objects will collide next frame

large time step / high velocity

static

static

no collision during steps

collision during steps

bounding volumes

bounding sphere

• defined by radius and center
• point in sphere?
• line passes sphere?
• sphere vs sphere?
• sphere vs plane?
• calculate perfect sphere for object?
• contact calculation?
• transform sphere?

calculate bounding sphere

• we want (radius,center_xyz)
• needs O(n) for N vertices for "bad" solution
• needs O(n²) for better solution
• methods
• calculate AABB and use "outer sphere"
• use "linear programming"

sphere: trivial operations

• point p in sphere s(center,radius)
• |p-s.center| < s.radius
• sphere s1 vs sphere s2
• |s1.center - s2.center| < s1.radius+s2.radius
• sphere s vs plane p(normal,distance)
• distance of s.center from plane < s.radius

​

sphere vs ray

• construct origin to center vector V
• project V onto ray R = V_R
• construct V_c = center-V_R

​

​

​

​

​

hit on |V_C| < radius !

center

origin

V

R₁

R₂

V_R

V_C

contact point between two spheres

n = c₁-c₂

p = (c₁ + r₁n + c₂-r₂n) * 0.5

axis-aligned bounding boxes AABB

• axis alignment make tests very simple
• if p.x < x ...
• store as location and extent
• center and radius(x,y,z,...)
• min and max location

point in AABB

given point P (x,y,z) and AABB A(min,max)

​

P inside A if

​

P>A.min and P<A.max

​

construct AABB

vector min = (flt_max,...)

vector max = (flt_min,...)

for each p in object :

if p.x < min.x : min.x = p.x

if p.y < min.y ...

...

if p.x > max.x : max.x = p.x

...

construct AABB (2)

• AABB is in world coordinates
• construction is expensive
• if object moves...
• apply move to AABB.center
• if object scales...
• apply scale to AABB.radius
• if object rotates...
• rotate both center and radius

​

so AABB reconstruction can be avoided during scenegraph transform

merging two AABBs

• if you group multiple objects (scenegraph) you need a combined AABB
• merging is fast
• get new min/max from all extents
• place new center

​

vector min = min( a.min, b.min )

vector max = max( a.max, b.max )

vector center = min+(max-min)/2

intersection of two AABB

min₁

min₁

min₂

min₂

max₂

max₁

max₂

max₁

intersection of two AABB (2)

​

bool aabb3_intersects( ABB3 a, AABB3 b )

{

var r[3];

r[0] = a.r[0] + b.r[0];

r[1] = a.r[1] + b.r[1];

r[2] = a.r[2] + b.r[2];

return

((a.c[0]-b.c[0] + r[0]) < 2*r[0]) &&

((a.c[1]-b.c[1] + r[1]) < 2*r[1]) &&

((a.c[2]-b.c[2] + r[2]) < 2*r[2]);

}

​

oriented bounding boxes

OBB

• may fit much tighter
• bounding box rotates with�object coordinate system
• => OBB transform cheap
• calculating perfect OBB is�expensive, only needed if mesh is�changing

​

intersection of two OBB

a on p

b on p

c on p

c

a

a

plane p

project box onto a plane and make 1D check

capsules

• cylinder with caps (two half-spheres)
• ideal for many things
• point in capsule
• check in cylinder and two spheres
• transform capsule points only

r

calculate capsules

• calculate bounding box
• use longest axis for length
• next longest axis as radius
• store line segment (two points) and radius

capsule-capsule collision

• compute distance d between lines
• collision if d < r₁+r₂
• faster than OBB!

r₁

r₂

déjà-vu

time step problem

• capsule is a swept-sphere
• sweeping good to avoid lost collisions

r₁

r₂

collision detection

• N objects have to be tested against each other = O(n²)
• use spatial partitioning to reduce tests
• pair (a,b) same as pair (b,a) (mark pairs)

​

bounding volume hierarchies (BVH)

• binary space partitioning (BSP)
• divide space in equal halfs, supports binary search
• quad- and oc-tree for 2D and 3D
• problem of unbalanced trees
• balanced trees
• kd-tree

collision detection process

• find all colliding pairs
• reduce: a.collides(b) == b.collides(a) !
• test only moving objects (not static vs static)
• use hierarchy of bounding volumes
• analyse pairs with collision
• contacts
• normal of contacts
• penetration (needs to be fixed by integrator)

compound bounding volumes

• if object doesn't fit into basic BV primitive, try to make composite of basic primities

​

rigid bodies

adding mass and rotation

rigid body dynamics

• so far our objects were particles
• ignored mass
• no extent in space
• integrated only position, not orientation
• rigid body means
• equally distributed mass
• non-deformable (does not break, bend, ...)
• we're now going to
• add mass
• add rotation
• update our simulation

​

difference particle and rigid body

• particle
• has three degrees of freedom
• only laws of linear motion
• rigid body
• has linear and angular motion
• makes 3+3=6 degrees of freedom
• can rotate around itself
• Chasles' theorem
• you can split motion into linear and angular part
• rigid body motion is translation along a line and rotation about that line (in and order)

mass, density, volume

F = ma <=> m = F/a

1kg = 1 N/m/s²

​

weight ≠ mass !

​

density = mass per volume

​

density: δ = dm/dv <=>

dm = δ dv => mass = ∫dm = ∫δ dv

​

numerical method for finding mass

• you have a rigid body with unknown volume
• create bounding volume with known closed volume formula (ie sphere, box, ..)
• sample a lot of random points in known volume and test if point is in rigid body
• accumulate all positives

m₁

m₃

m₂

m₄

F₁

F₃

F₂

F₄

(length of force vectors is acceleration)

system of masses

m = Σ_i=1..∞ m_i = m₁+m₂+m₃...

Newton's 2nd law: F=ma <=> m=F/a

​

Σ_i=1..∞ m_i = Σ_i=1..∞F_i/a_i

velocity and acceleration

velocity (v) is

1st derivative of position over time in m/s

acceleration (a) is

1st derivative of velocity over time

2nd derivative of position over time over time in m/s^2�

​

center of mass / barycenter

• weighted average location of all masses in a body
• C = 1/M Σ_i=1..∞ p_im_i , where�M = mass of body�p_i = i-th position of volume element�m_i = i-th mass of volume element
• for rigid bodies
• C is a fixed position
• C must not lie inside of body (!)�

finding center of mass

• find weighted average of mass centers:
• take n samples of body in a bounding volume
• if inside: add sample location with 1/M
• if outside: add sample location with -1/M

applying force to a rigid body

if you apply a force to a point of the body, which is incident to the center of mass:

• you need to overcome inertia (newton's 1st law)
• the body will start to move into that direction (newton's 2nd law)
• you may want to distribute the force to the collider (newton's 3rd law)

​

=> so far we dealt only with linear motion

linear and angular motion

• position
• linear velocity
• linear acceleration
• linear momentum
• ...

• rotation (a,axis)
• angular velocity
• angular acceleration
• angular momentum
• ...

angular velocity

w(t) is called angular velocity

w

P

Origin

w x r

r

angular acceleration

a(t) = dw(t) / dt

situation so far

center of mass

angular velocity

linear velocity

rigid body simulation skeleton

• integrate position
• position += velocity * dt
• integrate linear velocity
• velocity += acceleration * dt
• integrate angular position
• angular position += angular velocity * dt ?
• integrate angular velocity
• angular velocity += angular acceleration * dt ?
• find collision contacts
• correct contacts (apply constraints etc)

integration

• if body "fixed", make infinite mass
• update linear velocity only for finite masses
• newton's first law: state resting
• gravity force needs to be integrated every tick
• angular velocity only at collision time
• updating position and orientation
• position += velocity * dt
• and for orientation...?

integrate angular position (1)

"angular position += angular velocity * dt"

• "angular position" means the orientation
• can be expressed as matrix or quaternion
• we choose quaternion (much simpler)
• "angular velocity" is first derivative of rotation (quaternion or matrix)
• so angular velocity is a vector quantity (delta of rotations per axis: w(t) = (x,y,z)
• we need to differentiate d q(t)/dt (quaternions) or d R(t)/dt (for rotation matrices)
• we need to convert angular velocity into quaternion/matrix

derivative of dq(t)/dt (quaternions)

dq(t)/dt = 0.5 w(t) q(t) where

​

w(t) is angular velocity as quaternion (over time

q(t) is our current (t) rotation as quaternion

​

we could call this "spin"

integrate angular position (2)

• represent angular velocity as quaternion w (with scalar part w)
• convert angular position to quaternion q ie from axis-angle or rotation matrix

q' = 0.5 w * q * dt

We skipped a lot (!) of math, proofs and background, but the result is obviously easy to implement!

(Note: w should be normalized after this)

integrate angular velocity (1)

"angular velocity += angular acceleration * dt"

​

angular acceleration a:

​

a = dw(t)/dt (first derivative of ang. velocity)

= d²q(t)/dt² (second derivative or rotation)

​

...but following this path will get more and more complex (no explanation here).

​

​

angular acceleration (2) in R³

Newton's second law (the one about inertia)

​

F = ma (for linear motion), where

F=force, m=mass, a=lin. acceleration

​

L = I w (for angular motion), where

L = angular moment, I = inertia, w=ang. velocity. Another equation says

torque = I α, where α is ang. acceleration.

​

​

mass properties of a rigid body

• we know the material density
• we have to calculate/measure the volume
• mass is density * volume
• we have to calculate the center of mass
• we need the so called "inertia" or moment of inertia tensor for the rotational part of motion

moment of inertia tensor

• tensors are in our case scalars, vectors or matrices.
• the moment of inertia (newton's 1st law!) describes the resistance of a body against angular change
• only a 3x3 matrix is sufficient to describe this with enough information
• moment of inertia has the same role like mass in linear motion

calculating moment of inertia

• m_k are N random point masses of the rigid body

changing rotation

Rotation won't change unless changed.

(Newton's 1st law, the one about inertia)

​

It only changes, if we apply an impulse to it (through collision response)!

Applying an impulse to the system

center of mass r

linear velocity v

ang. velocity w

Impulse J

contact c

body has mass m and moment of inertia I

impulse: position

• we use a very (!) simplified model of an impulse
• our impulse will just change velocities
• change in velocity is change in momentum
• momentum = velocity * mass

​

So, we have a mass m and a force J

​

=> velocity += J / m (F=ma !)

impulse: rotation

• moment of inertia for rotations has the same role as mass for linear motion
• the new rotation axis�is defined by the�cross product�between J and r

​

​

​

(this is 2D, correct axis is screen z)

r

J

J ⊗ r

impulse: rotation (2)

w = angular velocity

​

new w = old w + (r ⊗ J) * I^-1

​

​

summary

what have we learned

• newton's laws of motion
• basic simulation of time
• numerical integration of ODE�(ordinary differential equations)
• particle systems with springs
• rigid body basics

what's still missing

• loosing energy
• damping, friction, resistance, dragging, etc.
• more collision detection
• contact graph
• shock propagation
• ...
• more joint constraints
• hinge, spring, fixed

​

and much much more... have fun learning!

I hope you enjoyed it.