Download Discrete Elastic Rods: A Geometric Model for Thin Flexible Rods and more Lecture notes Discrete Mathematics in PDF only on Docsity!
Discrete Elastic Rods
Mikl´os Bergou
Columbia University
Max Wardetzky
Freie Universit¨at Berlin
Stephen Robinson
Columbia University
Basile Audoly
CNRS / UPMC Univ Paris 06
Eitan Grinspun
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
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- pling between bending and twisting modes. We validate the discrete rod model via quantitative buckling, stability, and coupled-mode experiments, and via qualitative knot-tying comparisons.
CR Categories: I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism—Animation
Keywords: rods, strands, discrete holonomy, discrete differential geometry
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].
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.
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.
Elegant model of elastic rods We build on a representation 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.
Efficient quasistatic treatment of material frame We addition- 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 introduce
Figure 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.
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 approach in § 8 handles both types of constraints, emphasizing phys- 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.
Validation We show how to model and simulate inextensi- 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 ).
2 Related work
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; Lin and Schwetlick 2004] with applications spanning medicine [Lenoir et al. 2006], biology [Yang et al. 1993; Klapper 1996], the study of knots [Phillips et al. 2002; Brown et al. 2004], and computer 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 ].
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-
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.
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.
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-
Algorithm 1 Discrete elastic rod simulation Require: u^0 // Bishop frame vector in frame { t^0 , u^0 , v^0 } at edge 0 Require: x 0... x n + 1 // position of centerline in rest state Require: ( x 0 , ˙ x 0 )... ( x n + 1 , x ˙ n + 1 ) // initial position/velocity of centerline Require: boundary conditions // free, clamped or body-coupled ends 1: precompute ωωω (^) ij using ( 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) // [Smith 2008] 6: compute forces on centerline (§7.1) 7: integrate centerline (§7.2) // [Hairer et al. 2006] 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
t i e i
t i-^1 e i-^1
π−φ i
x i- 1
x i+ 1
x i
Figure 4: Notation Angles and vectors used in our discrete model.
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.
Curve-angle representation The (twist-free) Bishop 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 ):
m 1 = cos θ · u + sin θ · v , m 2 = − sin θ · u + cos θ · 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 = ( m 1 )′^ · m 2 and u ′^ · v = 0. Hence, we write the twist energy as
E twist(Γ) =
∫ β ( θ ′)^2 d s.
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.
4.2 Discrete Kirchhoff rods
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.
Primal vs. dual We distinguish between primal quantities , as- 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:
e^0 e^1 e^2 e n
x 0
x 1 x 2
x 3 x n x n+
Discrete framed curves A discrete framed curve, Γ, consists of a centerline comprised of ( n + 2 ) vertices x 0 ,... , x n + 1 and ( n + 1 ) straight edges e^0 ,... , e n^ such that e i^ = x i + 1 − x i , together with an assignment of material frames M i^ = { t i , m i 1 , m i 2 } per edge (see
Fig. 4 ). We consider adapted frames, where t i^ = e i /| e i | is the unit tangent vector per edge. Notice that tangent vectors of polygonal curves are uniquely defined at edges, whereas their definition at
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.
Pointwise vs. integrated quantities In the smooth setting, we 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].
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 D i as the Voronoi region associated to each vertex, having length |D i | = li /2, where li = | e i −^1 | + | e i |.
4.2.1 Discrete bending energy
Recall that the key player for formulating elastic bending energy in the smooth case was the curvature of the centerline.
Discrete curvature Since each edge is straight, it follows that 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
κ i = 2 tan
φ 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.
Discrete curvature binormal We define the curvature binormal at a vertex (an integrated quantity ) by
( κ b ) i =
2 e i −^1 × e i | e i −^1 || e i | + e i −^1 · e i^
Observe that the vector ( κ b ) i is orthogonal to the osculating plane passing through two consecutive edges and has magnitude 2 tan( φ i / 2 ). In particular, this quantity is well-defined in the case of collinear edges, even though the binormal itself is not.
Discrete bending energy We now have all the pieces to assem- ble the bending energy of a discrete naturally straight, isotropic rod:
E bend(Γ) =
n
i = 1
α
κ b i li / 2
li 2
n
i = 1
α( κ b i )^2 li
The division by li /2 is due to converting the integrated quantity κ b i into a pointwise one (see above), whereas the multiplication by li / 2 takes care of the measure of the domain of integration, D i.
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 ) i as defined by ( 1 ), we express the material curvatures as
ωωω (^) ij =
( κ b ) i · m 2 j , −( κ b ) i · m 1 j
) T
for j ∈ { i − 1 , i }. (2)
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.
With these definitions, the bending energy of a discrete rod is
E bend(Γ) =
n
i = 1
2 li
i
j = i − 1
ωωω (^) ij − ωωω (^) ij
) T
B
j
ωωω (^) ij − ωωω (^) ij
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 B j = α Id 2 × 2 and ωωω j i =^00 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.
Discrete parallel transport and Bishop frame We define discrete parallel transport in analogy to the smooth case, i.e. , as a rotation Pi about the curvature binormal,
Pi
t i −^1
= t i , Pi
t i −^1 × t i
= t i −^1 × t i^.
By convention, Pi is the identity if t i −^1 = t i , whereas Pi is not defined if t i −^1 = − t i. Par- allel transport is a key notion that allows us to generalize the notion of twist from the smooth to the discrete setting.
In order to define Bishop frames , we draw upon the smooth case by transporting a unit vector u^0 , which is orthogonal to t^0 , along the rod using our rotation operators Pi , i.e. , we iteratively define
u i^ = Pi
u i −^1
The third axis of each Bishop frame is then v i^ = t i^ × u i.
Bishop frame update The Bishop frame is by definition adapted to the centerline. This requires that u^0 ⊥ t^0 , 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 u^0 ⊥ t^0 may be violated (the exception is when t^0 is 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.
Material frame representation Since both the material frame and Bishop frame are defined at edges, let θ i^ denote the angle needed to rotate the Bishop frame into the material frame at edge e i (see Fig. 4 ). Then the material frame vectors at edge i are
m i 1 = cos θ i^ · u i^ + sin θ i^ · v i^ , and m i 2 = − sin θ i^ · u i^ + cos θ i^ · v i^.
Discrete twisting energy In perfect analogy to the curve-angle representation in the smooth case, we adopt the twisting energy
E twist(Γ) =
n
i = 1
β
( θ i^ − θ i −^1 )^2 li
n
i = 1
β m^2 i li
where the scalar mi = θ i^ − θ i −^1 is the ( integrated ) discrete twist. The variable mi measures the angle between two adapted frames: the result of parallel transporting the material frame from edge e i −^1 to e i^ and the material frame at edge e i^ itself.
5 Quasistatic material frame postulation
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.
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.
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:
∂ E (Γ) ∂ θ j^
Note that the quasistatic treatment of the material frame eliminates the θ j^ variables 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.
Boundary conditions The value of θ j^ for an edge may be pre- 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 θ j^ variables 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.
5.1 Quasistatic update
During simulation, we enforce the quasistatic nature of the material frame by ensuring that ( 4 ) is satisfied before forces are computed.
Special case of naturally straight, isotropic rods For an isotropic rod with a naturally straight undeformed configuration, ( 4 ) implies that a clamped rod has uniform twist,
mi li
θ n^ − θ 0 2 L
= constant , (5)
where 2 L = (^) ∑ ni = 1 li , and the pointwise twist mi / li (recall §4.2.2) 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
E (Γ) =
n
i = 1
α( κ b )^2 i li
β
θ n^ − θ 0
2 L
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.
General case For the general case, we use Newton’s method to minimize the elastic energy E (Γ) with respect to the material
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:
F^0 ( ε)
F^1 ( ε)
· · · F j −^1 ( ε)
F j ( ε)
ψ 1 ( ε) ψ 2 ( ε) ψ (^) j ( ε)
F^0 ( 0 )
Ψ j
Here, the Bishop frames F i^ = { t i , u i , v i } 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 ( F i −^1 ) to F i ; the first vertical arrow is also labeled with 0 because we ensure that F^0 is always updated via parallel transport (see §4.2.2). We are interested in the angle of rotation, Ψ j , required to align ˜ P j ( ε)
F j ( 0 )
to F j ( ε).
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^ = (^) ∑ (^) ij = 1 ψ i. For computing forces, we require the gradient of this angle with respect to vertex positions, which is given by
∇ i Ψ j^ =
j
k = 1
∇ i ψ k. (10)
By ( 9 ), this sum will have at most three non-zero terms.
Discussion In hindsight, we view the relation of holonomy to 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.
7 Equations of motion
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.
7.1 Forces on centerline
In deriving the forces on the rod’s centerline, we must consider how the energy depends on the x i variables—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 x i is therefore given by
d E (Γ) d x i
∂ E (Γ)
∂ x i
n
j = 0
∂ E (Γ)
∂ θ j
∂ θ j ∂ x i
Here, the total derivative d E /d x i takes into account the implicit dependence of potential energy on centerline positions, whereas the partial derivative only takes into account explicit dependence.
Recall from ( 4 ) that ∂ E / ∂ θ j^ vanishes for all edges for which θ j is not prescribed by a boundary condition. Therefore, for stressfree
boundary conditions, the sum is zero. For clamped boundaries, only one component of this sum could be non-zero, allowing us to write the force as
−
∂ E
∂ x i
∂ E
∂ θ n
∂ θ n ∂ x i
To evaluate ∂ θ n / ∂ x i , recall from §6.1 that varying the centerline rotates the Bishop frame at edge e n^ by angle Ψ n^ about the tangent t n. Therefore, to keep the material frame aligned to the boundary condition, we must subtract this angle from θ n^ to obtain
∂ E
∂ x i
∂ E
∂ θ n
n
j = 1
∂ ψ (^) j ∂ x i
For clamped boundaries, we give the components required to eval- uate this force for both the special and general case.
Special case For naturally straight, isotropic rods, the forces on vertex x i are given by up to three contributions:
2 α l (^) j
∇ i ( κ b ) (^) j
) T
( κ b ) (^) j +
β ( θ n^ − θ 0 ) L
∇ i ψ (^) j i − 1 ≤ j ≤ i + 1 ,
with the gradient of holonomy given by ( 9 ) and the gradient of the curvature binormal given by
∇ i − 1 ( κ b ) i =
2 [ e i ] + ( κ b i )( e i ) T | e i −^1 || e i | + e i −^1 · e i^
∇ i + 1 ( κ b ) i =
2 [ e i −^1 ] − ( κ b i )( e i −^1 ) T | e i −^1 || e i | + e i −^1 · e i^
∇ 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.
General case For anisotropic rods with a curved rest shape, the required expressions for the forces are given by:
∂ E ∂ θ n^
ln
( ωωω nn ) T^ JB n ( ωωω nn − ωωω nn ) +
2 β mn ln
and
∂ E
∂ x i
n
k = 1
lk
k
j = k − 1
∇ i ωωω (^) kj
) T
B
j
ωωω (^) kj − ωωω (^) kj
where the gradient of the material-frame curvature ωωω j k is given by
∇ i ωωω j k =
( m j 2 )
T
−( m j 1 )
T
∇ i ( κ b ) k − J ωωω j k (∇ i Ψ^
j ) T (^). (11)
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.
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
M ¨ x = −
d E (Γ) d x
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.
8 Constraints
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.
Inextensibility constraints For each edge of the rod, we use the quadratic constraint equations e i^ · e i^ − e i^ · e i^ = 0, where (the pre- computed) e i^ refers to the undeformed configuration.
Rigid-body coupling constraints The welding of the body to 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:
q · q − 1 = 0 , q x 0 q ∗^ + r − x 0 = 0 , q x 1 q ∗^ + r − x 1 = 0 ,
where r ∈ R^3 is 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 qx 0 q ∗^ is a rotation of x 0 that 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.
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].
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.
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).
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.
Fast manifold projection For the enforcement step, we adopt 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
both positional as well as rotational degrees of freedom by consid- ering the metric induced by the kinetic energy 12 y M ˜ y T^ , where the ( 3 n + 12 ) × ( 3 n + 12 ) generalized mass matrix and the generalized velocity y ∈ R^3 n +^12 are defined respectively by
M ˜ =
4 · I
M · Id 3 × 3 M
(^) and y = ( q −^1 q ˙, r ˙, ˙ x ).
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) moment of inertia tensor expressed as a 3×3 matrix in the reference coordinates. Since q is a unit quaternion, q −^1 q ˙ corresponds to a vector in R^3 [Hanson 2006]. The kinetic energy has contributions from the body rotation, body translation, and centerline translation.
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 M ˜ y T^ − C T λλλ , 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−^8 in satisfying the constraint C = 0.
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 ) −
h
( 2 q 0 q −^1 , r 0 − r , x 0 − x ).
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 m^0 1 = qm^0 1 q ∗.
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 m^0 1 , to the coupled rigid-body. Note that the rigid-body’s relatively large moment of inertia ensures that the explicit coupling is non-stiff.
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, τ = | τ| t^0 , 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 e^0. Note that δ θ 0 also af- fects the body’s orientation, due to the rod-body constraint. The
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!
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.
Plectonemes A rod generally forms entangled structures called 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.
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
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.
Asymmetry of twist The combination of anisotropy and non- 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
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.
Helical perversion: Perversion is a classical pattern that can be 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 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 = 2 Id 2 × 2 and β = 2, with each link in the chain having a radius of approximately 1 unit.
Hanging chain: free boundaries In Figure 10 , we show a chain 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.
9.4 Rigid-body coupling
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.
Wilberforce pendulum This experiment reveals fascinating as- 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
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.
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.
Torque transducer Rods can be used to couple rigid-bodies, 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.
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.
Rigid-bodies for rigid couples Coupling rods with rigid-bodies opens the possibility for complex simulations—in particular, for
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. 5 49 1.0 0.025 0.018 0.021 0.12 0. 6 200 2.0 0.072 0.09 0.064 0.81 1. 9 275 1.0 0.13 0.19 0.17 0.42 0. 7 75 2.0 0.043 0.1 0.2 0. 8 67 50.0 0.039 0.096 0.28 0. 10 31 1.0 0.016 0.13 0.13 0. 11 20 1.0 0.0036 0.0043 0.019 0. 12 2158 0.1 0.87 0.64 21.0 22.
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.
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.
10 Conclusion
Limitations and future work Our use of the Fast Projection 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.
Acknowledgments We would like to thank Michael Herrmann for his indispensable insights on unifying the isotropic and anisotropic picture, S´ebastien Neukirch for his help in reviewing the 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-
- and the DFG Research Center MATHEON in Berlin, as well as by Autodesk, Elsevier, mental images, NVIDIA, and Walt Dis- ney Animation Studios.
