David Baraff
email: baraff@cs.cmu.edu

Curved Surfaces and Coherence for Non-penetrating Rigid Body Simulation

Keywords: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling; I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism dynamics, constraints, simulation

A formulation for the contact forces between curved surfaces in resting (non-colliding) contact is presented. In contrast to previous formulations, constraints on the allowable tangential movement between contacting surfaces are not required. Surfaces are restricted to be twice-differentiable surfaces without boundary. Only finitely many contact points between surfaces are allowed; however, the surfaces need not be convex. The formulation yields the contact forces between curved surfaces and polyhedra as well. Algorithms for performing collision detection during simulation on bodies composed of both polyhedra and strictly convex curved surfaces are also presented. The collision detection algorithms exploit the geometric coherence between successive time steps of the simulation to achieve efficient running times.

Introduction

One of the most difficult behaviors to simulate in rigid body dynamics is the non-penetration constraint between solid bodies. The two problems involved in simulating non-penetrating rigid bodies are [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF] detecting collisions and contact between pairs of bodies and (2) determining the contact forces present between contacting bodies.

The force determination problem can be solved by analytical and non-analytical methods. Both Barzel and Barr [START_REF] Barzel | A modeling system based on dynamic constraints[END_REF] and Baraff [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF] give motivations for preferring analytical methods over non-analytical methods in rigid body simulation. Analytical formulations for the contact forces that arise between polyhedral bodies have been presented in [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF][START_REF] Featherstone | Robot Dynamics Algorithms[END_REF][START_REF] Lo ¨tstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF][START_REF] Lo ¨tstedt | Numerical simulation of time-dependent contact friction problems in rigid body mechanics[END_REF]. These formulations are the most general possible in that they express the contact forces between bodies that are completely unconstrained in their tangential (sliding) movement. For curved surfaces, formulations for the contact forces that arise to prevent inter-penetration have only been realized for certain cases of constrained tangential movement. For example, if two curved surfaces are restricted to roll without slipping, the contact force between them is easily determined [START_REF] Goyal | Second order kinematic constraint between two bodies rolling, twisting and slipping against each other while maintaining point contact[END_REF][START_REF] Neȋmark | Dynamics of Nonholonomic Systems[END_REF]. Similar restrictions such as rolling with a specified slip velocity also have simple analytical solutions. The general case of non-penetration between curved surfaces without any constraint on the tangential movement poses a much more difficult problem; the extension from polyhedra to curved surfaces is not straightforward. We present a formulation of the contact forces between curved surfaces that are completely unconstrained in their tangential movement. We have not encountered a formulation for this problem in any previous literature. We restrict ourselves to the case of twice-differentiable curved surfaces without boundary that contact at only finitely many points. Configurations that result in one-or two-dimensional contact regions are not dealt with. The surfaces need not be convex in the neighborhood of a contact point and may be defined by either implicit or parametric equations. A formulation of the contact forces between curved surfaces and polyhedra is derived as a special case of contact between two curved surfaces.

The collision detection problem has an extensive literary background in the fields of computational geometry and robotics. Computational geometry focuses on problems posed in terms of a static environment. Robotics generally focuses on problems posed in terms of a dynamic environment; the movement of bodies is known in terms of some function of time. In both cases, the emphasis is on developing the best algorithm for solving a single problem, either static or dynamic. In contrast, dynamic simulation involves the solution of a sequence of static problems, one per time step. Although each problem of the sequence can be solved separately using previous collision detection methods, algorithms specifically designed to solve a sequence of related problems are more efficient. We present efficient collision detection algorithms for polyhedra and convex closed curved surfaces by exploiting the geometric coherence between successive collision detection problems of the simulation. Surfaces can be defined either implicitly or parametrically.

Overview

Simulation of non-penetrating rigid bodies by analytical methods involves the basic flow of control shown in figure 1. At every time step, bodies are examined pairwise for possible interpenetration. If two bodies are found to inter-penetrate, the simulator backtracks to the point in time immediately before the interpenetration occurred. Once a configuration without interpenetration is achieved, the contact points between all the bodies are found. Finally, a system of constraint equations based on the contact points yields the analytically correct contact forces and impulses at every contact point. After the contact forces and impulses are applied to the bodies, a new time step is begun. We group this series of steps into two phases: contact determination (steps A-D of figure 1) and force determination (step E of figure 1) For curved surfaces, both the inter-penetration check and the backtracking steps of the collision determination phase make use of the same mathematical derivations used in the force determination phase. Accordingly, we will deal first with extending the force determination model to curved surfaces.

Analytical Force Determination

Featherstone [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF] gives a complete derivation of the system of constraint equations used to find contact forces between rigid bodies. We present a brief review of the mathematical structure of the problem.

At some time t 0 we are given a collection of nonpenetrating perfectly rigid bodies that contact at some number of points and asked to calculate the forces between the objects that would naturally arise to prevent inter-penetration. Contact points at which bodies are colliding give rise to contact impulses. Methods for calculating the contact impulses between bodies of any geometry are given in [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF][START_REF] Featherstone | Robot Dynamics Algorithms[END_REF][START_REF] Goldstein | Classical Mechanics[END_REF]. Contact impulses are calculated and applied prior to considering contact forces. Accordingly it is assumed that the configuration of bodies being analyzed has no colliding contacts. Consider a contact point p c between two bodies A and B at time t 0 (figure 2). Let the unit surface normal at the contact point be n ˆ(see Baraff[1] for the case when n ˆis not well defined). In the absence of friction, the (as yet) unknown contact force F is written F = f n ˆwith f the unknown (scalar) contact force magnitude at time t 0 .

n F = f A B n p c Figure

Contact between two objects.

The primary consideration for calculating the unknown contact force magnitude f at each contact point lies in the nonpenetration constraint between bodies. At every contact point p c between two bodies A and B at time t 0 , a geometrical constraint that prevents A and B from inter-penetrating near p c is con-structed. Each constraint is converted by a differentiation operation into a constraint on the contact force magnitude f at p c . The geometrical constraint is expressed by a constraint function (or characteristic function ) χ(t). A characteristic function χ(t) is a function of time that characterizes the geometric relation of A and B close to p c at times near t 0 .1 χ takes on values as follows: Given χ, we can express the constraint that A and B not interpenetrate near p c as

χ(t) =
χ(t) ≥ 0. ( 2
)
For example, consider a potential vertex-face contact between two planar polygons (figure 3).

n pb B χ(t) > 0 n pb B χ(t) = 0 n pb B χ(t) < 0 p a A p a A p a A Figure 3. Vertex-face contact.
Let p a (t) be the position of the vertex of A, p b (t) be the position of any fixed point on the contact face of B and let n ˆ(t) be the outwards unit normal of the contact face. The characteristic function for this situation is

χ(t) = n ˆ(t) . (p a (t) -p b (t)).
(

) 3
From figure 3, χ(t) is positive, zero, or negative according to whether p a (t) lies above, on, or below the contact face of B at time t.

Intuitively, χ may be regarded as a measure of the distance between A and B, with negative distance indicating interpenetration. Similarly, χ . may be regarded as a measure of the relative velocity between A and B. If χ . (t 0 ) < 0 then χ is decreas- ing and the bodies are colliding; however, it was assumed that the configuration had no colliding contact points. Likewise, if χ . (t 0 ) > 0 then χ is increasing and the bodies are separating; in this case the contact force is automatically zero and the contact point may be disregarded. Thus, the only contact points considered are those for which χ

. (t 0 ) = 0.
In order to convert the geometric constraint χ(t) ≥ 0 (equation (2)) into a constraint on the contact forces, we take the second derivative of χ with respect to time and require that χ

.. (t 0 ) ≥ 0.

(4)

Informally, we have constrained the relative "acceleration" χ .. , between A and B, to be non-negative. (Strictly speaking however, χ .. is not a physical measure of acceleration). Analytically, since

χ(t 0 ) = χ
. (t 0 ) = 0, χ .. (t 0 ) < 0 would make χ a decreasing function at t 0 . This would violate the constraint of equation [START_REF] Baraff | Determining frictional inconsistency for rigid bodies is NP-complete[END_REF].

How does χ .. (t 0 ) ≥ 0 constrain the contact forces? While χ(t 0 ) and χ . (t 0 ) are independent of any internal or external forces, χ .. (t 0 ) is a linear expression of the contact forces at time t 0 . Intui- tively, the contact forces must be "strong enough" to satisfy equation (4), and thus prevent A and B from accelerating towards each other at p c . Appendices A and C show that equation ( 4) is a linear inequality constraint on the contact force magnitudes.

In addition to the geometrically motivated constraint χ .. (t 0 ) ≥ 0, there is an additional relationship between χ .. (t 0 ) and f that must be satisfied. If χ .. (t 0 ) > 0, then χ is an increasing func- tion at time t 0 and A and B are separating at p c . In this case the contact force is zero. However, if χ .. (t 0 ) = 0 then A and B are not separating and f need not be zero. The relationship between f and χ .. (t 0 ) is known as a complementarity condition; it is written as

f χ .. (t 0 ) = 0 (5)
to express the fact that either f or χ .. is zero. If we impose the res- triction that f be non-negative, so that objects can "push" but not "pull" on each other, then a configuration with N contact points must satisfy the system of equations χ ..

i (t 0 ) ≥ 0, f i χ .. i (t 0 ) = 0, f i ≥ 0 (1≤i≤N) (6)
where f i and χ i are the contact force and constraint function for the ith contact point. In a previous paper, we proposed a heuristic method for solving equation [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF]. For the case of frictionless contact, equation ( 6) forms what is known as a positive semidefinite (PSD) linear complementarity problem. Equation ( 6) can also be viewed as a PSD quadratic programming problem. Efficient numerical algorithms exist that solve PSD linear complementarity problems 2 and PSD quadratic programs [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF], and we advocate their use over the heuristic solution method. However, in the presence of friction, it is known that equation ( 6) is no longer necessarily PSD. Finding the solution of a non-PSD linear complementarity problem or quadratic program is NP-hard [START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF]. A recent result [START_REF] Baraff | Determining frictional inconsistency for rigid bodies is NP-complete[END_REF] shows that finding a solution to equation [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF] in the presence of friction is also NP-hard. Thus, heuristic solution methods may indeed be necessary for practical simulations. For both PSD and non-PSD systems, coherence based methods can be exploited [START_REF] Lo ¨tstedt | Numerical simulation of time-dependent contact friction problems in rigid body mechanics[END_REF][START_REF] Murty | Linear Complementarity, Linear and Nonlinear Programming[END_REF] to reduce the computational expense of solving equation [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF]. See Lo ¨tstedt [START_REF] Lo ¨tstedt | Mechanical systems of rigid bodies subject to unilateral constraints[END_REF][START_REF] Lo ¨tstedt | Numerical simulation of time-dependent contact friction problems in rigid body mechanics[END_REF] and Featherstone [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF] for further discussions on the properties of this constraint system and methods for solving it.

Analytical Forces between Curved Surfaces

In the case of polyhedral objects, the characteristic function χ and its second derivative χ .. are readily available. Additionally, contact between polyhedra may result in line segments or polygons of contact. Although this results in an infinity of contact points, constraint functions need be formulated only for the finitely many vertices of the convex hull polygon of the contact line or area [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF]. For curved surfaces however, the convex hull of the contact area may not be a polygon. For example, a cylinder standing upright on a plane has a circular area of contact points. The convex hull of this contact region is a circle and cannot be described by a finite number of vertices. We have not developed a constraint for contact regions of dimension one or higher. Although discretization of the boundary of the contact area is one possibility, we would rather deal with an analytical formulation over the entire boundary. For curved surfaces, we will restrict our attention to situations in which the number of contact points is finite. We will construct a characteristic function χ for each contact point.

The difficulty in formulating a geometric constraint function for curved surfaces in contact is the need to construct a formula specific enough to be differentiable. How can we formulate the geometric constraint that all points on surface A near a point p c remain on or outside surface B? Furthermore, how can we write this as a scalar-valued differentiable function that is positive, zero, or negative according to whether A and B are disjoint, contacting, or inter-penetrating? One possible start is to let χ be the minimum distance between A and B near p c , and require that χ always be non-negative. The minimum distance is positive when A and B are separate and zero when they are contacting. However, as figure 4 shows, the minimum distance is not negative when A and B initially inter-penetrate; it is zero. 3 Additionally, the minimum distance is not differentiable at the time that A and B first contact at p c .

Closely related to the minimum distance however is the concept of an extreme distance. The extreme distance between A and B near p c is defined as follows. If A and B are disjoint near p c then the extreme distance between A and B is just the normal minimum distance between A and B (near p c ). If A and B are in contact at p c , then the extreme distance is zero. If A and B have inter-penetrated near p c , then the extreme distance is the maximum distance between A and B (near p c ). The extremal points of A and B are the two points p a and p b on A and B that realize the extremal distance (figure 4). Given the above definitions, we can construct a constraint function by letting χ(t) be a positive multiple of the extreme distance when A and B are disjoint or contacting. When A and B have inter-penetrated, we will let χ(t) be a negative multiple of the extreme distance. In the next section we will show how p a and p b can be used to construct such a formula. The restriction to situations where only finitely many contact points arise guarantees that the extreme distance (and the extremal points) of A and B sufficiently near p c will be unique.

Deriving χ ..

We still must construct an explicit formula for χ so that it may be differentiated to find χ .. . Although the χ we develop in this section yields an impractical result (computationally speaking), we feel its presentation is necessary to clearly understand the final form of χ and χ .. derived in section 6. The derivation of χ and χ .. will assume implicit definitions of the curved surfaces; however, the end result depends solely on the derivatives of the surfaces at the contact points. As a result, parametric definitions of the surfaces can be used as easily as implicit definitions for determining contact forces; see appendix D for details.

We model the two curved surfaces of A and B as implicit time-varying functions F (p,t) and G (p,t) where p is a point in world space. At time t, a point p is on the surface of A iff F(p,t) = 0. (We will refer to F and G as both functions and surfaces). Furthermore, if F (p,t) < 0 then p is inside A, and if F (p,t) > 0 then p is outside A. The function G similarly defines the shape of B. We will use the notation

∂p ∂ F (p, t ) = F ′ (p, t ) = ∇F(p,t) T , (7)
where ∇F(p,t) is a column vector (and F ′(p,t) a row vector). If F (p,t) = 0, then ∇F(p,t) is the outwards directed surface normal of F at point p at time t (figure 5). Given these definitions, we can express χ as

χ(t) = ∇G(p b ,t) . (p a (t) -p b (t)) (8)
where p a and p b are the two extremal points between A and B at time t. From figure 6, we see that ∇G(p b ,t) is colinear with the vector p ap b . The derivative of equation ( 8) requires the derivatives of p a and p b . Although the curved surface characteristic function appears similar to the polyhedral characteristic function (equation (3)), the latter is easily differentiated while the former is not. In the polyhedral characteristic function, the points p a and p b denote fixed points of A and B. The derivatives of a fixed point of a rigid body are simply expressed in terms of the motion of the body [START_REF] Barzel | A modeling system based on dynamic constraints[END_REF][START_REF] Goldstein | Classical Mechanics[END_REF]. However, in the case of curved surfaces, p a and p b are not fixed points of A and B. p a and p b change positions in two ways. First, p a and p b move according to the rigid body motion of A and B. Second, p a and p b change positions in the body space of A and B (figure 7). In order to differentiate p a and p b , they must be redefined in such a way that they can be differentiated. We define the extremal points p a and p b at time t as the (unique) pair of points near p c that satisfy the conditions

χ(t) > 0 χ(t) = 0 χ(t) < 0
E 4 : (p b -p a ) + λ 1 ∇G(p b ,t) = 0 → . E 3 : G(p b ,t) = 0 E 2 : F(p a ,t) = 0 E 1 : ∇F(p a ,t) + λ 2 ∇G(p b ,t) = 0 → (9)
λ 1 and λ 2 are unconstrained scalar values. Condition E 1 guarantees that the surface normals at the extremal points are colinear. Conditions E 2 and E 3 guarantee that p a and p b are points on A and B, and condition E 4 guarantees that p b 's displacement from p a is colinear to the surface normals. A formal justification of conditions E 1 thru E 4 as a definition of the extremal points may be found in any advanced calculus text; see for example Taylor and Mann [START_REF] Taylor | Advanced Calculus[END_REF]. Figure 6 shows the geometric intuition behind equation [START_REF] Goldstein | Classical Mechanics[END_REF]. Equation ( 9) is closely related to the Lagrange multiplier formulation for constrained minimization (hence our choice of λ as a symbol for the multipliers of ∇G).

For the case of contact between a polyhedron A and a curved surface B, equation ( 9) is modified depending on whether p c is coincident with a face, edge or vertex of the polyhedron. If p c lies in a face, then F is the implicit function of the plane embedding the face. Otherwise, if p c lies on an edge and ∇G(p b ,t 0 ) is not perpendicular to any adjoining face, p a is constrained to be the extremal point on the edge. If p c lies on a vertex, and ∇G(p b ,t 0 ) is not perpendicular to any adjoining edge (and hence any adjoining face), p a is defined to be coincident with the vertex. In the above two cases, a system of equations similar to ( 9) is formed and used in place of equation [START_REF] Goldstein | Classical Mechanics[END_REF].

Conditions E 1 thru E 4 are used to derive expressions for the derivatives of p a and p b . (We will consider numerical solutions of equation ( 9) to find p a and p b when we deal with collision determination in section 7.) Given extremal points p a and p b , we can find p . a and p .

b by making use of the implicit function theorem for simultaneous equations from calculus [START_REF] Taylor | Advanced Calculus[END_REF]. This theorem asserts that under proper conditions, p a and p b may be regarded as functions of time; the theorem also gives an analytic expression for the derivatives of p a and p b (with respect to time). We will write the Jacobian determinant of a set of vector functions

H 1 (x → ) thru H n (x → ) as ∂(x 1 , . . . ,x n ) ∂(H 1 , . . . ,H n ) = ∂x 1 ∂H n . . . ∂x 1 ∂H 1 . . . . . . . . . ∂x n ∂H n . . . ∂x n ∂H 1 . ( 10
)
If we are dealing with 3-space, then p a and p b are each three (scalar) functions of time: p a x (t), p a y (t) and p a z (t) and similarly for p b . The implicit function theorem, applied to equation (9) yields

p . a x (t) = - J K ( 11
)
where

K = ∂(t,p a y ,p a z ,p b x ,p b y ,p b z ,λ 1 ,λ 2 ) ∂(E 1 ,E 2 ,E 3 ,E 4 ) (12)
and .. b ; these in turn require derivatives of K and J. Unfortunately, a symbolic expression for the determinant of K or J is impractical. Although K and J have considerable block structure, block structure cannot be exploited in computing determinants. In its present form, the Jacobian determinant contains more than 1,000 terms of seven factors each; the derivative would contain far more terms. A transformation of coordinate systems and functions is presented in the next section that yields more tractable expressions.

J = ∂(p a x ,p a y ,p a z ,p b x ,p b y ,p b z ,λ 1 ,λ 2 ) ∂(E 1 ,E 2 ,E 3 ,E 4 ) . ( 13

Coordinate Transformations

The formulation for p . a and p . b in the last section involved derivatives of determinants of 8 × 8 matrices. By choosing an appropriate coordinate system and transforming the representation of the surface functions F and G, we can find a tractable representation for p . a and p . b . First, we assume a rotated coordinate system in which ∇F(p a ,t 0 ) and ∇G(p b ,t 0 ) are colinear with the z axis at time t 0 , with ∇G(p b ,t 0 ) directed positively along z. (We will employ the standard right-handed coordinate system used to depict functions z=h(x,y), with z the vertical axis). The effect of this rotation on derivatives of F and G is discussed in appendix B. Next, we explicitly model A near p a as a time-varying scalar function f of x and y. Where F was a function F (x,y,z,t), f is a function f (x,y,t). The justification for the existence of f is the implicit function theorem of calculus. The formal definition of f is

F(x,y, f(x,y,t),t) = 0. ( 14
)
Thus, the point (x,y, f(x,y,t)) is a point on F at time t. A function g is chosen similarly for G. In this new coordinate system, the extremal points p a and p b share the same x and y coordinates at time t 0 . At times t near t 0 , A and B do not penetrate as long as the z value of p a is greater or equal to the z value of p b (figure 8). This can be expressed as

χ(t) = f(p a x ,p a y ,t) -g(p b x ,p b y ,t). (15)
The condition χ .. ≥ 0 is simply

dt 2 d 2 f(p a x ,p a y ,t 0 ) - dt 2 d 2 g(p b x ,p b y ,t 0 ) ≥ 0. ( 16
)
The extremal points are then expressed in terms of f and g, without the use of the multipliers λ 1 and λ 2 (see appendix C,

z = f (x,y) z = g (x,

Contact Determination

Given the vast literature on collision detection in the computational geometry and robotics fields, it is with some trepidation that we present new algorithms for collision detection and determination of contact points. The algorithms presented are very basic; nonetheless we have found them to be extremely efficient for our simulations. In a dynamic simulation, a series of collision detection problems is encountered. Each problem is similar to the collision detection problem posed and solved during the previous time step. The focus of the algorithms presented in this section is using information from the previous time step to solve the collision detection problem for the current time step. The collision detection problems addressed in computational geometry and robotics are of a different nature.

In the field of computational geometry, collision detection algorithms are by and large restricted to static geometrical configurations. Algorithms are typically developed to solve a single instance of a problem involving fixed objects in the smallest asymptotic time complexity. The use of algorithms of this nature for collision detection during simulation essentially ignores any information discovered in previous time steps. Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three space[END_REF] present an algorithm that efficiently computes the minimum distance between convex objects. Additionally, the algorithm can use information from previous time steps for fast initialization, and would appear to be an attractive candidate for detecting collisions. However, we are not interested in the value of the minimum distance per se. and we feel that the algorithms presented below are better suited to our simulation environment. For the simpler problem of determining the disjointness of two convex polyhedra with a total of n vertices, algorithms with asymptotic time complexities of O(nlogn) and even O (n) have been achieved; however, it is not clear that these algorithms are practically useful for reasonable values of n [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three space[END_REF].

In the field of robotics, algorithms are developed for geometric problems involving objects with specified continuous motions over some time period; for example, determining the first collision between polyhedral objects [START_REF] Canny | Collision detection for moving polyhedra[END_REF]. These algorithms presuppose known paths for the objects and detect collisions over the length of the paths. In our case, the paths of the objects are not known.

We have found that the geometric relationship between objects in our simulations does not change very much between successive time steps. As a result, collision detection algorithms for our simulation environment should be structured to take maximum advantage of the geometric coherence between successive time steps. Our simulator uses contact determination algorithms that assume a high degree of coherence between time steps. If the relative displacements of objects between successive time steps are large, this assumption breaks down. However, our experience has been that numerical considerations in solving the differential equations of motion limit the size of the time step. We choose not to treat the case where an object passes completely through another object in one time step; one solution to this problem is 4D space-time swept volume algorithms to detect collisions [START_REF] Moore | Collision detection and response for computer animation[END_REF].

In this paper we limit our objects to the union of convex polyhedra and strictly convex closed surfaces. We will refer to these polyhedra and curved surfaces as primitives . The first geometric problem we consider is the pairwise comparison of primitives to determine inter-penetration. 4 Our primary mechanism for exploiting coherence will be through the use of witnesses to the decision problem of inter-penetration. A witness is some piece of information that can be used to quickly answer a decision problem. We will utilize coherence by caching witnesses from one time step to the next; hopefully a witness from the previous time step will be a witness during the current time step.

Since primitives are convex, a pair of primitives do not inter-penetrate if and only if a separating plane between them exists. 5 A separating plane between two objects is a plane such that each object lies in a different half-space of the plane. A given plane can be verified to be a separating plane between two convex polyhedra in time O (n) where n is the total number of vertices of the two polyhedra. Cundall [START_REF] Cundall | Formulation of a three-dimensional distinct element model -Part I. A scheme to represent contacts in a system composed of many polyhedral blocks[END_REF] performed pairwise collision detection during simulation by initially finding separating planes that were approximately equidistant from each polyhedron of the pair. At later time steps, numerical techniques were used to quickly update the separating planes so that they maintained their equidistant relationship between pairs of polyhedra. A simpler solution exists however; for disjoint or contacting convex polyhedra, it can be shown that a separating plane exists which either embeds at least one face of one of the polyhedra or embeds an edge from each polyhedron. If a pair of polyhedra are shown not to inter-penetrate by one of these separating planes, we cache the face or two edges embedded in the separating plane as a witness. At the next time step, we use the cached face or edges to form a new separating plane. In this manner the new separating plane is obtained simply from the old one without the need for any numerical computations. Even better, by caching the nearest face or edge to the separating plane from each polyhedron, disjointness can usually be verified in sublinear time (figure 9).

If two polyhedra are inter-penetrating, it is almost always the case that either a vertex of one polyhedron is inside the other, or an edge of one polyhedron has intersected a face of the other. In this case, the inter-penetrating vertex, or intersecting edge and face are cached as a witness to the inter-penetration. In subse- quent comparisons, the witness is used to quickly check for interpenetration. When it is necessary to initially find a witness, a sophisticated computational geometry algorithm such as Gilbert et al. [START_REF] Gilbert | A fast procedure for computing the distance between complex objects in three space[END_REF] might be employed. Currently we use exhaustive search to initially find a witness; we have found the added expense, amortized over the length of the simulation, to be negligible.

The use of separating planes also makes the contact determination step simple. Contact points between a pair of polyhedra separated by a plane P can only occur on the plane P. Given the separating plane, the contact points are quickly and efficiently determined by comparing only those faces, edges or vertices coincident with the separating plane. This determination can itself be performed quickly by caching information from the previous time step.

For comparing two curved surfaces, we must employ numerical techniques to determine disjointness. We make use of the concept of extremal points from section 5 in determining disjointness. Given two disjoint convex curved surfaces, the extremal points are the points of minimum distance and are a witness to the disjointness of the surfaces. Otherwise, the extremal points near the intersection of the surfaces are a witness that the surfaces do inter-penetrate (figure 6). To find these extremal points, a non-linear equation solver may be used to solve equation [START_REF] Goldstein | Classical Mechanics[END_REF] for p a and p b 6 ; see Forsythe et al. [START_REF] Forsythe | Computer Methods for Mathematical Computations[END_REF]. However, equation ( 9) admits multiple solutions of p a and p b ; we are interested in the solution that globally minimizes ||p ap b ||. 7 Non-linear equation solvers proceed from some initial estimate of the solution to an exact solution (within numerical tolerances). If we are initially finding the extremal points between the surfaces, we require some rough estimate of p a and p b , so that the solver will converge to the proper solution. For implicit surfaces, we can initially estimate p a and p b by intersecting the two surfaces with the line connecting the centroids of the two surfaces. For parametric surfaces, we can generate and store a coarse mesh of surface points for each surface and use the parametric coordinates of the minimum distance pair of points as an initial estimate. Using the initial estimate as a starting point, the solver converges to the proper solution of the extremal points.

Once the extremal points p a and p b are determined, they are cached for the next time step. In subsequent time steps, the cached extremal points are used as an initial estimate and the solver converges in a few iterations to the new extremal points. Furthermore, the accuracy of this initial estimate can be significantly improved as follows. If in addition to caching the extremal points p a and p b at time t 0 we cache p . a and p . b , we can estimate p a and p b at time t 0 + ∆t by p a (t 0 + ∆t) = p a (t 0 ) + ∆t p . a (t 0 ) [START_REF] Tomlin | Robust implementation of Lemke's method for the linear complementarity problem[END_REF] and similarly for p b . (See appendix D for the parametric case). Improving the accuracy of the initial estimates of p a and p b adds to the speed and robustness of the algorithm. The derivatives of p . a (t 0 ) and p . b (t 0 ) can be calculated in terms of determinants of 4 × 4 matrices as described in appendix C. Once two curved surfaces have been found not to inter-penetrate, the contact determination process consists of merely comparing the distance between the extremal points to some numerical threshold.

Comparisons between a polyhedral primitive and a curved surface primitive are handled in an analogous manner to contact between a polyhedron and a curved surface (section 5).

Determining the Collision Time

The last problem in the collision detection phase is backtracking to the point of a collision (step E of figure 1). During a simulation, it may happen that two objects A and B come into colliding contact at some time t c . Suppose that t 0 was the time of the previous time step when the objects had not yet collided, and at time t 1 (t 0 < t c < t 1 ) it is found that A and B have interpenetrated. When this occurs, the simulator makes a prediction t p of the time t c at which the initial collision occurred, backs up to time t 0 and then moves forward to time t p . If the simulator finds that A and B have not yet collided at t p , it assumes that t p < t c and makes a larger prediction for t c . Conversely, if the simulator finds that A and B have inter-penetrated at t p , it assumes that t c < t p , and makes a smaller prediction for t c . Otherwise, the simulator has found t c to within numerically accepted tolerances and may proceed.

Conceptually, this can be viewed as a root finding problem. Previous papers have solved this root finding problem by using bisection [START_REF] Moore | Collision detection and response for computer animation[END_REF] or regula falsa [START_REF] Baraff | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF]. The bisection method is extremely robust, simple to implement, and independent of geometry; however, it has relatively slow convergence, especially where great accuracy is required. The regula falsa method linearly interpolates the distance between A and B at t 0 and the amount of inter-penetration at t 1 to predict t c . The regula falsa method handles any geometry as long as measures of separation and inter-penetration are available; additionally, the method converges faster than bisection. regula falsa is not as robust as bisection but a hybrid bisection-regula falsa algorithm [START_REF] Forsythe | Computer Methods for Mathematical Computations[END_REF] works well in practice.

An alternative to regula falsa is Newton's method. Newton's method requires both a measure of separation between A and B at t 0 and the relative approach velocity at time t 0 . These quantities however are exactly modeled by χ and χ . . Newton's method for solving h (t) = 0 near the point t 0 is based on the Taylor series expansion

h (t) = h (t 0 ) + n =1 Σ ∞ n ! h (n) (t 0 ) (t -t 0 ) n . ( 18
)
By throwing out terms for n > 1 and replacing h with χ, we obtain the linear approximation

t = t 0 -χ . (t 0 ) χ(t 0 ) . ( 19
)
While Newton's method gives better convergence than the regula falsa method, both regula falsa and Newton's method consistently either under-estimate or over-estimate t c for constant (non-zero) acceleration. Since constant acceleration occurs frequently, it makes sense to predict t c by using a quadratic model as opposed to the linear model used by Newton's method. The quadratic model requires χ, χ

. and χ .. ; section 5 and appendix C show how χ, χ . and χ .. are calculated for polyhedral and curved objects. Fol- lowing the derivation of Newton's algorithm, we predict t c by

χ(t) = χ(t 0 ) + χ . (t 0 )(t -t 0 ) + 2 χ .. (t 0 ) (t -t 0 ) 2 . ( 20
)
Solving for t we obtain

t = t 0 + χ .. (t 0 ) -χ . (t 0 ) ± √ χ . (t 0 ) 2 -2χ .. (t 0 )χ(t 0 ) . ( 21
)
We set t p to be the smallest real root of equation ( 21) greater than t 0 ; if no such root exists, or χ .. (t 0 ) is zero (within numerical toler- ance), we use Newton's method. The method is made robust by incorporating a bisection step whenever convergence is slow [START_REF] Forsythe | Computer Methods for Mathematical Computations[END_REF].

For constant acceleration, equation ( 21) gives an exact result as long as χ is a linear measure of the distance. 8 For non-constant acceleration, the quadratic model still converges faster than Newton's method, close to t c .

Conclusion

(t) = v (t), R . (t) = ω * (t)R(t) (25)
where ω * (t) is the dual [START_REF] Barzel | A modeling system based on dynamic constraints[END_REF][START_REF] Goldstein | Classical Mechanics[END_REF] of ω(t). 9 Define the point-velocity function V (p,t) as

V(p,t) = v (t) + ω(t) × (p -c (t)). (26)
Then

T ˜. (p,t) = ω * (t)R(t) T (p -c (t)) -R (t) T v (t) = -R (t) T ω(t) × (p-c (t)) + v (t) = -R (t) T V(p,t). ( 27
)
We can represent the movement of F r by the rigid body motion T (p,t) by defining

F(p,t) = F r (T ˜(p,t)). ( 28
)
If q is some point on F r then at time t the point T (q,t) is a point of F since

F(T (q,t),t) = F r (T ˜(T(q,t),t)) = F r (q) = 0. (29)
Using the above definitions and the relations

∇F(p,t) = R(t)∇F r (T ˜(p,t)) (see appendix B) and F r ′ T = ∇F r , F . (p,t) = F r ′(T ˜(p,t))T ˜. (p,t) = F r ′(T ˜(p,t))(-R (t) T V (p,t)) = -R (t)F r ′(T ˜(p,t)) T T V (p,t) = -∇F(p,t) . V (p,

Appendix D: Parametrically Defined Surfaces

Assume that at time t 0 , the surface normal of S at the point S(u 0 ,v 0 ,t 0 ) is colinear with the z axis. As in appendix C, we let f(x,y,t) explicitly describe S near S(u 0 ,v 0 ,t 0 ). The definition of f is f(S x (u,v,t),S y (u,v,t),t) = S z (u,v,t).

(53)

The existence of f follows from the implicit function theorem. Thus, the first partials of f may be expressed in terms of partial derivatives of the parametric function S by solving a simple linear system. The second partials of f are obtained by differentiating the system of equations (54) and solving another linear system. A normal vector S N (u 0 ,v 0 ,t 0 ) to S(u 0 ,v 0 ,t 0 ) at time t 0 is given by S N (u 0 ,v 0 ,t 0 ) = ∂u ∂S (u 0 , v 0 , t 0 ) × ∂v ∂S (u 0 , v 0 , t 0 ). (55)

Suppose that T(u,v,t) is the parametric function for a second surface, with T N defined accordingly. Let the parametric coordinates of the extremal points p a and p b (where p a is on S and p b is on T) be (u a ,v a ) and (u b ,v b ). The analogue to equation ( 9) is then 9)) to find (u a ,v a ) and (u b ,v b ). The extremal points p a and p b are then cached in terms of their parametric coordinates

Figure 1 .

 1 Figure 1. Simulator control flow.

a < 0

 0 if A and B are inter-penetrating near p c . a = 0 if A and B are touching near p c or a > 0 if A and B are separate near p c or (1)

Figure 4 .

 4 Figure 4. The extremal distance and extremal points.

Figure 5 .

 5 Figure 5. Implicit surface description of A and B in terms of F and G.

  b ,t) ∇G(p b ,t)∇G(p b ,t)

Figure 6 .

 6 Figure 6. χ χ expressed in terms of the extremal points. When A and B are disjoint, ∇G(p b ,t) is pointed in the same direction as p ap b ; hence χ(t) is a positive multiple of the distance || p ap b ||. Similarly, when A and B inter-penetrate, ∇G(p b ,t) points in the opposite direction of p ap b and χ(t) is a negative multiple of the distance || p ap b ||. Thus, | χ(t) | is the extreme distance at time t scaled by ||∇G(p b ,t)|| and equation (8) defines a valid characteristic function χ.

Figure 7 .

 7 Figure 7. Movement of p a and p b in both world and body space.

)

  Similar expressions give the derivatives for p a y , p a z , p b x , p b y and p b z . Appendix E discusses possible ill-conditioning of the Jacobian matrix J. If we are given the extremal points p a and p b , λ 1 and λ 2 are easily determined and p . a and p . b are easily calculated (assuming the needed derivatives of F and G are at hand). However, an expression for χ .. involves the symbolic computation of p .. a and p

Figure 9 .

 9 Figure 9. Sublinear time verification of disjointness based on cached witnesses.

E 2 :E 1 :J

 21 (T (u b ,v b ,t) -S (u a ,v a ,t)) + λ 1 T N (u b ,v b ,t) = 0 → . S N (u a ,v a ,t) + λ 2 T N (u b ,v b ,t) = 0 → (56)As in section 5, we may regard u a ,v a ,u b and v b as functions of time; derivatives are then given byu . a = -J ∂(t,v a ,u b ,v b ) a ,t,u b ,v b ) = ∂(u a ,v a ,u b ,v b ) ∂(E 1 ,E 2

  y) (p a x ,p a y , f(p a x ,p a y ,t)) (p b x ,p b y ,g(p b x ,p b y ,t)) This change of coordinate systems and functions reduces the number of variables from 8 to 4; this allows a formulation in terms of 4 × 4 matrix determinants as opposed to the 8 × 8 matrices of the previous section. While this is an improvement, calculating the derivatives of the determinants constructed from equation (42) is still a formidable challenge; however, they are no longer needed. The real advantage of the new formulation is that the second derivatives of p a x , p a y , p b x and p b y are not required when computing χ .. . Equation (47) of appendix C shows how the second derivatives of p a and p b drop out of the expression for

	z	y
		x
	Figure 8. Side view of the implicit surfaces F and G expressed
	explicitly by f and g.
	equation (42)). dt 2 d 2 f and d 2 dt 2 expression for χ g. The final result is a fairly simple symbolic .. .

Table 1

 1 gives a rough indication of the running time of two simulations (figures 10 and 11) on a Hewlett Packard 835 workstation. A "cache miss" means that a new witness was computed from scratch, while a "cache hit" means that a previously cached witness was successfully updated to a witness for the current time step. The first simulation had 97 polygons and 6 curved surfaces and encountered 60 discontinuities while time stepping. The second simulation had 89 polygons and 102 curved surfaces and encountered 343 discontinuities while time stepping.

	Figure	Total No.	Cache	Cache	CPU
		Time Steps	Hits	Misses Minutes
	10. Jack	1,475	5,243	42	2.1
	11. Dice	4,162	345,793	1,384	78.6

Table 1 . Running times and caching effectiveness.

 1 If the linear and angular velocities of the rigid body motion T at time t are v (t) and ω(t) then c .

	so that
	T ˜(T (p,t),t) = T (T ˜(p,t),t) = p.	(24)

  : ∇f(p a x ,p a y ,t) -∇g(p b x ,p b y ,t) = 0 Since D 1 and D 2 do not involve p a z and p b z , let p a and p b denote (only for the remainder of this appendix) the column vectors (p a x ,p a y ) T and (p b x ,p b y ) T at time t. Similarly, let p ... , we must doubly differentiate f and g with respect to time. Using the chain rule,dt d f(p a ,t) = f ′(p a ,t)p .From equation (36) and the fact that ∇F(p a ,t 0 ) is colinear with the z axis, ∇f(p a ,t 0 ) T = (0, 0) T . ,t 0 ) which is in turn a linear function of the contact forces. Thus, χ .. (t 0 ) is itself a linear function of the contact forces.

			∂(D 1 ,D 2 )					∂(D 1 ,D 2 )
	p . a x = -	J ∂(t,p a y ,p b x ,p b y )	, p .	a y = -	J ∂(p a x ,t,p b x ,p b y )	(43)
	(and similarly for p . b x and p .	b y ) where
	∂(D 1 ,D 2 ) J = ∂(p a x ,p a y ,p b x ,p b y )	.				.		(44)
	The derivatives of p a z and p b z are simply f g . (p b x ,p b y ,t).	(p a x ,p a y ,t) and
	a = (p . a x ,p . derived	a y ) T and p .	b = (p . b x ,p .	b y ) T at time t. In section 6, we	∂x ∂f = -	∂z ∂F ∂x ∂F	, ∂y ∂f = -	∂z ∂F ∂y ∂F	and ∂t ∂f = -	∂F ∂F ∂t ∂z	(36)
	χ(t) = f(p a ,t) -g(p b ,t).						(45)	and similarly for g. Second derivatives of f and g are obtained by repeated differentiation of equation (35).
	To obtain χ									Using f and g, condition E 4 of equation (9) may be written
												componentwise as
	Then dt 2 d 2	. (p a ,t) f(p ..	(p a ,t). (46) (47)	g(p b x ,p b y ,t) p b y p b x from which we obtain -f(p a x ,p a y ,t) p a y p a x	+ λ 1	∂z ∂G (p b ,t) ∂y ∂G (p b ,t) ∂x ∂G (p b ,t)	=	0 0 0	(37)
		= p .										λ 1 =	∂z ∂G (p b ,t) f(p a (38)
												Multiplying equation (37) by -1 and using equation (36) allows
											(48)	us to rewrite E 4 as
	t) = -F′′(p,t)V (p,t) + ω(t) × ∇F (p,t). Differentiating equation (30) with respect to time, F .. (p,t) = -∇F . (p,t) . V (p,t) + ∇F(p,t) . V . (p,t) . Since V (p,t) is the point velocity of p in its rigid body frame, (31) (32) V . (p,t) is the point acceleration of p and is a linear function of the forces in the system. Thus, F .. is also a linear function of force. This yields dt 2 d 2 f(p a ,t 0 ) = p . a T f ′′(p a ,t 0 )p . a + 2∇f . (p a ,t 0 ) T p . a + f .. (p a ,t 0 ). (49) A similar expression is obtained for g. Thus, neither p .. a nor p .. b is required for the symbolic computation of χ .. (t 0 ). Equation (49) is .. .. a linear function of F by its last term, f (p a	p a y -p b y p a -	∂y ∂g (p b x , p b y , t ) ∂x ∂g (p b x , p b y , t )	= 0 → .	(40)
												Conditions E 2 and E 3 are no longer required since f and g give
												explicit definitions of the surfaces. Using the notation
												∇f T = ( ∂x ∂f , ∂y ∂ f ) = f ′	(41)
												p a x , p a y , p b x and p b y may be defined as the solution to
												p a x -p b →	(42)
												D 2 : p a y -p b y
												The implicit function theorem for simultaneous equations gives
												the result

t). (30) By differentiating equation (30) with respect to p and using ∇V(p,t) = ω * (t) ∇F . (p,t) = -F′′(p,t)V (p,t) -∇F (p,t) . ∇V (p,t) = -F′′(p,t)V (p,t) -∇V (p,t) T ∇F (p,t) = -F′′(p,t)V (p,t) -ω * (t) T ∇F (p,y ,p a y ,t)g(p b x ,p b y ,t) . xp b x + (f(p a x ,p a y ,t)g(p b x ,p b y ,t)) ∂y ∂g (p b x , p b y , t ) ∂x ∂g (p b x , p b y , t ) = 0 → . (39) This new condition has one less equation than E 4 because of the reduction of variables from F to f. In a similar fashion, λ 2 of condition E 1 is eliminated and condition E 1 is rewritten as ∂y ∂f (p a x , p a y , t ) ∂x ∂ f (p a x , p a y , t ) x + (f(p a x ,p a y ,t)-g(p b x ,p b y ,t))∇g(p b x ,p b y ,t) = 0 → . D 1 a + f . (p a ,t) = ∇f(p a ,t) T p . a + f a ,t) = f ′′(p a ,t)p . a + ∇f . (p a ,t) T p . a + ∇f(p a ,t) T p .. a + f . ′ (p a ,t)p . a + f a T f ′′(p a ,t)p . a + 2∇f . (p a ,t) T p . a + ∇f(p a ,t) T p .. a + f .. (p a ,t).

  Differentiating equation (53) with respect to u, v, and t, we obtain the system

	∂x ∂f	∂u ∂S x + ∂y ∂f	∂u ∂S y	= ∂u ∂S z
	∂x ∂f	∂v ∂S x + ∂y ∂f	∂v ∂S y	= ∂v ∂S z	.	(54)
	∂x ∂f	∂t ∂S x + ∂y ∂f	∂t ∂S y + ∂t ∂f = ∂t ∂S z

  Given the derivatives of the parametric coordinates of the extremal points, we can calculate the needed derivatives of p a x , p a y , p b x and p b y (for equation (49)). From equation (50) and since p a = S(u a ,v a ,t), y .For contact determination, equation (56) is solved (in place of equation (

			)	.	(58)
	p . a x =	∂S	∂u

x (u a ,v a ,t) u . a + ∂v ∂S x (u a ,v a ,t) v . a + ∂t ∂S x (u a ,v a ,t)

(59

) and similarly for p . a y , p . b x and p . b

We stress the fact that χ is a local function and need only be valid for an arbitrarily small open neighborhood of t 0 .

Our simulator uses an implementation of Lemke's algorithm described in[START_REF] Tomlin | Robust implementation of Lemke's method for the linear complementarity problem[END_REF].

If one shape lies completely inside the other, the minimum distance between them is positive, but intersection between the surfaces must occur first.

Initial techniques such as hierarchical bounding volumes and spatial subdivision can be used to limit the number of pairs of primitives considered by this step. This initial step also benefits greatly by exploiting coherence.

Note that convexity is crucial as this argument does not hold for concave objects. The collision detection problem for concave objects is a considerably more difficult problem than for convex objects.

For parametrically defined surfaces, we replace equation (9) with a non-linear equation that defines the extremal points in terms of parametric coordinates; see appendix D for details.

[START_REF] Forsythe | Computer Methods for Mathematical Computations[END_REF] The points p a and p b that maximize the distance between the surfaces satisfy equation (9), but they are not the solution we are interested in. Parametric surface have multiple solutions of equation (56) for p a and p b in terms of their parametric coordinates.

Acknowledgements

This research was funded by an AT&T Bell Laboratories PhD Fellowship and two NSF grants (#DCR8203979 and #ASC8715478). Simulations were performed on equipment generously donated by the Hewlett Packard Corporation. Some displays were computed using an AT&T Pixel Machine, donated by AT&T.

Appendix A: Rigid Body Motion of Surfaces

The formulation for the constraint function χ between curved surfaces models the surfaces as time-varying implicit functions F (p,t), F : R 3 × R → R. A point p is on the surface of F at time t if F (p,t) = 0. We will represent the shape F (p,t) in terms of a rest function F r (p) and a rigid body transformation T(p,t).

Let F r (p) be a time-invariant function from R 3 to R; F r defines a rest shape by the equation F r (p) = 0. Let a rigid body motion be defined by the affine transformation T : R 3

where R (t) is a 3 × 3 rotation matrix and c (t) is a point in R 3 .

Define T ˜as

(23) 8 For polyhedral contact, χ of equation ( 3) is a linear measure of distance, and equation (21) converges in one step to t c for constant accelerations. For curved surfaces, χ of equation ( 8) is weighted by ||∇G||. For reasonably scaled functions, ||∇G|| does not vary much over the range of the prediction.

Appendix B: Rotation of Coordinate Systems

In section 6, the coordinate axes were rotated so that ∇G(p a ,t 0 ) would be colinear with z ˆ, the unit z axis vector. Let R be the change of basis matrix; R is a rotation matrix such that

Note that R is constant with respect to time. Let ∇G 0 and G 0 ′′ be the derivatives of G in the original coordinate system. In the rotated coordinate system, the derivatives become

The derivatives G . and G .. are invariant under rotation. See Goldstein [START_REF] Goldstein | Classical Mechanics[END_REF] for further discussion. 9 Given a vector a ∈ R 3 , a * is the 3 × 3 (anti-symmetric) matrix such that for any

Appendix C: Derivation of χ ..

In section 6, a change of functions is introduced by writing the implicit functions F and G in terms of explicit functions f and g near the extremal points p a and p b . This change of functions is made in a coordinate system where both ∇F(p a ,t 0 ) and ∇G(p b ,t 0 ) are colinear with the z axis. The explicit functions f and g are defined by F(x,y, f(x,y,t),t) = 0 and G(x,y,g(x,y,t),t) = 0.

(

The existence of f and g is seen by the implicit function theorem. By differentiating equation ( 35) we obtain Appendix C derived χ .. by defining explicit surfaces f and g that modeled the implicit surfaces F and G near extremal points. χ

.. was then written in terms of the derivatives of the extremal points and the derivatives of f and g; in turn, the derivatives of f and g were expressed in terms of the implicit functions F and G. Given two parametric surfaces, S and T, the same thing can be done. We will first show how to express the derivatives of f and g in terms of the parametric functions S and T. We will then show how the extremal points and their derivatives are defined for parametric functions. This enables χ .. to be computed as shown in appendix C.

Let a time-dependent parametric function S(u,v,t), S : R 2 × R →R 3 , define a surface. We will write S in terms of three component functions S x , S y and S z as

(50) As in appendix A, S is assumed to be the rigid body motion of some time-invariant base shape S 0 ;

S(u,v,t) = R(t)S 0 (u,v) + c(t).

(51)

Differentiating with respect to t,

(u a ,v a ) and (u b ,v b ). Equation ( 17) is replaced by

and similarly for v a , u b and v b .

Appendix E: Ill-conditioned Jacobian Determinants

Given the formulation of p . b diverge. Thus, the overall derivation for forces between pairs of curved objects is meant to be applied to pairs where at least one of the objects really is curved at the contact point. In regions where both curvatures approach but do not attain zero, the contact force can change rapidly. If both curvatures achieve zero simultaneously, the contact force may be unbounded (although its integral is not). The former case can be handled (inefficiently) by taking sufficiently small time steps. The latter case would seem to require some sort of analytic integral of the contact force; this is a topic for further research.