A Geometrically Exact Contact Model for Polytopes in Multirigid-Body Simulation Free

There is a wide spectrum of applications for simulation, but an understandable bias toward real-time interactivity has bred a culture that leaves concerns of physical accuracy in a dark corner. As a result, applications which require some level of physical fidelity from simulation are left with a relatively underdeveloped set of tools. The major focus of the work presented herein is on improving the simulation accuracy, and consequently stability, without significantly increasing the computational cost.

The tradeoff of speed versus accuracy in simulation is a classic one, and many applications simply do not gain from improved physical accuracy. The gaming community is naturally concerned with methods that are fast and reliable, but not necessarily accurate [1–3], inspiring the term “game physics.” Many games, in fact, intentionally incorporate nonrealistic physics into their play. The graphics community, which is constantly pushing the boundaries of realistic-looking visual effects, does not offer many improvements to physically accurate simulation. Although there are many examples in a long history of “physics-based” methods in use in animation [4–8], there is really no need for such methods to do more than approximate some convenient alternative version of reality. Indeed, the goal in many graphics applications is to present the viewer with something that “feels” correct.

Particularly in the field of robotics, physical fidelity is a pressing concern and is especially challenging when simulating systems experiencing intermittent frictional contact, which is the case when robots work in unstructured, cluttered environments. Robotics experiments tend to be exceedingly costly and time consuming, so that the ability to perform realistic experiments in simulation would be a great benefit. Although work has been done regarding validation of simulation [9–12], the problem of measuring and quantifying the physical accuracy of simulation is an open issue that is challenging to solve. Comparison of a particular simulator's results with those obtained from a commercial simulator is meaningless without proof of the commercial simulator's accuracy. Even validation of a simulator with a physical experiment can be misinterpreted when it takes the tuning of several parameters to match a single or small number of experiments. Such an attempt at validation says little, if anything, of the simulator's ability to predict a different experiment.

Simulation can be a powerfully useful tool if shown to be accurate. A simulation requires accuracy in order to make any claims as to its ability to predict behavior in the physical world, yet the growing number of variants on old simulation methods all have at least this in common: they sacrifice physical fidelity at multiple stages in order to achieve computational speed-up. The ideal simulation tool would be both accurate and fast, but when we sacrifice accuracy, we are necessarily sacrificing the usefulness of that tool.

The polytopal exact geometry (PEG) contact model presented herein is a step toward improving the physical fidelity of simulations of polytopal rigid bodies in intermittent contact. In most popular simulation tools, in order to make it easy to formulate the nonpenetration constraints, the shapes of the bodies are effectively altered. Specifically, facets of polytopes are represented as hyperplanes (e.g., faces of polyhedra in 3D are represented as planes). As will be illustrated below, the price for this convenience is very high, as it leads to nonphysical simulation artifacts, such as huge constraint forces that eject bodies from the scene, or worse, cause the simulation's solvers to fail completely. By contrast, PEG exactly represents the boundaries of polytopal bodies, including their finite extent, and thus avoids these problems. To achieve this, the PEG contact model is more complicated than the standard contact model. With the same time-step sizes, simulation with PEG takes approximately twice as long as simulation with the standard contact model. However, its improved accuracy avoids nonphysical artifacts, including penetration-induced instabilities that can crash a simulation. In our opinion, this speed reduction is a penalty worth paying.

In the remainder of this section, we review multibody dynamics and some of the approaches used in the simulation. In Sec. 2, we detail the deficiencies of standard contact modeling practices and begin to derive the PEG contact model by introducing new approaches to cases that standard contact models fail to accurately represent. In Sec. 3, we define the general formulations of three classes of nonpenetration constraints and provide abstractions of them in order to model contact in a clean way. In Sec. 4, we precisely define contact constraint applicability and feasibility, and establish geometric tests to determine which contacts to include in PEG. In Sec. 5, we utilize the abstractions from Sec. 4 in order to define four fundamental contact constraints for spatial problems, which reduce to two when specialized to the planar case. We then demonstrate how to formulate a time-stepping subproblem that incorporates PEG. In Sec. 6, we compare the performance of multibody simulations using a popular time-stepping method with the standard contact model against the same time stepper with the PEG contact model. Finally, Sec. 7 discusses our results and possibilities for future work.

The two main classes of simulation methods for multibody systems are known as event-driven and time-stepping. In event-driven methods, the simulation is monitored for impact and separation events as well as changes between sticking and slipping (see Ref. [13]). For systems with small numbers of contacts, this method is preferred, because it uses an adaptive step size to accurately pinpoint and process the contact state transitions. Its downfall is that in systems with many bodies and contact events, the size of the time step becomes so small that computation times become impractically large. On the other hand, this is precisely the circumstance in which time-stepping methods have their greatest value. Time-stepping methods use large time steps that can advance the simulation across a number of contact events in one step, such that at the end of the time step, Newton's laws are satisfied and bodies do not overlap [14].

Time-stepping simulation methods utilizing the complementarity problem (defined formally below) were introduced by Moreau and Panagiotopoulos [15]. In these methods, the data needed to advance the simulation by one time step are obtained as the solution of a complementarity problem, the so-called time-stepping subproblem. Stewart and Trinkle introduced a modification of this problem that incorporated stabilization of nonpenetration constraints [16], and they provided a proof that body trajectories approximated by the solutions to a sequence of these subproblems converged to a solution of the instantaneous-time formulation as the step size goes to zero [17]. Anitescu and Potra [18] included bilateral constraints and showed that the subproblem with linearized friction cones is always solvable.

Physics engines such as DART [19] use semi-implicit, velocity-based time-stepping subproblems, formulated as linear complementarity problems with nonpenetration constraint stabilization. However, because it uses standard contact modeling, it can suffer from inaccuracy and instability as discussed above. Chrono::Engine [20] is built on top of the Bullet collision detection engine [21], which reports penetration depth and gap distance, implying that penetration is not avoided. Across the gamut of available physics engines, simulations with bodies modeled as spheres or unions of spheres are popular for benchmarking, especially when there are large numbers of contacts. This is because when simulating polytopes in contact using the standard contact model, simulations can fail due to nonphysical, overconstraint imposed on the polygons. These overconstrained situations, henceforth referred to as standard-model traps (SM traps), can cause bodies to obtain huge nonphysical velocities or result in complete failure of the simulation.

In this section, we will introduce an instantaneous-time model of polyhedral bodies in intermittent contact, and then discretize it to obtain the discrete-time counterpart. To simulate the motions of the bodies, the discrete-time model (i.e., the time-stepping subproblem) must be formulated and solved repeatedly, each time advancing the simulation forward in time by the chosen time step. Note that the development herein excludes contact friction in order to maintain a sharp focus on the new contact model, which depends only on body geometry and not friction. Friction can, however, be easily incorporated by augmenting the time-stepping subproblem with one's favorite friction model, several of which are detailed in Refs. [22,23].

Before formulating PEG and a time-stepping subproblem that employs it, it is important to introduce the complementarity problem [24].

The complementarity problem offers a convenient way to model contact because it naturally represents disjunctive conditions. The two key cases in multibody dynamics are: two bodies are in contact OR not in contact and an existing contact is sticking OR sliding. There is a vast body of research on complementarity problems and methods for their solution [25,26].

where $⊥$ represents the orthogonality of condition 3 in Eq. (1). It is common to refer to all three of these conditions jointly as a complementarity condition.

There are many general-purpose solvers available for complementarity problems [27] and in particular for linear complementarity problems (LCPs) [24,28]. In an LCP, $f(z)$ is a linear function of $z$⁠, i.e., $f(z)=Az+b$⁠, where $A$ is a known constant square matrix, and $b$ is a known constant vector of compatible size. Much theory for LCPs has been developed and used to design reliable solvers. It is known, for example, that if $A$ is positive definite, then the LCP is equivalent to a convex quadratic program, and therefore has a unique solution and can be solved efficiently. If $A$ is positive semidefinite, then a solution can be found efficiently also, but some unknowns may not have unique solutions. In multibody dynamics, this case arises only for frictionless systems, and for these, the body motions have unique solutions, but the contact forces do not.

If friction is included in the model, these nice properties are no longer present. Time steppers using polyhedral friction cones still yield LCPs, but $A$ becomes indefinite, and correspondingly, the solutions are not unique; they might not exist and when they do, there can be infinitely many. In the special case when all contacts are sticking and Coulomb friction cones are approximated by four-sided pyramids, the solutions are known to exist [29]. However, the problems are more challenging to solve if the friction model is nonlinear, as is the case for quadratic Coulomb friction, where no existence results are available. Nonetheless, some solvers have success solving the nonlinear complementarity problems (NCPs) that result. We utilize the general-purpose path solver for many NCP-based simulations. Work has been done to create specialized solvers for multibody dynamics [30].

Consider two bodies $A$ and $B$ as depicted in Fig. 1. When the signed gap distance ψ_n is positive (as in the figure), the bodies are separated. When ψ_n is zero or negative, they are in contact or interpenetrating. Ideally, the latter is not permitted, that is, $ψn≥0$⁠. Generally, a contact between two bodies may be represented as a five-tuple $C={Aid,Bid,pa,n̂,ψn}$⁠, where $Aid$ and $Bid$ are body indices or identifiers for $A$ and $B$⁠, respectively, $pa$ is the contact point on $A$ in world coordinates, and $n̂$ is the contact normal vector in the direction of $A$ onto $B$ by convention. The point of contact $pb$ on $B$ is available by $pb=pa+ψnn̂$⁠. The vectors $ra$ and $rb$ span from the center of mass of each body to the point of contact on that body. While not required, we use the convention that the origin of a body's body-fixed reference frame is coincident with its center of mass. Vectors $t̂$ and $ô$ form an orthonormal basis with $n̂$⁠.

where $rij$ is the vector from the center of mass of body i to the point of application of the force due to the jth contact as in Fig. 1 above. The generalized body configurations are represented by $q=[xTqT]T$ for position $x$ and unit quaternion q [26].

for a contact where $ψ(q,t)$ is the signed gap distance which is dependent on the configuration of the bodies $q$ and time t, and λ_n is the magnitude of the force applied in the normal direction $n̂$ of the contact. Just as we will do for $G$ and $M$ below, we drop the functional dependence $(·)$ to simplify notation.

The orthogonality condition of Eq. (4) enforces for each contact that if the force being applied is positive (compressive), then the gap distance must be zero, and if the gap distance is positive, then the force must be zero. The inequalities imply that the force cannot be tensile and the bodies cannot overlap. In the frictionless case, the continuous-time dynamic model of rigid bodies in contact is composed of Eqs. (3) and (4).

where $p$ denotes an impulse over the time step, i.e., $p=λh$⁠.

Equations (8) and (9) define the time-stepping subproblem of Stewart and Trinkle [16] without friction. The first row of Eq. (8) is the discrete-time Newton–Euler equation. The second row is merely a definition used to shorten the statement of the normal complementarity condition (9). Importantly, Eq. (9) defines the standard model of unilateral contact used widely in multibody simulators. The term $ψn/h$ appearing in the standard model acts as a constraint stabilization term, which if removed, allows interpenetrations to occur and then remain indefinitely. Regardless of whether term $ψn/h$ is removed or retained, the inequality on the left-hand side of condition (9) is the source of SM traps that we will see PEG eliminates.

Since Eq. (10) does not preserve length, each quaternion should be renormalized after each time step.

One approach to multibody simulation is to allow interpenetrations (and other constraint violations) to occur and then apply whatever hacks are necessary to catch poorly defined geometric configurations that inevitably result. Our approach attempts to prevent nonphysical configurations from occurring. In order to achieve this, we must detect potential contacts before interpenetration occurs.

Contact constraints have been traditionally modeled as fundamentally unilateral, however, representing features as infinite half-spaces at points of potential contact generates nonphysical behavior at corners where multiple finite features meet and terminate. Consider a vertex v near a corner as depicted in Fig. 2. When the collision detection routine identifies both $C1=C(v,e1)$ and $C2=C(v,e2)$⁠, the vertex becomes erroneously trapped in the convex subspace of the free space surrounding the corner if both contacts are enforced. We refer to this as the standard-model trap (SM trap).

Figure 3 shows a simple situation in which two SM traps collude to generate a time-stepping subproblem with no solution. A common heuristic in collision detection is to choose contacts with deepest penetration. These contacts are then used to formulate a dynamics problem to move the simulation forward in time. Since the standard contact model requires that by the end of the next time step, both of the identified contacts must be non-negative, an impossible problem is generated. Poor contact choices lead to poorly defined dynamics problems, resulting in the all too familiar “popping” or “exploding” bodies in a simulation. Note that allowing bodies such as those in Fig. 3 to first interpenetrate before identifying contacts does not avoid this erroneous contact problem, as the same heuristic results in similarly poorly defined dynamics. Fundamental errors such as these are what lead to common ad hoc “corrections” to these problems, such as limiting the maximum magnitudes of contact forces.

allowing v to penetrate either half-space corresponding to e₁ or e₂, but not both.

Lemma 1.1. Given $a,b∈ℝ,max(a,b)=a+max(b−a,0)≥0$⁠.

Proof.

• Case $a≥b⇒a+max(b−a,0)=a+0=a$⁠.

• Case $a≤b⇒a+max(b−a,0)=a+(b−a)=b$⁠.◻

Lemma 1.2.Given$a,b∈ℝ,b=max(a,0)⇔0≤b−a⊥b≥0$⁠.

Proof.

• Case $a≤0$ : $max(a,0)=0⇒b=0⇒0≤b−a⊥b≥0$⁠.

• Case a > 0: $max(a,0)=a⇒b=a⇒0≤b−a⊥b≥0$⁠.

• Case $b−a=0,b≥0$ : $b=max(a,0)=a$⁠.

• Case $b=0,b−a≥0$ : $a≤0⇒max(a,0)=0=b$⁠.◻

To give some insight into the slack variable c and its relationship with ψ₁ and ψ₂, consider the figure below (Fig. 4) and the value of c on each side of the bisector where $ψ1=ψ2$⁠. A vertex in the free space right of the bisector has $ψ1>ψ2⇒ψ1−ψ2>0$⁠, which requires c = 0 to satisfy the complementarity condition of Eq. (12). A vertex to the left of the bisector has $ψ1<ψ2⇒ψ1−ψ2<0$⁠, which requires $c=ψ2−ψ1$⁠.

We now wish to include potential forces λ₁ and λ₂ while ensuring that only one is nonzero since allowing both to be simultaneously nonzero would be nonphysical.

Lemma 1.3.Given$a,b∈ℝ, b=|min(a,0)|⇔0≤b+a⊥b≥0$⁠.

Proof.

• Case $a≤0$ : $|min(a,0)|=−a=b⇒b+a=0⇒0≤b+a⊥b≥0$⁠.

• Case a > 0: $|min(a,0)|=0=b⇒(a+b)b=0⇒0≤b+a⊥b≥0$⁠.

• Case $b+a=0,b≥0$ : $a≤0⇒|min(a,0)|=−a=b$⁠.

• Case $b=0,b+a≥0$ : $a>0⇒|min(a,0)|=0=b$⁠.◻

Lemma 1.4. Given$a,b∈ℝ, a=0,b≤0⇔(max(a,b)≥0)∧(max(a,b)+|min(a,0)|=0)$⁠.

Proof. The only way to satisfy both $max(a,b)≥0$ and $|min(a,0)|≥0$ is for both to reach equality. When $|min(a,0)|=0$⁠, it is possible that $a≥0$⁠, so in order to have $max(a,b)=0$⁠, we must have $a=0,b≤0$⁠. The argument for the reverse of a and b is similar.◻

Essentially, we have taken a set of contacts, in this case C₁ and C₂, and replaced a logical AND condition of $C1∧C2$ with a NAND condition of $C1∧¯C2$⁠. Geometrically, this corresponds to replacing line representations of a body with line segments. At this point, Eqs. (12)–(14) correspond to the locally nonconvex (LNC) model proposed by Nguyen [32], which models the case for a single vertex against multiple edges. Although this model has been shown to contribute some improvement in physical fidelity over other methods [33], stopping here would correct the standard-model trap at the expense of permitting interpenetrations [34,35]. To illustrate this error, let us continue with a less trivial case.

Consider two bodies approaching as depicted in Fig. 5. In this configuration, our collision detection routine identifies four potential contacts. These are $C1=C(va,eb1)$⁠, $C2=C(va,eb2)$⁠, $C3=C(vb,ea1)$⁠, and $C4=C(vb,ea2)$⁠. We should not simply constrain all of these contacts, for this would result in a dual standard-model trap. It is possible to make a guess as to which contacts to enforce, but this inserts the likely possibility of generating contact forces where there are none. Of course, we could also ignore all four contacts and allow interpenetration, but this too is clearly a nonphysical approach and leads to instabilities.

but we make the important observation that this model erroneously allows for ψ₁ and ψ₄ or ψ₂ and ψ₃ to be simultaneously violated. An example of such a consequence is depicted in Fig. 6. The model of Eq. (15) can be completed by including constraints to prevent such interpenetrations. That is, we additionally constrain $max(ψ1,ψ4)≥0$ and $max(ψ2,ψ3)≥0$ with

Equations (15) and (16) together represent our first example of a geometrically accurate contact model. We will refer to the constraint type in Eq. (15) as intercontact constraint and the type in Eq. (16) as cross-contact constraint. These, along with unilateral constraint, form the set of three contact constraints necessary for our model.

In order to resolve the issues discussed in Sec. 2, we had to change the way we formulate contact constraints. In this section, we define three fundamental contact constraints in their general form and begin to describe how they can be applied.

where λ_i is the force corresponding to contact C_i.

An intercontact constraint (⁠$I$ -constraint) is that which was demonstrated in Eq. (15) for preventing the standard-model trap. An $I$ -constraint deals with contacts of a single feature of one body against multiple features of another body. Here, we will define the general formulation of $I$ -constraint for a vertex of $A$ against multiple features of $B$ where a subset of these potential contacts may feasibly result in a contact force at the end of the current time step.

Consider two sets of contacts: $Cι$ which is composed of all contacts which could feasibly result in a contact force, and $Cκ$ which contains all contacts that could not feasibly result in a contact force but are necessary to accurately represent the contact geometry. An $I$ -constraint seeks to

• constrain at least one gap distance $ψ{ι,κ}$ of ${Cι,Cκ}$ to be non-negative

• permit a nonzero contact force $λιi>0$ for at most one contact in $Cι$

• only permit $λιi>0$ when $ψιi=0$ and all other $ψι,κ$ are nonpositive

Equations (20) and (21) enforce that at least one contact in ${Cι,Cκ}$ is non-negative. Equations (22) and (23) allow at most one contact force to be nonzero, and only permit this force to be nonzero when all other gap distances are nonpositive. Solving such a set of constraints is possible as a linear program with complementarity constraints (LPCCs) [36], or by introducing an auxiliary variable to convert inequalities (21) into a linear complementarity conditions before using an LCP solver [32].

This greatly reduces the number of constraints and removes the inequality of Eq. (21) to leave a pure LCP. We see now the incentive for choosing our contact set with a single primary contact $Cι$⁠, possibly using heuristics like those described in Sec. 4.

In contrast to $I$ -constraint which deals with a single feature of $A$ against multiple features of $B$⁠, a cross-contact constraint (⁠$X$ -constraint) deals with sets of contacts that may not depend on the same features, but could result in interpenetration if simultaneously violated, e.g., Eq. (16).

Given the sets of contacts $Ca$ and $Cb$⁠, a $X$ -constraint attempts to constrain either

• at least one non-negative ψ_a in $Ca$

  OR

• at least one non-negative ψ_b in $Cb$

and can write the LCP in multiple forms. Conveniently, $X$ -constraints can be written in the form of Eqs. (20) and (21). That is, $X(Ca,Cb)$ is similar to $I(Ca,Cb)$⁠, but excludes Eqs. (22) and (23) which involve contact forces.

Although a $X$ -constraint necessarily deals with contacts between the same pair of bodies, the contacts need not share features. That is, the geometric interdependencies of identified potential contacts, regarding which are to be enforced as defined by an $X$ -constraint, may be between two otherwise seemingly independent contacts, e.g., opposite corners of the face of a box. Subsequently, it is conceivable for two potential contacts, neither of which have a potential force associated with them, i.e., contacts in two different sets of $Cκ$⁠, to be $X$ -constrained. In such cases, the dynamics problem requires sufficient $U$ - or $I$ -constraints so as to allow generation of forces which may contribute to achieving configurations that will satisfy all constraints. This issue, however, while an important one when employing a heuristic approach to choosing primary contacts in $I$ -constraints, is beyond the scope of this paper.

In this section, we briefly discuss the nature of interpolytope contact and define a set of geometric tests that are useful for identifying and making heuristic choices regarding potential contacts.

We define a convex polyhedron in terms of features (Fig. 7). We may refer to a vertex either as an object v or as the vector $v$ representing the vertex's position. Similarly, an edge object $e=(vt,vh)$ corresponds to a vector $e=vh−vt$⁠. A body is composed of a set of vertices V, edges E, and faces F defined by three vertices in counterclockwise order. The outward pointing normal vector $η̂$ for a given face $f∈F$ is perpendicular to f. We define two vectors $t1$ and $t2$ for each edge that are planar with f₁ and f₂, respectively, and perpendicular to $e$⁠. These vectors are defined as

for convex bodies.

As we develop these tests, we keep in mind that a contact can occur anywhere on a polytope, but the possible normal vectors are restricted by the contacted feature. Figures 8–12 depict the possible positions and orientations of contact normal vectors for each feature in two and three dimensions. Since a vertex does not have a uniquely defined normal, we use the normal from the edges that connected to the vertex. It is important to realize that when using a time-stepping method, we hope to identify contacts which will reflect these feasible contact normals at the end of the time step. Yet, many methods of contact identification use approaches that assume continuous or static body positions. This is why we geometrically relax the following set of tests that will prove useful for identifying and making heuristic choices regarding potential contacts.

Consider the simple example in Fig. 13 where the two vertices, $va1$ and $va2$⁠, of body $A$ are near an edge e_b of body $B$⁠. In order to achieve stability in such a configuration, it is clear that both potential contacts $C1=C(va1,eb)$ and $C2=C(va2,eb)$ must be considered by the dynamics formulation. If only one were included, e.g., C₁, then it is possible for the other to penetrate e_b given a large enough time step or velocity. We observe that the second contact normal n₂ is outside of the normal cone region of $va2$ and does not correspond to a physically reasonable contact normal. However, this is because n₂ is determined at the beginning of the time step; yet, we wish to prevent interpenetration for the end of the time step. In order to avoid interpenetrations, we must detect and consider such contacts as C₂, even when they are not physically feasible at the current time step. This necessity requires that we geometrically relax the normal cone region or, equivalently, relax the definition of applicability.

The time required to solve a time-stepping subproblem, and therefore the time required to compute a simulation, is strongly dependent on the number of contacts n_c. While it may seem that the number of contacts is determined by the simulation, taking a close look at an individual discrete time step, as we do in Sec. 1.2.3, reveals that to enhance simulation accuracy and reduce instability, it is beneficial to include “potential contacts” even though they increase the problem size. We slightly modify classical applicability conditions to derive a method to choose a small set of extra contacts.

where $η̂b$ is the outward unit normal of e_b, and the superscript ve denotes that the test is for a vertex and an edge.

Notice that the applicability test is dependent on $q$⁠; more specifically, it is dependent on the relative orientation of the polygons, and not their relative translation. Later, when considering potential contacts to include in a time-stepping subproblem, we will introduce a relaxed form of inequality (28) to allow for possible rotation in a time step. In the case of our example, if the top polygon rotated enough in the clockwise direction in the time-step under consideration, contact could occur first at vertex $νa2$ rather than vertex $νa1$⁠.

Whereas classical applicability is a Boolean function, our function returns a value as we may wish to use the value of applicability as a heuristic when comparing contacts.

Figure 14 depicts two edges, e_a and e_b, in contact. Just as there are vertex–face configurations that have proximity but could not result in a force, so too are there edge–edge configurations that should not be associated with potential contact forces. First, we should note that given an edge–edge contact normal $n=ea×eb$⁠, we cannot immediately say if we have the correct sign for $n$⁠. However, we can still determine edge–edge applicability in the following manner. Let

Conceptually similar to Eq. (30), Eq. (31) measures the maximal extent by which the adjacent faces of an edge are pointed “in” on another body.

Feasibility is similar to applicability, but is dependent on feature positions as opposed to orientations. The region of feasibility is essentially the region in which we anticipate a contact between a vertex v and a facet f could occur. We define vertex–edge and vertex–face feasibility similarly. Given an angle $ϕ$ and a small value δ greater than machine precision, a vertex with contact $C(va,fb)$ has feasibility if

• $va∈VR(f)$ and $ψ≥−δ$

  OR

• $va∈VR(f)$ and v_a within region defined by $ϕ$ below the point $ϵ$ below the nearest point on f

  where VR(f) refers to the Voronoi region of the facet f [41]. An example of this region is depicted in Fig. 15.

We can also write an edge–edge feasibility in 3D, which checks only to see if the nearest two points p_a on e_a and p_b on e_b are within some small distance $ϵ$⁠.

From the set of fundamental constraints presented in Sec. 3, which are written in terms of complementarity problems, abstract out another step and look at pairs of bodies by cases of feature pairs.

This abstraction layering is depicted in Fig. 16. At the lowest level is the complementarity problem which is utilized first as a model of unilateral contact with force, but further as a tool for generating logical NAND conditions on subsets of potential contacts. Sets of these low-level constraints can be entirely represented by abstracting to the three fundamental constraints of unilateral, intercontact, and cross-contact. In turn, the fundamental constraints are formulaically generated by the highest abstraction layer which looks only at what pairs of features between two bodies are near each other. In Sec. 5.2, we will see how precisely the same abstraction can be developed in 3D.

We consider two contact configurations in 2D: vertex–edge and vertex–vertex. We do not consider edge–edge because this case is redundantly covered by the other two configurations [38]. To be clear, edge–edge penetration in 2D is impossible if vertex–edge and vertex–vertex constraints are properly enforced, simply because edges defined and composed of vertices, and for two edges to come near each other implies that the vertices of at least one of those edges are nearing the other edge. We can formulaically generate the set of constraints for the two 2D configurations as given below.

In this case, the vertex is within a reasonably small distance $ϵ$ of the edge, but farther than $ϵ$ from either of the edge's vertices.

where $Cι1$ contains at least one of C₁ or C₂, $Cκ1$ may contain one of C₁ or C₂, and similarly $Cι2$ contains at least one of C₃ or C₄, and $Cκ2$ may contain one of C₃ or C₄. As is the case for all $I$ -constraints, in which contacts are included in $Cι$ or $Cκ$ these are determined using the geometric heuristics defined in Sec. 4.

There are four contact configurations in 3D: vertex–face, edge–edge, vertex–edge, and vertex–vertex. Again, we exclude certain feature pairs just as we did in 2D, namely, edge–face and face–face, because they become redundant when the other configurations are properly enforced. Edge–face is redundant since edges are composed of vertices, and for an edge to approach a face requires that either the vertices of that edge approach that face, or the edge approaches the edges at the perimeter of that face (possibly adjacent edges), reducing to either vertex–face or edge–edge, respectively, and edge–vertex in the case that the edges of the second body are adjacent. By a similar argument, face–face becomes redundant; a face, being composed of edges in turn composed of vertices, nearing another face implies that edges and vertices are approaching the other face. In all possible configurations, this results in one of the four feature combinations considered.

The vertex–face case is analogous to the 2D vertex–edge case. Given a vertex v_a near a face f_b as depicted in Fig. 17, we include

in the set of constraints.

Consider Fig. 18 where a vertex v_a of body $A$ approaches the edge e_b of body $B$⁠. There are an arbitrary number of edges (greater or equal to three) that could connect to v_a; however, we can limit the edge–edge contacts considered to a subset of those with relatively large magnitudes of edge orientation $|ea·Oeb|$⁠. Given $C1=C(va,f1)$ and $C2=C(va,f2)$⁠, the set of constraints for vertex–edge is given by

where C₁ and C₂ are distributed appropriately into $Cι$ and $Cκ$⁠. Each C_ei is an edge with positive edge orientation and each C_ej is an edge with negative edge orientation, and both sets of edges have applicability and feasibility with e_b.

In the case of a vertex v_a near another body's vertex v_b as depicted in Fig. 19, we wish to avoid trapping v_a against the faces of $B$⁠, and similarly for v_b against $A$⁠. This first requires the set of $I$ -constraints

where f_a has edges $ea1 and ea2$ connected to v_a, and f_b has edges $eb1$ and $eb2$ connected to v_b. These $X$ -constraints are enforced only for edges e_ar and e_bs, which have edge–edge applicability and feasibility. The constraint in Eq. (38) is similar to the $X$ -constraints in Eq. (36) for the vertex–edge case, but permits a vertex to pass on either side of a face instead of either side of an edge.