**Discrete Elastic Rods**Mikl´ os Bergou Max Wardetzky Stephen Robinson Basile Audoly Eitan Grinspun Freie Universit¨ at Berlin Columbia University Columbia University CNRS / UPMC Univ Paris 06 Columbia University Figure 1: Experiment and simulation: A simple (trefoil) knot tied on an elastic rope can be turned into a number of fascinating shapes when twisted. Starting with a twist-free knot (left), we observe both continuous and discontinuous changes in the shape, for both directions of twist. Using our model of Discrete Elastic Rods, we are able to reproduce experiments with high accuracy. Abstract Notably lacking is the application of DDG to physical modeling of elastic rods—curve-like elastic bodies that have one dimension (“length”) much larger than the others (“cross-section”). Rods have many interesting potential applications in animating knots, sutures, plants, and even kinematic skeletons. They are ideal for model- ing deformations characterized by stretching, bending, and twist- ing. Stretching and bending are captured by the deformation of a curve called the centerline, while twisting is captured by the rota- tion of a material frame associated to each point on the centerline. We present a discrete treatment of adapted framed curves, paral- lel transport, and holonomy, thus establishing the language for a discrete geometric model of thin flexible rods with arbitrary cross section and undeformed configuration. Our approach differs from existing simulation techniques in the graphics and mechanics lit- erature both in the kinematic description—we represent the mate- rial frame by its angular deviation from the natural Bishop frame— as well as in the dynamical treatment—we treat the centerline as dynamic and the material frame as quasistatic. Additionally, we describe a manifold projection method for coupling rods to rigid- bodies and simultaneously enforcing rod inextensibility. The use of quasistatics and constraints provides an efficient treatment for stiff twisting and stretching modes; at the same time, we retain the dy- namic bending of the centerline and accurately reproduce the cou- plingbetweenbendingandtwistingmodes. Wevalidatethediscrete rod model via quantitative buckling, stability, and coupled-mode experiments, and via qualitative knot-tying comparisons. 1.1 Goals and contributions Our goal is to develop a principled model that is (a) simple to im- plement and efficient to execute and (b) easy to validate and test for convergence, in the sense that solutions to static problems and trajectories of dynamic problems in the discrete setup approach the solutions of the corresponding smooth problem. In pursuing this goal, this paper advances our understanding of discrete differential geometry, physical modeling, and physical simulation. We build on a representation Elegant model of elastic rods of elastic rods introduced for purposes of analysis by Langer and Singer [1996], arriving at a reduced coordinate formulation with a minimal number of degrees of freedom for extensible rods that rep- resents the centerline of the rod explicitly and represents the mate- rial frame using only a scalar variable (§4.2). Like other reduced coordinate models, this avoids the need for stiff constraints that couple the material frame to the centerline, yet unlike other (e.g., curvature-based) reduced coordinate models, the explicit centerline representation facilitates collision handling and rendering. CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Animation Keywords: rods, strands, discreteholonomy, discretedifferentialgeometry 1 Introduction Recent activity in the field of discrete differential geometry (DDG) has fueled the development of simple, robust, and efficient tools for geometry processing and physical simulation. The DDG approach to simulation begins with the laying out of a physical model that is discrete from the ground up; the primary directive in designing this model is a focus on the preservation of key geometric structures that characterize the actual (smooth) physical system [Grinspun 2006]. We addition- Efficient quasistatic treatment of material frame ally emphasize that the speed of sound in elastic rods is much faster for twisting waves than for bending waves. While this has long been established, to the best of our knowledge it has not been used to simulate general elastic rods. Since in most applications the slower waves are of interest, we treat the material frame quasistat- ically (§5). When we combine this assumption with our reduced coordinate representation, the resulting equations of motion (§7) become very straightforward to implement and efficient to execute. Geometry of discrete framed curves and their connections Because our derivation is based on the concepts of DDG, our dis- crete model retains very distinctly the geometric structure of the smooth setting—in particular, that of parallel transport and the forces induced by the variation of holonomy (§6). We introduceFigure 2: Helical perversion in experiment and simulation: starting from the natural shape of a Slinky R ? (top) we first remove writhe by pinching the flat cross-section between two fingers and traveling from one end to the other, arriving at a fully stretched, almost flat ribbon (middle). By bringing the ends together, a persistent perversion forms (bottom), whose shape is surprisingly different from the natural shape. simple algebraic tools and mnemonic diagrams that make it possi- ble to carry out in a methodical manner derivations involving the discrete connection induced by parallel transport. These tools are the central building blocks for deriving forces associated to dis- placements of the rod’s centerline and twist of the material frame. plifies both the analytical formulation and the numerical implemen- tation of the dynamics of symmetric rods. In graphics, Terzopoulos et al. [1987] introduced tensorial treat- ments of elastica, and Pai [2002] applied a discretization of the Cosserat rod model to simulate a strand. Bertails et al. [2006] used a piecewise helical discretization to produce compelling an- imations of curly hair using few elements per strand. Hadap [2006] considered a serial multi-body chain and used differential algebraic equations to treat the attendant numerical stiffness. These recent works used an implicit centerline representation based on reduced (curvature) coordinates. Inextensibility and encapsulated coupling with rigid-bodies When simulating inextensible rods, it becomes profitable to enforce the rod’s inextensibility via constraints rather than stiff penalty forces. Also, in general simulation applications, and in computer animation and gaming, it is often useful to enforce boundary con- ditions by coupling the rod to another physical system—most natu- rally, to a rigid-body, which carries exactly the degrees of freedom as any point on an elastic rod (i.e., position and orientation). Our approachin§8handlesbothtypesofconstraints, emphasizingphys- ical correctness as well as an important software reuse principle: the method of coupling allows our rod simulation to interact with any existing rigid-body simulation without internal modification of either the rod or rigid-body codes. By contrast, Choe et al. [2005] represented the centerline explicitly by a sequence of edges connected by linear and torsional springs. Gr´ egoire and Sch¨ omer [2007] proposed an explicit centerline dis- cretization coupled to a quaternionic material frame representa- tion. Spillmann and Teschner [2007; 2008] built on this idea in their development of algorithms for dynamic contact. Theetten et al. [2006] presented a geometric spline-based approach that can model large-displacement plastic deformations. These authors ar- gue that an explicit (displacement) representation of the centerline facilitates the simulation of complex contact scenarios and looping phenomena. We share this view. We show how to model and simulate inextensi- Validation ble elastic rods with arbitrary curved undeformed centerline and anisotropic bending response (§7), and we show the simplifica- tions that occur for naturally straight rods with isotropic bending response. We validate our model against known analytic solutions and present empirical evidence supporting the good convergence behavior of our discrete model to its smooth counterpart (§9). 3 Overview In many applications, the rod tends to bend or twist rather than to stretch; therefore, the case of inextensible rods is prevalent and im- 2 Related work Algorithm 1 Discrete elastic rod simulation Require: u0 Require: x0...xn+1 Require: (x0,˙ x0)...(xn+1,˙ xn+1) Require: boundary conditions 1: precompute ω ω ωj iusing (2) 2: set quasistatic material frame (§5.1) 3: while simulating do 4: apply torque to rigid-body (§8.2) 5: integrate rigid-body (external library) 6: compute forces on centerline (§7.1) 7: integrate centerline (§7.2) 8: enforce inextensibility and rigid-body coupling (§8) 9: collision detection and response // [Spillmann and Teschner 2008] 10: update Bishop frame (§4.2.2) 11: update quasistatic material frame (§5.1) 12: end while Mathematical analysis, modeling, and simulation of elastic rods is an active field in mechanics [van der Heijden et al. 2003; Goyal et al. 2007], numerical analysis [Falk and Xu 1995], and geometry [Bobenko and Suris 1999; with applications spanning biology [Yang et al. 1993; Klapper 1996], knots [Phillips et al. 2002; Brown et al. 2004], graphics [Chang et al. 2007]. It is impossible to survey the many works in this area in a brief space, so we discuss only the most closely-related works. For a broader starting point, refer to the expositions of Rubin [2000] and Maddocks [1984; 1994]. For the state of the art in strand and hair simulation, refer to the survey by Ward et al. [2007] and the course notes of Hadap et al. [2007]. // Bishop frame vector in frame {t0,u0,v0} at edge 0 // position of centerline in rest state Lin and Schwetlick 2004] medicine [Lenoir et al. 2006], the // initial position/velocity of centerline // free, clamped or body-coupled ends study of and computer // [Smith 2008] // [Hairer et al. 2006] In mechanical engineering, elastic rods are typically treated with a finite difference [Klapper 1996] or finite element [Yang et al. 1993; Goyal et al. 2007] discretization of the smooth equations. Gold- stein and Langer [1995] observed that using the Bishop frame sim-

portant. We visualize an inextensible rod as an infinitesimally thin, fixed centerline that can bend but not stretch and that is surrounded by a finite, but thin, elastic material. We consider the general case ofnaturallycurvedrodswithanisotropicbendingresponseandnote thesimplificationsthatoccurinthespecialcaseofnaturallystraight rods with isotropic bending response. We also present a simple and efficient method for attaching a rod to a rigid-body (§8). The dif- ferent components of our method are summarized in Algorithm 1. Barred quantities are precomputed and subsequently held fixed. t t γ(s) 4 Kirchhoff rods Figure 3: Adapted framed curve (Left) The configuration of an elastic rod is represented by a curve γ γ γ(s) and a material frame {t(s),m1(s),m2(s)}. (Right) The material frame is encoded by an angle of rotation θ relative to the natural Bishop frame {t(s),u(s),v(s)}. 4.1 Smooth setting We describe the configuration of Framed-curve representation a rod by an adapted framed curve Γ = {γ γ γ ;t,m1,m2} (see Fig. 3). Hereγ γ γ(s) is an arc length parameterized curve in R3describing the rod’s centerline; the assignment of an orthonormal material frame {t(s),m1(s),m2(s)} to each point on the centerline contains the requisite information for measuring twist. The material frame sat- isfies t(s) = γ γ γ′(s), i.e., it is adapted to the centerline such that the first material axis is tangent to the curve. We will refer to κ κ κ = t′as the centerline’s curvature (normal) vector. energy density. We do not assume that B is diagonal—in fact, we determine B by requiring that the undeformed material frame is a Bishop frame (see below). We also generalize to naturally curved rods by subtracting away the undeformed centerline curvature ω ω ω, i.e., using (ω ω ω −ω ω ω) in place of ω ω ω above. Putting all this together, we have The Kirchhoff theory of elastic rods assigns an Elastic energy elastic energy, E(Γ), to any adapted framed curve Γ. This energy is assembled from three scalar functions that measure strain—given by the change of the orthonormal frame {t(s),m1(s),m2(s)} ex- pressed in its own coordinates: Ebend(Γ) =1 Z (ω ω ω −ω ω ω)TB(ω ω ω −ω ω ω)ds , 2 where barred quantities refer to the undeformed configuration. The particular case of an isotropic, naturally straight rod is recovered by taking B = αId2×2and ω ω ω =0 0 0. ω1= t′·m1, ω2= t′·m2, m = m′ 1·m2. and Notice that since t′= κ κ κ, the first two of the above terms, ω1and ω2, represent the rod’s curvature vector expressed in material coor- dinates and measure the bending of material frame. The last term, m, refers to the twist of the material frame around the tangent. Ac- cordingly, total elastic energy contains bending and twisting contri- butions: 4.1.2 Twisting energy Letting m = m′ the centerline, the twisting energy is given by 1·m2denote the twist of the material frame about Etwist(Γ) =1 Z βm2ds . 2 E(Γ) = Ebend(Γ)+Etwist(Γ) . The formula m = m′ of material vectors immersed in ambient space. We now seek an equivalent expression that exposes a reduced set of coordinates. 1·m2gives an expression for the twist in terms The from this type of energy using Lagrangian mechanics (see, e.g., [Audoly and Pomeau 2008]). The goal of the present paper is to derive a discrete form of these equations. classical Kirchhoff equations for rods are obtained Given a Parallel transport and the Bishop (natural) frame fixed centerline, consider the task of assigning to it the geometri- cally most natural (and physically the most relaxed) frame. The Bishop frame {t(s),u(s),v(s)} for a given centerline is an adapted frame with zero twist uniformly, i.e., u′·v = −v′·u = 0. The as- signment of an adapted frame to one point on the curve uniquely fixes the Bishop frame throughout the curve. Our convention is to assign the Bishop frame at the first endpoint (s = 0) of the curve. Our assumption that s is an arc length parameterization implies that the rod is inextensible; therefore, we do not include a stretching en- ergy. Instead, we enforce inextensibility via an auxiliary constraint (§8). It is straightforward to drop this assumption by also including a stretching term. 4.1.1 Bending energy Consider traversing the centerline from one end to the other at unit speed. The evolution of the Bishop frame (and any other orthonor- mal frame) can be described in terms of its Darboux vector, Ω(s): When the rod’s undeformed configuration is straight (as opposed to curved) and the bending response is isotropic (as opposed to giv- ing preference to some bending directions over others), the bending energy takes the simple form t′= Ω×t , u′= Ω×u , v′= Ω×v . and Ebend(Γ) =1 αω ω ω2ds =1 Z Z ακ κ κ2ds , Since by the definition of the Bishop frame, u′· v = 0, it fol- lows from the second equation that Ω has no tangential component (along t). This, together with the first equation implies that Ω=κb, where κb = t×κ κ κ is the curvature binormal along the centerline. The Darboux vector of the Bishop frame serves to define paral- lel transport, a concept that plays a central role in the remainder of this paper. We parallel transport a vector x from one point on the centerline to another by integrating the equation x′= κb×x. 2 2 where the 2-vector ω ω ω = (ω1,ω2)Trepresents the centerline curva- ture vector expressed in the material frame coordinates and α is the rod’s bending modulus. We generalize this to anisotropic bending response by replacing the (isotropic) dot product with a general quadratic form B (a sym- metric positive definite 2×2 matrix), so that ω ω ωTBω ω ω is the bending

vertices would be ambiguous. Since tangent vectors belong to the triad that makes up frames on curves, we naturally assign frames to edges, rather than to vertices. In the smooth setting, we Pointwise vs. integrated quantities consider certain quantities (i.e., curvature and twist) at each point on the rod, and we define the energy by integrating over the length of the rod. By contrast, in the discrete setting, certain quantities often emerge naturally in an integrated rather than a pointwise way. For example, relating discrete curvature to the turning angle be- tween incident edges corresponds to an integration over the curve’s Gauss image [Grinspun 2006]. ei-1 ti-1 xi π−φi ti xi-1 ei xi+1 Figure 4: Notation Angles and vectors used in our discrete model. We refer to an integrated quantity as a measure that associates a number to a domain D of the curve, corresponding to an integral of a function over that domain. To convert an integrated quantity back to a pointwise one, we simply divide by the length, |D|, of the domain of integration. In the discrete case, we define Dias the Voronoi region associated to each vertex, having length |Di|=li/2, where li= |ei−1|+|ei|. Thus, infinitesimally, parallel transport corresponds to a rotation about the binormal—a concept that we will use again in our dis- crete model. Parallel transport keeps the tangential component of x tangential, and evolves the cross-sectional component of x via a tangential velocity, in particular without rotating the cross-section about the centerline. Observe that the three Bishop axes evolve un- der parallel transport. 4.2.1 Discrete bending energy The (twist-free) Bishop Curve-angle frame allows for a simple parameterization of the material frame [Langer and Singer 1996]. Let θ(s) be the scalar function that measures the rotation about the tangent of the material frame relative to the Bishop frame (see Fig. 3): representation Recall that the key player for formulating elastic bending energy in the smooth case was the curvature of the centerline. Since each edge is straight, it follows that Discrete curvature discrete curvature is naturally associated with vertices. Letting φi denote the turning angle between two consecutive edges (see Fig. 4) we define (integrated) curvature by m1= cosθ ·u+sinθ ·v , The key observation is that twist, m, can be expressed in terms of θ by m(s) = θ′(s). This relation is easily derived using the facts that m = (m1)′·m2and u′·v = 0. Hence, we write the twist energy as Etwist(Γ) =1 2 m2= −sinθ ·u+cosθ ·v . κi= 2tanφi 2. We will see in §6 that this definition of discrete curvature emerges naturally from an analysis of the geometry of discrete rods, and to be consistent throughout, we will use this definition in the bending energy. Note that κi→ ∞ as incident edges bend toward each other, so this measure of curvature will penalize sharp kinks in the rod. Z β(θ′)2ds . Observe that we have expressed the elastic energy of Kirchhoff rods by two dominant players: the position of the centerline, γ γ γ, and the angle of rotation, θ, between the Bishop and the material frame. This reduced coordinate representation is the cornerstone of our discrete theory. We define the curvature binormal Discrete curvature binormal at a vertex (an integrated quantity) by 2ei−1×ei |ei−1||ei|+ei−1·ei. (κb)i= (1) 4.2 Discrete Kirchhoff rods Observe that the vector (κb)i is orthogonal to the osculating plane passing through two consecutive edges and has magnitude 2tan(φi/2). In particular, this quantity is well-defined in the case of collinear edges, even though the binormal itself is not. When painting the discrete picture, our guiding principle is to seek a viewpoint that builds on the same geometric principles as the cor- responding smooth theory. We distinguish between primal quantities, as- Primal vs. dual sociated with vertices and decorated with lower indices, and dual quantities, associated with edges and decorated with upper indices. This diagram summarizes our indexing and notation conventions: We now have all the pieces to assem- Discrete bending energy ble the bending energy of a discrete naturally straight, isotropic rod: ?2li n n α(κbi)2 li ?κbi Ebend(Γ) =1 ∑ i=1 ∑ i=1 α 2= . 2 li/2 x1 x3 xn xn+1 x0 x2 The division by li/2 is due to converting the integrated quantity κbi into a pointwise one (see above), whereas the multiplication by li/2 takes care of the measure of the domain of integration, Di. e0 e1 e2 en A discrete framed curve, Γ, consists of Discrete framed curves a centerline comprised of (n+2) vertices x0,...,xn+1and (n+1) straight edges e0,...,ensuch that ei= xi+1−xi, together with an assignment of material frames Mi= {ti,mi Fig. 4). We consider adapted frames, where ti= ei/|ei| is the unit tangent vector per edge. Notice that tangent vectors of polygonal curves are uniquely defined at edges, whereas their definition at Recall that for an anisotropic bending response, we require the ma- terial curvatures given by inner products between curvature vectors and material frame vectors. Using the curvature binormal (κb)ias defined by (1), we express the material curvatures as 1,mi 2} per edge (see ?T ? ω ω ωj (κb)i·mj 2,−(κb)i·mj i= j ∈ {i−1,i} . for (2) 1

5 Quasistatic material frame postulation Here, following our primal-dual notation, the double indices— upper and lower—correspond to an edge-vertex pair. Barred quan- tities we denote evaluation on the undeformed configuration. Our goal is to arrive at equations describing the time-evolution of the centerline of the rod. We first pave the way with an important simplifying assumption. With these definitions, the bending energy of a discrete rod is n i ?TBj? 1 ? ? ω ω ωj i−ω ω ωj ω ω ωj i−ω ω ωj The speed of a twist wave increases as the rotational inertia of the cross section vanishes. By contrast, bending waves are disper- sive [Sommerfeld 1964]—their speed depends on their wavelength, λ—with velocity vh/λ, where h is the rod radius, and v is solid ma- terial’s speed of sound. Consequently, twist waves travel faster than bending waves, for bending waves with large wavelengths λ ≫ h, i.e., much larger than the rod radius. For the problems we consider, the relevant dynamics are in (large-wavelength) bending modes, and it is wasteful of computation to resolve the temporal evolution of the twist waves. ∑ i=1 ∑ j=i−1 Ebend(Γ) = . (3) i i 2li This formulation allows for a general (anisotropic) bending re- sponse and a general (curved) undeformed configuration. As in the smooth case, the simple expression of the special case is recovered by taking Bj= αId2×2and ω ω ωj i=0 0 0. 4.2.2 Discrete twisting energy To formulate the discrete twisting energy, we follow the same pro- gression as in the smooth setting, first laying out the notions of parallel transport and the Bishop frame, and then introducing the angle of rotation to the material frame. Consider the limit of a vanishing cross-sectional rotational inertia: twist waves propagate instantly. At any instant in time, the material frames are the minimizer of elastic energy, subject to any given boundary conditions on the material frame: We define Discrete parallel transport and Bishop frame discrete parallel transport in analogy to the smooth case, i.e., as a rotation Piabout the curvature binormal, ∂E(Γ) ∂θj = 0 , (4) Note that the quasistatic treatment of the material frame eliminates the θjvariables as degrees of freedom from the system: given the rod’s centerline and the boundary conditions on the material frame, (4) is used to determine what the material frame must be. ? ti−1? ? ti−1×ti? = ti−1×ti. = ti, Pi Pi By convention, Piis the identity if ti−1= ti, whereas Piis not defined if ti−1= −ti. Par- allel transport is a key notion that allows us to generalize the notion of twist from the smooth to the discrete setting. The value of θjfor an edge may be pre- Boundary conditions scribed by a given boundary condition on the material frame. In practice, we consider stressfree ends (i.e., no boundary conditions) or clamped ends (i.e., assigning the material frame for j = 0 and j = n). We use (4) for all θjvariables whose value is not pre- scribed by a boundary condition, ensuring that the number of equa- tions matches the number of unknowns for the angles θj. In order to define Bishop frames, we draw upon the smooth case by transporting a unit vector u0, which is orthogonal to t0, along the rod using our rotation operators Pi, i.e., we iteratively define ? ui−1? ui= Pi . 5.1 Quasistatic update The third axis of each Bishop frame is then vi= ti×ui. During simulation, we enforce the quasistatic nature of the material frame by ensuring that (4) is satisfied before forces are computed. The Bishop frame is by definition adapted Bishop frame update to the centerline. This requires that u0⊥ t0, a condition that must be maintained during simulation. Each simulation step moves the centerline to a new position, so that in general after Algorithm step 9 the requirement u0⊥ t0may be violated (the exception is when t0is clamped). We re-adapt the frame by parallel transporting the Bishop vector (in time)—in analogy to the parallel transport along the centerline—by computing a rotation operation. For an Special case of naturally straight, isotropic rods isotropicrodwithanaturallystraightundeformedconfiguration,(4) implies that a clamped rod has uniform twist, =θn−θ0 mi li = constant , (5) 2L Since both the material frame Material frame representation and Bishop frame are defined at edges, let θidenote the angle needed to rotate the Bishop frame into the material frame at edge ei (see Fig. 4). Then the material frame vectors at edge i are where 2L = ∑n depends only on the boundary conditions: the difference of angles of the end edges, θn−θ0. Substituting the above into the formula for E(Γ) gives the simplified expression i=1li, and the pointwise twist mi/li(recall §4.2.2) 1= cosθi·ui+sinθi·vi, and mi 2= −sinθi·ui+cosθi·vi. mi +β?θn−θ0?2 n α(κb)2 i ∑ i=1 E(Γ) = . (6) In perfect analogy to the curve-angle Discrete twisting energy representation in the smooth case, we adopt the twisting energy 2L li Observe that the twist energy depends only on the angle between material and Bishop frame at the boundary edges of the rod. For the naturally straight, isotropic case, the above implies that a rod with stressfree ends has no twist. β(θi−θi−1)2 n n βm2 i ∑ i=1 ∑ i=1 Etwist(Γ) = = , li li where the scalar mi= θi−θi−1is the (integrated) discrete twist. The variable mimeasures the angle between two adapted frames: the result of parallel transporting the material frame from edge ei−1 to eiand the material frame at edge eiitself. For the general case, we use Newton’s method General case to minimize the elastic energy E(Γ) with respect to the material

6 Discrete holonomy Holonomy is a classical concept from differential geometry: it mea- sures the deficit of geometric data (frames) to close up when paral- lel transported around a closed loop. In our case (of discrete rods), holonomy depends—just like parallel transport itself—only on the centerline; it measures the (scalar) angle when parallel transporting an adapted frame along a closed loop of discrete edges. Indeed, the fact that holonomy can be expressed as a scalar is the reason why it is such a useful concept for computing forces. Consider the variation of two consecutive edges ei−1(ε) and ei(ε) with tangents ti−1(ε) and ti(ε), and parallel transport Pi(ε), where ε denotes the variation parameter. We denote by˜Pi(ε) the parallel transport from ti(0) to the deformed configuration, ti(ε), i.e., the rotation that satisfies˜Pi(0) = Id and Figure 5: Asymmetry of twist: anisotropy of the cross-sections and natural curvature of the centerline “conspire” to produce non- uniform twist distribution. From top to bottom: reference configu- ration, moderate twist, large twist. ? ? ? ? ˜Pi(ε) ti(0) = ti(ε) ˜Pi(ε) ti(0)×ti(ε) = ti(0)×ti(ε) . and Now consider the concatenation of parallel transports given by Ri−1(ε) =?˜Pi−1(ε)?T◦?Pi(ε)?T◦˜Pi(ε)◦Pi(0) , whichwedepictbythefollowingmnemonicdiagram, wherearrows represent parallel transport: (8) frames, which requires both the gradient and Hessian of the energy with respect to the θjvariables. The gradient is given by ! -ti(ε) ∂E(Γ) ∂θj ∂ −mj+1 lj+1 mj lj ti−1(ε) ?Wj+Wj+1 where ?+2β (ω ω ωj = , (7) ∂θj 6 6 ∂ ∂θjWi=1 ψi(ε) i)TJBj(ω ω ωj i−ω ω ωj i) . li -ti(0) ti−1(0) To derive the latter identity, note that ∂mj ∂mj two dimensional vectors by counter-clockwise π/2 rotation. As ex- pected, (5) is a special case of (7). 1/∂θj= mj i, where J acts on 2and 2/∂θj= −mj 1to obtain ∂ω ω ωj i/∂θj= −Jω ω ωj Observe that Ri−1(ε) corresponds to traversing this diagram in counter-clockwise order, starting at ti−1(0). It follows from the def- initions that Ri−1(ε) maps the tangent ti−1(0) to itself and is there- fore a rotation by angle ψi(ε) about axis ti−1(0). In the language of differential geometry, ψi(ε) is the holonomy of the connection induced by parallel transport around the depicted closed loop. The Hessian, H, is obtained by differentiating (7) with respect to θj−1, θj, and θj+1. The resulting components of the Hessian are: The ingredient that we require is the gradient of ψifor 1 ≤ i ≤ n. Building on the literature of a related quantity known as the writhe of polygonal curves [de Vries 2005], we obtain the variation of ψi with respect to a centerline displacement δx: Hj,j−1= −2β Hj,j+1= −2β , , lj lj+1 ∂2 Hj,j=2β +2β lj+1 ?Wj+Wj+1 i)−1 li ?, + (∂θj)2 lj ?1 ? δxi−δxi−1 |ei−1| δxi+1−δxi |ei| δψi=−2ti−1×ti 1+ti−1·ti It can be shown that the corresponding expression in the smooth setting is δ(Rψds)=−Rκb·(δx)sds where the subscript s denotes differentiation with respect to s. Compare the discrete expression to the integrand of the smooth case; the second factor of the discrete expression is a finite-difference approximation of (δx)s. In analogy to the smooth form, we take the first factor to be the integrated curvature binormal. Note that this definition coincides with (1), and allows us to rewrite the discrete result as +1 · . ∂2 (∂θj)2Wi=1 i)TJTBjJ(ω ω ωj i)TBj(ω ω ωj (ω ω ωj (ω ω ωj i−ω ω ωj 2 2 i) . where li Note that the Hessian is tridiagonal and so it can be factored in O(n) time. Additionally, in order to accelerate the convergence of New- ton’s method, we combine it with a line search along the Newton direction [Press et al. 2002]. Finally, we note that Newton’s method converges quickly: the θjvariables do not change significantly be- tween time steps, so that previous values are a good initial guess. (κb)i 2|ei−1|, ∇i+1ψi= −(κb)i We presented the bending and twisting energy of a Forward smooth Kirchhoff rod and constructed by direct analogy the corre- spondingenergiesforadiscreteKirchhoffrod. Thechallengeahead is to derive the elastic forces induced by the gradients of these ener- gies. Theseforceswillincludetermsthatarisefromthedependence of the vectors ujand vjon centerline positions, xi. Since these Bishop vectors are defined via parallel transport, we must consider the variation of parallel transport with respect to moving the center- line. For this, we turn to the concept of holonomy. ∇i−1ψi= , (9) 2|ei| and ∇iψi= −(∇i−1+∇i+1)ψi. 6.1 Variation of parallel transport With our new tool in hand, we can derive the change in the Bishop frame when varying the rod’s centerline. Since the Bishop frame is

by requirement adapted to the curve, we are interested in the angle that the frame is rotated by about the tangent. This situation is again depicted in the following diagram: boundaryconditions, thesumiszero. Forclampedboundaries, only one component of this sum could be non-zero, allowing us to write the force as ∂θn ∂xi −∂E ∂xi −∂E ∂θn 0-F1(ε) 0-··· Fj−1(ε) 0-Fj(ε) . F0(ε) 6 6 6 6 To evaluate ∂θn/∂xi, recall from §6.1 that varying the centerline rotates the Bishop frame at edge enby angle Ψnabout the tangent tn. Therefore, to keep the material frame aligned to the boundary condition, we must subtract this angle from θnto obtain Ψj ψ1(ε) ψ2(ε) ψj(ε) 0 0-F1(0) -··· Fj−1(0) 0-Fj(0) F0(0) 0 n ∂ψj ∂xi −∂E ∂xi +∂E ∂θn ∑ j=1 . Here, the Bishop frames Fi= {ti,ui,vi} are displayed in place of the tangents to show that we are parallel-transporting these frames from one edge to the next. Additionally, each labeled arrow indi- cates the angle needed to align the result of parallel transporting the frame at the tail of the arrow with the frame at the head. Thus, the horizontal arrows are labeled with 0 to indicate that no twist is re- quired to align Pi(Fi−1) to Fi; the first vertical arrow is also labeled with 0 because we ensure that F0is always updated via parallel transport (see §4.2.2). We are interested in the angle of rotation, Ψj, required to align˜Pj(ε)?Fj(0)?to Fj(ε). Since parallel transport commutes with twist and holonomy is ad- ditive under concatenation of loops, it follows from traversing this diagram in counter-clockwise order that the resulting angle is Ψj=∑j i=1ψi. For computing forces, we require the gradient of this angle with respect to vertex positions, which is given by For clamped boundaries, we give the components required to eval- uate this force for both the special and general case. For naturally straight, isotropic rods, the forces on Special case vertex xiare given by up to three contributions: ?T(κb)j+β(θn−θ0) −2α ?∇i(κb)j ∇iψj i−1 ≤ j ≤ i+1 , L lj with the gradient of holonomy given by (9) and the gradient of the curvature binormal given by ∇i−1(κb)i=2[ei]+(κbi)(ei)T ∇i+1(κb)i=2[ei−1]−(κbi)(ei−1)T |ei−1||ei|+ei−1·ei ∇i(κb)i= −(∇i−1+∇i+1)(κb)i. Here [e] is a skew-symmetric 3×3 matrix acting on 3-vectors x by [e]·x = e×x. In deriving these gradients we have used the assumption of inextensibility of the rod’s edges. |ei−1||ei|+ei−1·ei, j , ∇iΨj= ∑ k=1 ∇iψk. (10) By (9), this sum will have at most three non-zero terms. In hindsight, we view the relation of holonomy to Discussion the profile of the centerline as an extension of the celebrated C˘ alug˘ areanu-White-Fuller theorem [1971] for discrete curves. Fuller [1978] noted that this theorem is applicable to the equilibria of closed elastic rods with isotropic cross-sections. Our develop- ment extends to the general case of open boundaries and even to anisotropic cross-sections. For anisotropic rods with a curved rest shape, the General case required expressions for the forces are given by: ∂E ∂θn=1 n)+2βmn n)TJBn(ω ω ωn n−ω ω ωn (ω ω ωn and ln ln n k ∂E ∂xi ?TBj? 1 7 Equations of motion ? ? ∇iω ω ωj ω ω ωj k−ω ω ωj ∑ k=1 ∑ = , k k lk j=k−1 We are now ready to write the equations governing the time- evolution of the rod’s centerline. Since the material frames depend on the rod’s centerline and are not independent degrees of freedom (see §5), we must consider this dependence as well when comput- ing the centerline forces. where the gradient of the material-frame curvature ω ω ωj kis given by ! (mj −(mj 2)T 1)T ∇iω ω ωj ∇i(κb)k−Jω ω ωj k(∇iΨj)T. k= (11) 7.1 Forces on centerline This last result is readily apparent from the definition of material- frame curvature in (2) and the variation of the Bishop frame in §6.1. In deriving the forces on the rod’s centerline, we must consider how the energy depends on the xivariables—both directly and indirectly by considering the effect that moving a vertex has on the rod’s ma- terial frame. The force acting on vertex xiis therefore given by 7.2 Integrating the centerline Since the material frame is always updated to be in quasistatic equi- librium (see §5.1), we only need to update the centerline based on the forces derived above. The equations of motion are n ∂θj ∂xi −dE(Γ) dxi = −∂E(Γ) ∂E(Γ) ∂θj ∑ j=0 − . ∂xi M¨ x = −dE(Γ) , Here, the total derivative dE/dxitakes into account the implicit dependence of potential energy on centerline positions, whereas the partial derivative only takes into account explicit dependence. dx where M is a 3(n+2)×3(n+2) (diagonal) mass matrix associ- ated to centerline positions. We discretize this equation using the symplectic Euler method [Hairer et al. 2006]. We use a manifold- projection method to enforce the inextensibility constraint. Recall from (4) that ∂E/∂θjvanishes for all edges for which θj is not prescribed by a boundary condition. Therefore, for stressfree

8 Constraints both positional as well as rotational degrees of freedom by consid- ering the metric induced by the kinetic energy1 (3n+12)×(3n+12) generalized mass matrix and the generalized velocity y ∈ R3n+12are defined respectively by Here M is the rod’s (diagonal) 3(n+2)×3(n+2) mass matrix, M is the body’s total (scalar) mass, and I is the body’s (symmetric) momentofinertiatensorexpressedasa3×3matrixinthereference coordinates. Since q is a unit quaternion, q−1˙ q corresponds to a vector in R3[Hanson 2006]. The kinetic energy has contributions from the body rotation, body translation, and centerline translation. 2y˜ MyT, where the For simulating extensible rods, it would suffice to include a stretch- ing component to the energy; however, simulating inextensible rods using stretch forces would lead to unnecessarily stiff equations. In addition, many interesting physical systems can be modeled by the coupling of elastic rods to rigid-bodies. Canonical examples include the torsional and Wilberforce pendulums, comprised of a rigid-body suspended under gravity by a thin elastic rod. To avoid numerical stiffness associated with maintaining inextensibility and rigid-body coupling, we add auxiliary constraints to our system. 4·I y = (q−1˙ q,˙ r,˙ x) . ˜ M = M ·Id3×3 and M For each edge of the rod, we use Inextensibility constraints the quadratic constraint equations ei·ei−ei·ei= 0, where (the pre- computed) eirefers to the undeformed configuration. We initialize the first iteration of fast projection with the results of an unconstrained time step. Each step of fast projection improves on the guess by computing the first Newton iteration for the min- imization of the functional y˜ MyT−CTλ λ λ, where (using Golden- thal’s notation) C is the vector of constraint equations, and λ λ λ is the vector of corresponding Lagrange multipliers. Thus, by for- mulating the kinetic energy using a generalized mass matrix in a high-dimensional configuration space, we are able to directly apply fast projection iterations to rigid-body coupling. In all of our exam- ples, we found that after 3–5 steps the method converged to within a tolerance of 10−8in satisfying the constraint C = 0. The welding of the body to Rigid-body coupling constraints the rod requires the body’s position and orientation, and the rod edge’s position and material frame, to be in one-to-one correspon- dence. Attaching the rigid-body to the first edge of the rod gives three constraints: qx0q∗+r−x0= 0 , qx1q∗+r−x1= 0 , q·q−1 = 0 , where r ∈ R3is the translation vector and q ∈ H is the unit quater- nion [Hanson 2006] mapping the body’s center of mass and orien- tation, respectively, from the reference to the current configuration, and q∗is the conjugate of q. The first equation ensures that q is unit length, so qx0q∗is a rotation of x0that is summed with r us- ing vector addition. The second and third equations ensure that the rod’s first edge and the rigid-body transform identically. After the iterations converge, fast projection requires a velocity up- date [Goldenthal et al. 2007]. Using our generalized coordinates, the appropriate update is (˙ q,˙ r,˙ x) ← (˙ q,˙ r,˙ x)−1 h(2q0q−1,r0−r,x0−x) . 8.1 Enforcing constraints Numerous approaches are available in the literature for maintaining constraints acting on a mechanical system [Shabana 2001]. Perhaps the most simple and straightforward of these is to use the penalty method; alas, as mentioned above, the penalty forces required to closely enforce the constraints are unacceptably stiff and require small time steps to ensure stability [Goldenthal et al. 2007]. Our discussion above immediately generalizes to the case of mul- tiple rigid-bodies attached to multiple rods. As a corollary, we can couple a rigid-body to any edge on the rod simply by splitting the rod at that edge. However, care must be taken to understand the induced boundary conditions: each rigid-body’s position and ori- entation serves to clamp the position and orientation of the center- line and material frame for each coupled rod end-edge. After fast projection, the material frame vector induced by the rigid body ori- entation is m0 1= qm0 An alternative is to employ an augmented Lagrangian formulation. Such formulations in general introduce additional variables (i.e., Lagrange multipliers) to enforce the constraints during the time in- tegration step of the algorithm. The exception are manifold pro- jection methods [Hairer et al. 2006], which perform constraint en- forcement as a post-integration step. We choose to maintain the constraints using this approach. 1q∗. 8.2 Torque transfer We consider a methodical categorization of those interface forces that are automatically transferred between the rod and rigid-body via projection and those that remain to be transferred explicitly. Observe that forces that correspond to gradients of the independent variables, q,r,x, are automatically transferred between the rod and the body by the projection step. The material frame is not an inde- pendent variable—it is a function of the centerline position and of the boundary conditions (see §5). Therefore, in Algorithm step 4, we explicitly transfer the torque acting on the material frame m0 the coupled rigid-body. Note that the rigid-body’s relatively large moment of inertia ensures that the explicit coupling is non-stiff. A manifold projection method integrates a mechanical system by alternating between an unconstrained time integration step and a constraint enforcement step. For the unconstrained step, we use an explicit symplectic Euler integrator (Algorithm step 7), and we call the ODE [Smith 2008] library to time step the rigid-body (Algo- rithm step 5). 1, to Manifold projection allows us to reuse an existing rigid-body code without modifying its time integration, collision response, or other internal structures. However, any other method for enforcing these constraints could be used, and we include a discussion of our ap- proach for completeness. The force acting on the material frame corresponds to the gradi- ent of energy with respect to the dependent angles, θi. By the quasistatic material frame assumption, the torque, τ = |τ|t0, ex- erted by the rod on the body is equal and opposite to the torque exerted by the body on the rod. To derive the magnitude of the torque we turn to the principle of virtual work [Lanczos 1986]. Holding the centerline fixed, consider a virtual displacement, δθ0, which varies the material frame about e0. Note that δθ0also af- fects the body’s orientation, due to the rod-body constraint. The For the enforcement step, we adopt Fast manifold projection and extend the fast projection method of Goldenthal et al. [2007]. This method takes an unconstrained configuration and finds a nearby constrained configuration. The notion of “nearby” is made precise by the natural metric on the configuration manifold. We ex- tend the application of fast projection to coupled systems involving

throughtheclamps, ϕ(s)=cos−1(t·ex)theangulardeviation of the tangent away from this axis, and ϕ0=maxsϕ(s) the maximal devi- ation at the center of the pattern. The envelope of the helix is given by f(ϕ(s))=tanh2?s given by s/s∗= 2α 1+cosϕ0 cosϕ−cosϕ0 1−cosϕ0 a fixed loading geometry, we obtain convergence towards this ana- lytical solution under refinement (see Fig. 7). ?, where s/s∗is the dimensionless arc length s, which simplifies to f(ϕ) = s∗ q1−cosϕ0 ?βm ? (see e.g. [van der Heijden and Thompson 2000]). With Figure 6: Torque transducer: two rigid-bodies coupled by a rod. The total energy, broken up into contributions from the bodies and the rod, is plotted. (rad) work performed by the body on the rod equates to the change in the rod’s stored energy: |τ|δθ0= (dE/dθ0)δθ0. Since this holds for any δθ0, it follows that |τ| = dE/dθ0. In the isotropic case, this yields |τ| = β(θn−θ0)/L. In the anisotropic case, we expand the full derivative as dE/dθ0= ∑i(∂E/∂θi)(∂θi/∂θ0). By the qua- sistatic assumption, most terms in this summation are zero, leading to |τ| = ∂E/∂θ0, which can be evaluated using (7). 25 20 15 10 5 9 Validation 0 0.4 0.6 0.8 1.0 1.2 The interaction between twisting and bending modes produces sur- prising instability phenomena, such as spontaneous buckling of rods. These instabilities are important for a variety of applications, such as for modeling DNA. When applied to some of these phe- nomena, our model not only reproduces the correct behavior quali- tatively, but in fact shows good convergence to analytical solutions (for those models where an analytical treatment is available in the literature). In the following we describe several example problems; refer to Table 1 for corresponding performance measurements. Figure 8: Michell’s buckling instability of an elastic ring with imposed internal twist θn. Left: above a critical value of θn, the planar, circular shape loses stability and buckles to a non-planar shape. Right: domain of stability of the circular shape with radius R = 1: simulations (dots) compared to theoretical threshold (black curve). Each dot corresponds to a simulation run with particular values of β and θn(α = 1 and n = 50 are fixed), initialized with a slightly perturbed circular shape; dots are colored in light blue when the amplitude of the perturbation decreases in time (stable) and in purple when it increases (unstable). 9.1 Convergence Another famous example of a rod becoming Michell’s instability unstable when subjected to twist is Michell’s instability, also known as Zajac’s instability (see Fig. 8). When the ends of a naturally straight, isotropic rod are closed into a ring in a twist-free manner, the resulting shape is circular. When these ends are twisted by an angle θnbefore they are closed up, and when this angle is pro- gressively increased starting from zero, the ring suddenly starts to writhe(buckle)outofplane, atawell-definedcriticalvalueoftwist. This effect is a perfect illustration of the coupling of the twist with the shape of the centerline. This critical twist can be computed an- alytically [Goriely 2006] as θn we study the stability of the circular shape numerically for different values of the rod parameters and compare to the analytical predic- tion. Our model displays excellent agreement with the predictions of the smooth theory. f(ϕ) f(ϕ) 1 1 analytical n = 180 n = 140 n = 110 n = 80 n = 60 n = 40 0.5 0.5 s/s* s/s* -3 3 -3 3 uniform non-uniform Figure 7: Localized buckling of a naturally straight, isotropic rod with an imposed angle of twist. The rod buckles into a helix with modulated amplitude (inset). Convergence of the envelope of this helix toward the analytical solution in the smooth case is observed for both a uniformly sampled rod (left) and a rod where one half is twice as refined as the other (right). The parameters for the simu- lation are: rod length L = 9.29, bending modulus α = 1.345, twist modulus β = .789, clamp rotated by 27 turns, imposed axial short- ening of 0.3 units, corresponding to a theoretical maximal deviation ϕ0= 0.919. c= 2π√3/(β/α). In Figure 8–right, 9.2 Special case Although knots have been extensively stud- Instability of knots ied as mathematical objects, the understanding of their mechanical behavior is far less advanced. Very recently, an analytical solution hasbeenobtaineddescribingtheequilibriumshapeofaloosetrefoil knot tied on a naturally straight, isotropic rod [Audoly et al. 2007] in the case where no twist is applied. We ran a simulation of a knotted rod with both ends clamped and were able to reproduce the typical equilibrium shape of the knot, which is made of a large, al- most circular loop connected to two flat tails by a braid region (see Fig. 1–left). We next twisted the ends of the rod and found that the shape of the knot first changes smoothly and then jumps at a critical value of the twist. This jump happens both for positive and negative applied twist, although the final shapes of the knot are qualitatively When a naturally straight, isotropic Localized helical buckling rod is clamped at its two ends with one end rotated by a given an- gle, the rod buckles as the two ends are quasi-statically translated towards one another. In the smooth case, the buckling of a long rod under twist is described by an exact solution of the equations of equilibrium. This analytical solution describes nonlinear (finite amplitude)bucklingawayfromthestraightconfigurationandmixes bending and twisting effects. We studied convergence of the equi- libria of our discrete model towards this solution in the limit of refinement. The geometry is shown in the inset of Figure 7. We de- notebysthearclength, extheunitvectorparalleltotheaxispassing

different (see Fig. 1). This fascinating behavior will be studied in detail in a future paper [Audoly et al. 2008]—it can be easily repro- duced using a small tube of an elastic material (a silicone wire in our experiments), and we encourage the reader to try it! occurs with the semi-circular shape, but it is less marked. In either case, the twist concentrates near the end that is rotated, although the rod properties are uniform along its length; this symmetry-breaking points to the fact that that the equations governing the quasi-static twist are nonlinear. Perversion is a classical pattern that can be Helical perversion: observed in naturally curved rods. It consists of a junction between two helices with opposite chiralities. It has been described by Dar- win in the tendrils of climbing plants and can be often observed in tangled phone cords [Goriely and Tabor 1998]. Perversions can be produced by first flattening a helical spring such as a Slinky R ? into a flat, straight ribbon by pulling on both of its ends (see Fig. 2– middle); next, by progressively releasing its ends from this straight configuration, one obtains a shape made of two mirror-symmetric helices (see Fig. 2–bottom). The final shape is surprisingly dif- ferent from the natural one as the chirality has been reversed in a half of the spring, as revealed by comparison with Figure 2–top. The existence of perversions is due to the nonlinear behavior of the rod—perversions belong to a general class of solutions that can be derived in nonlinear analysis by connecting two competing equi- librium configurations, which are here the right-handed and left- handed helices in the presence of natural curvature. Figure 9: Plectoneme formation: When the ends of a hanging elastic rod are twisted, it takes on structures known as plectonemes. The formation of plectonemes is governed by physical parameters, such as the twist rate, viscosity of the ambient fluid, and gravity. A rod generally forms entangled structures called Plectonemes plectonemes when subjected to large twist. Typically, the pattern formed by a twist-driven instability such as Michell’s instability grows from a weakly helical shape at short times, to fully developed plectonemes at long time. Such structures have been well-studied, for instance in the context of DNA supercoiling [Yang et al. 1993]. In this experiment, we started with a naturally straight rod and de- formed it into a parabola in a twist-free manner. Fixing the posi- tions of the endpoints of the rod and progressively increasing twist, we find that the rod starts to writhe out of the plane. Letting the instability develop fully, we observe plectonemes (see Fig. 9). De- pending on the viscous drag and gravity, we found that one or many plectonemes can be obtained. Single plectonemes have been de- scribed both analytically [van der Heijden et al. 2003] and numeri- cally [Goyal et al. 2007] in several papers; we observed an interest- ing phenomenon of plectonemes merging at long times, which has seemingly not yet been studied. Figure 10: Hanging chain: A chain consisting of curved elastic rods as links hangs under the influence of gravity. The material parameters for each rod are B = 2Id2×2and β = 2, with each link in the chain having a radius of approximately 1 unit. In Figure 10, we show a chain Hanging chain: free boundaries hanging from its two ends under the influence of gravity. Each link is an elastic rod with curved undeformed configuration with stressfree boundary conditions on the material frames. As shown in the accompanying animation, the two ends are pulled apart until the chain breaks and the links fly apart. Even though we are not applying a twist to any of the links in the chain, we still need the material frame to represent the undeformed curvature of the rod— recall that ω ω ω is defined as a 2-vector in material-frame coordinates. The simulation of plectonemes requires the treatment of rod self- contact, which is outside the scope of this paper; for the state of the art, refer to the treatment by Spillmann and Teschner [2008]. 9.3 General case 9.4 Rigid-body coupling The coupling of twisting and bending modes of naturally curved or anisotropic rods is fairly more complex than in the isotropic case, as demonstrated by in following experiments. When coupling rods with rigid-bodies, the transfer of energy be- tween these systems results in complex motions—in particular, through the interplay between bending and (non-uniform) twisting of the rod and the translational and rotational moment of the at- tached mass. In order to validate our model, we have in particular tested its ability to faithfully reproduce real-world experiments. The combination of anisotropy and non- Asymmetry of twist zero rest curvature can result in a non-uniform distribution of twist along the rod. In Figure 5, we show this effect for a rod with anisotropic cross section, whose centerline is either V-shaped or semi-circular in the reference configuration. We clamp one end of each of these rods and twist the other. In the V-shaped case, we observe that twist is first confined on one half of the rod—it takes a significant amount of twist before twist manages to “jump over” the kink in the middle of the rod. A similar phenomenon This experiment reveals fascinating as- Wilberforce pendulum pects of how bending and twisting interact and thereby affect the motion of an attached rigid-body. Consider a helicoidal spring (an isotropic rod whose reference state is a helix) whose upper end is

fig. no. verts. time-step forces contact quas. update fast proj. total 1 60 4.0 0.026 0.05 0.021 0.21 0.31 5 49 1.0 0.025 0.018 0.021 0.12 0.18 6 200 2.0 0.072 0.09 0.064 0.81 1.0 9 275 1.0 0.13 0.19 0.17 0.42 0.92 7 75 2.0 0.043 0.1 0.2 0.34 8 67 50.0 0.039 0.096 0.28 0.42 10 31 1.0 0.016 0.13 0.13 0.2 11 20 1.0 0.0036 0.0043 0.019 0.027 12 2158 0.1 0.87 0.64 21.0 22.0 Table 1: Performance evaluation: This table summarizes timing information (in milliseconds per simulation step) for examples de- picted in the figures, as measured on a single-threaded application running on a 2.66GHz Core 2 Duo. For examples run without col- lision detection, collision timings are omitted. Figure 11: Wilberforce pendulum: Due to a weak coupling be- tween the bending and twisting modes of a stretched spring, the en- ergy of the system is transferred back and forth between the trans- lational (red) and angular (blue) modes of oscillations. simulating tree-like one-dimensional configurations with several T- junctions. In Figure 12, we show an example of such coupling, where we additionally used external forces (wind and gravity) to increase the dramatic effect. Here the mass of the rigid-bodies (but not necessarily their spatial extension) can be tuned to achieve a variety of dramatic effects. This example shows the power of cou- pling rods with rigid-bodies, indicating the attractiveness of this approach to be used in a wide range of graphics applications, such as for simulating plant motion, or the skeletal dynamics of rigged characters. fixed, and attach a rigid-body to its lower end. The weight of the body stretches the spring when the whole system is at rest. Mov- ing the mass slightly upwards from this equilibrium and releas- ing it, first leads to vertical oscillations. Progressively, the system switches to twist oscillations and the vertical ones are extinguished; later, the twist oscillations start to decrease and the vertical ones reappear, and so on (see Fig. 11). The nonlinear behavior of the spring, which is captured accurately by our model, is responsible for this energy transfer between the two modes [Sommerfeld 1964]: stretching a spring affects its eigenmodes. Note the presence of sta- tionary waves in the spring, which are also clearly visible in our simulation movie. 10 Conclusion Rods can be used to couple rigid-bodies, Torque transducer with the rod effectively acting as a “torque transducer” between them. In Figure 6, we plot the kinetic and potential energy of the two rigid-bodies and the rod (one of which has much higher mass than the other), observing the energy transfer among the stored po- tential energy of the rod and the rotational energy of the two rigid bodies. Our use of the Fast Projection Limitations and future work algorithm leads to energy being dissipated during the constraint- enforcement step. In many applications, the energy behavior of the system is of interest, so we would like to explore alternate meth- ods for enforcing constraints that are both efficient and have good energy behavior. The formulation of discrete rods allows us to impose boundary con- ditions on the material frame at any edge along the rod. We have presented boundary conditions that arise due to explicitly clamping certain edges or coupling them to rigid-bodies. We are currently considering the effect of frictional contact on the material frame and would like to have an implementation that also allows contact constraints to dictate boundary conditions for the material frame. We are also interested in employing adaptivity in our current model and believe that the ideas of Spillmann and Teschner [2008] pave the way for this. In addition, an interesting related topic would be in providing higher-order methods for simulating elastic rods. The main issues involved would be in defining convergent discrete operators, such as parallel transport and curvature, on higher-order elements, and in ensuring the resulting order of accuracy. We would like to thank Michael Herrmann Acknowledgments for his indispensable insights on unifying the isotropic and anisotropicpicture, S´ ebastienNeukirchforhishelpinreviewingthe literature on elastic rods, Jonas Spillmann and Matthias Teschner for sharing their code for CORDE with us, David Harmon for his help with collision detection and response, Kori Valz for her artistic contributions, and Aner Ben-Artzi, Mathieu Desbrun, and Jerrold E. Marsden for early discussions regarding this project. This work was supported in part by the NSF (MSPA Award No. IIS-05-28402, CSR Award No. CNS-06-14770, CAREER Award No. CCF-06- 43268) and the DFG Research Center MATHEON in Berlin, as well as by Autodesk, Elsevier, mental images, NVIDIA, and Walt Dis- ney Animation Studios. Figure 12: Tree in wind: rigid-bodies connecting multiple rods together at a single point (reference configuration shown in inlay). Using this method, we can simulate a tree bending and twisting in response to a strong arctic wind. Note that without resistance to twist, the tree would start spinning due to vortices in the wind. Coupling rods with rigid-bodies Rigid-bodies for rigid couples opens the possibility for complex simulations—in particular, for

References HAIRER, E., LUBICH, C., AND WANNER, G. 2006. Geometric Numerical Integration: Structure-Preserving Algorithms for Or- dinary Differential Equations. Series in Comp. Math. Springer. AUDOLY, B., AND POMEAU, Y. 2008. Elasticity and geometry: from hair curls to the nonlinear response of shells. Oxford Uni- versity Press. To appear. HANSON, A. J. 2006. Visualizing Quaternions. MK/Elsevier. KLAPPER, I. 1996. Biological applications of the dynamics of twisted elastic rods. J. Comp. Phys. 125, 2, 325–337. AUDOLY, B., CLAUVELIN, N., AND NEUKIRCH, S. 2007. Elastic knots. Physical Review Letters 99, 16, 164301. LANCZOS, C. 1986. The Variational Principles of Mechanics, 4th ed. Dover. AUDOLY, B., CLAUVELIN, N., AND NEUKIRCH, S. 2008. Insta- bilities of elastic knots under twist. In preparation. LANGER, J., AND SINGER, D. A. 1996. Lagrangian aspects of the Kirchhoff elastic rod. SIAM Review 38, 4, 605–618. BERTAILS, F., AUDOLY, B., CANI, M.-P., QUERLEUX, B., LEROY, F., AND L´ EVˆ EQUE, J.-L. 2006. Super-helices for pre- dicting the dynamics of natural hair. In ACM TOG, 1180–1187. LENOIR, J., COTIN, S., DURIEZ, C., AND NEUMANN, P. 2006. Interactive physically-based simulation of catheter and guidewire. Computer & Graphics 30, 3, 417–423. BOBENKO, A. I., AND SURIS, Y. B. 1999. Discrete time la- grangian mechanics on lie groups, with an application to the la- grange top. Comm. in Mathematical Physics 204, 1, 147–188. LIN, C.-C., AND SCHWETLICK, H. R. 2004. On the geometric flow of Kirchhoff elastic rods. SJAM 65, 2, 720–736. BROWN, J., LATOMBE, J.-C., AND MONTGOMERY, K. 2004. Real-time knot-tying simulation. Vis. Comp. V20, 2, 165–179. MADDOCKS, J. H., AND DICHMANN, D. J. 1994. Conservation laws in the dynamics of rods. Journal of elasticity 34, 83–96. CHANG, J., SHEPHERD, D. X., AND ZHANG, J. J. Cosserat-beam-based dynamic response modelling. Computer Animation and Virtual Worlds 18, 4–5, 429–436. 2007. MADDOCKS, J. H. 1984. Stability of nonlinearly elastic rods. Archive for Rational Mechanics and Analysis 85, 4, 311–354. PAI, D. K. 2002. STRANDS: interactive simulation of thin solids using Cosserat models. CGF 21, 3. CHOE, B., CHOI, M. G., AND KO, H.-S. 2005. Simulating com- plex hair with robust collision handling. In SCA ’05, 153–160. PHILLIPS, J., LADD, A., AND KAVRAKI, L. 2002. Simulated knot tying. In Proceedings of ICRA 2002, vol. 1, 841–846. DE VRIES, R. 2005. Evaluating changes of writhe in computer simulations of supercoiled DNA. J. Chem. Phys. 122, 064905. PRESS, W. H., VETTERLING, W. T., TEUKOLSKY, S. A., AND FLANNERY, B. P. 2002. Numerical Recipes in C++: the art of scientific computing. Cambridge University Press. FALK, R. S., AND XU, J.-M. 1995. Convergence of a second- order scheme for the nonlinear dynamical equations of elastic rods. SJNA 32, 4, 1185–1209. RUBIN, M. B. 2000. Cosserat Theories: Shells, Rods and Points. Kluwer Academic Publishers. FULLER, F. B. 1971. The writhing number of a space curve. PNAS 68, 4, 815–819. SHABANA, A. 2001. Computational Dynamics, 2nd ed. John Wiley & Sons, New York. FULLER, F. B. 1978. Decomposition of the linking number of a closed ribbon: a problem from molecular biology. PNAS 75, 8, 3557–3561. SMITH, R., 2008. Open dynamics engine. http://www.ode.org/. SOMMERFELD, A. 1964. Mechanics of deformable bodies, vol. 2 of Lectures on theoretical physics. Academic Press. GOLDENTHAL, R., HARMON, D., FATTAL, R., BERCOVIER, M., AND GRINSPUN, E. 2007. Efficient simulation of inextensible cloth. In ACM TOG, 49. SPILLMANN, J., AND TESCHNER, M. 2007. Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects. In SCA ’07, 63–72. GOLDSTEIN, R. E., AND LANGER, S. A. 1995. Nonlinear dynam- ics of stiff polymers. Phys. Rev. Lett. 75, 6 (Aug), 1094–1097. SPILLMANN, J., AND TESCHNER, M. 2008. An adaptative contact model for the robust simulation of knots. CGF 27. GORIELY, A., AND TABOR, M. 1998. Spontaneous helix hand re- versal and tendril perversion in climbing plants. Physical Review Letters 80, 7 (Feb), 1564–1567. TERZOPOULOS, D., PLATT, J., BARR, A., AND FLEISCHER, K. 1987. Elastically deformable models. In Proceedings of SIG- GRAPH ’87, 205–214. GORIELY, A. 2006. Twisted elastic rings and the rediscoveries of Michell’s instability. Journal of Elasticity 84, 281–299. THEETTEN, A., GRISONI, L., ANDRIOT, C., AND BARSKY, B. 2006. Geometrically exact dynamic splines. Tech. rep., INRIA. GOYAL, S., PERKINS, N. C., AND LEE, C. L. 2007. Non-linear dynamic intertwining of rods with self-contact. International Journal of Non-Linear Mechanics In Press, Corrected Proof. VAN DER HEIJDEN, G. H. M., AND THOMPSON, J. M. T. 2000. Helical and localised buckling in twisted rods: A unified analysis of the symmetric case. Nonlinear Dynamics 21, 1, 71–99. GR´ EGOIRE, M., AND SCH¨ OMER, E. 2007. Interactive simulation of one-dimensional flexible parts. CAD 39, 8, 694–707. VAN DER HEIJDEN, G. H. M., NEUKIRCH, S., GOSS, V. G. A., AND THOMPSON, J. M. T. 2003. Instability and self-contact phenomena in the writhing of clamped rods. Int. J. Mech. Sci. 45, 161–196. GRINSPUN, E., Ed. 2006. Discrete differential geometry: an ap- plied introduction. ACM SIGGRAPH 2006 Courses. HADAP, S., CANI, M.-P., LIN, M., KIM, T.-Y., BERTAILS, F., MARSCHNER, S., WARD, K., AND KAˇ CI´ C-ALESI´ C, Z. 2007. Strands and hair: modeling, animation, and rendering. ACM SIGGRAPH 2007 Courses, 1–150. WARD, K., BERTAILS, F., KIM, T.-Y., MARSCHNER, S. R., CANI, M.-P., AND LIN, M. C. 2007. A survey on hair model- ing: Styling, simulation, and rendering. IEEE TVCG 13, 2 (Mar./ Apr.), 213–234. HADAP, S. 2006. Oriented strands: dynamics of stiff multi-body system. In SCA ’06, 91–100. YANG, Y., TOBIAS, I., AND OLSON, W. K. 1993. Finite element analysis of DNA supercoiling. J. Chem. Phys. 98, 2, 1673–1686.