Developable surfaces are ubiquitous in our daily life. Although their continuous properties have been well understood , 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  or non-conforming triangle meshes . 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  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 , 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.
2 Related Works
Rod Modelling The theory of elasticity for 1D rod has been established in . Later on, various discretization scheme for this model has been developed and applied in robotics , virtual surgery  and computer animation [13, 8, 14, 15, 9].
These discretization schemes fall in two categories:  and  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  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.
Finite Element Shells Finite element method is another promising alternative for modelling thin shells, see  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  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 . 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 , 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  by enforcing the length constraints on a non-conforming mesh. However,  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,  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 . Like , developability in their method depends on a set of nonlinear constraints to be satisfied. Finally, our representation is most closely related to , 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,  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 , 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.
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 , 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 vectorand 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. :
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 . The mean curvature measure of is which follows from the tube theory  and the mean curvature is then approximated in a area averaged manner. Although  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 , 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.
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.
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  as our underlying solver, which is wrapped into an Augmented-Lagrangian framework  to handle nonlinear 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.
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.
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  and globally coupled methods such as non-rigid impact zone . 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 , 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.
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.
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.
|Helix Unrolling (h=0.01)|
|Helix Unrolling (h=0.05)|
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.
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.
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.
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.
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 .
-  M. P. Do Carmo and M. P. Do Carmo, Differential geometry of curves and surfaces. Prentice-hall Englewood Cliffs, 1976, vol. 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: http://dx.doi.org/10.1111/j.1467-8659.2007.01059.x
-  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: http://dx.doi.org/10.1111/j.1467-8659.2012.03162.x
-  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.
-  E. English and R. Bridson, “Animating developable surfaces using nonconforming elements,” ACM Transaction on Graphics, 2008.
-  ——, “Animating developable surfaces using nonconforming elements,” in ACM Transactions on Graphics (TOG), vol. 27, no. 3. ACM, 2008, p. 66.
-  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.
-  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.
-  R. Casati and F. Bertails-Descoubes, “Super space clothoids,” ACM Transactions on Graphics (TOG), vol. 32, no. 4, p. 48, 2013.
-  E. Cosserat and F. Cosserat, “Théorie des corps déformables,” Paris, 1909.
-  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.
-  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: http://dx.doi.org/10.1111/1467-8659.00594
-  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.
-  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.
-  F. Bertails, “Linear time super-helices,” in Computer graphics forum, vol. 28, no. 2. Wiley Online Library, 2009, pp. 417–426.
-  P. G. Ciarlet, Theory of shells. Elsevier, 2000.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  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.
-  J. Nocedal and S. J. Wright, “Numerical optimization, second edition,” Numerical optimization, pp. 497–528, 2006.
-  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.
-  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.
-  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:
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 .