Second-Order Finite Elements for Deformable Surfaces

We present a computational framework for simulating deformable surfaces from planar rest shape with second-order triangular finite elements. Our method develops numerical schemes for discretizing stretching, shearing, and bending energies of deformable surfaces in a second-order finite-element setting. In particular, we introduce a novel discretization scheme for approximating mean curvatures on a curved triangle mesh. Our framework also integrates a virtual-node finite-element scheme that supports two-way coupling between cut-cell rods without expensive remeshing. We compare our approach with traditional simulation methods using linear and higher-order finite elements and demonstrate its advantages in several challenging settings, such as low-resolution meshes, anisotropic triangulation, and stiff materials. Finally, we showcase several applications of our framework, including cloth simulation, mixed Origami and Kirigami, and biologically-inspired soft wing simulation.


INTRODUCTION
Linear triangles have served as a fundamental computing primitive, accommodating various cloth and shell simulations, ever since its first success in simulating large-timestep cloth a quarter of a century ago [Baraff and Witkin 1998].By discretizing a deformed surface with planar triangles that share edges and vertices and solving elasticity using the finite element method (FEM), these triangle meshes can produce a broad range of visually appealing surface phenomena, such as cloth wrinkles and paper ceases.However, linear triangles contain several inherent weaknesses.First, in terms of its geometric expressiveness, the planar assumption of each triangle requires a large number of elements to approximate a complex, curved surface.Bending deformation can only be expressed with the dihedral angle between triangles that share an edge.Remeshing operations, sometimes requiring highly anisotropic elements, must be performed frequently to characterize bending behaviors that are not aligned with mesh edges (e.g., see Narain et al. [2012]).Second, in terms of dynamics, linear triangles assume a piecewise constant deformation gradient within each element, resulting in a demanding prerequisite on the mesh quality (e.g., see Shewchuk [2002]).The simulation accuracy will be negatively affected if a mesh contains low-quality triangles.Third, the linear geometry and constant deformation gradient further cause kinematic locking artifacts on a triangle mesh, especially when its resolution is low, which produces unnatural bending and overly damped cloth dynamics.In short, simulations on linear triangles are highly mesh-dependent.High-quality meshes (in terms of both mesh resolution and element aspect ratio) are essential for accurate simulations.
In recent years, high-order triangles have sparked interest in computer graphics research.These methods stem from the highorder finite element analysis in mechanical engineering [Johnson 2012], by placing multiple degrees of freedom on edges or interiors of finite elements.In this way, a triangle is no longer planar but rather has its shape described by a high-order polynomial and approximated with high-order basis functions.As pointed out by Schneider et al. [2018], the simulation accuracy can be decoupled from the mesh quality when employing high-order FEM schemes.Following this thread, a natural next step is to explore the efficacy of high-order methods in simulating deformable surfaces, with a particular focus on addressing the various mesh-dependent bending, locking, and coupling artifacts in cloth or thin-shell simulation.Nonetheless, devising a high-order cloth and shell simulator to tackle these highly dynamic and nonlinear physical problems is difficult.One of the main challenges is to appropriately define the stretching, shearing, and bending energy and their discretization that fits a high-order discretization.To begin with, a curved triangle can bend due to its intrinsic non-planar shape.In addition, two curved triangles connect on a curved edge, which voids the possibility of characterizing bending with dihedral angles.To the best of our knowledge, there exists no practical algorithm that can support cloth and shell simulations with high-order finite elements for graphics applications.
Filling such a gap, our paper proposes a second-order deformable surface simulation algorithm on curved triangle elements.At the heart of our approach is a novel elastic energy definition in conjunction with its discretion on non-planar triangles to accommodate the second-order simulation of elastic surfaces.Furthermore, we augment our simulator with a novel scheme to support secondorder cut-cell finite elements, enabling native two-way coupling with elastic rods in our simulation.We evaluate our simulation on different test scenarios to demonstrate its merits in producing dynamic, smooth, and deformable surface simulations.
Our main contributions include the following: • We present a discretization scheme for mean curvature and its induced bending energy on curved triangulation of smooth surfaces; • We present a novel virtual-node finite element scheme to couple deformable surfaces with cut-cell curves, enabling creases and coupling with elastic rods at flexible locations; • We validate the efficacy of the proposed framework on lowresolution, low-quality triangulation and compare it with an existing open-source finite-element implementations, • We demonstrate the application of our method in cloth simulation, mixed Origami and Kirigami, and biologically-inspired soft-wing simulation.

RELATED WORK
Deformable-surface simulation.Deformable-surface simulation in computer graphics originated from the early endeavors in cloth animation [Baraff and Witkin 1998;Breen et al. 1994;Eberhardt et al. 1996;Terzopoulos et al. 1987], which model cloth as an elastic body.To counteract the numerical instability evoked by the high stiffness necessary for stretch resistance, Baraff and Witkin [1998] introduced an implicit temporal-integration scheme to enable large time steps with high stability, at the cost of significant numerical damping.To reduce such artifacts, several implicit-explicit solving schemes [Boxerman 2003;Eberhardt et al. 2000;Stern and Grinspun 2009] are developed by treating stiff and non-stiff forces differently.Recently, a line of projective dynamics methods [Bouaziz et al. 2014;Ly et al. 2020] has also been devised to solve the implicit integration problem from the variational perspective efficiently.To discretely model the deformation and bending mechanics of thin shells, a number of convenient hinge-based models have been developed [Bridson et al. 2003;Gingold et al. 2004;Grinspun et al. 2003]; we adopt the paradigm of Grinspun et al. [2003] that computes the bending energy by comparing squared mean curvatures.To date, a large gamut of thin shell phenomena have been successfully captured, including plastic folding and crumpling [Narain et al. 2013], large-scale tearing and fracture [Busaryev et al. 2013;Pfaff et al. 2014], interaction with ambient fluid [Chen et al. 2013;Sifakis et al. 2008], reaction to environmental stimuli [Chen et al. 2018], etc.In addition, the accurate and affordable simulation of thin shells has facilitated a range of inverse design tasks, such as the natural shape of shells [Ly et al. 2018], FlexMaps [Malomo et al. 2018], deployable X-shells [Panetta et al. 2019], and inflatable surfaces [Panetta et al. 2021].
Higher-order methods.High-order discretization schemes and algorithms offer a competitive alternative to their low-order counterparts in resolving physical problems at higher fidelity at a comparable computational cost [Schneider et al. 2022;Vincent and Jameson 2011].For example, higher-order continuous Galerkin methods [Vos et al. 2010;Zhou and Lu 2005] and discontinuous Galerkin methods [Huynh et al. 2014;Jameson 2011;Wang et al. 2013] are developed to improve accuracy and efficiency in computational structural or fluid mechanics.Interested readers can refer to Johnson [2012] for a comprehensive discussion.Similarly, higher-order schemes have been bridged with the finite volume [Barth and Frederickson 1990;Harten et al. 1997] and finite difference [Shu 2003;Vasilyev 2000] paradigms.As identified by Vincent and Jameson [2011], a hindering challenge for adopting these high-order schemes is the generation of unstructured, high-order (curved) meshes.A number of works take on this challenge by directly generating meshes from CAD models [Dey et al. 1999] or by post-processing linear meshes [Ims et al. 2015].In hybrid Lagrangian-Eulerian simulation, methods that transfer higher-order velocity modes have been shown to provide better energy and momentum conservation for both fluid [Jiang et al. 2015] and solid materials [Hu et al. 2018].Higherorder level-set methods have also been developed to enhance mass conservation in multi-phase simulation [Moghadam et al. 2016], promote smoothness in surface construction [Bajaj et al. 2007], and preserve sub-grid structures [Nave et al. 2010].
Higher-order shells.Exploration of higher-order methods for shell simulation keeps being an important and major topic especially for structural analysis in mechanical and civil engineering.Previous works primarily target at solving (quasi-)static problems with small deformation.Many previous works [Batoz et al. 1980;Belytschko et al. 1984;Macneal 1982;Phaal and Calladine 1992;Sze and Zhu 1999;Tessler and Hughes 1985;Zhongnian 1992] proposed shell elements to achieve an accurate and consistent solution with the assumption of geometric linearity.To accurately present the curved mid-surface of a shell, lots of conventional shell elements introduced rotational degrees of freedom [Bathe and Ho 1981;Batoz et al. 1980;Belytschko et al. 1984;Campello et al. 2003;Macneal 1982;Tessler and Hughes 1985;Zhongnian 1992].This category also includes lots of existing shell algorithms in commercial software, e.g., S4 and S4R in ABAQUS [Abaqus 2011].In computer graphics, subdivision-surface shells are of particular interest to researchers [Cirak et al. 2000].A series of works have been performed to solve static and dynamic problems using different subdivision strategies, e.g.Catmull-Clark subdivision [Wawrzinek et al. 2011].
Previous higher-order shell elements have explored widely the methods of improving accuracy and consistency for simulating loads and displacements for (quasi-)static shells with nonzero thickness.However, scenes in graphic application usually require methods supporting stable simulation for dynamic large-scale deformation and convenient coupling with other structures, which we believe cannot be easily realized with the existing approaches.
The assumption of material or geometric linearity often leads to a visually implausible shape when simulating large deformations.Also, rotational DoFs suffer from potentially slow convergence and make crease handling considerably more expensive [Hu et al. 2021].Finally, handling flexible coupling with creases or rods in some of these higher-order FEMs is typically challenging.For example, Wawrzinek et al. [2011] and Cirak et al. [2000] assumed C1 continuity on thin shells, so handling creases is not straightforward.Laurent and Rio [2001] used only 9 DoFs to model a curved triangle, which would lead to incompatible boundary conditions if coupled with cut-cell creases or rods.

CONTINUOUS MODEL
Kinematics.We consider a deformable surface in R 3 with zero thickness and a planar rest shape Ω ⊂ R 2 .We use X = (, ) ∈ Ω to denote a material point in the rest shape.A spatial-temporal function x : Ω × R + → R 3 defines the motion of each material point X(, ) at time  as x(X(, ), ) or x(X, ) for brevity.We use x() to denote the surface at time , i.e., x() = x(Ω, ).The velocity and acceleration of each material point are denoted as v = x and a = x.The Jacobian of a surface point w.r.t.its (, ) coordinates is denoted as F with its subscript indicating the surface, e.g., F x( ) : Ω → R 3×2 returns the Jacobian of the material point X(, ) at time .
where  is the material density.We use f int and f ext to denote the internal and external force densities exerted on the rest shape.We assume that f int consists of the stretching and shearing force f ss and the bending force f bd , both of which are conservative forces derived from the corresponding potential energy  ss and  bd .We complement Eqn. ( 1) with time-varying Dirichlet boundary conditions x(X, ) = x  (X, ) on Γ  ⊆ Ω where x  prescribes the kinematic motion on X(Γ  ).
Stretching and shearing energy.We follow previous works [Kim and Eberle 2020] to assume a hyperelastic material model.Therefore, the strain energy for stretching and shearing can written as  ss = ∫ Ω Ψ ss (, ), where the energy density Ψ ss is defined as a function of the Green strain tensor  = 0.5(F ⊤ x F x − I).In our implementation, we applied St. Venant-Kirchhoff model [Sifakis and Barbic 2012].
Bending energy.We consider the bending energy as an integral  bd = ∫ Ω Ψ bd (, ) of the squared mean curvature [Grinspun et al. 2003]: where  x( ) : Ω → R stands for the mean curvature and  the bending stiffness.

NUMERICAL DISCRETIZATION 4.1 Variational Form
We discretize Enq.(1) with the backward Euler time integration, leading to the following energy minimization problem [Kane et al. 2000]: (3) subject to the Dirichlet boundary conditions.Here, the superscript  indexes the time step with ℎ its step size.

Higher-Order Finite Elements
We use the finite element method to discretize Eqn.(3).We assume Ω is a polygonal domain and partition it into triangles {  }.These triangles together contain  basis functions Depending on the order of finite elements, these nodes may include triangle vertices, edge midpoints, or triangle centers.The basis functions can be linear, quadratic, or higher-order polynomials on each triangle (Fig. 2) and can be precomputed.For any function  on Ω, we define f , an approximation to  , as the sum of all basis functions   weighted by predefined nodal values   , i.e., and it is common to use   =  (  ,   ) whenever possible.f is a continuous function on Ω, allowing us to approximate  at any (, ) ∈ Ω with finite degrees of freedom  1 ,  2 , • • • ,   (Fig. 2).Note that a triangle becomes curved in a second-order finite element.
The integral of a function  can be computed by quadrature rules: where {q  ∈ Ω} are the quadrature points and   their weights.We discretize X and x at triangle nodes to obtain X  , x  ∈ R 3× , on which we use the second-order finite element scheme to discretize the energies in Eqn.
(3).The main challenge is to discretize the bending energy due to curved triangles, which we detail below.

Energy Discretization
Previous works [Baraff and Witkin 1998;Gingold et al. 2004;Grinspun et al. 2003] have studied discrete bending energy on triangular surfaces extensively and formulated it as a function of dihedral angles.Generalizing this idea to second-order elements is non-trivial for two reasons: the discretized surface now consists of curved triangles with nonzero curvatures, and these triangles now meet at quadratic edges with varying dihedral angles.To resolve them in a principled way, we recall the discrete mean curvature measure from Cohen-Steiner and Morvan [2003] and generalize it to curved triangular surfaces: Definition 4.1 (Discrete mean curvature measure).Consider a discretized surface Π built from a second-order finite element partition {  } of its (, )-parametrization in Ω.For every (Borel) set  ⊂ R 3 , we define the discrete mean curvature measure  Π () as where  Π (p) is the mean curvature of Π at a surface point p,    is the edge shared by the triangles   and   , and ( ∩ {Π(   )}) computes the signed area swept by the arc spanned by the outward unit normals of Π(  ) and Π(  ) along the edge  ∩ {Π(   )}.The sign of  is positive if    is convex and negative if concave.
Fig. 3 shows a geometric interpretation of the second term in  Π (), referred to as the "tube formula" in Cohen-Steiner and Morvan [2003].The definition states that the contribution of an edge to the mean curvature equals the (incomplete) signed surface area of the tube.Our  Π () is a natural extension of the discrete mean curvature measure in first-order elements, in which case the first term in  Π () vanishes, and the second term becomes the dihedral angle times the edge length (Fig. 3

left).
We now return to discretizing the bending energy with the discrete mean curvature measure.We aim to provide a discretization scheme that can be computed with quadrature on elements.To achieve this, it is sufficient to explain how to obtain an estimate of  x+1 at each quadrature point q  , where the hat symbol indicates the surface is discretized with FEM in Eqn.(4).We first subdivide each   into three triangles    using the center of   , where  denotes the triangle edge index.Next, for each edge    , we consider a ) (Fig. 3).Eq.6 suggests the following mean curvature estimate for any (, ) ∈ where the denominator computes the surface area of the two deformed subdivided triangles.The intuition is that {   ∩ x+1 } constitutes a partition of the surface x+1 .Therefore, if we uniformly distribute the edge contribution  onto every point in the two adjacent triangles, integrating H on the whole surface x+1 equals  x+1 (R 3 ) as desired.We detail our implementation of the bending energy in the supplementary material.

COUPLING
We now extend our numerical model in Sec. 4 to couple surfaces with embedded curves with potential applications in simulating creases, wrinkles, or elastic rods.We skip the analysis on the continuous model and present our results in the discretized form.For simplicity, we will illustrate the basic idea using one curve on the planar rest surface.
Geometric representation.We recall the domain Ω with a triangulation {  } and consider a curve embedded in Ω with both ends on its boundary.We discretize the curve by a series of line segments {  } connecting the intersections of the curve with triangle edges, where the subscript  is the index of the triangle that contains   (Fig. 4).Each   now cuts   into two pieces    and    that can deform in their own ways as long as they remain coincident on   , which can also stretch or bend along with the surface deformation.
Cut-cell creases.An immediate solution here is to sub-divide every quadrilateral    and    into two triangles and apply Sec. 4 to the new triangulation.However, this conceptually simple method suffers from the complexity of adaptive remeshing.Instead, we ground our method on the virtual node algorithm [Bedrossian et al. 2010] to avoid the expensive remeshing of Ω.Consider a specific   that splits   into    and    .Each of them obtains copies of nodes in   which control its deformation (Fig. 4

middle and right).
To ensure    and    remain coincident on   after deformation, it is sufficient to check they agree on the positions of a finite number of nodes in   .For example, when   is a first-order finite element,   remains a line regardless of the deformation, whose location can be fully determined by checking its two endpoints.Noting that the (, )-coordinates of these endpoints are known beforehand, we conclude from Eqn. (4) that their deformed positions are linear on the degrees of freedom at virtual nodes.Therefore, such constraints are linear, which forms a large linear equality constraint Ax = b in Eqn. 3. Here, the notation has flattened x into a column vector, and A is sparse and typically has much fewer rows than columns.From counting the rows and columns of A, we notice that this combination of virtual node schemes with second-order instead of first-order finite elements is not arbitrary.Interestingly, coupling it with first-order elements leads to a loss of degrees of freedom on the curve, which we detail in the supplementary material.
In summary, the full numerical method in Sec. 4 remains valid but requires three modifications: First, the quadratures need to be sampled on    and    instead of on the full   .Second, we need to add extra bending energy along   .Finally, the numerical optimization now includes slightly more linear equality constraints, which requires us to solve a symmetric indefinite KKT system at each Newton iteration [Boyd et al. 2004].
Dynamic curves.Handling dynamic curves, e.g., elastic rods with their own density and elasticity, requires small algorithmic updates on the method stated before.First, we add conventional spring energy between neighboring vertices on a rod to our previously defined potential energy.Second, we add to  bd the bending energy from the curve introduced by Bergou et al. [2008], which is defined by integrating the squared curve curvature along the curve.

Implementation Details
We implement our method in C++ on a personal workstation with moderate hardware specifications (12 CPU cores and 16 GB memory).We implemented Newton's method with the standard line search and Hessian definiteness fix techniques [Nocedal and Wright 2006] and used SuiteSparse [Davis 2023] to solve sparse linear systems.We used the geometry library from Sharp et al. [2019] for geometry processing like triangulation and subdivision.Contact and collision were handled with our implementation of IPC [Li et al. 2020].Simulating one timestep in our largest experiment (80K second-order elements with 480K degrees of freedom) took 30 seconds to finish.The total time for running each experiment ranges from several minutes to at most two hours.

Evaluation
We compare our approach with first-order elements and analyze their differences in handling low-quality triangulation, kinematic locking, and computational cost.Additionally, we compare our implementation with other higher-order methods implemented in OOFEM, an open-source finite element implementation.
Low-quality meshing.For the finite element method, the choice of basis functions and the mesh quality are two critical factors that jointly affect its numerical performance [Johnson 2012].Because quadratic basis functions are more expressive than linear ones, we speculate that our approach is more tolerant of low-quality meshes than traditional first-order methods.An immediate implication is that we anticipate second-order elements to generate similarquality simulation results as first-order elements but with a coarser triangulation.To validate this hypothesis, we consider a motivating example of a piece of cloth falling on a static, rigid sphere.We deliberately used an extremely coarse triangulation and visualized the simulation results from first-order elements (32 triangles) and second-order elements (8 triangles) in Fig. 5.The planar triangles from first-order elements essentially quantized the underlying smooth motion, and noticeable artifacts from the piecewise-linear shapes are evident under such a low resolution.On the other hand, the curved triangles captured the motion more convincingly and exhibited a smoother surface, thanks to the expressiveness of their quadratic basis.This experiment indicates that our method can be a promising approach when high-resolution meshes are unaffordable.We include more experiments about the performance of our finite element methods on low-quality triangulations in the supplementary material.
Kinematic locking.Strong stretching coefficients rigidify firstorder triangular elements due to their planar nature and results in the kinematic locking issue, i.e., a multi-facet surface that can only bend along prescribed edges.On the contrary, our model uses bendable curved elements that alleviate the strong coupling between stretching and bending.We provide results of two experiments to illustrate this property.In each experiment, the mesh resolutions were chosen so that linear and quadratic elements created a similar number of nonzeros in the energy Hessian, resulting in a similar computational time in total (Table.1).
In the first experiment, we simulated a scene Saddle with four groups of test cases (Fig. 6).Each group compares the quadratic elements on a square with linear ones.The resolutions of quadratic and linear cases are 10 × 10 × 2 and 30 × 30 × 2 triangles.We vary the stretching stiffness (Young's modulus) across groups.We fix the locations of a pair of diagonal corners and simulate the motion of the square under gravity with an expectation of a saddle surface.The results are shown in Fig. 6.We adjust the bending stiffness such that when the diagonal edges are aligned with the bending direction, the quadratic and linear cases produce similar results.However, when the diagonal edges are orthogonal to the bending direction, the linear elements clearly show resistance to bending.This result shows the bending ability of linear elements is dependent on the triangulation of the surface.In contrast, the quadratic elements are better able to bend despite the same direction of triangulation.This observation suggests that our second-order approach remedies the kinematic locking issue.
In the second experiment, we simulated a scene Curtain with two groups of comparisons at different resolutions.Each group simulated a swinging curtain using linear and quadratic elements.Fig. 7 shows that under a similar computational time budget, our method produced a smoother surface at each frame that contrasts the polygonal artifacts in the linear elements caused by kinematic locking.The result suggests that our method is practical for alleviating kinematic locking without performance degradation.Table 1: Computational cost of the linear ("L") and quadratic ("Q") elements in Fig. 6 Saddle and Fig. 7 Curtain."NNZ" stands for the number of nonzeros in the resultant energy Hessian.The number after "Q-" and "L-" indicates the mesh resolution.
Comparison with other higher-order elements.As mentioned in the related works, lots of previous higher-order shell elements assumed geometric linearity, which often leads to a visually implausible shape with large deformation.To illustrate this, we compare our method with triangular shell elements in OOFEM [Patzák 2012], an open-source finite element implementation.We compare with three linear elements tr_shell02, tr_shell01 and RerShell which are not suitable for large deformation due to their assumption on geometric linearity and their lack of coupling between stretching and bending.Given a plate square mesh (80 × 80 × 2 triangles) on the XY-plane, we fix one boundary to have it swing under gravity.All three types of elements showed artifacts with large deformation (Fig. 8 and supplementary material).
For nonlinear shell elements, rotational DoFs often leads to convergence issues [Hu et al. 2021].We experimented Tr2Shell7, the only triangular element in OOFEM that supports large deformation, which showed slow convergence in a cantilever beam test with moderate loads.On a 5 × 5 × 2 (50 elements, 6 nodes per element, 7 DoFs per node) shell, with 0.03/0.04/0.044dimensionless loads per node on free edge, it requires 25/84/586 iterations to converge with Newton's method and failed to converge with larger loads.In comparison, our method converged within 5 iterations in a similar problem.

Applications
We present seven examples in cloth simulation, mixed Origami and Kirigami structures, and their coupling with elastic rods.Large curtain: This example simulates a curtain on 200 × 200 × 2 triangular elements (Fig. 9).The curtain swings under gravity while five evenly-spaced rings on top of the curtain move along a rod.The simulation shows an elastic, highly vibrant curtain with detailed wrinkles.
Lanterns: This example simulates two lanterns (Fig. 1) with a cylindrical mesh on a resolution of 200 × 200 × 2 triangles.There are 20 rods traversing the cylinder and no edge bending energy on edges where rods reside.The bottom of the cylinder is set to be fixed, and the top moves downward to squeeze the whole structure.With different stretching stiffness, the cylinders deform to different shapes : a smooth ellipsoid and a shape with curved creases.
Bat: This example simulates a bat wing on a mesh with 3600 elements.Elastic rods are inserted in the mesh by virtual nodes to simulate arms and fingers (Fig. 12).External periodic forces are applied at joints to simulate flying movements of the bat.
Tori: This example simulates a piece of cloth with 459 elements falling on two tori from above (Fig. 10).The cloth eventually rests on both tori and shows a smooth surface capturing the underlying tori shapes with the small number of curved triangles.
Origami and Kirigami: We present three more examples that simulate the dynamics of mixed Origami and Kirigami structures.Fig. 14 shows a paper-art work whose structure consists of an array of slanted paper tapes with Origami creases at linearly interpolated locations.The Kirigami design is simulated with only 572 triangles, yet it manifests smooth, large-scale elastic behaviors with visually appealing deformation and vibration.Fig. 13 shows a second Kirigami example on a 32 × 32 square mesh under gravity.Finally, we simulated another Kirigami design of a paper discretized into 600 triangular elements (Fig. 11).Folding the paper along its major crease produces a stair structure resulting from the spatially varying crease location at each paper tape.One of the four boundaries of the paper is fixed to a static rod so that when spinning the paper dynamically, it automatically folds to generate the final Kirigami.

CONCLUSIONS
This paper presented a novel second-order deformable-surface simulation algorithm.The framework discretizes a deformable surface using finite elements on curved triangles, which overcomes the prolonged limitations of linear element methods in tackling problems of anisotropic bending, locking, and low-resolution dynamics.We compared our approach with other state-of-the-art finite element methods and demonstrated its efficacy in terms of simulation details and smoothness.Our method facilitates simulations on low-resolution or low-quality meshes.This feature opens up new possibilities to incorporate our method into the various interactive design, visualization, and geometry processing pipelines.It provides high-fidelity simulation results based on imperfect mesh data.Thanks to the second-order nature of our discretization, our method is particularly suited for producing authentic simulations on Origami and Kirigami designs, which were challenging for traditional thin-shell methods due to their mixed creases, wrinkles, low-resolution meshes, and anisotropic mesh elements.

Figure 2 :
Figure 2: Top: Approximating a saddle surface with firstorder (left) and second-order (right) finite elements whose triangulation is shown in the insets.Bottom: The basis function landscape (blue surfaces) and its location (black dots) of a single first-order or second-order element.

Figure 3 :
Figure 3: A visual explanation of the discrete mean curvature measure (Def.4.1) on a surface discretized with two planar (left) and curved (right) triangles.The triangulation in the (, ) parameter space is shown in the middle inset, and the dashed lines further subdivide each triangle to distribute average edge curvature measure.A Borel set (not shown) in R 3 intersects the surface onto two of such subdivided triangles (darker blue) and fully encloses the edge of interest, resulting in the area of the gray surface to be counted toward the discrete mean curvature measure.

Figure 4 :
Figure 4: Left: a sample triangulation of a domain Ω with a crease (red curve) splitting multiple triangular elements.Middle and right: One such triangle is split into two linear (middle) and curved (right) triangles accompanied by extra virtual nodes (black dots, respectively).The crease moves passively with the deformation from triangles.Borel set    whose intersection with x+1 is exactly x+1 (  (   )  )

Figure 5 :
Figure 5: Cloth falling on a sphere.Top row are obtained by first-order elements with respective resolutions of 4 × 4, 4 × 16, and 16 × 4 (left to right).Bottom row are obtained by second-order elements with respective resolutions of 2 × 2, 2 × 8, and 8 × 2, matching the number of DoFs of its counterpart.

Figure 6 :
Figure 6: Saddle: This figure demonstrated the ability of quadratic elements to alleviate kinematic locking with various young's modulus (E).The first two rows are linear elements with diagonal edges aligned with and orthogonal to the bending direction.The next two rows are quadratic elements with such edges.The linear elements are on a resolution of 30 × 30 × 2 triangles, and quadratic elements are of 10 × 10 × 2 triangles.The triangulations of the surfaces are visualized by red/white/grey triangles.

Figure 7 :
Figure 7: Curtain: Quadratic elements at 30 × 30 (60 × 60) resolution are compared with linear elements at 90 × 90 (180 × 180) resolution.The curtain swings under gravity with 6 evenly-spaced positions on the top moving horizontally.The close-ups of location where locking is apparent are shown in blue boxes.

Figure 8 :
Figure 8: Comparison with other higher-order elements: Comparison of our method (top) with RerShell (middle) and tr_shell01 (bottom).A square swings under gravity with one side fixed.Only our method presents reasonable result.The other two elements deform incorrectly at the bottom of the surface.

Figure 9 :
Figure 9: Large curtain: The left figure shows an intermediate frame from simulating the curtain with 200 × 200 × 2 triangles with textures.The right figure displays the same frame without texture to showcase the rich wrinkles.

Figure 10 :
Figure 10: Tori: A piece of cloth simulated with 459 elements is in contact with two solid tori (top-right inset) and gently slips.

Figure 11 :
Figure 11: Stairs: A Kirigami design simulated with 600 elements.Starting from a planer rest shape (top left), the paper was pulled by its creases with non-zero rest angles (middle and bottom left) to form a stair-like kirigami (right).

Figure 12 :
Figure 12: Bat: A bat wing composed of elastic rods and shells.The motion is controlled by rods representing "arms".See the video for more details.

Figure 13 :
Figure 13: Kirigami: A Kirigami design simulated with 32 × 32 elements.Starting from a planer rest shape (top right), the shell drooped under gravity.

Figure 14 :
Figure 14: Paper art: A paper-art example assembled by ten paper strips with creases is simulated using 572 elements.The left ends of the strips are fixed, and the right ends move leftward to force the strips to pop upward.The right rows are middle frames of the simulation.