Class myphysicslab.lab.engine2D.ContactSim
Provided By  

Extends  
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 math and physics underlying RigidBodySim, ImpulseSim and ContactSim are described on the myPhysicsLab website.

ContactSim Math has more details about the math.
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 2334. 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

ParameterString named
RigidBodySim.en.COLLISION_HANDLING
see#setCollisionHandling

ParameterString named
RigidBodySim.en.EXTRA_ACCEL
see#setExtraAccel

ParameterBoolean named
RigidBodySim.en.SHOW_FORCES
see#setShowForces

ParameterNumber named
RigidBodySim.en.DISTANCE_TOL
see#setDistanceTol

ParameterNumber named
RigidBodySim.en.VELOCITY_TOL
see#setVelocityTol

ParameterNumber named
RigidBodySim.en.COLLISION_ACCURACY
see#setCollisionAccuracy
Find External Forces
Within evaluate()
we first let the superclass 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
where
a =
vector of accelerationsA =
matrix describing how thej
th contact force affects the acceleration of thei
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.
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 nonzero 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 substeps 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 )
Parameters 

