14. Contact: distance and collision

Chapter 11 measured the gap between two bodies at one instant with distance_to. A moving-body code advances its bodies by a time step and needs more: the instant inside the step at which two of them first touch. This chapter throws a small cube, the pebble, at the dragon. Four runs build the answer step by step:

1. sample the gap over the step and see what sampling misses;
2. advance the pebble by conservative steps, each certified by the gap it starts from, to the time of first contact;
3. find the earliest contact among several bodies at once;
4. do the same in 2-D around the NACA 0012 section, then slide a particle along a wall, where the plain bound creeps, the directional bound does not, and the default switches from one to the other.

use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle, pose_interpolate
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ez_R8P

implicit none
real(R8P), parameter     :: PI = acos(-1._R8P)
type(surface_stl_object) :: dragon, pebble
type(rigid_pose_t)       :: p0, p1
type(vector_R8P)         :: bmin, bmax
real(R8P)                :: scale, d, s, first, last
integer(I4P)             :: k, touching
logical                  :: hit, was

call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)          ! the dragon, placed as in chapter 3
bmin = dragon%get_bmin() ; bmax = dragon%get_bmax()
call dragon%translate(delta=-0.5_R8P * (bmin + bmax))
call dragon%rotate(axis=ez_R8P, angle=PI / 2._R8P)
call dragon%analyze
bmin = dragon%get_bmin() ; bmax = dragon%get_bmax()
scale = 1._R8P / maxval([bmax%x - bmin%x, bmax%y - bmin%y, bmax%z - bmin%z])
call dragon%resize(x=scale, y=scale, z=scale)
call dragon%analyze

call pebble%load_from_file(file_name='cube-01.stl', guess_format=.true.)         ! a pebble: a cube of side 0.04
call pebble%translate(x=-0.5_R8P, y=-0.5_R8P, z=-0.5_R8P)                       ! centred on its own origin
call pebble%resize(x=0.04_R8P, y=0.04_R8P, z=0.04_R8P)
call pebble%analyze
! one step: thrown along x through the dragon's body, turning a quarter turn about y
p0 = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=0._R8P, translation=[-0.6_R8P, -0.3_R8P, -0.2_R8P])
p1 = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=PI / 2, translation=[0.6_R8P, -0.3_R8P, -0.2_R8P])

• The dragon is placed as in chapter 3, without the repair: distance_to and the collision queries take any triangulated surface, open or closed. The pebble is cube-01.stl scaled to side 0.04 and centred on its own origin. A body turns about the origin of its own frame, so a pebble centred there spins in place.
• The step takes the pebble from the pose p0 at s = 0 to p1 at s = 1: 1.2 along x, straight through the dragon's body, while it turns a quarter turn about y. Every run of the chapter uses this setup.

Sampling misses the contact

do k = 0, 4                                                                      ! the gap at five instants
   s = 0.25_R8P * k
   call dragon%distance_to(pebble, distance=d, intersecting=hit, pose_other=pose_interpolate(p0, p1, s))
   print '(A,F5.2,A,F8.5,A,L1)', 's = ', s, ': gap ', d, ', touching ', hit
enddo

• pose_interpolate(p0, p1, s) is the pose at s. All the step queries use this motion: the rotation turns about the axis through the body origin by s θ, and the translation is linear in s (see continuous collision detection).
• distance_to(..., intersecting=, pose_other=) gives the gap at that pose with the trees as built (chapter 11).

touching = 0 ; was = .false.                                                     ! 200 instants, 1/199 apart
do k = 0, 199
   s = k / 199._R8P
   call dragon%distance_to(pebble, distance=d, intersecting=hit, pose_other=pose_interpolate(p0, p1, s))
   if (hit) touching = touching + 1
   if (hit .and. .not. was) first = s                                            ! an interval of contact starts
   if (was .and. .not. hit) print '(A,F7.4,A,F7.4)', '  touching from s = ', first, ' to s = ', last
   if (hit) last = s
   was = hit
enddo
print '(A,I0,A)', '200 samples: touching at ', touching, ' of them'

$ tut14a_sampling
s =  0.00: gap  0.41853, touching F
s =  0.25: gap  0.12466, touching F
s =  0.50: gap  0.05901, touching F
s =  0.75: gap  0.07438, touching F
s =  1.00: gap  0.35909, touching F
  touching from s =  0.3568 to s =  0.4171
  touching from s =  0.5578 to s =  0.6080
200 samples: touching at 24 of them

• The five samples print gaps from 0.059 to 0.419 and no contact. At s = 0.5 the pebble is inside the dragon, 0.059 from its surface. It entered the dragon between s = 0.25 and 0.75 and left again, and no sample saw it.
• 200 samples, 1/199 apart, find two intervals of contact: from s = 0.3568 to 0.4171 and from 0.5578 to 0.6080. Each tells only that the contact began within the last spacing. A body thinner than its travel between two samples, here 0.006, can still pass through unseen. More samples shrink the gap in time, but they never close it.

Conservative advancement

The guarantee comes from a speed bound. No point of the pebble moves faster than μ (distance per unit of s), so from a gap d the two bodies cannot touch before s + d/μ. Conservative advancement steps there, measures the gap again and repeats until the gap falls below a tolerance.

n = 3 * pebble%get_facets_number()
allocate(x(n), y(n), z(n))
do i = 1, pebble%get_facets_number()                                              ! the pebble's vertices (body frame)
   f => pebble%facet_at(i)
   x(3 * i - 2:3 * i) = f%vertex%x ; y(3 * i - 2:3 * i) = f%vertex%y ; z(3 * i - 2:3 * i) = f%vertex%z
enddo
mu = motion_bound_3d(p0, p1, x, y, z)                                             ! |dt| + |theta| max rho
da = dragon%get_bmax() - dragon%get_bmin() ; db = pebble%get_bmax() - pebble%get_bmin()
tol = 1.e-6_R8P * max(da%normL2(), db%normL2())                                   ! the default tolerance
print '(A,F7.4,A,F7.4,A,ES10.3)', 'pebble speed bound mu = ', mu, ' (translation 1.2, spin ', mu - 1.2_R8P, &
                                  '), tolerance ', tol

• motion_bound_3d(p0, p1, x, y, z) is the bound of the motion of pose_interpolate: μ = |Δt| + |θ| max ρ, with ρ the distance of a vertex (body frame) from the rotation axis. The pebble's μ is 1.2444: 1.2 of translation and 0.0444 of spin, a quarter turn at ρ = 0.028. A body at rest has μ = 0, and the two bounds add.
• The tolerance defaults to 1e-6 of the larger bounding-box diagonal of the two bodies, here 1.303e-6.

call dragon%distance_to(pebble, distance=d, pose_other=p0)                        ! s = 0
print '(A,F9.6,A,ES10.3,A,ES10.3)', 'step 0: s = ', 0._R8P, ', gap ', d, ', next step (d - tol/2) / mu = ', &
                                    (d - 0.5_R8P * tol) / mu
do k = 1, 5                                                                       ! stop after k steps: s_k, d(s_k)
   call dragon%time_of_impact(pebble, pose_identity(), pose_identity(), p0, p1, s, status, distance=d, &
                              max_iterations=k)
   if (d > tol) then
      print '(A,I0,A,F9.6,A,ES10.3,A,ES10.3)', 'step ', k, ': s = ', s, ', gap ', d, ', next step (d - tol/2) / mu = ', &
                                               (d - 0.5_R8P * tol) / mu
   else
      print '(A,I0,A,F9.6,A,ES10.3,A)', 'step ', k, ': s = ', s, ', gap ', d, ' <= tolerance: contact'
   endif
enddo

• max_iterations=k stops time_of_impact after k steps. Its toi is then the k-th iterate s_k (a certified lower bound on the contact) and distance the gap there. The loop prints the iterates one by one.

$ tut14b_toi
pebble speed bound mu =  1.2444 (translation 1.2, spin  0.0444), tolerance  1.303E-06
step 0: s =  0.000000, gap  4.185E-01, next step (d - tol/2) / mu =  3.363E-01
step 1: s =  0.336323, gap  2.208E-02, next step (d - tol/2) / mu =  1.774E-02
step 2: s =  0.354068, gap  1.362E-03, next step (d - tol/2) / mu =  1.094E-03
step 3: s =  0.355161, gap  8.522E-05, next step (d - tol/2) / mu =  6.796E-05
step 4: s =  0.355229, gap  5.909E-06, next step (d - tol/2) / mu =  4.225E-06
step 5: s =  0.355234, gap  9.784E-07 <= tolerance: contact
contact T at s =  0.355234 (5 steps)
  gap there  9.784E-07, dragon facet 1085, pebble facet 4
  on the dragon -0.1462 -0.2800 -0.1936, on the pebble -0.1462 -0.2800 -0.1936
directional: s =  0.355234 (4 steps, 789114 exact pair evaluations), same contact to 1e-6: T

• The first step covers 0.336 of the 0.355 that separate the pebble from the contact: the pebble flies almost straight at the dragon, so its closing speed is almost μ. After that each step leaves 5 to 7 % of the gap, and step 5 ends at a gap of 9.784e-7, below the tolerance.
• Each step stays within the gap it starts from, so the pebble cannot tunnel through, whatever its speed. The price is that the step shrinks with the gap: a body moving along a wall, not towards it, takes many short steps (Sliding contact below).

The time of first contact

call dragon%time_of_impact(pebble, pose_identity(), pose_identity(), p0, p1, toi, status, distance=d, point_self=pa, &
                           point_other=pb, facet_self=fa, facet_other=fb, iterations=steps)
print '(A,L1,A,F9.6,A,I0,A)', 'contact ', status == CCD_CONTACT, ' at s = ', toi, ' (', steps, ' steps)'
print '(A,ES10.3,A,I0,A,I0)', '  gap there ', d, ', dragon facet ', fa, ', pebble facet ', fb
print '(A,3F8.4,A,3F8.4)', '  on the dragon', pa, ', on the pebble', pb
call dragon%time_of_impact(pebble, pose_identity(), pose_identity(), p0, p1, toi_dir, status, iterations=steps, &
                           directional=.true., pairs=pairs)
print '(A,F9.6,A,I0,A,I0,A,L1)', 'directional: s = ', toi_dir, ' (', steps, ' steps, ', pairs, &
                                 ' exact pair evaluations), same contact to 1e-6: ', abs(toi_dir - toi) < 1.e-6_R8P

• time_of_impact(other, pose_self0, pose_self1, pose_other0, pose_other1, toi, status, ...) runs the loop above to its end. The dragon is at rest, so both of its poses are pose_identity(). The step is (d − tol/2)/(μ_A + μ_B), and the loop stops at a gap d(toi) ≤ tol.
• status is CCD_CONTACT (with d(toi) ≤ tol), CCD_SEPARATE (toi = 1), CCD_INITIAL_CONTACT (toi = 0) or CCD_NOT_CONVERGED (the cap of max_iterations, default 1000, was hit, and toi is a certified lower bound).
• distance, point_self, point_other, facet_self, facet_other are the gap, the witness points (world frame) and the closest facets at toi; iterations counts the steps. The run prints contact at s = 0.355234 in 5 steps, on dragon facet 1085 and pebble facet 4, the two witness points equal to 4 decimals.
• directional=.true. adds the directional bound of C2A: per pair of triangles, the closing speed projected on the pair's separating direction. The bound holds for each pair of triangles, which is convex, and its minimum over the pairs holds for the whole body. The loop takes the larger of the two certified steps. On this head-on throw it saves one step of five and costs 789114 exact pair evaluations (pairs).
• Without directional the loop is adaptive (issue #75): plain steps, and the directional bound only after eight steps in a row that stall (each advancing s by less than 1 % of what is left and closing less than 1 % of the margin d − tol/2). This throw never stalls, so the default is the plain run above, bit for bit. directional=.false. forces the plain bound only.

Several bodies

Now the dragon is at rest and three pebbles fly over the same step: pebble 2 on the throw above, pebble 3 falling along z through the dragon, pebble 4 flying along x across pebble 3's path, above the dragon.

bodies(1) = dragon ; poses0(1) = pose_identity() ; poses1(1) = poses0(1)       ! the dragon at rest
bodies(2:4) = pebble                                                             ! three pebbles, three throws
poses0(2) = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=0._R8P, translation=[-0.6_R8P, -0.3_R8P, -0.2_R8P])
poses1(2) = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=PI / 2, translation=[0.6_R8P, -0.3_R8P, -0.2_R8P])
poses0(3) = pose_from_axis_angle(axis=[1._R8P, 0._R8P, 0._R8P], angle=0._R8P, translation=[0._R8P, 0.3_R8P, 0.6_R8P])
poses1(3) = pose_from_axis_angle(axis=[1._R8P, 0._R8P, 0._R8P], angle=PI, translation=[0._R8P, 0.3_R8P, -0.6_R8P])
poses0(4) = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=0._R8P, translation=[-0.15_R8P, 0.3_R8P, 0.45_R8P])
poses1(4) = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=PI, translation=[1.05_R8P, 0.3_R8P, 0.45_R8P])
call first_impact(bodies, poses0, poses1, toi, pair, status, tolerance=1.e-6_R8P, tested=tested)
print '(A,L1,A,F9.6,A,2I2,A,I0,A)', 'first contact ', status == CCD_CONTACT, ' at s = ', toi, ', bodies', pair, &
                                    ' (', tested, ' of 6 pairs tested)'

• first_impact(surfaces, poses0, poses1, toi, pair, status, tolerance=, directional=, tested=) gives the earliest contact among N bodies over the step and the pair [i, j] (i < j) that makes it. Body i moves from poses0(i) to poses1(i); a body at rest has equal poses.
• Broad phase: the box of each body at s = 0, grown by its μ and by tol/2, contains the body over the whole step, so two bodies with disjoint grown boxes never come within the tolerance. The other pairs are ranked by a lower bound on their time of impact and run time_of_impact in that order. The search stops when the next bound is later than the best contact found, and tested counts the pairs that ran.

toi_all = 1._R8P                                                                 ! the same by every pair
do i = 1, 3
   do j = i + 1, 4
      call bodies(i)%time_of_impact(bodies(j), poses0(i), poses1(i), poses0(j), poses1(j), s, status, &
                                     tolerance=1.e-6_R8P)
      if (status == CCD_CONTACT) then
         toi_all = min(toi_all, s)
         print '(A,I0,A,I0,A,F9.6)', '  pair ', i, '-', j, ': contact at s = ', s
      else
         print '(A,I0,A,I0,A)', '  pair ', i, '-', j, ': separate'
      endif
   enddo
enddo
print '(A,L1)', 'equal to the loop over all pairs, bit for bit: ', toi_all == toi

• Ties go to the smallest pair, so the result equals the loop over all pairs, bit for bit, when both use the same tolerance. The default tolerance of first_impact is 1e-6 of the largest diagonal of all the bodies, and that of time_of_impact 1e-6 of the larger of the two: here both get tolerance=1.e-6_R8P.

do i = 1, 4                                                                      ! the bodies at the first contact
   at_toi(i) = bodies(i)
   call at_toi(i)%move(pose_interpolate(poses0(i), poses1(i), toi))
enddo
call proximity_pairs_of(at_toi, cutoff=0.3_R8P, pairs=near, distances=gaps)
do i = 1, size(gaps)
   print '(A,I0,A,I0,A,ES10.3)', 'within 0.3 at s = toi: bodies ', near(1, i), ' and ', near(2, i), ', gap ', gaps(i)
enddo

• proximity_pairs_of(surfaces, cutoff, pairs, distances) is the instantaneous screen: the pairs of surfaces closer than cutoff where they stand, with their gaps, found by a box broad phase and distance_to. Here it runs on the bodies moved to the first contact. proximity_pairs(cutoff, pairs, distances) does the same for the bodies (facet_body) of one surface. See body–body distance.

$ tut14c_several
first contact T at s =  0.087759, bodies 3 4 (2 of 6 pairs tested)
  pair 1-2: contact at s =  0.355234
  pair 1-3: contact at s =  0.406194
  pair 1-4: separate
  pair 2-3: separate
  pair 2-4: separate
  pair 3-4: contact at s =  0.087759
equal to the loop over all pairs, bit for bit: T
within 0.3 at s = toi: bodies 1 and 4, gap  2.580E-01
within 0.3 at s = toi: bodies 3 and 4, gap  6.599E-07

• Pebbles 3 and 4 touch first, at s = 0.087759, before either reaches the dragon. first_impact ran 2 of the 6 pairs. The loop over all of them finds three contacts (1-2 at 0.355234, the throw of run 2; 1-3 at 0.406194; 3-4) and the same earliest toi, bit for bit.
• At that instant two pairs are within 0.3: pebbles 3 and 4 at 6.6e-7 (the contact, below the tolerance) and the dragon and pebble 4 at 0.258.

The 2-D twins

use fossil, only : surface_stl_object, curve_object, rigid_pose_2d_t, pose_2d, first_impact_2d, slice_surface, &
                   CCD_CONTACT, CCD_SEPARATE, CCD_NOT_CONVERGED
use penf,   only : I4P, R8P

implicit none
real(R8P), parameter     :: PI = acos(-1._R8P)
type(surface_stl_object) :: naca
type(curve_object)       :: wing, wall, disc, curves(4)
type(rigid_pose_2d_t)    :: poses0(4), poses1(4), q0, q1
real(R8P)                :: d, toi, qa(2), gap
integer(I4P)             :: k, steps, status, pair(2), tested, sa, sb

call naca%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call slice_surface(naca, 3_I4P, 0.5_R8P, wing)                                  ! the section, chord 1, in (x, y)
call make_disc(0.03_R8P, disc)                                                   ! a particle of radius 0.03
curves(1) = wing ; curves(2:4) = disc
poses0(1) = pose_2d(angle=0._R8P) ; poses1(1) = poses0(1)                        ! the wing at rest
poses0(2) = pose_2d(angle=0._R8P, translation=[-0.2_R8P, 0.25_R8P])              ! falls onto the upper side
poses1(2) = pose_2d(angle=PI, translation=[0.4_R8P, -0.25_R8P])
poses0(3) = pose_2d(angle=0._R8P, translation=[-0.4_R8P, -0.12_R8P])             ! passes below
poses1(3) = pose_2d(angle=PI, translation=[0.9_R8P, -0.12_R8P])
poses0(4) = pose_2d(angle=0._R8P, translation=[-0.45_R8P, 0.2_R8P])              ! passes ahead
poses1(4) = pose_2d(angle=PI, translation=[-0.45_R8P, -0.2_R8P])
call wing%distance_to(disc, distance=d, point_self=qa, segment_self=sa, segment_other=sb, pose_other=poses0(2))
print '(A,F8.5,A,I0,A,2F8.4)', 'gap at s = 0 ', d, ', wing segment ', sa, ' at', qa
call wing%time_of_impact(disc, poses0(1), poses1(1), poses0(2), poses1(2), toi, status, iterations=steps)
print '(A,L1,A,F9.6,A,I0,A)', 'particle 2: contact ', status == CCD_CONTACT, ' at s = ', toi, ' (', steps, ' steps)'
call first_impact_2d(curves, poses0, poses1, toi, pair, status, tested=tested)
print '(A,L1,A,F9.6,A,2I2,A,I0,A)', 'first contact ', status == CCD_CONTACT, ' at s = ', toi, ', curves', pair, &
                                    ' (', tested, ' of 6 pairs tested)'

• curve%distance_to(other, distance=, point_self=, point_other=, segment_self=, segment_other=, intersecting=, pose_self=, pose_other=) is the 2-D gap, with rigid_pose_2d_t poses from pose_2d(angle=, translation=). Curves are small, so every segment pair is evaluated behind a box test, and contact is decided by the exact orient2d.
• curve%time_of_impact and first_impact_2d(curves, poses0, poses1, ...) take the arguments of the 3-D calls, with segment_self, segment_other for the facets (and motion_bound_2d for the bound). A body is a curve_object: here the section of naca0012-ascii.stl at mid-span (slice_surface) and three discs built with adopt_loops.

$ tut14d_curves
gap at s = 0  0.17935, wing segment 47 at -0.1608  0.0444
particle 2: contact T at s =  0.322097 (14 steps)
first contact T at s =  0.322097, curves 1 2 (5 of 6 pairs tested)
gap  2.0E-02, plain:       CCD_SEPARATE      at s =  1.000000 (41 steps)
gap  2.0E-02, default:     CCD_SEPARATE      at s =  1.000000 (41 steps)
gap  2.0E-02, directional: CCD_SEPARATE      at s =  1.000000 (1 steps)
gap  4.0E-06, plain:       CCD_NOT_CONVERGED at s =  0.004375 (1000 steps)
gap  4.0E-06, default:     CCD_SEPARATE      at s =  1.000000 (9 steps)
gap  4.0E-06, directional: CCD_SEPARATE      at s =  1.000000 (1 steps)

• At s = 0 particle 2 is 0.179 above the section, closest to segment 47. It falls onto the upper side and touches it at s = 0.322097 after 14 steps. That is the first contact of the four curves; first_impact_2d tested 5 of the 6 pairs.

Sliding contact

call wall%adopt_loops(u=[0._R8P, 1._R8P, 1._R8P, 0._R8P], v=[-0.1_R8P, -0.1_R8P, 0._R8P, 0._R8P], loop_sizes=[4_I4P])
do k = 1, 2                                                                      ! sliding 0.8 along the wall
   gap = merge(0.02_R8P, 4.e-6_R8P, k == 1)                                      ! at a gap of 0.02, then 4 tolerances
   q0 = pose_2d(angle=0._R8P, translation=[0.1_R8P, 0.03_R8P + gap])
   q1 = pose_2d(angle=0._R8P, translation=[0.9_R8P, 0.03_R8P + gap])
   call wall%time_of_impact(disc, pose_2d(angle=0._R8P), pose_2d(angle=0._R8P), q0, q1, toi, status, &
                            tolerance=1.e-6_R8P, iterations=steps, directional=.false.)
   print '(A,ES8.1,A,A,A,F9.6,A,I0,A)', 'gap ', gap, ', plain:       ', status_name(status), ' at s = ', toi, &
                                      ' (', steps, ' steps)'
   call wall%time_of_impact(disc, pose_2d(angle=0._R8P), pose_2d(angle=0._R8P), q0, q1, toi, status, &
                            tolerance=1.e-6_R8P, iterations=steps)
   print '(A,ES8.1,A,A,A,F9.6,A,I0,A)', 'gap ', gap, ', default:     ', status_name(status), ' at s = ', toi, &
                                      ' (', steps, ' steps)'
   call wall%time_of_impact(disc, pose_2d(angle=0._R8P), pose_2d(angle=0._R8P), q0, q1, toi, status, &
                            tolerance=1.e-6_R8P, iterations=steps, directional=.true.)
   print '(A,ES8.1,A,A,A,F9.6,A,I0,A)', 'gap ', gap, ', directional: ', status_name(status), ' at s = ', toi, &
                                      ' (', steps, ' steps)'
enddo

• The plain bound cannot tell motion along a wall from motion towards it: a particle moving at speed μ is assumed to close the gap at μ. At a gap of 0.02 it advances by about the gap per step and needs 41 steps to cross a step in which it never comes closer (last six lines of the run above, directional=.false. for the plain bound alone).
• At a gap of 4e-6, four tolerances, each step is about 3.5e-6 long: the plain loop hits the cap of 1000 steps at s = 0.004375 and returns CCD_NOT_CONVERGED with that certified lower bound.
• The directional bound sees that the particle does not approach the wall and crosses the whole step in one, at either gap (directional=.true.).
• The default arguments switch between the two (issue #75). At a gap of 0.02 each plain step covers a few per cent of what is left, so the loop never stalls: the same 41 steps as the plain bound. At 4e-6 the plain steps stall, and after eight of them the directional step crosses the rest: CCD_SEPARATE in 9 steps. Resting contact, particles in a channel, a valve against its seat converge with the default arguments.

What you learned

A gap sampled in time can miss a contact. motion_bound_3d bounds the speed of a body over a step, and time_of_impact advances by the steps that bound certifies, to the first contact, and switches to the directional bound by itself when bodies slide along each other. first_impact and first_impact_2d give the earliest contact among several bodies, equal to the loop over all pairs, and proximity_pairs_of the pairs close at one instant. Reference: Continuous collision detection, Body–body distance.

Next: 15. Body-fitted grids.