Class ComputeForces

Computes forces at contact points between RigidBodys, or impulses at collision points between RigidBodys. The compute_forces method is an implementation of the algorithm given in the paper

• Fast contact force computation for nonpenetrating rigid bodies by David Baraff, Computer Graphics Proceedings, Annual Conference Series: 23-34, 1994. 12 pages.

More info at:

• 2D Physics Engine Overview

• The Math and Physics Underlying ContactSim

• ContactSim

• ImpulseSim

This documentation is written assuming that contact forces and resulting accelerations are being calculated, but everything applies equally when calculating multiple simultaneous collision impulses and resulting velocities.

Terminology

This documentation uses several terms from the Baraff paper, such as C, NC, Acc. See that paper for precise definitions. Roughly these are:

• C is the set of contacts that have some force applied and zero acceleration

• NC is the set of contacts that have no force applied because they are separating (they have negative acceleration)

• Acc is the subset of the A matrix corresponding to just the the set of contacts C

The algorithm starts with both C and NC being empty. We then examine one contact at a time, moving it into C or NC, and possibly moving existing contacts between C and NC as necessary.

Constraints

The acceleration of the gap distance at each contact point is given by

a = A f + b
a = acceleration of the gap (positive means gap is tending to widen)
A = n x n matrix where A[i,j] = change in acceleration at contact i resulting from
  force of 1 being applied at contact j
f = force applied at each contact (vector, length n)
b = external and inertial forces in the system (vector, length n)
n = number of contacts
Copy

We find forces such that

0 <= A f + b
f . a = 0
Copy

The first condition ensures that acceleration is non-negative, so that bodies are not accelerating into each other: contact forces can only push, never pull. The second condition means that either

• the force at a contact is zero and acceleration is positive (contacts are separating), or
• the force at a contact is positive and acceleration is zero

Joints

Joints are contact points that can both push and pull, and which never break their contact. Regular contact points can only push, and the contact is broken if the objects move apart. Joints are called 'bilateral constraints' in the Baraff paper.

Joints have different constraints: the force can be be positive or negative, but the acceleration is always exactly zero.

Return Value

The return value from compute_forces is -1 if successful, meaning that a set of forces were found so that the acceleration satisfies the above constraints. If not successful, the caller can check the set of forces that were calculated to see if the resulting acceleration at each contact point is still acceptable.

See the method checkForceAccel for how to check the accelerations. For example, the following code calculates and checks the acceleration from the calculated forces.

const error = computeForces.compute_forces(A, f, b, joint, false, time);
if (error != -1) {
  let accel = UtilEngine.matrixMultiply(A, f);
  accel = UtilEngine.vectorAdd(accel, b);
  if (!computeForces.checkForceAccel(1E-8, f, accel, joint)) {
    throw '';
  }
}
Copy

Redundant Contacts

It can often happen that the set of contacts is redundant in that pushing on one contact leads the algorithm to increase forces at other contacts in such a way that the push at that contact is negated, leading to no change in acceleration at that contact. This means that we cannot independently set the acceleration at that contact.

This shows up mathematically as the Acc matrix being 'poorly conditioned' or singular if we were to add that contact to C, with some row being a linear combination of other rows. (Condition number is a measure of how close a matrix is to being singular). Each row of the matrix corresponds to a contact point. Another way to say the same thing is that there is a linear combination of rows that equals the zero vector. Then any one of those rows is 'redundant' (because it can be expressed as a linear combination of other rows).

Simulations where redundant contacts show up include: the Pile simulation where a pile of rectangular blocks is resting on each other in a corner on the ground; and the DoNothingGrinder where the shuttle blocks are wedged between immoveable blocks.

Deferred (Rejected) Contacts

We avoid adding redundant contacts to C, to keep the Acc matrix non-singular as long as possible. When starting to drive a contact to zero acceleration, we first check to see if adding this contact to C will make Acc singular. If so, we defer (also called reject) the contact, so the contact goes on a list of rejects (the R array) which are handled only after all other non-deferred contacts have been treated, and then only if the acceleration at the deferred contact is large enough to worry about.

For the normal non-deferred contacts we have a limit on the acceptable acceleration of SMALL_POSITIVE; but for deferred contacts, we have a different larger limit on the acceptable acceleration of 100 times that amount. Here 'acceptable' means whether the overall solution is acceptable, so that the compute_forces method can indicate success in its return value.

There are some other ways for a contact to be 'deferred' (or 'rejected'). One is when moving a contact from NC to C. When that happens, we do the same kind of check of whether the contact will make Acc singular, and if so we defer the contact.

Another way that a contact can be deferred is if we notice that a 'zero step' was made repeatedly at the same contact. Often when a contact moves from NC to C, the step size is close to zero because the contact had zero acceleration. This is OK, but if we then notice that the contact moves back from C to NC with another zero size step, we defer that contact because this can lead to an infinite loop. Note that there can be intervening zero steps of other contacts; for example, contact A, then B, and then C all move from NC to C with zero size steps, then A moves back from C to NC with a zero step -- we would defer contact A. But if any other the intervening steps (for B and C) were non-zero size then we would not defer contact A.

Note that we can only defer a contact when it has zero force on it in the solution as calculated to date. This is because all contacts in C usually have some non-zero force, and if you removed a contact from C without first reducing its force to zero then the solution would no longer be balanced and the acceleration at other contacts in C would no longer be zero and contacts in NC might have negative acceleration. Therefore we can defer a contact d before starting to drive it to zero acceleration, because it is not yet in C and has no force. But as soon as you start driving to zero, you have committed to putting the contact d into C because each step increases the force at d. We can defer any contact that is currently in NC because it has no force. In the 'zero step' case, we can defer the contact that is in C only if it has zero force on it.

Order of Treating Contacts

The order in which we handle (or 'treat') contacts is important and can affect what solution is found. The policy is set via the setNextContactPolicy method. The default policy is NextContactPolicy.HYBRID which first treats Joints in random order, and then non-Joints in the order defined by which has the most negative acceleration.

There are three other contact order policies: NextContactPolicy.MIN_ACCEL, NextContactPolicy.RANDOM, NextContactPolicy.PRE_ORDERED. Some of these are used for testing.

Infinite Loop Detection

There is a mechanism to detect infinite loops where a series of contacts keeps being rejected over and over. Part of the mechanism looks at whether any progress was made in the latest step by seeing if the acceleration at the contacts has changed.

The details of the infinite loop detection are as follows: There is a second 'reRejects' list which contains twice-rejected contacts. If we try to treat a reject, but then reject it again, it goes into the reRejects list and is removed from the rejects list. When any progress is made, the reRejects go back to the rejects list to be treated again. If the rejects list is exhausted without making any progress, then an infinite loop is detected, then we abandon the entire process, returning an error code. It is then up to the caller to decide if the resulting solution is adequate or not.

Sometimes Acc Becomes Singular

Despite the effort to keep Acc non-singular, we sometimes need to treat a contact that will make Acc singular because the contact has acceleration that is unacceptably large. In most cases this 'unacceptably large' acceleration is actually very small, like 1E-8 where the limit is 1E-10.

This algorithm is able to still find a solution when Acc is singular, but then the forces can become unreasonably large in order to drive the acceleration to a small value. What seems to often happen is the following: we are driving contact d to zero even though it makes Acc singular (if d were added to C) – this happens when we are treating a previously deferred contact, and is towards the end of the process when all non-deferred contacts are in the solution. What usually happens is that some other contact immediately moves from C to NC and then the Acc+d matrix becomes non-singular, which is a good result.

This algorithm is able to find a solution as long as the b vector is in the column space of the A matrix. This shows up in two places: first, we use a method of solving the matrix problem A x = b that can deal with a singular matrix like this. Second, we will see that when trying to drive a 'redundant' contact to zero acceleration that the delta_a (the change in acceleration from applying force at the contact) is zero; normally this means that we cannot drive that contact to zero acceleration and would fail; but instead it typically is the case that the total acceleration at that contact is already zero (or close to zero) because we have driven the other contacts to zero, and the redundant contact is dependent on those.

Will Not Find Minimal Forces

This algorithm is not guaranteed to find the minimum set of forces that will satisfy the constraints. Rather, the solution found (the set of forces) is sensitive to the order in which contacts are treated.

See myphysicslab.lab.engine2D.test.UtilEngine_test for unit tests that use random contact orderings; those tests show that the maximum force and the length of the force vector depends on the ordering, and also on the criteria for when a matrix is poorly conditioned (which affects when we defer treating a contact that would make the Acc matrix poorly conditioned).

Performance Tweaks

ComputeForces keeps a matrix allocated that is reused, to avoid re-allocating a large matrix which seems to be a performance bottleneck. Some of the matrix algorithms were modified so that the matrix can be over-sized. Also we reuse several vectors between calls to compute_forces, to avoid reallocation.

This resulted in an 11% reduction in running time for the pile_20_random_blocks performance test.

Future Improvements

See Future Improvements in 2D Physics Engine Overview.

Remaining Mysteries

There are two remaining 'math mysteries' in compute_forces:

1. When we go to drive a redundant contact d to zero, it pushes a contact from C to NC, so that d can be added and Acc stays non-singular. Is that guaranteed somehow?

2. When unreasonably large forces are calculated, it looks like it’s usually because there is a pair of opposed contacts and somehow choosing the order so that these are treated close together (in the sequence of which contact to treat when) is what causes a large force to occur. Is there a way to recognize this and avoid it? Perhaps the two contacts are close to linearly dependent? Or maybe adding the second makes the condition of Acc bad?