Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014), pp. 1–10
Vladlen Koltun and Eftychios Sifakis (Editors)

Position-based Elastic Rods
Nobuyuki Umetani

Ryan Schmidt

Jos Stam

Autodesk Research

Abstract
We present a novel method to simulate complex bending and twisting of elastic rods. Elastic rods are commonly
simulated using force based methods, such as the finite element method. These methods are accurate, but do
not directly fit into the more efficient position-based dynamics framework, since the definition of material frames
are not entirely based on positions. We introduce ghost points, which are additional points defined on edges, tonaturally endow continuous material frames on discretized rods. We achieve robustness by a novel discretizationof the Cosserat theory. The method supports coupling with a frame, a triangle, and a rigid body at the rod’s end
point. Our formulation is highly efficient, capable of simulating hundreds of strands in real-time.
Categories and Subject Descriptors (according to ACM CCS):
modeling—Animation

I.6.8 [Computer Graphics]: Simulation and

1. Introduction
Position-based dynamics (PBD) has been widely acceptedin the field of computer animation due to its efficiency, robustness and simplicity. The goal of the PBD is not to simulate physics as accurately as possible, but rather to sacrifice some quantitative accuracy to generate visually plausible simulation results very quickly. To this end, PBD has
broadly been applied in many game engines and visual effects, where speed and controllability is crucial. PBD has
primarily been used to simulate various physical phenomena
associated with solid and thin-shell (i.e., clothing) deformations [MHHR07]. We present a new application of PBD toelastic rod simulation, which is essential for animating thinstrands such as hair, fur, ropes, and so on.
Because of the existence of twist, rod simulations are
more complicated than solid or thin-shell simulations, making it hard to implement them in the PBD framework. PBD
requires an object’s deformations be characterized with discrete positions of points. However twist of a rod cannot
be directly specified with only positions of the rod center line, because twist requires angular information –rather
than positional– which describes how much the material
is twisted around the center line. Since complex nonlinear
deformation, such as out-of-plane buckling under torsionalstrain (e.g., Fig. 2), is caused by an interaction between twist
and bending, handling of twist is critical in nonlinear deformation of rods.
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

Figure 1: A squishy ball hits a wall. Tentacles of the squishy
ball was modeled with our position-based elastic rods.
Representation of twist in elastic rods is actively researched, but so far no previous work has successfully
implemented it in a PBD framework. While bending andstretching can be characterized by the position of the centerline, the twist is represented using material frames definedon the centerline. Constructing a material frame is difficult
because its orientation is coupled with the centerline. Various approaches have been proposed. For example, Spillman
et al. [ST07] use the penalty method to orient the material
frame to rods’ tangential direction. [BWR∗ 08] defined a material frame by rotating a Bishop frame in the direction of
the rod’s tangential direction. However, theses approaches
are not suitable for PBD in which all deformations are described via discrete point positions, not frame orientations.
Our goal is to incorporate the twisting and bending
physics of elastic rods into the PBD framework. With proper
handling of twist, we can produce (i) visually plausible sim-

2

N.Umetani et al. / Position-based Elastic Rods

ulation of rod deformations, capable of reproducing various phenomena, which are (ii) fully controllable by physically meaningful parameters, (iii) the PBD framework ensures robustness, speed and simple implementation. Our
position-based elastic rod model is derived from Cosserat
theory [CC07], which describes the physical principles ofelastic rods. Even though our model may not quantitatively
agree with real-world experiments, it reproduces plausible
deformation such as looping phenomena.
Our key contribution is to introduce ghost points, which
are placed on edges to represents material points distributed
around the edge. Using ghost points, we can represent thematerial frame so that its orientation naturally follows therods’ tangential direction, leading to an efficient formulation. The dynamics of these ghost points are treated straightforwardly in the framework of PBD, hence our implementation is simple and fast.
We propose a novel discretization of the Cosserat theory,
which further enhances robustness and simplicity of the simulation. Qith this enhancement, we can enforce the twisting
and bending forces without computing expensive trigonometric functions or their inverses. Moreover, our formula
prevents the rod from flipping in high-torsion situations. We
extend our discretization to also specify boundary conditions, in which the rod is connected to a rigid body or atriangle with one- or two-way coupling.
We find that we can interpret the traditional constraint enforcement procedure of PBD in a variational formulation.
This allows us to extend PBD so that it can handle multiple constraints simultaneously, resulting in a significant
speed-up. Furthermore, because the material’s deformation
is usually written in a vector or tensor form, our variational
approach enables us to manipulate the multi-dimensional
anisotropic stiffness parameters of the material.
Our formulation is very efficient. For example, on a modern workstation we can simulate 200 strands of hair, discretized with 20 edges, at an interactive rate. Our rods arealso highly stable, to the point of being able to recover froma random initial deformation without kinks or instability. We
demonstrate the flexibility of our approach by showing various examples coupling rods with other rigid and deformableobjects. Our contributions include:• Ghost points to discretize material frames on a rod which
are directly compatible with position-based dynamics.
• New robust discrete elastic rod formulation based on
Cosserat theory which also can specify boundary conditions of rods.
• New variational interpretation of the PBD constraint enforcement, which enables fusing multiple constraints.
2. Related Work
Since Pai et al. [Pai02] first introduced the Cosserat theoryof elastic rod simulation into graphics, it has remained an

Figure 2: Rods show dynamic looping phenomenon known
as Plectoneme when the both ends of a rod is clamped andtwisted. This nonlinear behavior cannot be reproduced witha model without coupling of twist and bending.
active research topic. Most elastic rod models are based on
Cosserat theory, which describes nonlinear deformation ofrods. A primary difference between formulations is whether
the centerline is represented implicitly or explicitly.
Implicit discretization of Cosserat theory. This approach
uses curvatures to represent curved rod shape and recovers
the rods’ explicit centerline shape by integrating frames fromone end. The first model from Pai et al. [Pai02] falls into
this category. Later, this approach was further extended forhelical elastic rods [BAC∗ 06] to simulate curled hairs. Recently, a super-clothoid model has been proposed [CBD13]
that uses power-series to integrate curvature to find the center line of a clothoid with a stability guarantee. Implicit models achieve high-order smoothness of the rod with fewer
DoFs and are suitable for hair simulations where one end
is clamped and another end is free. However, since this approach solves for derivative information of material frames,
it cannot be combined with position based approaches.
Explicit discretization of Cossarat theory. Our methods
are greatly inspired by the CORDE model proposed by Spillman et al. [ST07]. We both begin with an explicit representation of the centerline and define material frames to apply
Cosserart theory. However, Spillman et al. define material
frames independently from the centerline position, and thus
require penalty-based methods to make the frames follow thecenterline. This approach does not directly fit into the PBD
framework because the variables have frame orientation, in
addition to the position of the centerline.
Formulation based on Bishop frame. In the explicit centerline model, the definition of a material frame that tracks
the rods’ tangential direction is not straightforward. The discrete elastic rod proposed by Bergou et al. [BWR∗ 08] wasa breakthrough because their material frames, constructedfrom the Bishop frame, automatically follows the deformation of the rod center line. The Bishop frame construction
was later transferred from the spatial domain to the temporal
domain to enhance efficiency [BAV∗ 10, KJM10]. However
as transported frames are not straightforward to represent ina point-based formulation, formulating constraints for PBD
is a challenge with Bishop-based approaches.
Mass-spring model. Several studies [SLF08, MSW∗ 09]
modeled hair deformation with network of masses connectedsubmitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

3

N.Umetani et al. / Position-based Elastic Rodsby springs. The mass-spring model is simple and robust, but
it is difficult to separate the bending stiffness from twisting stiffness. In contrast, our model is highly controllable;
the properties of rods can be manipulated with independent,
physically meaningful parameters.
Position-based dynamics. PBD provides fast and robust
simulation tailored for computer animations [MHHR07,
BMOT13]. The core ideas, such as Verlet integration andstiffness based on constraints, are presented in [Jak01].
While PBD is not quantitatively accurate, in that it cannot reproduce real-world experimental data, PBD animation is sufficiently plausible for computer animations. PBD has been
applied in game engines such as Nvidia PhysXTM and visual effect engines such as Nucleus [Sta09]. The numerous
games and films produced with these tools are a testament to
the importance and utility of PBD. Recently, application of
PBD has been extended to fluid dynamics [MM13]. Muller
et al. [MKC12] presented a fast method to apply inextensible
constraints in the PBD for hair and fur without the twisting
and bending forces. To our knowledge, twisting and bendingof elastic rods have yet to be solved in the PBD framework.
Shape matching method [MHTG05] is the most similar to
PBD, and in that domain several studies tackle simulation ofelastic rods. The model in [RKM∗ 12] represents the elastic
rod by appending twisting force, computed by method using Bishop frame, to the shape matching deformation. The
twisting force is modeled as an external force hence thismethod doesn’t benefit from shape matching techniques unconditional stability. Oriented particles [MC11] consider rotational inertia at the nodes and perform rod simulation stably, but do not handle the coupling of bending and twisting. Shape matching method treats elasticity as a force to getback to the rigid goal shape, thus cannot easily handle coupling of twisting and bending where the rigid shape cannot
be a goal shape.
3. Discrete formulation of rod
Like other PBD formulations, our position-based elastic rod
handles elasticity in the form of constraints. Similar to theGauss-Seidel method, the particle locations are iteratively
updated, so that linearized constraints are locally satisfied.
Hence, how to construct constraint for twisting and bendingis essential for elastic rod simulation.
We our technique is based on the Cosserat theory, which
we briefly review (§ 3.1). To define differentiable continuous
material frames on a rod in a discrete setting, we introduce
ghost points (§ 3.2). Next, we describe a novel discretization of Cosserat theory that makes simulation simpler andfaster (§ 3.3). We extend the types of constraints that PBD
can handle by introducing a variational interpretation (§ 4).
Using this extension, constraints for our elastic rod are described (§ 5). Finally, we describe end point constraints that
fix the position and orientation of a rod at its end (§ 5.1).
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

3.1. Cosserat theory
We briefly review elastic rod simulation based on the
Cosserat theory. For a more detailed discussion, please refer to [Ant05]. We assume that a rod is parameterized byarclength s ∈ [0, L̄], where L̄ is the initial rod length. The
rod’s one dimensional configuration is described by the centerline, where p(s) ∈ R3 is the position of the rod’s centroid
at the location s.
The twist of a rod cannot be represented by the centerline
positions. Thus, at every point s on the rod, we define anorthonormal material bases d1 (s), d2 (s), d3 (s). A material
frame is a 3×3 orthogonal matrix D(s) = [d1 , d2 , d3 ] which
has material bases as columns. In differential geometry, theDarboux vector describes how the material frame evolves
along the curve, and is a key component of Cosserat theory.
The Darboux vector is an axial vector of frame rotation withrespect to change of s. More specifically, using the bases ofthe material frame, the Darboux vector can be written as3

ω (s) =

1X k
d (s) × dk (s)0 ,
2

(1)

k=1

where the prime ( 0 ) denotes derivative with respect to s.
Material coordinates of the Darboux vector at a location s
become ωi (s) = ω (s) · di (s). This ωi encodes how the frame
rotates around the material frame basis di . In Cosserat theory, the rod’s bending and twisting energy is defined using
the material coordinates of the Darboux vector as
Z L̄ X3 X
31V=(ωi − ω̄i )Ki j (ω j − ω̄ j )dt,(2)
0 2
i=1 j=1

where ω̄i is the material coordinate of the Darboux vector
in an initial zero-stress configuration, and Ki j ∈ R, i, j ∈{1, 2, 3} is a positive definite symmetric stiffness matrix.
The initial material frame D(s) can be chosen arbitrarily.
In practice, the material coordinates are chosen such that oneof the basis vectors, such as d3 , corresponds the rod’s tangential direction so that p(s)0 become parallel to d3 ( thisframe choice is called the adapted frame [BWR∗ 08]). Next,
d1 and d2 are chosen as the principal axes of the second moment of area computed for the rod’s cross-section. In thisparticular setting, assuming isotropic homogeneous material, K becomes diagonal and the K11 and K22 correspondto stiffness of bending around the axes d1 and d2 , and K33
specifies the twisting stiffness.
Material frames must track the motion of the centerline
over time. More specifically, when a rod changes its centerline, the material frames should move so that the tangential
direction always corresponds to the basis d3 . This constraint
on the material frame complicates the simulation of rods. In
contrast to the penalty method used in the CORDE model,
the construction of the material frame that we describe in
the following section naturally satisfies the above condition.

4

N.Umetani et al. / Position-based Elastic Rods

3.2. Continuous material frames on a discrete rod
We explicitly model the rod’s centerline as a polyline with
N + 1 vertices {p0 , p1 , . . . , pN }. Edge e ∈ {1, . . . , N} has end
points pe−1 and pe , and masses {m0 , . . . , mN } are associatedwith the vertices. As shown in Fig. 3, we assume each edgege has a ghost point pe ∈ R3 , which is placed on the perggpendicular bisector of each edge. Masses {m1 , . . . , mN } areassociated with the ghost points .

The vector Θ parameterizes the rotation, hence we canscale rotation by scaling this vector with the ratio inside oriΘ]r), where r = (s − se )/l¯e .
entation element as R(s) = exp([ΘHere, l¯e is the initial length of the orientation element e andse is the arclength parameter s at the midpoint of the edgee. With this uniform scaling, bases of the frame are interpolated with s asΘ]r)die ,di (s) = exp([Θ(5)

iFigure 3: Configuration of edges, points (red) and ghost
points (cyan).
Within the PBD framework, our ghost points are treatedin the same manner as any other points. The distance between the ghost point and the edge is set to a fixed value l¯g ,which we set to the average edge length. Ghost points represent material distributed around the edges, thus inertia forg
the rotation around an edge becomes l¯g2 me . As a result, we
can control rotational inertia by changing l¯g .
Material frame at the center of an edge. Using the ghost
points, we define frame De at the center of edge e as
 3
 de = (pe − pe−1 ) / |pe − pe−1 |,ggd2 = d3e × (pe − pe−1 ) / |d3e × (pe − pe−1 )|, (3) e1de = d2e × d3e .
For the rest of the rod, the frame is continuously interpolated from these edge centers. Following [ST07], we define
an orientation element e ( e ∈ {1, . . . , N − 1} ) as a segment
starting at the midpoint of edge e and ending at the midpointof edge e + 1. Within the orientation element, the frame isinterpolated from the two end frames De = [d1e , d2e , d3e ] andDe+1 = [d1e+1 , d2e+1 , d3e+1 ].
Interpolation of frames. In order to compute the Darboux
vector and its material coordinates for Cosserat theory, weinterpolate between two frames De and De+1 inside the orientation element e (see Fig. 4-left). We assume that thematerial frame changes uniformly from De into De+1 (i.e.,
spherical linear interpolation). This uniform interpolation
between two frames can be easily computed using an axisangle formulation. The rotation matrix between two frames
R = De+1 (De )T is represented with a 3-vector Θ = θn,where θ ∈ [0, π) is a rotation angle and n is an unit axis vector (see Fig. 4-middle). Using the matrix exponential, therotation matrix can then be written asΘ]),
R = exp([Θ(4)

where “[·]” denotes the cross product matrix (i.e., [a] is a 3x3
skew-symmetric matrix as [a]b = a × b, a, b ∈ R3 ).

The differentiation of basis d (s) with respect to s can be¯ i (s). Thus, from (1), the Darboux
Θ]/l)dwritten as di (s)0 = ([Θvector can be written as ω (s) = Θ /l,¯ which is constant overs. Note that the constant Darboux vector means this orientation element is modeled as a helix. However our discretization is different from super-helix [BAC∗ 06] in that helices
are an internal representation to model bending and twisting; they are not used to represent outer shape of the rod. In
fact, our helices are not connected to each other.
Because the basis di is rotated around the axis vector n,
the angle between the axis vector and the basis does notchange. Hence, the inner product of the two vectors n anddi (s) is constant with s (see Fig. 4-right). Thus, the material
coordinate of the Darboux vector is also constant as11
ωi (s) = ¯ Θ · die = ¯ Θ · die+1 , i = 1, 2, 3.
(6)
l
l

θr
θ(1-r)

Figure 4: Left: configuration of points and frames. Middle:
the rotation between neighboring frames. Right: rotation ofbasis around an axis.
3.3. Modified frame twist representation
So far, we described a straightforward way to obtain the
Darboux vector using axis-angle parameterization (R → Θ ).
This parameterization involves inverse trigonometric functions and is computationally intensive for PBD, where constraint application may be iterated hundreds of times within
a single time step. In this section, we further simplify ourmodel - while increasing robustness - by modifying the axisangle parameterization Θ of the rotational matrix R.
To obtain a better parameterization of Θ , we start from
Rodrigues’ rotation formula, which describes the relationship between the rotation matrix and axis-angle parameter:
R = I + [n] sin θ + (1 − cos θ)(nnT − I),(7)

where I is an identity matrix. To obtain cos θ, we take thesubmitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

5

N.Umetani et al. / Position-based Elastic Rods

trace of both sides of (7), giving
1 + 2 cos θ = trR =3
X
dne · dne+1 .

(8)

n=1

Here we use the relationship tr(XY) = tr(XY), which isapplicable for square matrices X and Y. Let “vect(·)”
be an operator R3×3 7→ R3 which maps a matrix to an
axis vector which is the anti-symmetric part of the matrix(i.e., vect(A) = 1/2 (A32 − A23 , A13 − A31 , A21 − A12 )). Using Rodrigues’ parameter [MLSS94], rotation R can be parameterized with 3-vector Ψ asP
2 3k=1 dke × dke+1
4vect (R)θ
Ψ = 2n tan =
=
. (9)
P
2
1 + trR
1 + 3n=1 dne · dne+1
Here we use a trigonometric identity tan θ/2 = sin θ/(1 +cos θ). Note that Rodrigues’ parameter Ψ have the same orientation n as Θ , but the parameter has a different magnitude
Ψ| = 2 tan(θ/2).
|Ψ
Material coordinates of the Darboux vector. In our
Cosserat-based formulation, we interchange axis angle vector (nθ) with Rodrigues’ parameter Ψ in the Darboux vector
as ω (s) ' Ψ /l¯ ≡ Ω (s). We call this approximated Darboux
vector the modified discrete Darboux vector, denoted by Ω .
This modification may affect the quantitative accuracy of ourposition-based elastic rod, but brings us closer to our goal of
simplicity and robustness.
There are four motivations for using this approximation.
First, Rodrigues’ parameter Ψ is a good approximation ofΘ if the rotation angle θ is small. Second, the Ψ is easy
to compute from the coordinate bases without using any
trigonometric functions or their inverse. Third, as mentioned
in [BWR∗ 08], the material coordinate of the modified discrete Darboux vector Ωk goes to infinity when the rod iskinked completely (θ = π). In such a case, the energy also
goes to infinity, as it tries to push the rod back to the original
shape. Finally, since the magnitude 2 tan(θ/2) is a monotonically increasing convex function, this makes the relativerotation at V = 0 unique. In other words, the Rodrigues’ rotation parameter is a bijective map to the rotation matrix.
This ensures that there are no kinked rest shapes other than
the initial shape.
From (9), the local coordinate of the modified discrete
Darboux vector Ωi (s) = Ω(s)·di (s), (i = 1, 2, 3) can be written as   j kjk
4 vect(Q)i2 de · de+1 − de · de+1
Ωi = ¯= ¯, (10)
P3
1 + trQll 1 + n=1 dne · dn
e+1
where Q = (De )T De+1 and indexes i, j, k permute as{i, j, k} = {1, 2, 3}, {2, 3, 1}, {3, 1, 2}. Note that R and Q areboth rotation matrices from De to De+1 , but they are different. While R is rotation in the global coordinate, Q is relative rotation describing rotation to De+1 observed from thesubmitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

frame De . Since the coordinate of the modified discrete Darboux vector is also constant over an orientation element, thebending and twisting energy on the orientation element e canbe written as
Ve =

3
l¯e X
Kii (Ωi − Ω̄i )2 .
2

(11)

i=1

Note that in our formula, the rods’ bending and twisting energy are computed exactly for both bending and twisting
components. Hence, the bending and the twisting stiffness
parameters are completely decoupled. In other words, we
can make a rod that has zero twisting resistance with highbending stiffness or vice versa. This gives the computer animation artist fine-grained control over the rods’ behavior.

4. Variational interpretation of constraint enforcement
In the PBD framework, elasticity is formulated via constraints. These constraints take zero values for the rest shapeand increase if strain is applied. In each time step, the particle is updated in the gradient direction of the constraint, representing the internal force pushing object back to the rest
shape. Let each constraint be written in the form C(p) = 0,
where p is a concatenation of the positions of points involvein this constraint. We denote A as the set of points involving in this constraint, and we assume the constraint is a Mdimensional function. While the constraint function is scalarin the original PBD formula, here we want generalize it for
multi-dimensional function. For this purpose, we introduce
a variational formula for the constraint enforcement.
In the classical mechanics, the points trajectory underconstraints is known to comply the Gauss’s principle of least
constraint [GJS01]. This principle states that the trajectories of constrained
points under gravity g minimize a quanPtity Z = a∈A ma |p̈a − g|2 , where p̈ is a points’ accelerations. Because the p̈a − g means how much the constraintchanged the accelerations, the principle means constrained
motion takes a path that minimizes the sum of acceleration
change caused by constraints. Batty et al. [BBB07] leverage this variational principle for incompressible constraint
around coupling boundary of fluid against rigid bodies.
Let us explain how the Gauss’s principle unfolds for the
PBD time integration scheme. Let pt and vt are positions andvelocities of points at time t and ∆t is the time step. In the
PBD scheme [MHHR07], the positions in the next time step
can be computed as: pt+∆t = pt + ∆t(vt + ∆tg) + ∆p, wherethe ∆p is the position updates by constraint enforcement.
Hence, the velocity of points becomes vt+∆t = (pt+∆t −pt )/∆t = vt + ∆tg + ∆p/∆t. Finally, the acceleration of apoint can be written asp̈ = (vt+∆t − vt )/∆t = ∆p/∆t 2 + g.

(12)

By plugging in (12) in the Gauss’s principle, we obtain the6

N.Umetani et al. / Position-based Elastic Rods

following variational formula for the position update:X∆p = arg minma |∆pa |2 , where C(p + ∆p) = 0. (13)∆psquare norm is quadratic, the linearization of the constraint
is not a good approximation. Thus, more iterations are required to achieve a rod with the same stiffness.

a∈A

This means the position update minimizes weighted sum of
its squares while satisfying the constraints. This QP problem
can be easily solved by introducing a Lagrange multiplierT
λ ∈ RM . First, we add λ C to the minimization target, which
has to have zero derivative at the solution:XTma |∆pa |2 + λ C = 0 (∀b ∈ A).
(14)
∇b 
a∈A

Here ∇b means derivative with respect to the position ofpoint b. Then, with a little calculation exercise, we get the
point update using the Lagrange multiplier as:T∆pa = −0.5wa λ ∇a C,(15)

where wa is a inverse of the point weight. Nest, we put
this (15) into linearized constraint
C + ∇C∆p = 0,(16)

to solve the Lagrange multiplier. Finally, we obtain the final
position update:

−1
X
T
∆pa = (wa ∇a C) wb ∇b C ∇b C C.
(17)
b∈A

Note that when the mass of all points are equal, this update is equivalent to multiplying the C with Moore-Penrose
matrix inverse of the matrix ∇C. This formula is equal tothe original one when the constraint is a scalar function(M = 1) and thus, is the natural extension of the original
formula for solving multiple constrains at the same time. In
other words, the original PBD approach solves constraint in
the Gauss-Seidel manor, but we solved in the block GaussSeidel manor, which gives better convergence if there is a
block structure in the coefficient matrix.
This position update preserves the linear and angular momentum. The constraints are chosen to be invariant to therigid transformation of points (i.e., its value doesn’t change
when the points are translated and rotated). Hence, the gradient of constraint ∇C is perpendicular to the null space ofthe constraint. As long as we update the positions in the gradient space as (15), each constraint enforcement conservestranslational and rotational momentum, and consequently,
the sequence of constraint enforcements also consequently
do.
We have also tried a different approach to handle themulti-dimensional constraint. Inspired by the penalty-based
method, we create a new scalar constraint by taking square
norm of the multi-dimensional constraint (C = CT C). Thisapproach was also functional, because if the constraint iszero, the square norm becomes also zero. However, since the5. Constraints for elastic rod simulation
Based on the previous sections, we construct constraints
suitable for elastic rod simulation. Please refer to the accompanying supplemental document for differentiations of these
constraint functions.
Edge constraints. We consider three constraints for eachedge. The first one is an edge length constraint
CeL (pe−1 , pe ) = |pe−1 − pe | − L̄e ,(18)

where L̄e is the initial length of edge e. This constraint keeps
the length of the edge similar to the rest length. The second
constraint enforces that the ghost point lies on the edge’s
perpendicular bisectorgg

CeP (pe−1 , pe , pe ) = (pe − pe−1/2 )T (pe − pe−1 ),(19)

where pe−1/2 = (pe−1 + pe )/2 is the midpoint of the edge e.
g
The third constraint enforces a constant distance L̄e between
the ghost point and the edgeggg

CeD (pe−1 , pe , pe ) = pe−1/2 − pe − L̄e .

(20)

The rods’ diameter change is usually very small, hence weenforce this constraint to prevent it.
Anisotropic Bending-twisting constraint. With thesethree constraints, CeL , CeP , and CeD , we can define coordinates
at the center of edge. Next, we define bending and twisting
constraints based on the material frames defined on the orientation element. The constraint tries to minimize the difference of modified Darboux vector (defined in (10)) between
rest shape and current shape.
We specify a rod’s anisotropic stiffness via parametersαi ∈ [0, 1] (i = 1, 2, 3). More specifically, α1 denotes theresistance to bending around d1 (the direction from midpoint to ghost point), α2 denotes the resistance to bending around d2 , and α3 denotes twisting stiffness. To handle
this anisotropic stiffness, we define the goal Darboux vectorΩGoal= (1 − αi )Ωi + αi Ω̄i , with αi ∈ [0, 1]. This is simply
ian interpolation between the current Darboux vector and thatof the rest shape. The bending and twisting constraint can be
written as:Ω1 − ΩGoal
1g g
BT . (21)Ce (pe−1 , pe , pe+1 , pe , pe+1 ) =  Ω2 − ΩGoal
2
GoalΩ3 − Ω3
Note that if all the αi are zero (zero stiffness), this constraint
C takes value zero and has no effect on the constraint enforcement. Similarly, if all αi = 1, the constraint tries tomatch the Darboux vector to the one from the rest shape.
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

7

N.Umetani et al. / Position-based Elastic Rods

5.1. Endpoint constraints
So far, we have described constraints for bending and twisting deformation of rods. With a few extensions to this formulation, we can constrain the rod’s orientation at its endpoints.
We consider three cases where a rod is attached to a frame,
to a triangle, and to a rigid body.
initialdeformedframe rotation

Figure 5: Configurations of frames on a midpoint of an edge
D1 and on a frame E. The rod tip’s orientation is firmly constrained to the frame and it moves together with frame.
Frame attachment. Assume a rod’s end point p0 is attached to a frame E. The frame attachment can be interpreted
as the rod being attached to an infinitely heavy rigid body.
More specifically, this attachment is one-way coupling; theframe rotation is given and is not affected by the rod. Asgdescribed in § 3.2, the three points {p0 , p1 , p1 } adjacent to
the end make a frame D1 . We assume the rod’s tip direction
is fixed to the frame. Hence, when the frame rotates fromĒ to E, the direction of the tip rotates from D̄1 to EĒT D̄1 .
In the deformed state, the rod’s frame changes from EĒT D̄1
to D1 in the length of 0.5l¯1 . We compute coordinates of themodified Darboux vector Ωi between the frame E and D1 as
vect(P)i4
i = 1, 2, 3,(22)Ωi (s) =¯
0.5l1 1 + trP
where P = (EĒT D̄1 )T D1 , which encodes the amount of
twisting at the tip in the material coordinate of the rod. Notethat Ωi is zero in the initial rest configuration. Thus, the constraint for keeping the relative orientation between the triangle and rod identical can be written asg

COf (p0 , p1 , p1 ) = (α1 Ω1 , α2 Ω2 , α3 Ω3 )T .

The position where the end point p0 is attached to a triangle can be written as qt = w0 pt0 + w1 pt1 + w2 pt2 , where
w0 , w1 , w2 are the barycentric coordinates for p0 in the triangle in its initial configuration. By plugging in Et and qt
into the frame attachment constraint, we can convert it to atriangle attachment constraint. Note that as this is two-way
coupling constraint; three triangle vertices {pt0 , pt1 , pt2 } andgthree rod vertices {p0 , p1 , p1 } are involved.
Rigid body attachment. With the triangle constraint, we
can attach the end of a rod to an object. To represent networks of rods, such as trees and nets, we put a small rigid
body at the joint where multiple rod join together. Each rod
is connected to this rigid body to transmit torque between
rods. We represent the frame of this rigid body attachedat the tip of the rod using two additional ghost points. Atgg
the endpoint of the rod p0 , we put ghost points {p01 , p02 },g2
g1
so that the two vectors p0 − p0 and p0 − p0 are perpendicular. With this ghost point, the frame of the rigid body
Er = {e1r , e2r , e3r } at the end of a rod can be written as
 1
g1
g1

 er = (p0 − p0 ) / p0 − p0 ,
e3r




e2r

g

g

=

e1r × (p02 − p0 ) / e1r × (p02 − p0 ) ,=e3r × e1r .

(26)

By plugging Er into the frame attachment constraint, we canattach a rigid body at the tip of the rod. The orthogonalityggand length constraints for p01 − p0 and p02 − p0 are enforced
as we did for ghost point on edges in (19) and (20).

(23)

where αi are the same parameter as the rod’s stiffness. Notethat to represent the initial rest configuration, a rotation matrix ĒT D̄1 is stored for the tip constraints, while three local
coordinates of the Darboux vector Ω̄i are stored for twisting
and bending constraints. The position constraint for a tip ofa rod p0 can be simply written as
CTf (pt0 ) = |q − p0 | ,the triangle’s movement is influenced by the rod. We assume the rod tip edge is attached to a triangle t with vertices{pt0 , pt1 , pt2 } as shown in Fig. 6. We define coordinate bases
Et = {et1 , et2 , et3 } as
 1
t
t
t
t

 et = (p1 − p0 ) / p1 − p0 ,et3 = et1 × (pt2 − pt0 ) / et1 × (pt2 − pt0 ) ,
(25)

 2
3
1
et = et × et .

(24)

where the q ∈ R3 is the position where the tip is attached.
Triangle attachment. The attachment to a triangle can behandled similarly to the case of attachment to a frame. Here,we consider two ways coupling between rod and triangle;
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

Figure 6: Left: configuration of points and frames in a triangle attachment. Upper-right: configuration of points in asmall rigid-body attachment. Lower-right: Rigid body attachment is used for a network of rods.
6. Implementation details
Resolving order-induced instability. In the FEM approach, multiple constraints are simultaneously enforced by

8

N.Umetani et al. / Position-based Elastic Rods

solving a linear system. In contrast, position-based frameworks sequentially enforce constraints one-by-one, and as aresult, the simulation can be affected by the order of application. Since a rod is one dimensional, a natural ordering isnaïve sequential order (i.e., enforcing constraints from one
end to the other). However, we found that with this ordering
the rod vibrates on its own, even without any external forces
and starting the simulation from the rest shape.
In position-based frameworks, each enforcement of a constraint takes the form of internal force. Enforcing the rod’s
twist constraint in sequential order is similar to an elastic
wave propagating inside the rod. The rod has many waving modes where the elastic wave goes back and forth in
the strand and the sequential constraint enforcement can excite such vibrations. Note that this artifact related to the constraint order is specific to rod simulation, because in solids
or shells the order is generally randomly distributed.
3

1
10

5
8

7
6

9
4

2

Figure 7: Bilateral interleaving ordering.
Our solution to prevent order-dependent instability is
to utilize a bidirectional interleaving order to apply constraints. As shown in Fig. 7, we interleave the order for both
directions of the rod. Applying constraints from two directions at the same time cancels the vibration modes, producing stable results without excessive energy gain.
Applying constraints in a random order also gives a goodresult that is stable without energy gain. However, we observed a tiny flicker in the high-torsional strain arising fromrandom permutations, especially for curved rods. Generally,
the interleaving permutation gives better results.
Initial position of ghost points. The ghost points on edges
are placed so that the frame axes d1 , d2 correspond to theprincipal axis of the cross section. For example, if we model
bending of a thin planar rod, such as a ribbon, one of the axes
is chosen to align with the normal direction of the plate. This
will give the expected anisotropic behavior, where the rod is
stiff in one direction and soft in the other.
In the case of an isotropic rod (i.e., a rod with a circularcross section), the material frame can be chosen arbitrarily.
But ideally, the initial frame should be chosen such that the
Darboux vector becomes as small as possible. If the frame
rotation between two edges is 180◦ , the discrete Darboux
vector become singular. Hence, we initialize the initial ghost
position using parallel transport [BWR∗ 08], so that the material frames become twist-free. If there are loops in the rod
topology, we relax the material frame after parallel transport,so that there are no excessive jumps between frames. For thisrelaxation, we average two material frames at the two edgesof an orientation element iteratively for every orientation element inside the rod.
Gravity on ghost points. We assume that material is distributed uniformly around an edge. Hence, the center ofgravity of an edge is located at its mid-point and the gravity
force should be applied there. However, in our discretization,
the ghost point is located off the axis of an edge, and so weneed to cancel the spurious rotational moment that would be
added to the edge.
We address this problem with a simple heuristic. We
tweak the amount of gravitational force on a ghost pointsuch that the gravitational acceleration becomes equal tothat of the mid point. We compute acceleration at the midpoint by forward differentiation of its velocities: am = (vtm −t
vt−∆tm )/∆t. Here vm is the velocity of the mid-point at time t,
computed as the average velocity of the end points. After enforcing constraints and updating velocities, we compute theratio r = (am · g)/ |g|2 , which indicates how much the midpoint was accelerated by gravity. Finally, we re-distribute thegravitational forces of the ghost points to the edge points
such that the total gravitational force applied to the ghost
point is in the ratio r against gravity. Please refer to the accompanying supplemental document for more detailed algorithm and the effect of our gravity force modification compared to other approaches.
7. Results
To demonstrate our twist model, we shows an elastic rod
pendulum in Fig. 8-left. The rod is connected to a fixed
frame on one side (the frame attachment), and is connected
to a cube modeled as a triangle mesh on another side (the triangle attachment). After the cube hits another object, the rod
is twisted and then recovers. The twist resistance is achievedby the constraints in equations (21) and (23), which are
based on Cosserat theory. Fig. 8-right demonstrates stability and resistance to flipping with large twist angles. Pleaserefer to the accompanying video for more detail.

Figure 8: Left: elastic pendulum. Right: rod under large
twist.
In Fig. 9, inspired by [BWR∗ 08], we create a ring andapply twist to the ring. As expected, the ring buckles andtwists into a figure-8 pattern, similar to the results shownin [BWR∗ 08]. In this experiment we also tested anisotropic
rods, where the cross section is significantly larger in one
direction than the other (e.g., a ribbon). As twist is applied tosubmitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

9

N.Umetani et al. / Position-based Elastic Rods

the ring-shaped ribbon, it first buckles similar to the isotropic
ring, but then because it can bend only in one direction, it
shows a different deformation result.

twist=0

twist=0.5π

twist=0.8π

twist=1.2π

Figure 9: Twisting ring with (upper) isotropic cross-sectionand (lower) anisotoropic cross-section.
Comparison with previous work. We compare our rod formulation with two previous PBD rod simulation techniques.
The first is oriented particles [MC11], where the torsion resistance is modeled by averaging neighboring oriented particles, and boundary twist is given by explicitly specifying rotation matrices at the end points. The second is Chain ShapeMatching (CSM) [RKM∗ 12], which applies a twisting force
as an external force, on the top of a shape matching deformation. Using these three methods, we simulate twist deformations of a straight and spring rod (discretized with 30
edges) with both ends clamped. In our PBD-based rod, twist
and bending constraints are formulated using two neighboring edges. To compare fairly, we choose two neighboringedges as a cluster of shape matching operations in our CSM
and oriented particle implementations. We apply the same
number of constraint enforcement iterations (ten iterations)for these three methods. Since we wish to compare the deformations originating purely from elasticity, we did not include self-contact and damping in the comparisons.
Fig. 10 shows the comparison result. The oriented particles method gives stable deformation but it lacks a twisting effect because the method doesn’t model coupling of
twisting and bending. The CSM shows twisting effects but
produces unstable deformations and doesn’t converge into astatic rest shape. This is because the twisting force is applied
as an external force, not as a constraint. Generally, the external forces in shape matching work only for simple potential
forces such as gravity, and not for elastic forces. Furthermore, CSM is significantly softer than other models becauseof the singularity in the momentum matrix, which is solvedin the oriented particle model by introducing rotational inertia. Our rod deforms stably and plausibly in both examples.
We also note that our PBD approach is faster than previous
shape matching approaches because PBD doesn’t require thecostly polar decompositions.
Robustness. Fig. 11 demonstrates the robustness of our rod
model. Starting with an configuration with random vertex
placement, the simulation quickly recovers, and within asmall number of frames, has returned to the rest shape.
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)

oriented particles

chain shape matching

our method

stableno-twisting effect
4.88ms/frame

unstabletwisting effect4.96 ms/frame

stabletwisting effect1.06 ms/frame

Figure 10: Twisting (upper) straight and (lower) spring rod
simulated with different techniques. The bottom row gives
computational time per frame for each method.

Fig. 1 and 12 demonstrate some complex examples created with our rod, implemented in a widely-used visual
effects framework. These examples include large numbersof rods interacting and colliding with each other. Please
see the accompanying video to observe the stability of ourmodel in these complex situations. In the tentacle-ball example (Fig. 1), the sticky behavior is achieved with large
friction forces. The core ball itself uses a dynamic positionbased cloth model, which is coupled with our rods using
the triangle attachment constraint. In the wire-mesh example (Fig. 12), the wire mesh is created by weaving rods inan over/under configuration. With our bending and twisting
constraint, the mesh is stiff enough to be prevented from separating or unraveling, maintaining the woven configuration.

0th frame

20th frame

40th frame

100th frame

Figure 11: Starting from random initial configuration, therod quickly converges to the original shape stably.
Interactive hair design. Since our PBD simulation is fastand robust, it can be used inside interactive modeling tools.
We implemented our PBD model with a simple collisionmodel in a 3D modeling tool, to allow the user to create
hair volumes. Because our model is robust, we can allow theuser to interactively manipulate rest shape of the hair, suchas length and twisting of strands.
In addition to designing hair for virtual applications,this tool is useful to create solid hair for 3D printing purposes. Fig. 13 shows a simple example where we designed
a hairstyle for the bunny model, and the 3D-printed result.
Please see the video for more examples, including a demon-

10

N.Umetani et al. / Position-based Elastic Rods
[BAC∗ 06] B ERTAILS F., AUDOLY B., C ANI M.-P., Q UERLEUX
B., L EROY F., L ÉVÊQUE J.-L.: Super-helices for predicting the
dynamics of natural hair. ACM TOG 25, 3 (2006). 2, 4
[BAV∗ 10] B ERGOU M., AUDOLY B., VOUGA E., WARDETZKY
M., G RINSPUN E.: Discrete viscous threads. ACM TOG 29, 4(2010). 2

Figure 12: A wire mesh modeled by our position-based elastic rods. The wire mesh goes through a ring pushed by asphere. All sequence was solved robust and stably.

[BBB07] BATTY C., B ERTAILS F., B RIDSON R.: A fast variational framework for accurate solid-fluid coupling. ACM TOG
26, 3 (2007). 5[BMOT13] B ENDER J., M ÜLLER M., OTADUY M. A.,
T ESCHNER M.: Position-based methods for the simulation ofsolid objects in computer graphics. In EG STARs (2013). 3

stration showing over 200 hairs, each is modeled with 20edges, being simulated at interactive rates.

[BWR∗ 08] B ERGOU M., WARDETZKY M., ROBINSON S., AU DOLY B., G RINSPUN E.: Discrete elastic rods. ACM TOG 27, 3(2008). 1, 2, 3, 5, 8

8. Discussion and Limitations[CBD13] C ASATI R., B ERTAILS -D ESCOUBES F.: Super space
clothoids. ACM TOG 32, 4 (2013). 2

We present position based elastic rods, which can simulate arod’s bending and twisting deformation efficiently using theframework of position based dynamics.
We remove the order-induced instability of constraint enforcement for rods by introducing a special ordering. Furthermore, we remove artificial rotational momentum induced by the gravity on ghost points. While these methods
works perfectly for various cases, there are no theoretical
proofs. Guaranteeing stability is one of our future works.
Our variation interpretation of constraint enforcement enables a handling multi-dimensional constraint that is specifically useful to modeling deformation of anisotropic materials. We are planning to apply this concept for anisotropic
solid and shell deformation simulation. Furthermore, handling wider variety of constraints such as ones shown in[SJLP11] is an interesting future work.
Although our bending and twisting constraints prevent
over 180◦ twisting via the Rodriguez parameter, under hightorsional circumstances the rod may still exceed this angle
limit. This is because we apply partial constraints sequentially. In future work, we will explore additional constraints
to prevent this flip completely.

3D printing[CC07] C OSSERAT E., C OSSERAT F.: Sur la statique de la ligne
déformable. Académie des Sciences 145 (1907), 1409–1412. 2[GJS01] G OLDSTEIN H., J R . C. P. P., S AFKO J. L.: ClassicalMechanics (3rd Edition), 3 ed. Addison-Wesley, 6 2001. 5[Jak01] JAKOBSEN T.: Advanced character physics. In Game
Developers Conference (2001). 3[KJM10] K ALDOR J. M., JAMES D. L., M ARSCHNER S.: Efficient yarn-based cloth with adaptive contact linearization. ACM
TOG 29, 4 (2010). 2[MC11] M ÜLLER M., C HENTANEZ N.: Solid simulation withoriented particles. ACM TOG 30, 4 (2011). 3, 9
[MHHR07] M ÜLLER M., H EIDELBERGER B., H ENNIX M.,
R ATCLIFF J.: Position based dynamics. Visual Communication
and Image Representation 18, 2 (2007). 1, 3, 5
[MHTG05] M ÜLLER M., H EIDELBERGER B., T ESCHNER M.,
G ROSS M.: Meshless deformations based on shape matching.
ACM TOG 24, 3 (2005). 3[MKC12] M ÜLLER M., K IM T.-Y., C HENTANEZ N.: Fast simulation of inextensible hair and fur. In VRIPHYS (2012). 3[MLSS94] M URRAY R. M., L I Z., S ASTRY S. S., S ASTRY S. S.:
A Mathematical Introduction to Robotic Manipulation, 2nd ed.
CRC Press, 3 1994. 5[MM13] M ACKLIN M., M ÜLLER M.:
ACM TOG 32, 4 (2013). 3

Position based fluids.

[MSW∗ 09] M C A DAMS A., S ELLE A., WARD K., S IFAKIS E.,
T ERAN J.: Detail preserving continuum simulation of straight
hair. ACM TOG 28, 3 (2009). 2[Pai02] PAI D. K.: Strands: Interactive simulation of thin solids
using cosserat models. CGF 21, 3 (2002). 2[RKM∗ 12] RUNGJIRATANANONW.,
K ANAMORIY.,M ETAAPHANON N., BANDO Y., C HEN B.-Y., N ISHITA
T.: Animating strings with twisting, tearing and flicking effects.
Computer Animation and Virtual Worlds 23, 2 (2012). 3, 9

Figure 13: Left: Designed bunny hair. Right: 3D print.

[SJLP11] S UEDA S., J ONES G. L., L EVIN D. I. W., PAI D. K.:
Large-scale dynamic simulation of highly constrained strands.
ACM TOG 30, 4 (2011). 10

Acknowledgement. We thank Duncan Brinsmead for helping us preparing animations.

[SLF08] S ELLE A., L ENTINE M., F EDKIW R.: A mass springmodel for hair simulation. ACM TOG 27, 3 (2008). 2

References[ST07] S PILLMANN J., T ESCHNER M.: CoRdE: Cosserat rod
elements for the dynamic simulation of one-dimensional elastic
objects. In ACM SCA (2007). 1, 2, 4

[Ant05] A NTMAN S. S.: Nonlinear problems of elasticity,
vol. 107. Springer, 2005. 3[Sta09] S TAM J.: Nucleus: Towards a unified dynamics solver forcomputer graphics. In CAD/Graphics 2009 (2009), pp. 1–11. 3
submitted to Eurographics/ ACM SIGGRAPH Symposium on Computer Animation (2014)