Efficient and robust strain limiting and treatment of simultaneous collisions with semidefinite programming

We present an efficient and robust method which performs well for both strain limiting and treatment of simultaneous collisions. Our method formulates strain constraints and collision constraints as a serial of linear matrix inequalities (LMIs) and linear polynomial inequalities (LPIs), and solves an optimization problem with standard convex semidefinite programming solvers. When performing strain limiting, our method acts on strain tensors to constrain the singular values of the deformation gradient matrix in a specified interval. Our method can be applied to both triangular surface meshes and tetrahedral volume meshes. Compared with prior strain limiting methods, our method converges much faster and guarantees triangle flipping does not occur when applied to a triangular mesh. When performing treatment of simultaneous collisions, our method eliminates all detected collisions during each iteration, leading to higher efficiency and faster convergence than prior collision treatment methods.


Strain limiting
In the real world, many materials, such as cloth and animal tissues, can only deform to a limited degree.Although compliant to small deformations, they are highly resistant to deformations larger than some threshold.For cloth, its structure of fibers and threads easily accommodates small-tomoderate amounts of stretching, but once the structural slack has been taken up, resistance to further stretching increases dramatically [1].When developing characteristic dynamics models to simulate such materials, many researchers have taken account of this biphasic property as a key to developing wrinkle patterns observed in many fabrics.Similarly, animal tissues, such as skin or relaxed muscles, are also compliant to small strains but very tough and resistant to larger ones.
Unfortunately, most simulation methods perform poorly for materials with highly non-compliant constitutive regimes.Standard finite-element methods, including mass-spring systems, model strong resistance using large material coefficients, leading to integration problems over time [1].Explicit integration methods require very small time steps to stay stable and avoid divergence.Implicit methods can maintain stability using relatively large time steps, but may converge slowly, and suffer from high residuals or excessive damping.To prohibit excessive extensibility in simulation, most dynamic models use projection to enforce a hard limit on large strains, i.e., strain limiting.
Many strain limiting methods have been proposed.Anisotropic strain limiting methods [2] in cloth simulation limit strains of edges.Isotropic methods [1,3] act on the strain tensors and limit the singular values of the deformation gradient matrices of triangles.Our method also performs isotropic strain limiting for both triangular surface meshes and tetrahedral volume meshes.Compared to prior methods, our method has faster convergence and produces better strain limited triangular meshes without any flipped triangles.

Treatment of collisons
Collision handling is another vital research area in computer animation.Deformable bodies naturally bring about large numbers of collisions varying in strength from resting contact to high speed impact.Collisions between soft and thin sheets, such as cloth and paper, are especially difficult to handle.There are two phases in collision handling, collision detection and collision response.Collision detection is an important component of physically based simulation systems targeting cloth and hair, and finite-element simulation.Even a single undetected collision can lead to simulation failure.Continuous collision detection (CCD) algorithms are widely used in cloth simulation to detect collisions between cloths, cloths and rigid bodies, and self-collisions of the same cloth.Robust collision response is also vital for cloth and shell simulation.Not only must collisions be prevented, but the response must be physically plausible.Unsuitable treatment of collisions may lead to divergence in collision handling and simulation failure.While the handling of individual collisions is well understood, simultaneous collisions can halt existing methods.When facing a cluster of interacting simultaneous collisions, a rigid impact zone (RIZ) approach as proposed by Bridson et al. [4] may be adopted to prevent collisions, but this also eliminates all relative tangential velocity in a physically implausible way.Another approach based on inelastic projection impact zone (IIZ), proposed by Harmon et al. [5], formulates collisions as linear constraints and then solves a nearest projection optimization problem.Our method is similar to IIZ, but we impose an additional determinant inequality constraint on each collision which ensures we can eliminate every detected collision in only one iteration.

Semidefinite programming
The notation A 0 implies that A is a symmetric, positive semidefinite (PSD) matrix.
Such an expression is a linear matrix inequality (LMI) [6].In the same manner, A B implies that A − B is PSD and thus for a scalar c ∈ R, the equation S cI implies that the eigenvalues of S are greater than or equal to c.
A semidefinite program (SDP) is a convex optimization problem formulated with LMI constraints and a linear objective.Semidefinite programming unifies several standard problems.It is easy to formulate any linear program, convex quadratic program, and second-order cone program as an SDP.SDP solvers are still not as mature as more classical optimization methods, and have higher time complexity.However, they are already efficient enough for many applications in computer graphics.Semidefinite programs are as easy to solve as linear programs.Most interior-point methods have been generalized to semidefinite programs.

Main results
In this paper, we present an efficient and robust method to perform both strain limiting and treatment of simultaneous collisions.
• We perform strain limiting of triangular meshes in the fashion of strain limiting of tetrahedral meshes.
The benefit is that it can prevent triangle flipping.• We impose an additional determinant constraint on collision response, which can ensure the detected collisions to be eliminated by only one iteration.• We formulate strain constraints and collision constraints as a serial of LMIs and LPIs, then solve a projection optimization problem with semidefinite programming.

Organization
The rest of the paper is organized as follows.
We survey prior work on strain limiting and collision handling in Section 2. We present a position based projection optimization problem with strain and collision constraints in Section 3. We present our semidefinite programming solution to our optimization problem in Section 4; implementation is considered in Section 5. We highlight the results and performance of our method in Section 6.

Related work
Provot [7] introduced strain limiting as a technique for stably modeling stiff springs by imposing constraints on the maximum and minimum allowed strain of each link.Subsequently, many extensions of this technique have been developed.Bridson et al. [4] applied momentum-conserving impulses to the velocities to ensure that all springs are deformed by a maximum of 10% from their rest lengths at the end of each time step.Goldenthal et al. [8] proposed an efficient constrained Lagrangian method for modeling inextensible spring networks.For triangle meshes, English and Bridson [9] used nonconforming elements to allow more degrees of freedom of a strain limited triangle to model inextensible cloth.Thomaszewski et al. [10] presented a continuumbased technique that independently constrains the three components of the strain tensor of a triangle.A technique for isotropic strain limiting was proposed by Wang et al. [1], who also introduced a multiresolution approach for enforcing these constraints.Narain et al. [3] posed strain limiting as a nonlinear optimization problem and used an augmented Lagrangian method to solve it.Collision handling plays a significant role in physically based simulation.Continuous collision detection (CCD) methods are widely used to detect collisions and intersections in cloth, hair, and thin sheet simulations.Brochu et al. [11] proposed a volume-based geometric predictor, and Tang et al. [12] proposed a Bernstein sign classification (BSC) based predictor, both of which can provide exact collision results by taking advantage of exact geometric arithmetic.Wang [13] introduced error analysis into a traditional cubic solver for CCD to achieve conservative but acceptable collision results.Detected collisions also need to be handled properly.Bridson et al. [4] applied repulsion impulses to collision elements to remove collisions.To handle clusters of interacting simultaneous collisions that repulsion impulses are not able to deal with, Provot [14] and Bridson et al. [4] used a rigid impact zone (RIZ) approach to prevent collisions.Huh et al. [15] divided the impact zone into clusters to partially alleviate the rigidification.Tsiknis [16] considered shearing modes of the impact zone.Harmon et al. [5] proposed a new inelastic projection impact zone (IIZ) method to better deal with simultaneous collisions.Tang et al. [17] presented a method to compute continuous penalty forces to determine collision responses between rigid and deformable models bounded by triangle meshes.

Position based projection
Many approaches to the simulation of dynamic systems in computer graphics are force based methods.In physically based simulation, strain limiting and collision response are used as remedies when excessive deformations or collisions appear after numerical time integration.Force based methods [4,5,18] act on velocities and then use the updated velocities to find final positions meeting with various kinds of constraints, such as strain constraints, collision constraints, and position constraints (e.g., fixed points).It usually needs many iterations to determine final good-quality positions, and requires relatively small time steps to keep the simulation system stable.In contrast, Müller et al. [19] proposed position based dynamics (PBD), which acts directly on positions to get a well constrained simulation configuration.Many constraints, such as strain and collision constraints, are very easy to handle by projecting points to valid locations in PBD.It is also very stable and allows simulations to take relatively large time steps.
Strain limiting and collision response can both be viewed as finding the closest projections to meshes which are correctly strain limited and collision-free.Given a mesh with m vertices, we can get its candidate positions [q 1 , . . ., q m ] T = q ∈ R 3m after numerical time integration.Let q ∈ R 3m be the well constrained positions of the mesh.The objective of the position based projection optimization problem is defined as following: where M is the mass matrix and ||.|| F denotes the Frobenius norm.Strain constraints and collision constraints are imposed respectively when performing strain limiting and collision response.We adopt semidefinite programming to solve the optimization problem stably and efficiently, as presented in Section 4.

Strain constraints
A well strain-limited surface mesh is one in which the strains of all faces lie within the user-specified strain limits [s min , s max ].We use [u 1 , . . ., u m ] T = u ∈ R 2m to denote the material coordinates of a surface mesh in material space, with u i ∈ R 2 .The deformation gradient of a triangle face is F = D q D −1 u , where We diagonalize F into F = U SV T using singular value decomposition (SVD); U and V are orthogonal matrices and S is a 3 × 2 matrix with nonnegative values on the diagonal, which are the two principle strains (s 1 , s 2 ) of the triangle face.So the strain constraints for the triangle face are s min s(q i , q j , q k ) = (s 1 , s 2 ) s max (2) i.e., s min s 1 , s 2 s max .
Strain limiting for a volume mesh is performed in a similar way.The deformation gradient of a tetrahedron [q i , q j , q k , q l ] is also written as F .Unlike the surface mesh case, now D q = [q j −q i , q k − q i , q l − q i ] and Finding the SVD of F gives three nonzero singular values (s 1 , s 2 , s 3 ), so the strain constraint for a tetrahedron is s min s(q i , q j , q k , q l ) = (s 1 , s 2 , s 3 ) s max (3) To unify strain limiting for surface meshes and volume meshes, we extend the 2D material coordinates of surface meshes to 3D by simply setting the third value to zero, i.e., u i = [u 1 , u 2 , 0] T .In addition, we introduce an auxiliary vertex p l for each triangle [q i , q j , q k ] of the surface mesh, as shown in Fig. 1(a).Its material coordinate is u p l = u i + [0, 0, 1] T , and its initial position in world space is p l = q i + n where n is the unit normal vector of the triangle.This transforms a triangular surface into a tetrahedral volume mesh, as shown in Fig. 1(b).This allows Eq. ( 2) to be updated to s min s(q i , q j , q k , p l ) = (s 1 , s 2 , s 3 ) s max (4) Moreover, we add a matrix determinant inequality constraint as below to make sure that strain limited triangles do not flip; a comparison with the previous method without this constraint is shown in Fig. 2. d(q i , q j , q k , p l ) = det([q j −q i , q k −q i , p l −q i ]) > 0 (5) This indicates that the matrix is orientation preserving and the volume of the tetrahedron  constructed by a triangle and its auxiliary vertex is always positive.The strain constraints for multiple triangles are formulated as follows: s min s(q, p) s max (6a) d(q,p) > 0 (6b) If there are w triangles, then [p 1 , ..., p w ] T = p ∈ R 3w , s ∈ R 3w and d ∈ R w .

Collision constraints
When performing treatment of collisions, two elementary kinds of primitive, vertex-face (VF) pairs and edge-edge (EE) pairs, need to be tackled properly.Each pair consists of four vertices, [q i , q j , q k , q l ].For a VF pair, q i is a vertex and [q j , q k , q l ] represents a triangle face.For an EE pair, [q i , q j ] and [q k , q l ] represent the two edges.Two distance constraint functions for a VF and an EE pair are defined as respectively in Ref. [5], as follows: c VF (q i , q j , q k , q l ) = n•[α 0 q i −(α 1 q j +α 2 q k +α 3 q l )] c EE (q i , q j , q k , q l ) = n•[(α 0 q i +α 1 q j )−(α 2 q k +α 3 q l )] (7) where n is a normal vector and α 0 , α 1 , α 2 , α 3 ∈ [0, 1] are parameters.For a VF pair, n is the normal of the triangle face, α 0 = 1, α 1 + α 2 + α 3 = 1.For an EE pair, n is the cross product of the two edges, A pair is collision-free provided that c > 0. At an intermediate time in a time step, a pair is colliding if c 0, as shown in Fig. 3; the projection of vector q c q i onto the normal vector n is negative.So the collision constraint for a pair is c(q i , q j , q k , q l ) > 0 (8) Because the normal n dependends on the positions of the four vertices, the collision constraints are nonlinear.Harmon et al. [5] just sets n to be the normal at the time of collision, so the constraint function c becomes linear.This corresponds to removing collisions by just modifying the positions along the normal vector, which essentially conforms to the laws of physics.Although the linear constraint in Eq. ( 8) is met, it cannot guarantee the elimination of collisions at the end of a time step because the vertex may still be on the negative side of the triangle face for a VF pair, like in the case shown in Fig. 3(b).
To simplify c(q), we also make the same choice as Harmon et al. [5].In the meantime, we impose a similar constraint to Eq. ( 5) on collision: which indicates the volume of the tetrahedron consisting of the four vertices of a collision pair is positive, i.e., a vertex is always on the positive side of its opposite triangle in the tetrahedron.Combining Inequality (9) with Inequality (8) simplifies the processing of collisions and ensures that collisions are eliminated.Figure 3(c) shows a result that imposes constraints with both Inequalities ( 8) and ( 9), ensuring that the vertex is on the positive side of the Fig. 3 Collision response for a VF pair.A VF collision pair consists of vertex qi and triangle ∆qj q k q l .n is the normal vector at collision time and qc = α1qj + α2q k + α3q l is the collision point, where α1, α2, α3 are barycentric coordinates.N is the normal vector of the triangle after the collision response.

Positional constraints
Positional constraints, such as fixed-points and gluing constraints, are very common in cloth simulation.

Fixed points
A fixed-point constraint for vertex i can be formulated as in which x i is a fixed point in world space.The fixed-point constraint function can be reformulated to enforce multiple fixed-point constraints using: [x 1 , ..., x m ] T = x ∈ R 3m is a column vector and P is an m × m diagonal block matrix, where each block is a 3×3 matrix.If vertex i is constrained, the ith entry x i is a user-defined vector in world space and the ith entry on the diagonal of P is an identity matrix, i.e., P ii = I.The other entries of x and P are zero.

Glue
A gluing constraint for two vertices i and j can be formulated as d min g(q i , q j ) = n T • (q i − q j ) d max (13) where n ∈ R 3 is a unit normal and [d min , d max ] is a distance interval.The glue constraint function can also be reformulated to enforce multiple glue constraints using: If vertices i and j are glued together, then n j = −n i and the other entries are zero.

Solution by semidefinite programming
Strain constraints (i.e., Inequalities (4) and ( 5)) and collision constraints (i.e., Inequality ( 9)) are nonlinear and non-convex, which makes the projection optimization problem complicated and difficult to solve.Narain et al. [3] adapted an augmented Lagrangian method to solve the problem with strain constraints.However, it converges slowly and is not robust in practice, as shown in Figs. 2 and 4. Instead, to solve our problem, we adopt the method proposed by Kovalsky et al. [20], which can control the singular values of a square matrix to lie within a positive interval [γ, Γ ] by use of semidefinite programming (SDP).
A meta-problem is defined as follows to control the singular values of an arbitrary square matrix.
where A is an n × n matrix, s min (A) and s max (A) are the minimum and maximum singular values respectively, and f (A, s min (A), s max (A)) is a convex objective function.Obviously, the feasible set of the meta-problem is non-linear and non-convex.It is difficult to solve it using traditional optimization methods.Recently, Kovalsky et al. [20] proposed a simple and efficient method to solve the metaproblem, which reformulates the non-linear and nonconvex constraints in Inequalties (15b) and (15c) as two linear matrix inequalities respectively.Inequality (15b) can be replaced equivalently by the following LMI: Furthermore, a family of maximal convex subsets is found to cover the entire set defined by Inequality (15c).Each maximal convex subset is defined by an LMI of the form: To find a global solution, the procedure rotates the current maximal convex subset to the next iteratively.
Positional constraints such as fixed-point and glue constraints are easy to handle because they are linear.We just detail how to deal with strain constraints and collision constraints.The problem of strain limiting for a single triangle or tetrahedron, or collision for a single pair, is very similar to the meta-problem.Strain and collision constraints can be reformulated as LMIs and LPIs, allowing us to take advantage of standard convex SDP solvers to solve strain limiting and collision problems in physically based simulation.

Strain limiting
Strain limiting for a single triangle can be defined via the following projection optimization problem with LMI constraints: iterations respectively.The table at top right shows the number of iterations taken by these two methods to converge.It is clear that our method converges faster (in one iteration) than Narain et al.'s method [3].The disadvantage of our method is that each iteration takes longer to perform.
where F is the deformation gradient matrix of the triangle defined in Section 3.1.After transforming the triangle into a tetrahedron, F is a square matrix.
Because F is just a linear transformation of q, we can easily apply the method to our strain limiting problem.Having shown how strain limiting for a single triangle or tetrahedron is done, it is easy to extend Eq. ( 18) to the case for multiple triangles and tetrahedra.

Collision response
For a single pair involved in a collision, the problem can be defined as a projection optimization problem with both LMI and LPI constraints: where q = [q i , q j , q k , q l ] represents the collision pair, and q∆ = [q j −q i , q k −q i , q l −q i ] is a 3×3 matrix which is a linear transformation of q.Inequality (19b) is linear so that we can also apply this method to the collision problem.Furthermore, scenarios with complex simultaneous collisions are also very easy to handle.

Local strain limiting
There are two manners in which we can perform strain limiting for a deformable mesh, global and local.In global strain limiting, the candidate positions of all vertices in a mesh are constrained simultaneously.Global strain limiting is very simple, but relatively slow.
In local strain limiting, we detect triangles which violate strain limits.Correlated triangles which share vertices are put into an SL zone which represents the regions where strain limiting is needed, just like an impact zone which is widely used in collision response methods to deal with complex scenes with simultaneous collisions.We take each SL zone as a unit when dealing with regions which violate strain constraints.For each SL zone, a projection optimization problem is solved to project the region to the nearest location which satisfies the strain constraints.Local strain limiting needs to detect triangles violating strain limits and handle them iteratively until no triangles are detected.Each detected triangle is a smallest SL zone.Correlated SL zones are merged with each other.The extreme case is that the entire mesh is covered by a sinlge SL zone, which is equivalent to global strain limiting.
Handling an SL zone may make correctly strain limited triangles become no longer strain limited, necessitating another iteration.The worst case is that no longer strain limited triangles may appear one by one, causing slow local strain limiting convergence, and taking much time.To accelerate convergence, we extend each SL zone by merging it with its one ring of neighbouring triangles, as shown in Fig. 5.This approach makes SL zones enlarge more quickly, contributing to faster convergence of local strain limiting.
When applying global strain limiting in physically based simulations, the internal energies stored in meshes propagate faster and better than in local strain limiting.The disadvantage is that it takes relatively more time to solve a big optimization problem, especially when meshes are generally already correctly strain limited.In contrast, local strain limiting is very fast when meshes are correctly strain limited.As SL zones expand, it takes more and more iterations to ensure the meshes retain good qualities.Additionally, detecting no longer strain limited triangles in each iteration is expensive.Thus, global strain limiting is more suitable for simulations using large time steps where large deformations Fig. 5 Local strain limiting.The region is an SL zone.Before performing strain limiting for this SL zone, we extend it by merging with its one ring neighbour triangles, i.e., the light blue region in the left.Then we perform strain limiting for the extended bigger SL zone.This helps local strain limiting converge faster.
appear frequently, while local strain limiting is more suitable to simulations using small time steps.

SDP optimization
In our method, we have to solve the optimization problem defined in Eq. ( 20), in which the objective is a quadratic function: subject to LMI and LPI constraints (20) To solve the problem with an SDP solver, we reformulate the problem as the following equivalent optimization problem with linear objective: min t (21a) LMI and LPI constraints (21d) in which z and t are auxiliary variables.Inequality (21b) is a convex conic quadratic constraint and Inequality (21c) is a linear constraint.The new optimization problem is easy to solve with a standard convex SDP solver.

Results
We reformulate strain constraints and collision constraints as a series of linear matrix inequalities and linear polynomial inequalities in the projection optimization problem.The transformed problem is easy to solve using standard convex semidefinite programming solvers; we use the one in Mosek [21].
Compared to Narain's strain limiting method [3], our method converges faster and can prevent triangle flipping when performing strain limiting for triangular meshes.Figure 6 shows that our method performs well in strain limiting for tetrahedral meshes.However, Narain et al.'s method [3] takes Fig. 6 Strain limiting for a tetrahedral volume mesh.Bottom left: original un-deformed mesh.Top: the mesh is stretched to twice its original length and compressed to half its original height.Bottom right: the mesh produced by our method has well limited strain and is closest to the deformed mesh.
less time in each strain limiting iteration than our method.
Compared with Harmon et al.'s collision response method [5], our method takes fewer iterations and less time to converge, making the collision handling process faster.

Performance
We evaluated the performance of our method with various cloth simulation benchmarks, as shown in Figs. 7, 8, and 9, using a standard PC (3.4 GHz Intel i7-4770 CPU, 8 GB RAM, 64-bit Windows 7, NVIDIA GeForce GTX 780 GPU).This includes a C++ implementation of strain limiting and collision response based on Mosek's semidefinite programming solver.Figure 9 illustrates results of our method when performing strain limiting on cloth simulations, as well as a comparison with a prior strain limiting method.Table 1 highlights the performance of our method when computing collision responses on different benchmarks, as shown in Fig. 7 A moving ball hits a hanging cloth.The cloth mesh has 8.2k triangles; the strains of each triangle are limited to lie in [0.95, 1.05].This benchmark uses our method to perform strain limiting and collision response.Fig. 9 Strain limiting benchmarks.We use two cloth meshes with different resolutions (2k and 8.2k triangles respectively) to demonstrate the difference between our strain limiting method and Narain et al.'s method during cloth simulation.We limit the deformation of these meshes using different strain limits, [0.9, 1.1] allowing 10% deformation at most compared with the rest mesh, [0.95, 1.05] allowing 5% deformation at most, and [0.99, 1.01] allowing 1% deformation at most, respectively.The cloth exhibits more detailed wrinkles when using a higher resolution mesh.Comparing our method with Narain et al.'s method using the same meshes and the same strain limits, it is clear that meshes generated by our method are better strain limited.In contrast, meshes generated by Narain et al.'s method are more loose only after 10 iterations.After 35 iterationes, Narain et al.'s meshes in the third row become tighter and closer to our meshes in the first row.In our method, we replace the linear equations by linear inequalities and add additional determinant inequality constraints.This ensures that detected collisions are eliminated after only one iteration, as shown in Fig. 3.

Analysis
Strain limiting and treatment of collisions are two important processes in physically based simulation, particularly cloth and hair simulations.We have presented an efficient and robust method which can deal with both of them well.There are several advantages of our method.When performing strain limiting, we transform a triangle into a tetrahedron; our method applies to both triangular surface meshes and tetrahedral volume meshes.Additionally, our method can ensure the volume of a tetrahedron to be positive preventing triangle flipping during strain limiting.Compared with prior strain limiting methods, our method converges faster, although our method takes more time in each iteration when performing global strain limiting.Strain limiting for many triangles produces many low-dimensional LMI constraints, many more than the number of variables.Thus, standard SDP solvers may be non-optimal and need more time to find a optimal solution.When handling simultaneous collisions, our method eliminates all detected collisions during every iteration, which contributes to higher efficiency and faster convergence than prior collision handling methods.

Limitations, conclusions, and future work
We have presented an efficient and robust method which performs well both for strain limiting and handling simultaneous collisions.In our method, strain constraints and collision constraints are reformulated as a seriesl of linear polynomial inequalities (LPIs) and linear matrix inequalities (LMIs).Our method solves a projection optimization problem with Mosek's standard semidefinite programming solver.We have tested our method on some cloth simulation benchmarks and highlighted its benefits compared to prior strain limiting methods and collision response methods.
Our method has a few limitations.When performing strain limiting for a high-resolution mesh in global fashion, our method takes more time than Narain et al.'s method [3].When combining strain constraints with collisions constraints, our method may be unstable when many collisions occur simultaneously.In the future, it is very possible and indeed essential to optimize our method to make it faster when performing strain limiting for a highresolution mesh.A more efficient SDP solver may also help to solve the global strain limiting problem faster.
Fund of Ministry of Education of China (No. 20130101110133).Ruofeng Tong is partly supported by National Natural Science Foundation of China (No. 61572424).

Fig. 1
Fig.1Transforming a triangular mesh into a tetrahedral mesh.(a) An auxiliary vertex p is added so that the vector q0p is the unit normal vector of the triangle in both world space and material space.(b) After adding an auxiliary vertex for each triangle, a triangular mesh is transformed into a tetrahedral mesh.Only some of the auxiliary vertices are shown in (b).

Fig. 2
Fig.2Strain limiting for a triangular surface mesh.The top mesh is stretched lengthwise by a factor of 1/2 from the original undeformed mesh.The strain of the stretched mesh is limited to lie within [0.99, 1.01].The bottom-left mesh is produced by Narain et al.'s strain limiting method[3]; some triangles have flipped as highlighted in red.The bottom-right mesh is produced by our method; the strain is well limited, without flipping triangle.

Fig. 4
Fig.4 Comparison of strain limiting methods.The triangular mesh at top right (a) has 9600 triangles; it is stretched twice in length and compressed by half in height compared to the original un-deformed mesh.Mesh (b) was generated by our method.As the picture shows, our method generates the best result.The meshes in (c) were generated by Narain et al.'s strain limiting method [3] using 10, 33, and 987 iterations respectively.The table at top right shows the number of iterations taken by these two methods to converge.It is clear that our method converges faster (in one iteration) than Narain et al.'s method[3].The disadvantage of our method is that each iteration takes longer to perform.

Fig. 8
Fig. 8 Collision response benchmarks: (a) a ball hits three layers of cloth; (b) three layers of cloth fall on a rotating ball and are twisted by it; (c) a falling ball pushes three layers of cloth through a funnel.All of (a), (b), and (c) produce many simultaneous complex collisions which may lead to cloth simulation failure.(d) is a clothed dancing human, in which less complex collisions occur.

Table 1
[5]parison of collision response methods.Note the advantages of our collision response method over Harmon et al.'s inelastic projection impact zone approach[5].Using the same collision detection method, our method takes fewer iterations in each time step to handle simultaneous collisions.Furthermore, our method takes less time in each iteration to deal with multiple simultaneous collisions.Both contribute to a faster collision handling process (including collision detection and collision response) Harmon et al.'s collision response method[5]using inelastic projection.Inelastic projection is actually a velocity filter because it acts on velocities of meshes.It solves a projection optimization problem with linear equation constraints.