Class myphysicslab.lab.engine2D.ContactSim

Provided By
All Implemented Interfaces

Physics engine for rigid bodies with contact forces to allow resting contact. The contact forces prevent the bodies from interpenetrating when they are in resting contact. Resting contact means the bodies are not colliding, but have edges and corners that are in continuous contact and exerting force on each other.

The overall idea is to calculate the exact amount of force needed to just barely prevent the objects from penetrating. These contact forces are calculated in the #evaluate method, which is called by the myphysicslab.lab.model.DiffEqSolver at the request of the myphysicslab.lab.model.AdvanceStrategy to advance the state of the simulation.

Background and References

See explanations at:

The algorithm used here is based on these papers:

David Baraff, [Fast Contact Force Computation for Nonpenetrating Rigid Bodies.](Baraff_Fast_Contact_Force_94.pdf) Computer Graphics Proceedings, Annual Conference Series, 1994; pages 23-34. David Baraff, [An Introduction to Physically Based Modeling: Rigid Body Simulation II—Nonpenetration Constraints.](Baraff_Siggraph_97_Course_Notes.pdf) Siggraph '97 Course Notes.

See also the list of David Baraff's papers.

See the paper Curved Edge Physics paper by Erik Neumann for modifications to contact forces when curved edges are involved.

Parameters Created

Find External Forces

Within evaluate() we first let the super-class apply the external forces to the RigidBody objects. The external forces include things like gravity, thrust, springs, damping. The external forces result in accelerations of the bodies, so at this point many of the bodies would start to penetrate into each other if we did not find contact forces to prevent that.

Find Contacts

Next in evaluate() we call findCollisions() to find all the contact points, and possibly collisions as well. If we find any actual (penetrating) collisions, then evaluate() returns the set of collisions found, which should then be handled before trying to step forward again, see myphysicslab.lab.engine2D.ImpulseSim.

The criteria for finding a contact is:

  • the corner (or edge) of one body must be very close to the edge of the other body, as specified by getDistanceTol().

  • the bodies must be moving very slowly relative to each other (the normal velocity) at the contact point, as specified by getVelocityTol().

The evaluate() method optionally finds independent subsets of collisions, because that can make the compute_forces algorithm, which is O(n^4), run faster in some cases. See the flag #SUBSET_COLLISIONS. Collisions are independent when there is no chain of moveable (finite mass) bodies in common between the collisions.

The Matrix Equation for Contact Forces

Now we have the set of contacts and we know the accelerations of the bodies due to external forces. We set up a matrix equation

a = A f + b


  • a = vector of accelerations
  • A = matrix describing how the j-th contact force affects the acceleration of the i-th contact distance (the separation between the bodies)
  • f = vector of contact forces (to be found)
  • b = external forces (gravity, thrust, rubber band, damping)

Set Up the A Matrix

Here is how to set up the A matrix: For each contact distance d_i, find how the acceleration of that contact distance d_i'' is related to the force at the j-th contact point. The force at the j-th contact point is f_j N_j, where f_j is a scalar and N_j is the vector normal. The a_ij entry in the A matrix tells what that relationship is between f_j and d_i''. This a_ij value is dependent only on the current geometry of how the objects are oriented and touching each other.

Baraff Figure 26

For example, consider Figure 26 of Baraffs Siggraph 97 Course Notes shown above. The figure shows two bodies (B,C) resting on the ground, and a third body (A) resting on top of the other two. There are 5 points of contact among the bodies and the ground. Here is how the matrix equation would look for this situation:

a1     a11  a12  a13   0    0     f1     b1
a2     a21  a22  a23   0    0     f2     b2
a3  =  a31  a32  a33  a34   0  *  f3  +  b3
a4      0    0   a43  a44  a45    f4     b4
a5      0    0    0   a54  a55    f5     b5

Consider the first contact at p1 which is between the ground and body B. The acceleration of the contact distance at p1 is affected by the forces at p1, p2, and p3 because all of those forces affect the movement of body B. But the forces at p4 and p5 have no effect, so their entries are zero in the first row of the A matrix.

The first row of the matrix equation can be written out as

a1 = a11*f1 + a12*f2 + a13*f3 + 0*f4 + 0*f5 + b1

That equation says the acceleration of contact distance at point p1 is equal to a certain linear combination of the contact forces f1, f2, f3, plus the acceleration due to the external forces which is b1.

The particular values for a11, a12, a13 are dependent on the geometry of the situation: where the forces f1, f2, f3 are applied on the body; in what direction do the forces act; where is the center of mass of the body; how the force causes the body to accelerate and spin; where is the point p1 is in relation to the center of mass; etc.

The third contact at p3 is more complicated because it is affected by any forces acting on body A or body B. These include all the forces except f5. Therefore the third row of the A matrix has four non-zero entries corresponding to the four forces that affect the acceleration at p3.

Constraints On the Solution

Assume we now have the b vector of external forces and the A matrix of dependencies between contact forces and resulting accelerations of different bodies. We solve for the f vector of contact forces subject to the following constraints:

a >= 0
f >= 0
a.f = 0

The constraints in words are:

  • a >= 0 we require all the relative normal accelerations (between bodies at contact points) to be either zero (remain in contact) or positive (separating).

  • f >= 0 We require contact forces to be zero or positive because they can only push, not pull.

  • a.f = 0 The third constraint is a math way of saying that if there is a force, then acceleration is zero OR if there is acceleration (separation), then there is no force.

Note that for Joints the constraints are different:

  • a == 0 a Joint always remains in resting contact, it never separates

  • No constraint on f because it can push or pull.

Solving For the Contact Forces

The Baraff 'Fast Contact Force' paper goes into full detail about the algorithm to solve this constraint problem. Here is a quick summary:

  • Start by setting all the forces to zero

  • Add in one force at a time, just enough to maintain the constraints on the points that have been considered so far (ignoring the others), and readjust the other forces as necessary.

  • Continue adding one force at a time, ignoring points that haven't been considered yet.

Suppose that we are working on adding in the third force after already finding the forces for points 1 and 2. The trick is that we only look at the constraints on the first 3 contact points, we ignore the other contact points and other forces. Once we've found the 3rd force (and rebalanced forces 1 and 2 as needed) we then move on to consider the 4th force.

The last step is to apply the contact forces to the bodies, which gives the final set of accelerations that the evaluate() method then returns to the differential equation solver.

Extra Acceleration

The contact forces are calculated so that there is zero acceleration at contact points; but this does not immediately affect the remaining small velocity at a contact point. As a result, objects that are in resting contact will often have some undesirable jittery motion.

One way to deal with this is to request a small amount of additional acceleration which will eliminate that velocity over a few time steps. If the objects are moving towards each other (a small negative velocity) we request a little more acceleration which leads to a little more force being applied there. If the objects are moving apart (a small positive velocity) we request a little less acceleration.

The extra acceleration is added to the b vector in the private method calculate_b_vector. See myphysicslab.lab.engine2D.ExtraAccel enum for explanations of the various options. See #setExtraAccel for how to specify the desired ExtraAccel option.

Intermediate Steps During evaluate

A DiffEqSolver works by 'averaging' several calculated states for each time step; for example ModifiedEuler averages 2 states, and RungeKutta averages 4 states. The collision detection done in evaluate() is based on those intermediate states within a step of the DiffEqSolver, so it is arguable whether that state actually ever occurs.

This point of view argues for only using the collisions detected between full steps of the DiffEqSolver. However, contacts can come and go during these sub-steps so it seems in practice to be more accurate to find the set of contacts anew in each call to evaluate(). Also if a penetrating collision is detected we need to stop the process and handle that collision instead, so it is important to do collision detection for that as well.

Prior to February 2012, there was an experimental option specified by the flag REUSE_COLLISIONS for which we did not find the collisions anew in the evaluate() method, instead we used the collisions found outside the evaluate() method during AdvanceStrategy.advance() which are based on a complete ODE step. See the git archive.

new ContactSim( opt_name )


name of this Subject

Instance Methods

Instance Properties

Static Properties