Xavier Merlhiot
email: xavier.merlhiot@cea.fr

A ROBUST, EFFICIENT AND TIME-STEPPING COMPATIBLE COLLISION DETECTION METHOD FOR NON-SMOOTH CONTACT BETWEEN RIGID BODIES OF ARBITRARY SHAPE

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

This paper proposes an efficient collision detection method which is compatible with time-stepping methods in the sense that it enables the robust simulation 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.

INTRODUCTION

On one hand, time-stepping methods such as those proposed in Refs. [START_REF] Anitescu | Time-stepping for the three dimensional rigid body dynamics[END_REF], [START_REF] Stewart | An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Inelastic Collisions and Coulomb Friction[END_REF] and [START_REF] Jean | The Non-Smooth Contact Dynamics Method[END_REF], 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. [START_REF] Redon | A Fast Method for Local Penetration Depth Computation[END_REF] or [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF]. 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 timestepping 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. 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.

ON THE GEOMETRICAL NEEDS OF THE MODELING OF NON-SMOOTH CON-TACT 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.

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) ⊂ R 3 called its shape, satisfying S i (t) = H i (q(t), t) S i (t 0 ), where H i is an at least twice continuously differentiable function from M × R into SE(3) and t 0 a reference time. We assume that 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 C(t) ⊂ M defined by the non-penetration of the body shapes may receive the following description :

C(t) = {q ∈ M, ∀ (i, j), i = j ⇒ H i (q, t) (S i ) • ∩ H j (q, t) (S j ) • = ∅} .
An effective contact situation between two bodies then corresponds to the existence of points shared between the boundaries of their shapes.

If 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. [START_REF] Brogliato | Nonsmooth Mechanics[END_REF]) that the unilateral constraint forces f (t) must lie in the opposite of the normal cone to C(t) at q(t):

f (t) ∈ -N C(t) (q(t)).
A case of major practical importance is when the admissible domain C(t) can be explicitly described by the non-negativity of a finite collection of p constraint functions (g i ) i∈{1,...,p} , each g i being a real function defined over M × R:

C(t) = {q ∈ M, ∀ i ∈ {1, ..., p}, g i (q, t) ≥ 0} .
If the constraint functions are differentiable with respect to q and satisfy a weak qualification hypothesis, then C(t) is tangentially regular and the normal cone to C(t) at q can be described in the following way :

N C(t) (q) = - p i=1 α i ∇g i (q, t), α ≥ 0, α ⊥ g(q, t) ,
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. [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF], Ref. [START_REF] Jean | The Non-Smooth Contact Dynamics Method[END_REF] or Ref. [START_REF] Renouf | 3D frictional contact and impact multibody dynamics: A comparison of algorithms suitable for real-time applications[END_REF]) 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.

Rigid body shapes regularity and smooth contact kinematics

In a first but restrictive definition attempt, we could consider a relative rigid body displacement h ∈ 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 ∂S i (t 0 ) ∩ h ∂S j (t 0 ). Let P be such a point for h = h * . Intuitively, the associated gap function g p : h ∈ SE(3) → 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 × 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. [START_REF] Montana | The Kinematics of Contact and Grasp[END_REF]) 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 vertexvertex 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. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF] 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. [START_REF] Park | Second-Order Contact Kinematics for Regular Contacts[END_REF], while in Ref. [START_REF] Glocker | Impacts with Global Dissipation Index at Reentrant Corners[END_REF] Glocker proposes the extension of impact laws that apply in those cases. 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. [START_REF] Pfeiffer | Relative Kinematics of Multibody Contacts[END_REF], and with a slightly different formalism in Refs. [START_REF] Duindam | Modeling the kinematics and dynamics of compliant contact[END_REF] and [START_REF] Visser | Screw Bondgraph Contact Dynamics[END_REF], is a prerequisite to a wide range of non-smooth and smooth contact models between rigid bodies (see for example Refs. [START_REF] Kraus | Compliant Contact Models For Rigid Body Collisions[END_REF], [START_REF] Song | Distributed Compliant Model for Efficient Dynamic Simulation of Systems with Frictional Contacts[END_REF] or [START_REF] Barbagli | Simulating human fingers: a Soft Finger Proxy Model and Algorithm[END_REF]).

Defining point contacts and non-negative gaps through LMDs

We claim that a natural and less restrictive definition of point contacts between two nonpenetrating 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 ) × 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 ) ∈ S i (t 0 ) × 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 ∂S i (t 0 )×h ∂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. [START_REF] Pfeiffer | Relative Kinematics of Multibody Contacts[END_REF] and [START_REF] Duindam | Modeling the kinematics and dynamics of compliant contact[END_REF]. Each resulting gap function g is hence defined at least on a non-empty subset of SE [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF]. 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).

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:

π(A, B) = inf τ , τ ∈ R 3 , (A + τ ) ∩ B = ∅. (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), π(h S j (t 0 ), S i (t 0 )) is attained for a unique translation vector τ * (h).

• If S i (t 0 ) and h S j (t 0 ) interpenetrate for some h, then

S i (t 0 ) ∩ (h S j (t 0 ) + τ * (h)) is an isolated point of ∂S i (t 0 ) ∩ h ∂S j (t 0 ) + τ * (h), denoted a(h).
With these conditions, the couple of points (a(h), a(h)τ * (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. [START_REF] Renouf | 3D frictional contact and impact multibody dynamics: A comparison of algorithms suitable for real-time applications[END_REF]). For convex polytopes, the inf bound in (1) is a min bound (see Ref. [START_REF] Agarwal | Penetration depth of two convex polytopes in 3D[END_REF]) but it is attained for some non-unique τ * 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-know 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. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF], [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF], [START_REF] Heidelberger | Consistent Penetration Depth Estimation for Deformable Collision Response[END_REF] and [START_REF] Redon | A Fast Method for Local Penetration Depth Computation[END_REF]). Their main weak point is that the use of history introduces hysteresis in the definition. For instance, Baraff proposes in Refs. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF], [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF] 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. [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF] and [START_REF] Johnson | 6-DOF Haptic Rendering Using Spatialized Normal Cone Search[END_REF]). The one found in Ref. [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF] closely resembles a local penetration depth definition and demands both local convexity and regularity assumptions. The one found in Ref. [START_REF] Johnson | 6-DOF Haptic Rendering Using Spatialized Normal Cone Search[END_REF] relies on discrete considerations.

Conforming contact situations and point contacts selection strategies

Now let us take again h ∈ 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 ∂S i (t 0 ) ∩ h ∂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. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF], and remarks in Ref. [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF] 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. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF], 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.

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 be can 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 ?

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. [START_REF] Lin | Collision detection between geometric models: A survey[END_REF], [START_REF] Jiménez | 3D collision detection: A survey[END_REF], [START_REF] Teschner | Collision Detection for Deformable Objects[END_REF], [START_REF] Klein | Efficient Collision Detection for Point and Polygon Based Models[END_REF] and the book Ref. [START_REF] Ericson | Real Time Collision Detection[END_REF] 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. [START_REF] Ponamgi | Incremental Algorithms for Collision Detection Between Polygonal Models[END_REF], [START_REF] Ehmann | Accurate and Fast Proximity Queries Between Polyhedra Using Convex Surface Decomposition[END_REF] and [START_REF] Kim | Six-Degree-of-Freedom Haptic Display Using Localized Contact Computations[END_REF]), interference volumes (see Refs. [START_REF] Heidelberger | Real-Time Volumetric Intersections of Deforming Objects[END_REF] and [START_REF] Hasegawa | Real-time Rigid Body Simulation for Haptic Interactions Based on Contact Volume of Polygonal Objects[END_REF]) or even more exotic approaches like intersection contours like in Ref. [START_REF] Volino | Resolving Surface Collisions through Intersection Contour Minimization[END_REF]. Some methods are dedicated to particular regularized contact models, like the approach found in Ref. [START_REF] Pauly | Quasi-Rigid Objects in Contact[END_REF] which is based on point cloud representations in the context of computer animation, or the one explained in Ref. [START_REF] Hippmann | An Algorithm for Compliant Contact Between Complexly Shaped Surfaces in Multibody Dynamics[END_REF] 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. [START_REF] Redon | Fast continuous collision detection between rigid bodies[END_REF]. The latter approach is closely to related to the determination of the instant of first contact, as exploited in event-driven methods (see again Baraff's excellent descriptions in Refs. [START_REF] Baraff | Analytical Methods for Dynamic Simulation of Non-penetrating Rigid Bodies[END_REF] and [START_REF] Baraff | Curved Surfaces And Coherence For Non-penetrating Rigid Body Simulation[END_REF]).

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).

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. [START_REF] Johnson | A Framework for Efficient Minimum Distance Computations[END_REF], [START_REF] Johnson | Unified Distance Queries in a Heterogeneous Model Environment[END_REF], [START_REF] Sohn | Computing the distance between two surfaces via line geometry[END_REF] and [START_REF] Hartmann | CAD-Based Contact Modelling For Multibody Systems[END_REF]). 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. [START_REF] Gottschalk | OBB-Tree: A hierarchical structure for rapid interference detection[END_REF] and the surveys cited in the preceding subsection) to attain tractable efficiency. The control of geometric approximation will be discussed in section 4.

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. [START_REF] Stewart | An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Inelastic Collisions and Coulomb Friction[END_REF], [START_REF] Anitescu | Time-stepping for the three dimensional rigid body dynamics[END_REF], [START_REF] Brogliato | Nonsmooth Mechanics[END_REF], [START_REF] Stewart | Rigid-Body Dynamics with Friction and Impact[END_REF] and [START_REF] Renouf | 3D frictional contact and impact multibody dynamics: A comparison of algorithms suitable for real-time applications[END_REF] 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.

SIMPLICIAL COMPLEXES AS VERSATILE DISCRETE SHAPE REPRESENTA-TIONS

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 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.

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. [START_REF] Munkres | Elements of Algebraic Topology[END_REF]), we recall that a simplex (or k-simplex) S in 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 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. [START_REF] Munkres | Elements of Algebraic Topology[END_REF] 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 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 Bd S. The simplicial interior of S is defined as S \ Bd S. We also recall that the relative interior (resp. the the relative boundary) of a subset A of 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 R d generated by A. The relative interior (resp. boundary) of a simplex coincides with its simplicial interior (resp. with its simplicial boundary). From a numerical point of view, simplicial complexes appear as generalized meshes, built upon 1-simplices of the complexes. The points in 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 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.

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 dependant 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 format 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 dilated1 simplicial complexes for representing geometric shapes of arbitrary topological dimension, as in Fig. 3. This is coherent with our preliminary assump-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.

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 (generalized) 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 visualization 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. Besides industrial tessellators exhibit fairly diverse output quality. However, high-quality tessellators 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 tetrehedral 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 multibody 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 performing the union all closed balls of radius ε having their centers in A. The resulting set is the Minkowski sum of A and the closed ball of radius ε centered at the origin. 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.

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.

A retraction 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 ε i -dilation of K i and ε j -dilation K j by deciding that the couple of points defining the generalized point contact kinematics is (a i + ε i n i , a j + ε j n j ), where n i and n j are defined by:

n i = -n j = 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 nonsmooth contact laws imposing non-negative gaps, although small positive and negative gap values are accessible for internal time-stepping computations2 . 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 surfacesor thin beams or cables. A comparable dilation method applied to convex polyhedra is described in Ref. [START_REF] Van Den Bergen | Proximity queries and penetration depth computation on 3d game objects[END_REF], but we have not found any reference which exploits it in the case of non-convex complex shapes.

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 timestepping 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 material 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.

Admissible domain boundary covering conditions

With non-convex simplicial complexes, constraint functions based on LMD gaps are defined only locally in the configuration manifold 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 (•, t) is defined and continuously differentiable, by V(X) the set of all neighborhoods of any set X ⊂ 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 (•, 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):

∀ t, ∀ i ∈ {1, ..., p}, ∃ Ω i ∈ V (g i (•, t)) -1 ({0}) ∩ ∂C(t) . Ω i ⊂ D i (t). (2)
Otherwise, "flip-flop" instabilities may occur in the corners of ∂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:

∀ t, ∀ x ∈ ∂C(t), ∃ Ω x ∈ V({x}), ∀ i ∈ {1, ..., p}, g i (x, t) = 0 ⇒ Ω x ⊂ 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 M. Such small displacements may cause both relative translations and rotations of the contacting rigid bodies, which explains the stabilizing effect of quasi-LMDs.

EFFICIENT QUASI-LMD COMPUTATIONS BASED ON AUGMENTED BOUND-ING 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. [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF] and Ref. [START_REF] Johnson | 6-DOF Haptic Rendering Using Spatialized Normal Cone Search[END_REF]). However the definition of normal cones given in Ref. [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF] 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.

Tangent and normal cones, Gauss map

We follow the terminology found for example in Ref. [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF], and adopt the following definitions for the tangent and normal cones.

Let X be a subset of R 3 and x a vector in X. A vector y ∈ R 3 is said to be a tangent of X at x if either y = 0 or there exists a sequence (x k ) k∈N such that x k ∈ X \ x for all k, lim k→∞ x k = x and lim k→∞ (x kx)/ x kx = 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 ∈ R 3 is said to be a normal of X at x is there exists sequences (x k ) k∈N and (z k ) k∈N , with respective limits x and z and such that for all k, x k ∈ X and z k ∈ 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 ) * ∩ S 2 , where S 2 is the unit sphere in R 3 . Now, let S be a d-simplex in R 3 with d > 0, σ a face of S of dimension d -1 and v the vertex of S opposite to σ. Let a be the orthogonal projection of v on the affine hull on of σ. The points a and v are distinct. The unit vector (va)/ av will be denoted u S (σ). Now let K be a finite simplicial complex in 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 σ of K, the function • For each i, let n K i (σ i ) ⊂ S 2 be a minimal generator of the polyhedral cone U K i (σ i ) * , and let C K i (σ i ) be the smallest semi-infinite revolution cone, with axis co-linear with the mean of the vectors in n K i (σ i ), which contains U K i (σ i ) * . The vector x jx i should be in the revolution cone obtained by augmenting of ε i (c K 1 (σ 1 ), c K 2 (σ 2 )) the half angle of C K i (σ i ).

x → T |K| (x) * dim(σ) dim T |K| (Int σ) * St K σ ∩ ∂|K| is a plane c K (σ) 0 3 • 0 0 2 • 1 1 2 • 2 0 1 false 3 0 1 true 4 1 1 • 5 2 1 • 6 • 0 • 7

Bounding volume hierarchies for efficient quasi-LMD computations

We propose to use a sphere-cone bounding volume hierarchy (BVH) similar to introduce in Ref. [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF], 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. [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF], but with the half-angles of the node cones augmented of the regularization parameter θ. 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. [START_REF] Fischer | The Smallest Enclosing Ball of Balls: Combinatorial Structure and Algorithms[END_REF]). For the bounding cones, we adopt a nested hierarchy, starting from the leafs and using the cones C K (σ).

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).

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. [START_REF] Jean | The Non-Smooth Contact Dynamics Method[END_REF]. This section shows experimental results demonstrating the robustness and efficiency of our approach.

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.

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 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.

Time per request (ms) (s, l) = (2, 10), "no opt." (s, l) = (2, 10) (s, l) = ( 1 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. [START_REF] Johnson | Spatialized Normal Cone Hierarchies[END_REF]. Other columns: our full 5D+1 BVH method.

"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.

CONCLUSION AND FUTURE WORK

Our collision detection method enables time-stepping schemes to be used robustly for nonsmooth 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.

Figure 1 :

 1 Figure 1: Left: complex shapes from industrial CAD models. Right: a conforming contact situation.

Figure 2 :

 2 Figure 2: Left: the usual vision of the gap values associated to a point contact P . Center: LMDs between nonpenetrating compact shapes. Left: typical issue concerning the locality of penetrations between non-convex shapes.

Figure 3 :

 3 Figure 3: Left: a simplicial complex. Right: A point contact between a plane and adilated shell.

Figure 4 :

 4 Figure 4: Left: condition (2) in configuration manifold M. Center: gaps between dilated simplicial complexes. Right: a LMD and a quasi-LMD (at the right hand side in green).

Figure 5 :

 5 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.

Figure 6 :

 6 Figure 6: Left: a highly constrained industrial part mounting scenario. Right: a peg-in-hole situation.

Table 1 :

 1 Classification of the simplices σ of the complex K. A dot indicates all possible values.

	ε(c 1 , c 2 )	0	1	2	3	4	5	6	7
	0	(0, 0) (0, 0) (θ, 0) (θ, 0) (θ, 0) (θ, 0) (θ, 0) n.d.
	1	(0, 0) (θ, θ) (θ, θ) (θ, 0) n.d.	n.d.	n.d. n.d.
	2	(0, θ) (θ, θ) (θ, θ) (θ, 0) n.d.	n.d.	n.d. n.d.
	3	(0, θ) (0, θ) (0, θ) n.d.	n.d.	n.d.	n.d. n.d.
	4	(0, θ) n.d.	n.d.	n.d.	n.d.	n.d.	n.d. n.d.
	5	(0, θ) n.d.	n.d.	n.d.	n.d.	n.d.	n.d. n.d.
	6	(0, θ) n.d.	n.d.	n.d.	n.d.	n.d.	n.d. n.d.
	7	n.d.	n.d.	n.d.	n.d.	n.d.	n.d.	n.d. n.d.

Table 2 :

 2 Definition of the values of the regularization function ε, the scalar θ 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.

The ε-dilation of a subset A in R 3 is an operation defined for any non-negative ε value which consists in

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.

ACKNOWLEDGEMENTS

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

takes a constant value over Int σ, which will be denoted T |K| (Int σ) * , and is described by:

where Dir Aff σ denoted the direction of the affine hull of σ and U K (σ) denotes the set

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

Characterization of LMDs between simplicial complexes

Let K 1 and K 2 be two finite and disjoint simplicial complexes in R

and only if the following condition holds:

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

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 ) ∈ X 1 × X 1 which locally minimizes d |K 1 |×|K 2 | is that the following condition holds:

(6)

Definition of quasi-LMDs

Let K 1 and K 2 be two finite and disjoint simplicial complexes in R 3 . The definition quasi-LMDs comes with a special classification of the faces σ of each complex K i , depending on the topological dimensions of σ and T |K| (Int σ) * , which is summarized in table 1. We define eight classes numbered from 0 to 7, and the class of σ in K is denoted c K (σ). Table 2 then defines a symmetric binary compatibility relation between the classes as well as a regularization function ε 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 σ 1 and σ 2 be faces of K 1 and K 2 respectively. A quasi-LMDs (x 1 , x 2 ) is supported by σ 1 and σ 2 if the following conditions hold:

• The classes of σ 1 and σ 2 are compatible, i.e. the value of ε(c K 1 (σ 1 ), c K 2 (σ 2 )) is defined.

• d Aff σ 1 ×Aff σ 2 has a unique minimizer which is (x 1 , x 2 ).

• (x 1 , x 2 ) is in Int σ 1 × Int σ 2 .

• For each i {1, 2}, for all σ ′ ∈ U K i (σ i ), x j -x i / x j -x i , u σ ′ (σ) < sin(ε i (c K 1 (σ 1 ), c K 2 (σ 2 )), with j ∈ {1, 2} \ i.


Matthieu Hamel
email: matthieu.hamel@cea.fr

Michel Vazeux

Mihaela Gulea

Diastereoselective oxidation of menthyl arenesulfenates to sulfinates and access to enantioenriched aryl methyl sulfoxides

Keywords: Sulfenate ester, sulfinate ester, sulfoxide, diastereoselective oxidation

Introduction

The interest of tri-coordinated sulfur derivatives and especially sulfoxides has successfully mutated from purely academic research to synthetic applications these last decades [START_REF] Wojaczyńska | Enantioselective Synthesis of Sulfoxides: 2000-2009[END_REF][START_REF] Otocka | Chiral Organosulfur Ligands/Catalysts with a Stereogenic Sulfur Atom: Applications in Asymmetric Synthesis[END_REF][START_REF] Wojaczyńska | Modern Stereoselective Synthesis of Chiral Sulfinyl Compounds[END_REF]. This success is mainly due to the discovery in the 1980s of the use of sulfoxides as efficient chiral auxiliaries for asymmetric synthesis. For example, enantiopure sulfoxides have found many applications in chiral induction [START_REF] Fernández | Recent Developments in the Synthesis and Utilization of Chiral Sulfoxides[END_REF], asymmetric catalysis [START_REF] Pellissier | Use of chiral sulfoxides in asymmetric synthesis[END_REF][START_REF] Han | Chiral sulfoxides: advances in asymmetric synthesis and problems with the accurate determination of the stereochemical outcome[END_REF], and drug design, with the wellknown Esomeprazole or its derivatives [START_REF] Cotton | Asymmetric synthesis of esomeprazole[END_REF][START_REF] Bentley | Role of sulfur chirality in the chemical processes of biology[END_REF][START_REF] O'mahony | Synthesis of enantioenriched sulfoxides[END_REF], Sulforaphane and its analogs (Kiełbasiński et al. 2014). Some years ago, we developed an innovative access to enantiopure ortho-phosphorylated phenylsulfoxides, by using, as the key step, an asymmetric oxidation of sulfenates to sulfinates [START_REF] Hamel | Asymmetric oxidation of sulfenates to sulfinates as a new route to optically active orthophosphorylated phenyl sulfoxides[END_REF][START_REF] Abrunhosa | New chiral ortho-P,Sdifunctionnalized aromatic compounds[END_REF]. ortho-Phosphorylated phenylsulfinate esters have been extensively studied by the groups of Drabowicz and Mikołajczyk, and such structures have found applications as mono-or bidentate ligands for the preparation of cisplatin derivatives [START_REF] Laforgia | Synthesis and Characterization of Platinum(II) Complexes of Diethyl [(Methylsulfinyl)methyl]phosphonate: Potential Drugs against Bone Tumors[END_REF][START_REF] Laforgia | Platinum Complexes with the Novel Ligand Diethyl [(Methylsulfinyl)-methyl]phosphonate (SMP): Solid-State Characterization of Potassium Trichloro(SMP)platinum(II) which, in Solution, Gives Dichloro(SMP)-platinum(II) and Potassium Chloride[END_REF][START_REF] Mauger | Reaction between the diaqua form of cisplatin and 2methylsulfanylphenylphosphonic acid yields a dinuclear phosphonato-bridged complex via NH3 elimination[END_REF][START_REF] Hamel | Ortho-(methylsulfanyl)phenylphosphonates and derivatives: Synthesis and applications as mono-or bidentate ligands for the preparation of platinum complexes[END_REF]. In this work, we demonstrated the potential of this approach by its extension to other arenesulfenate substrates, bearing an electron-withdrawing group (i.e. nitro or carboxylic ester) in the ortho or para position of the aryl moiety.

Materials and Methods

Reactions were carried out under nitrogen atmosphere. THF was purified with a PURESOLV TM apparatus developed by Innovative Technology Inc. Triethylamine (Et3N) was distilled over calcium hydride (CaH2). Carbon tetrachloride CCl4 was distilled over phosphorus pentoxide (P2O5) and stored over CaCl2 under a nitrogen atmosphere. The synthesis of methyl thiosalicylate was performed according to a published procedure (yield 83%; lit. 97% as the crude material [START_REF] Hamel | Asymmetric oxidation of sulfenates to sulfinates as a new route to optically active orthophosphorylated phenyl sulfoxides[END_REF]).

NMR spectra were recorded in deuterated chloroform (CDCl3), unless otherwise stated on either a 250 MHz or 400 MHz NMR spectrometer. Chemical shifts δ were indicated in ppm using tetramethylsilane (TMS) as the internal standard. 13 C { 1 H} NMR spectra were recorded at 62.9 or 100.6 MHz. Mass spectra were obtained on a GC/MS Saturn 2000 or on a Waters QTOF micro. Optical rotation values were measured on a Perkin-Elmer 241 automatic polarimeter; specific rotation and concentration are given in deg cm 3 g -1 dm -1 and g cm -3 , respectively. IR spectra were recorded with a Perkin-Elmer 16 PC FT-IR instrument. Satisfactory analytical data were obtained using a THERMOQUEST NA 2500 instrument. Thin layer chromatography (TLC) was performed on silica gel 60 F254, with UV detection or potassium permanganate solution (5% aqueous) staining. Flash chromatography was achieved on silica gel columns (40 -63 µm) using nitrogen pressure.

General procedure for the synthesis of sulfenate esters 1a-c. Abstract: Three aryl methyl sulfoxides were prepared via the following sequence: synthesis of the menthyl arenesulfenate from the corresponding thiol, diastereoselective oxidation leading to the sulfinate, then reaction with a Grignard reagent. The aryl moiety was substituted with a methyl ester on the ortho position or a nitro group on the ortho or the para position. The sulfenate esters were obtained in fair to excellent yield (42-82%).

Four different oxidizing agents were tested to obtain the corresponding sulfinate esters with diastereomeric excesses (de) ranging from 10 to 48%. After separation of the diastereomers and reaction with methyl Grignard reagent, two enantiopure sulfoxides, and one enantioenriched sulfoxide (32% ee) were obtained.

Sulfuryl chloride (1.3 equiv, 9.5 mmol) was slowly added to a previously stirred solution of arenethiol (1 equiv, 7.3 mmol) in dichloromethane (20mL). The resulting solution was stirred for 30 min at room temperature. After evacuation of the solvent and gases, tetrahydrofuran (THF, 25 mL) was added followed by (-)menthol (1 equiv, 7.3 mmol) and trifluoroacetic anhydride (0.1 equiv, 0.73 mmol). After cooling to -78 °C, triethylamine (2 equiv., 14.6 mmol) was added. The reaction mixture was allowed to reach room temperature, stirred for one hour and then hydrolyzed with a solution of sulphuric acid (H2SO4, 5% aqueous). The aqueous layer was extracted with diethyl ether. The combined extracts were washed with a 5% solution of potassium carbonate (K2CO3) followed by brine, dried over magnesium sulfate (MgSO4) and concentrated. The resulting oil was purified by silica gel chromatography.

Typical oxidation procedure for the synthesis of sulfinate esters 2a-c.

With mCPBA

Sulfenate ester 1 (1 0.5 mmol) was oxidized by the addition of meta-chloroperbenzoic acid (mCPBA, 1 equiv, 0.5 mmol) in CH2Cl2 at -78 °C, after which the reaction was allowed to reach room temperature. The solvent was evaporated and the crude sulfinate 2 directly purified by silica gel chromatography.

With NBS/H2O

N-Bromosuccinimide (NBS, 1 equiv., 1 mmol) was added portionwise to a stirred solution of menthyl sulfenate 1 (1 equiv., 1 mmol) in a mixture of acetonitrile and water 2:1. The concentration of sulfenate ester was ca. 50 mM. After addition of all the NBS, stirring was continued for 2 h and the mixture was quenched by adding a solution of H2SO4 (5% aqueous). The aqueous layer was extracted three times with dichloromethane. The organic layer was successively washed out with a solution of potassium carbonate (5% aqueous) and sodium thiosulfate (5% aqueous) in water. After drying and evaporation of the organic solvent, the crude extract, whose purity was usually satisfactory as shown by the NMR spectrum, was purified on a silica gel column with pentane-ethyl acetate as the eluent to obtain menthyl sulfinates 2.

With (+)-or (-)-(8,8-dichlorocamphorylsulfonyl)oxaziridine

The chiral (8,8-dichlorocamphorylsulfonyl)oxaziridine (1 equiv, 0.5 mmol) was added portionwise to a stirred solution of sulfenate ester 1 (1 equiv, 0.5 mmol) in carbon tetrachloride under a nitrogen atmosphere. After complete conversion, the solvent was evaporated and the residue dissolved in pentane. The mixture was filtered and washed once with pentane. The sulfinate was further purified on silica gel chromatography when necessary.

Grignard reaction with sulfinate esters

Methyl magnesium bromide (MeMgBr, 1.1 equiv, 1.21 mmol) was added dropwise to a stirred solution of menthyl sulfinate 2a-c (1 equiv, 1 mmol) in 10 mL of diethyl ether at 0 °C. The resulting mixture was stirred for 1 h at room temperature, after which 10 mL of a 5% aqueous solution of sulphuric acid was added and the aqueous layer extracted with diethyl ether. The organic extracts were washed with a solution of potassium carbonate (5% aqueous) followed by brine, dried over MgSO4, filtered, evaporated and the mixture was further purified on silica gel chromatography, yielding 3a-c. No scale-up experiment was performed in this study.

Numeral identifications of carbon atoms for compounds listed below are arbitrary and given only for helping the description of NMR spectra, e.g. with compound 1b.

(-)-Menthyl 2-(methoxycarbonyl)benzenesulfenate 1a. 78% yield; Rf 0.63 (pentane/ethyl acetate 90:10). White solid; mp 69 °C;

[α]D 20 -53.9 (c 0.9 in CHCl3); 1 H NMR (250 MHz, δ): 7.97 (dd, 1H, 3 J = 7.8, 4 J = 1.3, H6), 7.79 (d, 1H, 3 J = 8.2, H3), 7.54 (dt, 1H, 3 J = 7.2, 4 J = 1.4, H4), 7.14 (dt, 1H, 3 J = 7.8, 4 J = 1.0, H5), 3.92 (s, 3H, CH3), 3.54 (dt, 3 J = 10.6, 3 J = 4.1, H1'), 2.42 -2.48 (m, 1H, H6'), 2.27 -2.41 (m, 1H, H3'), 1.61 -1.72 (m, 2H, H3', H5'), 1.41 -1.58 (m, 1H, H4'), 1.23 -1.40 (m, 2H, H4', H6'), 1.01 -1.14 (m, 1H, H8'), 0.97 (d, 3H, 3 J = 7.0, H7'), 0.90 -0.97 (m, 1H, H3'), 0.88 (d, 6H, 3 J = 6.4, H9'); 13 C NMR (62.9 MHz, δ): 167. (-)-Menthyl 2-nitrobenzenesulfenate 1b; 82% yield; Rf 0.74 (pentane/ethyl acetate 90:10). Yellow solid; mp 91 °C; [α]D 20 -17.0 (c 1 in CHCl3); 1 H NMR (400 MHz, δ): 8.28 (dd, 1H, 3 J = 8.4, 4 J = 1.3, H3), 7.90 (dd, 1H, 3 J = 8.2, 4 J = 1.0, H6), 7.66 (dt, 1H, 3 J = 7.1, 4 J = 1.3, H5), 7.26 (dt, 1H, 3 J = 7.3, 4 J = 1.1, H4), 3.61 (dt, 1H, 3 J = 10.7, 3 J = 4.2, H1'), 2.40 (ddd, 1H, 3 J = 7.0, 3 J = 6.9, 4 J = 2.7, H6'), 2.25 -2.28 (m, 1H, H3'), 1.64 -1.72 (m, 2H, H3', H5'), 1.46 -1.52 (m, 2H, H2', H4'), 1.30 -1.42 (m, 2H, H4', H6'), 1.06 -1.14 (m, 1H, H8'), 0.99 (d, 3H, 3 J = 7.0, H7'), 0.90 (d, 3H, 3 J = 6.5, H9'), 0.88 (d, 3H, 3 J = 6.9, H9'); 13 ; 1 H NMR (400 MHz, δ): 8.18 (dt, 2H, 3 J = 9.0, 4 J = 2.0, H3), 7.31 (dt, 2H, 3 J = 9.0, 4 J = 1.9, H2), 3.53 (dt, 1H, 3 J = 10.7, 3 J = 4.3, H1'), 2.22 -2.37 (m, 2H, H3', H6'), 1.62 -1.71 (m, 2H, H4', H5'), 1.30 -1.49 (m, 2H, H2', H3'), 0.97 (d, 3H, 3 J = 7.0, H7'), 0.94 -1.07 (m, 2H, H6', H8'), 0.91 (d, 3H, 3 J = 6.6, H9'), 0.80 -0.87 (m, 1H, H4'), 0.82 (d, 3H, 3 J = 7.0, H9'); 13 C NMR (62.9 MHz, δ): 151.9 (C4), 144.1 (C1), 123.0 (C6), 119.4 (C5), 87.9 (C1'), 47.8 (C2'), 39.9 (C6'), 33.0 (C4'), 30.8 (C5'), 24.5 (C8'), 22.2 (C3'), 21.1 (C7'), 20.0 (C9'), 15.0 (C9'); IR (KBr) ν = 1576, 1507 (NO2), 1337; HRMS (ESI, m/z): [M + Na] + calcd for C16H23NO3NaS: calcd 332.1296, found 332.1303.

(-)-Menthyl (S)-2-(methoxycarbonyl)benzenesulfinate (S)-2a; 68% yield; Rf 0.3 (pentane/diethyl ether 70:30). White solid; mp 146 °C; [α]D 20 -360.0 (c 0.6, CHCl3); Lit. [START_REF] Klunder | Convenient synthesis of sulfinate esters from sulfonyl chlorides[END_REF]: -348.6 (c 1.18, acetone); 1 H NMR (400 MHz, δ): 8.26 (dd, 1H, 3 J = 8.0, 4 J = 1.2, H6), 7.99 (dd, 1H, 3 J = 7.6, 4 J = 1.2, H3), 7.71 (dt, 1H, 3 J = 7.0, 4 J = 1.2, H4), 7.54 (dt, 1H, 3 J = 7.1, 4 J = 1.2, H5), 4.08 (dt, 1H, 3 J = 10.8, 3 J = 4.4, H1'), 3.91 (s, 3H, CH3), 2.28 -2.31 (m, 1H, H6'), 2.03 -2.07 (m, 1H, H3'), 1.60 -1.64 (m, 2H, H3', H5'), 1.33 -1.47 (m, 1H), 1.23 -1.28 (m, 1H), 0.94 -1.14 (m, 2H), 0.91 (d, 3H, 3 J = 6.4, H7'), 0.80 (d, 3H, 3 J = 6.8, H9'), 0.76 -0.84 (m, 1H, H3'), 0.67 (d, 3H, 3 J = 7.2, H9'); 13 C NMR (62.9 MHz, δ): 166.0 (C=O), 148.7 (C2), 133.5 (C4), 131.7 (C6), 131.1 (C3), 128.4 (C5), 124.9 (C1), 80.4 (C1'), 52.9 (CH3O), 48.3 (C2'), 42.5 (C6'), 34.4 (C4'), 32.1 (C5'), 25.6 (C8'), 23.6 (C3'), 22.5 (C7'), 21.2 (C9'), 15.9 (C9').

(-)-Menthyl (S)-2-nitrobenzenesulfinate (S)-2b; 78% yield (for both diastereomers); Rf 0.47 (pentane/diethyl ether 60:40). Clear yellow solid; [α]D 20 -413.9 (c 1.95, CHCl3); Lit. [START_REF] Hajipour | Highly diastereoselective synthesis and easy method for synthesis of optically active sulfinate esters from aromatic disulfides[END_REF]: -447.3 (c 1.5, acetone); 1 H NMR (400 MHz, δ): 8.36 (dd, 1H, 3 J = 7.6, 4 J = 1.2, H6), 8.23 (dd, 1H, 3 J = 8.0, 4 J = 1.2, H3), 7.92 (dt, 1H, 3 J = 7.6, 4 J = 1.2, H4), 7.73 (dt, 1H, 3 J = 7.6, 4 J = 1.6, H5), 4.16 (dt, 1H, 3 J = 10.8, 3 J = 4.4, H1'), 2.34 -2.39 (m, 1H, H6'), 2.03 -2.06 (m, 1H, H3'), 1.65 -1.71 (m, 2H, H2', H6'), 1.45 -1.55 (m, 2H, H4'), 1.04 -1.17 (m, 2H, H5', H8'), 0.98 (d, 3H, 3 J = 6.8, H7'), 0.82 (d, 3H, 3 J = 7.2, H9'), 0.82 -0.88 (m, 1H, H3'), 0.76 (d, 3H, 3 J = 6.8, H9'); 13 C NMR (100.6 MHz, δ): 145.7 (C2), 142.8 (C1), 134.9 (C4), 132.4 (C5), 126. (-)-Menthyl (R)-2-nitrobenzenesulfinate (R)-2b; Rf 0.38 (pentane/diethyl ether 60:40). Clear yellow solid [α]D 20 +266.2 (c 1.42 in CHCl3); 1 H NMR (400 MHz, δ): 8.32 (dd, 1H, 3 J = 8.0, 4 J = 1.6, H6), 8.23 (dd, 1H, 3 J = 8.0, 4 J = 1.2, H3), 7.90 (dt, 1H, 3 J = 7.6, 4 J = 1.2, H4), 7.72 (dt, 1H, 3 J = 8.0, 4 J = 1.6, H5), 4.26 (dt, 1H, 3 J = 10.8, 3 J = 4.4, H1'), 2.22 -2.28 (m, 1H, H6'), 1.78 -1.85 (m, 1H, H3'), 1.62 -1.66 (m, 2H, H2', H6'), 1.41 -1.51 (m, 1H, H4'), 1.18 -1.33 (m, 2H, H5', H8'), 0.99 -1.10 (m, 1H, H4'), 0.91 (d, 3H, 3 J = 6.8, H7'), 0.76 -0.84 (m, 1H, H3'), 0.78 (d, 3H, 3 J = 6.8, H9'), 0.74 (d, 3H, 3 J = 6.8, H9'); 13 C NMR (100.6 MHz, δ): 145.3 (C2), 142.9 (C1), 135.0 (C4), 132.5 (C5), 126.1 (C6), 125.0 (C3), 83.1 (C1'), 48.4 (C2'), 43.4 (C6'), 33.9 (C4'), 31.8 (C5'), 25.0 (C8'), 23.1 (C3'), 22.1 (C7'), 21.0 (C9'), 15.4 (C9').

(-)-Menthyl 4-nitrobenzenesulfinate 2c (mixture of diastereomers); 75% yield; Rf 0.2 (pentane/ethyl acetate 90:10). White solid; mp 75 °C; [α]D 20 -100.7 (c 0.4, CHCl3); for a de 32%, dr 34/66; Lit. [START_REF] Hajipour | Highly diastereoselective synthesis and easy method for synthesis of optically active sulfinate esters from aromatic disulfides[END_REF]: -438 (c 1.2, acetone); 1 H NMR (400 MHz, δ): 8.39 (t, 2H, 3 J = 8.8, H3), 7.92 (d, 2H, 3 J = 8.8, H2), 4.23 & 4.28 (dt, 1H, 3 J = 10.7, 3 J = 4.5, H1'), 0.72 -2.36 (m, 20H); 13 C NMR (62.9 MHz, δ): 152.4, 150.2, 126.6, 126.1, 124.4, 124.3, 83.7, 81.7, 48.8, 48.0, 43.7, 42.9, 34.0, 33.9, 32.0, 31.9, 25.8, 25.5, 23.3, 23.2, 22.2, 22.0, 21.0, 20.9, 15.8, 15.6;IR (KBr) ν = 1529, 1346, 1141, 763; EIMS (m/z (%)): 326 (1) [M + ], 170 (7), 139 (33), 83 (100), 69 (19), 55 (48).

(R)-(+)-Methyl 2-methylsulfinylbenzoate (R)-3a; Yield 25%; Rf 0.2 (ethyl acetate/pentane 80:20). White solid; [α]D 20 +139.0 (c 1, CHCl3); Lit. [START_REF] Pitchen | An efficient asymmetric oxidation of sulfides to sulfoxides[END_REF], for 60% ee: +84.0 (c 0.4, acetone); 1 H NMR (400 MHz, δ): 8.32 (dd, 1H, 3 J = 7.9, 4 J = 1.2), 8.09 (dd, 1H, 3 J = 7.8, 4 J = 1.2), 7.84 (dt, 1H, 3 J = 8.0, 4 J = 1.3), 7.58 (dt, 1H, 3 J = 7.1, 4 J = 1.2), 3.96 (s, 3H, CH3O), 2.85 (s, 3H, CH3).

(R)-(+)-Methyl 2-nitrophenyl sulfoxide (R)-3b; Yield 62%; Rf 0.2 (ethyl acetate). White solid; [α]D 20 +164.6 (c 1, CHCl3); Lit. [START_REF] Brunel | Catalytic enantioselective oxidation of sulfides with a chiral titanium complex[END_REF], for 75% ee: +118.5. 1 H NMR (400 MHz, δ): 8.03 (d, 1H, 3 J = 8.0), 6.92 -7.01 (m, 2H), 6.62 (t, 1H, 3 J = 7.1), 2.22 (s, 3H, CH3).

Methyl 4-nitrophenyl sulfoxide 3c; Yield 59%; Rf 0.2 (ethyl acetate). White solid; [α]D 20 +50.0 (c 1 in CHCl3); Lit. [START_REF] Brunel | Highly Enantioselective Oxidation of Sulfides Mediated by a Chiral Titanium Complex[END_REF], for 99.3% ee: +156.9 (c 0.75 in CHCl3); 1 H NMR (400 MHz, δ): 8.36 (d, 2H, 3 J = 8.2), 7.80 (d, 2H, 3 J = 7.9), 2.28 (s, 3H, CH3).

Results and Discussion

In order to demonstrate its potential, the sequence we described earlier [START_REF] Hamel | Asymmetric oxidation of sulfenates to sulfinates as a new route to optically active orthophosphorylated phenyl sulfoxides[END_REF]) has been extended to other substrates. Thus, aryl methyl sulfoxides where the aromatic ring is functionalized with either a nitro (in ortho or para position) or an ortho-methoxycarbonyl groups were prepared following this strategy. The overall reaction sequence is described in Scheme 1. Some of the diastereomerically pure sulfinate esters 2a, b and c have previously been reported by other groups. 2a was prepared starting from the corresponding sulfinyl chloride which reacts with menthol, with a low diastereoselectivity [START_REF] Klunder | Convenient synthesis of sulfinate esters from sulfonyl chlorides[END_REF]. Compounds 2b and 2c were obtained from the reaction of ortho-or para-nitrophenyldisulfide with menthol in the presence of lead tetraacetate, with a moderate diastereomeric excess [START_REF] Hajipour | Highly diastereoselective synthesis and easy method for synthesis of optically active sulfinate esters from aromatic disulfides[END_REF]. Another method consists in reacting the sulfonyl chloride and the alcohol in the presence of trimethylphosphite in dichloromethane [START_REF] Lee | Chiral acetylenic sulfoxide in enantioselective synthesis of tetrahydroisoquinoline and tetrahydro-β-carboline alkaloids. Total synthesis of (R)-(+)-Carnegine and (R)-(+)-tetrahydroharman[END_REF] or under solid state conditions [START_REF] Hajipour | A convenient and highly diastereoselective synthesis of optically active sulfinate esters from sulfonyl chlorides under solid state conditions[END_REF]. Crystal structure of (S)-2a has even been assigned, with a hydrogen bond linking the oxygen of the sulfinate ester with the hydrogen located in the para position of the methoxycarbonyl group [START_REF] Altamura | 1R,2S,5R)-(-)-Menthyl (S)-2-(methoxycarbonyl)benzenesulfinate[END_REF].

Scheme 1. Synthesis of sulfoxides 3a-c from the corresponding arylthiols.

In the literature the three sulfoxides 3a-c have been prepared in an (almost) enantiopure form. Some generic methods for their preparation involve either the enantioselective oxidation of the corresponding sulfide, catalyzed for instance with an organometallic chiral complex, with an enzyme, or via the synthesis of chiral sulfinate esters. Table 1 gives some examples of the oxidation of the prochiral sulfide leading to the sulfoxides 3a-c. These examples were sorted so that both enantiomeric excess and chemical yield were optimum.

According to the description of the synthesis given in Scheme 1, we started from commercial aromatic thiols. Methyl thiosalicylate was prepared according to a published reaction [START_REF] Sasaki | Polycyclic Nhetero compounds. XXXVIII. Synthesis and evaluation of antidepressive activity of [1]benzothieno[2',3':4,5]furo[3,2d]pirimidines and their precursors[END_REF]. The thiols were first converted into sulfenyl chlorides, then reacted with (-)-menthol in the presence of triethylamine and a catalytic amount of trifluoroacetic acid anhydride (TFAA). Enantiopure menthyl sulfenate esters were thus obtained in good yields, from 50 up to 78%. The second step consists in the diastereoselective oxidation to afford the menthyl sulfinates. Various methods can be used for the oxidation of sulfenate esters to their corresponding sulfinates. We used both achiral and chiral oxidants, the two chiral compounds being the (+)-and the (-)-isomers of (8,8dichlorocamphorylsulfonyl)oxaziridine (hereafter, designated as (+)-Ox and (-)-Ox, for the sake of clarity), whereas simple oxidants were either mCPBA or the couple NBS/H2O. Oxidations of 1b and 1c were performed only with achiral oxidants. Diastereoselectivity was measured following the 1 H NMR chemical shift from the CH-O of the menthyl moiety. The results are summarized in Table 2. .0 (c 1, CHCl3) a absolute configuration at sulfur atom b dr after separation of diastereomers All these oxidations performed smoothly and led to good conversions, but took an unexpected long time when performed with the oxaziridine (up to 15 days). As previously stated, oxidation with NBS/H2O couple did not favor the same diastereomer as when performing the oxidation with mCPBA [START_REF] Hamel | Asymmetric oxidation of sulfenates to sulfinates as a new route to optically active orthophosphorylated phenyl sulfoxides[END_REF]. For the three sulfinate esters 2, the observed diastereomeric excess was in the range of 10-48%, so slightly lower than those obtained with ortho-phosphorylated menthyl benzenesulfinate [START_REF] Hamel | Asymmetric oxidation of sulfenates to sulfinates as a new route to optically active orthophosphorylated phenyl sulfoxides[END_REF]. Moreover, the best diastereomeric ratio values were obtained when oxidizing compounds 1 with the NBS/H2O couple. Surprisingly, the best de was observed with 1c, suggesting that the position of the substituent had almost no influence on the steric hindrance, the latter being therefore driven by the menthyl moiety. Therefore, one or the other diastereomer can be favored either by using mCPBA or NBS/H2O or by using (-)-or (+)-menthol. Starting from the diastereomeric mixtures obtained respectively in entries 2 and 6 from Table 2, the diastereomerically pure major sulfinates 2a and 2b have been successfully separated using silica gel chromatography or by separative crystallization in acetone. In the case of 2c (entry 8, Table 2), the diastereomer separation was not possible. Specific rotation for optically pure compounds 2a and 2b are negative with respective values of -360.0 and -453.2 (c 0.6 and 1 in chloroform). These are slightly more negative than reported in the literature with respective [α]D 20 of -348.6 and -447.7 (c 1.18 [START_REF] Klunder | Convenient synthesis of sulfinate esters from sulfonyl chlorides[END_REF] and 1.5 in acetone [START_REF] Hajipour | Highly diastereoselective synthesis and easy method for synthesis of optically active sulfinate esters from aromatic disulfides[END_REF]). According to the literature, these sulfinate esters have the (S) configuration at the sulfur atom. From the two optically pure sulfinate esters (SS)-2a and (SS)-2b and the optically enriched compound 2c, the preparation of the corresponding methyl aryl sulfoxides 3a-c has been performed. As methylmagnesium bromide is able to attack both the electrophilic sulfur atom and the ester function of a sulfinate, to limit the second undesired process, the reaction was performed at 0 °C with a very slight excess of Grignard reagent. However, a low yield was obtained for 3a (25%), whereas 3b and 3c were isolated with satisfactory yields of 62% and 59%, respectively (Table 3). The comparison of the specific rotation values of the obtained aryl methyl sulfoxides with data from the literature (see in the experimental part), allowed us to assign the absolute (R) configuration at sulfur in the three cases. This was expected to be obtained after full inversion of configuration at the sulfur atom from the starting corresponding (Ss)-sulfinates 2a-c leading to sulfoxides (R)-3a and (R)-3b in an optically pure form (entries 1 and 2, Table 3), and (R)-3c with 32% ee (entry 3, Table 3). Despite the fact that this method was exemplified on a small number of substrates, the obtained results are encouraging for application on other arylsulfides bearing an electron withdrawing substituent on the aromatic ring. The oxidation of sulfenate esters was performed with simple oxidants such as mCPBA or NBS/H2O and the asymmetric induction was provided by the cheap enantiopure menthol. Since the chiral induction is provided by the alcohol, one can expect that this strategy should work not only for arylsulfenates but also for alkylsulfenates as well.

Conclusion

In conclusion, an efficient method for the access to three ortho and para substituted aryl methyl sulfoxides is described herein.

The diastereoselective oxidation of sulfenate to sulfinate ester is the key step, leading to diastereomeric excesses ranging from 10 to 48%. After separation of the diastereomers and reaction with methyl Grignard, two enantiopure sulfoxides, the (Rs)-(2methoxycarbonylphenyl) methyl sulfoxide and the (Rs)-2nitrophenyl methyl sulfoxide, and one enantioenriched (32% ee) were obtained.

  4 (C=O), 151.2 (C2), 132.7 (C4), 130.3 (C6), 123.3 (C3), 121.3 (C5), 121.0 (C1), 86.6 (C1'), 52.4 (CH3O), 49.0 (C2'), 41.0 (C6'), 34.3 (C4'), 31.9 (C5'), 25.6 (C8'), 23.3 (C3'), 22.2 (C7'), 21.1 (C9'), 16.1 (C9'); IR (KBr) ν = 1690 (C=O st), 757; EIMS (m/z (%)): 322 (1) [M + ]; 184 (85), 167 (26), 152 (100), 83 (33), 69 (15), 55 (30), 41 (15), Anal. calcd. for C18H26O3S: C 67.04,H 8.13, S 9.94, found: C 66.47, H 8.20, S 10.47.

  C NMR (62.9 MHz, δ): 148.0 (C2), 134.1 (C1), 134.0 (C4), 125.1 (C5), 124.3 (C6), 123.1 (C3), 87.0 (C1'), 49.0 (C2'), 41.0 (C6'), 34.2 (C4'), 31.8 (C5'), 25.7 (C8'), 23.3 (C3'), 22.1 (C7'), 21.0 (C9'), 16.1 (C9'); IR (KBr) ν = 1504 (NO2), 734; EIMS (m/z (%)): 332 (100) [M + Na] + , 194 (32); HRMS (ESI, m/z): [M + Na] + calcd for C16H23NO3NaS: 332.1296, found 332.1315; Anal. calcd. for C16H23NO3S: C 62.11, H 7.50; found C 62.03, H 7.99. (-)-Menthyl 4-nitrobenzenesulfenate 1c; 42% yield; Rf 0.78 (pentane/ethyl acetate 90:10). Yellow solid; mp 70 °C; [α]D 20 -162.3 (c 1.15 in CHCl3)

  3 (C6), 125.0 (C3), 81.1 (C1'), 48.0 (C2'), 41.4 (C6'), 34.0 (C4'), 31.9 (C5'), 25.3 (C8'), 23.3 (C3'), 22.2 (C7'), 20.9 (C9'), 15.6 (C9'); IR (KBr) ν = 1531 (NO2), 1456, 1349, 1137, 949, 757; EIMS (m/z (%)): 348 (11) [M + Na], 228 (9), 210 (100); HRMS (ESI, m/z): [M + Na] + calcd for C16H23NO4NaS: 348.1245, found 348.1247.

Table 1 . Some examples of selective sulfide oxidation leading to enantiopure sulfoxides 3a-c.

 1

	Sulfoxide	Config. Yield (%)	ee (%) Method	Ref.
	3a	(R)	53	99	Ti(Oi-Pr)4/(R,R)-DET 1:4	De Sio et al. 2010
	3a	(R)	95	80	Mn(OTf)2/porphyrin deriv./H2O2	Dai et al. 2014
	3a	(S)	n.d.	82	Myoglobin active site mutants L29H/H64L	Ozaki et al. 1999
	3b	(R)	90	92	H2O2, chiral amino alcohol, VO(acac)2	Aydin 2013
	3b	(S)	84	99	H2O2, Al salalen	Matsumoto et al. 2008
	3c	(R)	95	99.5	H2O2, MgSO4, chiral Brønsted acid	Liao et al. 2012
	3c	(S)	92	92	H2O2, Ti Salen	Saito et al. 2001
	3c	(S)	66	96	N,N'-dioxide-iron(III) complex	Wang et al. 2020

Table 2 . Diastereoselective oxidation of arenesulfenate esters 1

 2

	Entry	Sulfenate	Oxidant	Conditions	Sulfinate	Yield (%)	dr	de (%)
	1	1a	NBS/H2O	CH3CN/H2O, 20 °C, 2 h 2a	68	69:31	38
	2	1a	mCPBA	CH2Cl2, -78 → 20 °C	2a	90	45:55	10
	3	1a	(+)-Ox	CCl4, 20 °C, 15 d	2a	56	39:61	22
	4	1a	(-)-Ox	CCl4, 20 °C, 15 d	2a	63	65:35	30
	5	1b	NBS/H2O	CH3CN/H2O, 20 °C, 2 h 2b	6	72:28	44
	6	1b	mCPBA	CH2Cl2, -78 → 20 °C	2b	78	42:58	16
	7	1c	NBS/H2O	CH3CN/H2O, 20 °C, 2 h 2c	69	74:26	48
	8	1c	mCPBA	CH2Cl2, -78 → 20 °C	2c	75	34:66	32

Table 3 . Transformation of arenesulfinate esters 2 into enantioenriched sulfoxides 3

 3

	Entry Sulfinate dr (config.) a	Sulfoxide	Yield (%)	ee (%) (config.) a Measured [α]D 20
	1	2a	100:0 b (S)	3a	25	100 (R)	+84.0 (c 1, CHCl3)
	2	2b	100:0 b (S)	3b	62	100 (R)	+164.6 (c 1, CHCl3)
	3	2c	34:66 (S)	3c	59	32 (R)	+50

Acknowledgements

The Authors gratefully acknowledge financial support from the Ministère de la Recherche for the grant allowed to MH for his PhD thesis. The Reviewers and the Editor are also gratefully acknowledged for their valuable comments.

Conflict of interest statement

The Authors declare no conflict of interest.