Pascal Volino

Nadia Magnenat-Thalmann

Miralab

Francois Faure Ljk

Simple, yet Accurate Nonlinear Tensile Stiffness

Keywords: I.6 [Computing Methodologies]: Simulation and Modeling, I.3.5 [Computing Methodologies]: Computer Graphics -Computational Geometry and Object Modeling -Physically-Based Modeling, J.6 [Computer Applications]: Computer-Aided Engineering. General Terms: Mechanical Simulation, Cloth Simulation Particle Systems, Finite Elements

Recent Particle System models have evolved toward accurate representation of elastic stiffness based on continuum mechanics, converging to formulations that make them quite analogous to fast Finite Element methods. These formulations usually involve the linearization of tensors that help their formulation in the context of linear elasticity. Toward our objective of simulating the nonlinear properties of cloth accurately, we show through this work that this linearization can indeed be suppressed and replaced by adapted strain-stress laws relating precisely the nonlinear behavior of the material. This leads to very streamlined computations that are particularly efficient for simulating the nonlinear anisotropic tensile elasticity of deformable surfaces. Through a simple and efficient implementation using the Particle System formalism, we demonstrate the efficiency of this method with examples related to garment simulation.

INTRODUCTION

The mechanical properties of cloth materials are highly isotropic and nonlinear: The internal forces in the material are not at all proportional to the deformations, and they furthermore vary much with their orientation relatively to the thread directions. This anisotropy and nonlinearity still creates significant challenges when it comes to define a mechanical simulation system that would accurately reproduce these effects on virtual objects, which would also be applicable for actual applications such as accurate CAD systems (Figure 1). At the meantime, interactive and Virtual Reality applications require very efficient simulation models that compute quickly enough for offering good reactivity to user interaction.

Combining these requirements is still one of the main challenges of cloth simulation, and we will contribute to address this by showing that nonlinearity is indeed not necessarily synonym of complex implementation and slow computation, and that the way-too-often unconsidered hyperelastic St.Venant-Kirchhoff materials can indeed be key actors in accurate interactive cloth simulation.

Existing Simulation Schemes

Particle Systems have always been of large interest in the field of cloth simulation, and more generally in the field of interactive mechanical simulation, as they offer a simple, intuitive and flexible way to model mechanical systems. Furthermore, they can be combined with a large range of numerical integration schemes, according to the important features of the simulation context (dynamic accuracy, convergence speed, fast and approximate simulation, robustness...).

Early Particle Systems were grid-based [START_REF] Breen | Predicting the Drape of Woven Cloth using Interacting Particles[END_REF] [START_REF] Eberhardt | A Fast, Flexible, Particle-System Model for Cloth Draping, Computer Graphics in Textiles and Apparel[END_REF], and already featured simulation of nonlinear behavior curves through formulations that made them quite analogous to continuummechanics models. They were however fairly limited accuracy for large deformations, and required quite long computation times. Faster models, based on spring-mass grids, have once again become popular since implicit integration methods have started to be used for Particle Systems [Baraff and Witkin, 1998], because they allowed a simple expression of the Jacobian of the particle forces while requiring only simple computations [START_REF] Desbrun | Interactive Animation of Structured Deformable Objects[END_REF]] [START_REF] Meyer | Interactive Animation of Cloth-like Objects in Virtual Reality[END_REF] [START_REF] Choi | Stable but Responsive Cloth, Computer Graphics[END_REF]. Combined with advanced implicit integration methods [START_REF] Eberhardt | Implicit-Explicit Schemes for Fast Animation with Particles Systems[END_REF] [START_REF] Hauth | A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods[END_REF] [Volino and Magnenat-Thalmann, 2005], these simulation schemes have particularly become popular for real-time and interactive applications. Unfortunately, spring-mass systems are quite unable to model surface elasticity accurately. Although some techniques have been developed to identify their parameters to those of the simulated material, they do not allow full discrimination between deformation modes [START_REF] Bianchi | Simultaneous Topology and Stiffness Identification for Mass-Spring Models Based on FEM Reference Deformations[END_REF], and they remain particularly inaccurate for anisotropic and nonlinear models. Particle Systems, as a whole, have inherited this reputation of inaccuracy from them.

On the other hand, Finite Elements have now acquired a good maturity for mechanical simulation. Their traditional field of application is elastic solid or shell modeling for mechanical engineering purposes, a context where linear elasticity and small deformations are the rules. These formulations are not so well adapted to very deformable objects such as cloth, and early attempts to model cloth using high-order elements [START_REF] Eischen | Finite Element Modeling and Control of Flexible Fabric Parts, Computer Graphics in Textiles and Apparel[END_REF] led to impractically high computation times. Some recent developments have attempted to speed up the computation times required for Finite Elements.

These have particularly been used for simulating deformable volumes in the context of interactive simulation for virtual surgery systems [START_REF] Debunne | Dynamic Real-Time Deformations Using Space & Time Adaptive Sampling[END_REF] [START_REF] Hauth | Interactive Physically-Based Solid Dynamics[END_REF]. For instance, pre-inverting the linear system matrix (as done by [START_REF] Desbrun | Interactive Animation of Structured Deformable Objects[END_REF] for Particle Systems) may speed up the computation [Bro-Nielsen and [START_REF] Bro-Nielsen | Real-Time Volumetric Deformable Models for Surgery Simulation using Finite Elements and Condensation[END_REF], but restricts the application field to linear models and very small mechanical systems [START_REF] James | ArtDefo -Accurate Real-Time Deformable Objects[END_REF]. "Explicit Finite Elements" [O' Brien and Hodgins, 1999] come close to good computational charge compromise by locally approximating the solution of each element [START_REF] Cotin | Real-Time Elastic Deformations of Soft Tissues for Surgery Simulation[END_REF]. However, these models rely on simple linear elasticity which is inappropriate for the large deformations or displacements encountered in cloth simulation. Alleviating this problem through the use of a linearized form of the Green Lagrange strain tensor has recently led to more accurate models [START_REF] Debunne | Dynamic Real-Time Deformations Using Space & Time Adaptive Sampling[END_REF]] [Etzmuss et al, 2003]. At the meantime, Particle Systems have also evolved toward accuracy, and there is indeed a significant convergence in the formulation of accurate Particle System models based on continuum mechanics [Etzmuss et al, 2003] with their Finite Element counterparts, as both describe accurately surface elasticity through representations of the Green-Lagrange strain tensor with various degrees of linearization, and these intermediate approaches allow very efficient simulation [START_REF] Teschner | A Versatile and Robust Model for Geometrically Complex Deformable Solids[END_REF].

One major drawback remains with these accurate schemes: They are only applicable in the context of linear elasticity. In most models, this limitation essentially results from the linearization process commonly used to turn the Green-Lagrange tensor into a linear Cauchy tensor. In order to preserve accurately for large deformations, this is carried out by rotating the coordinates system of parameter space of the elements so as to match the current eigendirections of the strain tensor, expressing it as a shear-less deformation [START_REF] Hauth | Corotational Simulation of Deformable Solids[END_REF][START_REF] Muller | Interactive Virtual Materials, Proceedings of Graphics Interface[END_REF][START_REF] Muller | Stable Real-Time Deformations[END_REF]] [Etzmuss et al, 2003[START_REF] Nesme | Efficient, Physically Plausible Finite Elements[END_REF]. Besides significant additional computation, these corotational schemes make it complicated to express anisotropic mechanical behaviors, particularly if they are nonlinear. Attempts to avoid the corotational formulation [Volino and Thalmann, 2005] also lead to inaccurate approximations of the shear strain and stress for large deformations.

St.Venant-Kirchhoff Materials

Indeed, there is a way to simply describe accurately the strain and stress of materials: Using the Green-Lagrange tensor without any linearization. This is the basis of hyperelastic and St.Venant-Kirchhoff materials [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF][START_REF] Zhuang | Haptic Interaction with Global Deformations[END_REF]] [Picinbobo et al, 2003] [ [START_REF] Barbic | Real-Time Subspace Integration for St.Venant-Kirchhoff Deformable Models[END_REF]]. This approach is less "intuitive", through the nonlinearity of their strain and stress tensors. Hence, their strain is not proportional to the deformation of the material, and the exerted force is not proportional to the stress. Actually, with stress linearly related to strain, the tensile force-deformation curve of such a material is rather cubic (Fig. 2). However, their mathematical definition is indeed the most mathematically "natural" way of expressing strain and stress, and indeed the simplest, despite the nonlinearity. However, the nonlinear native behavior of St.Venant-Kirchhoff materials has a significant adverse effect: Under large compression deformation, the tensile compression forces stop increasing, and even revert back to zero (Fig. 2 left). Hence, these materials are prone to "collapse" under large compression, which makes them unsuitable for volume simulation [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]Wood, 1997] [Irving et al, 2004]. This is indeed why these materials are usually overlooked in mechanical simulation systems, preferring linearized formulations instead.

Fortunately, this collapse behavior is usually not an issue in the context of cloth simulation: A cloth surface is almost never subject to large compression forces, since weak surface bending forces quickly allow compression to relax into buckling. Furthermore, the rapidly increasing tensile forces during extension are a much more realistic representation of the behavior of real fabric materials than a simple linear behavior.

Along this, the most attracting perspective of these materials is the simplicity of their formulation, ensuring efficient computation. Furthermore, they can easily be generalized for supporting arbitrary nonlinear strain-stress laws, for reproducing accurately the complex behavior of cloth materials through the weft, warp and shear deformation modes, a feature that most cloth simulation models lack, and that we intend to provide through this work.

Our Goal

Our ultimate goal is to be able to simulate very accurately complex cloth objects, such as complete garments on animated characters, with precise reproduction of the nonlinear mechanical behavior of cloth. Unlike what is possible to do with most existing cloth simulation models, we intend to reproduce accurately the nonlinear and anisotropic strain-stress behaviors, as discussed in Section 2.1, so as to match as closely as possible the actual behavior of cloth materials, particularly when it comes to large deformations. Our model should address elasticity as well as viscosity, making the model suitable not only for draping applications, but also for dynamic motion computations which require mechanical damping. We also want to perform the computation at a reasonable speed, possibly compatible with real-time or interactivity. Hence, we cannot afford models much more complicated than the simple spring-mass models which are typically used in this context.

Our new computation scheme should also support arbitrary triangle meshes which are typically generated from Delaunay triangulation. On these elements, we express strains and stresses according to the nonlinear Green-Lagrange tensors, also considering nonlinear strain-stress behaviors for representing accurately the nonlinear anisotropic behavior of cloth. Indeed, such use of continuum mechanics makes our scheme equivalent to classic tensile Finite Element schemes. However, we take advantage of its simplicity for formulating it in a very simple and direct Particle System approach, expressing particle forces directly from their positions. Furthermore, we formulate an accurate computation of the Jacobian, allowing high-performance simulation through efficient implicit integration methods or relaxation schemes which remain stable and robust even with very large deformations. All these developments are expressed without excessive abstract formalisms through "ready-to-implement" expressions allowing straightforward integration into computation algorithms. After some definitions on the mechanical behavior of cloth in Section 2.1, the computation scheme is detailed in Sections 2.2 and 2.3. We demonstrate the accuracy and efficiency of this scheme in Section 4.

The presented scheme only models tensile elasticity, which deals with in-plane deformations. Meanwhile, bending elasticity deals with out-of-plane deformations (surface curvature), and its main visible effect is to limit fold curvature and wrinkle size. However, the resulting internal forces are comparatively very weak, and very often, wrinkle size is rather limited by the size of the surface elements, which often bring more bending rigidity than the actual physical bending stiffness. Therefore, it can often be ignored in the context of interactive cloth simulation using large elements and soft cloth materials. However, for high-accuracy simulations, the presented tensile model can easily be complemented by a bending model using the schemes defined by [START_REF] Grinspun | Discrete Shells[END_REF] or [START_REF] Volino | Simple Linear Bending Stiffness in Particle Systems[END_REF]].

THE SIMULATION SCHEME

Tensile Viscoelasticity

From the theory of elasticity [START_REF] Timoshenko | Theory of elasticity[END_REF]Goodier, 1970] [Gould, 1993], the internal tensile deformations of a surface are characterized by its strain, measured through a strain tensor, represented by three independent values ε uu , ε vv , ε uv to the coordinate system of the material. In dynamic systems, their rate is measured though their time derivatives ε uu ', ε vv ', ε uv '. Meanwhile, the internal tensile forces are characterized by its stress, measured accordingly though a stress tensor represented by three independent values σ uu , σ vv , σ uv . The strain and stress values are related through the current energy per surface unit w of the material by the following relationships, for any deformation mode m among (uu,vv,uv):

m m w ε σ ∂ ∂ = (1)
The relationship between strain and stress defines the mechanical behavior of the material. In the most general context, this is expressed through the following strain-stress relationship:

( ) ( ) ( ) uv vv uu uv vv uu uv uv vv uu uv vv uu vv uv vv uu uv vv uu uu ε ε ε ε ε ε σ ε ε ε ε ε ε σ ε ε ε ε ε ε σ ′ ′ ′ ′ ′ ′ ′ ′ ′ , , , , , , , , , , , , , , , (2)
An isotropic material behaves exactly identically whatever its orientation, and its strain-stress relationship does not depend on the orientation of the considered coordinate system.

In the case of linear viscoelasticity, the strain-stress relationship can be expressed as a linear expression, the elastic and viscous stiffness of the material being represented as symmetric matrices E and E':

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ′ ′ ′ ′ + ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ uv vv uu uv vv uu uv vv uu E E ε ε ε ε ε ε σ σ σ (3)
In the particular case of isotropic linear elasticity, the behavior of the material is only described with two parameters: The Young modulus e relates the stiffness of the material while the Poisson coefficient ν relates its transverse contraction upon extension (Fig. 3). The corresponding matrix is the following: Many simulation models exist for dealing with linear materials [Etzmuss et al, 2003]. However, through the complexity of its structures, cloth materials have very nonlinear behaviors: The internal deformations of the materials produce internal forces which are not at all proportional to them. Linear strain-stress models are a very rough approximation of the actual behavior of cloth materials, and they quickly show their limits for large deformations (as for instance, the "superelasticity" artifact observed by [START_REF] Provot | Deformation Constraints in a Mass-Spring Model to Describe Rigid Cloth Behavior[END_REF]) (Fig. 5). Therefore, accurate cloth simulation requires nonlinear models for fitting closely the simulated cloth behavior to the actual one, usually measured through specific tensile tests (see Appendix). Our main interest is however to address the case of nonlinear materials through their most general strain-stress behavior (2), which is needed for representing accurately the large-scale behavior of structures which have high small-scale complexity, such as cloth.

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - = 2 1 0 0 0 1 0 1 1 2 ν ν ν ν e E (4)
Figure 5: A large square piece of cloth, hanging attached at two corners. The material is isotropic, with linear (left) (red curve far left) and nonlinear (right) (green curve far left) elongation strain-stress behavior. Colorscale shows average tensile stress (back) and strain (front). In both cases, there is a high concentration of stress around the attachment points which support the weight of the whole cloth. With a linear strain-stress behavior, this creates locally high strain, and therefore unrealistically large deformations in these areas. Meanwhile, a nonlinear strain-stress behavior smoothes out the strain and the deformations are much more evenly distributed over the cloth.

While elastic forces oppose deformation, viscous forces oppose deformation rate. Viscosity does not affect the rest geometry of the cloth in static contexts (draping), but it contributes to the damping of the cloth motion through energy dissipation. Meanwhile, plasticity results from hysteresis in the stain-stress behavior of the material. It is quite significant for cloth materials, and it mainly results from the friction between the textile fibers. Like viscosity, the main noticeable effect of plasticity is also the cloth motion damping caused by energy dissipation.

Our model can be extended for modeling plasticity through an adequate processing of the strain-stress behavior.

The easiest solution, derived from [START_REF] Breen | Predicting the Drape of Woven Cloth using Interacting Particles[END_REF] is to modify the rest strain of the strain-stress laws according to the current strain-stress state of the system, simulating the solid friction between the fibers. More advanced models, which take into account time-dependent evolution of the rest strain, could take advantage on models based on Prony series [START_REF] Soussou | Application of Prony Series to Linear Viscoelasticity[END_REF]. However, a precise modeling of plasticity requires an accurate knowledge of the possibly time-dependent strain-stress behavior inside the hysteresis envelope (not only the envelope itself) and this requires measuring the tensile properties of the cloth materials with more sophisticated procedures than the single load-unload cycle usually performed (Fig. 6). The complete study of plasticity and its simulation is a topic on its own, far beyond the scope of this work. Furthermore, the numerical simulation is itself quite unpractical and very expensive due to the numerical problems caused by the high nonlinearity of hysteresis. This is why not many cloth simulation systems do actually support accurate simulation of plasticity. While unnecessary in the context of cloth draping applications, the accurate dynamic reproduction of cloth motion along time requires significant damping for avoiding unrealistically "bouncing" materials. Whereas viscosity cannot be explicitly measured in cloth materials through simple and standard test procedures, it is nevertheless a simple solution for reproducing the global dissipative behavior caused by plasticity along viscosity. In this context, some techniques attempt to identify viscous parameters of cloth materials by experimentally matching the simulated damping of its motion to actual samples [START_REF] Bhat | Estimating Cloth Simulation Parameters from Video[END_REF]].

In the following, we formulate our model for the context of general nonlinear anisotropic tensile viscoelasticity

(2). Adequate simplifications may then be carried out according to the actual simulation context.

Computational Details

Our algorithm processes triangle elements of the mesh describing the surface (Fig. 7 left). The mesh does not need to have grid regularity, and may be obtained on arbitrary surface regions by processes such as Delaunay triangulation. Each element is described by its 2D parametric coordinates (u a ,v a ),(u b ,v b ),(u c ,v c ) of its vertices, referring to an orthonormal parametric coordinate system (in the context of cloth simulation, aligned to the weft and warp fiber directions). The current position of the deformed element is defined by the 3D world coordinates P a ,P b ,P c of its vertices, and possibly velocity coordinates P' a ,P' b ,P' c . The weft and warp vectors are expressed in 3D world coordinates as U and V, which are not necessarily orthonormal anymore because of material deformation (fig. 7 right). In the following, these vectors will be used for measuring the deformation state of the element, as well as expressing any vector value related to the element in world coordinates. Solving these linear systems leads to the following weights, to be precomputed:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) b a c a c b c b a a b vc c a vb b c va b a uc a c ub c b ua v v u v v u v v u d u u d r u u d r u u d r v v d r v v d r v v d r - + - + - = - = - = - = - = - = - = - - - - - - 1 1 1 1 1 1 (6)
During the simulation, these values are the weights for computing the current 3D vectors U and V directly as a weighted sum of the current vertex positions P i as follows:

∑ ∑ ∈ ∈ = = ) , , ( ) , , ( c b a i i vi c b a i i ui P r V P r U (7)
When viscosity has to be considered in the context of dynamic simulations, the current evolution rates of the coordinate vectors U' and V' can be computed as well from the current vertex velocities P' i :

∑ ∑ ∈ ∈ ′ = ′ ′ = ′ ) , , ( ) , , ( c b a i i vi c b a i i ui P r V P r U (7')
Our model is based on the Green-Lagrange strain tensor, which allows the rotation-invariant description of internal surface strain in the context of large displacements. This symmetric tensor G is defined by the coordinate vectors as follows, I being the identity matrix representing the rest state:

[

] [ ] ( )

I V U V U G T - = 2 1 (8)
From this tensor, we can extract the weft warp and shear strain values, which respectively measure the elongation deformations along weft warp directions and the shear deformation between them, as follows:

( ) ( ) ( ) U V V U V V U U T T uv T vv T uu + = - = - = 2 1 2 1 2 1 1 1 ε ε ε (9)
At this point, we can note that the strain is a quadratic expression of the vertex positions. We do not attempt to perform any linearization of these expressions.

Similarly, if viscosity is considered, the strain rate values are computed accordingly:

( ) ( ) ( ) U V V U V V U U T T uv T vv T uu ′ + ′ = ′ ′ = ′ ′ = ′ ε ε ε (9')
Having computed the strain state of the triangle surface, we can now obtain the stress state by using the strainstress relationship (2) that characterizes the material of the surface.

The Green-Lagrange strain tensor is associated to the second Piola-Kirchhoff stress tensor through (1).

Therefore, the forces F j exerted on the vertex j is computed by derivation of the weft, warp and shear components of the total elastic energy W of the triangle relatively to the particle position P j . Since we assume linear deformation of the triangle element, we have constant strain and stress and therefore surfacic energy w over its surface of area |d|/2 (6), and we have, for any j among (a,b,c): With explicit expression of the derivatives of the Green-Lagrange strain values (9) using (7), we obtain:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ - = ∂ ∂ - = ∂ ∂ - = ∑ ∈ ) , , ( 2
( ) ( ) ( ) ( ) U r V r V r U r d F vj uj uv vj vv uj uu j + + + - = σ σ σ 2 (11)
Direct implementation of ( 7), ( 9), ( 2), ( 11) can be the basis of an accurate Particle System model integrated with explicit numerical methods, such as Runge-Kutta [START_REF] Eberhardt | A Fast, Flexible, Particle-System Model for Cloth Draping, Computer Graphics in Textiles and Apparel[END_REF].

Computing the Jacobian

The Jacobian of the forces is necessary for implementing implicit numerical integration methods [Baraff and Witkin, 1998] [START_REF] Eberhardt | Implicit-Explicit Schemes for Fast Animation with Particles Systems[END_REF][START_REF] Hauth | A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods[END_REF]] [Volino et al, 2005], and any other relaxation schemes based on Newton iterations. Its accurate evaluation is necessary for good convergence and stability of state-of-theart simulation systems. This Jacobian is computed from (10) and (2) as follows, for any i and j among (a,b,c): Considering viscosity, the following contribution has to be considered as well: With explicit expression of the derivatives of the Green-Lagrange strain values (9) using ( 7), we obtain (for simplicity, here written without the contributions of the cross-dependencies between the deformation modes):

⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂ - = ∂ ∂ ∑ ∑ ∈ ∈ ∈ ) , , ( ) , , ( ) , ,
⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ′ ∂ ∂ - = ′ ∂ ∂ ∑ ∈ ∈ ) , , (
( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ + + + + + ∂ ∂ + ∂ ∂ + ∂ ∂ - = ∂ ∂ I r r r r r r r r U V r r V U r r V V r r U U r
( ) ( ) ( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ′ ∂ ∂ + ′ ∂ ∂ + ′ ∂ ∂ - = ′ ∂ ∂ T ui vj T vi uj uv uv T vi vj vv vv T ui uj uu uu i j U V r r V U r r V V r r U U r r d P F ε σ ε σ ε σ 2 (13')
This Jacobian, which is globally symmetric thanks to the symmetric strain-stress laws (2), has a fairly simple structure, and the numerical solving algorithm should take advantage of this for numerical solving or sparse storage.

The expression of the elasticity Jacobian (12) (13) has two components: * A stiffness component, which depends on the strain-stress stiffness ∂σ/∂ε. * A geometric component, which depends on the stress value σ.

The expression of the viscosity Jacobian (12') (13') has only a stiffness component.

It can be noted that in the context of linear viscoelasticity (3), the terms ∂σ/∂ε and ∂σ/∂ε' and are simply the elements of the elastic stiffness matrices E. and E'.

The stiffness component represents how the strain-stress stiffness of the material acts in the relationship between particle forces and positions. It is usually the most important component of the Jacobian. Meanwhile, the geometric component represents the particle force changes which result only from the evolution of the element geometry. For instance, it relates the rotation of the particle forces as the element rotates.

In the context of small deformations, it would make sense to ignore the geometry component of the Jacobian, assuming that the vectors U and V do not evolve far from their initial state. This is typically the approach done in approaches which use a linear approximation of the Green-Lagrange tensor, possibly addressing larger deformations through the corotational approach [Etzmuss et al, 2003]. However, as pointed out by [START_REF] Choi | Stable but Responsive Cloth, Computer Graphics[END_REF] in the context of spring-mass systems, the full evaluation of the Jacobian is important for the stability of the simulation in the context of large deformations such as cloth simulation, where elements may exhibit large orientation changes within few iterations. Fortunately, through the use of the simple, non-linearized expression of the Green-Lagrange tensor, the expression of the geometric component is quite simple, and furthermore isotropic (it does not directly depend on the current values of U and V). Therefore, this approach is the best for performing stable simulations of highly deformable nonlinear materials through the use of a Jacobian accurately matching the actual deformed state of the mechanical system, whatever the amount of deformation (Fig. 8). A potential problem related to this expression results from the possible negative eigenvalues that the isotropic Jacobian component might introduce in numerical linear system matrix of the integration scheme. While this is not an problem in the context of the accurate linear system solving schemes used in most Finite Element models, this could potentially bring trouble if the Conjugate Gradient method is used, such as with most implicit Particle System integration methods. There is no much worry to have about this concerning the stiffness component, where positive eigenvalues are ensured by the physical plausibility of the material (where usually forces monotonously oppose deformations). However, the geometric component is more likely to generate negative eigenvalues, particularly when material compression state creates negative eigenvalues in the stress tensor. While, in the context of cloth simulation, this is not much an issue as surface buckling quickly absorbs large tensile compression, this still needs to be addressed for obtaining an instability-proof simulation system (transitionally highly compressed elements can always occur during simulation). While filtering out negative eigenvalues from the global Jacobian is not a computationally practical option, a formal, fairly conservative solution would be to filter out the negative eigenvalues of the stress tensor through adequate projection of the σ uu , σ vv , σ uv values to be used in the computation of the geometric component (similarly to what [START_REF] Choi | Stable but Responsive Cloth, Computer Graphics[END_REF] does with spring-mass systems). In practice, we have experimented that a sufficient approximation for ensuring stability is to use null values of σ uu or σ vv if they are negative, thus avoiding any eigensystem computation.

Algorithm Summary

The simple computational process can directly be formulated in a ready-to-implement manner through the Particle System formalism, summarized as follows: Precomputation * On each element, evaluate and store the vertex distribution factors using (6): Computation * On each element, evaluate the current material referential U and V from P a , P b , P c using (7).

* From this, compute the strain values ε uu , ε vv , ε uv using (9). * And compute the stress values σ uu , σ vv , σ uv using the mechanical behavior of the material (2). * Then compute the corresponding particle force contributions F a , F b , F c using (11). * If needed, compute the Jacobian contribution using (13) with the partial derivatives of (2).

CONSIDERATIONS

St.Venant-Kirchhoff Materials

Considering a linear strain-stress relationship (3), the material is a hyperelastic St.Venant-Kirchhoff material, which has a cubic tensile force-deformation relationship (Fig. 2). This nonlinear behavior is indeed much more representative of the nonlinear tensile behavior of usual cloth materials than a linear tensile force-deformation behavior.

The largest interest of such a model is the simplicity of its computation, and the computation of ( 7), ( 9), ( 11), (13) can also be carried out with additions and multiplications only, without the need of square roots (no vector normalizations) or trigonometry (no coordinate rotations). Hence, it is possible to implement a very simple simulation system suited for real-time or interactive applications. Furthermore, there are also no discontinuities or singular configurations, and this makes such model very appropriate for robust interactive systems.

Mass Lumping

In our Particle System scheme, the mass given to each particle should be exactly one third of the mass of all triangle elements adjacent to that particle in the mesh. Obviously, the mass of an element is the surfacic mass of the cloth material times the area |d|/2 of the triangle.

A significant difference between Finite Elements and Particle Systems is that while the former considers accurately the mass of the mechanical system distributed over the surface of the elements, the latter condensates the mass of the elements on the discrete particle locations. Hence, the mass of a uniform triangle element is distributed equally on its three vertices. As a side effect, the rotational inertia of the element around its mass center becomes four times larger. In most contexts however, this is not a significant issue since the rotational inertia of small elements is usually not significant compared to their translational inertia. As a major benefit of this approximation, the inertia matrix only contains diagonal elements, and this drastically simplifies the mechanical computations. This is why this approximation, known as mass lumping, is also often considered in usual Finite Element models.

PERFORMANCE

This model has been integrated in a cloth simulation system that allows interactive editing of cloth objects and collision processing, for the creation and simulation of complete garments.

In our implementation, the mechanical model accepts independent weft, warp and shear elastic strain-stress curves modeled as polynomial splines of any order. Possible cross-dependencies among these modes are also accepted by our system, for instance modeling transverse contraction. Viscosity is also modeled in the same way.

External forces include gravity, anisotropic viscous aerodynamic drag (wind) and collision effects (friction).

The Particle System is integrated with Backward Euler, Implicit Midpoint [Volino and Thalmann, 2005] or alternatively BDF-2 [START_REF] Hauth | A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods[END_REF] numerical integration methods for performing dynamic cloth simulations.

Meanwhile, for cloth relaxation and draping applications, an iterative Newton resolution scheme which finds the particle positions that minimize the particle forces through the Jacobian was also implemented. In any of these implicit schemes, our implementation uses, for each iteration, the Jacobian of the forces corresponding to the actual current state of the system. This variable Jacobian scheme, which captures accurately all the nonlinearities of the system, is made possible through on-the-fly evaluation of the Jacobian directly inside the numerical solving of the computation iteration, which is performed using the Conjugate Gradient method.

The computation code is written in standard C++ using double precision floating-point, and the tests are performed on a 3GHz Pentium4 PC.

Accuracy of the Model

Our accuracy test consists in simulating a virtual material described by its weft, warp and shear nonlinear strainstress curves, modeled from force-displacement curves of a tensile test. Through simulation, we perform a virtual tensile test of the material, comparing the resulting force-displacement curves to the experimental ones.

We have performed this test using a 150g/m 2 wool gabardine fabric, which is a fairly nonlinear material, highly anisotropic with weak shear stiffness.

We have first measured the tensile weft, warp and shear (along weft) force-elongation curves using a tensile tester on the normalized 20cm×5cm sample. The averages of the load and unload curves, once converted into strainstress curves according to the formula given in the Appendix, are then modeled with polynomial splines (3 segments, 3rd-order for elongation and 2nd-order for shear).

The virtual simulation of the experimental setup has then been carried out on a virtual sample of 2500 triangles:

The cloth rectangle is attached along its two longest edges, and the total attachment force along one edge is measured according to its displacement. No other external force is considered and gravity is set to 0. Between each state change, the cloth equilibrium was computed using the Newton relaxation method.

The following force-displacement curves give the comparison between the polynomial model modeling the experimental data and the virtual test data (strain-stress curves converted back to force-displacement curves) (Fig. 9). From these curves, we can see that the force-displacement behavior our simulation system duplicates very precisely the force-displacement model resulting strain-stress polynomial spline model (converted back to forcedisplacement). The error remains always below 0.01%, with force differences below 0.01 N whereas total forces may exceed 100 N.

Tensile Weft Force-Displacement Curve

This illustrates the very high accuracy that our model can provide even in the context of nonlinear anisotropic materials in the context of large deformations. Such accuracy is indeed expected, since our model represents accurately the mechanical constitutive laws of the material model through continuum mechanics without any approximation. Furthermore, since the computation of the Jacobian is also exact, the model converges fairly quickly to the theoretical value described by the curve model (in less than 10 iterations to the error specified above, and only limited by the numerical floating-point accuracy of the computations).

Real-World Accuracy

In order to verify the validity of our model for simulating cloth in a more "real-world" context, we have set-up an experiment for comparing the behavior of the real cloth with the simulated virtual cloth.

A piece cloth is clamped horizontally on a circular frame (30cm diameter), with as little initial tension as possible. Then, a vertical force is exerted on a disk (6cm diameter) on the middle of the cloth. The vertical displacement of the disk is measured according to that force (Fig. 10). We have chosen a fairly soft terry cloth, so as to observe fairly large displacements that can be accurately measured.

For the real tests, we have performed an average of several measurement sequences, which allow us to assess the accuracy of these measurements through the computation of the standard deviation. Observed measurement variations (around 10%) would mainly result from variations of the initial tension of the cloth on the frame, as well as some possible material plasticity effects of the material between each measurement. For the simulation, strainstress curves have been measured accurately using tensile tests (as described in the Appendix), modeled as three independent (weft, warp, shear) polynomial splines ( 14) averaging the hysteresis loop, and implemented in our simulation system. The measured curve (Fig. 10 bottom right) shows that our simulation can reproduce the real cloth with a fairly good accuracy. Along inaccuracy causes mentioned above, another could be the approximate modeling of the material, not taking into account possible dependencies between weft warp and shear deformation modes, as discussed in the Appendix.

Computation Time

Computational timings were performed on the dynamic simulation of a 1m×1m cloth square (Fig. 11), initially horizontal, attached along one of its edges. Numerical integration was performed using the Inverse Euler integration, with constant simulation timesteps of 10 milliseconds.

For comparison, we have implemented in our simulation system the linearized corotational scheme, obtained by rotation of the local element coordinates to the eigendirections of the strain tensor and linearization of the strain and stress expressions. It can be noted that the expression of the Jacobian is more complex (because of the linearization)

and requires approximations (because of the state-dependent rotation).

We have compared the computation time of our scheme and the linearized corotational scheme for computing 100 timesteps of the computation with several mesh resolutions. The considered material is a simple isotropic Non-Linearized

Spring-Mass

Figure 11: Computation times required for the iterating the computation of 1m 2 of cloth, using various mesh resolutions, using our simulation scheme, the linearized corotational scheme and a simple spring-mass scheme.

From these curves, we can see that our implementation of the model is able to iterate more than 17 500 elements per second. Meanwhile, the linearized corotational scheme is about two times slower. This slowdown results not only from the computation of the eigensystem for evaluating the rotation, but also from the additional coordinate rotations done in the computation. Furthermore, using a linearized scheme also requires additional vector normalisations (square roots), as well as a more complex, non-isotropic geometric component in the exact expression of the Jacobian (13). It should be noted that the numerical solving process (involving usually between 10 and 15 iterations of the Conjugate Gradient algorithm for each timestep, depending on the stiffness) also accounts for more than half of the total computation time.

The interesting point is that the simple spring-mass system is only marginally faster compared to our scheme, with only 20 000 elements iterated per second. Hence, our simulation scheme offers a huge increase of simulation accuracy for only about 15% additional computation.

Tested in the same conditions with our simulation scheme, the nonlinear anisotropic strain-stress behavior defined is the previous section requires roughly 10% additional computation.

Application: Virtual Prototyping

The presented model has been implemented in a virtual prototyping system which allows the design of complex garments on mannequins. The efficiency of the model allows the simulation of the accurate mechanical properties of cloth on high-resolution meshes as the mannequin moves, allowing the garment designer to assess the stretch forces for particular postures (Fig. 12). The design of the garment patterns may then be corrected accordingly, with interactive mechanical feedback on the garment drape. The Particle System approach allows complex multilayer garments to be processed through an adequate handling of complex collision situations. Through the use of various implicit integration schemes, the model can efficiently compute both static drapes and dynamic animations.

CONCLUSION

Through the presented computation scheme, we have obtained an accurate simulation system which offers highaccuracy representation of cloth deformation based on continuum mechanics. Through the use of the non-linearized Green-Lagrange tensor, this system performs the simulation of hyperelastic St.Venant-Kirchhoff materials, which are particularly adapted for representing the behavior of cloth materials. Yet, this scheme is also extended for supporting accurately nonlinear and anisotropic representations of weft, warp and shear strain-stress curves, which can be obtained from force-deformation curves measured with traditional tensile tests on actual fabric samples.

The presented scheme is one of the simplest available featuring the accuracy of continuum mechanics. It remains accurate for large deformations without the need of complicating the computation process with additional transformations, such as corotational formulations. Furthermore, it can be expressed as a classical Particle System (with particle forces directly computed from particle positions), offering a high versatility of applications in numerous simulation contexts (draping, dynamic simulation, collision response and other geometric constraints, etc). The exact Jacobian of the model can also be efficiently computed, offering efficient and stable simulations through the use of implicit integration methods (Inverse Euler, BDF-2, etc).

Through these qualities, this simulation scheme is particularly well adapted for applications that need to combine good mechanical accuracy with computation times compatible with interactive applications. It is therefore a very good candidate for garment prototyping applications.

The highly streamlined computational process also opens the doors for parallelization and hardware implementation, which, through the rise of dedicated chips, represent the future of high-performance mechanical simulation.

Figure 1 :

 1 Figure 1: Accurate simulation of nonlinear anisotropic cloth materials is required for garment prototyping applications.

Figure 2 :

 2 Figure 2: Tensile elongation force-deformation curve of a St.Venant-Kirchhoff material (red) compared to a linear material (blue), global view (left) and close-up (right). On the St.Venant-Kirchhoff material, we can note the "collapse" of the compression force for large compression deformation.

Figure 3 :

 3 Figure 3: A square piece of cloth attached at two opposite edges (far left), extended at 150% of its initial length. The material is linear isotropic (4), with a Poisson coefficient ν of 0 (left), 0.25 (center), 0.50 (right). Colorscale shows transverse compression strain.

Figure 4 :

 4 Figure 4: A large square piece of cloth cut along various fiber directions, hanging attached at one edge. The material is linear anisotropic (3), with 2 E uu,uu = 1 E vv,vv = 8 E uv,uv (radial elongation stiffness diagram far left). This represents typical anisotropy of cloth materials, which usually have a much lower shear/elongation stiffness ratio than isotropic materials. Strips show weft fiber direction. Colorscale shows average tensile strain.

Figure 6 :

 6 Figure 6: Weft (left) and warp (right) elongation tensile test of a cloth material with several load-unload cycles of various amplitude. The precise evaluation of material plasticity requires a complex study of strain-stress hysteresis, far beyond the evaluation of the hysteresis envelope only.

Figure 7 :

 7 Figure 7: A triangle element is defined in 2D parametric coordinates (left), and deformed in 3D world coordinates (right).

Figure 8 :

 8 Figure 8: An accurate Jacobian is necessary for stable simulation of very large deformations without compromising convergence speed.

Figure 9 :

 9 Figure 9: Comparing simulated weft, warp and shear force-displacement curves with their counterparts defined in the model. With less than 0.01% difference, the curves are practically identical. The virtual tensile setup is shown at bottom right.

ForceFigure 10 :

 10 Figure 10: The experimental setup (top), and the virtual simulation (bottom left). Depending on the force exerted by the weight on the middle of the cloth, the vertical displacement is plotted (bottom right), comparing the simulation (green curve) with the average of several experimental measurements (red dots with error bars showing average deviation).

  material (3) (4) of Young modulus e = 1000 N/m, null Poisson coefficient ν = 0. We have also compared the computation time of a simple spring-mass system (springs defined along mesh edges) having an equivalent stiffness.

Figure 12 :

 12 Figure 12: Virtual prototyping applications require an accurate representation of cloth material behavior for evaluating precisely the stretch forces (color scale) on the garment along particular postures of the character (top). An efficient simulation model is also required for the dynamic simulation of high-resolution complex garments involving several layers of cloth.

ACKNOWLEDGEMENTS

We are grateful to all MIRALab team who participated to this research through their contributions, advices, testing and creative artwork. Special thanks for Christiane Luible for her expertise in garment design and animation. This research is funded by the European FP6 project LEAPFROG IP (NMP2- CT-2005-515810).

BIBLIOGRAPHY

APPENDIX: TENSILE TESTS FOR MEASURING STRAIN-STRESS CURVES OF CLOTH

The tensile behavior of cloth materials is usually characterized by its tensile elasticity strain-stress relationship, measured through adequate tensile tests. Standard procedures exist for characterizing cloth properties. For instance, the Kawabata Evaluation System [START_REF] Kawabata | The Standardization and Analysis of Hand Evaluation[END_REF] defines normalized hardware and tests for measuring weft warp elongation and shear strain-stress curves. It is important to note the material is only considered through independent curves relating weft and warp elongation and shear, with no consideration of any possible cross-dependencies between deformation modes, and with also no consideration of the deformation rate (which is set to a single standard value). Then, the parameters of the cloth are evaluated as a set of characteristic values quantifying the shape of these curves. Since our concern is only related to the strain-stress curves of tensile elasticity, we can exploit weft warp elongation and shear curves measured by the Kawabata hardware (Fig. 13 left), or use other tensile tests that may capture more complex behaviors in wider range of deformations.

According to the Green-Lagrange tensor (8), The simplest and most natural way to integrate measured curves into our computation framework is to consider a simplified relationship between strain and stress (2) expressed as three decoupled (weft, warp and shear) curves, as follows:

However, through the nonlinearity of the Green-Lagrange tensor (8), the weft, warp and shear strain and stress values cannot be directly identified to the measured force-deformation curves through a simple linear transformation. The purpose of this section is to establish the adequate conversion formulas (Fig. 14), which should remain valid for large deformations.

For this, we consider a general tensile test (Fig. 13 Using ( 9) with the deformed values of U and V, we obtain the conversion formula from displacements to strain values:

We also obtain the conversion formula from stress values to forces:

The forces are converted back to stress values as follows: It can be noted that such tensile tests are unable to measure compression stiffness of cloth materials, since these buckle and exhibit wrinkle patterns upon compression because of their low bending stiffness. Nevertheless, cloth materials do have some compression stiffness, which is indeed necessary for their buckling behavior. Therefore, realistic cloth simulation requires compression stiffness to be modeled. High-accuracy models are however not required, since buckling allows relaxation and limits the actual cloth compression to low values. In practice, as compression stiffness is generally not measured, we typically extrapolate elongation strain-stress curves to negative values using antisymmetric functions.

It can also be noted that elongation and shear measurements can be carried out both along weft and warp directions (as done in the Kawabata standard), leading to a total of four curves. However, the mathematical definition of tensile deformation only allow its description as three independent values, that we have chosen in ( 14) to be weft, warp and shear, as defined by the Green-Lagrange tensor (8) (9). Assuming independence of the two elongation modes, the two shear force-deformation curves cannot be unrelated, and they should indeed produce identical shear strain-stress curves after conversion (this should also be true in any case in the context of linear elasticity). Therefore, a mismatch between them indicates that there is a significant nonlinear dependency between the deformation modes.

More generally, some behaviors of the cloth material cannot be evaluated using curves from the Kawabata tests only. These include those resulting from the coupling of deformation modes, as for example transverse shrinking (which, in the context of linear elasticity (4), is represented by the Poisson coefficient) (see Fig. 3 in Section 2.1).

Also, non-symmetric cloth behaviors (possibly caused by non-symmetric yarn patterns) need to be modeled with coupling between elongation and shear modes. It would be possible to evaluate with better accuracy the nonlinear strain-stress behavior of a cloth through more comprehensive force-deformation tests than what is proposed in the Kawabata standard, by combining simultaneously various values of elongation and shear deformations. Adequate interpolation functions would then model this data as a general strain-stress relationship (2) to be used as input in the proposed simulation scheme.