# A robust, efficient and time-stepping compatible collision detection method for non-smooth contact between rigid bodies of arbitrary shape

Xavier Merlhiot

► To cite this version:

Xavier Merlhiot. A robust, efficient and time-stepping compatible collision detection method for non-smooth contact between rigid bodies of arbitrary shape. MULTIBODY DYNAMICS 2007, ECCOMAS Thematic Conference, Jul 2007, Milan, Italy. [hal-03334104](https://hal.science/hal-03334104)

HAL Id: hal-03334104

https://hal.science/hal-03334104

Submitted on 3 Sep 2021

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A ROBUST, EFFICIENT AND TIME-STEPPING COMPATIBLE COLLISION DETECTION METHOD FOR NON-SMOOTH CONTACT BETWEEN RIGID BODIES OF ARBITRARY SHAPE

Xavier Merlhiot*

* CEA-LIST

Interactive Simulation Laboratory

18 route du Panorama, 92260 Fontenay-aux-Roses, France

e-mail: xavier.merlhiot@cea.fr

**Keywords:** Collision detection, contact dynamics, interactive-time multibody simulation.

**Abstract.** This paper proposes an efficient collision detection method which is compatible with time-stepping methods in the sense that it enables the robust simulation of non-smooth contact between rigid bodies with complex shapes, including industrial CAD models of various topology and in presence of conforming contact situations. It introduces a discrete representation of rigid body shapes based on dilated simplicial complexes, which generalizes the notion of triangulation to domains of arbitrary topological dimension. It also defines finite collections of point contacts between those shapes thanks to quasi-LMDs, which are defined as an extension of local minimum distances with respect to small relative rotations, between the base complexes. Smooth gap functions associated to these point contacts are defined, as well as complete and smooth generalized contact kinematics, enabling the use of non-smooth contact laws like Signorini or Coulomb. Quasi-LMDs also lead to the stable treatment of conforming contact cases. An efficient method based on 5D+1 bounding volume hierarchies for computing quasi-LMDs is presented. Finally, robustness and performance benchmarks show that our method combined with a fast time-stepping-based solver allows interactive-time simulations of complex and possibly conforming contact situations.
# 1 INTRODUCTION

On one hand, time-stepping methods such as those proposed in Refs. [2], [41] and [20], are known to be powerful computational mechanics tools for multibody dynamics or quasi-statics with many contacts and degrees of freedom. However their capacities have been mostly illustrated with applications involving simple geometric shapes, like collections of spheres, boxes and planes, or convex polyhedral models of modest complexity. On the other hand, computer graphics and real-time applications like virtual reality and haptics have motivated the design of collision detection algorithms that run efficiently on complex polyhedral models such as in Refs. [37] or [23]. But the geometric information they compute is often specialized for the use of computationally inexpensive (but less mechanically correct) contact models like penalty methods based on penetrations depth or repulsive potentials based on separation distances. In the next two sections, we will give arguments showing how their direct combination with time-stepping schemes and more mechanically correct or non-smooth contact models suffers limitations and generally leads to severe robustness issues, therefore introducing the ingredients of our method as natural and efficient alternatives.

<figure><img src="image_2.png" /><figcaption>Figure 1: Left: complex shapes from industrial CAD models. Right: a conforming contact situation.</figcaption></figure>

For computational efficiency reasons, our method relies on a discrete representation of the rigid body shapes which will be described in section 4. Since it is based on simplicial complexes, it has the capacity to represent geometric shapes of arbitrary topological dimension. This versatility is particularly useful when working with combinations of volumes, shells and beams for instance. Then the key elements that give consistency and robustness to our approach will be presented in section 5 with the introduction of the notion of *quasi-LMD*. Finally the efficient computational method proposed in section 6 will be illustrated on industrial benchmarks in section 7 featuring complex shapes and conforming contacts as shown in Fig. 1.

# 2 ON THE GEOMETRICAL NEEDS OF THE MODELING OF NON-SMOOTH CONTACT BETWEEN RIGID BODIES

In this section are reported and discussed the main geometric hypotheses that appear in the definition of the most common non-smooth contact models between rigid bodies. In order to introduce the terminology and notations used in this paper, we will begin with restating the definitions of a few basic concepts.

## 2.1 Non-interference of rigid bodies and admissible configuration subset representations

Consider a system of *N* rigid bodies in dimension three, with generalized coordinates $q(t)$ in the configuration manifold $M$. The *i*-th rigid body occupies a time-dependent spatial domain
$S_i(t) \subset \mathbb{R}^3$ called its shape, satisfying $S_i(t) = \mathcal{H}_i(q(t), t) S_i(t_0)$, where $\mathcal{H}_i$ is an at least twice continuously differentiable function from $M \times \mathbb{R}$ into $SE(3)$ and $t_0$ a reference time. We assume that $\mathcal{H}_i(q(t_0), t_0) = I$. $S_i(t_0)$ will be called the reference configuration shape of the $i$-th rigid body. It is assumed to be a connected compact volume (i.e. having topological dimension three). Considering shapes of lower topological dimension appears to us as a degenerate case in the mechanical modeling of solid bodies (in dimension three, even thin shells models include a positive thickness parameter), which may result in consistency issues in the definition of local contact models (see subsection 3.1 for examples).

Let us first examine the case of perfect unilateral constraints. With the preceding notations, the geometrically admissible configuration set $\mathcal{C}(t) \subset M$ defined by the non-penetration of the body shapes may receive the following description :

$$ \mathcal{C}(t) = \{q \in M, \forall (i, j), i \neq j \Rightarrow \mathcal{H}_i(q, t) (S_i)^\circ \cap \mathcal{H}_j(q, t) (S_j)^\circ = \emptyset\}. $$

An effective contact situation between two bodies then corresponds to the existence of points shared between the boundaries of their shapes.

If $\mathcal{C}(t)$ is tangentially regular (see section 6.1 for a definition), an abstract formulation of the perfect unilateral constraints can be (see for example Ref. [6]) that the unilateral constraint forces $f(t)$ must lie in the opposite of the normal cone to $\mathcal{C}(t)$ at $q(t)$:

$$ f(t) \in -N_{\mathcal{C}(t)}(q(t)). $$

A case of major practical importance is when the admissible domain $\mathcal{C}(t)$ can be explicitly described by the non-negativity of a finite collection of $p$ constraint functions $(g_i)_{i \in \{1, ..., p\}}$, each $g_i$ being a real function defined over $M \times \mathbb{R}$:

$$ \mathcal{C}(t) = \{q \in M, \forall i \in \{1, ..., p\}, g_i(q, t) \geq 0\}. $$

If the constraint functions are differentiable with respect to $q$ and satisfy a weak qualification hypothesis, then $\mathcal{C}(t)$ is tangentially regular and the normal cone to $\mathcal{C}(t)$ at $q$ can be described in the following way :

$$ N_{\mathcal{C}(t)}(q) = \left\{ -\sum_{i=1}^{p} \alpha_i \nabla g_i(q, t), \alpha \geq 0, \alpha \perp g(q, t) \right\}, $$

with $g = (g_1, ..., g_p)$, thus leading to gradient-type complementarity formulations of perfect unilateral constraints between rigid bodies rather than more abstract differential inclusions.

Many authors implicitly or explicitly (see for example Ref. [4], Ref. [20] or Ref. [38]) suppose that those constraint functions can be defined thanks to the notion of *gaps* associated with a finite number of *point contacts* between rigid bodies. In the following subsections, we will discuss reasonable geometric hypotheses under which such situations may arise.

## 2.2 Rigid body shapes regularity and smooth contact kinematics

In a first but restrictive definition attempt, we could consider a relative rigid body displacement $h \in SE(3)$ such that the rigid body shapes $S_i(t_0)$ and $h S_j(t_0)$ have disjoint interiors, and think of point contacts between the rigid bodies as being isolated points of $\partial S_i(t_0) \cap h \partial S_j(t_0)$. Let $P$ be such a point for $h = h^*$. Intuitively, the associated *gap function* $g_p: h \in SE(3) \to \mathbb{R}$ should be based on some kind of “signed local distance” between the shapes $S_i(t_0)$ and $h S_j(t_0)$,
and be defined in a whole vicinity of $h^*$: the gap value must be positive in case of separation (i.e locally disjoint shapes), negative in case of interpenetration (i.e. non-empty intersection of the shape interiors “at the vicinity $P$”), and zero otherwise (i.e. in case of exact contact), like in the typical situation of Fig. 2. This definition and terminology agrees with most of the existing literature. If $g_P$ is differentiable on a vicinity of $h^*$, then, thanks to the smoothness of $H_i$ and $H_j$, it can be used to define a differentiable constraint function, locally in an open subset of $M \times \mathbb{R}$.

If we suppose that the shape boundaries are twice continuously differentiable in the vicinity of a point contact $P$ obtained at relative pose $h^*$, then they must be strictly relatively convex (see for example Ref. [30]) in the vicinity of $P$. In such a situation, $g_p(h)$ may receive an adequate definition, to which we will come back in the two next subsections, based on the local geometry of the shape boundaries, that makes it a smooth function at the vicinity of $h^*$.

In contrast, point contacts located on non-smooth areas of shape boundaries, like vertex-vertex contacts between polyhedral shapes, are known to be the source of reentrant corners in $C(t)$. In those cases, $C(t)$ is not tangentially regular, hence if gap functions are ever defined for these point contacts, they cannot be smooth functions. Those situations have been identified by numerous authors: Baraff in Ref. [3] calls them “degenerate point contact cases” between polyhedral shapes and, for convenience, converts them into vertex-plane contacts by arbitrarily choosing a normal direction, Park et al. categorizes them as “singular contacts” in Ref. [32], while in Ref. [13] Glocker proposes the extension of impact laws that apply in those cases.

<figure><img src="image_3.png" /><figcaption>Figure 2: Left: the usual vision of the gap values associated to a point contact $P$. Center: LMDs between non-penetrating compact shapes. Left: typical issue concerning the locality of penetrations between non-convex shapes.</figcaption></figure>

Anyway, in the cases where a normal direction to contact does not receive a reasonable and univocal geometric definition, the use of frictional contact laws (like non-smooth Coulomb's law), or even more sophisticated or regularized contact laws, is severely compromised. More precisely, the availability of smooth point contact kinematics, including smooth gap functions, as derived in Ref. [34], and with a slightly different formalism in Refs. [9] and [44], is a prerequisite to a wide range of non-smooth and smooth contact models between rigid bodies (see for example Refs. [28], [40] or [5]).

## 2.3 Defining point contacts and non-negative gaps through LMDs

We claim that a natural and less restrictive definition of point contacts between two non-penetrating rigid bodies of non-necessarily convex shapes $S_i(t_0)$ and $h S_j(t_0)$, which are recalled to be compact volumes, should be based (see Fig. 2) on the strict local minima of the Euclidean distance function restricted to $S_i(t_0) \times h S_j(t_0)$. Since shapes are compact, such a minimum will be attained at a minimizing couple of points $(a_i, a_j) \in S_i(t_0) \times h S_j(t_0)$, called a *local minimum distance*, or simply *LMD*, between $S_i(t_0)$ and $h S_j(t_0)$. It is easy to see LMDs
between $S_i(t_0)$ and $h S_j(t_0)$ lie in $\partial S_i(t_0) \times h \partial S_j(t_0)$. LMDs consistently extends the definition of point contacts given in the preceding subsection to all the relative configurations $h$ where the shape interiors are disjoint. It also defines associated non-negative gaps without ambiguity, and is coherent with the definition of generalized contact kinematic given for convex shapes in Refs. [34] and [9]. Each resulting gap function $g$ is hence defined at least on a non-empty subset of $SE(3)$. If the boundary surfaces are twice differentiable at the vicinity of a LMD, then the resulting generalized contact kinematics are smooth (e.g. the gap function and normal direction to contact are differentiable with respect to $h$).

## 2.4 Negative gaps definition issues

The problem of defining negative gap values that extend the gap functions defined in the preceding subsection is more difficult. When two shapes are convex, with at least one them being strictly convex, then the non-intersection of their interiors implies that there exists exactly one LMD between them. Otherwise, if the shapes interiors intersect, then one could try to invoke the notion of penetration depth between convex sets $A$ and $B$, defined as:

$$ \pi(A, B) = \inf \|\tau\|, \tau \in \mathbb{R}^3, (A + \tau) \cap B = \emptyset. \quad (1) $$

Suppose now that $S_i(t_0)$ and $S_j(t_0)$ meet sufficient supplementary conditions for the following assertions to hold:

* For any $h$ in $SE(3)$, $\pi(h S_j(t_0), S_i(t_0))$ is attained for a unique translation vector $\tau^*(h)$.

* If $S_i(t_0)$ and $h S_j(t_0)$ interpenetrate for some $h$, then $S_i(t_0) \cap (h S_j(t_0) + \tau^*(h))$ is an isolated point of $\partial S_i(t_0) \cap h \partial S_j(t_0) + \tau^*(h)$, denoted $a(h)$.

With these conditions, the couple of points $(a(h), a(h) - \tau^*(h))$ may serve as a basis for defining generalized contact kinematics in interpenetration situations. If the shapes boundaries are sufficiently smooth, one can hope for the smoothness of these contact kinematics. The most typical and widely used example of shapes that satisfy all the conditions listed above is the one of balls (see for example Ref. [38]). For convex polytopes, the inf bound in (1) is a min bound (see Ref. [1]) but it is attained for some non-unique $\tau^*$ translation vector.

For non-convex shapes, defining a “local penetration depth” is difficult. Strong local convexity assumptions are necessary, and the resulting definition is necessarily limited to “small” penetrations due to the lost of locality of penetration (see Fig. 2). The domain of validity of such a local definition is difficult to estimate in practice on complex non-convex shapes, causing well-known consistency problem. Circumventing solutions that use causality principles in the definition of local penetration depth have been proposed in several approaches (see for example Refs. [3], [4], [17] and [37]). Their main weak point is that the use of history introduces hysteresis in the definition. For instance, Baraff proposes in Refs. [3], [4] to go back in time in a situation of exact contact and use polyhedral contact models defined in this configuration (see also the next subsection).

Another encountered solution for defining generalized contact kinematics in interpenetration situations relies on the notion of extreme distance between surfaces, which in general is not precisely defined (see for example Refs. [4] and [25]). The one found in Ref. [4] closely resembles a local penetration depth definition and demands both local convexity and regularity assumptions. The one found in Ref. [25] relies on discrete considerations.
## 2.5 Conforming contact situations and point contacts selection strategies

Now let us take again $h \in SE(3)$ such that the rigid body shapes $S_i(t_0)$ and $h S_j(t_0)$ have disjoint interiors. The connected components of $\partial S_i(t_0) \cap h \partial S_j(t_0)$ which have topological dimension one or two will be called conforming contacts. Even in the case of twice differentiable boundaries, an immense diversity of cases may arise, and we will restrict our discussion to planar contact regions.

Under the assumption of perfect unilateral constraints, one could be tempted to extend the notion of point contacts to non-isolated points, provided that the associate gap and contact kinematics can still receive an acceptable definition, and decide that every extreme point of the convex hull of a planar contact is a point contact, hoping that this construction, that we call selection of point contacts, will providing a satisfying set of constraint functions to describe $C(t)$. Baraff calls this procedure “restriction of contact points” in Ref. [3], and remarks in Ref. [4] that such points may exist in infinite number, like in the example of a cylinder resting on a plane along one of its bases. The same author also remarks that contacts between polyhedral shapes give conforming contact regions which are planar polyhedra, and hence yield to a finite number of selected point contacts.

However, since these points are not isolated in the intersection of the shape boundaries, it is generally not possible, except maybe in special cases, to define associated gap functions based the considerations of the preceding subsection, especially in the case of polyhedral shapes (take for example a cube resting on a plane along one of its faces), since arbitrarily small relative displacement of the contacting bodies may cause penetrations which do not remain local to each selected point contacts. Baraff proposes a solution in Ref. [3], which consists in referring to an exact contact situation and applying a contact model (vertex-plane or edge-edge) for each selected point contact. We think that the main drawback of Baraff’s approach is that it implies a penetration-withdraw operation to define the exact contact reference configuration. This operation introduces hysteresis in the definition of the point contacts and may pose compatibility problems with time-stepping methods (see 3.3 for a more detailed discussion).

We finally remark that the frequently encountered combination of such selection strategies with a non-smooth frictional point contact model, like Coulomb’s, seems to us to be an arbitrary modeling choice of frictional conforming contact between rigid bodies that has no particular mechanical justification.

## 3 EFFICIENT COLLISION DETECTION METHODS

From the point of view of computational contact mechanics, collision detection can be seen as the effective computation of sufficiently complete geometric information, as detailed in section 2, for the chosen contact models to be used in a simulation framework. Hence, apart from the computational efficiency of collision detection methods, which is of major importance for real-time or interactive simulations, the first question that can be asked is about the exact nature of the output of such algorithms. Its variable definition explain the various understanding of the term “collision detection”. The second one is about its quality : in other words, does the computed geometric information satisfy sufficient hypotheses for the computational methods used for solving the possibly non-smooth contact dynamics or quasi-statics to exhibit satisfactory numerical behavior, and hopefully convergence ?
## 3.1 Diversity of existing methods

Although determination of contact between simulated rigid or flexible bodies is a common issue to several research areas and technical applications, including computational mechanics, computer graphics and real-time mechanical simulations (for computer games, virtual reality, haptics and robotics), the two latter communities seem to have been the most productive in the field of computationally efficient methods, which are commonly referred to as *collision detection* algorithms. Due to the immense literature on the subject, we have preferred selecting some key references cited in the next paragraph. From a more complete panorama, we orient the reader in direction of several good surveys that can be found in Refs. [29], [21], [43], [27] and the book Ref. [11] which also details many important implementation issues.

Indeed, the problem of collision detection have received numerous formulations corresponding to different interpretations and levels of richness of the computed geometric information (from mere boolean intersection queries to spatial and/or temporal localization of contact), with diverse geometric definitions of contact situations: either based on interference or proximity of shapes, global separation distance, local penetration depths and vectors with possible use of convex decompositions and incremental computations (see Refs. [35], [10] and [26]), interference volumes (see Refs. [18] and [16]) or even more exotic approaches like intersection contours like in Ref. [45]. Some methods are dedicated to particular regularized contact models, like the approach found in Ref. [33] which is based on point cloud representations in the context of computer animation, or the one explained in Ref. [19] which explicitly computes polyhedral intersections. Others are especially designed for the use of mechanical solvers which are not based on time-stepping algorithms, like “continuous collision detection” as defined in Ref. [36]. The latter approach is closely related to the determination of the instant of first contact, as exploited in event-driven methods (see again Baraff’s excellent descriptions in Refs. [3] and [4]).

However we retain that collision detection methods can be considered as belonging to the field of computational geometry and may show a relatively high level of algorithmic complexity. Besides, when involving polyhedral shapes, these methods very often demand strong conformity hypotheses like 2-manifold triangulations of shape boundaries (see section 4.3).

## 3.2 Necessity of discrete shape representation

Another recurrent aspect of discussions about collision detection methods is the necessity of a tradeoff between efficiency and the quality of the computed geometrical information. As body shapes defined by CAD models are mostly made of piecewise smooth surfaces and curves, some methods, seeking geometric accuracy, have been proposed in order to compute distances or penetration depths between curved surfaces (see Refs. [22], [24], [39] and [15]). All of them exhibit computation times that currently make them incompatible with interactive multibody contact simulations with complex shapes, since they mostly rely on iterative optimization algorithms with possibly slow convergence rates or involving high-degree polynomial root-finding subroutines. In contrast, methods based on approximate discrete shape representations can benefit from the use of bounding volume hierarchies of diverse types (see Ref. [14] and the surveys cited in the preceding subsection) to attain tractable efficiency. The control of geometric approximation will be discussed in section 4.

## 3.3 Time-stepping compatibility issues

In section 2, we have already pointed out that, when relying on set of point contacts for defining the admissible domain $C(t)$, smooth gap functions taking both positive, zero and negative
values have to be geometrically defined. In their practical resolution, time-stepping methods as those presented in Refs. [41], [2], [6], [42] and [38] effectively need gap functions to be able to take both positive and – at least small – negative values.

The first element that may interfere with the convergence of these methods is the introduction, that we already mentioned in subsections 2.4 and 2.5, of hysteresis in some definitions of negative gaps. Another issue comes from the fact that in some methods (e.g. those based on penetration depths only), numerical access to positive gaps is possible only after exploring a penetration configuration. In these two cases, dangerous “accumulation of constraints” heuristics may be necessary, for example during iterations of Josephy-Newton or other Newton-like or Gauss-Seidel-like algorithms used to solve nonlinear complementarity problems or discretized differential inclusions, leading to robustness issues in regions where $C(t)$ is not convex. Note that the use of popular fast time-stepping solvers, relying on a single execution of the collision detection method and linearizing contact equations in order to solve only one LCP per time step, may then give unstable numerical behaviors and/or large interpenetrations, unless the contact unilateral constraints are heuristically accumulated from the preceding time steps.

The second class of frequently encountered issues comes from the fact that the computed geometric information may not be adapted to the use of a non-smooth contact model like Signorini or Coulomb, because of a lack a regularity in the computed contact kinematics. If those are known to be non-smooth for modeling reasons, then non-smooth numerical methods that are able to treat those non-regularities have to be invoked. Otherwise, if the smoothness of the contact kinematics is an hypothesis which is exploited by the time-stepping method, then its numerical violation may result in unpredictable numerical behavior.

# 4 SIMPLICIAL COMPLEXES AS VERSATILE DISCRETE SHAPE REPRESENTATIONS

In the preceding section, we have seen that efficient collision detection methods applying to objects with complex shapes currently make mandatory the use of discrete representations. Moreover, working with industrial CAD models implies being able to represent shapes of arbitrary topology and to control all geometric approximations. We thus propose a unified representation of shapes using simplicial complexes generated by meshing tools.

Roughly speaking, a 3D simplicial complex can be seen as generalized triangulation made of simplices of $\mathbb{R}^3$. In this section the reader will first find a rigorous definition of a simplicial complex and may get familiarized with the associated terminology and the notations we will use in the rest of the paper. Secondly, we motivate and explain the use of simplicial complexes as rigid body shape representations, especially in the case of industrial CAD models, and give guidelines about their generation by meshing tools.

## 4.1 A short introduction to simplicial complexes

Simplicial complexes generalize the notion of triangulation to domains of arbitrary topological dimension (see for example Fig. 3 at the right hand side). Particular cases include tetrahedral meshes, triangulated surfaces, polylines and sets of isolated points.

Formally (see for example Ref. [31]), we recall that a simplex (or $k$-simplex) $S$ in $\mathbb{R}^d$ is defined as the convex hull of $k + 1$ affinely independent points $\{p_0, ..., p_k\}$. The simplices defined as the convex hulls of subsets of $\{p_0, ..., p_k\}$ are called the *faces* of $S$.

A simplicial complex $K$ in $\mathbb{R}^d$ is a collection of simplices, called its *faces*, such that:

* every face of a simplex of $K$ is in $K$;
• the intersection of any two simplices of $K$ is a face of each of them.

In the terminology of algebraic topology, the union of all the faces of a simplicial complex $K$, equipped with the smallest topology which is compatible with the ones of its faces (see Ref. [31] for details), is called the *polytope* of $K$ and denoted $|K|$. A topological space which the polytope of a finite simplicial complex will be called a *polyhedron*.

We also give a few more definitions that will be useful in the next section. Let $S$ be a simplex in $\mathbb{R}^d$. The faces of $S$ which are different from $S$ itself are the *proper faces* of $S$. Their union is called the *simplicial boundary* of $S$ and denoted $\text{Bd } S$. The *simplicial interior* of $S$ is defined as $S \setminus \text{Bd } S$. We also recall that the *relative interior* (resp. the *relative boundary*) of a subset $A$ of $\mathbb{R}^d$ if the interior (resp. the boundary) of $A$ in the natural topology (i.e. the one which comes from the Euclidean structure) of the affine subspace of $\mathbb{R}^d$ generated by $A$. The relative interior (resp. boundary) of a simplex coincides with its simplicial interior (resp. with its simplicial boundary).

<figure><img src="image_4.png" /><figcaption>Figure 3: Left: a simplicial complex. Right: A point contact between a plane and adilated shell.</figcaption></figure>

From a numerical point of view, simplicial complexes appear as generalized meshes, built upon 1-simplices of the complexes. The points in $\mathbb{R}^d$ which belong to a 1-simplex of a complex $S$ are called the *vertices* of $S$, but we will also use the term *vertex* for a 1-simplex itself. Hence a natural way to represent a simplicial complex in a data structure is to first give a list of its vertices as numbered points in $\mathbb{R}^d$, and then represent each of its $k$-simplices by a list of $k + 1$ vertex numbers. A necessary conformity condition is that each for each represented $k$-simplex, the $k + 1$ listed points are affinely independent. Otherwise, we will speak of a *degenerate* simplex. Examples in 3D include flat triangles and tetrahedra and zero-length edges. A complex containing such simplices will also be called degenerate.

## 4.2 Working with complex shapes from industrial CAD models

Another reason for the use of discrete shape representations instead of continuous models is an industrial motivation, since discrete representations are not dependent on a particular CAD kernel native format. As industrial constraints impose a complete control geometric errors, including during collision detection and contact simulation, one could think of reducing the total errors by directly working on native formats for collision detection. Unfortunately, not all of these formats are open, and the transfer formats containing continuous representations and made available in industrial CAD software products, like IGES for instance, are also based on approximations of the native models.

We actually propose to use *dilated*<fnref n="1" /> simplicial complexes for representing *geometric shapes* of arbitrary topological dimension, as in Fig. 3. This is coherent with our preliminary assumptions

<footnote n="1">¹The $\varepsilon$-dilation of a subset $A$ in $\mathbb{R}^3$ is an operation defined for any non-negative $\varepsilon$ value which consists in</footnote>
tion stating that the *body shapes* should be compact volumes (see section 2), even if they are
obtained by dilating geometric shapes of lower (and even not uniform) topological dimension.

We propose to represent two-dimensional models like thin shells by triangulations of the medial surfaces dilated by half the thicknesses of the shells. One-dimensional models like thin beams or cables will be preferably represented by dilated polylines (chained straight line segments) coming from the discretization of their medial fiber. For volumes, two options exists: use a complete tetrahedral triangulation of the volume, or simply an oriented 2-manifold triangulation of the boundary. We will see in the next section that our method has the drawback of requesting a positive dilation of complexes. To satisfy this constraint, the use of a small negative offset before or during the boundary surface meshing of volume models can compensate the necessary dilation.

4.3 Meshing recommendations, mesh flaws and geometric error control

The geometric approximations made during the meshing operations are generally controlled
via error criteria, like the maximum chordal deviation. Other meshing constraints may imposed
in order satisfy geometric conformity criteria. Industrially available tools for generating (gen-
eralized) triangulations may be classified in two main categories, which are *tessellators* and
*meshers*.

On one hand, tessellators are generally first aimed at generating discretized representation
of lines (or “polylines”) and triangulations of surfaces (including volume boundaries) for vi-
sualization purposes. Therefore they must be capable of treating every native model of almost
arbitrary quality with an extremely reduced reject rate, but the conformity criteria they enforce
will be very modest, focalizing on the optimization of the number of generated triangles. Be-
sides industrial tessellators exhibit fairly diverse output quality. However, high-quality tessella-
tors are able to generate 2-manifold triangulations of “clean” CAD surface models and volume
boundaries, with a control over the maximum edge length and sometimes over the maximum
dihedral angle between neighbor triangles. Such tessellators are currently available in a strong
proportion of industrial CAD software suites.

On the other hand, meshers are expected to give the possibility to control more various
geometric criteria for the needs of numerical methods like finite elements. A volume tetrahedral
mesher is also the only solution for the triangulation of the interior of volumes. Unfortunately,
meshers are known to demand input CAD models of relatively high quality, and sometimes
manual corrections on their output, which may represent large amounts of pre-processing time.

In industrial contexts where rigid multi-body simulations would be useful, large pre-processing delays and costs are not always acceptable. In the early conception stages of a product, they could even destroy the potential advantages of digital mock-ups over traditional physical ones. In those cases, a high-quality tessellator may be the tool of choice. Still, if the input CAD models are of poor quality and suffer flaws like cracks, holes, or more serious topological aberrations, the output tessellation is likely to exhibit comparable defects. Therefore multi-body simulators relying on collision detection algorithms which absolutely require flawless input triangulations are practically unusable in these contexts.

Our method only requests (see section 6) that the input simplicial complexes are non-degenerate, which is most of the time attained by the use a high-quality tessellator followed by very basic mesh untangling algorithms when degenerate simplices are detected. If cracks, holes, or

<footnote>performing the union all closed balls of radius $\varepsilon$ having their centers in $A$. The resulting set is the Minkowski sum of $A$ and the closed ball of radius $\varepsilon$ centered at the origin.</footnote>
non-manifold edges are detected in the triangulation of a volume boundary, then we can choose to treat it as a shell from the point of view of collision detection.

We conclude this section by signaling that the control of geometric errors induced by the successive negative offset (for volumes only), meshing, and dilation operations can be completely controlled, provided that the tools used for the two first steps enforce at least a maximum chordal deviation criterion.

# 5 TIME-STEPPING COMPATIBLE COLLISION DETECTION BETWEEN DILATED SIMPLICIAL COMPLEXES

We follow our claim of subsection 2.3, and we propose to define point contacts between dilated simplicial complexes thanks to LMDs, using the dilation for defining negative gaps. We will define actually the slightly weaker notion of quasi-LMDs which enables the stable treatment of conforming contact cases.

## 5.1 A retraction method for defining negative gaps and smooth contact kinematics

Consider two dilated simplicial complexes and suppose that the non-dilated complexes are disjoint. If the shapes obtained by dilation are disjoint too, we remark that each LMD between them can be obtained by “shortening” a LMD between the complexes as shown in Fig. 4 at the center. Moreover, the non-penetration of the dilated simplicial complexes is equivalent to the condition that all the LMDs between the base complexes give distances greater than the sum of the dilations.

Starting from these central remarks, our approach consists in using the LMDs between the non-dilated simplicial complexes for defining point contacts between the dilated complexes. Let $(a_i, a_j)$ be a LMD between the disjoint simplicial complexes $K_i$ and $K_j$. For each such LMD, we define a point contact between the shapes obtained by $\varepsilon_i$-dilation of $K_i$ and $\varepsilon_j$-dilation $K_j$ by deciding that the couple of points defining the generalized point contact kinematics is $(a_i + \varepsilon_i n_i, a_j + \varepsilon_j n_j)$, where $n_i$ and $n_j$ are defined by:

$$n_i = -n_j = \frac{a_j - a_i}{\|a_j - a_i\|}.$$

$n_i$ defines the normal direction to contact.

Maintaining the complexes separated during simulation is then ensured by the use of non-smooth contact laws imposing non-negative gaps, although small positive and negative gap values are accessible for internal time-stepping computations<fnref n="2" />. Note that the resulting generalized contact kinematics are smooth, hence enabling the introduction of frictional contact laws. Moreover all computations are reported on the non-dilated simplicial complexes. Consequently this purely implicit “retraction method” appears like an advantageous alternative to local penetration depths or any other difficult construction mentioned in subsection 2.4. It makes the definition of generalized contact kinematics as easy as if the shapes were simply two balls, and it is not surprising since actually they are unions (a potentially infinite number) of balls.

Besides, our method is purely geometric and does not invoke the history of relative movement between the shapes. Thanks the positive dilation, it also avoids the well-known consistency issues in defining contacts between pure surfaces or thin beams or cables. A comparable

<footnote n="2">²Incremental control of displacements, resembling incomplete descent strategies, can be used to prevent the possible intermediary computation of time-stepping solvers to explore configuration where the base complexes interpenetrate.</footnote>
dilation method applied to convex polyhedra is described in Ref. [7], but we have not found any
reference which exploits it in the case of non-convex complex shapes.

<figure><img src="image_5.png" /><figcaption>Figure 4: Left: condition (2) in configuration manifold $\mathcal{M}$. Center: gaps between dilated simplicial complexes. Right: a LMD and a quasi-LMD (at the right hand side in green).</figcaption></figure>

## 5.2 Quasi-LMDs : stable and generic treatment of conforming contact cases

First remark that conforming contact regions, as defined in subsection 2.5, between two
shapes obtained by the dilation of simplicial complexes are planar polyhedra. To treat those
cases, we would like our method to reproduce a contact selection strategy, but without relying
on history.

Now let us examine the simple example of a dilated cube resting a plane, the contact region
being a square. Since it exists no strict local minimum of the distance function restricted the
cartesian product of the shapes, all LMDs vanish in this relative configuration. Ideally we would
like to define point contacts based on the four contacting corners of the cube, by it is easy to see
that, for any relative pose, at most two LMDs exist between the cube and the plane.

If for defining LMDs we accept all minima of the distance function, then an infinity of them
can be defined for a conforming contact region. Further restricting those minima to the corners
of the cube is only a partial solution, since at least two of the point contacts defined in this
manner vanish under infinitely small rotations of the cube, impeding convergence of most time-
stepping scheme and giving rise to unstable numerical behaviors.

Our approach is actually based on the definition of *quasi-LMDs*, which corresponds to a
prolongation of the definition of particular LMDs under small rotations (see Fig. 4 at the right
hand side). If the LMDs supported by the corners of the cube receive such a regularization, then
four point contacts are defined for every sufficiently small perturbation of the exact conforming
contact situation. This regularization of constraints with respect to small rotations yields a stable
treatment of conforming contact cases with a typical time-stepping scheme, for example a fast
solver based on a single LCP, as will be demonstrated in section 7.

The general definition of quasi-LMDs necessitates the introduction of supplementary mate-
rial and will be done in section 6, but the following subsection gives an immediate analysis of
the regularization of constraints giving by LMD gaps with respect to small rotations.

## 5.3 Admissible domain boundary covering conditions

With non-convex simplicial complexes, constraint functions based on LMD gaps are defined
only locally in the configuration manifold $\mathcal{M}$, having the convention that the constraint defined
by a gap function is simply ignored outside of its definition domain. Let us denote by $D_i(t)$
the domain where the $i$-th constraint function $g_i(\cdot, t)$ is defined and continuously differentiable,
by $\mathcal{V}(X)$ the set of all neighborhoods of any set $X \subset \mathbb{R}^n$. $C(t)$ is the domain where all gap functions are non-negative. The prolongation of LMDs into quasi-LMDs actually corresponds the enlarging the $D_i(t)$ domains without modifying $C(t)$.

Actually typical time-stepping schemes need the constraints functions $g_i(\cdot, t)$ to be defined on a whole vicinity of the portion of its zero level which belongs to the boundary of $C(t)$, which can be formalized by the following qualification condition (see Fig. 4 at the left hand side):

$$ \forall t, \forall i \in \{1, ..., p\}, \exists \Omega_i \in \mathcal{V}((g_i(\cdot, t))^{-1}(\{0\}) \cap \partial C(t)). \Omega_i \subset D_i(t). \quad (2) $$

Otherwise, “flip-flop” instabilities may occur in the corners of $\partial C(t)$, the constraint solver used by time stepping scheme retaining only a subset of the constraints forming the corner, either during its successive iterations in the same time step or from one time step to another. Indeed condition (2) is equivalent to the following one:

$$ \forall t, \forall x \in \partial C(t), \exists \Omega_x \in \mathcal{V}(\{x\}), \forall i \in \{1, ..., p\}, g_i(x, t) = 0 \Rightarrow \Omega_x \subset D_i(t), $$

which means that all constraints which are active at a corner of $C(t)$ (corresponding for example to a conforming contact situation) should still be defined under sufficiently small displacements in $\mathcal{M}$. Such small displacements may cause both relative translations and rotations of the contacting rigid bodies, which explains the stabilizing effect of quasi-LMDs.

# 6 EFFICIENT QUASI-LMD COMPUTATIONS BASED ON AUGMENTED BOUNDING SPHERE-CONE HIERARCHIES

Our method computes LMDs and quasi-LMDs efficiently thanks to a bounding volume hierarchy (BVH) inspired from Johnson and Cohen's (see Ref. [23] and Ref. [25]). However the definition of normal cones given in Ref. [23] is derived from normals to volumes with smooth boundaries. The adequate tools for giving a rigorous characterization of a LMD between volumes with polyhedral boundaries are actually provided by convex analysis and apply to general simplicial complexes. The cones that have to be considered are the polar cones of the tangent cones to the complexes, which are slightly different from the normal cones.

## 6.1 Tangent and normal cones, Gauss map

We follow the terminology found for example in Ref. [8], and adopt the following definitions for the tangent and normal cones.

Let $X$ be a subset of $\mathbb{R}^3$ and $x$ a vector in $X$. A vector $y \in \mathbb{R}^3$ is said to be a tangent of $X$ at $x$ if either $y = 0$ or there exists a sequence $(x_k)_{k \in \mathbb{N}}$ such that $x_k \in X \setminus x$ for all $k$, $\lim_{k \to \infty} x_k = x$ and $\lim_{k \to \infty} (x_k - x)/\|x_k - x\| = y/\|y\|$. The set of all tangents of $X$ at $x$ is called the tangent cone of $X$ at $x$, and is denoted $T_X(x)$. A vector $z \in \mathbb{R}^3$ is said to be a normal of $X$ at $x$ is there exists sequences $(x_k)_{k \in \mathbb{N}}$ and $(z_k)_{k \in \mathbb{N}}$, with respective limits $x$ and $z$ and such that for all $k$, $x_k \in X$ and $z_k \in T_X(x_k)^*$. The set of all normals of $X$ at $x$ is called the normal cone of $X$ at $x$ and is denoted $N_X(x)$. If $N_X(x) = T_X(x)^*$, $X$ is said to be regular, or tangentially regular, at $x$. We also define the Gauss map of $X$, denoted $G_X$, to be the set-valued function defined by $G_X(x) = T_X(x_k)^* \cap S^2$, where $S^2$ is the unit sphere in $\mathbb{R}^3$.

Now, let $S$ be a $d$-simplex in $\mathbb{R}^3$ with $d > 0$, $\sigma$ a face of $S$ of dimension $d-1$ and $v$ the vertex of $S$ opposite to $\sigma$. Let $a$ be the orthogonal projection of $v$ on the affine hull on of $\sigma$. The points $a$ and $v$ are distinct. The unit vector $(v - a)/\|v - a\|$ will be denoted $u_S(\sigma)$.

Now let $K$ be a finite simplicial complex in $\mathbb{R}^3$. The relative interiors of the faces of $K$ form a partition of $|K|$. It is not difficult to show that, for each face $\sigma$ of $K$, the function $x \mapsto T_{|K|}(x)^*$
takes a constant value over Int $\sigma$, which will be denoted $T_{|K|}(\text{Int } \sigma)^*$, and is described by:

$$T_{|K|}(\text{Int } \sigma)^* = (\text{Dir Aff } \sigma)^\perp \cap U_K(\sigma)^*, \quad (3)$$

where Dir Aff $\sigma$ denoted the direction of the affine hull of $\sigma$ and $U_K(\sigma)$ denotes the set $\{u_{\sigma'}(\sigma), \sigma' \in K, \sigma \subset \sigma', \dim(\sigma') = \dim(\sigma) + 1\}$.

Let $x$ be a point in $|K|$, and denote $\sigma_x$ be the unique face $K$ which contains $x$ in its relative interior. We have:

$$T_{|K|}(x)^* = T_{|K|}(\text{Int } \sigma_x)^*. \quad (4)$$

## 6.2 Characterization of LMDs between simplicial complexes

Let $K_1$ and $K_2$ be two finite and disjoint simplicial complexes in $\mathbb{R}^3$, and denote $d_{|K_1| \times |K_2|}$ the restriction of the Euclidean distance to $|K_1| \times |K_2|$. Actually we will characterize all local minima, strict or not, of the function $d_{|K_1| \times |K_2|}$. A couple $(x_1, x_2) \in |K_1| \times |K_2|$ is a local minimizing point of $d_{|K_1| \times |K_2|}$ if and only if the following condition holds:

$$x_2 - x_1 \in T_{|K_1|}(x_1)^* \cap -T_{|K_2|}(x_2)^*. \quad (5)$$

Here is an outline of the proof: assume that (5) holds. Thanks to (3) and (4), we can show that for each $i$, there exists $\varepsilon_i > 0$ such that $T_{|K_i|}(x_i)^* = T_{\nu_i}(x_i)^*$, where $\nu_i$ denotes the convex hull of $|K_i| \cap B(x_i, \varepsilon_i)$. Besides $T_{|K_1| \times |K_2|}(x_1, x_2)^*$ is equal to $T_{|K_1|}(x_1)^* \times T_{|K_2|}(x_2)^*$ and (5) implies that $-\nabla d_{|K_1| \times |K_2|}$ is in $T_{|K_1| \times |K_2|}(x_1, x_2)^*$. Hence $-\nabla d_{|K_1| \times |K_2|}$ is in $T_{\nu_1 \times \nu_2}(x_1, x_2)^*$, and since $d_{|K_1| \times |K_2|}$ and $\nu_1 \times \nu_2$ are convex, $(x_1, x_2)$ locally minimizes $d_{|K_1| \times |K_2|}$ in $\nu_1 \times \nu_2 \subset |K_1| \times |K_2|$. The converse implication is straightforward.

Now let $X_1$ and $X_2$ be subsets of $|K_1|$ and $|K_2|$ respectively. A necessary condition for the existence of a couple $(x_1, x_2) \in X_1 \times X_1$ which locally minimizes $d_{|K_1| \times |K_2|}$ is that the following condition holds:

$$\bigcup_{x \in X_1} T_{|K_1|}(x_1)^* \cap \bigcup_{x \in X_2} -T_{|K_2|}(x_2)^* \cap X_2 - X_1 \neq \emptyset. \quad (6)$$

## 6.3 Definition of quasi-LMDs

Let $K_1$ and $K_2$ be two finite and disjoint simplicial complexes in $\mathbb{R}^3$. The definition quasi-LMDs comes with a special classification of the faces $\sigma$ of each complex $K_i$, depending on the topological dimensions of $\sigma$ and $T_{|K|}(\text{Int } \sigma)^*$, which is summarized in table 1. We define eight classes numbered from 0 to 7, and the class of $\sigma$ in $K$ is denoted $c_K(\sigma)$. Table 2 then defines a symmetric binary compatibility relation between the classes as well as a regularization function $\varepsilon$ which determines which type of regularization of LMDs with respect to orientation is to be made.

Precisely, we define the quasi-LMDs between $K_1$ and $K_2$ as follows. Let $\sigma_1$ and $\sigma_2$ be faces of $K_1$ and $K_2$ respectively. A quasi-LMDs $(x_1, x_2)$ is supported by $\sigma_1$ and $\sigma_2$ if the following conditions hold:

* The classes of $\sigma_1$ and $\sigma_2$ are compatible, i.e. the value of $\varepsilon(c_{K_1}(\sigma_1), c_{K_2}(\sigma_2))$ is defined.

* $d_{\text{Aff } \sigma_1 \times \text{Aff } \sigma_2}$ has a unique minimizer which is $(x_1, x_2)$.

* $(x_1, x_2)$ is in $\text{Int } \sigma_1 \times \text{Int } \sigma_2$.

* For each $i \{1, 2\}$, for all $\sigma' \in U_{K_i}(\sigma_i)$, $\langle x_j - x_i / \|x_j - x_i\|, u_{\sigma'}(\sigma) \rangle < \sin(\varepsilon_i(c_{K_1}(\sigma_1), c_{K_2}(\sigma_2)))$, with $j \in \{1, 2\} \setminus i$.
<table><caption>Table 1: Classification of the simplices $\sigma$ of the complex $K$. A dot indicates all possible values.</caption><thead><tr><th>dim(σ)</th><th>dim (T<sub>|K|</sub>(Int σ)<sup>*</sup>)</th><th>St<sub>K</sub> σ ∩ ∂|K|] is a plane</th><th>c<sub>K</sub>(σ)</th></tr></thead><tbody><tr><td>0</td><td>3</td><td>.</td><td>0</td></tr><tr><td>0</td><td>2</td><td>.</td><td>1</td></tr><tr><td>1</td><td>2</td><td>.</td><td>2</td></tr><tr><td>0</td><td>1</td><td>false</td><td>3</td></tr><tr><td>0</td><td>1</td><td>true</td><td>4</td></tr><tr><td>1</td><td>1</td><td>.</td><td>5</td></tr><tr><td>2</td><td>1</td><td>.</td><td>6</td></tr><tr><td>.</td><td>0</td><td>.</td><td>7</td></tr></tbody></table>

<table><caption>Table 2: Definition of the values of the regularization function $\varepsilon$, the scalar $\theta$ being a fixed small positive angle which is a parameter of our method. “n.d.” indicates an undefined value, i.e. incompatibility between the classes.</caption><thead><tr><th>ε(c<sub>1</sub>, c<sub>2</sub>)</th><th>0</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><th>0</th><td>(0, 0)</td><td>(0, 0)</td><td>(θ, 0)</td><td>(θ, 0)</td><td>(θ, 0)</td><td>(θ, 0)</td><td>(θ, 0)</td><td>n.d.</td></tr><tr><th>1</th><td>(0, 0)</td><td>(θ, θ)</td><td>(θ, θ)</td><td>(θ, 0)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>2</th><td>(0, θ)</td><td>(θ, θ)</td><td>(θ, θ)</td><td>(θ, 0)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>3</th><td>(0, θ)</td><td>(0, θ)</td><td>(0, θ)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>4</th><td>(0, θ)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>5</th><td>(0, θ)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>6</th><td>(0, θ)</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr><tr><th>7</th><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td><td>n.d.</td></tr></tbody></table>

* For each $i$, let $n_{K_i}(\sigma_i) \subset S^2$ be a minimal generator of the polyhedral cone $U_{K_i}(\sigma_i)^*$, and let $C_{K_i}(\sigma_i)$ be the smallest semi-infinite revolution cone, with axis co-linear with the mean of the vectors in $n_{K_i}(\sigma_i)$, which contains $U_{K_i}(\sigma_i)^*$. The vector $x_j - x_i$ should be in the revolution cone obtained by augmenting of $\varepsilon_i(c_{K_1}(\sigma_1), c_{K_2}(\sigma_2))$ the half angle of $C_{K_i}(\sigma_i)$.

## 6.4 Bounding volume hierarchies for efficient quasi-LMD computations

We propose to use a sphere-cone bounding volume hierarchy (BVH) similar to introduce in Ref. [23], and generalize it by deciding that the leafs of the hierarchy contains the simplicial interiors of the faces of a simplicial complex $K$, which we recall to be a partition of $||K||$. We also augment this 5D BVH to a 5D+1 by addition the discrete dimension of the simplex classes, i.e. we add eight boolean value to each sphere-cone bounding volume indicating whether the node of the hierarchy contains simplices of each class. Our node-node test is based on the compatibility relation defined by table 2 on the necessary condition (6) just like in Ref. [23], but with the half-angles of the node cones augmented of the regularization parameter $\theta$. The leaf-leaf test immediately follows the preceding subsection.

Moreover, for the top-down creation of the hierarchy, we propose to adopt a dimensional splitting strategy. Recursive subdivision is done by alternating between relying on the position and relying on the value of the Gauss map of the primitives (we perform a median cut along the first principal axis). For the computation of sphere bounding volumes, we use a layered hierarchy and compute the minimal enclosing ball of the simplices (see Ref. [12]). For the bounding cones, we adopt a nested hierarchy, starting from the leafs and using the cones $C_K(\sigma)$.
For the special case of volumes, the quasi-LMDs are located on the boundaries which are two-dimensional *oriented* simplicial complexes. In this particular case, our method, like Johnson's, can take as input this 2D complex and rely on its orientation instead of using a triangulation of the whole volume (which would also contains the filling tetrahedrons).

# 7 EFFICIENCY AND ROBUSTNESS BENCHMARKS ON INDUSTRIAL DIGITAL MOCK-UP CASES

Our method has been integrated in an interactive multibody simulation engine developed at CEA and based on a complementarity formulation of non-smooth unilateral contact solved by a semi-implicit time-stepping scheme based on Ref. [20]. This section shows experimental results demonstrating the robustness and efficiency of our approach.

## 7.1 Complex conforming contacts

Tests of robustness and time-stepping compatibility have been performed on various models ranging from the simple cube example introduced in section 5.2 (see Fig. 5) to complex shapes coming from industrial CAD models (see Fig. 1 and Fig. 5). Stable conforming contacts are obtained thanks to quasi-LMDs, even in geometrically complex situations.

<figure><img src="image_6.png" /><figcaption>Figure 5: Left: exact conforming contact between a cube and plane, showing four quasi-LMDs. Center: quasi-LMDs do not vanish under small rotations of the cube. Right: a complex conforming contact situation.</figcaption></figure>

## 7.2 Efficiency and robustness evaluations on industrial models

Industrial benchmarks featuring volumes and shells (see Fig. 1 and Fig. 6) have been used for the quantitative evaluation of efficiency and robustness of our method. The evaluation has been done in the case of a fast solver, with explicit collision detection (only one quasi-LMD query per time step, LCP formulation of perfect unilateral constraints). Table 3 summarizes statistics made on whole trajectory which corresponds to the mounting of a winds shield wiper's actuating mechanism (rigidified during the mounting phase). All other mechanical parts are rigid bodies which are fixed to the global inertial frame. The trajectory starts outside of the car structure, goes through a narrow passage between the base of the windshield and the superior edge of the hood, and finishes with a peg-in-hole insertion into the supporting parts (see Fig. 6). This is real scenario which was proposed to us by PSA-Peugeot-Citroën.

Several meshing parameters have used for the generation of the triangulations and shells and volume boundaries, which were obtained using CATIA V5 tessellator. The maximum chordal deviation *s* was chosen between 1 and 2 millimeters, which is coherent with level of geometrical accuracy demanded in the scenario. The maximum edge length *l* was varied between 3 and 20 millimeters, and in table 3 we have retained only two combinations of these parameters which are (*s*, *l*) = (2, 10) and (*s*, *l*) = (1, 5). The corresponding tessellations had respectively 521, 965
<figure><img src="image_7.png" /><figcaption>Figure 6: Left: a highly constrained industrial part mounting scenario. Right: a peg-in-hole situation.</figcaption></figure>

and 1, 357, 649 triangles for the whole scene, which lead to a total of respectively 1, 577, 453 and 4, 091, 987 simplices, for a total number of 26 rigid body shapes.

The performance measurements were done a Hewlett-Packard Proliant DL 145 G2 server, equipped with two dual-core Opteron 275 processors. Our implementation is multithreaded at the level of the 25 pairs of shapes on which quasi-LMD computations are requested. The query times performed by our method have enabled the interactive simulation of the whole mounting process using a space-mouse for 6-degree-of-freedom manipulation of the moving body and visual feedback about the contact constraints indicated by colored arrows for each quasi-LMD.

With other models of slightly more modest complexity, real-time simulation rates were obtained, enabling high-quality haptic interaction. Those results will be presented in a future communication.

<table><caption>Table 3: First column: 5D BVH, no use of orientation in the recursive subdivisions when building the BVH, which is very close to Johnson's original method found in Ref. [23]. Other columns: our full 5D+1 BVH method.</caption><thead><tr><th>Time per request (ms)</th><th>(s, l) = (2, 10), “no opt.”</th><th>(s, l) = (2, 10)</th><th>(s, l) = (1, 5)</th></tr></thead><tbody><tr><th>Average</th><td>15.42</td><td>14.00</td><td>10.06</td></tr><tr><th>Min</th><td>0.16</td><td>0.30</td><td>0.41</td></tr><tr><th>Max</th><td>142.58</td><td>45.04</td><td>29.18</td></tr></tbody></table>

## 7.3 “Less simplices implies faster queries” is false

Efficiency of bounding volume hierarchies is known to greatly depend on the properties of the discretized models, classical approaches consisting in minimizing the number of triangles in surface meshes. We conclude this section by signaling that experiments show (see third column of table 3) that optimal efficiency of our method is attained for a certain size of simplex, controlled here by imposing by a maximum edge length during the mesh generation. Inhomogeneous triangle sizes, as found in visual tessellations, are also to be avoided.

## 8 CONCLUSION AND FUTURE WORK

Our collision detection method enables time-stepping schemes to be used robustly for non-smooth contact between rigid bodies of arbitrary topology and geometric complexity, including industrial CAD models, with a high level of efficiency, an all-time control of geometric approximations, and a stable and generic treatment of conforming contact situations. The computed contact kinematics are sufficiently smooth and complete for classical non-smooth contact models to be applied.

We see no obstacle in using our method in the context of event-driven schemes, or with some compliant or regularized point contact models. In addition, a highly optimized parallel
implementation for distributed memory architectures is in progress, as well as an adaptation to
the case of deformable bodies. Future work also includes quasi-LMD computations between
shapes represented by medial axis transforms in the rigid case, remarking that those shapes are
also defined by the union of an infinite number of balls, but with variable radius.

9 ACKNOWLEDGEMENTS

We thank PSA-Peugeot-Citroën, Renault and Dassault Aviation for providing industrial CAD
models and scenarios for our software benchmarks.

REFERENCES

[1] P. Agarwal, L. J. Guibas, S. Har-Peled, A. Rabinovitch, and M. Sharir. Penetration depth of two convex polytopes in 3D. *Nordic J. Computing*, **7**, 227–240, 2000.

[2] M. Anitescu, F.A. Potra, and D.E. Stewart. Time-stepping for the three dimensional rigid body dynamics. *Comp. Meth. Appl. Mech. Eng.*, **177**, 183–197, 1999.

[3] D. Baraff. Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies. In Proc. *ACM Siggraph '89*, published as *Computer Graphics*, **23:3**, 223–232, 1989.

[4] D. Baraff. Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation. In Proc. *ACM Siggraph '90*, published as *Computer Graphics*, **24:4**, 19–28, 1990.

[5] F. Barbagli, A. Frisoli, K. Salisbury, M. Bergamasco. Simulating human fingers: a Soft Finger Proxy Model and Algorithm. In Proc. *Int. Symp. on Haptic Interfaces for Virtual Environment and Teleoperator Systems*, 9–17, 2004.

[6] B. Brogliato. *Nonsmooth Mechanics*. Springer Verlag, Second Edition, 1999.

[7] G. van den Bergen. Proximity queries and penetration depth computation on 3d game objects. In Proc. *Game Developers Conference*, 2001.

[8] D. P. Bertsekas. *Convex Analysis and Optimization*. Athena Scientific, 2003.

[9] V. Duindam and S. Stramigioli. Modeling the kinematics and dynamics of compliant contact. In Proc. *2003 IEEE Int. Conf. on Robotics and Automation*, **3**, 4029–4034, 2003.

[10] S. A. Ehmann and M. C. Lin. Accurate and Fast Proximity Queries Between Polyhedra Using Convex Surface Decomposition. In Proc. *Eurographics'2001*, published as *Computer Graphics Forum*, **20:3**, 500–510, 2001.

[11] C. Ericson. *Real Time Collision Detection*. Morgan Kaufmann Publishers, Elsevier, 2005.

[12] K. Fischer and B. Gärtner. The Smallest Enclosing Ball of Balls: Combinatorial Structure and Algorithms. *Int. J. of Computational Geometry and Applications*, **14**, 341–378, 2004.

[13] C. Glocker. Impacts with Global Dissipation Index at Reentrant Corners. In Proc. *3rd Contact Mechanics Int. Symp.*, 45–52, 2001.

[14] S. Gottschalk, M. Lin, and D. Manocha. OBB-Tree: A hierarchical structure for rapid interference detection. In Proc. *ACM Siggraph'96*, published as *Computer Graphics*, **30**, 71–180, 1996.
[15] S. Hartmann and M. Hiller and A. Jennert and M. Sabinarz. CAD-Based Contact Mod-
elling For Multibody Systems. In Proc. *IUTAM Symposium on Unilateral Multibody Dy-
namics*, 1998.

[16] S. Hasegawa and M. Sato. Real-time Rigid Body Simulation for Haptic Interactions Based on Contact Volume of Polygonal Objects. In Proc. *Eurographics'2004*, published as *Computer Graphics Forum*, **23:3**, 529-538, 2004.

[17] B. Heidelberger, M. Teschner, R. Keiser, M. Mueller and M. Gross. Consistent Penetration Depth Estimation for Deformable Collision Response. In Proc. *Vision, Modeling, Visualization VMV'04*, 339-346, 2004.

[18] B. Heidelberger, M. Teschner and M. Gross. Real-Time Volumetric Intersections of De-
forming Objects. In Proc. *Vision, Modeling, Visualization VMV'03*, 461–468, 2003.

[19] G. Hippmann. An Algorithm for Compliant Contact Between Complexly Shaped Surfaces in Multibody Dynamics. In Proc. Multibody Dynamics 2003, ECCOMAS Thematic Conference, 2003.

[20] M. Jean. The Non-Smooth Contact Dynamics Method. *Comp. Meth. Appl. Mech. Eng.*, **177**, 235–257, 1999.

[21] P. Jiménez, F. Thomas, and C. Torras. 3D collision detection: A survey. *Computers and Graphics*, **25:2**, 269–285, 2001.

[22] D. E. Johnson and E. Cohen. A Framework for Efficient Minimum Distance Computa-
tions. In Proc. *IEEE Int. Conf. on Robotics and Automation*, 3678–3684, 1998.

[23] D. E. Johnson and E. Cohen. Spatialized Normal Cone Hierarchies. In Proc. *2001 ACM Symposium on Interactive 3D Graphics*, 129–134, 2001.

[24] D. E. Johnson and E. Cohen. Unified Distance Queries in a Heterogeneous Model Envi-
ronment. In Proc. ASME 2004 Design Engineering Technical Conferences, 2004.

[25] D. E. Johnson, P. Willemsen E. and Cohen. 6-DOF Haptic Rendering Using Spatialized Normal Cone Search. *IEEE Trans. on Visualization and Comp. Graphics*, **11:6**, 661–670, 2005.

[26] Y. Kim, M. Otaduy, M. C. Lin and D. Manocha. Six-Degree-of-Freedom Haptic Display Using Localized Contact Computations. In Proc. *IEEE Symp. on Haptic Interfaces for Virtual Environment and Teleoperator Systems*, 209–216, 2002.

[27] J. Klein. *Efficient Collision Detection for Point and Polygon Based Models*. PhD Thesis,
University of Paderborn, Germany, 2005.

[28] P. R. Kraus and V. Kumar. Compliant Contact Models For Rigid Body Collisions. In Proc. *1997 IEEE Int. Conf. on Robotics and Automation*, **2**, 1382–1387, 1997.

[29] M. Lin and S. Gottschalk. Collision detection between geometric models: A survey. In Proc. *IMA Conference on Mathematics of Surfaces*, 33–52, 1998.

[30] D.J. Montana. The Kinematics of Contact and Grasp. *Int. J. Robotics Res.*, **7:3**, 17–32.
[31] J.R. Munkres. *Elements of Algebraic Topology*. Perseus Publishing, 1984.

[32] J. Park, W.-K. Chung and Y. Youm. Second-Order Contact Kinematics for Regular Contacts. In Proc. 2002 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS 2002), 2239-2244, 2002.

[33] M. Pauly, D.K. Pai, L.J. Guibas. Quasi-Rigid Objects in Contact. In Proc. *ACM SIGGRAPH/Eurographics Symposium on Computer Animation*, 2004.

[34] F. Pfeiffer and P. Wolfsteiner. Relative Kinematics of Multibody Contacts. In Proc. *Int. Mechanical Engineering Congress and Exposition*, 1997.

[35] M. K. Ponamgi, D. Manocha and M. C. Lin. Incremental Algorithms for Collision Detection Between Polygonal Models. *IEEE Trans. on Visualization and Comp. Graphics*, **3:1**, 51–64, 1997.

[36] S. Redon, A. Kheddar, and S. Coquillart. Fast continuous collision detection between rigid bodies. In Proc. *Eurographics'02*, 2002.

[37] S. Redon and M.C. Lin. A Fast Method for Local Penetration Depth Computation. *Journal of Graphics Tools*, **11:2**, 37–50, 2006.

[38] M. Renouf, V. Acary and G. Dumont. 3D frictional contact and impact multibody dynamics: A comparison of algorithms suitable for real-time applications. In Proc. Multibody Dynamics 2005, ECCOMAS Thematic Conference, 2005.

[39] K.-A. Sohn, B. Jüttler, M.-S. Kim, and W. Wang. Computing the distance between two surfaces via line geometry. In Proc. *IEEE Pacific Conf. on Comp. Graphics and Applications*, 236-245, 2002.

[40] P. Song and V. Kumar. Distributed Compliant Model for Efficient Dynamic Simulation of Systems with Frictional Contacts. In Proc. 2003 ASME Design Engineering Technical Conferences, 2003.

[41] D.E. Stewart and J.C. Trinkle. An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Inelastic Collisions and Coulomb Friction. *Int. Jour. of Num. Meth. Eng.*, **39**, 2673–2691, 1996.

[42] D.E. Stewart. Rigid-Body Dynamics with Friction and Impact. SIAM Review. **42:1**, 3–39, 2000.

[43] M. Teschner, S. Kimmerle, B. Heidelberger, G. Zachmann, L. Raghupathi, A. Fuhrmann, M.P. Cani, F. Faure, N. Magnenat-Thalmann, W. Strasser, P. Volino. Collision Detection for Deformable Objects. *Computer Graphics Forum*, **24:1**, 2005.

[44] M. Visser, S. Stramigioli, and C. Heemskerk. Screw Bondgraph Contact Dynamics. In Proc. 2002 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems, 2239-2244, 2002.

[45] P. Volino, N. Magnenat-Thalmann. Resolving Surface Collisions through Intersection Contour Minimization. *ACM Transactions on Graphics* (*SIGGRAPH 2006 proceedings*), **25:3**, 1154–1159, 2006.