01 Dec 2019
part 1
part 2
part 3
part 4
part 5
part 6

In order to handle collisions and resting contacts, it is necessary to determine contact points.
One can achieve this using the separating plane algorithm (also see Baraff).

## Separating Planes Algorithm

The separating plane algorithm only applies to convex polyhedra.
However non-convex objects can in many cases be approximated by a union of convex polyhedra.
A separating plane is a plane such that all points of the polyhedra A and B lie on opposite sides of the plane.
If two polyhedra are not intersecting, a separating plane can either be defined using one of the faces of one of the polyhedra or
be constructed from one edge of A and one edge of B.

The distance *d* of a point *x* from the plane *P(p,n)* defined by the point *p* and the normal vector *n* is:

In the case of using a face from polyhedron A, the point *p* is a corner of the face.
The normal vector *n* is the normalized cross product of two edges of the face and pointing outward from the polyhedron.
The point *x* of polyhedron B with the smallest distance *d* defines the distance of B to the face of A.
If this distance is positive, the face of polyhedron A is a separating plane.

In the case of using two edges, the point *p* is the corner of an edge from polyhedron A.
The normal vector *n* is the positive or negative normalized cross product of the two edge vectors.
Again the point *x* of polyhedron B with the smallest distance *d* defines the distance of B to the face of A.

From all the planes, the one with the largest distance *d* is selected.
If the distance is zero or negative, the two polyhedra are in contact.
One can place a plane in the "middle" of the contact by moving it by half the distance *d*, i.e.

In the next step all points from A and B, which have a distance less than *ε* to the plane, are selected.
Two vectors *u* and *v* orthogonal to the normal *n* and orthogonal to each other are determined.
The points are projected onto the plane and represented using the vectors *u* and *v*:
The convex hull of the each of the two polygons is determined.
Then the intersection of the resulting two 2D polygons is computed.
Finally the convex hull of the intersection is taken.
The points are projected back into 3D space.
These are the contact points used for simulating the colliding and resting contacts.

29 Nov 2019
part 1
part 2
part 3
part 4
part 5
part 6

The following article is partially based on Hubert Eichner's article on inequality constraints.

## Inequality Constraints for Colliding Contacts

At a colliding contact the relative speed of the two objects is negative.
I.e. a contact is a colliding (and not a resting) contact if the normal speed at the contact point is below a small negative threshold.
Colliding contacts have restitution, i.e. the collision is partially elastic and the two objects will separate after the collision.
The colliding contacts and affected joints are handled using a special time step of zero duration.
Again the inequality constraint is
The linear and angular components of *J* are the same as for a resting contact.
The correction term *b* however depends on the initial normal speed at the contact point *vn* and the restitution coefficient *ε*:

## Sequential Impulses for Collisions

Again the contact and joint impulses are estimated iteratively.
Note that the normal speed *vn* at each contact is determined beforehand and stays constant during the iterations.

- determine
*vn* for all colliding contacts
- for each iteration
- for each joint
- compute Jacobian
*J* and correction vector *b*
- predict speed
*u*
- compute
*λ*
- compute the impulse
*P*
- add
*P* to accumulated impulses of the two objects

- for each colliding contact
- subtract old impulses
*P* from previous iteration from accumulated impulses of the two objects
- compute Jacobian
*J* and correction vector *b*
- predict speed
*u*
- compute new
*λ* and clip it to be greater or equal to zero
- compute the impulse
*P*
- add
*P* to accumulated impulses of the two objects

- apply impulses to objects

The value *λ* is stored as *Pn* for limiting friction impulses lateron.
The time step here is zero!
Therefore external forces do not need to be considered while handling collisions.
An object falling to the floor will experience several collisions until the linear and angular speed has decreased sufficiently for the contacts to become resting contacts.

25 Nov 2019
part 1
part 2
part 3
part 4
part 5
part 6

The following article is based on Hubert Eichner's article on inequality constraints.

## Inequality Constraints for Resting Contacts

Contact points (resting contacts) are represented as inequality constraints.
In contrast to a joint, a resting contact can only create repellent forces and no attracting forces.
Similar to a joint constraint, the resting contact is represented using a function *C(y(t))*.
Here *C* is one-dimensional and it is basically the distance of the two objects at the contact point.
The inequality constraint of the resting contact is
with the matrix *J* having one row and twelve columns.
Instead of anchor points, one uses vectors *ri* and *rj* from the center of each object to the contact point.
Using the contact normal *n* the inequality constraint becomes:
The rotational component can be simplified as shown below:
Thus the linear components of *J* are
And the angular components of *J* are
The correction term depends on the distance *d*.
I.e. if the distance is negative, a correction is required so that the objects do not penetrate each other any more.

## Sequential Impulses (updated)

The impulses generated by contact constraints and joint constraints are accumulated in a loop.
The (cumulative) *λ* obtained for contact constraint furthermore is clipped to be greater or equal to zero.
Note that it has to be the cumulative value which is clipped.
The algorithm for solving the joint and contact constraints becomes:

- for each iteration
- for each joint
- compute Jacobian
*J* and correction vector *b*
- predict speed
*u*
- compute
*λ*
- compute the impulse
*P*
- add
*P* to accumulated impulses of the two objects

- for each resting contact
- subtract old impulses
*P* from previous iteration from accumulated impulses of the two objects
- compute Jacobian
*J* and correction vector *b*
- predict speed
*u*
- compute new
*λ* and clip it to be greater or equal to zero
- compute the impulse
*P*
- add
*P* to accumulated impulses of the two objects

- use impulses and external forces to perform Runge Kutta integration step

The sequential impulse iterations have to be performed four times when using the Runge-Kutta method.
The value *λ* is stored as *Pn* for limiting friction impulses lateron.

24 Oct 2019
part 1
part 2
part 3
part 4
part 5
part 6

Inspired by Hubert Eichner's articles on game physics I decided to write about **rigid body dynamics** as well.
My motivation is to create a space flight simulator.
Hopefully it will be possible to simulate a space ship with landing gear and docking adapter using multiple rigid bodies.

## The Newton-Euler equation

Each rigid body has a center with the 3D position *c*.
The orientation of the object can be described using a quaternion *q* which has four components.
The speed of the object is a 3D vector *v*.
Finally the 3D vector *ω* specifies the current rotational speed of the rigid body.

The time derivatives (indicated by a dot on top of the variable name) of the position and the orientation are (Sauer et al.):
The multiplication is a quaternion multiplication where the rotational speed *ω* is converted to a quaternion:

In a given time step *Δt* the position changes as follows

The orientation *q* changes as shown below
Note that *ω* again is converted to a quaternion.

Given the mass *m* and the sum of external forces *F* the speed changes like this
In other words this means

Finally given the rotated inertia matrix *I(t)* and the overall torque *τ* one can determine the change of rotational speed
Or written as a differential equation

These are the Newton-Euler equations.
The "x" is the vector cross product.
Note that even if the overall torque *τ* is zero, the rotation vector still can change over time if the inertial matrix *I* has different eigenvalues.

The rotated inertia matrix is obtained by converting the quaternion *q* to a rotation matrix *R* and using it to rotate *I₀*:

The three column vectors of *R* can be determined by performing quaternion rotations on the unit vectors.
Note that the vectors get converted to quaternion and back implicitely.

According to David Hammen, the Newton-Euler equation can be used unmodified in the world inertial frame.

## The Runge-Kutta Method

The different properties of the rigid body can be stacked in a state vector *y* as follows.
The equations above then can be brought into the following form
Using *h=Δt* the numerical integration for a time step can be performed using the Runge-Kutta method:

The Runge-Kutta algorithm can also be formulated using a function which instead of derivatives returns infitesimal changes

The Runge-Kutta formula then becomes

Using time-stepping with *F=0* and *τ=0* one can simulate an object tumbling in space.