An Eigenanalysis of Angle-Based Deformation Energies

Angle-based energies appear in numerous physics-based simulation models, including thin-shell bending and isotropic elastic strands. We present a generic analysis of these energies that allows us to analytically filter the negative eigenvalues of the second derivative (Hessian), which is critical for stable, implicit time integration. While these energies are usually formulated in terms of angles and positions, we propose an abstract edge stencil that succinctly parameterizes the edge deformation, and allows us to derive generic, closed-form analytical expressions for the energy eigensystems. The resultant eigenvectors have straightforward geometric interpretations. We demonstrate that our method is readily applicable to a variety of 2D and 3D angle-based elastic energies, including both cloth and strands, and is up to 7× faster than numerical eigendecomposition.


INTRODUCTION
Many elastic energies used for deformable objects are based on angles, which often exhibit complex nonlinear relationships to vertex positions, and incorporate geometric quantities in a reduced space. This makes the energies difficult to analyze, and in particular it becomes unclear when the Hessians associated with their force gradients become indefinite. This is unfortunate, because many implicit solvers [Baraff and Witkin 1998] require the underlying system to be semi-positive definite (SPD) in order to guarantee convergence.
Previous works have produced SPD systems by using Gauss-Newton approximation [Choi and Ko 2002] or by explicitly filtering the Hessian at each element [Teran et al. 2005]. Recent methods efficiently addressed this issue by constructing analytical eigendecompositions at each element.
While this approach has been applied to hyperelastic volume Smith et al. 2019] and membrane [Kim 2020] energies, no equivalent analysis has ever been devised for angle-based energies, in part because previous analyses rely on deformation invariants that do not generalize to angle-based energies.
In this paper, we address this challenge by devising a new low-dimensional geometric representation that allows us to analyze the eigensystems of this class of energies in closed form. Our analysis applies to a variety of different energies in both 2D and 3D. We first make use of an edge deformation matrix reminiscent of the deformation gradient for thin shells, which then allows us to derive generic, compact expressions for the energy gradient and Hessian. From there, we transform into a sparse, canonical coordinate system that decouples in-plane and out-of-plane movements. Within these reduced dimensions, we are then able to derive an analytical eigendecomposition.
The analytic expressions allow us to perform an in-depth investigation of angle-based energies. We are able to construct a geometric picture for each eigenvector, and directly examine the conditions under which indefiniteness arises. Finally, we show that our derivation can be applied to a variety of angle-based energies, including the bending energy of isotropic strands and thin shells, and results in improved speed and robustness. To summarize our contributions: • A generic formulation for angle-based energies based on an abstract, two-edge stencil.
• An intuitive, geometric interpretation of the eigenvectors.
• An application of our analysis to strand and cloth energies.

Strand Models
Since at least the 1990s, computer graphics papers have used angle-based models for the simulation of strands and hair [Anjyo et al. 1992] by employing multi-body [Hadap 2006], "angular" [Bourguignon and Cani 2000] or "bending" formulations [Daldegan et al. 1993;Ward and Lin 2003]. Some works [Iben et al. 2013;Selle et al. 2008] have approximated these angular forces with linear springs, but these are known to insert non-physical tangential forces that are often undesirable. Implementing angular energies is more involved than linear springs, but by laying out generic eigensystems here, we hope to reduce this complexity for future practitioners.
Cosserat-based strand models also use angle-based bending and twisting [Pai 2002], and have been extended to include anisotropic cross-sections [Bergou et al. 2010]. In this paper, we focus on isotropic cross-sections, which for the aforementioned models reduce to half-angle tangent energies [Gingold et al. 2004]. Quadratic variants have also been useful for modelling constrained strand assemblies [Sueda et al. 2011], which in particular includes knitted cloth [Cirio et al. 2015;Sánchez-Banderas et al. 2020]. Angle-based strand formulations have also been applied in fabrication [Pérez et al. 2015] and architecture [Panetta et al. 2019a]. Finally, concurrent recent work found an analytic eigensystem for unitary bending energies in the course of formulating a new strand model [Wu et al. 2023], but we instead analyze the fully general case here.

Shell Models
Similar to the case of strands, some works approximate angle-based bending energies with linear springs [Bridson et al. 2002;Choi and Ko 2002;Provot et al. 1995], but these again inject tangential forces into the simulation that introduce artificial stiffness in undesired directions.
Angle-based models have also been used for cloth and shell simulation since at least the 1990s [Breen et al. 1994;Volino et al. 1995]. Many different formulations have been proposed, including quadratic [Baraff and Witkin 1998;Grinspun et al. 2003], cos-based [Bridson et al. 2003], sin-based [Wardetzky et al. 2007], and tan-based [Gingold et al. 2004] models. Deriving the derivatives of the bending angle, especially the Hessian required by implicit integrators, can be cumbersome, so Tamstorf and Grinspun [2013] provide explicit expressions as well as architectural recommendations for efficient performance. In their supplement, they provide some intuition regarding sources of indefiniteness, and describe the rank-4 nullspace of the Hessian. Our analysis confirms their nullspace, but goes on uncover closed-form expressions for all the non-trivial eigenpairs.

ANGLE-BASED EIGENANALYSIS
In the following, we will denote matrices with bold uppercase (A), vectors with bold lowercase (a), and scalars with unbolded lowercase ( ). We use vec(·) to denote tensor vectorization [Golub and Van Loan 2013], and adopt the convention that a vectorized matrix is denoted with a lowercase version of the same symbol, e.g. vec(A) = a.

Abstract Edge Formulation
To begin, we assume that the energy is computed based on angles over some discrete object, and that each angle is formed by two abstract edges, e 0 and e 1 (Fig. 2). This can correspond to two adjacent edges in a strand, or two edges along one "flap" of a shell. Once our generic analysis is complete, we will show how to apply our results to these specific domains in §4.
We define an angle-based elastic energy Ψ( ), where is the angle between the two edge vectors: The range of should be [− , ] to fully describe all rotational configurations, but the range of acos is [0, ] and the sign of is ambiguous. To resolve this ambiguity, domain-specific information such as the direction of the bending edge is needed. We abstract this auxiliary geometric information into a sign function S. Under differentiation, S becomes a constant function, aside from a point singularity at 0 that can be addressed using a simple numerical guard. Thus, we can treat S as a constant that does not participate in the ensuing analysis. Without loss of generality, we define a deformation gradient for our edge stencils: Like the deformation gradient, F fully describes the scaling and rotation of the stencil related to , but subtracts off rigid body translation modes. Deformation gradients are usually defined as the linear map component of an affine transformation [Bonet et al. 2021], which maps the rest state to current state. For example, the deformation gradient for membrane stretching can be computed by multiplying the current edge matrix with the inverse of the material space coordinate matrix [Kim 2020]. In the case of a 1D edge stencil, the material coordinates become constant scaling factors that will self-cancel in Eqn. 3, so we can omit them here. Therefore, we can assemble the deformation gradient of our abstract stencil using only the current state edges. We can write the derivatives of Ψ( ) in terms of F where the Jacobian F x takes on different forms for strands and shells. We will later show in §4 that 2 F x 2 either resolves to zero or is of low rank. Thus, analyzing the Hessian 2 Ψ F 2 becomes the most significant task.

2D Eigenanalysis
We start with the 2D case of e 0 , e 1 ∈ R 2 and F ∈ R 2 . The results of this 2D analysis will then serve as a building block for 3D. We assume the form of Ψ( ) can be arbitrary, so we first write the derivatives of the intermediate variable .
In R 2 , the sign of can be assigned with respect to a canonical -axis. In particular, we can define S = sgn( 00 11 − 01 10 ) where is the th component of e . To avoid dealing with higher order notations, we will write the analysis in terms of the vectorized deformation gradient f = vec(F). The gradient and Hessian of in terms of the deformation gradient can be written as the blocks where t = e /∥e ∥ is a normalized version of each edge and C = 0 1 −1 0 is a /2 clockwise rotation ( Fig. 2). The ⊥ subscript is added to a vector whenever it is multiplied by C. e.g. Ct 0 = t 0⊥ . The gradient and Hessian of Ψ( ) can thus be written: where the subscript on H 2 indicates that it is for the 2D case. To reveal the eigenstructure of H 2 ∈ R 4×4 , we observe it has properties relative to the edges and their orthogonal vectors: This leads us to formulate a change-of-basis based on these vectors that then yields a sparsified version of H 2 where The eigensystem of a 4 × 4 matrix can be obtained analytically by solving for the roots of the 4th order characteristic function. Eqn. 14 shows that H 2 can be sparsified, which then yields simpler root expressions. The eigensystem of A −1 2 H 2 A 2 is then given by: where: The eigenvectors of H 2 in the original coordinate system can then be retrieved by re-applying A 2 where q * are written in unnormalized form. This particular change-of-basis does not alter the spectrum of the eigendecomposition, so the eigenvalue expressions remain the same. These expressions may look onerous at first, but most of the complexity comes from the edge length scaling factor , and they reduce to much simpler forms when ∥e 0 ∥ = ∥e 1 ∥: For simplicity, in the above expressions we assume that > 0 and removed some scaling factors. More importantly, we will see in §3.4 that they have clear geometric interpretations.

3D Eigenanalysis
The 3D case is similar to 2D, except that the osculating plane is not static. To track its orientation, we define the binormal b = e 0 × e 1 and its normalized form t = b /∥b∥. The in-plane orthogonal edge vectors are thus defined as t ⊥ = t × t for ∈ {0, 1}, where t again denotes a normalized e . The gradient and Hessian with respect to the vectorized deformation gradient vec(F) = f ∈ R 6 are: where , ℎ, g 0 , g 1 are defined the same as Eqns. 4 and 7. T ∈ R 3×3 where ∈ {0, 1, } is defined as the Levi-Civita tensor E right-multiplied by the vector t such that ∀ ∈ R 3 , × t = T . Relative to Eqns. 5 and 8, the 3D Hessian now includes terms that arise from differentiating b. Similar to the 2D case, we can multiply the edges and their orthogonal vectors onto H to get their projections onto the osculating plane. Interestingly, H displays the exact same characteristics as H 2 in Eqn. 9-12. If we denote the 3D version of A 2 as A ∈ R 6×4 and defined it analogously to Eqn. 13, we obtain HA = AA ⊤ 2 H 2 A 2 . This implies that H is block-diagonal in the 2D coordinate system and that the rank-4 in-plane dynamics is linearly independent from the rank-2 out-of-plane dynamics. Thus, we can write the first four eigenvectors of H in terms of the osculating plane where q ∈ R 6 . The remaining two eigenvectors live in the complementary orthogonal space spanned by the basis of A. If we define q = t ⊤ t ⊤ ⊤ , where is an arbitrary scalar, then q is an eigenvector of H if and only if for some scalar , the following holds: Solving for and , we obtain the last two eigenpairs: Similarly, the complexity here is also dominated by the scaling.
This concludes our eigenanalysis of Ψ( ) for our abstract stencil. Compared to previous positionbased formulations [Tamstorf and Grinspun 2013], our deformation gradient-based formulation and the subsequent transformation into a canonical coordinate system reveals the analytic structure of the Hessian. As our analysis was performed generically on , it can be applied to any isotropic angle-based energy.

Eigenmode Interpretation
The Hessian can be viewed as the Jacobian H = g f , where g = g ⊤ 0 g ⊤ 1 ⊤ is the vector of gradients, and essentially maps the differential df = de ⊤ 0 de ⊤ 1 ⊤ to dg. The eigenvectors then represents directions along which the change in gradient aligns with the motion, i.e. dg = df. We refer to these motions as eigenmodes.
By examining the signs of the entries in the eigenvectors, we can view the geometry of each eigenmodes when > 0 in Figure 3. The top two rows illustrate the in-plane motion, where the first row shows a rotation of the stencil (q 0 ) and a scaling of the edge lengths (q 1 ). The second row shows a compression (q 2 ) and expansion (q 3 ) of the hinge angle. The last row shows out-of-plane rotations about an axis in the osculating plane corresponding to q 4 and q 5 . q 0 q 1 q 2 q 3 q 4 q 5 Fig. 3. Geometry of the eigenvectors. The black vectors are current edges, light green are differentials corresponding to the eigenvectors, and orange the edges e * + de * after deforming along the eigenmode. The first four represents in-plane motion, and the last row describes out-of-plane rotation.
With geometric interpretations of the eigenvectors established, we can now investigate the signs of the eigenvalues. We first examine 0,1,2,3 from Eqns. 20 and 21 which depend on ∈ [0, √ 1 + 4 2 ). By examining the monotonicity of the eigenvalues in terms of , we can see that the signs of the eigenvalues depend on the sign of ℎ. When ℎ > 0, 0 and 2 are non-positive and 1 and 3 are non-negative. When ℎ < 0 the signs switch. Thus, H always contains at least a rank-2 negative-semi-definite subspace.
Taking the indefinite mode q 2 as an example, when the stencil is compressed, there should be increased force to resist the compression. Since the force is scaled with the reciprocal of the edge length, and the edges are stretched, it will be attenuated due to the scaling. Thus, the change in forces will encourage this kind of motion and make the stencil deviate further from the rest state, causing instability.
For the out-of-plane eigenpairs, the signs of 4 and 5 are determined by sin . If sin < 0, then 4 is negative and 5 is positive, while when sin > 0, the signs switch. This means that only one of 4 and 5 can be negative at a time, depending on the value of . Thus, the out-of-plane subspace of H always contains at least one negative eigenvalue, which corresponds to a small torque inducing a sudden out-of-plane rotation.
Using the above analysis, we can assemble an SPD system by ignoring all negative eigenpairs and assembling the Hessian using the outer products of the positive eigenpairs: H SPD = >0 q q ⊤ .

APPLICATION TO RODS AND SHELLS
We now apply our eigenanalysis to different discrete meshes and elastic energies. Our analysis assumes that the energy model being examined depends on an angle, which is computed with respect to two vectors. In order to apply our analysis to a specific energy, we must define a deformation gradient (Eqn. 2) such that an angle is computed using its columns (Eqn. 1). More specifically, we need to assign our abstract variables e 0 , t 1 , Ψ( ) (including its derived quantities and ℎ) and S to corresponding variables in a concrete discrete model, such that e 0 , e 1 and satisfy Eqn. 1. To integrate our analysis into an actual solver, we will then need the derivatives of the deformation gradient with respect to vertex positions.

Elastic Rod Bending Energies
The most straightforward application of our analysis is to the bending energy of elastic rods. This energy involves a stencil with three consecutive vertices (Fig. 4, left) that already closely corresponds to our abstract stencil. In the isotropic case of the Discrete Elastic Rods model [Bergou et al. 2010], the curvature-based energy can be written as an angle-based energy: Other models are instead quadratic in the bending angle [Cirio et al. 2014], where denotes the stiffness and¯is the average rest edge length. For these models, we can compute = S · acos y ⊤ 0 y 1 ∥y 0 ∥ ∥y 1 ∥ . In general, S cannot be determined independent of the twisting, so for simplicity we set S = 1. We assign the abstract edges to e 0 = y 0 , e 1 = y 1 from Fig. 4 for this particular model, which yields the deformation gradient: Then we have the Jacobian of the vectorized F: where I 3 is the 3 × 3 identity matrix and In this case, it becomes clear that 2 F x 2 in Eqn. 3 resolves to zero. Thus, projecting 2 Ψ F 2 to be PSD is sufficient to guarantee that 2 Ψ x 2 is also PSD. With e 0 , e 1 , , ℎ, and S all assigned, we can directly apply the eigenanalysis from §3 to rod bending energies. Specifically, we need to assign these values and vectors to the change-of-basis matrix from Eqn. 13, the coefficients and from Eqn. 22, and the out-of-plane rotation eigenpairs from Eqns. 30-32. After obtaining the F-based eigensystem, plugging the Jacobian of the deformation gradient into Eqn. 3 will generate forces and force gradients for any FEM-based solver.

Shell Bending Energies
Many shell models define their bending energy using a hinge angle, and our abstract stencil can compactly characterize many of the different Ψ cloth ( ) proposed in previous papers. Here we list the derivatives of some common bending energy models with respect to .
The quadratic Discrete Shells model [Grinspun et al. 2003] becomes: The cos-based model from [Bridson et al. 2003] resolves to: The sin-based energy from [Wardetzky et al. 2007] becomes: where the hinge angle is defined as the bending angle formed by two neighboring triangles with a shared "hinge" edge. The right side of Fig. 4 shows a hinge stencil on a triangle mesh surface, where the three edges characterizing the stencil are In other works, the hinge angle is computed using the normal vectors of the two triangles (see e.g. [Tamstorf and Grinspun 2013]). However, to more closely mirror our abstract stencil, we compute the hinge angle as = S − acos z 0 ·z 1 ∥z 0 ∥ ∥z 1 ∥ , where are versions of y 2 and y 3 projected to the orthogonal plane of y 1 . This is equivalent to the normalbased definition, but also captures the stencil scaling, which plays a significant role under differentiation. Thus, we assign the abstract edge vectors as e 0 = z 0 and e 1 = z 1 , and the deformation gradient then becomes The unit binormal can then be defined as t = z 0 ×z 1 ∥z 0 ×z 1 ∥ , which is parallel to the hinge edge y 1 . We can then assign S to a comparison of the binormal vector to the hinge vector: S = t · t 1 .
At this point, we have already assigned e 0 , e 1 , S, and ℎ so that we have everything we need to assemble the eigensystem of the Hessian. We can substitute t 0 and t 1 from Eqn. 4 with 0 = z 0/∥ z 0 ∥ and 1 = z 1/∥ z 1 ∥ to get the F-based gradient. Now we have everything we need to assemble a positive definite F-based Hessian based on the analysis in §3. To obtain the position-based derivatives in Eqn. 3, we will need the Jacobian of the deformation gradient: Because our novel definition of deformation gradient is based on projected vectors, we need to be particularly careful with the chain rule in Eqn. 3. Unlike many FEM-based elastic energy models, the first term of Eqn. 3 is nonzero. Instead, it becomes: where W ∈ R 12 . To guarantee the positive-definiteness of the position-based Hessian, we need to efficiently filter W, whose full analytic expression we give in Appendix A. Since W is a rank-4 matrix composed of outer products, its null space is easy to analyze. We can then find a basis in the orthogonal complement of the null space that can further sparsify W We will present the definitions of p and show that they encode motions that complement the 6-dimensional F-based eigensystem in Appendix A. We can write W in terms of the new basis: Now the eigenvalues of A −1 W WA W can be computed by solving a depressed quartic equation analytically, which has existing, mature, and efficient implementations [Khashin 2014]. For an eigenvalue , the corresponding unnormalized eigenvector can be computed by: Our experiments have shown that two out of the four eigenvalues of W are always non-negative. Therefore, to perform semi-positive definite projection on W, we find the two non-negative eigenvalues and assemble their corresponding outer products.

DISCUSSION AND RESULTS
Our analysis is applicable to many angle-based energies and can be integrated into an FEM solver. We have included in our supplemental materials Matlab/Octave implementations that verify the correctness of our derivation and expressions using finite difference tests. We show a variety of elastic rod and thin shell simulations to demonstrate the robustness of our method across different configurations, and under large time steps.

Implementation.
To test the performance of our model on strands, we simulated a stiff net of crossing wires, such as those used in architectural [Panetta et al. 2019b] applications. We inserted a bending term between each pair of intersecting rod segments, resulting in six bending stencils for each internal vertex in a checkerboard net (Fig. 5). We applied the stretching energy from Discrete Viscous Threads (DVT) [Bergou et al. 2010], but omitted the twisting energy because torsion is minimal in such constrained configurations. For bending, we tried both the quadratic model from Eqn. 34 and the tangent-based model from Eqn. 33. We use a backward Euler [Baraff and Witkin 1998] integrator with a large time step of = 1 /30 s, and solved the global system using a Cholesky solver. We converted to using the formula = 4 4 where = 0.1 cm.

Results.
We test the robustness of our model under three different stiffnesses: = 1 10 Pa, = 1 11 Pa, and = 1 12 Pa. We compared our filtered approach to an unfiltered method and a Gauss-Newton approximation [Choi and Ko 2002] that only includes the leading outer product term of the Hessian (Fig. 6). . Cloth with bending stiffnesses of = 1, 10, 1 2 , 1 3 , 1 4 and 1 5 . All are simulated using a quadratic hinge energy [Grinspun et al. 2003]. Unfiltered and Gauss-Newton approximations explode under a variety of scenarios, which are denoted with an ×. Our filtered Hessians produce smooth results in all cases.
As can be seen in the video, at = 1 10 Pa, the unfiltered, Gauss-Newton and filtered (ours) approaches are all relatively stable. The Gauss-Newton approach contains some jittering at the beginning (Fig. 7a), but stabilizes in later frames.
Under larger stiffness of = 1 11 Pa and = 1 12 Pa, bending forces start to dominate and both unfiltered and Gauss-Newton explode immediately after the simulation starts. Fig. 7 shows the frames right before explosion. For both methods, the instability emerges at the corners and boundaries, and propagates to the whole model very quickly. We also found that Gauss-Newton is generally less stable than the unfiltered method and explodes faster. In contrast, our method is totally stable and settles to an equilibrium for both the quadratic and tangent-based energies. Our method efficiently simulates stiff strand bending, even with large time steps.

Simulating Shells
5.2.1 Implementation. For shells, we constructed a row of V-shaped cloth strips with varying stiffness, and set the rest angle of the V to /4. One end of the cloth is attached to a wall, and collisions with the wall are handled with Baraff-Witkin filtering [Baraff and Witkin 1998]. We use quadratic bending (Eqn. 37), and ARAP-type membrane stretching [Kim 2020]. The scene is simulated with backward Euler with a time step of = 1 /30. For the linear system solver, we used both preconditioned conjugate gradients (PCG) and Cholesky method.

Results.
We again compare our method against unfiltered and Gauss-Newton strategies. Fig. 8 shows frames from all three methods using PCG, where a red cross appears when the simulation exploded. In the first row of unfiltered Hessians, we see that the simulation is unstable under high bending stiffnesses. As the strips droop under gravity, the bending energy increases, and local instabilities quickly cause the simulation to explode.
The second row shows the results using Gauss-Newton. Unlike the unfiltered case, Gauss-Newton is able to keep the system semi-positive definite under large stiffnesses, as can be seen in the rightmost = 1 5 strips. However, the middle strips show that Gauss-Newton is unstable . Comparison of the three approaches with bending stiffnesses of = 1, 10, 1 2 , 1 3 , 1 4 and 1 5 using the same setting as Fig. 8 except that the linear system is solved using the Cholesky method.
= 1 10 1 2 1 3 1 4 1 5 Avg.  (Fig. 8) timings. The three methods being compared are the unfiltered baseline, our method, and a brute-force numerical eigendecomposition. All timings are in milliseconds. The last columns shows average timing across all stiffnesses.
under mid-range stiffnesses that undergo larger deformations. Some strips explodes when they encounter a kinematic collision with the wall. The second strip from the right does not explode, but drifts to one side and non-physically breaks the symmetry of the system.
Our method is shown in the last row of Fig. 8, and is consistently stabler than both unfiltered and Gauss-Newton. Our filtered Hessian is able to fully describe the significant eigenmodes of the system while discarding unstable modes. Fig. 9 shows the simulation results using a Cholesky solver. Cholesky is usually more sensitive to instabilities when handling large linear systems, so both the Gauss-Newton and unfiltered methods explode much earlier than PCG. The second row shows that with Gauss-Newton, the rightmost strip with largest stiffness drifts to one side and hits the wall, eventually exploding. This is because Gauss-Newton is not able to adjust to sudden local displacements caused by kinematic collisions. The last row shows that our method also generates asymmetric motion, which we hypothesize is due to numerical noise pushing the simulation out of an energetic saddle point. Regardless, our method remains stable.

Crushed Cylinders
To test the robustness of our method on more complicated examples, we crushed stiff cylinders composed of both thin shells and wire nets (Fig. 1). We attach the top and bottom of the cylinder to two cubes that then crush the geometry. We simulated the thin shell using a bending stiffness of = 1 12 1 11 1 10 Avg.
Baseline 1764  1775 1752 1764  Ours  1898  1866 1879 1881  Numerical  12689  12920 13203 12937   Table 2. Net scene (Fig. 6) timings. All timings are in milliseconds and the Young's moduli are in Pascals. Methods are the same as in Table 1. = 1 6 , and the wire cylinder under the Young's modulus of = 1 12 Pa. Other settings are the same as §5.1 and §5.2. Non-trivial buckling patterns form and evolve during the simulation, and the dynamics remains stable even under sudden, large deformations. In the video, we show that compared to our method, both unfiltered and Gauss-Newton failed to produce the wrinkles.

Timing and Efficiency
We measure efficiency of our method by using unfiltered Hessian construction time as a baseline, and comparing it to our method, as well as a brute-force numerical call to SelfAdjointEigenSolver from the Eigen library [Guennebaud et al. 2010]. In Tables 1 and 2 we list the total time spent, in milliseconds, of computing the bending energy Hessian over 400 simulation frames. All the tests were implemented in C++ using the Eigen library, and were run on an Apple M1 Pro chip with 10-core CPU and 32 GB memory.
Our method is 6× to 7× faster than the brute-force numerical call. Compared to the unfiltered baseline, our method adds a negligible 14% to the shell scene and 7% to the net scene. Considering the robustness under large time step that our method adds, this trade-off is acceptable.

CONCLUSIONS AND FUTURE WORK
We present the first clean eigenanalysis of angle-based deformation energies based on an abstract stencil. Our analysis is physically interpretable and reveals the intrinsic geometric structure behind the energy Hessian. We provide a semi-quantitative interpretation of the eigenmodes of angle-based energies, and a fast scheme for assembling a semi-positive definite Hessian per stencil. We apply our method to elastic rod and shell simulation and have achieved consistent improvements to stability compared to popular existing methods, while also maintaining efficiency.
Per-element SPD projection can be suboptimal, because the finite element approximation can be sensitive to the discretization scheme. In some extreme cases involving large velocities, fixing the indefiniteness locally cannot prevent the instability from propagating to the global modes. Whether there is a projection scheme that can account for global low-frequency dynamics is future work. Even for a local per-element scenario, a fully analytical eigensystem is not guaranteed to exist for energy types not covered in this paper, such as anisotropic rod bending.
A general eigendecomposition that is applicable to this type of angle-based energy can have other potential applications, such as exerting controls and constraints over the simulation system in a physically valid way. Thus, generalizing our ideas to broader applications remains future work.

ACKNOWLEDGMENTS
We thank the reviewers for their comments. This work was supported by Adobe, the Bungie Foundation, Teng and Han Family Fund, and National Science Foundation (IIS-2132280). Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

A FULL EXPRESSIONS OF THE SHELL EIGENANALYSIS
We present the full expression of the first term in Eqn. 3. This is used to derive the eigensystem of W in §4.2. We also implemented this for experiments. First, we give the building blocks based on the block-wise formulation in Eqn. 43: Then the full expression of W ∈ R 12×12 is: We can see that W is sparse and each block is an outer product. This property enables fast matrix assembly, and leads us to analyze the structure of the null space and row space. Next, by investigating the structure of W, we found four vectors in the row space that has an interesting property: Wp 0 = − 0 p 0 + 1 p 2 Wp 1 = 0 p 1 + 1 p 3 where p 0 and p 1 are defined as and computing p 2 and p 3 using Eqn. 55 yields: where p 2 and p 3 are displacements along the hinge direction t 1 for each of the four vertices. The projected edges z 0 and z 1 is invariant in this case and these modes p 2 and p 3 are not captured in the F-based Hessian. The p 0 and p 1 directions represent a rotation of the hinge edge y 1 , and they are supplementary to the rotation modes of the F-based Hessian q 4 and q 5 . After further exploiting Eqn. 54 and Eqn. 55 using the symmetry of W, we are able to construct the matrix in Eqn. 47 in the chosen basis. Because assembling the 12 by 12 filtered W directly from two outer products of 12-dimensional vectors is relatively lightweight compared to the F-based term, the shell-based variation of our analysis remains efficient as shown in §5.4.