Efficiently simulating thin objects such as cloth, paper, or skin remains significantly more challenging than simulating their rigid or volumetric counterparts, due to their kinematic complexity and nonlinear physics. Specifically, numerical discretizations of thin plates must overcome two challenges: the dramatic stiffness disparity between stretching and bending material forces, and the tendency for discretizations of thin surfaces to lock.
Stiffness scale-separation. The behavior of thin plates is governed by a combination of membrane forces, resisting in-plane stretching and shear of the material, and bending forces. These forces have dramatically different stiffnesses: stretching stiffness scales like the thickness of the material, whereas bending scales like thickness cubed; cloth, for example, is approximately thick, so that the stretching forces are million times stronger.
For sufficiently thin objects, the scale separation between these stiffnesses is so vast that stretching of the material is imperceptible; in this setting paying the price in element or time step size in order to resolve the stretching modes does not make sense. Instead, one can look at isometric kinematics, where constraints enforcing zero in-plane strain replace force-level resistance to stretching. However, such a strategy requires discrete surface kinematics that supports cleanly decoupling bending from stretching modes.
Membrane locking. A smooth surface can bend into a general developable surface isometrically; the same is not true for a triangle mesh. Regular equilateral meshes can isometrically bend about their three axes of symmetry only; irregular meshes cannot bend at all without distorting some of the triangle edges. As a consequence, any formulation of stretching forces based on triangle edge lengths will suffer from membrane locking: situations, such as holding a piece of cloth up by two corners as shown in Figure 1, left, where deformations that are bending-dominated in the continuous regime are stretching-dominated in the discrete regime (Quaglino, 2016), due to failure of bending to kinematically decouple from stretching. Note that unlike other discretization errors, locking is independent of mesh resolution: an irregular mesh with inextensible edges will not bend no matter how fine.
On the other side of the same coin, insufficiently constraining discrete inextensible plates, so that stretching-dominated deformations in the continuous regime are instead bending-dominated (or completely uncontrolled), allows equally undesirable spurious deformation modes.
We present a method for simulating isometric thin plates by combining two simple ideas: replacing the very stiff membrane forces by hard constraints, and doing so on a formulation of the strain tensor that is averaged over surface patches via moving least squares, to avoiding membrane locking. We can then simulate the material by adding a bending energy and enforcing the hard inextensibility constraints using standard methods for constrained Lagrangian dynamics, such as the method of Fast Projections (Goldenthal et al., 2007).
The resulting numerical method is simple to implement, and because stretching forces are replaced with hard constraints, it can simulate infinitesimally thin, inextensible materials without requiring very fine meshes or small time steps. Crucially, the averaged inextensibility constraints do not lock, and do not exhibit spurious modes, so that even very coarse, efficient simulations of thin materials undergoing large amounts of bending and crumpling give good qualitative results.
2. Related Work
There is a deep body of work in both computer graphics and finite element analysis on the simulation of thin plates and shells. Membrane locking is a well-known phenomenon for simulations involving low-order elements, and has been studied extensively by Quaglino (2012; 2016), who proposed several tests characterizing locking and taxonomized the locking behavior of a variety of triangle-mesh-based kinematics. English and Bridson (2008) suggested gluing triangles at edge midpoints rather than at vertices, which avoids locking but unfortunately suffers from spurious modes. Popular workarounds to locking include adaptive refinement of the mesh, as in the popular ArcSim (Narain et al., 2012) code and its extensions; avoiding triangles entirely and using higher-order elements; or compensating for locking by tweaking stiffness parameters on an ad-hoc per-example basis.
Treating membrane strain with constraints rather than forces has proven to be a powerful technique for reproducing characteristic wrinkle patterns in thin shell materials. Early application of this idea was used to introduce fine wrinkles where contact introduces compression of the materials (Provot, 1995; Bridson et al., 2003). Goldenthal et al. (2007) proposed a framework for constraint-based limiting of strain in the warp and weft directions on a quadrilateral mesh, and demonstrated that constraint-based inextensibility is significantly more efficient than force-based simulation of membrane strain, even when using implicit methods. Chen and Tang (2010) enforce isometry in a least-squares sense while also respecting collision constraints. Other methods for enforcing strain limits have been proposed, with support for constraining all elements of the strain tensor (Thomaszewski et al., 2009) and for boosting performance by applying constraints in a hierarchical manner (Wang et al., 2010), though these methods treat strain based on triangle deformations and will lock unless the material is sufficiently compliant when under the strain limit. Also related is the method of position-based dynamics (Müller et al., 2007; Stumpp et al., 2008), a stable and fast alternative to traditional physical simulation that replaces all forces (even the soft bending forces) with constraints.
Tension Field Theory
In tension-dominated problems, wrinkles perpendicular to the direction of tension have no effect on the coarse material shape, and so can be neglected; this idea is at the heart of tension field theory (Mansfield, 1964; Steigmann, 1990), which has been discretized in graphics for predicting the inflated shape of balloons (Skouras et al., 2014). The theory has also been applied to cloth, either via inequality constraints on edge lengths (Jin et al., 2017) or incorporation of wrinkle parameters into the kinematics (Quaglino, 2016).
Although most cloth solvers are mesh-based, meshless methods have also been explored (Yuan et al., 2008). They are particularly appealing for problems involving fracture; peridynamics (Silling, 2000) has had significant success in computer graphics (Levine et al., 2014; Chen et al., 2018; He et al., 2018; Xu et al., 2018) and has been extended to shells (Chowdhury et al., 2016). Similar in spirit are Elastons (Martin et al., 2010), a quadrature scheme for deformation energy based on unstructured sample points which can be used to simulate elastic bodies of arbitrary codimension in a unified way. Also related are the “F-bar” methods in finite elements for simulating incompressible volumes, which avoids locking by imposing area/volume constraints on patches of elements rather than on single quadrilaterals (de Souza Neto et al., 1996) or triangles (Neto et al., 2005). This idea has been applied to simulation of incompressible volumes (Ortiz-Bernardin et al., 2015; Boroomand and Khalilian, 2004), including in graphics (Irving et al., 2007).
3. Meshless Kinematics
In this section, we describe the meshless kinematics we use, and our notation. The heart of our method will be in Section 4, where we formulate isometry constraints on neighborhoods of the material. We assume we are given a surface whose rest state is flat. We adopt a standard meshless discretization of the plate, and sample at points of , and associate to each point a neighborhood of other sample points. We require that ; the neighboring points can be chosen based on a threshold distance away from the sample points on , graph distance on a user-provided input triangle mesh, etc.
We assume that locally, each neighborhood can be isometrically parameterized by a region of the plane (this parameterization might be given, in the case of cloth sewing patterns, or precomputed from by plane-fitting). We will write for the position of sample point on ; note that a sample point likely belongs to multiple neighborhoods and might have different material coordinates in each such neighborhood. Let be the embedded position of sample point in 3D. The set of
are then the degrees of freedom of the simulation. Finally we associate to each neighborhood a lumped mass(based on a given material density and the barycentric area of its associated sample point).
4. Isometry Constraints
The most obvious way to enforce isometry of a surface is to discretize it as a triangle mesh, and constrain each edge length to remain constant. However, such constraints are doomed to lock. Consider for instance that for a mesh with vertices and edges, there are kinematic degrees of freedom and constraints, and yet from the Euler characteristic formula, ; the edge-based isometry constraints are therefore expected to remove almost all of the cloth’s degrees of freedom. Contrast this situation with the smooth setting, where the cloth can deform into a rich variety of developable surfaces. Note that replacing the triangle mesh with quadrilaterals does not in any way solve the problem: one can then enforce isometry for only the quadrilateral edges (as in Goldenthal et al. (2007)), which allows non-isometric shear and stretching in the diagonal direction (see Figure 2), or also constrain the lengths of the diagonals, which essentially triangulates the mesh.
We propose instead to enforce isometry on the vertices: we will use moving least-squares to formulate a strain on each of the neighborhoods , and constraint the principle strains to be zero. We will then have only , rather than , isometry constraints, leaving leftover degrees of freedom for isometric bending modes.
Near a sample point , the deformation gradient of the thin plate, with respect to the local parameterization of , must satisfy
Of course, for points with more than two neighbors, this relation is overconstrained; we instead work with an averaged deformation gradient in , based on satisfying Equation 1 in the moving least-squares sense,
The averaged deformation gradient can then be expressed in closed form as
where each column of is one , each column of is one current displacement , and
Given current deformation gradient for a neighborhood of point , we can formulate the strain tensor,
is the identity matrix. Typically for dynamics we would then derive forces by applying a constitutive law to this strain; since we instead want to enforce inextensibilty as a hard constraint, we need to write down constraint functions specifying vanishing of the strain. There are several sets of constraints we could choose; unfortunately, all are nonlinear. We use the pair
Notice in particular that both constraints have zero as a regular value: neither gradient vanishes when , a crucial requirement for the numerical stability of techniques like Fast Projections that are based on the method of Lagrange multipliers. Intuitively, the trace of the strain tensor corresponds to the sum of the squared principle stretches, which is in a sense the averaged squared distance to the neighboring points. The determinant of the strain tensor is the product of the squared principle stretches and measures a notion of squared neighborhood area. Both constraints together imply that the stress tensor, averaged in each neighborhood, is the zero matrix, equivalent to isometry of the neighborhood.
We mention in passing some alternative sets of constraints we tried, and abandoned: perhaps the most obvious is to directly constrain each entry of , which amounts to three constraints per neighborhood since strain is symmetric. However, this formulation bloats the set of constraints by 50%, while introducing rank deficiency in the constraint gradients, leading to poor performance. One might also consider constraining the trace and determinant of the strain tensor, rather than of —the problem with this approach is that the function has a critical point at the strain-free state, which as mentioned above makes the constraint unsuitable for projection.
Relation of Locking to Constraint Jacobian Rank
We stress that the locking phenomena described in the introduction, and the first paragraph of this section, are entirely unrelated to linear dependence or independence of the constraint gradients. Obviously, unsatisfiable constraints cannot be enforced. But note that constraining the length of each edge in a triangle mesh locks the mesh, despite these constraints being linearly independent. Linearly independent constraint gradients are insufficient to prevent locking, nor will redundant (dependent) constraint gradients cause locking.
Our constraints are always satisfiable (assuming that the initial configuration is an isometric embedding of the plate) but the constraint gradients are not guaranteed to be independent: for instance, when , the trace and determinant constraints have identical gradient. This linear dependence is expected, from the geometry of bending: a flat neighborhood has two isometric deformation modes, since it can bend in any direction, but once the neighborhood begins to bend in one particular direction and symmetry is broken, only one isometric mode remains.
5. Time Integration
We enforce the isometry constraints (2) at every time step using the method of Fast Projections (Goldenthal et al., 2007), with implicit time integration of bending and external forces. Any bending model can be used, though we note that for materials which are rest flat (including cloth, paper, etc), the simple and efficient quadratic bending model (Bergou et al., 2006) is quite attractive, since isometry of the material is precisely the assumption required for validity of the model:
where is the system mass matrix, the Laplace-Beltrami operator (computed from a meshless discretization (Petronetto et al., 2013), or from a triangulation of the input surface), and is a bending stiffness proportional to the Young’s modulus and cubed thickness of the material. Note that the quadratic bending energy, true to its name, has constant Hessian which can be prefactored, so that the performance bottleneck of each time integration step is enforcing the isometry constraints.
In most of our experiments, we used a time step size of seconds and a constraint tolerance of 0.1; Fast Projections typically converges in one to three iterations. However, Fast Projections does not perform a true projection onto the constraint manifold, and comes with no guarantees; we found that regularizing each projection by the geometric stiffness matrix of the constraints (Tournier et al., 2015) improves performance. In cases where Fast Projections fails to decrease the constraint residual after a few iterations, we switch to true projection using the augmented Lagrangian method; in practice this fallback was only needed when high-energy impacts cause large violations in the constraints within a single time step.
We perform a sequence of tests demonstrating that the constraints (2) neither lock nor exhibit spurious modes. In Figure 2, we pin a square sheet of cloth at two corners, and let the sheet relax under gravity. We simulate this experiment using our isometry constraints for different constraint tolerances, for different resolutions of sampling points, and for both ordered and disordered sampling of the cloth. Our simulations consistently give comparable results. Similarly, for tests that drape cloth on rigid bodies (Figure 3), our method produces qualitatively similar results, without locking, even for very coarse discretizations of the cloth. Statistics for these and other experiments described in this section can be found in Table 1.
Effect of Parameters
As we decrease the constraint tolerance as shown in the left column in Figure 2, the sagging in the middle of the cloth decreases as expected, since the isometry constraints enforce that the distance between the two pinned corners is identical to that of the flat rest shape.
Comparison to Related Work
We compare to force-based shell methods (Grinspun et al., 2006), using the same bending stiffness, on the third column of Figure 2, where enforcing isometry via high stiffness locks the cloth. Lastly, the same experiment was performed using the constraints proposed by Goldenthal et al. (2007). Due to the fact that the edge length constraints on a quad method do not control the length of each quad diagonal, the cloth sags significantly.
We also perform a test where the cloth deformation is driven by shear. We pin one side of a piece of cloth, and attach a heavy mass to a point on the other. The cloth stretches taut along diagonal lines of tension, but does not lock, and in contrast to the Goldenthal et al. constraints on the principal strain along the axes of a quadrilateral grid, obeys the isometry constraints, which we verify by plotting the change in extrinsic distance from the top left corner to every other points in the material; see Figure 4.
Our method, unlike thin plate formulations based on tension field theory, works equally well in problems involving pure compression, such as the crumpling of an elastic cylinder. We replicate the cascade of diamond buckling patterns observed by Martin et al. (2010), and although only our high-resolution simulation can resolve the fine-scale pattern at the beginning of the simulation, it correctly reproduces the highly-buckled final state.
Spinning tank top
We show the dynamics of a tank top draping over a spinning hanger in Figure 6. Our method resolves the wrinkles generated during the motion.
One important feature of our method is that the bending stiffness of the simulated cloth can be controlled, while preserving isometry. For most thin shell models, extensive tuning of the hyper-parameters is often required to achieve the balance between avoiding locking and maintaining isometry of the material. As shown in Figure 7, just by adjusting the bending stiffness, our method produces cloth with different wrinkling patterns without stretching.
Several papers (Rémillard and Kry, 2013; Li and Kry, 2014) describe how to simulate wrinkling of skin via one-way coupling of thin plate simulations to a volumetric elastic substrate. The coupling is enforced via a sparse set of average-position constraints sampled on the skin surface. This setup is ideal for our method, since (i) the behavior of the upper skin layers is governed by bending and by the coupling constraints, so that replacing in-plane strain of the skin with isometry of these fine layers is justified; (ii) the coupling constraints can be trivially incorporated by including them as extra terms during Fast Projections. Figure 8 shows the behavior of our simulation when the substrate is compressed in one and two directions.
- A quadratic bending model for inextensible surfaces. In Proceedings of the Fourth Eurographics Symposium on Geometry Processing, SGP ’06, pp. 227–230. External Links: Cited by: §5.
- On using linear elements in incompressible plane strain problems: a simple edge based approach for triangles. International Journal for Numerical Methods in Engineering 61 (10), pp. 1710–1740. Cited by: §2.
- Projective dynamics: fusing constraint projections for fast simulation. ACM Trans. Graph. 33, pp. 154:1–154:11. Cited by: §2.
- Simulation of clothing with folds and wrinkles. In Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’03, pp. 28–36. External Links: Cited by: §2.
- A fully geometric approach for developable cloth deformation simulation. The Visual Computer 26 (6), pp. 853–863. External Links: Cited by: §2.
- Peridynamics-based fracture animation for elastoplastic solids. Computer Graphics Forum 37 (1), pp. 112–124. Cited by: §2.
- A peridynamic theory for linear elastic shells. International Journal of Solids and Structures 84, pp. 110 – 132. External Links: Cited by: §2.
- Design of simple low order finite elements for large strain analysis of nearly incompressible solids. International Journal of Solids and Structures 33 (20), pp. 3277 – 3296. External Links: Cited by: §2.
- FEPR: fast energy projection for real-time simulation of deformable objects. ACM Transactions on Graphics (TOG) 37 (4), pp. 1–12. Cited by: §2.
- Animating developable surfaces using nonconforming elements. ACM Trans. Graph. 27 (3), pp. 66:1–66:5. External Links: Cited by: §2.
- Efficient simulation of inextensible cloth. ACM Transactions on Graphics (Proceedings of SIGGRAPH 2007) 26 (3), pp. to appear. Cited by: §1, Figure 2, §2, §2, §4, §5, §6.
- Computing discrete shape operators on general meshes. Computer Graphics Forum (Eurographics) 25 (3), pp. 547–556. External Links: Cited by: Figure 2, §6.
- Projective peridynamics for modeling versatile elastoplastic materials. IEEE Transactions on Visualization and Computer Graphics 24 (9), pp. 2589–2599. External Links: Cited by: §2.
- Volume conserving finite element simulations of deformable models. ACM Trans. Graph. 26 (3). External Links: Cited by: §2.
- Inequality cloth. In SCA ’17 Proceedings of the ACM SIGGRAPH / Eurographics Symposium on Computer Animation, Los Angeles, California, pp. 1–10. Cited by: §2.
- A peridynamic perspective on spring-mass fracture. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation, SCA ’14, pp. 47–55. Cited by: §2.
- Multi-layer skin simulation with adaptive constraints. In Proceedings of the Seventh International Conference on Motion in Games, MIG ’14, New York, NY, USA, pp. 171–176. External Links: Cited by: §6.
- The bending and stretching of plates. Vol. 6, Pergamon, Amsterdam, The Netherlands. Cited by: §2.
- Unified simulation of elastic rods, shells, and solids. In ACM SIGGRAPH 2010 Papers, SIGGRAPH ’10, New York, NY, USA, pp. 39:1–39:10. External Links: Cited by: §2, §6.
- Position based dynamics. Journal of Visual Communication and Image Representation 18 (2), pp. 109 – 118. External Links: Cited by: §2.
- Adaptive anisotropic remeshing for cloth simulation. ACM Trans. Graph. 31 (6), pp. 152:1–152:10. External Links: Cited by: §2.
- F-bar-based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. part i: formulation and benchmarking. International Journal for Numerical Methods in Engineering 62 (3), pp. 353–383. Cited by: §2.
- Improved robustness for nearly-incompressible large deformation meshfree simulations on Delaunay tessellations. Computer Methods in Applied Mechanics and Engineering 293, pp. 348 – 374. External Links: Cited by: §2.
- Mesh-free discrete laplace-beltrami operator. Comput. Graph. Forum 32 (6), pp. 214–226. External Links: Cited by: §5.
- Deformation constraints in a mass-spring model to describe rigid cloth behavior. Graphics Interface 23 (19), pp. 147–154. Cited by: §2.
- Membrane locking in discrete shell theories. Ph.D. Thesis, Georg-August-Universität Göttingen. Cited by: §2.
- A framework for creating low-order shell elements free of membrane locking. International Journal for Numerical Methods in Engineering 108 (1), pp. 55–75. Cited by: 2nd item, §2, §2.
- Embedded thin shells for wrinkle simulation. ACM Trans. Graph. 32, pp. 50:1–50:8. Cited by: §6.
- Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids 48 (1), pp. 175 – 209. External Links: Cited by: §2.
- Designing inflatable structures. ACM Trans. Graph. 33 (4), pp. 63:1–63:10. External Links: Cited by: §2.
- Tension-field theory. Proceedings of the Royal Society of London A 429, pp. 141–173. Cited by: §2.
- A Geometric Deformation Model for Stable Cloth Simulation. In Workshop in Virtual Reality Interactions and Physical Simulation ”VRIPHYS” (2008), F. Faure and M. Teschner (Eds.), Grenoble, France, pp. 1–8. External Links: Cited by: §2.
- Continuum-based strain limiting. Computer Graphics Forum 28 (2), pp. 569–576. Cited by: §2.
- Stable constrained dynamics. ACM Trans. Graph. 34 (4), pp. 132:1–132:10. External Links: Cited by: §5.
- Multi-resolution isotropic strain limiting. ACM Transactions on Graphics 29 (6), pp. 156:1–10. Note: Proceedings of ACM SIGGRAPH Asia 2010, Seoul, South Korea Cited by: §2.
- Reformulating hyperelastic materials with peridynamic modeling. Computer Graphics Forum 37, pp. 121–130. Cited by: §2.
- Application of meshless local petrov-galerkin (mlpg) method in cloth simulation. CMES - Computer Modeling in Engineering and Sciences 35, pp. 133–154. Cited by: §2.