Keywords: collision, contact, simulation, physically based animation

We propose a new approach for collision modeling in physically based animation. Contrary to most current approaches, our method can be used with all physical models rather than any specific class. At each time step, the geometry of the bodies is mapped to a sparse regular Eulerian grid. Each grid point carries a mass, a velocity and a spatial density. This grid acts as a common mechanical layer where detection, modeling and reaction to collision take place, without any assumption about the internal dynamics of the bodies in contact. Mappings are then used to propagate penalty-and constraint-based reactions back to the bodies. We show that mappings can be easily set up for the most commonly used physical models. Our approach greatly simplifies the implementation of collision modeling since we only have to consider each body's individual mapping to the Eulerian grid, rather than numerous model pair-specific methods. Moreover, it allows us to design and reuse efficient collision response strategies independently of the physical models. We demonstrate our method with a variety of models including rigid bodies, deformable solids and fluids.

Contact Eulérien pour le Traitement Unifié des Collisions

Résumé : Nous présentons une nouvelle approche pour la modélisation du contact en animation par modèles physiques. Contrairement à la plupart des autres approches, notre méthode peut être appliquée à tous les modèles physiques plutôt qu'à une classe particulière. À chaque pas de temps, la géométrie des corps est plaquée sur une grille eulérienne par un mapping. Chaque point de la grille contient une masse, une vitesse et une densité spatiale. Cette grille joue le rôle d'une couche mécanique commune dans laquelle la détection, la modélisation et la réaction aux collision se produisent, sans connaissance de la dynamique interne des objets en contact. Le mapping est alors utilisé pour rétropropager les forces et contraintes de contacts vers les objets. Nous montrons que les mappings sont simples à construire pour la plupart des objets. Notre approche simplifie considérablement l'implémentation des contacts car il suffit de détermnier le mapping de chaque corps vers la grille eulérienne, plutôt que traiter une multitude de paires de modèles au cas-par-cas. De plus, elle permet de concevoir et réutiliser des stratégies efficaces de réponses aux collisions, indépendamment des modèles physiques. Nous en faisons la démonstration desur de nombreux modèles dont des solides rigides, des objets viscoélastiques et des fluides.

Mots-clés : collision, contact, simulation, animation par modèles physiques 1 Introduction

The field of collision processing has been thoroughly investigated for the last twenty years. However, none of the approaches proposed so far is both generic and efficient. Most methods are dedicated to a single type of contact such as triangle-based meshes. This makes the development of versatile physically based animation systems difficult involving the implementation and the run-time combination of numerous algorithms and techniques to deal with the variety of geometries (spheres, boxes, cylinders, lines or polygons) and mechanical behaviors (rigid, deformable, fluid). Using a new kind of geometric primitive requires implementing its intersection detection algorithm with each geometric primitive already used. This motivates us to build a totally new collision modeling method in which arbitrary mechanical models are mapped to a common collision model designed to detect and model contact efficiently.

Overview

Our model treats contact using a hierarchy of clearly separated layers, as illustrated in figure 1. The top layer contains arbitrary Lagrangian bodies governed by their own internal laws. Any type of physical body can be used, provided that a spatial description of its volume is available. We call this layer the Lagrangian layer. The lower layers allow interaction between bodies the top layer.

The second layer models the Cartesian space occupied by the bodies. At each time step, each body is sampled over a set of Eulerian points representing voxels with an associated filling ratio, which we call density. We call this layer the Eulerian layer. The physical interactions between the bodies are modeled at this level. Each body is sampled over the same set of points, resulting in a finite number of possible contact points.

The third layer, called the collision layer, is used for collision detection. It represents the sum of the volumetric images of the bodies. Interaction forces are applied at the collision points in the direction of the gradient of density, and are then mapped to generalized forces applied to the bodies. This allows them to react to collisions.

Our contribution is to introduce a new contact modeling method based on Eulerian density grids, with the following properties:

• Any model can interact any other, provided a mapping from the body to the Eulerian layer is available.

• Density fields have the same versatility as distance fields for collision detection and modeling; however, they are much easier to resample at each time step because density is an additive, local property. Moreover, estimates of contact directions have better stability.

• Levels of detail are available for collision detection and modeling using coarse or fine versions of the density grids, allowing an easy trade off between geometric accuracy and speed.

• Collision response strategies, both penalty-and constraint-based, can be designed and re-used independently from the body models. This contrasts with traditional approaches which are targeted to specific models.

The remainder of this paper is organized as follows. Previous work is summarized in section 2. Section 3 discusses contact modeling and reaction more specifically. Section 4 shows how contact penalties and constraints are expressed in the Eulerian layer, and which operators are necessary to achieve the mechanical coupling with the bodies. Section 5 explains how to build a collision model for a variety of bodies at each time step. An architecture for versatile stable time integration is proposed in section 6. Results are presented in section 7, and a discussion is given in section 8.

Previous work

Collision detection and response techniques are well studied in the Computer Graphics community, and the contributions are too numerous to list here. Many different approaches have been proposed, depending on the type of objects in contact. While we provide a brief description below, please refer to [LG98, JTT01, TKH * 04] for excellent surveys.

Rigid-body collision processing has been largely studied in the last two decades. While early work detected collisions between convex polyhedra [GJK88, GASF94, CLMP95] For non-convex polyhedra, various approaches using hierarchical bounding volumes were developed [BO04, GLM96, vdB97, KHM * 98]. Alternatively, distance fields based methods have also been used [START_REF] Fuhrmann | Distance fields for rapid collision detection in physically based modeling[END_REF][START_REF] Funfzig | Hierarchical spherical distance fields for collision detection[END_REF]. Their advantage is that they do not have restrictions about the object topology and provide the penetration depth for collision response.

Deformable objects are not able to exploit as much precomputed data as rigid bodies, hence they require more complex collision detection algorithms. Many approaches use a dynamic hierarchy of bounding volumes [START_REF] Zachmann | Geometric data structures for computer graphics[END_REF], and significant work has been done about the update of the hierarchy as body shape changes [LAM01, [START_REF] Mezger | Hierarchical Techniques in Collision Detection for Cloth Animation[END_REF][START_REF] James | Bd-tree: output-sensitive collision detection for reduced deformable models[END_REF]. Other approaches use body-space hierarchies [START_REF] Volino | Collision and self-collision detection: Efficient and robust solutions for higly deformable surfaces[END_REF][START_REF] Provot | Collision and self-collision handling in cloth model dedicated to design garments[END_REF]. Recent work explores the use of volumetric data-structures, such INRIA as spatial hashing [THM * 03], GPU-accelerated voronoi diagrams [SGG * 06], Layered Depth Images [START_REF] Heidelberger | Real-time volumetric intersections of deforming objects[END_REF]. Deformed distance fields [START_REF] Fisher | Fast penetration depth estimation for elastic bodies using deformed distance fields[END_REF][START_REF] Marchal | Collision between deformable objects using fastmarching on tetrahedral models[END_REF] have been used, however, they work for relatively small deformations. Most of the overhead of these techniques resides in updating the volumetric representation.

While contact response between rigid bodies is generally based on constraints [MW88, Hah88, Bar93, MC95], contact response for soft bodies can be based on penalty forces [START_REF] Baraff | Large steps in cloth simulation[END_REF][START_REF] Volino | Implementing fast cloth simulation with collision response[END_REF] or constraints [START_REF] Provot | Collision and self-collision handling in cloth model dedicated to design garments[END_REF]. [START_REF] Bridson R | Robust treatment of collisions, contact and friction for cloth animation[END_REF] combines penalties with constraints in a collision processing framework independent of the internal dynamics of simulated cloth.

Interactions with fluids are very challenging, due to the variety of representations and solvers used. While most solid simulations compute the positions and forces of a set of Lagrangian points, fluid solvers handle velocity and pressure fields defined over a fixed Eulerian grid [START_REF] Stam | Stable Fluids[END_REF] or a set of Lagrangian particles [SAC * 99, MCG03, PTB * 03]. If the fluid is not present on the whole domain, its surface needs to be tracked using marker particles [START_REF] Foster | Realistic animation of liquids[END_REF] or levelsets [START_REF] Foster | Practical animation of liquids[END_REF]. Recent approaches include combining particles and levelsets [EMF02, ELF05], using an octree representation [START_REF] Losasso | Simulating water and smoke with an octree data structure[END_REF][START_REF]Discontinuous fluids[END_REF], and using a polygonal representation to enhance precision near the interface [START_REF] Bargteil | A semi-lagrangian contouring method for fluid simulation[END_REF].

A simple approach to compute two-ways interactions between fluids and other objects is to consider them as fluids, either by adding particles at the surface of solids in contact with a particlebased fluid [GHD03, MST * 04], or by considering rigid bodies as fluids with an additional rigidity constraint [START_REF] Carlson | Rigid fluid: animating the interplay between rigid bodies and fluid[END_REF]. More accurately, fluid-structure interactions are physically expressed by integrating the fluid pressure over the interface to provide the force applied to the solid, and using the solid's velocity as boundary conditions for the fluid [START_REF] Benson | Computational methods in lagrangian and eulerian hydrocodes[END_REF][START_REF] Guendelman | Coupling water and smoke to thin deformable and rigid shells[END_REF]. , or to accurately compute the coupling by remeshing around immersed objects [START_REF] Klingner | Fluid animation with dynamic meshes[END_REF][START_REF] Chentanez | Simultaneous coupling of fluids and deformable bodies[END_REF].

While not strictly related to fluid-solid interactions, the recent work of Losasso et al. [START_REF] Losasso | Multiple interacting liquids[END_REF] presents a new, impressive method for handling multiple interacting fluids. Unfortunately it is not directly applicable to model interactions with other objects.

Background

In this section, we discuss the main approaches to model contact and apply reactions and show that our density-based contact model overcomes the drawbacks of the distance fields.

Contact modeling approaches

There are different classes of contact modeling approaches. The first is based on surface proximity. In this approach, contact is modeled as pairs of geometrical primitives (such as vertices, edges, and polygons) in close proximity and undergoing repulsive forces or constraints. However, primitives which cross the opposite surface are then repelled in the wrong direction. Consequently, a distance offset is typically needed to avoid intersections, and robustness requires adaptive time steps, continuous-time collision detection, and sophisticated mixes of elastic forces and hard constraints [START_REF] Bridson R | Robust treatment of collisions, contact and friction for cloth animation[END_REF].

The second class is based on penetration depth. In this approach, sample points on the surface of one object are tested against the volume of the other, as illustrated in the left of figure 2. Contact is modeled as pairs of points and targets, such as (A, A ) in figure 2. This approach is more robust because the primitives which cross the opposite surface are then repelled outward, which allows intersection. However, it can not be applied to pairs of lines and surfaces because they have no volume. Moreover, the directions of repulsion can be poorly approximated, like pair (B, B ) in figure 2, or unstable, like pair (C,C ) which depends on which side of the bisecting line of the square corner point C is located. The target points are the projections of the sample points to the closest primitive of the opposite surface. They are easily computed when simple volumes such as spheres or tetrahedras [THM * 03] are used. Generally, finding targets involves the expensive computation and the storage of distance fields. Precomputation can be used for rigid bodies [START_REF] Guendelman | Nonconvex rigid bodies with stacking[END_REF], but sophisticated update strategies are necessary in the case of elastic bodies [START_REF] Fisher | Fast penetration depth estimation for elastic bodies using deformed distance fields[END_REF] and fluids [START_REF] Enright | Animation and rendering of complex water surfaces[END_REF].

Recently, an alternative approach based on the minimization of the intersection volume has been presented by [START_REF] Volino | Resolving surface collisions through intersection contour minimization[END_REF]. They model the intersection volumes using the intersection lines and apply their technique to cloth. As far as we know, this approach has not been applied to volumetric objects due to the difficulty of efficiently modeling the intersection volumes. Our method overcomes this problem using an Eulerian density grid, as illustrated in the right of figure 2, and explained in section 4.

Contact reaction

Once contact is modeled, we have to compute how to move the points along the contact directions. The reaction methods basically fall in two categories, or a mix of both. Force-based methods apply elastic attraction between the contact points and their targets. Intersection occurs, and stiff springs are often necessary to keep it small enough. This requires either very small time steps to avoid instabilities, or the use of implicit time integration. Implicit integration uses the partial derivative of the forces with repect to the positions (i.e. stiffness). A linear equation system with a positive symmetric definite (PSD) matrix must be solved at each time step. Such equation systems have a unique solution and can be solved iteratively and efficiently using the conjugrate gradient algorithm [START_REF] Baraff | Large steps in cloth simulation[END_REF]. The advantages of this approach are the stability, the numerical efficiency, and the good control of INRIA the computation time. However, handling Coulomb friction is not straighforward. Moreover, implicit integration is possible only when contact forces are differentiable, which prohibits the use of distance fields due to the instability of the target directions (point C). Consequently, distance fields are generally used with contraint-based reaction methods.

Constraint-based methods aim to totally avoid intersections, by requiring the contact points to meet their targets. This generally results in singular systems which may have an infinite number of solutions as well as no solution at all, and sophisticated LCP solutions may be necessary [START_REF] Baraff | Fast contact force computation for nonpenetrating rigid bodies[END_REF], with poor control of computation time. A popular alternative to solve such systems is to apply Gauss-Seidel iterations, where the constraints are processed sequentially until the residual falls under a given threshold, or a maximum number of computations is reached. The advantages of this approach are simplicity, versatility, easy Coulomb friction handling, and the control of the computation time. The drawbacks are a very slow converging rate when pressure waves propagate through numerous objects or through objects with highly different masses, as well as oscillations in case of conflicting constraints. Bridson et al [START_REF] Bridson R | Robust treatment of collisions, contact and friction for cloth animation[END_REF] present sophisticated methods to control proximities as well as strain using mixed force-and constraint-based strategies.

Force-and constraint-based reaction methods have different advantages and drawbacks, but they all rely on the target directions computed by the contact modeling method. Our new density-based contact modeling method has the same versatility as the distance fields, while overcoming its main drawbacks. Good, stable target directions are computed, allowing the implementation of robust elastic repulsion forces using implicit integration.

Eulerian contact

This section discusses the middle and bottom layers of figure 1. We show how to detect contact and how to set up the contact equations for force-and constraint-based responses.

Density fields

We call spatial density (or density, for short) the ratio of space locally occupied by matter. Density is a discontinuous function with value 1 inside a body, 0 outside, and 0.5 at the surface. However, we can get a continuous approximation of its value using finite differences in a regular grid. This creates as contact volume around the surface with the thickness of one step of the density grid. Near the surface, the density gradient gives an approximation of the normal, and bodies can repel each other along this direction. Though sharp edges can not be modeled exactly using density fields, we can get arbitrarily close approximations of them using appropriately fine sampling. Density fields can actually be seen as volumetric images of the bodies, and all the classical greymap processing techniques can be applied to control the level of detail. Figure 3 shows that smooth density fields can be modeled even with large downscaling factors. More sophisticated compression techniques such as wavelet bases could be used to increase computational efficiency. A key point of our mechanical framework is the very simple form of the collision detection criterion used at the bottom layer in figure 1: ∑ i ρ i > 1 where index i ranges over all the bodies present at a given point. In the remainder of this document, we only discuss pair-wise interactions. Interactions between more than two bodies at the same place are processed as pairs of interacting bodies.

The target directions are parallel to the gradient of density, as shown in sections 4.2 and 4.3. The gradient is null inside the collision volume (see point E in figure 2), and parallel to the normal near the surface (point D) or to a weighted sum of the normals (point F). The contact force is thus proportional to the external area of the contact volume, which is consistent with real-world physics. The density field can be sampled efficiently at each time step, as explained in section 5. Lines and surfaces are modeled as volumetric objects with thickness equal to the sampling step of the grid.

Contact penalties

When two bodies intersect, we can model an elastic penalty derived from a potential energy related to intersection. Using

ρ 1 + ρ 2 > 1 : ρ = ρ 1 + ρ 2 -1 else : ρ = 0
we define E = 1 2 kV ρ2 as potential energy where V is the volume associated with contact point as a function of the grid resolution, and k is a user-defined parameter related to material properties. The following elastic forces at a given intersection point r can be derived from the potential

f 1 = - ∂ E ∂ r = -kV ρ∇ρ 1 , f 2 = -∂ E ∂ r = -kV ρ∇ρ 2 . (1)
The density gradient ∇ρ = ∂ ρ ∂ r is computed using finite differences in the regular grid. We use Eulerian point r as a contact point between the bodies. Newton's law f 1 =f 2 requires that ∇ρ 1 = -∇ρ 2 which is only true when the bodies are tangent. In practical simulations, the bodies intersect and the density gradients have different directions and intensities.

To maintain physical consistency, we have to average the forces computed using equations 1

f 1 = -f 2 = - kV 2 ρ(∇ρ 1 -∇ρ 2 ) . (2)

INRIA

A force orthogonal to the repulsion direction can be applied to model friction.

Note that these contact forces are only applied near the boundaries of the objects, where ∇ρ is not null. The net force is thus a function of the intersection area rather than intersection volume. When desired, it is possible to obtain volume-based penalties by decreasing the resolution of the contact grid.

Stiff penalties are necessary to avoid deep intersections. To maintain numerical stability, we have either to apply small time steps, or to perform implicit time integration [START_REF] Baraff | Large steps in cloth simulation[END_REF]. The latter uses the stiffness operator, which relates a variation of force to a variation of body positions. The variations of force at points r and r + ∆r due to the displacement can be computed by differentiating the Eulerian force (eq. 2) as detailed in appendix A:

∂ f ∂ r = kV 4 (∇ρ 2 -∇ρ 1 )(∇ρ 2 -∇ρ 1 ) T + ρ ∂ ∇ρ 1 ∂ r + ∂ ∇ρ 2 ∂ r (3)
which is a symmetric matrix.

Contact constraints

When constraint ρ 1 + ρ 2 ≤ 1 is violated, we can solve it approximately using a first-order differentiation of the equation:

d(ρ 1 + ρ 2 ) dr ∆r = - ρ
which gives us the scalar equation

∇ρ T 1 ∇ρ T 2 ∆r 1 ∆r 2 = - ρ
where ∆r 1 and ∆r 2 are two displacement for a single point in the grid for the two bodies sampled at this point. This equation has 6 unknown scalars. We can constrain the displacement along the average gradient direction n given by n = ∇ρ 1 -∇ρ 2 ∇ρ 1 -∇ρ 2 which reduces the equation to

∇ρ T 1 n ∇ρ T 2 n a 1 a 2 = - ρ
which has 2 unknown scalars. Newton's laws give us the additional equation m 1 a 1 + m 2 a 2 = 0 where m 1 and m 2 are the masses of the volume elements located at the grid point. This allows us to compute displacements at the grid particles. Continuous-time collision constraints can be modeled similarly. Note that due to finite differences, the computed displacement can not be larger than the Eulerian grid step. Resolving a deeper intersection would require several iterations.

Mappings

In the previous sections, we have shown that non-intersection penalties and constraints can be expressed using simple formulas based on ρ and ∇ρ at each contact point. In order to use these relations for the mechanical coupling of the intersecting bodies, we need to write these expressions as functions of the mechanical degrees of freedom (DOF) of the bodies. This section discusses the mapping between the top and middle layers in figure 1.

Mechanical requirements

Let q and q be the values and the time derivatives of the independent (also called generalized, or Lagrangian) DOF of the bodies. The associated forces are denoted using Q. In this section, we show which mathematical operators are necessary to map the DOF of the bodies to the Eulerian layer. The derivation of these operators in various cases is presented in section 5. The interaction forces are computed at grid points, and we need to map these forces to generalized forces applied to the bodies. To do this in a physically consistent way, we have to ensure that the virtual power of the contact forces f is equal to the virtual power of their corresponding generalized forces Q. The virtual power of a contact force acting on a body is the dot product of the force by the virtual velocity ṙ of a particle attached to the body and located at the contact point. The velocities at the contact points are given by ṙ

= ∂ r ∂ q q .
We thus need to compute Q so that f T ∂ r ∂ q q = Q T q , where ∂ r ∂ q is the kinematic matrix relating the generalized velocities of the body to the velocity of the body particle located at the contact point. Since this relation must hold for any set q of generalized velocities, we get at each contact point

Q = ∂ r ∂ q T f .
When implicit time integration is needed, The stiffness matrix K of the contact force at each contact point is necessary. This operator encodes the variation of force in response to a small variation of q:

K i j = ∂ Q i ∂ q j = ∂ ∂ q j ∂ r ∂ q j T f = ∂ 2 r ∂ q 2 j T f + ∂ r ∂ q j T ∂ f ∂ r ∂ r ∂ q j
where ∂ f ∂ r is given in equation 3. The kinematic operators ∂ r ∂ q j and ∂ 2 r ∂ q 2 j make the connection between the contact layer and the bodies. The application of hard constraints requires the computation of dq given dr. Since the velocity field at a given point can depend on an arbitrary number of mechanical variables, the pseudo-inverse of ∂ r ∂ q i is required. The same is used to apply a velocity change. This allows the application of INRIA instantaneous displacement or velocity changes, as commonly used in constraint-based methods. As usual, conflicting constraints can be solved iteratively using Gauss-Seidel approaches. In summary, the necessary operators for the bodies to interact in the Eulerian layer are only ρ(r), ∂ r ∂ q and ∂ 2 r ∂ q 2 . The first is the simplest geometric definition of a volume, while the second is the kinematic matrix relating particle velocities ṙ to the generalized velocities q, and the third encodes the non-linearity of r(q). The other operators can be deduced using finite differences in the regular grid. In the remainder of this section, we explain how to build these operators for the most commonly used models.

Simple geometric primitives

Simple geometric primitives such as points, lines, triangles, tetrahedra and spheres can be sampled directly using techniques inspired from classical graphics rasterization. The sampling algorithm loops over the target points to evaluate the density. Points are sampled directly in the grid using antialiasing techniques. Lines can be sampled using a 3D variant of Bresenham's algorithm. Similarly, triangles and tetrahedra can be sampled by looping over the regularly spaced target points r using the appropriate variants of the traditional polygon-filling algorithms. Density ρ(r) may be binary or real. In a translating sphere, each particle moves with the center and the mapping is the identity. In the other primitives, the velocities ṙ are linearly interpolated like colors in graphical rendering. The interpolation coefficients are the non-null entries of the kinematic matrix ∂ r ∂ q . The second-order

∂ 2 r ∂ q 2
is null due to linearity.

Rigid bodies

When a shape is rigidly attached to a moving reference frame, the kinematic matrix encodes the standard rigid body velocity field ṙ = ȯ + ω × (ro), where o is the origin of the reference frame and ω is the angular velocity associated with a given kinematic variable. Its derivative ∂ 2 r ∂ q 2 = ω × (ω × (ro)) encodes the curvature of the trajectories. Density can be evaluated by projecting back the Eulerian points to a precomputed distance or density field, or by projecting raw data to the Eulerian grid as explained in the next section.

Complex shapes

In the case of volumes defined by raw data such as medical images, or embedded in nonlinear deformation fields, we propose to use an additional layer, as illustrated in figure 4. We call it the volume layer. It represents the volume occupied by the body, using sampling particles attached to the body, with a volume assigned to each of them. The particles of this layer are not independent, since their positions and velocities q , q are entirely defined by the particles of the top layer. The designer of the physical model is responsible for implementing the mapping to the volume layer. Fortunately, implementing the mapping associated with a complex-shaped physical body is as simple as depicting it as a set of spheres with centers, velocities (along with the kinematic operators discussed in section 5.1), radii and volumetric masses. INRIA However, the resulting contact behavior is different from what can be obtained using traditional sphere-based contact models, because in our framework the contact spheres are then mapped to the Eulerian layer. This results in a discrete number of possible contact points, and much smoother contact surfaces, as shown in the accompanying video.

The mapping from the volume layer to the Eulerian layer is computed automatically as follows. Each particle distributes its density to the closest points in the Eulerian grid, and we average the contributions of several particles in the same Eulerian voxel. Let us consider a moving volume particle with position q i , velocity q i and mass m i . The masses are used to compute physically consistent linear mapping coefficients. They can be straightforwardly deduced from the volume and the volumetric mass of the body. Physical consistency is achieved by distributing the mass and momentum of this volume particle to the corners of the Eulerian grid cell which contains it using its barycentric coordinates. The net mass and momentum distributed to an Eulerian grid point are the sum of the contributions of all the volume particles contained in its neighboring cells:

m j = ∑ i∈body particles w i j m i , (4)
(mṙ) j = ∑ i∈body particles w i j m i q i (5)
where w i j is the barycentric coordinate of volume particle i with respect to Eulerian point j. The density associated with a volume particle is its volume divided by the volume of the Eulerian cell.

Once the mass m j and momentum (mṙ) j associated with an Eulerian grid point are computed, its associated velocity and kinematic mapping coefficients are

ṙ j = (mṙ) j m j , ∂ r j ∂ q i = m i m j , (6)
The operator ∂ 2 r ∂ q 2 is null due to linearity. Note that the radius of influence of the particles in the volume layer can be larger than one cell of the Eulerian grid. In this case, the density, the mass and the velocities are distributed over more than one cell.

The volume layer can be used to depict the surface of meshless deformable bodies. In this model, radial basis functions are attached to moving particles, the weighted sum of the particle displacements defines the kinematic field of the body, which is used to model the internal deformations, which is in turn used to compute appropriate response forces. This is a special case where the degrees of freedom q of the volume layer are the same as the independent degrees of freedom q.

Implementation

While the fundamental equations have been discussed, implementing them efficiently requires several important considerations, detailed bellow.

While mathematically we can consider computing the density of all objects everywhere, in order to obtain good performances it is important to compute it only when necessary. In order to achieve this we construct a hierarchical representation of the Eulerian grid. Whenever possible, only the root level of the hierarchy is computed, and finer levels are created on demand when more than one object is present. This strategy prunes the computation of most densities, except close to potential contacts. However in some cases the object is able to rasterize itself very efficiently, such as fluids internally computing a levelset, in which case it is more efficient to construct the hierarchy bottom-up.

In the previous sections, we have shown how potentially all physical models can exert interaction forces on each other through our Eulerian contact layer. We can thus consider the interactions as given external forces and animate each body independently using its dedicated differential equation solver. We have also designed an implicit time integration solver able to animate the interacting bodies when stiff forces and large time steps are applied. Consistently with our contact processing approach, the solver has no direct access to the internal degrees of freedom of the bodies. Each body stores internally its own state and auxiliary vectors including position, velocity, force and displacements. The solver uses the conjugate gradient algorithm, which triggers force computation and vector operations without addressing the vectors directly. The only feedback from the bodies to the solver is the value of dot products between auxiliary vectors. The computation of the net force applied to each body as a function of its positions and velocities is done recursively through the mappings. Stiffness computation is detailed in appendix B.

This approach makes no assumption about the type of degrees of freedom of the interacting bodies. One could be a viscoelastic body while the other could be a hierarchy of articulated rigid bodies, provided that the interaction forces are applied consistently. The method is reusable for all available physical bodies. Moreover, simple constraints such as maintaining an independent particle at a fixed place are internally applied by the bodies as filters. Notice that this flexible solver is not specialized to our contact framework and can be applied to traditional systems of bodies in contact. It allows us to straightforwardly apply [START_REF] Baraff | Large steps in cloth simulation[END_REF]'s implicit time integration using a filtered conjugate gradient solution to our layered body framework, and we have been using it successfully.

We have also designed a generic constraint solver using the formulas given in section 4.3. Constraints are processed one after another, and each correction is immediately propagated, in a Gauss-Seidel manner. In our experiments, the stiff penalty approach gives more stable results. However, we plan to use it in future work to apply Coulomb friction using the robust method proposed by [START_REF] Bridson R | Robust treatment of collisions, contact and friction for cloth animation[END_REF].

Results

Our new framework allows us to animate scenes including arbitrary physical models. We implemented and tested many algorithms, including rigid bodies, deformable objects (using either masssprings or FEM), an Eulerian grid-based fluid solver, as well as a Lagrangian SPH fluid.

To test the robustness of our method, we designed several simulations involving up to three different types of objects. Each algorithm is only completed by the mappings to and from the Eulerian density grid, and has no knowledge nor special codes to handle interactions to another specific type of objects. Despite this lack of specialized interactions, we were able to reproduce previous experiments, such as collisions between a fluid and a thin cloth [START_REF] Guendelman | Coupling water and smoke to thin deformable and rigid shells[END_REF] (figure 5).

While interactions between fluids, rigids and deformables have previously been demonstrated, coupling all three in a single simulation has not yet been clearly demonstrated. This is the motivation behind the animation shown in figure 5(down), involving a light poolrack floating over a fluid and crushed by an heavy rigid object falling to the ground. In this case the fluid solver is using a very different integration method than the rest of the simulation (semi-lagrangian advection, then pressure projection), hence the fluid is not included in the implicit integration system and interacts with the other bodies using forces integrated explicitly. While this requires a rather small timestep (of the order of 1/100 th second), it is nonetheless able to produce correct interactions. Expressing forces interacting with the fluid implicitly, as shown by Chenatez et al. [START_REF] Chentanez | Simultaneous coupling of fluids and deformable bodies[END_REF], would allow for much larger timesteps.

The previously mentioned simulations were executed on a single processor at speeds varying between a few frames per seconds to less than 5 seconds per frame. More importantly, this speed was not significantly slower than each individual simulated object taken individually. In many cases, the fluid solver is the bottleneck. The scene shown in the second picture in the first page includes rigid and deformable volumetric bodies, as well as pieces of cloth. The total amount of particles is more than 21,000. As shown in 6, the collision detection and modeling takes less than 25% of the computation time, even in this complex stacking configuration.

Discussion

Our framework meets the objective of allowing contact between arbitrary physical models. Contrary to all previous work, our approach allows contact interactions without any assumption about object geometry (e.g. mesh, implicit function) or constitutive laws (e.g. rigid, visco-elastic, fluid). We demonstrate it using rigid, visco-elastic and fluid bodies with mesh-based and implicit geometries. The differentiability of the density-based contact forces allows the use of implicit time integration, resulting in stable simulations. Using appropriate penalty stiffness, contact errors are easily smaller than the size of one density voxel. Even in the complex case of interaction between fluid and deformable bodies, the important phenomena such as buoyancy and the meeting of boundary conditions are present.

Our method is more efficient than distance fields. Density fields are computed much more efficiently because density is a local, additive property. This allows us to simulate deformable objects, while methods based on distance fields are generally limited to rigid bodies with precomputed distance fields, or, in the case of fluid simulation, involve complex interface tracking methods.

We believe that the computations in the Eulerian grid are very well suited to the hardware accelerations available with modern graphics architectures.

The limitations of our method are mainly due to density rasterization. The geometrical accuracy of our method directly depends on the resolution of the grid. Nonetheless, as shown in the first accompanying video, a very coarse grid gives surprisingly good results, and a finer grid makes intersections unnoticeable. Aliasing artifacts are hardly noticeable due to smooth density mapping. Note that in case of deformable bodies, contact precision is typically less of an issue, allowing a trade-off between accuracy and speed when needed.

Our current implementation has the limitation that the grid resolution must fit the smallest object in the scene. This results in non-interactive frame rates when thin objects such as cloth are involved. In future work, we plan to use hierarchical grids, and interpolate the coarse grids locally near the contacts.

In future work, there are a number of issues that we would like to address. Self-collision could be handled by partitioning an object's shape in several complementary sub-shapes tagged with different identifiers but referring to the same object. When a solid object is wholly contained within another solid object, the object inside will not be repelled because the forces acting on its surface will cancel each other. To handle this, time continuous collision detection could be used. Locally, we can estimate the time derivative of the density based on density gradients and velocities. When distant objects with high velocities are present, we could fill the voxels swept during the whole estimated INRIA trajectory rather than only the current location, similarly with the expanded bounding boxes used by some authors.

Our approach does not compete with specialized methods when only one kind of object is involved. The benefit of our method is to allow any kinds of physical bodies to interact.

A Contact stiffness

When a body moves, its density field moves accordingly. We start from equation 2 and consider a small displacement ∆r of body 1, while body 2 remains fixed. Density ρ 1 (r) undergoes a change ∆ρ 1 -∇ρ 1 ∆r resulting in a new contact force at the starting point r:

f (r) = kV 2 ( ρ -∇ρ T 1 ∆r)(∇ρ 2 -∇ρ 1 + ∂ 2 ρ 1 ∂ r 2 ∆r)
The first-order variation of force is thus:

∆ f (
Considering an infinitely small ∆r and averaging the variations at the starting point (eq. 7) and at the end point of the displacement (eq. 8), we obtain the symmetric stiffness matrix of equation 3. The same result is obtained by considering a small displacement of body 2 while body 1 remains fixed.

B The stiffness of bodies in contact

We write the net generalized force as a function of the independent degrees of freedom as:

Q(q) = Q internal (q) + ∂ r ∂ q T f
where Q internal is the internal force and f the contact force. This can be computed in three steps The variation of the generalized forces due to a variation of q is:

∆Q = ∂ Q internal ∂ q ∆q + ∂ r ∂ q T ∂ f ∂ r ∂ r ∂ q ∆q + ∇ q ∂ r ∂ q T ( f )
where ∇ q ∂ r ∂ q T is a function of the geometric properties of the mapping r(q) and computed by the mapping. Given ∆q, the first term is computed internally by the bodies. The second term is computed in three steps, by first propagating ∆q to the Eulerian layer through the mapping, then applying contact stiffness, and finally mapping the result back to the bodies. The last term is computed by the mappings.

Figure 2 :

 2 Figure 2: Contact modeling. Left: distance field. Right: density field.

Figure 3 :

 3 Figure 3: Density Fields. The red line shows the 0.5 isovalue. (a): Original 256x256 image. (b),(c),(d): Same field reduced to the resolution of 16x16 and then rescaled to original resolution using subsampling, linear interpolation, and bi-cubic interpolation, respectively.

Figure 4 :

 4 Figure 4: The volume layer inserted between the Lagrangian and Eulerian layers. Left: deformable body with nonlinear deformation, right: rigid body with raw volume data.

Figure 5 :

 5 Figure 5: Top: contacts involving an SPH-based fluid and a thin deformable cloth. Down: Interactions between an Eulerian grid-based fluid, a co-rotational FEM deformable model, and a rigid body.

Figure 6 :

 6 Figure 6: Performance measurements. The white upper area represents the time devoted to the computation and integration of internal and repulsion forces. The rest corresponds to collision detection and modeling.

  At the same time, the force at the end point r + ∆r becomes:

		r)	kV 2	-∇ρ T 1 ∆r(∇ρ 2 -∇ρ 1 ) +	ρ ∂ 2 ρ 1 ∂ r 2 ∆r
	and we can write:	∆ f (r) ∆r		kV 2	(-∇ρ T 1 (∇ρ 2 -∇ρ 1 ) +	ρ ∂ ∇ρ 1 ∂ r	)	(7)
		f (r + ∆r) =	kV 2	( ρ + ∇ρ T 2 ∆r)(∇ρ 2 +	∂ 2 ρ 2 ∂ r 2 ∆r -∇ρ 1 )
	and we obtain:	∆ f (r) ∆r	kV 2	(∇ρ T 2 (∇ρ 2 -∇ρ 1 ) +	ρ ∂ ∇ρ 2 ∂ r	)

RR n°6203

Faure & Allard & Nesme

http://www.inria.fr ISSN 0249-6399