# 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.

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

We find forces such that

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

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);
if (!computeForces.checkForceAccel(1E-8, f, accel, joint)) {
throw '';
}
}

### 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?

## Constructors

• #### Parameters

• ##### name: string

for debugging, this distinguishes whether this is used for contact forces or collision impulses

• ##### pRNG: Random

pseudo random number generator, used to randomly decide order in which to calculate forces

## Properties

name_: string

name of this ComputeForces instance, for debugging only

nextContactPolicy: NextContactPolicy = NextContactPolicy.HYBRID

The next contact policy to use for deciding order in which to treat contacts.

order: number[] = []

Order in which contacts were treated; each entry is index of contact in A matrix, for debugging only.

pRNG: Random

pseudo random number generator, used to randomly decide order in which to calculate forces

preOrder: number[] = []

Order in which to treat contacts; each entry is index of contact in A matrix, when the NextContactPolicy.PRE_ORDERED policy is used.

## Methods

• Calculates the forces at each contact point of a multi-body contact situation. Can also be used to find impulse needed at each collision point in a multi-body collision.

#### Parameters

• ##### A: Float64Array[]

an n x n matrix giving change in acceleration for force at each contact

• ##### f: number[]

force at each contact (vector, length n), this is what is solved for and is returned via this array (this array is zeroed out at start).

• ##### b: number[]

external and inertial forces in the system (vector, length n)

• ##### joint: boolean[]

indicates which contacts are Joints (vector, length n)

• ##### debugCF: boolean

true shows debugging messages

• ##### time: number

the current simulation time, used only for debugging

• ##### Optionaltolerance: number

the tolerance to use for checking the results. If not provided, then no check is done.

#### Returns number

error code, -1 if successful otherwise an error occurred

• Returns order in which contacts were treated; each entry is index of contact in the A matrix.

#### Returns number[]

• Sets the policy for choosing which contact to treat next.

#### Parameters

• ##### OptionalpreOrder: number[]

Order in which to treat contacts; each entry is index of contact in A matrix. Used only with NextContactPolicy.PRE_ORDERED.

#### Returns void

• Returns true if the given force and accel vectors satisfy the constraints that if f != 0 then a = 0, or if f = 0 then a >= 0.

#### Parameters

• ##### tolerance: number

ignore deviations from constraints smaller than this

• ##### force: number[]

array of forces applied to a set of contacts

• ##### accel: number[]

array of accelerations at a set of contacts

• ##### joint: boolean[]

whether each contact is a joint

#### Returns boolean

true if the force and accel vectors satisfy the constraints

• Returns the maximum unwanted acceleration at all contacts.

#### Parameters

• ##### accel: number[]

acceleration at each contact

• ##### joint: boolean[]

true when contact is a Joint

• ##### n: number

number of contacts

#### Returns number

the maximum unwanted acceleration at all contacts

Generated using TypeDoc