Eitan Grinspun

Pierre Alart

Bruno Lévy

Fabrice Neyret

Sébastien Neukirch

Roberto Santoprete

Stéphane Redon

PhD Alexandre Derouet-Jourdan

PhD Romain Casati

engineer, then PhD Gilles Daviet

Alejandro Blumentals

Florent Cadoux

Joëlle Thollot

Franck Hétroy

Bernard Brogliato

Pierre Saramito

Arnaud Lazarus

Damien Rohmer

Jean-Sébastien Franco

Rahul Narain

Marie-Paule Cani

Keywords: 17 2, 3, 1 Piecewise uniform curvature: the super-helix model, 19 2, 3, 2 Piecewise linear curvature: the super-clothoid model, 20, 26 2, 4, 1 Recursive integration scheme, 27

My activity has also involved industrial transfer, and I would like to thank all the industrial partners, from Néomis Animation & Beelight to L'Oréal Research, I had the opportunity to work with, as well as the Inria valorisation cell for making such collaborations possible. I am grateful to Bernard Brogliato for having welcomed me in the BiPop group in 2007. Being the sole researcher from Computer Graphics among a set of scientists focussed on more fundamental topics such as Mechanics or Optimization was a challenging situation for me, but corresponded to the scientic environment I was looking for at that time. Beyond the scientic immersion into nonsmooth mechanics and optimization, working in BiPop brought me the opportunity to reach some scientic independence very soon. It also denitely helped me step back and get a fresh point of view on various domains, including my own community. I precisely take this opportunity to thank my home community, Computer Graphics, for the intense creativity I found in all the people working in this eld, and the overall benevolence ruling over people. Moreover, I still remain impressed by the professionalism of Siggraph and Siggraph Asia program committees, which, I believe, is truly exceptional, and had a strong impact on my evaluation practice in various other contexts. Thinking back to my rst years of research, I feel very grateful to my PhD advisors,

Chapter 1 Introduction J'ai pas mérité de jouer du rock'n'roll. Mes ghettos, mes idées ne sont pas homologués. J'ai pas le bon blouson, j'ai pas les bonnes bottes. Et en haut de mon bras je n'ai rien fait tatouer. J'ai donné aux curés du sauvetage collectif. J'ai joué les mêmes notes, swingué les mêmes ris. Peu à peu j'ai compris les données du débat. Que rien ne bouge et l'égalité par le bas. Et tant pis si la foule gronde. Si je tourne pas dans la ronde. Papa quand je serai grand, je sais ce que je veux faire : je veux être minoritaire.

Jean-Jacques Goldman, Minoritaire

My research activity deals with the mechanical modeling and numerical simulation of elastic slender structures, in particular bers, modelled as thin elastic rods, but also surfaces in some more recent work, elaborating on models for thin elastic plates or developable shells. I've also got interested in the capture of dry frictional contact among such thin structures, as well as within granular ows. All these dynamical systems are useful to model a large range of physical phenomena at dierent scales, from twisted DNA to the mechanics of twining plants, sand, human hair, and cloth (see Figure 1.1). Target applications encompass visual eects for the feature lm industry, as well as virtual prototyping for various domains such as cosmetology (hair), virtual try on (hair and cloth), or risk management (granular materials).

Figure 1.1 Our main simulation results: from the numerical modeling of thin elastic rods and shells (top) to the coupling of rigid bodies, rods and shells for capturing the dynamics of complex macroscopic phenomena involving dry frictional contact (bottom).

2

Chapter 1. Introduction

Simulating visually rich phenomena

Deformable slender structures such as hair bers, cloth, ribbons, paper, cables, tree branches or leaves, are ubiquitous around us. It is actually striking to note that such thin structures represent most of the deformable objects composing our direct environment, but also a wide range of physical, biological or manufactured phenomena occurring at very dierent scales, from macromolecules and carbon nanotubes at nanoscale to galactic laments at the universe scale (see Figure 1.2).

At the macroscopic scale (and even smaller scales), slender structures often feature an intricate natural shape, ranging from straight to curly. For small strains, their deformation is mostly elastic (with a preference for bending and twisting when no overly constrained), yet generating large displacements due to geometric nonlinearities, and instabilities such as buckling. The diversity and complexity of the resulting shapes and motions greatly contributes to the visual richness of the real world.

When multiple thin structures are coupled together with contact and friction, which is the predominant interaction mode at the macroscopic scale, the range of emerging phenomena is even more exacerbated, giving rise to stick-slip dynamical instabilities, entangling, or spontaneous collective behavior. Human hair, which is typically composed of 150,000 thin bers, beautifully depicts such complex mechanical behaviors when uttering in the wind. 1.1. Simulating visually rich phenomena 3

Many physicists and mathematicians have strived for centuries to understand the principles governing those complex mechanical phenomena, providing a number of continuous models for slender structures, granular matter, and frictional contact.

In the XX th century, industrial applications such as process automatization and new ways of transportation have boosted the elds of Mechanical Engineering and

Computer-Aided Design, where material strength, reliability of mechanisms, and safety, standed for the main priorities. Instead, large displacements of structures, buckling, tearing, or entanglement, and even dynamics, were long considered as undesirable behaviors, thus restraining the search for corresponding numerical models.

In contrast, Computer Graphics, which has emerged in the 60's with the advent of modern computers, was eager to capture such peculiar phenomena, with the sole aim to produce spectacular images and create astonishing stories.

At the origin, Computer Graphics thus drastically departed from other scientic elds. Everyday-life phenomena such as cloth buckling, paper tearing, or hair uttering in the wind, mostly ignored by other scientists at that time, became actual topics of interest, involving a large set of new research directions to be explored, both in terms of modelling and simulation.

Since a few decades, a new generation of physicists became interested again in the understanding of such visually fascinating phenomena, and started investigating the tight links between geometry and elasticity 1 . Common objects such as folded or torn paper, twined plants, coiled honey threads, or human hair [START_REF] Ben Amar | Crumpled paper[END_REF][START_REF] Goriely | Mechanics of climbing and attachment in twining plants[END_REF][START_REF] Brun | A numerical investigation of the uid mechanical sewing machine[END_REF][START_REF] Goldstein | Shape of a ponytail and the statistical physics of hair ber bundles[END_REF] have thus regained some popularity among the community in Nonlinear Physics 2 .

Yet, while the phenomena of interest have become remarkably close to those of Computer Graphics, the goals and employed methodologies still dier substantially from one community to the other, as explained in next section.

In parallel, the engineering industry has recently shown some new and growing interest into the modeling of dynamic phenomena prone to large displacements, contact and friction. For instance, the cosmetology industry is more and more interested in understanding the nonlinear deformation of hair and skin, with the help of simulation. Likewise, auto and aircraft manufacturers are facing new challenges involving buckling or entanglement of thin structures such as carbon or optical bers; they clearly lack predictive, robust and ecient numerical tools for simulating and optimizing their new manufacturing processes.

4

Chapter 1. Introduction

Computer Graphics as my home community

Computer Graphics is a very particular domain. From outside, one may hardly understand what this community is about if one does not know that it represents, above all, the gathering of people from various career paths (artists, industrials, scientists) who share a common passion for pictural arts and somehow an aesthetic vision of the world. Most people I have met in Computer Graphics, even among scientists and industrials, do have an accute sense for artistic creativity, many of them being keen on drawing (like myself ), on the making of movies and video games, or on any other creative activity. By construction, the scientic community in Computer Graphics is thus tightly linked to the digital entertainment industry, which encompasses all activities related to the making of animation feature lms, special eects, and video games. It is noteworthy that a large number of researchers who are active in the eld are employed by private companies such as Autodesk, Disney

Research, Pixar, Weta Digital, Microsoft Research, Adobe, etc. The most prestigious (and selective, see Figure 1.3) conference of the eld, ACM Siggraph, gathers every year tens thousands experts and practitioners including research scientists, engineers, and artists, fostering an extraordinary mixing between academia, private research centers and studios, and engineering and arts schools. Creating and directing visual digital worlds The primary and historic goal of Computer Graphics is to digitize the world into synthetic images targeted for a human audience. Computer Graphics thus has the ambition to represent visually all the phenomena we may perceive through our eyes, at all relevant scales. The dream application should allow a user to navigate in real-time within an innite virtual space, being able for instance to look up at the clouds swirling in a stormy sky, then viewing down a forest canopy shaking under the wind, next zooming into the details of the undulating tree branches and distinguishing every uttering leaf, even this small iridescent beetle climbing onto that hairy leaf... As if this were not 1.2. Computer Graphics as my home community 5 enough, the user may also interact with this virtual world as intuitively as possible, for example by grasping a wreath of owers out of the prairie, creating a new river down the hills or changing the shape of the mountains, the shininess of the ocean, or the course of the wind, just at her will.

From other scientic elds, such an ambition of creating and directing an entire digital world similar to the real one, with all possible levels of details, may look quite adventurous and even untenable. Indeed, the usual scientic methodology is rst to study simplied and rened phenomena, at a specic space and time scale, before daring to tackle more complex scenarios possibly integrating multiple bodies or multiple scales. The apparent audacity of Computer Graphics may however be explained by the fact that its goals, and thus its evaluation criteria, are drastically dierent from those of usual scientic elds. Unlike physics or biology, the goal here is not to explain observed phenomena, but to reproduce them on an array of pixels, as faithfully as possible. Unlike computational mechanics, the goal is not to compute accurate internal stresses inside a structure to anticipate rupture and strengthen the material at the right locations, the goal is to compute visible quantities such as deformation or velocity elds at the outer surface of the body, at a good resolution and free of any visual artifact; the goal is not to predict nely the parts of the building that will be destroyed by a re starting at a certain point, the goal is to create ames that will blu the audience and let them believe 3 they are surrounded by re; nally, the goal is not to characterize the exact shape and magnitude of a perfectly round drop of water falling into a perfectly at basin; the goal is to provide the user with intuitive tools to create the drop animation that will best serve the underlying storytelling and convey the right emotions, or say otherwise, that will satisfy the eye of the art director4 .

Having said that, creating outstanding visual eects still demands powerful computational models which cannot be built by hand, and instead need to be automatized. Following physical principles is an excellent way to achieve this. This is particularly true in the sub-eld of Computer Animation, which deals with the representation of moving phenomena.

3 The Computer Graphics community often uses the term plausibility to justify that its models are intended for human eyes which should be the sole judges of their validity. However I do not like this term which is often misused and may serve as a pretence to escape from validation, should it be physical, perceptual, or anything else (validation is a vast but still immature research topic in Computer Animation). I will thus attempt to dene my own evaluation criteria in Section 1.3, and will avoid using the term plausibility in this document. I however acknowledge the need for more and better validation of my own work, which I aim at developing in future work.

6

Chapter 1. Introduction Physics-based animation More specically, in the case of a dynamical scene involving passive objects such as hair, cloth, uids, etc., physics-based simulation has proven over the years to be a method of choice for synthesizing resulting visual eects. Unlike phenomenological methods which develop descriptive models for reproducing a given emerging phenomenon, physics-based methods provide generative models which explicitly integrate the physical causes of the phenomenon. From a set of initial conditions as well as a few physical parameters (e.g., the mass, stiness, natural shape), a physics-based simulator may generate not just a single eect, but a wide range of emerging phenomena revealing the whole complexity of the underlying physics.

My research work precisely lies within physics-based animation. The goal in our sub-community is to design numerical models which can reproduce convincingly very complex scenes, such as detailed breaking waves, a full garment folding on a walking character, or a full head of hair suddenly shaking. Although huge progress has been achieved in visual eects since the 90's, leading today to entirely synthetic animated outdoor scenes in most adventure movies, a large number of scenarios, such as the ones mentioned just above, still resists automatic digitization and requires many hours of manual work in production studios to be properly animated. Challenging scenarios typically involve a large number of degrees of freedom, strongly nonlinear dynamics, and many constraints in particular, contacts and friction. All these features are very dicult to capture numerically. A supplementary obstacle may be the lack of relevant (continuous) physical models in some cases. For instance, while the theory of Kirchho for thin elastic rods ts in well with our needs for modeling a bending rope or a tree branch at the right (human) scale, to the best of our knowledge there does not exist an integrated, continuum model for representing a full bulk of very thin bers, like hair. Physical modeling thus represents an important, pluridisciplinary and long-term challenge for bridging the gap between, on the one hand, available representations of a certain class of phenomena, and on the other hand, the will to capture more complex ones, at the right scales.

In the following I explain the methodology and guidelines I have set up for tackling some of the challenges listed above, in the particular cases of elastic slender objects (bers, cloth) and divided materials (granular and brous media).

Methodology and guidelines

When designing a physics-based simulator for computer graphics, I have inevitably in mind the four following criteria:

Realism: Ingredients which are necessary to capture relevant visual eects should be identied, incorporated in the model and translated numerically with as few quality loss as possible, at the right resolution, and without visual artifacts.

Robustness: The simulator should converge properly for a relevant range of parameters. This is all the more challenging as we are interested in objects prone 1.3. Methodology and guidelines 7 to large displacements, and subject to a potentially huge number of frictional contact constraints, which causes standard simulators to break down.

Eciency: The simulator should be fast enough for allowing complex scenes to be simulated in reasonable timings on a standard PC (in our case, a few days of computation for a given scene is an upper-bound).

Control: The simulator should provide the user with some handles to control the shape and motion of the object in an intuitive way, without having to tweak manually a signicant number of meaningless parameters.

Unlike the common belief according to which favoring one of these criteria is necessarily done at the price of sacricing others, I am striving to develop compact numerical models satisfying all four criteria at the same time.

Building dedicated models To achieve this, the streamline idea is to identify specic emerging phenomena to simulate, and then build dedicated numerical models taking advantage of an upfront reduced mathematical formulation. For instance, instead of treating bers and cloth as generic 3D continuum elastic models with a potentially small thickness (modeled with an ε parameter), one may directly consider equations for thin elastic rods, plates and shells, where reduction is undertaken upstream. Such dedicated models just retain the right ingredients one wishes to capture, without having to care about numerical issues inherent to the degeneracy of the model when ε becomes smaller and smaller. Similarly, so far we have considered only inextensible models for thin elastic rods (namely, the Kirchho theory for thin elastic rods), because when bending and twisting deformations are the main focus of interest, stretching and shearing then become irrelevant deformation modes, as proven by a simple scaling argument [Audoly and Pomeau, 2010, Section 3.7].

Finally, for modeling frictional contact as eectively as possible, we have been considering macroscopic models which directly cope with the emerging rigidity and nonsmoothness of contact and friction, instead of using locally compliant models which get stier and stier when approaching the right scale.

Targeting the full modeling pipeline Of course, on the down side, such a strategy requires that many dierent elds theory for thin elastic rods and shells, nonsmooth contact mechanics, etc. be investigated in deep, before tackling merely numerical problems, which themselves involve other scientic elds such as numerical analysis, optimization, and algorithmic. However, to design eective simulators, I believe it is essential to master the full modeling pipeline, from mathematical models to numerical schemes, including algorithmic choices and implementation issues.

From a personal point of view, I am also eager to investigate and learn new insights from all these dierent scientic elds (despite sometimes a certain frustration due to some obvious lack of time). Besides, at the scale of a person or a small group, I can testify it is quite satisfactory to build prototypes which are entirely home-made, and thus completely under control.

Chapter 1. Introduction Developing cutting-edge, but large usage codes Finally, one may argue that the benet of dealing with generic models is to build some generic code, which can deal with (almost) any problem, while dedicated models are limited to very specic codes, and thus to narrow applications. In the context of my research work, I partially disagree with this argument, for two reasons. First, a generic code may probably deal with many dierent situations (e.g., a thin rod but also a thick rod, elastic bodies but also uids, etc.), yet instead of treating optimally one of the possible scenarios (w.r.t. the criteria dened above), it is likely to do a rather poor job for all of them. Second, we have noticed, interestingly, that even though we have modeled each one of our phenomena specically, core portions of our codes show up some common features. For instance, it is striking to note that it is exactly the same nonsmooth solver, originally published in the context of dynamic hair simulation [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF], which has served to resolve our hair inverse design problem [Derouet-Jourdan et al., 2013a], and which has then been leveraged to simulate granular ows modelled as a continuum [Daviet and Bertails-Descoubes, 2016b,a]. The emergence of common canonical problems has thus allowed us to concentrate our eorts on their ecient implementation.

Aiming at a pluridisciplinary research Although my research has been primarily driven and motivated by Computer Graphics applications, which are the most inspiring to me, the methods that we have employed to address related scientic problems largely go beyond the sole eld of Computer Graphics, sometimes addressing questions of interest for other communities.

For one thing, as commonly done in Computer Animation, we have frequently borrowed tools and concepts from more fundamental elds such as optimization, geometry and mechanics. In the study of slender elastic structures, the originality of our approaches stemmed from the use of Cosserat elastic theories barely known in Computer Graphics before the 2000's and likewise barely popular in Mechanical Engineering, and which I started to familiarize with during my PhD thesis. In the eld of contact mechanics, some new lightnings could be brought by using the point of view of nonsmooth optimization which I beneted from due to my immersion within the BiPop team since 2007.

Besides, we have from time to time tackled scientic issues which surpass the sole eld of Computer Graphics. We have gradually realized that publishing some specic technical contributions to graphics may leave them buried under the application (the nal application being necessarily put forward in a graphics paper), whereas they could, in some cases, be relevant to other elds and benet from a more careful evaluation from experts outside graphics. In recent years, we have thus strived to better identify the nature of our contributions and attempted to conduct a more diversied publication policy. One consequence is that we have started to publish some of our work in other relevant elds such as

Structure of this document

The rst part of this document provides a summary of the three main topics of research I have been investigating since 2007 (date of my arrival at Inria as a Chargée de recherche), and is organized according to three dedicated chapters:

1. Modeling of dynamic bers (Chapter 2) In this rst research axis, we introduce a family of high-order, reduced models for discretizing the Kirchho equations for thin elastic rods in both a faithful and robust way [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF][START_REF] Bertails | Linear time super-helices[END_REF][START_REF] Bertails-Descoubes | Super-clothoids[END_REF][START_REF] Casati | Super space clothoids[END_REF]. Such models are particularly well-suited for simulating inextensible bers subject to bending and twisting, and featuring an arbitrary curly intrinsic geometry. We have recently started to extend such discrete 5 schemes to the numerical modeling of in- extensible developable shells [Blumentals et al., 2016a].

2. Modeling of frictional contact (Chapter 3) A second research axis deals with the modeling of contact and Coulomb friction for slender structures and granular ows within the nonsmooth contact dynamics framework. Beyond some recent analytical study of the problem of existence and uniqueness of solutions to the continuous contact problem [START_REF] Blumentals | The contact problem in Lagrangian systems subject to bilateral and unilateral constraints, with or without sliding coulomb's friction: a tutorial[END_REF], we have designed some numerical solvers to resolve the discrete contact problem in a both robust and ecient way.

Such algorithms were leveraged for robustly handling thousands packed bers at reasonable frame rates [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF], and more recently adapted to the simulation of frictional contact in cloth [START_REF] Daviet | Fast Cloth Simulation with Implicit Contact and Exact Coulomb Friction[END_REF], as well as the exact solving of the Drucker-Prager law for granular materials [Daviet and Bertails-Descoubes, 2016b,a].

3. Inverse static design of bers (Chapter 4) A last research axis, which will be the main starting point of my future work, is dedicated to nding new insights into the inverse modeling of static bers [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF], 2011, 2013b,a].

This work consisted in taking an arbitrary curve geometry as input and inferring corresponding geometric and physical parameters of the simulator such that the input geometry corresponds to a stable conguration at equilibrium.

12

Chapter 2. Discrete curvature-based models for dynamic bers perfectly inextensible. Such intrinsic inextensibility thus removes the burden of adding subsequent (sti ) inextensibility contraints when solving the dynamics. The price to pay, however, is that the geometry of the rod is not readily available but has to be computed iteratively from the curvatures. Moreover, the inertia matrix of such systems turns out to be dense, leading to a N 2 complexity when solving them at each time step, where N is a number of elements.

A rst strategy to resolve the latter issue consisted in exploring higher-order spatial discretizations, namely piecewise linear curvatures [START_REF] Bertails-Descoubes | Super-clothoids[END_REF][START_REF] Casati | Super space clothoids[END_REF], giving rise to the so-called super-clothoid model, so as to reduce the number N of elements needed for a given spatial accuracy. However, we were faced with the loss of a closed-form kinematics compared to the piecewise uniform case. Our main contribution was to devise an accurate integration scheme, based on power series expansion, which proved to be orders of magnitude faster compared to classical integration methods [START_REF] Casati | Super space clothoids[END_REF].

A second strategy was to propose a linear, recursive time-integration scheme for super-helices [START_REF] Bertails | Linear time super-helices[END_REF], inspired by the Featherstone linear-time algorithm for articulated chains of rigid bodies [START_REF] Featherstone | The calculation of robot dynamics using articulatedbody inertias[END_REF]. This approach allowed us to simulate long Kirchho rods made of hundreds elements at a much cheaper cost, and was also leveraged to design tree-like structures of rods in a very simple way.

Organization of the chapter Section 2.1 introduces prior work regarding the modeling and simulation of dynamic bers, and motivates our choice for reduced, curvature-based discrete models. Section 2.2 briey recalls equations for a Kirchho rod, with a particular stress on the structure of the kinematic problem known as the Darboux problem. Section 2.3 introduces the reduction of Kirchho equations to nite spatial dimension, and provides two dierent spatial schemes relying on a piecewise constant and piecewise linear shape function for the curvature, respectively. Section 2.4 gives the main ideas of our recursive integration algorithm, which may not only apply to the super-helix model, but to curvature-based discrete models of any order. Finally, Section 2.5 concludes this part and gives some insights into a possible extension of our work to the dynamic of surfaces.

Prior art and choice for reduced models

The scientic study of bers has a long history in various elds, tracing back to the rst continuous mechanics theories a few centuries ago to their further analysis in physics and mathematics, and their recent numerical treatment in Mechanical

Engineering and Computer Graphics. Motivation originates from a number of applications ranging from the understanding of DNA supercoiling [START_REF] Benham | DNA mechanics[END_REF] and climbing plants [START_REF] Goriely | Mechanics of climbing and attachment in twining plants[END_REF] to the simulation of submarine cables [START_REF] Goyal | Non-linear dynamic intertwining of rods with self-contact[END_REF], surgery threads and needles [START_REF] Pai | Strands: Interactive simulation of thin solids using Cosserat models[END_REF][START_REF] Chentanez | Interactive simulation of surgical needle insertion and steering[END_REF], or hair [START_REF] Ward | A survey on hair modeling: Styling, simulation, and rendering[END_REF]. Theories for thin elastic rods Various theories were proposed in mechanics to model the equilibria and the dynamics of bers, depending on the type of deformation considered. Our goal here is to capture the geometric richness of typical ber deformations such as waving hair, coiling cables, curled ribbons or twining plants. These phenomena are largely nonlinear, dominated by bending and twisting elastic deformations, while stretching and shearing can be neglected. To account for this regime properly, we consider inextensible bers with a vanishing rotation inertia, neglect shearing, and assume moment strains to remain small making use of an elastic constitutive model while large displacements, at the origin of the desired geometric nonlinearities, are allowed. The model is thus strictly subject to nite 2 rotations about the cross-section axes (bending) and about the tangent of the centerline (twisting). The corresponding governing equations a set of partial dierential equations together with boundary conditions were rst developed by Kirchho and Clebsch in their theory of thin elastic rods under nite displacements [START_REF] Dill | Kirchho's theory of rods[END_REF]. Within a more general framework on shells, rods and points, the Cosserat brothers [1909] later on proposed a clever mathematical representation of the rod geometry, relying on a space curve (the centerline) together with a material frame attached to the rod cross-section and continuously rotating along the centerline about the so-called Darboux rotation vector. A modern description of these theories can be found in [START_REF] Antman | Nonlinear Problems of Elasticity[END_REF][START_REF] Audoly | Elasticity and Geometry: from hair curls to the nonlinear response of shells[END_REF]. [START_REF] Pai | Strands: Interactive simulation of thin solids using Cosserat models[END_REF] was the rst to introduce them to the Computer Graphics community.

2 As opposed to innitesimal.

14

Chapter 2. Discrete curvature-based models for dynamic bers Discretizing material rods In Mechanical Engineering, both nite dierences and nite elements approaches were developed to discretize material rods in space and time. Though nite dierences schemes have in principle the advantage of being easy to set up, properly accounting for the rod boundary conditions (typically, a clamped rod with the other end free) generally requires the use of a shooting strategy, which implies the solving of multiple nonlinear problems. Moreover, the sti nature of the Kirchho equations, stemming from the presence of fourth-order spatial derivatives, imposes the use of overly small steps in time and space, or sophisticated implicit integrators [START_REF] Goyal | Non-linear dynamic intertwining of rods with self-contact[END_REF]. In contrast, a nite elements strategy allows one to single out spatial terms from time-evolving quantities, and provides a vast choice of elements to approximate them together with the boundary conditions. A popular method is the so-called geometrically exact beam approach [START_REF] Reissner | One one-dimensional large-displacement nite-strain beam theory[END_REF][START_REF] Simo | A three-dimensional nite-strain rod model. part ii: Computational aspects[END_REF], which derives an exact weak formulation for a generalized Kirchho rod with stretching and shearing, and nally discretizes the displacement and rotation elds with interpolating shape functions. One important issue of this approach, which spurred many subsequent works in the nite elements community, deals with the proper interpolation of rotations for preserving objectivity, i.e., invariance of the strain measures under rigid motion [START_REF] Criseld | Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its nite-element implementation[END_REF]].

Moreover, regarding our specic needs here, this method is not directly applicable to the handling of inextensible and unshearable rods.

In Computer Graphics, nite dierences schemes initially proposed by [START_REF] Pai | Strands: Interactive simulation of thin solids using Cosserat models[END_REF] to solve the statics of Kirchho rods were subsequently superseded with more robust schemes so as to deal with the full dynamic case. Reduced -coordinates models, based on a minimal parametrization of the system, were proposed to account for the exact kinematics of the rod, and especially to preserve inextensibility: this was the case for the articulated rigid body approach [START_REF] Hadap | Modeling dynamic hair as a continuum[END_REF][START_REF] Hadap | Oriented strands -dynamics of sti multi-body system[END_REF], parameterized by angular joints, and then by our super-models (detailed hereafter), parameterized by curvatures and twist. In contrast, further work focused on nodal models in order to get an explicit, point based representation of the centerline leading to a sparse mass matrix, at the price of adding external constraints to preserve the true kinematics. In the CoRde model [START_REF] Spillmann | CoRdE: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects[END_REF], both positions and orientations are considered as degrees of freedom. The Lagrange equations of motion are written for discrete approximations of kinetic and potential energies including a stretch term and orientations are coupled back to the centerline through soft constraints. Relying on the Bishop frame, [START_REF] Bergou | Discrete elastic rods[END_REF] use a curve-angle parameterization to reduce the number of redundant parameters and guarantee that the orientation frame naturally remains adapted to the centerline. Discrete equations of motion are then established by leveraging principles from discrete dierential geometry. Due to the choice of a nodal parameterization, inextensibility however needs to be explicitly enforced, e.g., through a fast projection scheme [START_REF] Bergou | Discrete elastic rods[END_REF] or a sti stretch term [START_REF] Bergou | Discrete viscous threads[END_REF]. Finally, to ensure proper stability at an acceptable computational cost, an implicit scheme based on Newton's method is advocated to discretize the nonlinear sti bending and stretching forces [START_REF] Bergou | Discrete viscous threads[END_REF].

Kirchho equations for thin elastic rods 15

Reduced Lagrangian dynamics One advantage of reduced dynamics is that, by directly considering moment strains as actual degrees of freedom (i.e., material curvatures and twist instead of positions), the rod kinematics is exactly preserved, without redundancy and without adding any further constraint. Models parameterized by curvatures 3 also benet from an inexpensive implicit handling of bend- ing forces, as those forces are linear in curvature. Finally, while multiple collision tricks such as the position alteration technique [START_REF] Bara | Large steps in cloth simulation[END_REF] were specically developed for nodal models in Computer Graphics, more sophisticated, constraint-based frictional contact solvers naturally cope with reduced Lagrangian models (see Chapter 3) without having to worry about getting intermingled with external kinematic constraints.

Kirchho equations for thin elastic rods

Notation In what follows, s denotes the space variable and t the time variable.

Space derivatives are represented by the prime symbol, so that a (s, t) = ∂a ∂s and time derivatives by the dot symbol, so that ȧ(s, t) = ∂a ∂t . For the sake of clarity, we may omit the time variable when describing the geometry of the rod. The special orthogonal group of dimension 3, denoted SO(3), collects nite rotations of R3 (represented as direct orthogonal matrices) and is a non commutative Lie group.

We consider an inextensible and unshearable material rod of length L, represented by a centerline r(s) together with a material frame R(s), both parameterized by arc length s ∈ [0, L]. At location s, the vector r(s) ∈ R 3 gives the 3D position of the centerline and the rotation matrix R(s) ∈ SO(3) encodes the (unitary) tangent vector 4 n 0 (s) = r (s) as well as the two normal vectors n 1 (s) and n 2 (s) attached to the cross-section of the rod.

s = 0 s = L n 0 (s) n 1 (s) n 2 (s) r(s)
For simplicity, we assume that the rod is clamped at s = 0, and that its clamped position r(0) = r cl and orientation R(0) = R cl are given. We also assume that the end s = L is free, that is, its position or velocity is not prescribed kinematically.

Note that these two assumptions hold in most real bers we wish to model, e.g., plants and hair. Otherwise, any of them could easily be dropped out; the former, by releasing r cl and R cl as degrees of freedom; the latter, by adding some constraint at the s = L end.

Chapter 2. Discrete curvature-based models for dynamic bers Kinematics From s = 0 to s = L, the material frame R(s) continuously evolves along the centerline r(s) through innitesimal rotations about the so-called Darboux vector Ω(s) which represents the instantaneous space rotation vector of the rod. This space evolution mathematically writes

R (s) = [Ω(s)] × R(s), (2.1)
where [u] × denotes the skew symmetric matrix corresponding to the vector cross product operator, i.e., [u]

× v = u × v.
R (s) = R(s) [κ(s)] × .
(2.2)

Finally, by compacting the centerline and the material frame into one single variable F(s) = {r(s); R(s)} and assuming κ(s) is xed, the full kinematics of the rod can be formulated as an explicit 5 linear rst-order Cauchy-Lipschitz problem, referred to as the Darboux problem (see, e.g., [START_REF] Ivanova | On one approach to solving the Darboux problem[END_REF]),

F (s) = n 0 (s) ; R(s) [κ(s)] × with F(0) = {r cl ; R cl } as initial conditions, (2.3)
which admits a unique solution. Note that the ambient space is not a vector space but rather a nonlinear dierentiable manifold, since the kinematic relationship for the material frame operates onto the non commutative Lie group SO(3). Due to non commutativity, the solution has no formal expression in the general case.

Dynamics Let ρ be the volumetric mass of the rod and S the surface area of its cross-section. We assume the rod is subject to external forces such as gravity or contact forces. Expressing the balance of linear and angular momentums on an innitesimal portion of the rod and neglecting inertial momentum due to the vanishing cross-section leads to the following dynamic equations for a Kirchho rod,

T (s) + p(s) = ρSr(s) M (s) + n 0 (s) × T(s) = 0 (2.4a) (2.4b)
where p is the the sum of external forces distributed per unit length (for instance, in the case of gravity, p = -ρSge z ) and T(s) and M(s)) are the internal force and moment respectively, transmitted from the free part of the rod through its cross-section at s. Elastic constitutive law Dynamic equations are completed with a constitutive law that expresses the ability of the rod to bend and twist elastically, M(s) = K 3 (κ(s) -κ(s)) in the local basis R(s), (2.5) where K 3 = diag 3 (K 0 , K 1 , K 2 ) is a diagonal 3 × 3 matrix collecting the twisting and bending stinesses K 0 , K 1 , and K 2 respectively, and where κ(s) ∈ R 3 collects the intrinsic twist and curvatures of the rod, used to model natural curliness (i.e., the shape the rod would have in the absence of external forces). As we are modeling a homogeneous material, both the linear mass density ρS and the stiness matrix K 3 are assumed to be constant with respect to space and time 6 . In contrast, the intrinsic curvature vector κ(s) may vary spatially to account for a wide range of natural shapes.

Boundary conditions As mentioned ealier, we consider the rod to be clamped at end s = 0 and free at the other end s = L. Corresponding boundary conditions read

         r(0) = r cl Enforced clamped position (2.6a) R(0) = R cl Enforced clamped frame (2.6b) T(L) = 0
No external force at free end (2.6c)

M(L) = 0

No external torque at free end (2.6d)

Numerical model Equations (2.32.5) together with the boundary conditions (2.6) at s = 0 and s = L form a nonlinear and sti boundary value problem, which has no explicit solution is the general case and is known to be dicult to solve numerically.

Discrete curvature-based models

Realizing that curvature plays a key role in both the kinematics and the dynamics of the rod, an interesting idea consists in approximating the curvature vector with a simple, polynomial expression that is function of s. The coecients of the polynomial are then taken as primary variables of the discrete model. One immediate consequence is that bending forces, which are linear in curvature, become linear in the discrete variables. Being sti in nature, those forces can thus be treated implicitly in a straightforward manner, without having to solve a nonlinear problem.

Furthermore, the kinematic (Darboux) problem becomes numerically tractable, and even solvable with a closed-form solution in the piecewise constant case.

Chapter 2. Discrete curvature-based models for dynamic bers Discrete kinematics In the following, we shall denote by q ∈ R m our state variable collecting the discrete curvatures and twists 7 , and q the vector of same size storing the discrete intrinsic curvatures and twists. For now, we assume that the centerline of the rod r can be computed as a function of s, q, r cl and R cl , by solving the Darboux problem (2.3) with an accurate numerical method. The latter point is specically addressed in Sections 2.3.1 and 2.3.2 in the case of the super-helix and the super-clothoid model, respectively. Formally dierentiating the centerline twice with respect to time leads to the following expression for acceleration, r(s, t) = r * (s, t) + q(t) ∂ 2 r ∂q 2 (s, t) q(t) + ∂r ∂q (s, t)q(t),

(2.7)

where r * (s, t) is the acceleration generated by the clamping motion, which can be dropped when the clamped end is static. Expression (2.7) puts in evidence the linear dependence of the centerline acceleration r with respect to the generalized acceleration q, a property that will be leveraged when designing a linear-time integration scheme in Section 2.4.

Discrete dynamics Discrete equations of motion result from a weak formulation of the strong Kirchho equations (2.4), where the trial functions are deduced from the constrained, piecewise polynomial kinematics. Consider an innitesimal virtual displacement δq of our discrete degrees of freedom. This translates into a perturbation δκ(s) in the curvatures and twist function, which causes an innitesimal rotation of the material frame about a virtual rotation vector δθ(s), such that δR(s) = [δθ(s)] × R(s), as well as an innitesimal displacement δr(s) of the centerline. Applying the principle of virtual work [START_REF] Reissner | One one-dimensional large-displacement nite-strain beam theory[END_REF] Finally, relating perturbed quantities to the virtual displacement δq and using the kinematic equation (2.7) yields the discrete dynamic equations M(q)q + K (qq) + G(q) + A(q, q) = 0

(2.8) 7 Discrete curvature and twists q corresponds to per-element curvatures and twist for the superhelix model, and per-joint curvatures and twist for the super-clothoid model.

where

M(q) = ρS L 0 ∂r ∂q
∂r ∂q ds

G(q) = -ρS L 0 ∂r ∂q g ds A(q, q) = ρS L 0 ∂r ∂q
q ∂ 2 r ∂q 2 q + r * ds,

(2.9)

and where the stiness matrix K is a constant function linear in K 3 , with a pattern depending on the order of discretization (see Sections 2.3.1 and 2.3.2).

The main challenge consists in evaluating vectors G and A and matrix M in a both accurate and fast way. Sections 2.3.1 and 2.3.2 propose two dierent spatial discretizations one piecewise constant and one piecewise linear with respect to s and show how each of them addresses this problem.

Semi-implicit time-integration Finally, once the discrete dynamic equation (2.8) is assembled, we may use a semi-implicit Euler scheme to nd q at next time step, by simply taking linear terms as implicit and others as explicit, i.e., such that the time-discrete dynamic equation reads M(q t )q t+1 + K q t+1 -q + G(q t ) + A(q t , qt ) = 0.

(2.10)

As the elastic term of the dynamics appears to be the stier term of the equation, this cheap integration scheme proves sucient to guarantee a good stability of the algorithm in most situations, even those involving some strong motion of the rod.

Piecewise uniform curvature: the super-helix model

The super-helix element We discretize the rod into N elements with arc lengths at joints denoted by s i , i ∈ {0..N }. On each element E i between two successive joints s i and s i+1 , the curvature vector κ(s) is assumed to be constant with respect to arc length s, i.e., κ(s) = κ i ∀s ∈ [s i , s i+1 ].

The generalized coordinate vector q of the system thus collects the 3 N curvatures and twist values κ i,j with j ∈ {0..3}. The stiness matrix is simply dened as

K =          0 K 3 0 • • • 0 0 1 K 3 . . . . . . 0 . . . . . . 0 . . . . . . N -1 K 3 0 0 • • • 0 N -1 K 3          = K 3 ⊗ diag N ( ) , (
R(s) = R(s i ) exp s [κ i ] × and r(s) = r(s i )+ s i+1 s i R(s) e x with e x = [1 0 0] ,
which can be computed in closed-form, i.e., using only basic and cheap unitary operations. In our case, the above expressions can be evaluated using only cosine and sine functions, and it is noteworthy that the centerline r(s) takes the form of a perfect circular helix on each element E i . Starting integration at s = 0 with clamping conditions (2.6a 2.6b) and then using continuity conditions for R(s) and r(s) at each joint s i , we can integrate the full kinematics of the rod from s = 0 to s = L. The super-helix model then features a C 1 -smooth piecewise helical centerline, and all the terms of the kinematics (material frame, centerline, but also deformation gradients), and nally all the terms of the dynamics (2.8)), can be computed in closed-form 8 .

Limitations A discrete rod with piecewise constant curvature may still represent a rather rough approximation of the continuous case, with an improper degree of continuity at the joints (see Figure 2.2). Instead of using an excessively rened primitive, one may think it would be worth designing a richer, higher-order element with linear curvatures, that would better stick to the actual curvature prole of real bers and guarantee visually pleasing smoothness (i.e., at least C 2 -continuity) of the centerline at any resolution. One becomes unfortunately faced with the loss of a formal expression for the kinematics. Yet, observing that the Darboux problem still possesses a lot of structure, we show in the following that such a space clothoid element can be conveniently derived. The key is to introduce a fast and accurate integration scheme based on an adaptive power series summation algorithm. This numerical algorithm is then used as a formal computation tool to evaluate the spatial terms of the dynamics at a high precision.

Piecewise linear curvature: the super-clothoid model

In the piecewise linear case, where each element takes the form of a so-called 3D clothoid 9 , we were able to build an accurate integration scheme which proved to be 8 Computing each term of the dynamics (2.8) amounts to evaluating the integral over [0, L] of cosine and sine functions, and products of cosine and sine functions.

9 A clothoid or Euler spiral is a planar curve characterized by a linear curvature w.r.t. the arc length. By extension to 3D, the term 3D Euler spiral was coined by [START_REF] Harary | 3D Euler spirals for 3D shape completion[END_REF] to name a space curve characterized by linearly geometric curvature and torsion. In our case we rather use the term space clothoid to name the centerline of a rod characterized by linear material curvatures and twist. The centerline of our rod element is actually more general than the so-called 3D Euler spiral, the entire class of so-called 3D Euler spirals being obtained by cancelling the rst material curvature. Note that this dierence will be discussed again in Chapter 4. orders of magnitude faster compared to classical integration methods [START_REF] Casati | Super space clothoids[END_REF]. The key of our approach was to leverage the form of the solution as a power series expansion, while avoiding the pitfall of catastrophic cancellation through an adaptive integration strategy. With this tool in hand, we were able to demonstrate that the super-clothoid model could capture intricate shapes both robustly and eciently, with better spatial accuracy and geometric fairness compared to state-of-the-art methods (see Figure 2.1).

The super-clothoid element As before we discretize the rod into N elements.

However, now the discrete curvature variables κi are located at joints s i , and the curvature vector q collects 3 (N + 1) degrees of freedom. On each element E i of length i , the curvature vector κ(s) is assumed to vary linearly with arc length s, so that its expression reads

κ(s) = 1 - s -s i i κi + s -s i i κi+1 ∀s ∈ [s i , s i+1 ].
The stiness matrix is not diagonal anymore but becomes tridiagonal,

K =          0 3 K 3 0 6 K 3 0 • • • 0 0 6 K 3 0 + 1 3 K 3 . . . . . . . . . 0 . . . . . . . . . 0 . . . . . . . . . N -2 + N -1 3 K 3 N -1 6 K 3 0 • • • 0 N -1 6 K 3 N -1 3 K 3          = K 3 ⊗ tridiag N ( ) ,
(2.12)

where tridiag N ( ) is a tridiagonal matrix of size N linearly depending on the lengths.

Chapter 2. Discrete curvature-based models for dynamic bers Kinematic integration using power series expansion We rst focus on a single element with index i = 0. When the curvature vector is polynomial of degree ≥ 1, the key idea is to formulate the solution of the Darboux problem (2.3) as a power series expansion (PSE). Indeed, invoking Cauchy's theorem on analytic solutions of linear ODEs with analytic coecients (see, e.g., [Poole, 1936, 2]), we deduce that the solution F of (2.3) is C ∞ and admits a power series expansion

F(s) = ∞ n=0 r n s n ; ∞ n=0
R n s n on R. In our particular case where the curvature is a polynomial of degree 1 (space clothoid element), the general term of the series

rn (s) = r n s n ; Rn (s) = R n s n is recursively dened as                    R0 = R cl R1 (s) = s R0 [κ 0 ] × (2.13a) Rn+2 (s) = s n + 2 Rn+1 (s) [κ 0 ] × + s Rn (s) [γ] × ∀n ∈ N r0 = r cl rn+1 = s n + 1 Rn e x ∀n ∈ N, (2.13b)
where 0 is the length of the element, κ0 and κ1 are its initial and nal curva- tures respectively, and γ is the curvature slope, γ = κ1 -κ 0 0 . Note that according to (2.13b), computing the centerline follows from that of the material frame, the latter involving the recursive sequence (2.13a) of second order.

Furthermore, one important property of the above power series expansion is that only the very rst terms of the series are relevant, the following ones rapidly decreasing in norm and falling below machine precision. This means that the solution F can be evaluated accurately through a truncated series containing a small number of terms (in practice, around 100 terms). This nice property is due to the simple structure of our kinematic problem (2.3), which formulates as an explicit linear ODE with polynomial coecients. In this case indeed, we can prove that the general term of the series decreases super-linearly to zero as n tends to innity [START_REF] Neher | An enclosure method for the solution of linear ODEs with polynomial coecients[END_REF].

Catastrophic cancellation issue Unfortunately, when numerically evaluating the sum of the relevant terms in nite precision, one is inevitably faced with roundo issues leading to huge approximation errors. Consider expression (1 + y) -y which should be equal to 1 whatever the value of y. In oating-point arithmetic, this equality only holds if y is close enough to 1. In double precision for example, take y = 10 16 and compute the expression above. The numerical result is 0.0, yielding a relative error of 100%. This error is the consequence of rst, an absorption phenomenon when computing the sum 1+10 16 , which, due to machine overow when aligning mantissa, is approximated to 10 16 . Then, a cancellation phenomenon when subtracting 10 16 . Such unfortunate combination of absorption and cancellation leads to erroneous results and for this reason is called catastrophic cancellation. Details on oating-point arithmetic can be found in, e.g., [START_REF] Goldberg | What every computer scientist should know about oatingpoint arithmetic[END_REF].

(a) (b) (c)
Figure 2.3 2.3a: Dramatic loss of precision when naively summing power series of the kinematics. 2.3b: In contrast, our piecewise summation algorithm guarantees high precision of the summation. 2.3c: A long and highly curved space clothoid integrated with our piecewise computation method, using 109 subdivisions.

In our case, catastrophic cancellation problems precisely occur when naively computing the sum of our power series Rn (s) for a long and/or curly rod (see Figure 2.3a). To explain this phenomenon, we have plotted in Figure 2.4 (in blue)

the norm of the general term Rn (s) as a function of n, for dierent values of s. The resulting hillock-like prole implies that when computing the sum of the series, one actually adds very small values together with very large ones in norm, the widest range being obtained when getting to the top of the hillock. Note that the larger s, the higher the top of the hillock. More precisely, we have shown in [START_REF] Casati | Super space clothoids[END_REF]] that the top of the hillock H(s) grows quasi-exponentially with the increasing function λ(s) dened as

λ(s) = 2s( κ0 ∞ + s γ ∞ ).
Moreover, looking back at Recursion (2.13a), one notes that entries of the matrices to be added are of alternating sign, due to the product with skew symmetric matrices.

This results in cancellation when computing the sum. All this combined together, it is then not surprising that we are faced with a catastrophic cancellation issue when λ(s) (and thus s) becomes too large. As λ(s) increases with s as well as with intrinsic curvatures, we now understand why numerical issues show up for a long and/or curly rod. Luckily enough, as shown in the following, we can avoid catastrophic cancellation by devising an adaptive summation algorithm. As expected [START_REF] Neher | An enclosure method for the solution of linear ODEs with polynomial coecients[END_REF], the decreasing towards 0 appears to be super-linear. In red: Evolution, function of s, of the upper-bound H(s), in log scale. Note that the plot of this upper-bound visually matches the maximum function max n Rn (s) ∞ , meaning that the top of the hillock grows quasi-exponentially with s.

Fast and Precise Power Series Summation

should be bounded by 10 d 2 so as to be able to safely cover additions on the range

[10 -d 2 , 10 d 2 ]. We have proved that a sucient upper-bound for M is M max n ∈ N s.t. (n + 1) n 10 d 2 n! .
(2.14)

One obtains M 19 for d = 16. In practice, we set M to 10 to maintain good precision across summation. This choice allowed us to reach high precision for all the summations we have computed.

Adaptive piecewise summation Our full summation method then relies on an automatic subdivision of the integration domain into subintervals, on which integration can be safely performed. More precisely, our adaptive summation algorithm consists in splitting the evaluation domain

[0, 0 ] into p adaptive subintervals [0, σ 1 ], [σ 1 , σ 2 ], • • • [σ p-1 , 0 ] such that σ i+1 = σ i + s max (σ i ), where s max (σ i ) is pro- vided by this upper-bound function, s max (σ) =        √ κ(σ) 2 ∞ +2M γ ∞ -κ(σ) ∞ 2 γ ∞ if γ = 0 M 2 κ(σ) ∞ else if κ(σ) = 0 +∞ otherwise.
(2.15) On each subinterval, summation is thus guaranteed to be performed with good accuracy (see an illustration in Figure 2.5).

[ We show in [START_REF] Casati | Super space clothoids[END_REF] than our adaptive piecewise summation algorithm reaches a good precision (up to the machine precision) orders of magnitude faster compared to classical integration schemes such as Euler or Runge Kutta of order 2 and 4.

0, σ 1 ] [σ 1 , σ 2 ] [σ 2 , ]
Finally, using continuity conditions at joints, the kinematics of the full rod can be accurately integrated from s = 0 to s = L.

Propagation to the terms of the dynamics To compute all the terms of the discrete dynamic equation (2.8), we need to perform several operations on our kinematic terms, such as linear combination, scalar product, dierentiation with respect to q, or integration with respect to s. All these operations are actually easy to perform using power series expansion; furthermore, we show that our precision guarantee resists all these operations, thus providing us with a fast and accurate method for computing all the terms of the dynamics.

Model comparison

We have compared the super-clothoid model against the super-helix model and the discrete elastic rod model [START_REF] Bergou | Discrete elastic rods[END_REF], by per-

Linear-time integration of the dynamics

One major advantage of the aforementioned curvature-based formulations is that the kinematics of the discretized rod remains, by construction, perfectly inextensible. Such intrinsic inextensibility thus removes the burden of adding subsequent (sti ) inextensibility contraints when solving the dynamics. The price to pay, however, is that the geometry of the rod is not readily available but has to be computed iteratively from the curvatures. Moreover, the inertia matrix of such systems turns out to be dense, leading to a N 2 complexity when solving them at each time step, where N is a number of elements. Exploring high-order discrete models such as the clothoidal element helped us reduce the number of elements required for a given 2.4. Linear-time integration of the dynamics precision; however, for a very long rod requiring high spatial precision, the N 2 computation would remain a bottleneck.

One strategy to deal with the latter issue was to propose a linear, recursive scheme for curvature-based discrete models, inspired by the Featherstone algorithm for articulated chains of rigid bodies [START_REF] Featherstone | The calculation of robot dynamics using articulatedbody inertias[END_REF]. We have investigated such a linear strategy in [START_REF] Bertails | Linear time super-helices[END_REF] in the case of a super-helix, but our algorithm could be similarly applied to higher-order discrete models such as the super-clothoid model. Below we explain the general principle of the method.

Recursive integration scheme

Non-local dynamics We consider again the time-and-space discrete dynamic equations (2.10) for a super-model, which are recalled here, with the two unknowns (implicit variables) marked at time t + 1,

M(q)q t+1 non-local + K q t+1 -q local + G(q) + A(q, q) = 0.
(2.16) Let E i be the i th element of the rod. While the i th block of the elastic term is local to the element and only depends on the i th block of generalized coordinates q i , the block of the inertial term depends on all generalized accelerations q, as M(q) is obviously dense from its expression in (2.9), knowing that the term ∂r ∂q evaluated at s ∈ E i depends on all previous coordinates q j,j≤i due to the recursive kinematics of super-models.

Building local dynamic equations from tip to root

The key of our approach, detailed in [START_REF] Bertails | Linear time super-helices[END_REF], is to show that Equation (2.16) can be rewritten locally on each element E i as

K i (q i -qi ) = A i qi + R i ri (0) + C i θ θ θ i (0) + b i , (2.17)
where ri (0) and θ θ θ i (0) are the linear and angular accelerations of the rst end of element E i , respectively, and the 3 × 3 matrices A i , R i and C i are called cumulate inertias, and the vector b i cumulate force of the rod.

As in [START_REF] Featherstone | The calculation of robot dynamics using articulatedbody inertias[END_REF], the cumulate inertias and the cumulate force of the rod can be computed recursively from the free end (i = N -1) to the clamped end (i = 0) of the rod. This can be proved by induction, using the kinematic property (2.7) according to which the linear (and also angular) acceleration of the rod is linear w.r.t. the generalized acceleration q.

Recursive solving from root to tip Finally, once each local dynamic equation (2.17) is built, one may simply integrate the dynamics recursively from root to tip, using continuity conditions at the joints between elements.

Chapter 2. Discrete curvature-based models for dynamic bers

Comparison with the composite approach

The recursive approach presented above brings several advantages compared to the so-called composite approach detailed in Section 2.3, in which full terms of the dynamics are assembled.

First, we have observed a computational gain of about ×20 for rods made of more than 20 elements. In practice, this allowed us to simulate very long rods in real-time (see Figure 2.7(a)). Then, by simply adapting boundary conditions at the joints between elements, we were able to extend the recursive approach quite straightforwardly to the handling of tree-like structures (see Figure 2.7(b)). Note that tree-like structures would have been much less trivial to model using the composite approach.

On the downside, the recursive approach may require a smaller timestep to be stable. This is due to the explicit computation of qi in the rst pass of our algorithm, which makes it slightly less stable than the composite method for an equivalent time step. In the case of articulated rigid bodies, the explicit solving of Featherstone's algorithm was pointed out by [START_REF] Hadap | Hair Simulation[END_REF] who proposed an interesting framework for enforcing an implicit solving [START_REF] Hadap | Oriented strands -dynamics of sti multi-body system[END_REF]. Note however that in our case, the dierence in stability between the two methods keeps on being small.

Finally, we shall not use the recursive formalism of super-models in the next chapter, which is dedicated to the introduction of frictional contact for which composite models are better suited for. However, we believe that the recursive model would be particularly useful for modeling botanical or biological structures. In the latter case, a promising extension would be to design an adaptive scheme as in [START_REF] Redon | Adaptive dynamics of articulated bodies[END_REF], by relying upon Featherstone's divide-and-conquer algorithm [Featherstone, 1999a,b]. So far, to the best of our knowledge, such adaptive algorithms were leveraged for articulated rigid bodies only. Adaptive super-models could be very useful for manipulating in real-time long and complex biological laments (e.g., macro-molecules) prone to bending and twisting, by automatically tuning their number of active degrees of freedom. We have investigated new discrete models for solving the Kirchho dynamic equations for thin elastic rods. All our models share a curvature-based spatial discretization, allowing them to capture inextensibility of the rod intrinsically, without the need for adding any kinematic constraint. Moreover, elastic forces boil down to linear terms in the dynamic equations, making them well-suited for implicit integration. We shall furthermore show in Chapter 4 that this property will be leveraged for solving inverse design problems.

Interestingly, our discretization methodology can be interpreted from two dierent points-of-view. From the nite-elements point-of-view, our strain-based discrete schemes can be seen as discontinuous Galerkin methods of zero and rst orders.

From the multibody system dynamics point-of-view, our discrete models can be interpreted as deformable Lagrangian systems in nite dimension, for which a dedicated community has started to grow recently [START_REF] Sugiyama | Flexible multibody dynamics Essential for accurate modeling in multibody system dynamics[END_REF]. We note that adopting the multibody system dynamics point-of-view helped us formulate a lineartime integration scheme, which had only be investigated in the case of multibody rigid bodies dynamics so far.

Towards a curvature-based shell model Recently, we have started to investigate similar high-order modeling strategies for surfaces, in particular for the case of inextensible and developable shells. In [Blumentals et al., 2016a] we have built The second core topic of my research activity deals with the modeling and simulation of dry (Coulomb) frictional contact within (potentially large) collections of thin structures, such as rigid grains, thin elastic rods or plates/shells. Target applications include the simulation of brous and granular materials (such as human hair and sand, respectively), and more recently, cloth. For all these dynamical systems, dry frictional contact plays a major role in their visual appearance (see Figure 3.1).

Summary of contributions Our initial motivation for conducting research in frictional contact algorithms was to improve the realism of hair simulations, where self-interactions were simply neglected in the worse case [START_REF] Anjyo | A simple method for extracting the natural beauty of hair[END_REF], and regularized in the best case [START_REF] Mcadams | Detail preserving continuum simulation of straight hair[END_REF], yielding an overly smooth appearance of hair. We have started our study by modeling hair as a collection of individual thin elastic rods interacting through contact and friction a method known as the discrete element method. However, instead of using common compliant and regularized models for contact and friction, we have followed the nonsmooth contact dynamics approach developed by Moreau and Jean a few decades ago [START_REF] Moreau | Unilateral contact and dry friction in nite freedom dynamics. Nonsmooth mechanics and applications[END_REF][START_REF] Moreau | Numerical treatment of contact and friction: the contact dynamics method[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF], in which emerging eects of 32 Chapter 3. Frictional contact models for large multibody systems We have rst studied the numerical counter-parts for each main formulation and evaluated how each of them performs in terms of eciency, robustness, and scalability [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in ber assemblies[END_REF]. Facing robustness and scalability issues when using state-of-the-arts approaches, we have then designed a new frictional contact solver by coupling an iterative Gauss-Seidel strategy together with an extremely robust one-contact solver relying upon a hybrid local solver [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. This new solver proved to converge well in scenarios involving thousands bers subject to tens thousands frictional contact points, and thus allowed us to enhance considerably the realism of hair simulations (see Figure 3.8). Recently, we have adapted this solver for coping with frictional contact in cloth modeled as a nodal system [START_REF] Daviet | Fast Cloth Simulation with Implicit Contact and Exact Coulomb Friction[END_REF] (see Figure 3.9).

In parallel to our work on the discrete element modeling of ber assemblies, we have started to investigate continuum modeling strategies so as to scale up simulations to hundreds thousands bers the actual order of magnitude of a full head of hair. Faced with the glaring lack of macroscopic models for drybrous materials in the literature, we have started this long-term investigation by considering the simpler case of granular matter, where each element is not a deformable ber but a rigid and isotropic grain. We have modeled a large assembly of 3.1. Prior art and choice for nonsmooth models 33 grains as a dilatable viscoplastic uid combined with a Drucker-Prager yield criterion (which stands for a generalization of the Coulomb friction law to a continuum) and a unilateral compressibility constraint, rst in the dense regime [START_REF] Daviet | Nonsmooth simulation of dense granular ows with pressure-dependent yield stress[END_REF], then accounting for a varying volume fraction [Daviet and Bertails-Descoubes, 2016a]. Thanks to a semi-implicit time-stepping scheme and a careful spatial discretization of our rheology, we were able to preserve at each time step the exact coupling between normal and tangential stresses, in a stable way. This contrasts with previous approaches that either regularize or linearize the yield criterion for implicit integration, leading to articial sliding motion or visible grid artifacts [START_REF] Chambon | Numerical simulations of granular free-surface ows using smoothed particle hydrodynamics[END_REF][START_REF] Narain | Free-owing granular materials with two-way solid coupling[END_REF]. Remarkably, our discrete problem turned out to be very similar to the discrete frictional contact problem (DFCP) that we had formerly encountered when dealing with discrete element modeling, which allowed us to leverage again our robust nonsmooth solver [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. Our method was qualitatively validated by successfully capturing typical macroscopic features of some classical experiments, such as the discharge of a silo or the collapse of a granular column.

Organization of the chapter Section 3.1 motivates our choice for modeling frictional contact using rigid and nonsmooth frictional contact laws, in the context of both multibody and continuum dynamics. Section 3.2 rst formulates the discrete element problem where a large number of deformable structures are interacting through non-penetration constraints combined with the nonsmooth Coulomb friction law, then spells out the specicity of our resulting discrete problem (DFCP), and discusses robust and ecient numerical solutions for solving it. Following a similar outline, Section 3.3 presents the nonsmooth Drucker-Prager law for pressuredependent yield-stress uids, combined with a non-compressibility condition, and provides insights into the semi-implicit discretization of the full dynamical problem, which turns out to share a structure similar to that of our DFCP.

Prior art and choice for nonsmooth models

Most popular approaches in Computer Graphics and Mechanical Engineering consist in assuming that the objects in contact are locally compliant, allowing them to slightly penetrate each other. This is the principle of penalty-based methods (or molecular dynamics), which consists in adding mutual repulsive forces of the form k f (δ), where δ is the penetration depth detected at current time step [START_REF] Cundall | A computer model for simulating progressive large scale movements of blocky rock systems. in proceedings of the symposium of the international society of rock mechanics[END_REF][START_REF] Moore | Collision detection and response for computer animation[END_REF]. Though simple to implement and computationally ecient, the penalty-based method often fails to prevent excessive penetration of the contacting objects, which may prove fatal in the case of thin objects as those may just end up traversing each other (see Figure 3.2(a)). One solution might be to set the stiness factor k to a large enough value, however this causes the introduction of parasitical high frequencies and calls for very small integration steps [START_REF] Bara | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF]. Penalty-based approaches are thus generally not satisfying for ensuring robust contact handling.

34 Chapter 3. Frictional contact models for large multibody systems In the same vein, the friction law between solid objects, or within a yield-stress uid (used to model foam, sand, or cement, which, unlike water, cannot ow beyond a certain threshold), is commonly modeled using a regularized friction law (sometimes even with simple viscous forces), for the sake of simplicity and numerical tractability (see e.g., [START_REF] Spillmann | An adaptive contact model for the robust simulation of knots[END_REF][START_REF] Frigaard | On the usage of viscosity regularisation methods for visco-plastic uid ow computation[END_REF]).

Such a model cannot capture the threshold eect that characterizes friction between contacting solids or within a yield-stress uid. The nonsmooth transition between sticking and sliding is however responsible for signicant visual features, such as the complex patterns resting on the outer surface of hair, the stable formation of sand piles, or typical stick-slip instabilities occurring during motion (see Figure 3.1).

After having employed such regularized strategies during my PhD, the search for a realistic, robust and stable frictional contact method encouraged me to depart from those, and instead to focus on rigid contact models coupled to the exact nonsmooth Coulomb law for friction (and respectively, to the exact nonsmooth Drucker-Prager law in the case of a uid), which better integrate the eects of frictional contact at the macroscopic scale. Those modeling choices are summarized in Figure 3.3.

For stability and consistency purposes, nonsmooth models require implicit integration schemes 1 when being simulated numerically. Resulting discrete problems then take the form of algebraic equations subject to complementarity set-valued constraints, which are better understood and manipulated using convex analysis tools, the latter having been mainly developed in the early 60's by Jean-Jacques Moreau and R. Tyrrel Rockafellar. Such nonsmooth methods being disregarded in Figure 3.3 Summary of our nonsmooth modeling choices for dealing with frictional contact (right column), and comparison with their most simple regularized counterparts (left column). Notation: h(q) is the gap function between the two contacting objects, r N and r T are the normal and tangential components of the contact force, respectively, u T is the tangential component of the relative velocity between the two contacting objects, σ T the tangential stress tensor of the continuum, and εT the strain rate tensor.

h(q) r N (a) Compliant contact model h(q) r N (b) Rigid contact model
most Computer Graphics groups 2 a decade ago (in favor of regularized or explicit 2 It is interesting however to note that nonsmooth contact mechanics have been introduced in Computer Graphics by David Bara as soon as in the late 80's [START_REF] Bara | Analytical methods for dynamic simulation of non-penetrating rigid bodies[END_REF][START_REF] Bara | Coping with friction for non-penetrating rigid body simulation[END_REF][START_REF] Bara | Issues in computing contact forces for non-penetrating rigid bodies[END_REF][START_REF] Bara | Fast contact force computation for nonpenetrating rigid bodies[END_REF].

But it has not been followed up by the community at that time, then considerably losing ground on the Mechanical Engineering community in this area. In Computer Graphics, nonsmooth mechanics has slowly regained popularity in the context of rigid body dynamics, notably with the work of Kauman and colleagues [2005;2008]. methods), their assessment was my main motivation for integrating in 2007 the Inria project-team BiPop, specialized in the implicit treatment of frictional contact within the nonsmooth mechanics framework [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF].

Nonsmooth contact dynamics for thin elastic objects

In this section we present a discrete element model for an assembly of thin elastic objects (typically, bers), where frictional contacts are modeled and time-discretized using the nonsmooth contact dynamics framework [START_REF] Moreau | Numerical treatment of contact and friction: the contact dynamics method[END_REF]. The resulting Discrete Frictional Contact Problem (DFCP) is specically solved for by a new robust and scalable solver in the case of a large assembly of rods (≈ 2000) modeled as super-helices [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF], and for which state-of-the-art nonsmooth solvers were shown to fail. We have also started to adapt our solver to the interaction handling between large nodal systems and rigid bodies, such as a full garment worn by a walking character.

Adding unilateral contact with Coulomb friction

As individual components, we consider discrete Lagrangian systems; that is, spatially discretized systems whose kinematics is fully described by the choice of a nite set of reduced (or generalized) coordinates, without the need for adding any constraint. This is typically the case of the discrete curvature-based models derived for inextensible rods in Chapter 2, but also of the nodal systems classically used for modelling extensible plates and shells in Computer Graphics [START_REF] Bara | Large steps in cloth simulation[END_REF][START_REF] Grinspun | Discrete Shells[END_REF].

Contact-free dynamics Let q ∈ R m be the m generalized coordinates of such a spatially discrete system and v = dq dt its generalized velocities. Equations of motion simply read as an ODE,

M(t, q) dv dt = f (t, q, v), (3.1)
where the inertia matrix M(t, q) is symmetric positive denite, and the generalized force f collects all forces applying onto the system, including internal elastic forces, nonlinear inertial forces, and external forces such as gravity or air friction. Now, when considering multiple such systems simultaneously by concatenating their respective contributions, the structure of Equation (3.1) is globally preserved.

If we look more into details this will be relevant when designing solving algorithms in Section 3.2.4 we note that the inertia matrix M naturally becomes sparse-block as several individual components are added. In the remainder of the text, we shall equally refer to Equation (3.1) for the unconstrained dynamics of an individual or multiple components, keeping in mind the sparse pattern of M in the case of multiple components.

Nonsmooth contact dynamics for thin elastic objects 37

Local contact basis When coupling one or various components together with contact and friction (including contact with an external object, mutual contacts between components, and self-contact), we make several important simplifying assumptions. First, we only consider a nite set of contact points. Such a simplication is debatable, as it is well-known that some contacting congurations, for instance a plectoneme (see Figure 3.2(b)) may exhibit a continuous line of contact beyond a certain number of applied turns, a so-called supercoiled conguration [START_REF] Neukirch | Extracting DNA twist rigidity from experimental supercoiling data[END_REF]. However, concentrating contact forces into an adequate 3 set of points often proves sucient to avoid penetration and capture the Coulomb friction law, while considerably simplifying the numerical settings. We further assume that one contact always occurs between exactly two bodies, and that the surface of contacting objects is suciently smooth so that we can dene a tangent surface S and a contact normal e, from which we build a local basis (see Figure 3.4). We can then express in this basis the relative velocity u between the two bodies, and the force r at the contact point applying from one body to the other. Non-penetration conditions The rst law that we wish to enforce is a strict non-penetration constraint between the two objects. Let h(q) be the gap function between the two bodies, which corresponds to a signed distance function. We assume impacts to be purely inelastic (that is, kinetic energy instantaneoulsy dissipates when the two objects collide), and thus enforce a vanishing post-impact normal relative velocity u N (t + ). When the two objects are in contact at time t (i.e., such that h(q) = 0), some velocity-level condition, named Signorini condition, can be formulated to express the non-penetration constraint at time t>0,

0 ≤ u N ⊥ r N ≥ 0, (3.2)
where the x ⊥ y notation means that the two scalar variables x and y should be orthogonal, i.e., x y = 0 (more generally, in R d , x ⊥ y means x y = 0). The complementarity condition (3.2) means that either the contact ceases (u N > 0) and necessarily the normal contact force vanishes (r N = 0), or the contact persists (u N = 0) and necessarily the normal contact force, if not vanishing, becomes (or remains) active (r N > 0) to prevent penetration.

3 The set of contact points is determined on the y using a collision detection process.

38 Chapter 3. Frictional contact models for large multibody systems Dry (solid) friction law The dry friction law, often attributed to [START_REF] Coulomb | Théorie des machines simples: en ayant égard au frottement de leurs parties et à la roideur des cordages[END_REF] (but which a number of other authors have contributed to, such as de Vinci, Amontons, and Euler), results from the macroscopic observation that the normal component r N and the tangential component r T of a contact force r between two solids satisfy the following law,

r T ≤ µ S r N if u T = 0 (static regime) r T = µ D r N if u T = 0 (dynamic regime), (3.3)
where µ S and µ D are the static and dynamic friction coecients, respectively. Although it has been observed that µ S is generally slightly higher than µ D , we shall assume in the following, for the sake of simplicity, that they are equal and correspond to the same parameter µ. The value of µ shall depend on the characteristics of the contacting surfaces (from µ i = 0 for perfect contact without friction to µ i = µ max for rough surfaces).

One may observe that the Coulomb law is a threshold model, since the tangential contact force should reach a certain magnitude before the objects are able to slide on top of each other. Moreover, this magnitude directly depends on the normal contact force, meaning that a heavy object requires some higher pulling force to be dragged on a surface compared to a light object. This normal-dependent threshold, which is key to realism 4 , is also responsible for the fact that Coulomb's law does not t in with the so-called associated rules or generalized standard materials [START_REF] De Saxcé | The bipotential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF]], which are easier to manipulate formally and to handle numerically.

As a practical consequence, it will not be possible to view the Coulomb friction constraint as the optimality condition of a convex minimization problem (see 3.2.3).

A full discussion about the inclusion of Coulomb's law in the family of implicit standard materials (non-associated ow rules) is provided in [START_REF] Daviet | Modèles et algorithmes pour la simulation du contact frottant dans les matériaux complexes : application aux milieux breux et granulaires[END_REF].

The full Signorini-Coulomb disjunctive law Combining non-penetration constraint (3.2) together with Coulomb friction (3.3) is commonly referred to as the Signorini-Coulomb law. It can be written in a compact and disjunctive manner, as

• Taking-o case: r = 0 and u N ≥ 0

• Sticking case: r ∈ K µ and u = 0 • Sliding case: r ∈ ∂K µ , u N = 0 and ∃α ∈ R * + , u T = -αr T ,
where K µ is the second-order cone (SOC), dened in dimension d as

K µ := x = [x N , x T ] ∈ R × R d-1 , µx N ≥ x T ,
4 Indeed, imagine a stack of ber layers, or a sand pile under gravitational eld. If the friction threshold is constant, sliding will occur at any arbitrary depth of the pile, leading to an unrealistic avalanching behavior. In contrast, if the threshold is correctly set proportional to the normal force, the top surface will be, as expected, more prone to sliding avalanches compared to deeper layers.

with x N the (scalar) normal component of x ∈ R d , and x T its tangential component.

The Signorini-Coulomb law, illustrated in Figure 3.5, can be interpreted as follows: if the rst case, the normal relative velocity is nonnegative and the contact breaks (taking-o case). Then, in the case when the two bodies are in resting contact, either they are sticking to each other and the contact force can lie anywhere in K µ , or they are sliding with respect to each other, then the contact force must belong to the boundary of K µ and the tangent force must be collinear to the relative velocity with the opposite direction (according to the maximum dissipation principle [START_REF] Moreau | Unilateral contact and dry friction in nite freedom dynamics. Nonsmooth mechanics and applications[END_REF]). This disjunctive formulation is intuitive, but not very practical because of its combinatorial nature (there are 3 n cases to check, if the system contains n contacts).

In Section 3.2.3, we give equivalent formulations of Signorini-Coulomb's law which prove to be more tractable numerically in the general case.

In the remainder of this document we shall denote by C µ the set of vectors (r, u) in R d × R d satisfying Signorini-Coulomb's law.

The Discrete Frictional Contact Problem (DFCP)

We now assume that our system (3.1) is subject to n frictional contacts, supposed to be purely inelastic. To simplify notation, the n relative velocities u i and contact forces r i are collected into u and r, respectively. We dene H(q) := ∂u ∂v (q) the deformation gradient matrix 5 of dimension (n d, m) relating the relative spatial velocities at contact points u to the generalized velocities v. Let u f be the value of u when v = 0, which can be non-zero in the case of forced motion 6 . The nonsmooth 5 For a super-model, the deformation gradient matrix reads H(q, s) = ∂r ∂q (s) and can be evaluated eciently at any contact point s = sc, either formally (super-helix) or using adaptive power series summation (super-clothoid). For a nodal system, the gradient deformation matrix simply consists of identity blocks.

6 Unlike nodal systems where u f = 0, for a super-model this term may be nonzero when the clamped end is moving; in that case an extra velocity ṙ * (s), independent of v = q, shows up, and we have for each contact block u f = ṙ * (sc), where sc is the location of contact. 40 Chapter 3. Frictional contact models for large multibody systems system including contact forces r can be written as

   M(q) dv dt = f (t, q, v) + H(q) r u = H(q)v + u f (t, q) ∀i = 1 . . . n, (r i , u i ) ∈ C µ i . (3.4)
The notation dv dt is actually misused as the inclusion of contacts implies that we might observe jumps in the generalized velocities v. To circumvent this diculty, we follow Moreau's time-stepping scheme and integrate the equations over an arbitrary timestep δt [START_REF] Moreau | Unilateral contact and dry friction in nite freedom dynamics. Nonsmooth mechanics and applications[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]. We thus get a formulation over mechanical impulses instead of forces, and a velocity jump instead of an acceleration.

Let us set r := δt r dt, and let v and u be the discrete approximations of the generalized velocity and of the relative velocity, respectively, at the end of the timestep. The constraint-free dynamics (3.1) is discretized using a θ-scheme and can be written at a given timestep as Mv = f . We further assume that H(q) and u f (t, q) are almost constant during the timestep, and approximate them with H and u f , respectively. The discretized velocity/impulse formulation of (3.4), with unknowns (v, u, r), then reads

   Mv -f = H r u = Hv + u f ∀i = 1 . . . n, r i , u i ∈ C µ i .
(3.5) Problem (3.5), referred to as the one-step problem in [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF], is the core problem we wish to solve in this chapter. In the remainder of this document, we shall call it Discrete Frictional Contact Problem (DFCP). Note that we may eliminate v in (3.5) by introducing the Delassus operator W := HM -1 H and by setting b := u f + HM -1 f , and nally obtain the following reduced DFCP,

Wr + b = u ∀i = 1 . . . n, r i , u i ∈ C µ i , (3.6)
which may again be compacted as a force (or impulse) based problem,

∀i = 1 . . . n, r i , (Wr) i + b i ∈ C µ i . (3.7)
There is a number of diculties when attempting to solve the DFCP:

• First, the inclusion of (r i , u i ) in the set C µ i makes the problem nonconvex and nonsmooth. Standard convex optimization algorithms, which have wellestablished convergence guarantees, can thus not be applied to our problem, and specic algorithms need to be developed.

• Then, the DFCP may not always possess a solution 7 . One typical example is the (discrete) Painlevé paradox 8 , described for instance in [START_REF] Bara | Coping with friction for non-penetrating rigid body simulation[END_REF][START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF]]. Yet, in our case where we wish to simulate objects that are not overly constrained (avoiding, e.g., the scenario of a hand strongly compressing a hair wisp), our problem is likely to possess a solution most of the time [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF]. When designing a splitting solving method however, the non-existence issue may show up again for each subproblem. Thus, robust schemes, able to nd an approximate solution in the absence of a theoretical one, have to be designed.

• Finally, our scenarios may rapidly involve thousands or even tens thousands of contacts, making the problem very large in size. Solving algorithms should thus be able to scale up well with the size of the problem, that is, maintain good convergence properties at acceptable timings. Typically, in our scenarios where the time step might be as small as 1 ms, we aim at solving the DFCP in a few seconds up to a few minutes, in order to keep a reasonable (i.e., about one hour-long) frame rate.

To face up all these challenges at the same time, we have been opting for a practical study, where we have been rst testing a wide range of nonsmooth solvers, and then designing new ones targeted at the scenario to be simulated. The result of our study is not a mathematical analysis of a given algorithm with some proof of convergence 9 . Instead, for each scenario of interest (hair, cloth), we have been char- acterizing the main features of the resulting DFCP and we have built a dedicated method or heuristics, leading to some fast and robust solver.

Multiple equivalent formulations of Signorini-Coulomb

We have rst been extensively playing with the many dierent equivalent formulations of the Signorini-Coulomb law, each one of them leading to various solving

strategies. An extensive study of these dierent strategies is proposed in Gilles Daviet's PhD thesis [START_REF] Daviet | Modèles et algorithmes pour la simulation du contact frottant dans les matériaux complexes : application aux milieux breux et granulaires[END_REF]. Here we only mention the formulations and methods we have been frequently making use in our dierent applications, or comparing against.

7 Non-existence typically occurs when rank H < nd, in which case W is not necessarily strongly positive-denite [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF]. Moreover, if the DFCP admits a solution, it might not be unique.

In this chapter, where we wish to solve the direct problem (i.e., to retrieve generalized positions and velocities from the dynamic equations, a set of parameters, and initial conditions), we will not care too much about non-uniqueness of the solution. However, this will become an actual concern when attempting to solve the inverse problem and to identify parameters, as discussed in Chapter 4.

Linearized friction cone One popular approach, which is actually not equivalent to the Signorini-Coulomb law but stands for a linear approximation of it, consists in approaching the second-order cone K µ with a convex polyhedron (a faceted cone).

This approach, initially proposed by [START_REF] Klarbring | Contact problems with friction by linear complementarity[END_REF], has been widely used in Mechanics, Robotics and Computer Graphics [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Miller | Implementation of multi-rigid-body dynamics within a robotic grasping simulator[END_REF][START_REF] Erleben | Velocity-based shock propagation for multibody dynamics animation[END_REF][START_REF] Kaufman | Staggered projections for frictional contact in multibody systems[END_REF][START_REF] Otaduy | Implicit contact handling for deformable objects[END_REF], as it advantageously yields a Linear Complementarity Problem 10 (LCP), for which many solvers have been proposed in the literature. While this approximate formulation captures both dynamic and static friction regimes properly, it introduces some anisotropy in the sliding direction (see Figure 3.6) unless a high number of facets is used [Acary and Brogliato, 2008, Sec. 13.3.7]. A trade-o thus has to be made between accuracy and computational cost. In our applications we have not found this approach to be more ecient compared to an exact treatment of the friction cone, even with a low number of facets, and thus have rapidly discarded it. All the formulations presented in the sequel consider the exact Coulomb cone K µ . Disjunctive formulation We have already mentioned that solving the disjunctive Signorini-Coulomb formulation by enumerating each case is intractable for an arbitrary number of contacts n. However, when n = 1 (and we shall see the interest of considering the one-contact problem when dealing with splitting algorithms), then an exact algorithm can be designed [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. Actually, the only dicult case is the sliding case, which requires computing the roots of a quartic polynomial. This analytic one-contact solver has the advantage of providing a certicate of (non)-existence of solutions, and when a solution exists, to give its value down to the machine precision.

10 A Linear Complementarity Problem (LCP) reads 0 ≤ x ⊥ y ≥ 0 with y = Ax + b.

Nonsmooth contact dynamics for thin elastic objects 43

Complementarity formulation The Signorini-Coulomb law implies that r i lies inside the cone K µ i and u i in the half-space R + × R d-1 . Actually, this can be expressed as a cone complementarity problem by mapping this half-space to the dual cone of K µ i , which is K 1 µ i

(see Figure 3.7, (a) and (b)). Such a change of variables has been described by De Saxcé and Feng [1998],

u i := u i + µ i u i T e, (3.8)
and leads to the following expression of the Coulomb law,

r i , u i ∈ C µ i ⇐⇒ K 1 µ i u i ⊥ r i ∈ K µ i .
(3.9)

Note that for frictionless contacts (µ i = 0), the complementarity problem (3.9)

simply boils down to an LCP on the normal parts of u i and r i .

(r i , u i ) ∈ C µ i ⇐⇒ f (r i , u i ) = 0, (3.10) where f is a nonsmooth merit function from R d × R d to R d .
One classical example of such a nonsmooth function was notably provided by [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], and dened as

f AC : R d × R d -→ R d (r i , u i ) -→ Π R + r i N -ξ N u i N -r i N Π B d-1 (µr N ) r i T -ξ T u i T -r i T , 44
f DS : R d × R d -→ R d (r i , u i ) -→ Π Kµ i r i -ξ u i -r i ,
where ξ is a positive real number.

Functional formulations may be solved for by using a generalized Newton method, in which the usual gradient is replaced with a generalized gradient at non-dierentiable points [START_REF] Qi | A nonsmooth version of Newton's method[END_REF].

In our work [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF], summarized in Section 3.2.4, we have been deriving yet another merit function, using the equivalence between a Second-Order

Cone Complementarity Problem (SOCCP) and a modied Fischer-Burmeister 11 functional formulation [START_REF] Fukushima | Smoothing functions for secondorder-cone complementarity problems[END_REF]. Due to better smoothness properties compared to the two other formulations above, this modied Fischer-Burmeister formulation proves in practice easier to solve using a generalized Newton solver, compared to the two formulations above.

However, in all cases, when attempting to solve our large DFCP using any of the functional method in a global way (that is, solving all the contacts simultaneously), our solver always failed to scale up to more than one hundred contacting bers. More precisely we have found on our range of problems that the convergence success of the generalized Newton method was directly related to the conditioning number ν = nd m [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in ber assemblies[END_REF]. For ν < 1 (loosely connected bers), the method would converge properly. However, for ν ≥ 1 (tight ber packing, i.e., our typical scenario of interest), convergence would be very slow and eventually the method would fail to converge. For this reason we have thus changed strategies and instead used a splitting method, where contacts are solved one by one in an iterative fashion (see 3.2.4).

Still, one advantage of functional formulations is that they give a natural stopping criterion for any iterative method attempting to solve the reduced DFCP (3.7).

Indeed, the norm

1 2 Φ(r) 2 where Φ(r) is the R d n -value function such that Φ i (r) = f (r i , (Wr) i + b i ),
gives the global error of the iterative method. In our splitting algorithm, we have evaluated the error provided by the global Alart-Curnier function for stopping iterations of our algorithm (see Section 3.2.4).

Fixed-point sequence of convex optimization problems Another class of methods attempts to formulate the DFCP as a sequence of convex optimization 11 For complementarity problems in R, a well-known merit function is the so-called Fischer-Burmeister function dened as fFB(x, y) = x 2 + y 2 -x -y. One may easily check that 0 ≤ x ⊥ y ≥ 0 ⇐⇒ fFB(x, y) = 0. [START_REF] Fukushima | Smoothing functions for secondorder-cone complementarity problems[END_REF] has extended this formulation to second-order conical complementarity constraints.

problems [START_REF] Haslinger | Approximation of the Signorini problem with friction, obeying the coulomb law[END_REF][START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF]. In particular, Florent Cadoux has shown in his PhD thesis [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF][START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF] that solving the DFCP is equivalent to performing a xed point loop on either a primal (i.e., velocity-based) or dual (i.e., impulse-based) convex optimization problem.

Let us give some insights to this interesting method. Starting from De Saxé and Feng's change of variable (3.8) and concatenating it for every contact, one gets

u := u + E s ∈ L * , where s = [ u T 1 , . . . , u T n ] = [ u T 1 , . . . , u T n ] , E = blockdiag(µ i e i ), and L * = i K * µ i . The DFCP (3.5) then reads        Mv -f = H r (a) u = Hv + u f + Es (b) L * u ⊥ r ∈ L (c) s = [ u T 1 , . . . , u T n ] . (d)
The key of the Cadoux approach is to note that if s is xed, then (a), (b), and (c) are exactly the optimality conditions of a convex optimization problem subject to conical constraints, which can be equivalently formulated in a primal form,

min 1 2 v M v -f v (quadratic, strict. convex) Hv + u f + E s ∈ L * (conical contraints) (3.11)
or in a dual form,

       min 1 2 r W r + b r (quadratic, convex) r ∈ L (conical constraints) W = H M -1 H b = H M -1 f + u f + E s (3.12)
Finally, the full DFCP can be solved thanks to a global xed-point loop,

F (s) = s with F i (s) := u i T (s) , (3.13)
where at each iteration, u i (s) is computed by solving either the primal problem (3.11),

or the dual problem (3.12).

From this algorithm, it becomes now clear that the DFCP cannot be simply cast into a convex optimization problem. Interestingly, previous approaches which have approximated the DFCP with a convex optimization problem [START_REF] Anitescu | Optimization-based simulation of nonsmooth rigid multibody dynamics[END_REF][START_REF] Mazhar | Using Nesterov's method to accelerate multibody dynamics with friction and contact[END_REF] can be interpreted as a single xed-point iteration of the procedure above. One shortcoming inherent to this approximation is that the relative velocity u, which should be tangent to the contacting surface when contact is active, is biased by the adding of a non-zero normal component, which in practice may cause articial jumps of sliding objects during simulation.

We have been testing the two variants of Cadoux's method on our large ber problems, using dierent solvers. For instance, we have experimented solving the 46 Chapter 3. Frictional contact models for large multibody systems dual problem with either a projected gradient method, or interior points. Solving the primal problem with interior points actually gave better results of convergence, due to the strongly convex objective in this case. However, although the latter variant turned out to perform better compared to global functional formulations, it did not allow us to handle more than a few hundreds contacting rods: beyond, the solver would fail to converge in a reasonable amount of time.

Splitting (Gauss-Seidel) formulation Instead of solving the DFCP globally, one may instead solve for each contact problem sequentially, and then iterate on the full set of contacts until convergence is reached. Such a sequential splitting approach, often called Gauss-Seidel due to its resemblance with the Gauss-Seidel algorithm to solve linear systems, is actually classically used for simulating large granular assemblies [START_REF] Dubois | The non smooth contact dynamic method: recent LMGC90 software developments and application[END_REF], due to its scalability properties. Convergence

to the solution at a high precision is however known to be slow in practice [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]. This was however not a problem in our case, where we only request to solver to converge to an acceptable precision, without the need to get beyond visual perception.

When adopting the Gauss-Seidel strategy, one is then left with the choice for the local solver. The key insight we have learnt from our multiple tests is that failing to solve one single contact properly may often cause the global loop to fail and introduce instabilities in the ber dynamics, ultimately leading to a simulation crash. We have thus strived to build a local solver as robust as possible, by devising a hybrid strategy.

A local hybrid solving strategy

Large ber assemblies Facing robustness and scalability issues when using stateof-the-arts approaches, we have designed a new frictional contact solver by coupling an iterative Gauss-Seidel strategy together with an extremely robust one-contact solver [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. The local solver combines two dierent methods for solving the 1-contact problem in contact force variables: on the one hand, it primarily relies upon the root nding of the modied Fischer-Burmeister function [START_REF] Fukushima | Smoothing functions for secondorder-cone complementarity problems[END_REF], thanks to a nonsmooth Newton method; on the other hand, in the rare case when the previous method fails to converge, the solver automatically switches to the more expensive but analytical solver which is garanteed to nd the solution to the 1-contact problem when it exists [START_REF] Bonnefon | Quartic formulation of Coulomb 3D frictional contact[END_REF]. The latter solver, called in 1% of the cases, thus acts as a failsafe of the method.

Our global solver proved to converge well in scenarios involving thousands bers subject to tens thousands frictional contact points, and thus allowed us to enhance considerably the realism of hair simulations. Since then, our method has been leveraged by academia and the special eects industry for simulating hair and fur realistically [START_REF] Kaufman | Adaptive nonlinearity for collisions in complex rod assemblies[END_REF]. In particular, our method has been used in production by Weta Digital, the leading company in special eects, and our accompanying software has been licensed to L'Oréal Research and AGT Digital. We have also recently adapted our solver for coping with frictional contact in cloth modeled as a nodal system [START_REF] Daviet | Fast Cloth Simulation with Implicit Contact and Exact Coulomb Friction[END_REF], as explained below. Cloth: towards a dedicated nodal solver In the ber assembly case, the matrix M is block-diagonal, so that the Delassus operator can be computed in an ecient way by leveraging sparse-block computations [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. This justies solving the reduced DFCP (3.6), where primary unknowns are forces r. For cloth however, where primal variables (nodal velocities of the cloth mesh) are all interconnected via elasticity through implicit forces, the method developed above is computationally inecient. Indeed, the matrix M (only block-sparse, but not blockdiagonal) is costly to invert for large systems and its inverse is dense. However, we may leverage the fact that generalized velocities of the system are 3D velocities, and slightly adapt the previous solver so as to recover eciency.

We assume body-cloth contacts occur at cloth vertices only, i.e., each contact involves only one vertex. A key observation is that each nonzero block of J related to a vertex i is simply a rotation matrix E i , corresponding to the local contact basis for the vertex i. Our idea is then to build the square block-diagonal matrix G with G i,i equal to E i if i is in contact, and to the 3 × 3 identity matrix otherwise. G can thus be inverted trivially as G -1 = G T . Augmenting u and r so that their size matches that of v, we can write (3.5) as

   GMG T u -Gf = r v = G T u u i , r i ∈ C µ i if i in contact and r i = 0 otherwise. (3.14)
Denoting W = GMG T and b = Gf , we obtain a system that is very close to the reduced DFCP (3.6), except that u and r have reversed roles. To retrieve symmetry, we apply De Saxcé's change of variable (3.8) u i := u i +s i (u) so that u i is orthogonal to r i . For a xed value of s, we identify (3.14) as the KKT conditions of the convex in size which typically distinguishes them from powders, made of much smaller grains. Considering this limit size, grain-grain interactions in granulars are mainly dictated by contact and dry friction, while air-grain interactions can be neglected. Cohesion between grains may furthermore be considered, typically in the case of wet materials.

Due to their discrete nature at macroscopic scale, granular materials are highly dissipative systems which exhibit various physical states, ranging from liquid (when the material continuously ows) to solid (when the material rigidies due to internal friction and undergoes plastic deformations), or even gaseous (when grains are dispersed in air and interact mainly through impacts). Such a rich visual behavior, very distinctive from that of Newtonian uids (e.g., water), explains why simulating granular matter has, beyond classical applications in Physics and Mechanical Engineering, gained increased attention in Computer Graphics for a few decades [START_REF] Zhu | Animating sand as a uid[END_REF][START_REF] Alduán | Simulation of high-resolution granular media[END_REF][START_REF] Narain | Free-owing granular materials with two-way solid coupling[END_REF].

However, unlike Newtonian uids or elastic materials for which robust physical models have been developed in the last centuries, granular matter still resists some comprehensive physical understanding: no universal equation exists yet for describing the collective behavior of such a peculiar material.

As a result, various methods have been developed for studying and simulating granular materials, following two main strategies. On the one hand, many approaches naturally rely upon a discrete model, where grains are represented individually as rigid bodies and grain-grain interactions are accounted for using a frictional contact model [START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF]. Because they explicitly model the material at the grain scale, those methods are able to capture the dierent states of the granular matter and their transitions quite accurately. However, considering that a simple teaspoon of sugar already consists of around 10 5 grains, one may easily imagine that such approaches rapidly suer from scalability issues.

On the other hand, continuum-based methods, where granular matter is modeled as a viscoplastic uid, have started to be investigated for a few decades. Though less accurate than discrete models, such methods have nevertheless demonstrated some good qualitative agreement with real scenarios in plastic and dense owing regimes.

In particular, a popular constitutive law is the µ(I) rheology [START_REF] Jop | A constitutive law for dense granular ows[END_REF], relying upon a pressure-dependent yielding model referred to as the Drucker-Prager yield criterion (detailed below) with a non-constant friction coecient.

For the sake of tractability, existing numerical models however make some important simplications, for instance assuming an incompressible ow, a uniform density, or a regularized or linearized Drucker-Prager yield criterion. Such approximations often result in a considerable loss of realism, for instance preventing a heap of sand from stabilizing after owing, or introducing some articial anisotropy in the ow.

Our main objective towards this work was to treat the Drucker-Prager rheology numerically with as little as possible downgrade. Typically, we addressed the nonsmooth rheology directly with no regularization and no linearization, and relaxed the incompressibility condition to incorporate a more realistic, unilateral compressibility condition. Altogether, these new settings allowed us to improve substantially 50 Chapter 3. Frictional contact models for large multibody systems the quality of simulations of large granular ows, while keeping reasonable computational timings.

Caveat Our work on the continuum modeling of granulars with a nonsmooth rheology is recent, and was the core contribution of Gilles Daviet's PhD thesis [START_REF] Daviet | Modèles et algorithmes pour la simulation du contact frottant dans les matériaux complexes : application aux milieux breux et granulaires[END_REF]. Here we only sketch the main ideas of our method, without going deeply in technical details. The interested reader may nd all technical issues in our publications [Daviet and Bertails-Descoubes, 2016b,a].

Notation We consider an Eulerian uid representing a granular medium. We use standard tensor notation for manipulating scalar and vector eld associated to this . The gradient of a scalar eld φ and of a vector eld v are written as ∇φ (vector) and ∇v (matrix), respectively, and the divergence of a vector v eld as ∇•v (scalar). Given a vector eld v, let D(v) := 1 2 (∇v +(∇v) T ) be the symmetric part of its gradient and W(v) := 1 2 (∇v -(∇v) T ) its skew-symmetric part. Finally, u denotes the velocity eld of the grains, and we use the shorter notation ε := D(u) for the so-called strain rate tensor. Note that Tr ε = ∇ • u.

Macroscopic yield criteria

Continuous models for granular materials are derived upon macroscopic yield criteria, i.e., inequalities that should be satised by the principal stresses the eigenvalues of the stress tensor σ for the material to remain stable. The most well-known is the Mohr-Coulomb criterion. For cohesionless 3D materials such as dry sand, and with σ 3 ≤ σ 2 ≤ σ 1 the principal stresses, it states

σ 1 -σ 3 ≤ sin ϕ (σ 1 + σ 3 ) ,
where ϕ is the so-called friction angle, which corresponds to the inclination of a stable heap of the granular material under gravity.

In the space of principal stresses, the Mohr-Coulomb criterion denes a cone with hexagonal basis, which is numerically unwieldy. It is thus often discarded in favor of the Drucker-Prager yield criterion on the second invariant J 2 = 1 2 Tr(Dev σ) 2 of the deviatoric stress tensor, which in 3D reads

J 2 ≤ -μ σ 1 + σ 2 + σ 3 3 where J 2 = 1 6 i =j (σ i -σ j ) 2 . (3.15)
This criterion denes a second-order cone in the space of principal stresses, and μ is called the friction coecient. Note that in 2D, the Mohr-Coulomb and Drucker-Prager criteria are equivalent.

3.3. Continuum modeling of granular materials 51

Nonsmooth rheology DP(µ)

The particularity of our uid is to be non-Newtonian: our rheology relies upon a viscoplastic constitutive law combining the Drucker-Prager yielding criterion together with an unilateral incompressibility constraint.

As done usually, we decompose the solid phase stress tensor as σ g := η ε + σ c , Moreover, the maximum dissipation principle states that in the yielded regime, friction should be saturated and the frictional stress tensor should be colinear to the deviatoric part of the strain rate. The deviatoric part of σ c should thus satisfy one of the two regimes,

     Dev σ c = (μp) Dev ε | Dev ε| if Dev ε = 0 (yielded) | Dev σ c | ≤ μp if Dev ε = 0 (unyielded).
(3.16)

Unilateral compressibility constraint Most continuum-based models for granulars consider the uid to be perfectly incompressible. In contrast, we want to take into account the typically asymmetric yielding behavior of granulars by allowing the uid to expand as much as desired, while strictly preventing compaction.

With a dense ow hypothesis, the unilateral compressibility constraint can be expressed simply as ∇ • u ≥ 0. We set the pressure p to enforce this inequality, i.e.,

p ≥ 0 if ∇ • u = 0 p = 0 if ∇ • u > 0
or, using an equivalent complementarity notation,

0 ≤ p ⊥ ∇ • u ≥ 0.
(3.17)

In our results, we have shown that relaxing the common incompressibility assumption ∇ • u = 0 prevents the arising of an ill-dened rheology in some typical scenarios such as the ow in the wake of an obstacle. For instance, in [START_REF] Chauchat | A three-dimensional numerical model for dense granular ows based on the mu(i) rheology[END_REF], where incompressibility is enforced, a negative pressure is observed behind the obstacle, which is outside the domain of validity of the well-accepted µ(I) rheology [START_REF] Jop | A constitutive law for dense granular ows[END_REF]. In contrast, our method guarantees that the pressure remains non-negative (see Figure 3.10). Moreover our complementarity constraint (3.17) naturally ts in with our numerical framework, without adding any computational cumbersomeness.

52 Chapter 3. Frictional contact models for large multibody systems Solution set DP(µ) Our full rheology (3.16 3.17) can be compactly rewritten as

               Dev σ c = - μ d Tr(σ c ) Dev( ε) | Dev( ε)| if Dev ε = 0 | Dev(σ c )| ≤ μ d Tr(σ c ) if Dev ε = 0 0 ≤ 1 d Tr(σ c ) ⊥ Tr( ε) ≥ 0, (3.18)
which is reminiscent of the Signorini-Coulomb law derived in Section 3.2.1, with ε playing the role of the relative velocity u, and σ c playing the role of the contact force r. More precisely, by leveraging the following isomorphism converting the vector space to the space of symmetric tensors S d ,

χ : R × R T → S d (a; b, c) → b c c -b + a I if d = 2 (a; b, c, d, e, f ) →    b -c √ 3 d e d -b -c √ 3 f e f 2c √ 3    + √ 2 √ 3 a I if d = 3,
we can prove that our rheology is similar to the Signorini-Coulomb law, albeit expressed in higher dimension (in dimension 3 for the 2D case, and in dimension 6 for the 3D case) [START_REF] Daviet | Nonsmooth simulation of dense granular ows with pressure-dependent yield stress[END_REF].

Discretization After discretizing the full dynamic equations of our uid using nite elements, together with an implicit discretization of the rheology (3.18), we build a discrete system which shares the exact same structure as our DFCP introduced in Section 3.2.2, though with higher dimension. We can thus leverage exactly the same tools as before: in particular, we use our ecient Gauss-Seidel solver featuring a modied Fischer-Burmeister (extended to higher dimensions) zero-nding solver as the local solver. Note however that the original enumerative solver was specically designed for 3D, and cannot be trivially extended to higher dimensions.

Some results

We consider a gravity-induced 2D granular medium inside a narrow channel (of width L), owing around a cylinder of diameter D = L 4 , with no-slip boundary conditions on the sides of the channel and on the cylinder. • First, it validates the benet of allowing dilation of the ow, as the pressure in the wake (right part) of the obstacle is indeed zero, and does not become strictly negative as in [START_REF] Chauchat | A three-dimensional numerical model for dense granular ows based on the mu(i) rheology[END_REF];

• The zone of highest pressure is not located at the very front of the obstacle.

Instead, we observe the formation of a high-pressure arch above this point;

• Above this arch, there exists a region where the pressure is lower than further upstream; moreover the transition between these two zones is very abrupt, as shown by the absence of two isosurfaces on the gure. Note that this 54 Chapter 3. Frictional contact models for large multibody systems phenomenon disappears when considering free-slip boundary conditions for the channel walls.

Extension to a varying volume fraction

Our previous model does not consider a volume fraction eld. Instead, the uid is considered to be dense that is, already at its maximal volume fraction φ max everywhere inside the simulation domain. The material is theoretically allowed to dilate, but not to compress, through the complementarity condition 0 ≤ ∇•u ⊥ p ≥ 0.

However, in this model, the use of a xed domain precludes an actual dilatation of the material through time. More fundamentally, even if the volume fraction eld were somehow tracked through time (e.g., with particles), the ∇ • u ≥ 0 condition would still prevent the material from recompacting after dilatating. We thus follow Narain and colleagues [2010] to account for a non-constant volume fraction φ, and instead enforce the complementarity condition

0 ≤ φ max -φ ⊥ p ≥ 0. (3.19)
where φ max ≤ 1 is the maximum packing fraction for the material, beyond which grains cannot be compacted anymore. For monodisperse spherical grains, φ max usually lies in the 0.550.65 range. We used φ max = 0.6 in all our simulations. Considering this non-constant volume fraction φ, we rederived in [Daviet and Bertails-Descoubes, 2016a] the dynamics of our uid together with the new Drucker-Prager condition. Using a semi-implicit time integration scheme together with a Material Point Method for spatial discretization, our nal problem boiled down, again, to a discrete problem with a similar structure compared to our initial DFCP.

This nice property allowed us to simulate free-owing sand in several minutes per frame only (see Figure 3.11).

Conclusion: Towards a continuum-based model for brous materials

In this chapter we have motivated our choice for nonsmooth frictional contact laws, in order to capture the typical emerging eects due to contact and friction between solids. We have also presented a range of dierent formulations of the Signorini- medium such as hair. One key ingredient that we have not been considering in our previous models is the inuence of air inside divided materials. Typically, air plays a considerable role in hair motion. To advance in that direction, we have started to look at a diphasic uid representation of granular matter, where a Newtonian uid and the solid phase are fully coupled, while the nonsmooth Drucker-Prager rheology for the solid phase is enforced implicitly [START_REF] Daviet | Simulation of DruckerPrager granular ows inside Newtonian uids[END_REF]. This rst approach could be a starting point for modeling immersed granulars in a liquid, or ash clouds, for instance.

There still remains a long way to go for taking into account long bers instead of isotropic grains in the solid phase. Coupling ber elasticity with our current formulation is clearly an open and challenging problem that we would like to deal with in the future.

Chapter 4

From geometry to mechanics: Inverse design of bers In a third part of our work, we have considered some new challenges regarding inverse static design. Indeed, while current simulators may succeed in reaching a good level of realism, they remain dicult to control in order to achieve a precise artistic goal or, more generally, to match a target observation. More precisely, to generate some desirable shapes and motions, one should be able to feed a simulator with the right parameters. Finding such parameters remains a very dicult task, which is often performed through a tedious trial and error process. To make this task fully automatic, we have started looking at inverse solutions in the case where a static shape is provided as input: the inverse model should be able to interpret automatically this shape as a stable equilibrium 1 of the simulator, under gravity and other external forces such as contact and friction.

58 Chapter 4. From geometry to mechanics: Inverse design of bers Summary of contributions In the case of an isolated ber, we have shown that inverting any of our super-model [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF][START_REF] Casati | Super space clothoids[END_REF]] boils down to two decoupled problems that are both easy to solve [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF], 2011;[START_REF] Bertails-Descoubes | Super-clothoids[END_REF]: rst, an equilibrium condition which appears to be linear in the intrinsic shape of the ber, thanks to the curvature-based parameterization of our ber models; second, a sucient stability condition that can be simply set by xing a lower-bound for the ratio of stiness over mass. Actually, the only remaining diculty is to solve a merely geometric tting problem converting a curve as a piecewise helix or clothoid.

In the case of helical tting, we have already brought some ecient and robust solutions [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF][START_REF] Derouet-Jourdan | Floating tangents for approximating spatial curves with G 1 piecewise helices[END_REF]].

In the presence of contact and friction, Coulomb sticking constraints have to be considered, which makes the overall inverse problem nonsmooth and ill-posed. We have shown that assuming known mass and stiness and a simplied inverse model, it is possible to recover a reasonable intrinsic shape as well as frictional contact forces at play, both of them satisfying exactly static equilibrium [Derouet-Jourdan et al., 2013a]. This work allowed us, for the rst time, to animate a few hair geometries stemming from recent hair captures, such as the one depicted in Figure 4.3(a).

Finally, we are currently looking at the inverse problem in the continuous case, that is on the strong form of the Kirchho static equations [START_REF] Bertails-Descoubes | About curves, frames, and material rods[END_REF]. This ongoing work allows us to better characterize the space of solutions of the inverse problem in the case of an isolated ber, through a generic curve-angle parameterization of the rod. This also allows us to understand how discrete models may automatically select a subspace of solutions, by imposing a certain form for the material frame of the rod. Such insights give us the hope to reach in the near future some actual identication process, i.e., a clear selection of a solution among the set of potential solutions.

Organization of the chapter Section 4.1 species our motivation for tackling the inverse design problem, and discards the generic nonlinear optimization approach in favor of a more eective approach dedicated to Kirchho rods. Sections 4.2 tackles the discrete inverse design problem relying upon super-models, in the case of an isolated ber under gravity only. Then Section 4.3 takes a look at the more complex case of a ber assembly subject to gravity and frictional contact. Finally, Section 4.4

elaborates on the problem of retrieving material curvatures from a merely geometric curve, rst in the continuous case, then in the discrete case. actually suggested a very interesting scientic problem to me: how to let the artist be completely free to design the hair shape she/he desires, and to have the hair simulator comply with this input shape and predict its motion?

Obviously, the artist would draw hair as perceived within a physical world, that is, under external forces such as gravity and contacts. Moreover, at the character design stage, it is very likely that the drawn hairstyles intend to correspond to a hair conguration at stable equilibrium. However, such a shape is not readily usable by a physics-based simulator, which in contrast requires, among its parameters, the intrinsic shape 2 of the object to be simulated, that is the shape it would have in the absence of external forces (see Figure 4.2).

More generally, the intended hair design may come from various processes, ranging from artistic hairstyling to automatic hair capture. In all cases, one is left with, on the one hand, an observed3 (or target) shape q obs which integrates the eect of all surrounding forces, and on the other hand, a simulator which has to be fed with the right parameters, in particular q, so as to match the input shape q obs at stable equilibrium.

To simulate an elastic object from such a target shape q obs , a naïve though common practice so far in Computer Graphics was to use q obs as the intrinsic shape of the object, that is, to set q = q obs . However, when launching the simulator, the object would sag under gravity, thus losing the initial shape and ruining the design 2 Recall that the intrinsic shape of a Kirchho rod was given in Chapter 2 by the intrinsic curvature function κ(s); in our nite-dimension Lagrangian settings, it corresponds to the vector q. The initial conguration q(t = 0) is set to a trivial intrinsic shape q here, a set of straight rods (a). Then, after feeding the simulator with other parameters (such as the mass and stiness of each rod), a static shape q eq can be computed under gravity and frictional contacts (b).

work made by the artist, or the 3D geometry patiently reconstructed from real hair capture (see Figure 4.3).

Inverse design problem To address this problem eectively, an inverse modeling process 4 thus has to be investigated. The inverse model should be able to interpret automatically the target shape q obs as a stable equilibrium of the simulator, under gravity and other external forces such as contact and friction.

As described in Section 2.2, the Kirchho rod model possesses two kinds of parameters. On the one hand, its intrinsic curvature κ, discretized as vector q, and which may vary spatially along the centerline. On the other hand, its material parameters, composed of the linear mass density ρS and the bending/twisting stinesses of the rod K 3 = diag 3 (K 0 , K 1 , K 2 ); material parameters are assumed to be constant with respect to space and time.

A generic and standard way to model our inverse problem is through least squares minimization [START_REF] Kern | Problèmes inverses[END_REF]. That is, given the input conguration q obs , one seeks to solve the following constrained minimization problem, min (q,q) F (q,q) ∈ K 1 2 qq obs 2 , (4.1)

4 In Mechanical Engineering, inverse problems are usually split in two categories [START_REF] Beck | Inverse problems and parameter estimation: integration of measurements and analysis[END_REF]]: (a) inverse design problems, where material parameters are known and one is looking for the intrinsic shape of the object, and (b) inverse measurement problems, where material properties are searched for. Our goal here is clearly to solve for an inverse design problem, but in the medium term, we would also like to extract some information about the material properties, from the input shape.

4.2. Case of an isolated ber 61 (a) Observed shape q obs (b) Sagged shape q eq under gravity where the constraint F (q, q) ∈ K expresses a generalized (stable) equilibrium condition, possibly incorporating frictional contact forces when the set K is not reduced to {0}. The above formulation is nonlinear, and even nonsmooth in the presence of frictional forces: it is thus challenging to solve. We have recently attemped to solve this problem in the case of nodal models for cloth using the adjoint method, and have reported our results in the research report [START_REF] Casati | Inverse elastic cloth design with contact and friction[END_REF]. In particular, we have shown how to extend the adjoint method classically used in the context of bilateral constraints to deal with frictional contact constraints (see details in the PhD thesis of R. Casati [START_REF] Casati | Quelques contributions à la modélisation numérique de structures élancées pour l'informatique graphique[END_REF]). Though promising, such a method is, even in the contactless case, prone to convergence issues inherent to nonlinear optimization, and much work still remains to be done to improve the robustness of our inverse process.

Inspecting Kirchho equations as well as our discrete curvature-based models (described in Chapter 2), we have fortunately noted that a much more eective inverse model could be built, which avoids nonlinear optimization. Result is not only the design of a simple, fast and robust inversion process, but also a better characterization of the whole space of solutions, which opens the way for a reliable parameter identication process in the short term.

Case of an isolated ber

Let us consider a discrete curvature-based model (super-model) clamped at one end, and subject to gravity only (see Figure 4.4(b)). We assume the conguration q 62 Chapter 4. From geometry to mechanics: Inverse design of bers to be given as input, and we are looking for the parameters q (intrinsic shape), ρS (linear mass density and K 3 = diag 3 (K 0 , K 1 , K 2 ) (bending and twisting stinesses) so that q corresponds to a stable equilibrium of the super-model. using our oating tangent algorithm (see Section 4.2.4), and nally automatically computes the parameters of a super-helix so that its centerline matches the input curve at stable equilibrium under gravity (c).

Necessary and sucient condition for equilibrium

Let us rst write down the equilibrium equation. From Equation (2.8), dropping all time-dependent terms and making parameter dependence explicit with brackets yields

K [K 3 ] (q -q) + G [ρS] (q) = 0, (4.2)
where

K [K 3 ] = K 3 ⊗ L N and G [ρS] (q) = -ρS L 0 ∂r ∂q g ds G(q)
, where L N is a symmetric positive-denite matrix of size N , which only depends on the length = { 0 , . . . , N -1 } of the rod elements. For a super-helix, L N is simply diagonal (see Expression (2.11)), whereas for a super-clothoid, L N is tridiagonal (see Expression (2.12)).

Note that if E p (q) is the potential energy of the rod, Equation (4.2) is equivalent to writing ∇E p [q,K 3 ,ρS] (q) = 0, meaning that we are looking for the parameters q, K 3 and ρS so that q is a critical point of energy of the rod.

Assuming the material parameters ρS and K 3 to be constant, Equilibrium (4.2) 4.2. Case of an isolated ber 63 is exactly satised for a unique family of solutions,

q = q + K -1 [K 3 ] G [ρS] (q) = q -ρS diag 3 K -1 0 , K -1 1 , K -1 2 ⊗ L -1 N G(q) (4.3)
using inverse of Kronecker product [Prasolov, 1994, Section 27.4],

where K 0 , K 1 , K 2 and ρS can span the entire positive space.

In practice, solving (4.3) amounts to choosing a positive set of parameters {K 0 , K 1 , K 2 , ρS}, and then computing q by solving a full rank linear system of ∼ 3 N equations, which turns out to be diagonal in the case of super-helices, and tridiagonal in the case of super-clothoids. The number of elements N being generally low (≤ 20), the latter problem is thus extremely fast to solve.

Sucient condition for stability

Evaluating the stability of the equilibrium requires the computation of the Hessian matrix ∇ 2 E p of potential energy. The equilibrium will be stable if ∇ 2 E p is a positivedenite matrix, i.e., if all its eigenvalues are positive.

Dierentiating the left-hand member of Equation (4.2) gives the expression for the Hessian of potential energy,

∇ 2 E p [K 3 ,ρS] = K 3 ⊗ L N + ρS S(q), (4.4)
where S is a real symmetric matrix of size 3 N which (nonlinearly) depends on the input conguration q, but remains independent of our unknown parameters. Note that ∇ 2 E p (q) thus only depends on the material parameters K 3 and ρS, in a linearly way, while dependence w.r.t. the intrinsic shape q has vanished. Using Expression (4.4), we can nd a sucient condition for the Hessian ∇ 2 E p to be positive-denite. Indeed, Horn's theorem [START_REF] Fulton | Eigenvalues, invariant factors, highest weights, and Schubert calculus[END_REF] implies that the minimum eigen value of a sum of real symmetric matrices is greater than or equal to the sum of the eigen values that are minimum for each matrix. Let λ m , µ m (q), τ m be the minimum eigen values of ∇ 2 E p , S, and L N , respectively. Let K = min{K 0 , K 1 , K 2 }.

From [Prasolov, 1994, Section 27.4], we deduce that the minimum eigen value of the Kronecker product K 3 ⊗L N is the product of minimum eigen values K τ m . Applying Horn's theorem then gives

λ m ≥ K τ m + ρS µ m (q) (4.5)
with K > 0 and τ m > 0. Inequality (4.5) provides a lower bound for the smallest eigenvalue of ∇ 2 E p . It thus yields a sucient condition on K 3 and ρS for guaranteeing that all the eigenvalues of ∇ 2 E p are positive,

K ρS > - µ m (q) τ m . (4.6)
64 Chapter 4. From geometry to mechanics: Inverse design of bers

Recall that the two eigen values µ m (q) and τ m can be directly computed from the input conguration q and the element lengths , using standard eigen value extraction algorithms 5 .

According to Inequality (4.6), by either increasing K or decreasing ρS, the stability of the rest shape under gravity can thus always be guaranteed, whatever the input conguration q. This is actually in line with intuition. Imagine an input ber oriented downwards. To be stable, this conguration will not require a high stiness since it is already pointing towards a direction that is preferred by gravity. Now, if the ber is oriented upwards, a large stiness (or equivalently, a small mass) will be required to guarantee that the ber can stably holds in this position and will not bend downwards. This is quite similar to applying gel on hair bers so that they comply to some arbitrary desired shapes, even if those are in strong contradiction with gravity forces.

Inverse design algorithm

Computing the parameters of the super-model so as to match an arbitrary conguration q at stable equilibrium is thus elementary in the case where only gravity is involved. The user may rst set the stiness K and the linear mass density ρS so that their ratio lies in the authorized halfspace given by (4.6), then compute the intrinsic curvature q satisfying equilibrium by solving the system of equations (4.3).

If unhappy with the motion of the ber which looks to sti or too soft, she can modify the stiness, check that stability is still guaranteed (if not, automatically retrieve the minimum stiness satisfying stability), and nally recompute the intrinsic curvature yielding equilibrium. This simple design process is summarized in Note that our stability criterion (4.6) only provides a sucient condition for sta-5 In the case of a super-helix, the minimum eigen value τ m boils down to the minimum length of elements, mini{ i}.

Case of an isolated ber 65

bility, meaning that equilibrium may still be stable when decreasing stiness beyond the suggested lower bound. In practice however, we have noticed that violating this condition always causes equilibrium to become unstable. This gives us the hope to establish in the short term a necessary condition for stability, useful this time not just for designing stable rods matching an input curve, but for characterizing precisely the range of potential stiness and mass parameters from a real observed rod under gravity (see Conclusion of this chapter).

From an arbitrary curve to a discrete material conguration

So far we have considered that a discrete material conguration q was provided as input. However, this is not exact. The true input data is an arbitrary curve Γ(s), represented with a nite number of degrees of freedom, for instance as a spline (see Figure 4.4(a)). Two steps then need to be completed for reconstructing an admissible material conguration q, before performing inverse statics.

The rst step is to convert the arbitrary input curve into a kinematically admissible centerline r(s) for a super-model, that is, into a C 1 -smooth piecewise helical curve if dealing with the super-helix model, or a C 1 -smooth piecewise clothoidal 6 curve if dealing with the super-clothoid model.

Once a piecewise arc curve r(s) has been reconstructed, a second step consists in reconstructing admissible material curvatures and twist q for the discrete rod, which is equivalent to reconstructing a kinematically admissible material frame R(s) for the rod. To anticipate on the following, we have shown that the only admissible material frame for a super-helix consists in building a continuous frame which shares, on each element, the same Darboux vector as the Frenet frame. We defer the corresponding proof to Section 4.4 where the general question of retrieving a material frame from a curve will be deeply examined, from the continuous (Kirchho ) case to the discrete case. Here, we focus only on the rst step, that is on the geometric problem aiming at converting an arbitrary curve into a piecewise arc curve.

Piecewise helical tting In our work we have tackled the rst step mentioned above in 3D in the case where the target super-model is a super-helix 7 . That is, we have addressed the geometric problem of converting an arbitrary input curve into a C 1 -smooth piecewise helix.

6 In 3D, the notion of a clothoidal curve is not clear given our denition of a space clothoid element in Chapter 2. Indeed, unlike what happens in the piecewise uniform curvature case (superhelix), assuming that the material curvatures and twist are piecewise linear does not imply that the geometric (or Frenet) curvature and torsion of the centerline are piecewise linear (the reader may take a look at Equations (4.16) to get convinced about it). As noted already in Chapter 2, the shape of the centerline of the super-clothoid model is thus more general than the so-called 3D

Euler spiral [START_REF] Harary | 3D Euler spirals for 3D shape completion[END_REF].

7 In 2D, approximation algorithms are much simpler to design due to the absence of torsion. Our 3D oating tangent algorithm directly translates to the 2D problem consisting in tting a curve to a smooth piecewise circular curve [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF]. Moreover, ecient methods for tting a 2D curve to a smooth piecewise clothoid curve have been designed, for instance [START_REF] Mccrae | Sketching piecewise clothoid curves[END_REF] that we leverage in [START_REF] Bertails-Descoubes | Super-clothoids[END_REF]. 66 Chapter 4. From geometry to mechanics: Inverse design of bers In a nutshell, our algorithm consists in segmenting the input rod into N pieces of equal length, and approximating each curve segment with a helix under the constraint that two neighboring helical segments should have a C 1 -smooth connection, i.e., tangents should remain continuous at the junction. To achieve this, we enforce tangents at junction points to match initial tangents, but since having both tangents and positions xed may not be compatible with an interpolating helix, we relax the initial positions of the junction points, and minimize the error between admissible and initial positions. The key of our approach relies upon the following theorem [START_REF] Ghosh | Geometric approximation of curves and singularities of secant maps. A dierential geometric approach[END_REF][START_REF] Derouet-Jourdan | Floating tangents for approximating spatial curves with G 1 piecewise helices[END_REF],

Theorem 4.1. Given two points p 0 and p 1 such that p 0 = p 1 and two tangents t 0 and t 1 such that t 1 = ±t 0 , there exists a unique short helix starting at p 0 with tangent t 0 and ending at p 1 with tangent t 1 if and only if p 1 -p 0 , t 1 -t 0 = 0, (4.7) which allows us to formulate admissible positions recursively as a simple linear function of the input tangents. Our algorithm, coined oating tangents algorithm, then simply amounts to solving a full-rank linear system of size N ; this obviously proves order of magnitude faster compared to a nonlinear least-squares optimization method (see Figure 4.7) Theorem 4.1, illustrated in Figure 4.6, was rst stated in Ghosh's PhD thesis [START_REF] Ghosh | Geometric approximation of curves and singularities of secant maps. A dierential geometric approach[END_REF], but proof was incomplete. We have completed the missing parts in [START_REF] Derouet-Jourdan | Floating tangents for approximating spatial curves with G 1 piecewise helices[END_REF], especially thanks to the PhD work of Alexandre Derouet-Jourdan [START_REF] Derouet-Jourdan | Inversion statique de bres: de la géométrie de courbes 3D à l'équilibre d'une assemblée de tiges mécaniques en contact frottant[END_REF]. Note that the proof of existence is constructive, and provides an ecient algorithm for building the unique short helix once Condition (4.7) is satised.

Case of a ber assembly subject to frictional contact

While inverse design of an isolated ber is interesting from a theoretical point of view, and may even apply to a few examples such as the tail of an animal or a plant stem, in practice most useful scenarios involve interactions with other objects, and typically contact and friction forces.

We have thus started to investigate the inverse design problem for a ber subject not only to gravity but also to frictional contact forces. More specically, we have considered a ber assembly subject to external contacts (with a body) as well as ber-ber contacts, in the presence of Coulomb friction. A typical example of such a system is a head of hair, usually modelled in Computer Graphics with a few hundreds to thousands thin elastic rods. With the recent advance in 3D hair reconstruction from real hair photographs, many geometric datasets, composed of 3D curves which resemble real hair shapes, are nowadays available (see Figure 4.8(a)). 4.3. Case of a ber assembly subject to frictional contact 67

x y z t 0 t 1 t(ϕ) ϕ ϕ * p 0 p 1 α(ϕ) α ⊥ (ϕ) Figure 4.6
In red, the unique short helix that interpolates p 0 and p 1 with tangents t 0 and t 1 when Equation (4.7) is satised. The helix is said to be short if ϕ * ≤ 2π.

Underdetermined conditions for equilibrium

Let us consider a dynamic system composed of a set of N F bers discretized as super-models. Let q be the vector containing the m generalized coordinates of the full system. We assume the system to be subject to n frictional contacts, which may occur between a ber and an external rigid object, or between two dierent bers, or within a single ber.

Using notations of Chapter 3, our previous equilibrium equation (4.2) gets transformed into the constrained problem

K [K 1 3 ,K 2 3 ,...,K N F 3 ] (q -q) + G [ρS 1 ,ρS 2 ,...,ρS N F ] (q) = H(q) r ∀i = 1 . . . n, r i ∈ K µ i , (4.8)
meaning that elastic and gravitational forces should balance contact forces, and all contacts should be in sticking mode (no sliding velocity), and thus each local contact force r i should belong to the interior of the Coulomb friction cone K µ i .

Compared to Section 4.2, new parameters enter the game: the n friction coecients µ i corresponding to the n contacts. Material parameters (stiness and linear mass density) are also augmented with the number of bers, and likewise the intrinsic curvature vector q, which is of same size as q. Finally, as in the direct problem, the n contact forces r i collected in vector r are part of the unknowns.

To simplify the problem, we have so far only considered q and r as unknowns. All the N F material parameters K j 3 and ρS j as well as the n friction coecients µ i are thus assumed to be xed in advance by the user (and explicit dependence to those denoted by brackets will be omitted from now on). Our problem amounts to searching for intrinsic shapes q and contact forces r which are compatible with Equilibrium (4.8),

q = q + K -1 G(q) -H(q) r r ∈ K µ with K µ = K µ 1 × K µ 2 × . . . × K µn .
(4.9)

One immediate observation is that Problem (4.9) is underdetermined, as any valid choice for r leads to a mathematically valid solution for q (but not the reverse).

In particular, it is possible to choose r = 0, meaning that contact forces are inactive; one is then left with the same solution for q as in the isolated case treated in Section 4.2. However, this choice may be fairly unrealistic in many situations, as illustrated in Figure 4.9, where it is perceptible that contact denitely plays a role in the observed conguration.

Figure 4.9 Two examples of input congurations q that should, intuitively, be accounted for by contact forces and not just by the intrinsic shape q of the rods.

A simple heuristics for a well-posed problem

Our idea to better pose the inverse problem (4.9) is to prevent the intrinsic shape q from taking any fancy value, and instead to restrict its value to be close to a good guess q0 . For instance, in Figure 4.9, it is unrealistic to think that the sudden curvature deviation of the curves is due to a brutal change in the intrinsic shape q, especially if we have in mind representing natural bers (hair, plant stem) which are rather homogeneous; instead, a more realistic assumption is to consider that q is probably vanishing, and thus, that the resulting curved shape is mainly explained by contact itself.

More generally, we may formalize our problem as nding the contact forces r which minimize the drift between an admissible intrinsic shape q, which satises 70 Chapter 4. From geometry to mechanics: Inverse design of bers Equilibrium (4.9), and an initial guess q0 . More precisely, our new problem reads min r 1 2 q q -K -1 H(q) r -G(q) -q 2 + γ r 2 s.t. with W = H(q)K -2 H(q) +γI and b = -H(q)K -2 G(q)(q) + K(q -q0 ) , where I is the identity matrix. Note that the objective function f (r) = 1 2 r W r + r b is quadratic, and has to be minimized under second-order conic constraints: our problem is a second-order cone quadratic program (SOCQP), and belongs to the family of the so-called quadratically constrained quadratic programs (QCQP) [START_REF] Boyd | Convex Optimization[END_REF]. Since W is symmetric positive-denite (deniteness comes from the regularization term), our problem is strictly convex and admits a unique solution r. Rening our inverse model by taking a priori information on the parameters we wish to retrieve has thus allowed us, unsurprisingly, to transform the initial underdetermined problem (4.9) into the well-posed convex problem (4.11).

Choosing a good guess For now we have assumed the intrinsic curliness of a given hairstyle could be guessed easily. If one adopts the simplistic assumption that the ber grows regularly with a uniform intrinsic shape, it is indeed possible to retrieve this value approximately. Boundary conditions for thin elastic rods subject to gravity tell us that at the free end of hair bers, the actual curvature equals the intrinsic curvature (see Equations (2.6d) and (2.5)). Provided no contact substantially deforms the tip of the ber, one may thus simply measure the actual curvature at the tip to get a good estimation of the intrinsic curvature.

However, real bers seldom feature a purely uniform intrinsic shape, which is in reality modulated by small defects due to a non-perfect growth process or history involving irreversible (plastic) changes. When dealing with input data coming from real hair captures, we found out that taking as an estimation the full actual conguration q yielded better results than taking a uniform intrinsic curvature matching q at the tip only. Looking back at our minimization problem (4.10), this choice can be interpreted mechanically. Taking q0 = q means that we are searching for contact forces that guarantee an exact equilibrium state for hair while minimizing hair internal elastic energy. This implies that, as far as possible, we rely on the contact forces to compensate for gravity. Of course, as the number of contact points is arbitrarily sparse and as the contact forces are bounded and constrained to belong to the 4.4. From a geometric curve to a material curve 71 friction cone, contact is unlikely to support any arbitrary input shape on its own.

In the (common) case where contacts are not sucient to compensate for gravity, then the elastic energy becomes active and also contributes to the equilibrium.

Ecient solving The optimality conditions of our SOCQP (4.11) are equivalent to the following complementarity condition,

∀i = 1 . . . n, K 1 µ (∇f (r)) i ⊥ r i ∈ K µ i .
(4.12)

Interestingly, this condition is mathematically similar to the complementarity condition (3.9) of the dynamic case, except that the dynamic variable ũi has been replaced with (∇f (r)) i , i.e., the i th 3 × 1 block of the gradient. Now, our solver derived in [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF] was designed to solve exactly this kind of conic complementarity, i.e., problems which can be formulated as

∀i = 1 . . . n, K 1 µ i (Wr r r + c) i ⊥ r r r i ∈ K µ i ,
where W is a symmetric positive (semi-)denite matrix of size (3 n, 3 n), and c a vector of size 3 n. In the direct (dynamic) case, W was identied to the sparse (discrete) Delassus operator HM -1 H [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]. In our static inversion problem, W is identied to our W matrix, which is sparse, symmetric and positivedenite.

We can thus apply our DFCP's Gauss-Seidel solver straightforwardly to solve our inverse problem (4.12), and thus our initial SOCQP (4.11). Unlike the dynamic problem (3.5), our inverse problem (4.12) is convex. It is thus likely to be easier to solve. In practice, simply using the primary Fischer-Burmeister solver (with no fail-safe) proved to converge well for all the problems we have tested.

4.4 From a geometric curve to a material curve

One important problem that we have put under the carpet so far is the following:

Given a piecewise helical (or clothoidal) curve, how is it possible to reconstruct the material frame R(s) or equivalently, the material curvatures and twist q of the corresponding super-model? This question is indeed fundamental, because it raises the issue of transforming a merely geometric curve, i.e., our input curve, into a material curve, i.e., the geometry of a (discrete) Kirchho rod, composed of a centerline r(s) and of a material frame R(s). While the input curve only carries geometric information, the material rod carries some mechanical information, such as how matter deforms along the centerline.

To investigate this question, we shall rst have a look at the continuous picture, taking as input a general curve Γ(s) with no prescribed shape. We shall establish the necessary and sucient conditions enforced by the (strong) Kirchho equations on the framing of Γ(s), both in the direct case (i.e., when the intrinsic curvature κ(s) is known), and in the inverse case (i.e., when the κ(s) is not known and looked for to guarantee the equilibrium of the rod). Then we shall examine the special case of our super-models, as well as that of other discrete models of the literature, in particular the discrete elastic rod model [START_REF] Bergou | Discrete elastic rods[END_REF][START_REF] Bergou | Discrete viscous threads[END_REF].

Note that this section partly covers unpublished material [START_REF] Bertails-Descoubes | About curves, frames, and material rods[END_REF], and we present here only the main conclusions of our study. Since the defense of this habilitation, more results on the inverse design problem have been obtained

and published in [START_REF] Bertails-Descoubes | Inverse design of an isotropic suspended kirchho rod: theoretical and numerical results on the uniqueness of the natural shape[END_REF].

4.4.1 From a geometric curve to a material curve Let Γ(s) be a C 2 -smooth curve8 of length L parameterized by its arclength s ∈ [0, L].

The idea is to equip Γ(s), which is a merely geometric curve, with a material frame R(s) = (n 0 (s), n 1 (s), n 2 (s)), supposed to be orthonormal and adapted to the curve Γ(s), i.e., such that n 0 (s) = Γ (s) and the two other axes n 1 and n 2 lie in the normal plane P ⊥ (s) of the curve (see Figure 4.10(a)). Then, we shall further assume that R(s) is C 1 -smooth and that the frame R(s) is transported along the curve Γ(s) through innitesimal rotation about a Darboux vector Ω Ω Ω(s), just similarly as in Section 2.2 of Chapter 2. This mathematically reads

R (s) = [Ω Ω Ω(s)] × R(s) = R(s) [κ(s)] × ,
where κ(s) = {κ 0 (s), κ 1 (s), κ 2 (s)} is the material curvature vector, which collects the local coordinates (the so-called material twist and curvatures) of the Darboux vector in the material frame R(s).

Note that at this stage, we have already converted our input curve Γ(s) into a material curve {r(s); R(s)} with only material twist and curvatures as degrees of freedom. This however does not mean it is an acceptable frame for a Kirchho rod, as Kirchho equations impose further assumptions about the material frame. We shall now examine such conditions, in two dierent settings:

Direct case: We shall study the law R(s) should obey when Γ(s) corresponds to the centerline of a Kirchho rod, with all parameters xed; Inverse static case: We shall study the law the intrinsic curvatures κ(s) should obey when Γ(s) corresponds to the centerline of a Kirchho rod at equilibrium.

Curve-angle parameterization

Let us rst reduce the dimension of our problem, by observing that R(s) can only rotate about the tangent vector n 0 (s), since it should remain orthonormal and adapted to Γ(s). This degree of freedom can be parameterized by an angle θ(s) representing the rotation angle about n 0 (s) between a reference material frame R ref (s), also orthonormal and adapted to Γ(s), and our material frame R(s). We further suppose that the reference frame R ref (s) can be fully determined by the curve Γ(s): this is typically the case for the Frenet frame and the Bishop frame 9 , which are explicitly taken as examples of reference frames hereafter.

Finally, we assume the reference frame to be suciently smooth for being interpreted as a material frame transported along the curve Γ(s) through innitesimal rotation about a Darboux vector Ω Ω Ω ref (s), similarly to R(s). This allows us to dene a reference material curvature κ κ κ ref (s), which gives the coordinates of Ω Ω Ω ref (s) in the local frame R ref (s).

After a few derivations, we can express our (unknown) curvature vector κ function of the (known) curvature κ κ κ ref , as (4.14) 9 Note that unlike the Frenet frame which is uniquely determined by the curve (provided some extra regularity assumptions for Γ(s) are met, as mentioned earlier), the Bishop frame is not unique but determined by the curve up to a rigid rotation: the initial frame at s = 0 may be arbitrarily chosen.

κ(s) = R -θ (s) κ κ κ ref (s) + θ (s)e x ,
74 Chapter 4. From geometry to mechanics: Inverse design of bers which can be expanded into Examples of reference frames The Frenet frame is a possible reference frame, which admits the Darboux vector Ω Ω Ω ref (s) = τ g (s) n f 0 (s) + κ g (s)n f 2 (s) and the material curvature vector κ κ κ ref (s) = (τ g (s), 0, κ g (s)), with n f 2 (s) the binormal of the curve Γ(s), and κ g (s) and τ g (s) the geometric curvature and torsion, respectively.

       κ 0 (s) = κ ref 0 (s) + θ (s) κ 1 (s) = cos θ(s) κ ref 1 (s) + sin θ(s) κ ref 2 (s) κ 2 (s) = -sin θ(s) κ ref 1 (s) + cos θ(s) κ ref 2 (s).
From (4.15) we retrieve the equations

     κ 0 (s) = τ g (s) + θ (s) κ 1 (s) = sin θ(s) κ g (s)
κ 2 (s) = cos θ(s) κ g (s). which were already derived directly in [Love, 1927, Section 253].

Another popular reference frame is the Bishop frame, which has the particularity to have a vanishing material twist. It has for Darboux vector Ω Ω Ω ref (s) = κ g (s)n f 2 (s) and for material curvature vector κ κ κ In particular, we nd that the material twist κ 0 (s) of the material frame R(s) is given by the derivative of the angle θ(s). This nice property made the curveangle (Γ, θ) parametrization of a material rod from the Bishop frame particularly appealing [START_REF] Langer | Lagrangian aspects of the Kirchho elastic rod[END_REF] and is now preferred by many authors to the classical Euler angles parametrization. This reduced curve-angle parametrization was in particular leveraged in Computer Graphics for designing the so-called discrete elastic rod model, for which primary variables are discrete node positions of the centerline [START_REF] Bergou | Discrete elastic rods[END_REF].

ref (s) = (0, κ b 1 (s), κ b 2 (s)) with κ b 1 (s) and κ b 2 (s) two functions satisfying κ b 1 (s) 2 + κ b 2 (s) 2 = κ 2 g (s). When R ref (s)

Direct problem: framing a Kirchho rod

Let us now state the condition on R(s) guaranteeing it is a suitable frame for a Kirchho rod whose centerline exactly coincides with Γ(s). As before, we assume the Kirchho rod to be clamped at s = 0, and free at s = L. We also suppose that all its parameters κ(s), K, and ρS are known, and that the rod is only subject to a known external force such as gravity. Then from the curve Γ(s), it is possible to recover an admissible material frame (and thus admissible material twist and curvatures) for the corresponding Kirchho rod, according to the following property,

A(s) = 2 (K 2 -K 1 ) κ ref 1 (s)κ ref 2 (s) B(s) = (K 2 -K 1 ) ((κ ref 2 (s)) 2 -(κ ref 1 (s)) 2 ) C(s) = K 1 κ ref 2 (s)κ 1 (s) -K 2 κ ref 1 (s)κ 2 (s) D(s) = -K 1 κ ref 1 (s)κ 1 (s) + K 2 κ ref 2 (s)κ 2 (s) E(s) = -(K 2 -K 1 ) κ ref 1 (s)κ ref 2 (s) + K 0 (κ ref 0 ) (s) -κ 0 (s) ,
subject to the boundary conditions (2.6).

(ii) θ(s) is a stationary point of the potential elastic energy of the rod,

E el = L 0 (κ(s) -κ(s)) K 3 (κ(s) -κ(s)) ds. (4.19)
Moreover, in the particular case of an isotropic rod (K 1 = K 2 ) with vanishing intrinsic curvatures and twist, then the material frame is necessarily the Bishop frame.

Proof. We only sketch here the main lines of the proof. Conditions on the material frame are given by the Kirchho angular momentum equation (2.4b), which can be expressed in the local frame R(s) as with boundary conditions (2.6) for θ(s).

[κ(s)] × K 3 (κ(s) -κ(s)) + K 3 κ (s) -κ (s) + e x × T loc (s) = 0 (4.20) Projecting (4.20) onto each vector of R(s) leads to        κ 1 κ 2 (K 2 -K 1 ) + K 1 κ 2 κ1 -K 2 κ 1 κ2 + K 0 κ 0 -κ 0 = 0 κ 0 κ 2 (K 0 -K 2 ) + K 2 κ 0 κ2 -K 0 κ 2 κ0 + K 1 κ 1 -κ 1 -T loc,2 = 0 κ 0 κ 1 (K 1 -K 0 ) + K 0 κ 1 κ0 -K 1 κ 0 κ1 + K 2 κ 2 -κ 2 + T loc,1 = 0 (4.
Then, we consider the total potential energy of the rod and examine how it gets perturbed as the material frame is perturbed through an innitesimal rotation δϕ(s) about the tangent vector to the centerline Γ(s). The dierence with the general analysis in [Audoly and Pomeau, 2010, Section 3.6] is that we restrict here δϕ(s) to be a cross-sectional perturbation, as the tangent vector of the material frame should not be perturbed since the centerline is xed. Thus, δϕ(s) may be written as δϕ(s) = δθ(s) n 0 (s), which represents an innitesimal rotation of angle δθ(s) about the tangent vector n 0 (s). One may then show that for a rod with free end, the perturbed potential energy boils down to the perturbed potential elastic energy, and also that it vanishes for any small perturbation δϕ(s) if and only if Equation (4.18) subject to the boundary conditions (2.6b2.6d) is satised.

Finally, by choosing the reference frame to be the Bishop frame, one may easily show that when K 1 = K 2 and κ = 0, then Equation (4.1) boils down to θ (s) = 0, meaning that the unique set of admissible material frames is given by the Bishop frame.

This quasistatic frame property, and in particular the interpretation of the material frame as a minimizer of potential elastic energy, has already been leveraged in [START_REF] Bergou | Discrete elastic rods[END_REF] for building the discrete elastic rod model (see Section 4.4.5).

To our knowledge, the optimality condition (4.1) has however not been explicitly derived before, furthermore with respect to a generic angular function θ(s).

Note that a material frame retrieved by Property 4.1 may be that of a Kirchho rod under motion, i.e., with a non-vanishing linear acceleration. Moreover, if the potential elastic energy of the rod is non convex, there may exist several admissible material frames, yielding dierent dynamical states for the rod. However, the constant feature is that whatever the dynamic state of the rod is, the material frame always obeys a quasistatic equation, due to the Kirchho assumption neglecting torsional inertia. Finally, it is interesting to note that Equation (4.1) takes some familiar forms in special cases. For instance, when the rod is isotropic with a piecewise helical centerline and uniform intrinsic curvatures and twist κ, then θ(s) satises the equation of the Euler elastica. We shall further study these special cases in future work.

Inverse static design problem

We now assume that the intrinsic curvatures and twist κ(s) are unknown too (other parameters K and ρS being still xed). From the input curve Γ(s), we aim at nding conditions on κ(s) and the material frame R(s) so that Γ(s) coincides with the centerline of a Kirchho rod at equilibrium under gravity, framed by R(s). We have the following property, which, to the best of our knowledge, has not been formulated before, Property 4.2. For any choice of function θ(s) (i.e., any choice of the material frame R(s)), Kirchho 's equilibrium equations are satised if and only if κ(s) satises the explicit linear ODE of rst order,

κ (s) + K -1 3 [κ(s)] × K 3 κ(s) = κ (s) + K -1 3 b(s), (4.22) with b(s) =   κ 1 (s) κ 2 (s) (K 2 -K 1 ) κ 0 (s) κ 2 (s) (K 0 -K 2 ) -T loc,2 κ 0 (s) κ 1 (s) (K 1 -K 0 ) + T loc,1   ,
where κ(s) is a function of θ(s) given by (4.14), κ(s) = R -θ (s) κ κ κ ref (s) + θ (s)e x and the tension T(s) = R(s) T loc (s) satises the linear momentum equation (2.4a) at equilibrium,

T (s) -ρSge z = 0, i.e., T (s) = - L s ρSge z = ρSg(s -L) e z .
Moreover, we have the boundary condition κ(L) = κ(L). Proof. Proof follows from Equations (4.21), which are reorganized function of the new unknown κ(s) = {κ 0 (s), κ1 (s), κ2 (s)} as

-K 3 κ (s) + A(s) κ(s) = -K 3 κ (s) -b(s), (4.24) with A(s) =   0 K 1 κ 2 (s) -K 2 κ 1 (s) -K 0 κ 2 (s) 0 K 2 κ 0 (s) K 0 κ 1 (s) -K 1 κ 0 (s) 0   = -[κ(s)] × K 3 ,
and then multiplied by -K -1 3 and simplied.

The above property allows us to characterize the set of solutions κ(s): it is the set of (unique) solutions to the Cauchy problem (4.24 4.23) parameterized by the function θ(s). At this stage however, one may feel that it is not possible to identify the intrinsic shape of the rod from the sole data of a curve Γ(s). But actually, after the defense of this habilitation, we have proved that although the set of admissible natural curvatures κ(s) is innite, the intrinsic shape of the curve Γ(s) is unique, and can be computed easily, starting from any framing of the input curve Γ(s) [START_REF] Bertails-Descoubes | Inverse design of an isotropic suspended kirchho rod: theoretical and numerical results on the uniqueness of the natural shape[END_REF].

Discrete picture

Having in mind the continuous picture for framing a Kirchho rod, it is now interesting to get back to discrete Kirchho models. We especially study the framing of two models built upon very dierent concepts: the discrete elastic rod model [START_REF] Bergou | Discrete elastic rods[END_REF][START_REF] Bergou | Discrete viscous threads[END_REF], and the super-helix model.

Framing a discrete elastic rod The quasistatic frame property 4.1 was a core motivation for building the discrete elastic rod model [START_REF] Bergou | Discrete elastic rods[END_REF]. In this model, discrete positions on the centerline are the main degrees of freedom of the dynamics, and discrete angles θ i , determining the material frame from the Bishop frame, are computed at each time step so as to minimize the potential elastic energy of the rod. Note that in [START_REF] Bergou | Discrete viscous threads[END_REF], the θ i may be free degrees of freedom in the case where torsional inertia is considered. In that case the quasistatic frame property is satised at equilibrium only (in the limit of a vanishing torsional inertia, the property is fully recovered). For the discrete elastic rod model, Property 4.1 is thus satised by construction, at any resolution. Moreover, as predicted by Property 4.2, an input curve might be arbitrarily framed, and a set of intrinsic curvatures (that should all be equivalent, given our latest results [START_REF] Bertails-Descoubes | Inverse design of an isotropic suspended kirchho rod: theoretical and numerical results on the uniqueness of the natural shape[END_REF]), can be found to satisfy the equilibrium of the corresponding rod.

Framing a super-helix In contrast, the super-helix relies on a strong assumption regarding the shape of its elements, as it imposes a piecewise uniform shape function for the material curvatures and twist. This assumption has a strong impact on the shape of the material frame. Indeed, from (4.16) we get that on each element, θ(s) should be a uniform function, and as a consequence τ g (s) and κ g (s) should also be uniform. We then retrieve the fact that a super-helix element represents a circular helix (since its Frenet curvature and torsion are uniform), plus the fact that the angle between the material frame and the Frenet frame is uniform. The latter condition implies that the two Darboux vectors are equal (but not the frames), that is,

Ω Ω Ω SH (s) = Ω Ω Ω f (s) ∀s ∈ [0, L]. (4.25)
We note that the θ(s) = C te condition breaks the continuous law (4.1) written in the case of a helical centerline. The super-helix model thus does not satisfy the quasistatic frame condition exactly, but only at the limit when the number of elements tends to innity.

When dealing with inverse design, once a piecewise helix has been reconstructed by some geometric approximation (see Section 4.2.4), then no much choice is left for the material frame: it may start arbitrarily, but then should be transported on each element like the Frenet frame, while remaining continuous at joints (thus shifting from the Frenet frame at each joint). Interestingly, the choice of an arbitrary shape function for the material strains thus acts as a selecting process for the set of admissible material frames. It would be interesting to investigate the case of super-clothoids, which in turn impose a piecewise linear shape function for material 4.5. Conclusion: Towards non-invasive identication 79 strains. Such an assumption is often compatible with real bers for which the growing process remains quite smooth, and would thus pre-lter admissible material frames in a more realistic way.

Conclusion: Towards non-invasive identication

In this chapter we have introduced the inverse static design problem for bers, and have proposed a full pipeline for nding some solutions in the case of the superhelix model, possibly subject to frictional contact. We have also brought some new insights on the full space of solutions, both in the discrete and continuous settings.

Since our rst work on ber inverse design [START_REF] Derouet-Jourdan | Stable inverse dynamic curves[END_REF], many other authors in Computer Graphics got interested in the problem, usually for more general deformable objects discretized with low-order nite elements [START_REF] Twigg | Optimization for sag-free simulations[END_REF][START_REF] Skouras | Computational design of rubber balloons[END_REF][START_REF] Zhao | Interactive authoring of simulation-ready plants[END_REF][START_REF] Chen | An asymptotic numerical method for inverse elastic shape design[END_REF][START_REF] Pérez | Design and fabrication of exible rod meshes[END_REF]. This growing interest has been especially boosted by the advent of additive fabrication, making it possible to fabricate the models with the predicted properties, and then check for their validity in practice.

However, the latter methods, heavily relying upon nonlinear optimization, are meant to converge to one particular solution of the inverse problem. As we have shown here, there can be many ways to explain a given static conguration. In contrast to concurrent methods, we would like to better characterize the exact space of solutions of the problem, in order to understand the precise mechanical information which can be extracted from the geometry of objects. Such a non-blind methodology has already proven useful for computing solutions to the inverse design problem eciently, without having to resort to nonlinear optimization. In the longer run, it would also serve to control the exact amount of input information needed (e.g., a minimal number of poses) for dealing with an accurate identication process. This dream for designing a merely non-invasive identication pipeline, that is, recovering all physical parameters of an object from the observation of its shape only, is the very topic of my research program, which is provided in next chapter. Acquisition technologies range from expensive structured light or laser scans to new cheap devices such as depth cameras [START_REF] Chang | Dynamic geometry processing[END_REF] which are often sucient for capturing a static pose precisely.

Yet, while more and more geometrical data is collected and shared among the communities, there is currently very little study about how to infer the underlying mechanical properties of the captured objects merely from their geometrical congurations. One can however suspect that the pure static shape of a physical object may already give some insights about the constitutive material of the object and the interplaying contacts: from the folding patterns of a tablecloth or a curtain, the human eye may perceive whether the fabric is made of rough cotton or silk, and identify zones of contacts. One may then have the dream that feeding a well-designed physics-based simulator with such easy-available initial data could help predict the deformations or even the dynamics of the physical objects of interest.

Material tests for measuring physical parameters In parallel, contactless measurement methods, which reconstruct full-displacement elds based on camera capture and digital image correlation, have recently gained much interest in Experimental Mechanics [START_REF] Avril | Overview of identication methods of mechanical parameters based on full-eld measurements[END_REF]. Indeed, unlike sensor-based capture, image-based capture does not interfere with the displacement eld being measured.

Combined with FEM-based inverse modeling, contactless measurement methods allow for a complete parameter identication of complex materials. They however request that a number of specic material tests (e.g., tensile and shear tests) be performed, which may often require some expensive material and time-costly measurement protocols, and sometimes may even be impracticable when objects are not directly manipulable. Moreover, although some recent developments in Computer

Graphics have extended the range of studies from small to moderate 3D deformations [START_REF] Miguel | Data-driven estimation of cloth simulation models[END_REF] and partly lightened the necessary amount of control in the experimental setup [START_REF] Wang | Data-driven elastic models for cloth: modeling and measurement[END_REF][START_REF] Miguel | Modeling and estimation of internal friction in cloth[END_REF], such methods remain 5.1. From geometry to mechanics, a broken pipe 83 limited to the study of contact-free objects. In contrast, using geometrical acquisitions of a minimum number of uncontrolled static poses would release the burden of material testing and provide a breakthrough in the non-invasive and fast measurement of many mechanical features including rest conguration, material parameters, self-contacting forces, and friction coecients at contact.

Inaccurate reconstruction of motion Besides, while geometry acquisition from a static pose has become a mature technology, a current major challenge in Computer Graphics and Vision deals with the accurate reconstruction of motion. More precisely, the non-rigid shape registration, which aims at nding consistent correspondences between successive acquisitions of the moving object and resolving ambiguities or occlusions, is still an active area of research [START_REF] Chang | Dynamic geometry processing[END_REF].

Physical predictions of the motion of the physical object could be of great help to resolve such inconsistencies. Conversely, the recovery of accurate dynamic geometrical poses would enrich mechanical interpretation and greatly help rene the identication of all physical parameters at work.

Giving a physical meaning to mere geometrical data would not only serve as an innovative parameters measurement method, but also as a powerful strategy to convert a purely descriptive approach into a generative one, able to predict an innite number of new and rich dynamic scenarios. Some sparse connection tentatives A few works have attempted to identify the physical parameters of slender objects such as bers, cloth or skin, from an arbitrary shape at static equilibrium [START_REF] Bhat | Estimating cloth simulation parameters from video[END_REF][START_REF] Twigg | Optimization for sag-free simulations[END_REF][START_REF] Bickel | Physical face cloning[END_REF]. However, because these approaches rely upon nodal mechanical models whose rest shape has to conform to some constraints (e.g., spring rest length should be nonnegative), they cannot guarantee that a static equilibrium conguration will match the input shape exactly. Moreover, as the functions to minimize are non convex and of large size, the cost of problem solving tends to get prohibitive, precluding an accurate treatment of contact and friction. Finally, such a straight numerical treatment provides only a single solution to the (under-determined) inverse problem and does not give insights about the general structure of the subspace of solutions.

Even sparser, inverse dynamic studies greatly suer from the limited quality of the 3D reconstruction of moving geometry. To complete the geometrical reconstruction consistently, physically-guided approaches have started to be explored. So far however, they have only proven to be successful in the case of smooth dynamic systems such as Newtonian uids [START_REF] Wang | Physically guided liquid surface modeling from videos[END_REF]. Recently, [START_REF] Bouman | Estimating the material properties of fabric from video[END_REF] have proposed to skip the 3D reconstruction step and to use instead statistics characterizing temporal textures in order to predict the material properties of fabric from gentle motion. However, their method is not tailored for making a perfect match between the real phenomenon and a (predictive) physic-based model of it, 84 Chapter 5. Research Perspectives: From Geometry to Mechanics nor for retrieving contact forces and friction coecients at play.

To the best of my knowledge, our preliminary work on the inversion of ber assemblies was the rst to take into account the presence of dry frictional contact for (static) inversion [Derouet-Jourdan et al., 2013a]. We however made some important assumptions, such as the a priori knowledge of the material parameters and an a priori warm start estimation of the rest shape.

Interpreting geometry as a mechanical state

The key challenge of my research project is the automatic connection between the geometrical shape of physical objects and their underlying mechanical properties.

More precisely, I intend to focus this study on complex deformable objects featuring detailed geometrical congurations. Typical objects of interest include slender deformable structures such as rods, plates and shells, all of them being widespread in our environment, from the macroscopic scale (e.g., tree branches and leaves, hair, cloth, skin, paper) to nanoscopic and molecular scales (e.g., carbon nanotubes, DNA). Such structures, which are prone to strongly nonlinear behaviors as well as to possibly prominent self-contacting causing knots, plectonemes, or folds, exhibit very rich geometrical congurations.

My claim in this research project, supported by the preliminary results we have gained on hair bers after years of research, is that these complex geometrical features reveal a lot about the underlying mechanical structures.

Scientic problems

To be able to extract such mechanical properties from a minimum set of geometrical shapes, my goal is to develop specic eective computational models addressing four major scientic problems: SP1: Design of well-suited discrete models for slender structures I believe that the quality of the upstream, reference physics-based model is essential to the effective connection between geometry and mechanics. Typically, such a model should properly account for the nonlinearities due to large displacements of the structures.

It should also be parameterized and discretized in such a way that inversion gets simplied mathematically, possibly avoiding the huge cost of large and nonconvex optimization. In that sense, unlike concurrent methods which impose inverse methods to be compatible with a generic physics-based model, I instead advocate the design of specic physics-based models which are tailored for the inversion process.

More precisely, from our experience on ber modeling, I believe that reduced Lagrangian models, based on a minimal set of coordinates and physical parameters (as opposed to maximal coordinates models such as mass-springs), are particularly well-suited for inversion and physical interpretation of geometrical data. Furthermore, choosing a high-order coordinate system (e.g., curvatures instead of angles) 5.3. Scientic problems 85 allows for a precise handling of curved boundaries and contact geometry, as well as the simplication of constitutive laws (which are transformed into a linear equation in the case of rods). SP2: Static inversion of physical objects from geometrical poses I believe that pure static inversion may by itself reveal many insights regarding a range of parameters such as the undeformed conguration of the object, some material parameters or contact forces.

The typical settings that I will consider will be composed of, on the one hand, a reference mechanical model of the object of interest (SP1), and on the other hand a single or a series of complete geometrical poses corresponding each to a static equilibrium under gravitational load. The core challenge will consist in analysing theoretically the amount of information that can be gained from one or several geometrical poses, and to understand how the fundamental under-determinacy of the inverse problem can be reduced, for each unknown quantity (parameter or force) at play. Both the equilibrium condition and the stability criterion of the equilibrium will be leveraged towards this goal. SP3: Dynamic inversion of physical objects from geometrical poses To rene the solution subspaces searched for in SP2 and estimate dynamic parameters (e.g., damping coecients), a dynamic inversion process accounting for the motion of the object of interest is necessary.

In contrast to the static case SP2 where we can aord to rely on exact geometrical poses, our analysis in the dynamic case will have to take into account the imperfect quality of input data with possible missing parts or outliers. One interesting challenge will be to combine our physics-based model from SP1 together with the acquisition process in order to rene both the parameter estimation and the geometrical acquisition. SP4: Experimental validation with respect to real data The goal will be to confront the theories developed in SP2 and SP3 to real experiments. Compared to the statics, the dynamic case will be particularly involving as it will be highly dependent on the quality of input data as well as the accuracy of the motion predicted by the physics-based model designed in SP1. Such experiments will not only serve to rene the models developed in SP1, SP2 and SP3, but will also be used to improve the 3D geometrical acquisition of moving objects. Besides, once validation will be performed, we shall work on the setting up of new non-invasive measurement protocols to acquire physical parameters of slender structures from a minimal amount of geometrical congurations. Scientic impact Fundamentally, this research project may considerably enhance our understanding of how geometry is tightly linked to mechanics. From the mere observation of shapes through imaging, it will allow scientists to better identify and even predict the physical behavior of many physical and biological slender structures ranging from soft tissues to DNA. Regarding Computer Science, our project may signicantly contribute to pave the way for robust 3D dynamic reconstruction of geometry for slender structures, a longstanding challenge in Computer Vision.

Applications of the future In the long run, I expect this project to contribute to new performant and non-invasive measurement systems for capturing automatically the physical properties of slender structures, from a few samples geometry only.

Beyond huge applications in the measurement of human features such as hair, skin, cloth and even internal thin organs from 3D capture, one interesting application would be the automatic retrieval of physical objects properties from 2D images such as photographs or even paintings. Furthermore, the ability to predict the motion of objects from a mere static geometry is highly desirable not only in the artistic design of dynamic scenes for movies or games, but also in the dynamic virtual try-on industry of garments and hairstyles. Ultimately, our contributions could serve for the fast, image-based creation of complete dynamic digital human doubles, desirable not only for the virtual acting in movies but also for patient-specic surgery.

Chapter 6

Conclusion

Now listen to the rule of the last inch. The realm of the last inch. The job is almost nished, the goal almost attained, everything possible seems to have been achieved, every diculty overcome -and yet the quality is just not there. The work needs more nish, perhaps further research. In that moment of weariness and self-satisfaction, the temptation is greatest to give up, not to strive for the peak of quality. That's the realm of the last inch -here, the work is very, very complex, but it's also particularly valuable because it's done with the most perfect means. The rule of the last inch is simply this -not to leave it undone. And not to put it o -because otherwise your mind loses touch with that realm. And not to mind how much time you spend on it, because the aim is not to nish the job quickly, but to reach perfection.

Alexander Solzhenitsyn, The First Circle

In this memoir we have presented our main ideas for designing realistic, robust and computationally ecient models for thin deformable objects and divided materials prone to frictional contact. The general philosophy we are tending to adopt is a pluridisciplinary approach, driven by applications but also inspired by theoretical concepts. To build a consistent and eective model, we strive to master the full modeling pipeline, from the continuous physical model to its discretization and implementation. From this transverse viewpoint, which led us to meet dierent scientic communities, we retain a number of clues and preferred paths for advancing our research work; these are listed below.

Upfront modeling Throughout this long-term work on the numerical modeling of bers and frictional contact, we have learnt that systematically concentrating the eorts on the upstream modeling and formulation of problems often pays o: even for very complex problems, the resulting numerics may be greatly simplied and thus solved more easily and robustly. Keeping in mind this key lesson, we are starting to investigate the case of 2D slender structures (plates and shells), for which many exciting challenges remain open, regarding both direct and inverse modeling.

Small data for big understanding 1 My current research trajectory is clearly an attempt to build stronger connections with physics and physicists, with the aim to better understand and master the macroscopic modeling of complex natural phenomena; indeed, there is still a huge gap between our actual need for relevant macroscopic models, and the set of available models of the literature. However, in a time where big data is considerably expanding and seems to prevail over standard 1 First coined by Pascal Barla while talking together at Siggraph 2016.

88

Chapter 6. Conclusion modeling techniques, one may argue that such a tentative, denitely upstreams, is vain. Indeed, big data analysis spreads the word that models may not need to be designed by humans anymore, but rather spontaneously emerge from the automatic analysis and learning of terabytes of collected data, performed by supercomputers.

Yet, in this deep change of modeling strategies where statistical and computational approaches may want to take the lead, I believe that the physical modeling point of view is more than ever desirable. First of all, the human mind feels like understanding the world, and not just reproducing it blindly. Then, the idea of feeding computers with unlimited and inexpensive data is probably illusory, and unsustainable in certain elds like human biology where medical imaging cannot be overused on a given patient. In contrast, I believe that compact models, able to give some physical meaning and to provide some predictive tools from a minimal amount of input data, are still extremely valuable. They should also complement nicely statistical approaches which treat huge amounts of data, by providing them with some useful hints for interpreting classied objects.

Physics and Computer Science: a promising alliance? Finally, as pointed out in the introduction of this document, there are some striking common points between Nonlinear Physics and Computer Graphics, as some phenomena of interest are becoming more and more similar. Links between the two elds are however quite loose, with only a handful of experts in the world working at the frontier between the two elds. Such a disconnection can be explained by higher education, which usually clearly separates Physics (and Mechanics) from Computer Science, the latter being more willingly attached to Applied Mathematics. In France, such a segmentation is clear at university or in engineering schools: computer scientists will not tackle Natural Sciences during their study, and physicists will probably learn that Computer Science is just a programming tool, readily accessible through black boxes such as Mathematica or Abaqus. In research labs, this segmentation persists, although more and more eorts are made to establish connections between digital science and other elds. For instance, the last decades have seen the emergence of Computational Biology, an actual pluridisciplinary eld where biologists, mathematicians and computer scientists have indeed combined their expertise to advance a modern vision of biology.

At my (much humbler) scale, I believe valuable connections could be established between Physics of Nonlinear Elasticity and Computer Animation, where scientists share the common goal to model complex and integrated mechanical phenomena at the macroscopic scale. Indeed, while computer scientists are eager to learn and understand new physical models, physicists get more and more interested in the numerical tools, in which they perceive not only a means to conrm predictions afterwards, but also a support for testing hypothesis and for getting insights into the search for analytic solutions, starting at the modeling stage. However, they may be limited by a blind usage of numerical black boxes, which may not be dedicated to their specic needs. According to me, promoting a science of modeling in numerical physics would thus be a promising and rich avenue for our two research elds. Numerical Modeling of elastic slender structures subject to contact and friction:

From dynamic simulation to inverse static design Abstract: Slender structures (laments or surfaces) are ubiquitous around us.

It is actually striking to note that such thin structures represent most of the deformable objects composing our direct environment, stemming either from natural biological phenomena plant stems and leaves, ower petals, hair, skin or manufactural processes cables, ribbons, paper, cloth. When multiple thin structures are coupled together with contact and friction, the range of emerging phenomena is even more exacerbated, giving rise to stick-slip dynamical instabilities, entangling, or spontaneous collective behavior.

This memoir presents an overview of our work on the simulation of slender elastic structures and entangled materials, such as hair, with a specic interest for virtual prototyping and computer graphics applications.

I rst introduce a family of high-order, reduced models for discretizing Kirchho 's equations for thin elastic rods in a both faithful and robust way. Such models are particularly well-suited for simulating inextensible bers subject to bending and twisting, and featuring an arbitrary curly resting geometry. Then I show how such models can be coupled to frictional contact using the nonsmooth contact dynamics framework, and I present a hybrid iterative solver suitable for robustly handling thousands packed bers at reasonable frame rates. A rst continuum-based approach. Finally, I give some insights into the inverse modeling of bers, consisting in taking an arbitrary curve geometry as input and inferring corresponding geometric and physical parameters of the simulator such that the input geometry corresponds to a stable conguration at equilibrium.

Figure 1 .

 1 Figure 1.2 Top: Slender structures such as hair, cables, plants, or cloth are widespread at the macroscopic scale, where their interactions are dominated by dry frictional contact. Bottom: it is interesting to observe that slender structures also appear at extreme scales, such as supercoiled DNA (left) and carbon nanotubes (middle) when zooming in some biological and manufactured matter, or even galactic laments (right) when zooming out up to the universe scale. Image credits: (c) S. Paris, T. Judd, 2007, (d) J. Bahnson's group, (e) L. Qu, U. Dayton, 2008, (f) A. Pontzen, F. Governato, 2014.

Figure 1 . 3

 13 Figure 1.3 The famous Siggraph bar. Nobody knows exactly where it stands, but everyone knows you should go over it to have your paper accepted. Be aware that after a rst success, you will get addicted to it.

Figure 2 . 1

 21 Figure 2.1 Many physical bers exhibit a smooth curled geometry with linear-like curvature prole, which is captured and deformed accurately thanks to our super-clothoid model [Casati and Bertails-Descoubes, 2013]. From left to right and top to bottom, three examples of real bers whose shapes are synthesized and virtually deformed in real-time using a very low number of 3D clothoidal elements: a vine tendril (4 elements), a hair ringlet (2 elements), and a curled paper ribbon (1 single element with an exaggerated at cross-section). Left photograph courtesy of Jon Sullivan, pdphoto.org.

0 M

 0 while considering an inextensible rod as well as the boundary conditions (2.6), leads to the following weak formulation, L (s) + n 0 (s) × T(s) • δθ(s) ds = 0, where T(s) = L s (p(s ) -ρSr(s )) ds . Integrating by parts and noting that δκ(s) = δθ (s) and [δθ(s)] × n 0 (s) = (δr) (s), we get L 0 M(s) • δκ(s) ds + L 0 p(s) • δr(s) ds = ρS L 0 r(s) • δr(s) ds.

  Figure 2.2 Comparison of fairness (or smoothness) between (a) the super-helix model and (b) the super-clothoid model. Whereas joints between elements are particularly visible for the super-helix model (5 elements here), the super-clothoid model generates a smooth, visually pleasing shape even at a very coarse resolution (2 elements here).

Figure 2

 2 Figure 2.4 Hillock-like prole of the general term Rn (s) in norm. In blue: Evolution of Rn (s) ∞ function of n (at xed s), in log scale.As expected[START_REF] Neher | An enclosure method for the solution of linear ODEs with polynomial coecients[END_REF], the decreasing towards 0 appears to be super-linear. In red: Evolution, function of s, of the upper-bound H(s), in log scale. Note that the plot of this upper-bound visually matches the maximum function max n Rn (s)

Figure 2 . 5

 25 Figure 2.5 Our piecewise summation algorithm applied to the rod's kinematics, for one clothoidal element. The length of each subinterval which guarantees a safe evaluation of the geometry is automatically provided by our method.

  Figure 2.6 Comparisons between the super-clothoid model (SSC ), the super-helix model (SH ), and the discrete elastic rod model (DER), in terms of accuracy and eciency.

  Figure 2.7 (a) Interactive manipulation of a long super-helix made of 30 helical elements. (b) Interactive animation of a weeping willow swept by wind, represented as a tree-like structure of super-helices (composed of 242 helical elements) with various stiness values.

  Figure 2.8 (a) Cylindrical geometry of our shell patch with uniform material curvatures. (b-d) Dynamical oscillations of our inextensible shell patch clamped at bottom right corner and subject to gravity, with an initially curved intrinsic shape.

  Figure 3.1 Taking into account dry frictional contact is of utmost importance for capturing signicant visual features such as a stable sand pile or complex hair patterns.

  Figure 3.2 Simulation of complex geometrical congurations involving self-contact, such as a knot (a) and a plectoneme (b). With a classical penalty-based method, tying the knot causes the contact to break, whereas contact resists arbitrarily large user forces when using an implicit constraint-based method, without having to tune an extra stiness parameter.

  Newtonian uid (e.g., water) εT σ T (f ) Yield-stress uid (e.g., foam)

Figure 3 . 4

 34 Figure 3.4 Local contact basis, with normal and tangent subspaces.

  Figure 3.5 Illustration in 3D of the three cases of the Signorini-Coulomb law.

  Figure 3.6 Exact vs. approximate model for Coulomb friction. (a) Periodic orbit of the free end of a rod resting on top of a rotating sphere. Frictional contact is simulated using our eight dierent solvers. Solvers that model exact Coulomb friction all reach the same orbit, no matter the choice of the error function. Other models (dashed lines) reach completely dierent ones. (b) Visual comparison of a full hair simulation with exact (left) and linearized (right) Coulomb friction, for equal computation time. With exact Coulomb friction, spontaneous hair clumping emerges in a more visible way.

  The interpretation of the Signorini-Coulomb law as a complementarity condition on dual cones (3.9) is at the origin of a myriad of other equivalent formulations, each one of them giving rise to dierent solving algorithms. In the following, we shall focus on two new formulations inspired by (3.9): the so-called De Saxcé functional formulation, and Cadoux's xed-point algorithm.

Figure 3 . 7

 37 Figure 3.7 Changes of variables applied onto (a) the relative velocity u and the contact force r, in order to exhibit (b) the complementarity of u and r proposed by De Saxcé and Feng [1998].

  Chapter 3. Frictional contact models for large multibody systems where Π C is the orthogonal projection operator on the convex space C, ξ T and ξ N are positive real numbers, and B d-1 (a) ⊂ R d-1 is the ball of radius a ≥ 0 centered at the origin.As suggested above, another well-known function satisfying (3.10) may be derived based upon De Saxcé and Feng's change of variable (3.8), by using the projection theorem (see, e.g.,[Hiriart-Urruty and Lemaréchal, 2001, Proposition A.5.3.3]),

Figure 3 . 8

 38 Figure 3.8 Comparison of the hair collective behavior between (top) real hair motion sequences and (bottom) our corresponding simulations, based on large assemblies of (up to 2,000) individual bers with massive self-contacts and Coulomb friction. Our model retains typical emerging eects such as transient coherent motions or stick-slip instabilities, at competitive timings (typically, a few seconds per time step, with dt = 1ms).

  continuum. Let S d be the space of symmetric d × d tensors (represented as matrices of dimension d), where d is the space dimension (d = 2 or d = 3). For σ ∈ S d , we denote by Tr σ its trace (normal part) and by Dev σ its deviatoric (tangential) part, Dev σ = σ -1 3 Tr σI, where I is the identity tensor. We use the norm | • | associated to the scalar product < σ, τ >= σ:τ 2 = 1 2 τ ij σ ij . Note that | • | amounts to the Frobenius norm scaled by 1 √ 2

d

  where the rst part corresponds to a standard Newtonian viscosity (dissipative term due to random collisions in the owing material), and σ c is the additional stress due to the Coulomb interactions between individual grains. Drucker-Prager rheology The contact normal stress is p = -1 Tr σ c , such that σ c = Dev σ c -pI. The second stress invariant can be computed through the identity J 2 = 1 2 Tr(Dev σ c 2 ) = | Dev σ c | 2 . The Drucker-Prager yield criterion (3.15) therefore boils down to | Dev σ c | ≤ μp.

  Figure 3.10 Simulation of a granular ow with our nonsmooth numerical method, with μ = 0.3. The velocity and stress elds are shown in the case where the granular ows from left to right in a narrow channel, around a cylindrical obstacle.

  Figure 3.11 A steel ball impacts a granular bed with varying volume fraction. Unlike previous approaches, our fully resolved (non-linearized) Drucker-Prager rheology allows us to retrieve a perfectly round crater. Observe also the solid, liquid and gaseous phases of the granular ow that are retrieved by our model.

4. 1

 1 Figure 4.1 Stylized hand-drawn hairstyle. The personality of the character is largely determined by the shape of her hair, hence such geometric features should be preserved during animation. Image courtesy of Néomis Animation Studio Jeroen Dejonckheere.

  Figure 4.2 Classical initialization of a physics-based simulator.The initial conguration q(t = 0) is set to a trivial intrinsic shape q here, a set of straight rods (a). Then, after feeding the simulator with other parameters (such as the mass and stiness of each rod), a static shape q eq can be computed under gravity and frictional contacts (b).

Figure 4 . 3

 43 Figure 4.3 When setting the intrinsic shape of the hair to the observed shape q obs (a), the simulator yields a new equilibrium conguration (b) which can be far from the observed shape. Acquired geometry (a) is data courtesy of Hao Li and colleagues.

  Figure 4.4 Example of ber design session.The user rst carves the tail of the cat geometrically as she desires (a), then converts this curve into a piecewise helical arc curve (b) using our oating tangent algorithm (see Section 4.2.4), and nally automatically computes the parameters of a super-helix so that its centerline matches the input curve at stable equilibrium under gravity (c).

Figure 4

 4 Figure 4.5, and an example of interactive design session is illustrated in Figure 4.4.

Figure 4 . 5

 45 Figure 4.5 Our inverse design process for an isolated ber.

  Figure 4.7 Comparison of helical approximations (with N = 15 elements) on a input curve dataset (a), between least-squares optimization (b) and our oating tangents algorithm (c).Our method proves not only more accurate, but also faster (10× on this example) than least-squares optimization.

  introduced a regularization term γ r 2 for preventing contact forces to reach overly big values. Problem (4.10) may be put under the canonical form min

  Figure 4.10 Framing an arbitrary curve Γ(s) with an adapted frame R(s), parameterized by the angle θ with respect to a reference frame R ref (s).

  The relationship between R(s) and R ref (s) is illustrated in Figure 4.10(b), and mathematically reads R(s) = R ref (s) R θ (s), θ(s) -sin θ(s) 0 sin θ(s) cos θ(s)   is the rotation matrix of axis e x and angle θ, which represents the change of basis matrix from R ref (s) to R(s).

  is the Bishop frame, System (4.15) boils down to s) = θ (s)κ 1 (s) = cos θ(s) κ b 1 (s) + sin θ(s) κ b 2 (s) κ 2 (s) = -sin θ(s) κ b 1 (s) + cos θ(s) κ b 2 (s).

  21a) represents the axial projection of the Kirchho angular momentum conservation (4.18) expressed in the material frame. Note that this equation is 76 Chapter 4. From geometry to mechanics: Inverse design of bers independent of the tension T(s), and thus of the acceleration of the centerline. It can actually be reformulated function of the single angular variable θ(s), all other quantities depending only of the centerline r(s) (thus being known). This can be simply achieved by replacing the curvatures with their angle parameterization given by formulas (4.15), leading to the second-order nonlinear dierential equation (4.1)

  24) and (4.23) form a Cauchy problem, which admits a unique solution.

  Figure5.1 Recent acquisition technologies make it now aordable to reconstruct the 3D geometry of many complex objects featuring intricate shapes, such as (a) cloth[START_REF] Bradley | Accurate multiview reconstruction using robust binocular stereo and surface meshing[END_REF], (b) skin[START_REF] Beeler | High-quality passive facial performance capture using anchor frames[END_REF], (c) hair[START_REF] Luo | Structure-aware hair capture[END_REF], and (d) trees[START_REF] Livny | Automatic reconstruction of tree skeletal structures from point clouds[END_REF]. However, the consistent connection of such a geometrical description to the mechanical characterization of the underlying object remains an open issue.

86

  Chapter 5. Research Perspectives: From Geometry to Mechanics 5.4 Expected Impact Converting geometry into physics promises a considerable amount of scientic impact and potential applications for a number of elds ranging from Computer Graphics to Physics and Biology. Thanks to the funding of this research project from the European Research Council (grant started in September 2015), I intend to create and lead a group internationally active on this pluridisciplinary topic.

Figure 6 . 1

 61 Figure 6.1 My daily activity, a constant search for the right balance. Between mechanics and graphics, between theory and applications, between paper work and programming, between supervision and learning, between publication and industrial transfer... and between work and private life.

Keywords:

  Thin elastic rods, frictional contact, numerical simulation, direct and inverse modeling Modélisation numérique de structures élancées en présence de contact frottant : De la simulation dynamique à la conception statique inverse Résumé : Les structures élancées (objets liformes ou surfaciques), font partie intégrante de notre environnement quotidien : les exemples biologiques naturels branches d'arbre, feuilles, pétales, cheveux, peau rivalisent en nombre avec les objets issus de l'industrie câbles, rubans, papier, vêtements. La diversité des formes et la non-linéarité des déformations de ces structures procurent une grande richesse visuelle, qui est davantage exacerbée lorsque des contacts et du frottement entrent en jeu, donnant lieu à des phénomènes d'instabilité de type glissement-adhérence, des comportements collectifs spontanés, ou de l'emmêlement. Synthétiser de manière réaliste la forme et le mouvement de ces phénomènes à l'échelle macroscopique intéresse grandement l'industrie du cinéma d'animation et des eets spéciaux, mais aussi, et de plus en plus, des secteurs industriels plus traditionnels en ingénierie mécanique, soucieux de raner et d'optimiser leurs procédés de conception et de fabrication grâce au prototypage virtuel. Ce mémoire présente une synthèse de dix années de recherche visant à capturer quelques-uns de ces phénomènes mécaniques complexes de manière à la fois réaliste, robuste, et ecace en temps de calcul. J'introduis tout d'abord une famille de modèles réduits et de haut degré en espace permettant de discrétiser les équations de Kirchho pour les tiges élastiques minces. Ces modèles sont particulièrement adaptés à la simulation de laments capables de se échir et de se tordre tout en restant parfaitement inextensibles. Après l'évocation de perspectives d'extension au cas surfacique, je montre comment ces structures minces peuvent être couplées entre elles via du contact et du frottement solide, en m'appuyant sur des principes de mécanique non-régulière. En particulier, j'introduis un solveur numérique hybride capable de résoudre ecacement le problème de contact frottant discret, pour des dizaines de milliers de contacts. Ce modèle numérique nous a notamment permis de simuler la dynamique de milliers de bres enchevêtrées, et de commencer à envisager la simulation dynamique de vêtements. Nous avons également pu étendre la portée de notre approche à la résolution du mouvement de millions de grains en contact frottant, modélisés comme un écoulement uide dont la loi de comportement (loi de Drucker-Prager), est non-régulière. Finalement, j'évoque nos travaux en conception statique inverse de bres, consistant à inférer les paramètres physiques de nos modèles à partir de la seule observation de leur forme sous gravité et éventuellement en présence de contact frottant. Au-delà de la perspective d'améliorer le contrôle utilisateur d'une simulation physique dans un contexte de création artistique, ces recherches laissent espérer à long terme la conception de protocoles d'identication paramétrique non-invasifs pour les structures élancées en contact frottant.

  , in particular in Physics and Mechanical Engineering. In parallel, we are considering diverse applications of our work, encompassing not only special eects

	for the movie industry and feature lm animation, but also virtual prototyping
	for the engineering industry (past collaborations with L'Oréal Research, starting
	collaboration with Safran).

Computer Geometric Design, Multibody System Dynamics, or Non-Newtonian Fluid Mechanics. Likewise, we have started to participate in conferences and scientic events held on various 1.4. Structure of this document domains

  It is noteworthy that the local coordinates of Ω in the material frame represent the material twist κ 0 and material curvatures

	κ 1 and κ 2 of the rod, i.e., Ω(s) = R(s) κ(s), where κ(s) = [κ 0 (s), κ 1 (s), κ 2 (s)] is
	called the material curvature vector (or simply curvature vector ) in the remainder of
	the document. By further using properties of rotation matrices, one can reformulate
	Equation (2.1) as

  Prior art and choice for nonsmooth models . . . . . . . . . . 33 3.2 Nonsmooth contact dynamics for thin elastic objects . . . . 36

	Chapter 3
	Frictional contact models for large
	multibody systems
	Contents
	3.1

3.2.1 Adding unilateral contact with Coulomb friction . . . . . . . 36 3.2.2 The Discrete Frictional Contact Problem (DFCP) . . . . . . 39 3.2.3 Multiple equivalent formulations of Signorini-Coulomb . . . . 41 3.2.4 A local hybrid solving strategy . . . . . . . . . . . . . . . . . 46 3.3 Continuum modeling of granular materials . . . . . . . . . . 48 3.3.1 Discrete vs. continuum models for granular matter . . . . . . 48 3.3.2 Macroscopic yield criteria . . . . . . . . . . . . . . . . . . . . 50 3.3.3 Nonsmooth rheology DP(µ) . . . . . . . . . . . . . . . . . . . 51 3.3.4 Extension to a varying volume fraction . . . . . . . . . . . . . 54 3.4 Conclusion: Towards a continuum-based model for brous materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

  Inverse design problem: motivation and goal . . . . . . . . . 58 4.2 Case of an isolated ber . . . . . . . . . . . . . . . . . . . . . 61

	Contents
	4.1

4.2.1 Necessary and sucient condition for equilibrium . . . . . . . 62 4.2.2 Sucient condition for stability . . . . . . . . . . . . . . . . . 63 4.2.3 Inverse design algorithm . . . . . . . . . . . . . . . . . . . . . 64 4.2.4 From an arbitrary curve to a discrete material conguration . 65 4.3 Case of a ber assembly subject to frictional contact . . . . 66 4.3.1 Underdetermined conditions for equilibrium . . . . . . . . . . 68 4.3.2 A simple heuristics for a well-posed problem . . . . . . . . . . 69 4.4 From a geometric curve to a material curve . . . . . . . . . 71 4.4.1 From a geometric curve to a material curve . . . . . . . . . . 72 4.4.2 Curve-angle parameterization . . . . . . . . . . . . . . . . . . 73 4.4.3 Direct problem: framing a Kirchho rod . . . . . . . . . . . . 74 4.4.4 Inverse static design problem . . . . . . . . . . . . . . . . . . 76 4.4.5 Discrete picture . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.5 Conclusion: Towards non-invasive identication . . . . . . . 79

  Property 4.1. {Γ(s), R(s)} is an admissible Kirchho rod conguration, i.e., it is compatible with dynamic Kirchho equations, if and only if one of the two equivalent conditions is satised, (i) θ(s) is solution of the axial projection of the angular momentum equation,

	M (s) • n 0 (s) = 0	(4.18)
	subject to the boundary conditions (2.6b2.6d), meaning that θ(s) is solution of the
	second-order dierential equation
	θ (s)+A(s) cos 2 θ(s)+B(s) cos θ(s) sin θ(s)+C(s) cos θ(s)+D(s) sin θ(s)+E(s) = 0
	with

In France this new trend was particularly stimulated by the work of Yves Pomeau, who con- vinced many young scientists to study the nonlinear physics of common objects such as paper, plants, or hair [

ESPCI, 2016].2 It is however amusing to observe that research in these areas is quite successful in obtaining the IGNobel prize[START_REF] Audoly | Fragmentation of rods by cascading cracks: Why spaghetti does not break in half[END_REF][START_REF] Goldstein | Shape of a ponytail and the statistical physics of hair ber bundles[END_REF], thus still being considered as an exotic research topic by physicists.

In some situations, strictly simulating the known equations of physics may not be eective in conveying the right emotions, even if the underlying model is realistic: the audience may simply neither be convinced nor impressed by the simulated phenomenon! In many cases, a beautication or magnication of the phenomenon is thus desired, in a way comparable to the exaggerated deformed balls drawn by cartoonists to mimick expressive rebound[START_REF] Williams | The animator's survival kit[END_REF]. In the same way expressive rendering has emerged in the 90's, expressive animation has recently gained some interest. Providing user control over a physics-based simulation through inverse design, as we do in the third part of this memoir (see Chapter 4), is one research direction among others of this emerging subeld.

In the whole document, the term discrete refers to nite-dimensional models in space (which are subsequently discretized in time). Our contributions mostly deal with this spatial reduction of the Kirchho equations, yet with some implications on the full space-time discretization scheme.

In this document we often use the shortcut curvatures to refer to material twist and curvatures.

The material frame R(s) is said to be adapted to the curve r(s).

This assumption is debatable in some cases, like for instance vegetal stems for which it has been shown that the stiness may vary spatially during the plant growing process[START_REF] Beusmans | Mechanical properties within the growth zone of corn roots investigated by bending experiments. ii. distributions of modulus and compliance in bending[END_REF]. Note however that such an assumption is not mandatory for building our discrete rod models in Section 2.3. Yet, it will be leveraged by our inversion process presented in Chapter 4.

This means that at each time step, contact forces r (and the stress tensor σ respectively) are unknowns of the problem, together with the relative velocities u of the problem (and the strain rate tensor ε respectively).

Note that the Painlevé paradox, as well as more general cases of nonexistence of solution, do not only occur in the discrete settings, but also in the time-continuous problem of Lagrangian systems subject to frictional contact. This was our very objet of study in[START_REF] Blumentals | The contact problem in Lagrangian systems subject to bilateral and unilateral constraints, with or without sliding coulomb's friction: a tutorial[END_REF].

In the last decades there has been many attempts to solve the DFCP, but the literature still remains very sparse when dealing with convergence guarantees[START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]

In all this document, equilibrium refers to the static equilibrium of the structure, meaning that the sum of forces and the sum of torques applying onto the object vanish.

In practice, our input data take the form of geometric curves in the case of rods, and surfaces in the case of plates/shells. Such primitives are represented with a nite number of degrees of freedom (for instance an ordered sequence of points, a spline, a triangulated mesh, etc.).

At this stage, required degree of smoothness varies upon the reference frame that will be chosen in the curve-angle parameterization (see Section 4.4.2). If the Frenet frame is chosen, C 3 -continuity is required to dene torsion as a continuous function of s (provided curvature does not vanish), allowing us to interpret the Frenet frame as a C 1 -smooth material frame.

Acknowledgments

quadratic optimization problem [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF][START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF] min ũ 1 2 ( u)

which can be solved using our hybrid solver, since this time u plays the exact same role as r (replacing K µ i with K 1 µ i

). Following [START_REF] Cadoux | Méthodes d'optimisation pour la dynamique non-régulière[END_REF], we compute the solution of the full problem (3.5) by iteratively updating s using the xed point algorithm (3.13). Unlike the Delassus operator W, our new operator W is easy to assemble and sparse. Solving our primal problem (3.14) thus turns out to be orders of magnitude faster compared to the reduced DFCP. We are currently trying to extend this simple method for dealing with cloth self-contact.

Figure 3.9 Our body-cloth contact solver captures exact Coulomb friction both eciently and robustly. On this dress example featuring 6000 cloth vertices and 1000 contact points on average, our solver converges at each time step (dt = 1ms) in a few hundred milliseconds..

Continuum modeling of granular materials

In this section we drastically change models for representing a collection of interacting objects. Instead of modeling individual objects explicitly as in Section 3.2, we adopt a continuum viewpoint where a fragment of this continuum represents an homogenized sample of individual elements. Such an approach is desirable when aiming at modeling hundreds thousands or even billions interacting elements, which would be intractable using a discrete element method as before 12 . However, to the best of our knowledge, continuum viewpoints are not well-developed enough for taking into account individual elements with complex geometries such as long and thin deformable bers. To start with a continuum viewpoint, we have thus restricted ourselves to a simpler, yet widely studied material: granular matter.

Discrete vs. continuum models for granular matter

Granular materials (see, e.g., [START_REF] Andreotti | Granular media: between uid and solid[END_REF] for a comprehensive description) commonly refer to a large collection of small solid grains larger than 100 µm 12 Note that another strategy to scale things up would be to take a look at high performance computation. However the solver we have designed in Section 3.2 is inherently sequential, and a completely new parallel algorithm would then have to be designed from scratch, requiring some specic expertise in parallel design. Motivated by physical insights instead, we have rather chosen, at least for now, to investigate the building of new, macroscopic physics-based models.

Chapter 5

Research Perspectives: From

Geometry to Mechanics

With the considerable advance of automatic image-based capture in Computer

Vision and Computer Graphics these latest years, it becomes now aordable to acquire quickly and precisely the full 3D geometry of many mechanical objects featuring intricate shapes. Yet, while more and more geometrical data get collected and shared among the communities, there are currently very few studies about how to infer the underlying mechanical properties of the captured objects merely from their geometrical congurations.

An important challenge consists in developing a non-invasive method for inferring the mechanical properties of complex objects from a minimal set of geometrical poses, in order to predict their dynamics. In contrast to classical inverse reconstruction methods, my research project is built upon the claim that 1/ the mere geometrical shape of physical objects reveals a lot about their underlying mechanical properties and 2/ this property can be fully leveraged for a wide range of objects featuring rich geometrical congurations, such as slender structures subject to contact and friction (e.g., folded cloth or twined laments).

To achieve this goal, we shall develop an original inverse modeling strategy based upon a/ the design of reduced and high-order discrete models for slender mechanical structures including rods, plates and shells, b/ a compact and wellposed mathematical formulation of our nonsmooth inverse problems, both in the static and dynamic cases, c/ the design of robust and ecient numerical tools for solving such complex problems, and d/ a thorough experimental validation of our methods relying on the most recent capturing tools.

In addition to signicant advances in fast image-based measurement of diverse mechanical materials stemming from physics, biology, or manufacturing, this research is expected in the long run to ease considerably the design of physically realistic virtual worlds, as well as to boost the creation of dynamic human doubles. Graphics these latest years, it becomes now aordable to acquire quickly and precisely the full 3D geometry of many mechanical objects featuring intricate shapes such as cloth and skin [START_REF] Miguel | Data-driven estimation of cloth simulation models[END_REF], or even hair bers [START_REF] Luo | Structure-aware hair capture[END_REF].