Modelling Developable Ribbons Using Ruling Bending Coordinates

03/13/2016 ∙ by Zherong Pan, et al. ∙ 0

This paper presents a new method for modelling the dynamic behaviour of developable ribbons, two dimensional strips with much smaller width than length. Instead of approximating such surface with a general triangle mesh, we characterize it by a set of creases and bending angles across them. This representation allows the developability to be satisfied everywhere while still leaves enough degree of freedom to represent salient global deformation. We show how the potential and kinetic energies can be properly discretized in this configuration space and time integrated in a fully implicit manner. The result is a dynamic simulator with several desirable features: We can model non-trivial deformation using much fewer elements than conventional FEM method. It is stable under extreme deformation, external force or large timestep size. And we can readily handle various user constraints in Euclidean space.



There are no comments yet.


page 2

page 3

page 4

page 5

page 6

page 7

page 8

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Developable surfaces are ubiquitous in our daily life. Although their continuous properties have been well understood [1], their accurate modelling and discretization is still an open problem. Recently, methods have been proposed in [2, 3] to model these surfaces statically. In this paper, we takes a step further to model the dynamic properties of developable ribbons, a special type of developable surfaces that can be isometrically mapped to two-dimensional strips with much smaller width (latitude dimension) than length (longitude dimension). Developable ribbons have seen a lot of applications for modelling stylish hairs, satin bows or films.

FEM solver is clearly a competitive solution for modelling developable surfaces, using either conforming [4] or non-conforming triangle meshes [5]. However, none of these methods are geometrically accurate in that their configuration space is not a subset of the true developable shape space. As a result, on conforming meshes large stiffness energies have to be introduced to limit the stretch, which in turn leads to the locking phenomena. On the other hand, the hard length constraints in [6] on non-conforming meshes greatly limit the timestep size. Moreover, the reconstructed conforming meshes are again not exactly developable and usually suffer from noisy perturbation.

Key to our dynamic ribbon simulator is a novel configuration space that parameterizes a subset of the developable shape space which is large enough to cover most non-trivial deformations. Specifically, we describe the shape of ribbon by a set of creases along the centerline and bending angles across them, see fig:Parameter. This method is in direct contrary to previous reduced models such as [7], where various constrains are introduced to pull the shape towards the true shape space. These constraints usually lead to stability issue or locking phenomena. Instead, our novel representation guarantees that the ribbon can be isometrically mapped to material space. As a result, no additional constraints are needed, making our method stable under large external force or timestep size. Moreover, the whole timestepping scheme can be formulated as a single optimization, which greatly simplifies implementation. We noticed that a similar idea has been exploited in [8, 9] for modelling helical rod.

Under this configuration space, we present a proper discretization of the kinetic and potential energies. Our discretization scheme bears several desirable features: First, material space remeshing and world space deformations are modelled uniformly; Energy gradients can be analytically evaluated allowing quasi-newton method to converge efficiently; The vertex positions are reintroduced as auxiliary variables so that conventional collision handlers can be trivially port to our new formulation. In conclusion, our contributions can be summarized as follows:

  • A novel parameterization of salient global deformations of developable ribbon.

  • A discrete timestepping scheme that can be efficiently integrated in a fully implicit manner.

  • An optimization-based framework for multi-ribbon simulation allowing flexible user constraints and collision resolution.

The rest of the paper is organized as follows. After briefly reviewing the related works, we first describe the transfer function between and our configuration space in sec:Transfer. We then present our discretization scheme for the kinetic and potential energies in sec:Discrete. Finally, in sec:Optimization, we go into some implementation details of our optimization strategy, constrain and collision handling before we conclusion our discussion.

Fig. 1: Illustration of our configuration space. The centerline of a ribbon is evenly segmented into elements separated by creases. The angles between these creases and the centerline are and the final shape of the ribbon in is reconstructed by bending along these creases. The corresponding bending angles are . Given these parameters and , the two ends of each crease, and , can then be derived analytically.

2 Related Works

Rod Modelling The theory of elasticity for 1D rod has been established in [10]. Later on, various discretization scheme for this model has been developed and applied in robotics [11], virtual surgery [12] and computer animation [13, 8, 14, 15, 9].

These discretization schemes fall in two categories: [13] and [14] adopted a hybrid representation. In their method, the centerline is discretized in with a frame attached to each segment. This configuration space is not a subset of the true developable shape space so that additional constrains are needed for inextensibility and consistency between the frames and the centerline. Our method is more closely related to [8, 14] and [9] where the configuration space is parameterized solely by the differentials of positions in . These methods share the advantage that no extra constraints are needed. But a reconstruction procedure is required to recover variables. Despite these similarities, none of them can be directly used to model developable ribbons because their configuration space has only a small intersection with the developable shape space. For example, one may extend a rod along its binormal directions to get a ribbon-like surface but its deformation away from the centerline is not isometric as illustrated in fig:compRod.

Fig. 2: Compared with our method (top), extending a twisted rod (red) along binormal directions (top) doesn’t give isometric deformation (area error=)

Finite Element Shells Finite element method is another promising alternative for modelling thin shells, see [16] for a description of their continuous model. The discrete counterpart has been introduced into the graphics community in [17, 18]. But these methods model elastic, instead of developable shells. Later, it is shown in [19] that isometric deformation can be approximated on a conforming mesh by enforcing hard length constraints in a post-projection step. Moreover, this method has the good property that bending energies become quadratic [20]. Although in this work we adopt the dihedral angle based formulates following [17, 18], our bending energies are also quadratic since we treat the bending angles as our generalized coordinates.

However, although a developable surface can be approximated using [19], a triangular conforming mesh has insufficient degrees of freedom to cover the developable shape space, leading to the so-called locking phenomena. This problem is resolved in [6] by enforcing the length constraints on a non-conforming mesh. However, [6] suffers from noisy vertex perturbation when a conforming mesh is reconstruction for rendering and collision resolution. Also, the fast-projection involved in [19, 6] greatly limits the timestep size and stability. Compared with these methods, our formulation allows the same or even higher accuracy with much less elements since no discretization along the latitude direction is needed.

Developable Surface Modelling Our method is also closed related to previous efforts towards static developable surfaces modelling. Among these works, [21] describes a rectifying developable surface by the envelop of its rectifying planes. However, their method cannot be directly used for dynamic modelling since the centerline is represented by an Bèzier curve, on which inextensible and other consistency constraints are hard to formulate. Developable surfaces have also been known in the community of architectural design as PQ meshes [22]. Like [13], developability in their method depends on a set of nonlinear constraints to be satisfied. Finally, our representation is most closely related to [3], where a developable surface is characterized explicitly by creases and bending angles. But we emphasize that, since the correct dynamic behaviour of a developable ribbon heavily depends on the material space crease direction changes, [3] cannot be directly used to this end because their crease directions are hard to parameterize.

3 Ruling-Bending Coordinates

One unique property of a developable surface is that the Gaussian Curvature is zero everywhere. As a result, each point on is attached to a ruling line or crease along which normal is constant and the ribbon is bended by angle . As is noted in [3], the world space shape of can be characterized by these crease directions and bending angles up to rigid transformation. In our formulation, these two sets of parameters define our configuration space. However, the crease directions for a general developable surface is hard to parameterize. Fortunately, we have observed that for with much larger longitude then latitude dimension, a large subset of salient deformations can be modelled with only creases that pass through the ribbon centerline, see fig:Deformation for an illustration. We can then parameterize our crease direction by , where is the angle between the crease and the centerline.

Fig. 3: Non-trivial global deformations (brown) can be modelled using only creases that cross the centerline. However, some local deformations (red) is excluded from our shape space.

Specifically, given with longitude dimension and latitude dimension , we first segment its centerline into elements . Then, between any two consecutive elements we introduce creases with directions and bending angles . Since we don’t allow singularity points, any two consecutive creases cannot intersection, leading to a set of crease constraints:


where and . In practice, we set to avoid degenerate triangles in the reconstructed mesh for collision handling. When these constraints are satisfied, our configuration space is thus parameterized by .

Like [8], a reconstruction procedure is needed to recover positions of bottom rim vertices and top rim vertices . The material space positions of these vertices are:

where we need to introduce boundary crease directions . Here we used homogeneous coordinates for convenience. Their world space positions and can then be derived by applying a series of crease transformations:


where is a rigid rotation by angle along crease . In this equation, we again need to introduce boundary value . For , if the first segment of the ribbon is fixed, we simply have as well. While if the ribbon is attached to a floating frame, is a global rigid transformation:

parameterized by rotation vector

and translation . If this is the case, our configuration space is parameterized by the set of variables . For numerical optimizaiton, our dynamic simulator heavily depends on an analytical formula for . We leave their derivations to appen:DerivRecon. can be found following the same procedure. Besides, when extra torsional forces are applied on the ribbon, we formulate them as additional constraints on the normal directions in world space. For element , its normal direction can be calculate by:


where . Its derivative can be found similarly.

4 Discrete Equation of Motion

The above configuration space naturally encodes the shape of a globally deformed developable ribbon. To find its motion in temporal domain, we have to discretize the equation of motion. To simplify our presentation, we start from temporal discretization using the Implicit Euler method. It has a simple variational form, which has been exploited in e.g. [23]:

where is position vector assembled from and . Here, denotes the internal or external potential energy terms. Since the ribbon mesh will deform in material space, the mass matrix derived from conventional FEM method is dependent on the crease direction . See appen:MassMat for more details.

In this work, since the configuration space is a subset of the true shape space, no stiffness energies are need to limit stretch and we are left with bending energies to be considered. Since we have the bending angles as an independent variable, bending energies based on dihedral angles [4, 18] become quadratic in in our case. Specifically, for each crease between , we introduce:

where is a diamond element between with area . This formula is derived in a similar way to [17]. The mean curvature measure of is which follows from the tube theory [24] and the mean curvature is then approximated in a area averaged manner. Although [3] adopted a more accurate form, taking the change of mean curvature along ruling into consideration, the difference is insignificant compared with our simplified form.

Other potential terms are discrete version of external forces or soft constraints. For example, the gravitational potential energies are simply: . In addition to these terms, we add a small regularization to resolve the ambiguity of crease directions for a flat ribbon,preferring ruling directions orthogonal to the centerline. The final is:

where is the bending stiffness coefficient and in all our examples. We also added an additional term for user controllability.

On the other hand, our compact configuration space poses great challenge on spatial discretization of the kinetic term. This is because encodes material space remeshing (by changing ) and world space deformation in a uniform manner. As a result, positions and may not correspond to the same points and in material space and cannot be subtracted directly. A common practice here is to temporarily fix to ensure material space consistency and perform remeshing regularly. Although this strategy has been successfully adopted in conventional FEM shell solver such as [25], it fails to work with our method because our configuration space is so compact that bending angles along cannot cover a large enough subset of the true developable shape space, leading to severe locking artifact, see fig:Locking.

Fig. 4: A torsional force is applied to a ribbon with one end fixed. Due to our compact configuration space, bending angles alone cannot cover large enough subset of the developable shape space (top). The naive method of updating crease direction once every frames (middle) still suffer from large deviation from the ground truth (bottom).

We thus need to introduce a resampling operation satisfying:


where is assumed to be piecewise linear in its first parameter. In fact, since all vertices lie on the top or bottom rim of the ribbon,

is a simple 1D-interpolation. We can then define a valid subtraction as:

, where means apply resampling on each block of . Note that in this way we have essentially encoded remeshing and deformation in a single optimization, and our final form of optimization becomes:

We want to emphasize that positions here is not required to be reconstructed from our generalized coordinates, allowing the ribbon to temporarily deviate from the developable shape space. This flexibility enables conventional collision handlers to be easily integrated into our framework.

Our solver for this nonlinear optimization is detailed in sec:Optimization, which heavily relies on an analytical formula for the energy gradient. A naive way of gradient evaluation may follow from appen:DerivRecon and chain rule. But we show in appen:AdjointMode that it could be greatly accelerated by evaluating in an adjoint mode.

5 The Solver Framework

In this section, we present our full-featured multi-ribbon solver framework, including various user constraints or external forces handling and collision resolution. Alg:Outline provides an outline of our two-step pipeline. In the first substep, the optimization problem is solved for each ribbon in parallel to predicate a desired new configuration where the vertex positions are reconstructed. The collision handler then finds a corrected collision free positions in the second substep, temporarily leaving the developable shape space.

1:a mesh with vertices
2:a set of user constraints
3:for each ribbon  do in parallel
4:     =optimize(,)
5:     =reconstruct()
6:end for
Algorithm 1 one timestep of ribbon solver

For the optimization in our first substep, since the analytical Hessian matrix of eq:Recon is too costly to evaluate, Newton-type solvers becomes largely unavailable. We thus choose the L-BFGS-B algorithm [26] as our underlying solver, which is wrapped into an Augmented-Lagrangian framework [27] to handle nonlinear constraints.

5.1 Constraints

Now we discuss several types of constraints supported by our framework. The most important is the non-intersecting ruling constraints. Since these constraints are large in number, we transform them into box constraints by a variable substitution as:

so that they can be handled by the L-BFGS-B algorithm. We are thus left with only two general linear constraints:

to be handled by the Augmented-Lagrangian framework.

Another common constraint is the loop constraint that requires the two ends of a ribbon to be connected:

for the orientable case or:

for the non-orientable case. These constraints are useful for modelling a ribbon chain or the Mobius band, see fig:Loop.

Fig. 5: Our algorithm can robustly handle nonlinear loop constraints. In this example, an non-orientable loop constraint is assigned to each of the 9 ribbons, forming a Mobius chain.

Finally, a lot of interesting deformations are resulted from torsional forces which is introduced into our framework as additional normal guiding energies: . Other kinds of common user constraints can be added as in conventional FEM methods. Now we can summarize our optimization substep in alg:Optimize.

2:for  do
4:     initial guess
7:     if  then
8:         return
9:     else
11:     end if
12:end for
Algorithm 2 optimize(,)

5.2 Collision Resolution

One advantage of our formulation is that existing collision detection and resolution methods can be trivially plugged into our framework as a post processor. After a new configuration is returned by our optimizer, a triangle mesh with vertices is reconstructed and passed to the collision handler, which in turn finds a closest collision free mesh with vertices: . Although this mesh may not lie in the developable shape space, its distance to our configuration space is very close in our experiments. Moreover, this developability error will not accumulate because a valid configuration is always recovered by our first substep at next frame.

To work the best with our method, one need to be careful in choosing of underlying continuous collision handler. There are generally two methods for resolving a large amount of continuous collisions: local methods based on randomized impulses [28] and globally coupled methods such as non-rigid impact zone [29]. Since the closeness to a valid configuration is important in our case, we choose to use globally coupled handlers. Specifically, we solve a quadratic programming of the following form to resolve a set of collision constraints are returned by the detector:

where is some closeness measure. In the original work [29], is simply . This measure leads to diagonal Hessian so that the QP problem can be solved efficiently in its dual form for a very large mesh, especially when the number of constraints are much smaller than number of vertices. Although this measure totally ignores the stiffness between vertices, this choice is appropriate for most cloth animation setups where a small timestep size is used.

But in our case, the situation is reversed. Due to our compact representation, the reconstructed mesh is orders of magnitude smaller for comparable results than FEM method. As a result, the number of potential constraints are usually comparable to the number of vertices. On the other hand, since we used fully implicit method, our solver is stable under large timestep. Unfortunately, such large timestep size also makes the collision force stiff. As a result, ignoring internal stiffness in would lead to large discrepancy between and . This may introduce large developability error and finally cause collision failure due to degenerated or flipped triangles which is illustrated in fig:CollStiff.

Fig. 6: A frame from the same animation as bottom fig:ribbonCollision. Compared with conventional method (left), the collision handler can be greatly stabilized with additional terms (right) under large timestep size.

Out of these considerations, we add an quadratic artificial stiffness term to . Specifically, for each edge of the reconstructed triangle mesh with vertices , we introduce additional energy terms: , where is an artificial stiffness coefficient set to in all our examples. The new QP problem:

is again solved in its dual form using the active set method, where the dual Hessian is calculated by pre-factorizing the sparse Hessian and solve the sparse right hand side . Since our mesh is rather small, the overhead of this solve is neglectable.

6 Results and Validations

The stability, accuracy and efficiency of our method is evaluated using several benchmark tests, the performance of our solver on all examples is summarized in table:performance. In a first set of tests illustrated in fig:LargeForce, we apply large torsional or dragging forces on a ribbon with or without loop constraint. Fortunately, our optimizer presents no performance degradation in both cases.

Fig. 7: Stability test on a ribbon with length= and width= segmented into elements. Top row: twisting the spherical ribbon by with strong normal constraints. Bottom row: a looped ribbon under large dragging force.

The stability is also validated in temporal domain. In this subsequent test, we compared two simulated helix unrolling sequences under different timestep size. Our solver finds a reasonably consistent result with up to , see fig:Timestep.

Fig. 8: Helix unrolling of a ribbon with length= and width= segmented into elements. Under different timestep size on the left and on the right, our solver finds consistent results.
Scene Ribbon Opt./[Coll.](sec) Inner/Outer
Extreme Torsion
Extreme Dragging
Torsion (ours)
Torsion (C-FEM) N/A
Torsion (NC-FEM) N/A
Helix Unrolling (h=0.01)
Helix Unrolling (h=0.05)
Ribbon Falling
Double Chain
Mobius Chain
TABLE I: Time cost for one step of our solver. From left to right: number of ribbon segments in longitude/latitude dimension, average time cost for optimization/collision detection if applicable and number of inner/outer iterations taken by the optimizer. All tests are done on a single desktop computer with dual Intel Xeon E5-2630 CPU and 128Gb memory.

To demonstrate our advantage over previous methods, we compared our solver with conventional FEM methods using (non-)conforming mesh. For the conforming solver, we have to introduce an artificial stiffness term which is set to . According to fig:Comparison and the accompanying video, the conforming solver would suffer from noisy spatial error due to locking phenomena and the non-conforming solver instead suffers from noisy temporal error due to the conforming reconstruction procedure and additional boundary constraints. By contrast, our method on a much smaller mesh faithfully regenerates the zigzag ruling pattern without any instability in both spatial and temporal domain.

Fig. 9: Twisting a ribbon with length= and width= segmented into longitude elements and latitude elements. From top to bottom: our method, FEM solver on a conforming mesh and FEM solver on a non-conforming mesh.

In addition, we investigated the accuracy of our kinetic energy formulation by showing the consistency with conventional FEM methods again. However, as is shown in fig:Comparison, FEM solver cannot serve as a valid groundtruth. To resolve this contradiction, we observe that our method is equivalent to an conventional FEM solver with dynamic remeshing respecting the ruling direction. In view of this, we work on an animation sequence generated using our method, and for any two consecutive frames and , we insert all vertices to get a combined mesh with material space vertex positions . From this mesh with world space vertex positions , we perform time integration using FEM method on conforming mesh to predicate a new combined mesh with world space vertex positions , which is then compared with . This essentially compares the per-frame discrepancy between our method and FEM method with remeshing factored out. The result is visualized in fig:ComparisonKinetic, which validates the accuracy of our formulation.

Fig. 10: The animation of free helix unrolling with . Visualization of discrepancy between (red) and (brown) at , from left to right. For reference, is also shown in green.

Unfortunately, such accuracy is achieved at the cost of an additional resampling operator . There is however a naive simplification to our kinetic term by lumping all the mass to the centerline, giving a simplified objective energy:

which just takes the average of velocities and positions along the ruling direction for the centerline. Here is the constant lumped mass for the th centerline vertex. Such approximation works well in some cases such as fig:ComparisonKinetic but may lead to severe artifact elsewhere. An extreme example is shown in fig:Rotate.

Fig. 11: A ribbon rotating along the centerline (no gravity). If lumped kinetic energy is used, the ribbon won’t even deform since is always achieved at the rest state (top). While our formulation correctly captures the triangular ruling pattern (middle), which finally leads to large deformation due to centrifugal force (bottom).

One major disadvantage of our algorithm is that the performance degenerates with the number of segments according to fig:perfSeg. The bottleneck of our method is the evaluation of and which is quadratic in , which can be largely removed using the adjoint method. But longer ribbon won’t affect the stability our method. In the falling ribbon example of fig:ribbonCollision, we used a long ribbon with 150 segments. In this case the overhead of optimization dominates our solver pipeline.

Fig. 12: Our method handles collision and contact robustly. Top: A long ribbon falling on the ground. Bottom: impact of two ribbon chains with 9 ribbons each.
Fig. 13: Time cost of finding the rest shape of a mobius band for ribbon of different length, using simple chain rule (red) or the adjoint method (blue) energy gradient evaluation.

The stability of collision and contact handling is illustrated in fig:ribbonCollision as well. For both of these examples we used . Under such large timestep size, our new closeness metric is an indispensable component especially for the double ribbon chain example. Due to the stiffness of collision forces, conventional non-rigid impact zone solver generates extremely distorted triangles and quickly fails after the first few segments of the chain get in contact.

7 Conclusion and Discussion

In conclusion, the paper presents an configuration space that covers a large subset of the shape space of developable ribbon. Based on this configuration space, we develop a optimization based ribbon simulator with flexible user controllability and robust collision handling. Thanks to the compactness of the configuration space, the solver is locking free compared with previous FEM based methods. This enables the solver high fidelity on a very small mesh. Moreover, it presents better stability in both temporal and spatial domain over conventional methods [17, 6].

We also noticed several drawbacks of the method. A clear problem is that our configuration parameters are densely related to vertex positions, which limits the scalability in terms of both time and memory for extremely long chain. The optimizer would also require more gradient evaluations for longer ribbon. To alleviate this problem, it is worth exploring acceleration techniques such as Newton-type optimizer with an approximate Hessian or massive parallelism. Another major drawback is that, as the width of ribbon increase, our configuration space represents a smaller subspace of the true shape space. And sometimes user may want to recover some local deformation. In these cases, locally deformed patches can be reintroduced by coupling the solver to conventional FEM methods or model reduction techniques such as [30].


  • [1] M. P. Do Carmo and M. P. Do Carmo, Differential geometry of curves and surfaces.   Prentice-hall Englewood Cliffs, 1976, vol. 2.
  • [2] P. Bo and W. Wang, “Geodesic-controlled developable surfaces for modeling paper bending,” Computer Graphics Forum, vol. 26, no. 3, pp. 365–374, 2007. [Online]. Available:
  • [3] J. Solomon, E. Vouga, M. Wardetzky, and E. Grinspun, “Flexible developable surfaces,” Computer Graphics Forum, vol. 31, no. 5, pp. 1567–1576, 2012. [Online]. Available:
  • [4] E. Grinspun, A. N. Hirani, M. Desbrun, and P. Schröder, “Discrete shells,” in Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation.   Eurographics Association, 2003, pp. 62–67.
  • [5] E. English and R. Bridson, “Animating developable surfaces using nonconforming elements,” ACM Transaction on Graphics, 2008.
  • [6] ——, “Animating developable surfaces using nonconforming elements,” in ACM Transactions on Graphics (TOG), vol. 27, no. 3.   ACM, 2008, p. 66.
  • [7] J. Spillmann and M. Teschner, “Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects,” in In Proc. ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007, pp. 63–72.
  • [8] F. Bertails, B. Audoly, M.-P. Cani, B. Querleux, F. Leroy, and J.-L. Lévêque, “Super-helices for predicting the dynamics of natural hair,” in ACM Transactions on Graphics (TOG), vol. 25, no. 3.   ACM, 2006, pp. 1180–1187.
  • [9] R. Casati and F. Bertails-Descoubes, “Super space clothoids,” ACM Transactions on Graphics (TOG), vol. 32, no. 4, p. 48, 2013.
  • [10] E. Cosserat and F. Cosserat, “Théorie des corps déformables,” Paris, 1909.
  • [11] S. Javdani, S. Tandon, J. Tang, J. F. O’Brien, and P. Abbeel, “Modeling and perception of deformable one-dimensional objects,” in Robotics and Automation (ICRA), 2011 IEEE International Conference on.   IEEE, 2011, pp. 1607–1614.
  • [12] D. K. Pai, “Strands: Interactive simulation of thin solids using cosserat models,” Computer Graphics Forum, vol. 21, no. 3, pp. 347–352, 2002. [Online]. Available:
  • [13] J. Spillmann and M. Teschner, “C o r d e: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects,” in Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation.   Eurographics Association, 2007, pp. 63–72.
  • [14] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun, “Discrete elastic rods,” in ACM Transactions on Graphics (TOG), vol. 27, no. 3.   ACM, 2008, p. 63.
  • [15] F. Bertails, “Linear time super-helices,” in Computer graphics forum, vol. 28, no. 2.   Wiley Online Library, 2009, pp. 417–426.
  • [16] P. G. Ciarlet, Theory of shells.   Elsevier, 2000.
  • [17] E. Grinspun, A. N. Hirani, M. Desbrun, and P. Schröder, “Discrete shells,” in Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation.   Eurographics Association, 2003, pp. 62–67.
  • [18] R. Bridson, S. Marino, and R. Fedkiw, “Simulation of clothing with folds and wrinkles,” in Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation.   Eurographics Association, 2003, pp. 28–36.
  • [19] R. Goldenthal, D. Harmon, R. Fattal, M. Bercovier, and E. Grinspun, “Efficient simulation of inextensible cloth,” ACM Transactions on Graphics (TOG), vol. 26, no. 3, p. 49, 2007.
  • [20] M. Wardetzky, M. Bergou, D. Harmon, D. Zorin, and E. Grinspun, “Discrete quadratic curvature energies,” Computer Aided Geometric Design, vol. 24, no. 8, pp. 499–518, 2007.
  • [21] P. Bo and W. Wang, “Geodesic-controlled developable surfaces for modeling paper bending,” in Computer Graphics Forum, vol. 26, no. 3.   Wiley Online Library, 2007, pp. 365–374.
  • [22] Y. Liu, H. Pottmann, J. Wallner, Y.-L. Yang, and W. Wang, “Geometric modeling with conical meshes and developable surfaces,” in ACM Transactions on Graphics (TOG), vol. 25, no. 3.   ACM, 2006, pp. 681–689.
  • [23] S. Martin, B. Thomaszewski, E. Grinspun, and M. Gross, “Example-based elastic materials,” in ACM Transactions on Graphics (TOG), vol. 30, no. 4.   ACM, 2011, p. 72.
  • [24] D. Cohen-Steiner and J.-M. Morvan, “Restricted delaunay triangulations and normal cycle,” in Proceedings of the nineteenth annual symposium on Computational geometry.   ACM, 2003, pp. 312–321.
  • [25] R. Narain, A. Samii, and J. F. O’Brien, “Adaptive anisotropic remeshing for cloth simulation,” ACM Transactions on Graphics (TOG), vol. 31, no. 6, p. 152, 2012.
  • [26] R. H. Byrd, P. Lu, J. Nocedal, and C. Zhu, “A limited memory algorithm for bound constrained optimization,” SIAM Journal on Scientific Computing, vol. 16, no. 5, pp. 1190–1208, 1995.
  • [27] J. Nocedal and S. J. Wright, “Numerical optimization, second edition,” Numerical optimization, pp. 497–528, 2006.
  • [28] R. Bridson, R. Fedkiw, and J. Anderson, “Robust treatment of collisions, contact and friction for cloth animation,” in ACM Transactions on Graphics (ToG), vol. 21, no. 3.   ACM, 2002, pp. 594–603.
  • [29] D. Harmon, E. Vouga, R. Tamstorf, and E. Grinspun, “Robust treatment of simultaneous collisions,” in ACM Transactions on Graphics (TOG), vol. 27, no. 3.   ACM, 2008, p. 23.
  • [30] D. Harmon and D. Zorin, “Subspace integration with local deformations,” ACM Transactions on Graphics (TOG), vol. 32, no. 4, p. 107, 2013.

Appendix A Explicit Formula for and

Since for is just rigid rotation, they have simple analytical form:

whose partial derivatives and can be found using a symbolic software. We omit these here to save space. And for the global rigid transformation, we have:

where can be found using the Rodriguez’s formula. Now can be founding using simply chain rule:

Appendix B The Mass Matrix

Without loss of generality, we use linear shape function for each quad element. In this case, element would contribute a mass block of the following form:

where .

Appendix C Adjoint Mode Gradient Evaluation

In order to evaluate , we start from the chain rule:

where the first two terms can then be evaluated in an adjoint mode by exploiting the special structure of eq:Recon and eq:ReconN, which is composed of two passes as illustrated in Alg:AdjointAlg. In this way, the algorithmic complexity is reduced from to .

3: forward pass
4:for  do
8:end for
10: backward pass
11: adjoint variable
12:for  do
19:end for
Algorithm 3 Adjoint Gradient Evaluation