Codimensional Incremental Potential Contact
CCodimensional Incremental Potential Contact
MINCHEN LI,
University of Pennsylvania & Adobe Research
DANNY M. KAUFMAN,
Adobe Research
CHENFANFU JIANG,
University of Pennsylvania (a) initial (b) 0.5m/s (c) 1m/s (d) 2m/s(e) 4m/s ( f ) 8m/s Fig. 1.
Table cloth pull:
Codimensional Incremental Potential Contact (C-IPC) enables high-rate time stepping of codimensional models with intersection-free,accurate, controllable strain-limiting, thickness modeling, and frictional contact. Here, stepping shell and volumetric FE models at โ = . ๐ , a square cloth islaid on a table with heavy, stiff volumetric dinnerware (a). C-IPC stably responds to rapid pulling at varying speeds (b-f) without stretching beyond limitswhile accurately resolving tight contacts between table edges and dinnnerware. With a slower pull in (b) of . ๐ / ๐ nearly all objects pull off due to friction. Aspull speed increase at ๐ / ๐ (c) and ๐ / ๐ (d), less dinnerware fall off with some just pulled a bit. In (e) and (f) pulling at higher speeds ๐ / ๐ or ๐ / ๐ leavesthe stetting on table with rapid acceleration overcoming friction. Here too, combined strain-limiting and resolution of elastodynamics generates detailedwrinkling behaviors in the cloth for these faster pulls as it flies off the table. We extend the incremental potential contact (IPC) model [Li et al. 2020a]for contacting elastodynamics to resolve systems composed of codimen-sional degrees-of-freedoms in arbitrary combination. This enables a unified,interpenetration-free, robust, and stable simulation framework that couplescodimension-0,1,2, and 3 geometries seamlessly with frictional contact. Ex-tending the IPC model to thin structures poses new challenges in computingstrain, modeling thickness and determining collisions. To address these chal-lenges we propose three corresponding contributions. First, we introducea ๐ถ constitutive barrier model that directly enforces strain limiting as anenergy potential while preserving rest state. This provides energetically-consistent strain limiting models (both isotropic and anisotropic) for cloththat enable strict satisfaction of strain-limit inequalities with direct couplingto both elastodynamics and contact via minimization of the incrementalpotential. Second, to capture the geometric thickness of codimensional do-mains we extend the IPC model to directly enforce distance offsets. Ourtreatment imposes a strict guarantee that mid-surfaces (respectively mid-lines) of shells (respectively rods) will not move closer than applied thickness Authorsโ addresses: Minchen Li, University of Pennsylvania & Adobe Research,[email protected]; Danny M. Kaufman, Adobe Research, [email protected]; Chenfanfu Jiang, University of Pennsylvania, [email protected]. values, even as these thicknesses become characteristically small. This en-ables us to account for thickness in the contact behavior of codimensionalstructures and so robustly capture challenging contacting geometries; anumber of which, to our knowledge, have not been simulated before. Third,codimensional models, especially with modeled thickness, mandate strictaccuracy requirements that pose a severe challenge to all existing continuouscollision detection (CCD) methods. To address these limitations we developa new, efficient, simple-to-implement additive CCD (ACCD) method that iter-atively refines a lower bound converging to time of impact. In combinationthese contributions enable codimensional
IPC (C-IPC). We perform extensivebenchmark experiments to validate the efficacy of our method in capturingintricate behaviors of thin-structure contact and resulting bulk effects. Inour experiments C-IPC obtains feasible, convergent, and so artifact-freesolutions for all time steps, across all tested examples โ producing robustsimulations. We test C-IPC across extreme deformations, large time steps,and exceedingly close contact over all possible pairings of codimensionaldomains. Finally, with our strain-limit model, we confirm C-IPC guaranteesnon-intersection and strain-limit satisfaction for all reasonable (and wellbelow โ verified down to 0 . โข Computing methodologies โ Physical simula-tion . Additional Key Words and Phrases:
Strain Limiting, Contact Mechanics,Mixed-Dimensional Elastodynamics, Constrained Optimization a r X i v : . [ c s . G R ] J a n โข Li et al. Reference Format:
Minchen Li, Danny M. Kaufman, and Chenfanfu Jiang. 2021. CodimensionalIncremental Potential Contact. arXiv (February 2021), 26 pages.
Cloth, hair and sand are everywhere. Simulating thin materials likethese has long remained a critical task in computational model-ing and animation. Each such solid is generally best modeled by acorresponding reduced codimensional formulation; respectively byshells, rods and particles. Doing so enables improved efficiency withmany less simulated degrees-of-freedom (DOF). Likewise it miti-gates the numerical conditioning and severe convergence challengesposed by locking that we would otherwise encounter if simulatingthin materials volumetrically [Yang et al. 2000]. But while we canwell-capture the elastic behavior of individual thin materials withcodimensional models, accurately and consistently simulating theircollective and contacting behavior on practical size meshes, withreal-world materials and conditions remains challenged. We enu-merate these challenges here.First and foremost simulations should remain intersection free atall times. Steps that permit even slight intersections of a codimen-sional geometry can and will produce tangles that simulations cannot recover from.Second, we require energetically consistent, controllable and ac-curate strain-limiting. While codimensional discretizations mitigatelocking they do not remove it. For example, shell models continueto suffer severe numerical stiffening in bending modes for clothmaterials unless midsurface meshes apply extremely high (and gen-erally impractical size) resolutions [Quaglino 2012]. Strain-limitingallows softer cloth materials coupled with tight limits to avoid theseartifacts while still recovering proper cloth response. However, ifstrain-limiting is not accurately and consistently applied, and care-fully coupled to elastodynamics, it often produces more numericalartifacts than it fixes.Third, geometric thickness must be captured in contact. Whilethin structures typically have small relative thicknesses, whose elas-tic behavior can be indirectly captured by codimensional DOFs, cor-rectly and directly modeling their thickness in contact is crucial forcapturing accurate geometriesand bulk behaviors. Consider, forexample, a deck of cards. Eachcard has a nearly imperceptiblethickness when considered indi-vidually. However, when stacked, their combined thickness is largeand clear, and can not be ignored. While contact-processing strate-gies often introduce thickness parameters, they are generally appliedas a nonphysical failsafe to mitigate collision-processing inaccura-cies. As analyzed in Li et al. [2020a], these thicknesses can not beconsistently enforced and moreover must be changed heuristicallyper scene and example (e.g., based on collision speeds) to avoidsimulation failures.Fourth, all codimensions (including volumes) should be seam-lessly simulated in a common framework without distinction norspecial casing. Accurate frictional contact forces should directlycouple interactions between all geometric types irrespective of howclose the contact or how thin the modeled thicknesses are.
Fig. 2.
Twisted cylinder.
Testing thickness modeling, contact resolutionand buckling under extreme and increasing stress we start twisting a ๐ -wide cloth cylinder with 88K nodes, a thickness offset of . ๐๐ and ห ๐ of ๐๐ . Time-stepped at โ = . ๐ , C-IPC quickly obtains interesting foldsand soon after forms a thick central cylinder of wound cloth supported by C-IPCโs finite thickness offset (upper left and bottom). In contrast, (upper right)without C-IPCโs thickness offset the correct geometry can not form (norcan it later capture the materialโs buckling behavior) further emphasizingthe importance of consistent thickness modeling for shells. After . ๐ offurther twisting, our offset-thickened cloth continues to support complexcontact-driven behaviors including the final buckled geometry; see oursupplemental video. Fifth, contact modeling between codimensional models with small but finite thicknesses challenges all collision-detection routines.Most critically we see unacceptable failures and/or inefficiencies inall available existing CCD methods and codes. We see this both forstandard floating point root-finding CCD methods as well as morerecent developments in exact CCD methods. These issues are attrib-utable to the well-known numerical sensitivity of the challengingunderlying root-finding problems posed. Here, most specifically, wesee that IPC is reliant on a bounded guarantee of non-intersectionfor all CCD evaluations. Large enough conservative bounds canhelp [Li et al. 2020a] account for these inaccuracies, ensuring non-intersection, but they do so at the cost of progress and so can evenstall convergence altogether when it comes to codimensional DOF.Finally, simulations should not snag in contact on sharp andcodimensional geometries and so should be able to robustly solvetime-step problems to convergence irrespective of step size, scene,materials, contact and boundary conditions.Here we propose a method that, to our knowledge, is the firstto directly resolve all six above challenges with a consistent andreliable solution for frictionally contacting dynamics. Our methodguarantees strict satisfaction of fully coupled strain-limits, withintersection-free trajectories, and controllable geometric thicknessresolution (at thin material scales) independent of time step size andseverity of contact conditions.
Specifically, to address these challenges we extend the incrementalpotential contact model (IPC) [Li et al. 2020a] for contacting elasto-dynamics to resolve systems composed of arbitrary combinationsof codimensional DOF with three critical components. odimensional Incremental Potential Contact โข 3
Fig. 3.
Hairs.
Left: two hair clusters are twisted to form a tightly wound braid without intersection (processing up to 1.5M contacts per time step); the bottomend is then released and the braid unwinds. Right: a hair cluster simulation example inspired by McAdams [2009] where one end of a hair cluster is fixed, andso falls under gravity upon another hair cluster with two ends fixed and then across a sphere.
Constitutive strain limiting.
We introduce a ๐ถ constitutivebarrier model that directly enforces strain limiting as anenergy potential while preserving rest state. This providesenergetically-consistent strain limiting models (both isotropicand anisotropic) for cloth that enable, for the first time, strictsatisfaction of strain-limit inequalities for all iterations andso for all time steps (verified down to 0 . IPC thickness model.
To capture the geometric thickness ofcodimensional domains in contact we extend the IPC modelto directly enforce distance offsets. Our treatment imposesa strict guarantee that mid-surfaces (respectively mid-linesand points) of shells (respectively rods and particles) willnot move closer than applied thickness values, even as thesethicknesses become characteristically small; see Sec. 5. Thisenables us to account for thickness in the contacting behaviorof codimensional structures and so capture challenging con-tacting effects; a number of which, to our knowledge, havenot been simulated before.
Additive CCD.
To provide the strict accuracy required of CCDto resolve codimensional models with thickness, we develop anew, efficient and simple-to-implement additive CCD (ACCD)method. ACCD iteratively accumulates a lower bound con-verging to time of impact and is stable for the exceedinglychallenging evaluations required. While we most immedi-ately focus here on ACCDโs application within our C-IPCframework, as we show in the following sections, ACCD isexceedingly simple and so easy to implement when comparedto all prior CCD routines (exact and floating point) with bothimproved performance, robustness and guarantees, and so issuitable for easy replacement in all applications where CCDmodules are employed.Together these form the core of codimensional IPC (C-IPC). C-IPC enables unified simulation of all codimensions including elasticvolumetric bodies, shells, rods and particles all coupled togethervia accurately solved contact and friction. Along with these core contributions we test, compare and analyze C-IPC across many newand pre-existing benchmarks. We confirm C-IPC provides tightlycontrollable geometric thickness behaviors and guarantees feasibil-ity at all time steps. Finally, C-IPC remains intersection-free andstrictly satisfies all set strain limits across all computed trajectories,likewise verified in our extensive benchmark testing.
Beginning with the pioneering work of Terzopoulos et al. [1987]the simulation of codimensional models, especially shells and rods,has remained a focus in computer graphics. Efficiently modeling thecomplex behaviors of cloth [Baraff and Witkin 1998; Bender et al.2013; Bridson et al. 2002; Grinspun et al. 2003; Harmon et al. 2009; Liet al. 2020b; Narain et al. 2012; Volino and Thalmann 2000] and hair[Choe et al. 2005; Daviet et al. 2011; Deul et al. 2018; Kaufman et al.2014; Kugelstadt and Schรถmer 2016; McAdams et al. 2009; Mรผlleret al. 2012; Selle et al. 2008; Ward and Lin 2003] are now particularlycritical tasks across applications.Pipelines for their simulation most commonly adopt implicit orelse linearly implicit time-integration methods [Baraff and Witkin1998; Bridson et al. 2002; Harmon et al. 2008; Kim 2020; Li et al.2020b, 2018a; Narain et al. 2012; Otaduy et al. 2009; Tang et al. 2016,2018] with a diverse array of collision filters and penalties appliedto help resolve contact processing.To accelerate performance, extensions to the GPU [Li et al. 2020b;Schmitt et al. 2013; Tang et al. 2013, 2016, 2018], fast projection[English and Bridson 2008; Goldenthal et al. 2007] and multigrid[Tamstorf et al. 2015; Wang et al. 2018; Xian et al. 2019] are all beingactively explored. While, to improve the fidelity of models, data-driven material estimation [Bhat et al. 2003; Clyde et al. 2017; Miguelet al. 2012], alternate spatial discretizations [Guo et al. 2018; Jianget al. 2017; Weidner et al. 2018] and even hybrid models combiningyarn and shell simulation [Casafranca et al. 2020; Sperl et al. 2020]are being investigated. Likewise, the differentiability of cloth simu-lation is now also becoming critical [Liang et al. 2019] for neuraltraining applications. โข Li et al. (a) (b)(c) (d)
Fig. 4.
Reef Knot.
To exercise extreme compression, contact and frictionwith large stresses we extend the challenge of the reef knot example fromHarmon et al. [2009] by simulating with 100K nodes ( ร that of the original),applying a strain limit of . , friction with ๐ = . , and pulling eventighter. Two ribbons are initially intertwined in (a) and then stretched in(b)-(c) to form a final tight knot as the ribbons are pulled to nearly theirstrain limit in (d). A key challenge then has been to reliably simulate shells withguaranteed results. Harmon et al. [2009] introduce an explicit timestep-ping method for simulating shells that offers a guarantee of non-intersection across all timesteps. This method then requires smalltimesteps and simplified friction. Here we provide a complementary,fully implicit, differentiable simulation method, unified for all codi-mensional models, with accurate frictional contact, guaranteeingnon-intersection, as well non-inversion and strain-limit satisfaction(when desired), that can be stepped at all reasonable time step sizes.
Strain-limiting methods seek to impose bounds on membrane de-formation. A wide range of models and algorithms have long beenproposed to provide strain limiting [Bridson et al. 2002; Provot et al.1995]. To understand recent methodsโ behaviors and limitations wecategorize them by: 1) constraint choice (limits imposed as eitherequality or inequality); 2) type of DOFs constrained (these are mostcommonly edge-based or singular values); 3) splitting model; and 4)solver applied to enforce these constraints.With this breakdown we see that equality constraints on length-based measures remain a consistent constraint choice. Goldenthal etal. [2007] apply bilateral constraints on quadmesh edge lengths andpropose a fast projection method to correct predicted displacements.English and Bridson [2008] adopt a non-conforming strategy, apply-ing equality constraints on triangle-edge midpoint distances, andextend fast projection to support a BDF-2 integrator. Thomaszewskiet al. [2009] apply lagged, corotational small strain on triangles todefine equality constraints and enforce limits by post-projection viaGauss-Seidel or Jacobi iterations. Chen and Tang [2010] likewisedefine equality constraints on triangle edge lengths. Improvementover edge-based constraints can often be obtained by constrainingsingular values of the deformation gradient tensor (generally pertriangle) to a finite range.However, irrespective of details, equality constraints always re-main active and so simple DOF counting generally explains themembrane locking artifacts often encountered when these strain constraints are enforced. Here, alternate constraint DOF choices[Chen et al. 2019; English and Bridson 2008] can help alleviate thisissue by reducing the ratio of constraints to DOFs but also can in-troduce new challenges, e.g. via non-conforming meshes. As analternative to bilateral enforcement inequality constraints on strain,in the form of upper and lower bounds [Hernandez et al. 2013;Narain et al. 2012; Wang et al. 2010, 2016] have long been applied;while, even more recently, applying just upper bounds [Jin et al.2017] has also been considered. Irrespective of bounds, unilateralconstraints are only activated when their strain measures are attheir limits and so the potential for an over-constrained system isreduced.While approaches thus clearly vary we observe that all meth-ods, with the exception of Chen and Tangโs [2010] least-squaresapproximation for frictionless elastostatics, currently employ time-step splitting. Here we see that two-step splits are often appliedto first solve elastodynamics before applying a constraint projec-tion to simultaneously resolve contact constraints and strain limits.Alternately, introducing even more errors, we also see three-stepsplits where sequential solves of elastodynamics, contact and strainlimiting are each decoupled and so resolved independently per timestep. Thus, for the latter three-step approach strain limits and con-tact constraints cannot be simultaneously satisfied, while for theformer two-step methods errors introduced by splitting elasticityand constraints also generate unacceptable artifacts; see Section 6.Finally, we note that, to our knowledge, no method (irrespec-tive of constraint-type or splitting choice) employs a constraintsolver that can guarantee enforcement of strain limits. As constraintenforcement errors then vary across simulation scenes, and evenindividual time steps, this leads to uncontrollable and inconsistentvariations in material behavior per scene and step. To enable con-sistent, effective strain-limiting C-IPC introduces a fully coupled,inequality based model solved with a strict guarantee of strain-limitsatisfaction.
Thin bodies modeled by codimensional geometries generally havethicknesses orders-of-magnitude smaller than other dimensions.It is thus tempting (and common) to ignore thickness in contact.However, in order to correctly capture thin material interactionswe must account for the geometric effects of thickness in contact.We see this everywhere, for example when playing cards stack in apile (Figure 8) or when noodles fill a bowl (Figure 10).A direct strategy for capturing thickness is then to model thinmaterials volumetrically with full translational DOF [Hauptmannand Schweizerhof 1998]. However, in doing so linear elements sufferfrom well-known shear-locking artifacts. And while higher-orderelements can alleviate these convergence issues, locking artifactsare still not easily avoided altogether, while computational costs in-crease significantly. To address these challenges solid shell methodscommonly employ reduced integration with hourglass mode stabi-lization [Reese 2007] or assumed natural strain methods [Cardosoet al. 2008] to mitigate locking with linear elements. However, thisadded complexity generally does not compete with codimensionalmodeling odimensional Incremental Potential Contact โข 5 (a) (b) (d) (e) ( f ) (g)(c) Fig. 5.
Garments.
A yellow dress with challenging folded and seamed knife pleats produced by FoldSketch [Li et al. 2018b] is staged (a) and then draped on amannequin to static equilibrium (b) with stitches seamed together and then simulated on a Rhumba dancing sequence (c). A multilayer skirt is draped in thesame way (d), and then used to simulate an input character animation sequence with large and fast motions (e,f,g). As the mannequin rapidly kicks, intricategarment details and high speed collisions are captured by C-IPC simulating with large time steps ( โ = . ๐ ). To resolve contact with codimensional models contact-processingstrategies often introduce thickness-like parameters to offset con-straints and so reduce collision-processing inaccuracies; see e.g.Narain et al. [2012] and Li et al. [2018a] for recent examples. How-ever, as analyzed in Li et al. [2020a], these thicknesses can notbe consistently enforced and moreover must be changed heuris-tically per scene and example (e.g., based on collision speeds) toavoid simulation failures from intersections and instabilities. In turn,while sometimes helpful to improve contact resolution, these tun-ing parameters can not reliably be applied to model consistent andchanging thickness behaviors; see Section 6.5.We extend the IPC model to capture the geometric thicknessof codimensional domains with offset barriers that guarantee arequested small minimal separation from mid-surface as well mid-line and point geometries. In turn this enables reliable and consistentmodeling of thickness behaviors in the interactions between thincontacting materials.
Continuous Collision Detection (CCD) methods have long beenemployed as a check against intersecting mesh boundaries. Provot[1997] formulated finding first times of impacts for linearly displac-ing point-triangle and edge-edge pairs by first testing for coplanarityvia solution of a cubic equation and then performing overlap checksto detect collision.This by now standard formulation of floating point CCD hasbeen broadly applied and fine tuned over time [Harmon et al. 2009;Tang et al. 2011]. However, while root-finding with floating-pointarithmetic can be made efficient, significant numerical errors willstill generate unacceptable, false results (both negative and positive).This is most commonly encountered when distances are small and/orconfigurations are unavoidably degenerate. Here, while switchingfrom floating-point to rationals can help avoid round-off issues, costcan unacceptably increase if special care is not taken.Recent methods address CCD accuracy [Brochu et al. 2012; Tanget al. 2014; Wang et al. 2015] by application of exact arithmeticfor improved robustness with efficiency. Alternately, others [Har-mon et al. 2011; Lu et al. 2019] perform conservative CCD with arequested small conservative separation distance to better avoid interpenetration. This latter strategy is employed in volumetric IPC[Li et al. 2020a] using a robust, floating-point CCD implementation as a base solver.Here, for modeling thickness in C-IPC, we require CCD queriesthat can maintain finite separation distances between mid-surface,mid-line and point elements. In turn, we rely on accurately pre-serving orders-of-magnitude smaller distances between the offsetsurfaces for computing accurate contact forces. This challenges theaccuracy and robustness of CCD queries and we see unacceptableerrors, resulting in failures, in all available existing CCD methodsand codes. We see this both for standard floating-point root-findingCCD methods as well as more recent exact CCD methods. In turn wesee that these failures can generate intersections, slow convergence,and often stop simulation progress altogether; see Section 6.6.For robust CCD evaluation we derive a useful lower bound ontime of impact and apply it to develop a new, simple to compute,numerically robust, floating-point, additive CCD (ACCD) algorithm.ACCD monotonically approaches times of impact without error-prone, direct root-finding. In Section 6 we confirm ACCD efficientlyand accurately succeeds on a wide range of challenging cases whereall other methods fail and find similar and often better performancein the cases where floating-point CCD methods can succeed. Finallywe verify that ACCD is also suitable for replacement in CCD mod-ules outside of the C-IPC framework with improved efficiency androbustness. With the diverse range of geometries commonly found in envi-ronments the unified simulation of all codimensions in one frame-work is critical for simulation efficiency and accuracy. Martin et al.[2010] focus on a unified elasticity model. They derive elastons โ ahigher-order integration rule to measure stretch, shear, bending andtwisting along all axes without distinction between codimensions.Elastons accurately capture a wide range of elastoplastic behaviors,while contact forces are determined by point-wise penalties, withobjects represented by a set of spheres. Chang et al. [2019] addressunification for mixed-dimensional elastic bodies by defining all https://github.com/evouga/collisiondetection โข Li et al. connections between domains of varying codimension via equalityconstraints. Here, Bridson et al.โs [2002] collision processing algo-rithm is then applied to resolve contact. The material point method(MPM) [Jiang et al. 2017] also offers general-purpose modeling ofall codomains and contact between them. MPM discretizes elasticityon Lagrangian particles while solving momentum balance on theDOFs of an Eulerian grid. Contact between codimensional objects isthen directly resolved via the Eulerian grid as a velocity flow. How-ever, sticking and gap errors in contact are well known artifactsin MPM and can be unacceptable if grid resolutions are too low.Position-based dynamics (PBD) [Bender et al. 2014, 2015; Mรผlleret al. 2007] also enables seamless, unified coupling of bodies withvarying codimensions. Here both constitutive model and contactsare resolved as constraints that are iteratively processed for efficienttime-integration at the cost of controllable accuracy as constraintcomplexity increases [Li et al. 2020a]. Extensions of PBD with XPBD[Macklin et al. 2016], Projective Dynamics [Bouaziz et al. 2014], andits further generalization via ADMM [Overby et al. 2017] all sim-ilarly offer platforms for co-simulating a diversity of codomains.Recent enhancements now provide improved friction [Ly et al. 2020]and increased efficiency [Daviet 2020]. However, the fundamentaltrade-off of accuracy and robustness (in the resolution of elasticityand contact) for efficiency in computation remains[Li et al. 2019; Lyet al. 2020]. C-IPC enables the direct simulation of all codimensionssimultaneously. Coupling is provided by the interaction betweenany and all codimensional pairings via accurate, intersection-freecontact. We focus on the combined solution of meshed, codimensional mod-els simulated jointly and so coupled arbitrarily via contact. To doso we generalize the IPC model to codimensional DOF and so mustaddress the challenges introduced by thin models both coupled by,and stressed via contact and large, imposed boundary conditions.
Fig. 6.
Cloth on rotating sphere.
A square cloth (85K-node 1st row, 246K-node 2nd row) is dropped upon a sphere set on the ground with friction( ๐ = . for both). As the sphere starts to rotate the cloth follows andwinds tightly while preserving strain limits throughout. Differing moduliare used for the ๐พ - and ๐พ -node cloth respectively to produce intricatewrinkling behaviors with differing frequency. Here from left to right weshow the 25th (right before rotation), 50th, and 75th frames. Here we first cover our generalization of IPC to mixed codimen-sional models and then, in the following sections, construct the keycomponents of our method that enable their simulation.
Elastodynamics with contact.
For elastodynamics, we performimplicit time stepping with optimization time integration [Gastet al. 2015; Kaufman and Pai 2012; Li et al. 2019; Liu et al. 2017;Overby et al. 2017; Wang et al. 2020] to minimize an IncrementalPotential (IP) [Kane et al. 2000] with line search, ensuring stabilityand global convergence. For a wide range of time steppers andassuming hyperelasticity, the IP to update from time step ๐ to ๐ + ๐ธ ( ๐ฅ ) = || ๐ฅ โ ห ๐ฅ ๐ || ๐ + ๐ผโ ฮจ ( ๐ฝ๐ฅ + ๐พ๐ฅ ๐ ) , (1)with the timestep update given by ๐ฅ ๐ + = argmin ๐ฅ ๐ธ ( ๐ฅ ) . Here โ is time step size, ๐ is the mass matrix, and ฮจ is an elastic energypotential. Scaling factors ๐ผ, ๐ฝ,๐พ โ [ , ] and explicit predictor ห ๐ฅ ๐ then determine the time step method applied. For example, herewe focus primarily on the graphics-standard implicit Euler with ๐ผ, ๐ฝ = ๐พ = ๐ฅ ๐ = ๐ฅ ๐ + โ๐ฃ ๐ + โ ๐ โ ๐ ext . Alternate implicittime stepping, e.g., with implicit midpoint or Newmark, are similarlymatched by applying alternate scalings and predictors [Kaufmanand Pai 2012].Incremental Potential Contact (IPC) then augments the IP withcontact and dissipative potentials [Li et al. 2020a]: ๐ธ ( ๐ฅ ) = || ๐ฅ โ ห ๐ฅ ๐ || ๐ + ๐ผโ ฮจ ( ๐ฝ๐ฅ + ๐พ๐ฅ ๐ ) + ๐ต ( ๐ฅ ) + ๐ท ( ๐ฅ ) . (2)Here ๐ต ( ๐ฅ ) and ๐ท ( ๐ฅ ) are respectively the IPC contact barrier andfriction potentials. The contact barrier, ๐ต , enforces strictly positivedistances between all primitive pairs while ๐ท provides correspond-ing friction forces. Following Li et al. [2020a] we apply a customNewton-type solver to minimize each IP with Continuous Colli-sion Detection (CCD) filtering executed in each line search of everyNewton iteration to ensure intersection-free trajectories throughout.Please see Li et al. [2020a] for details on the base IPC algorithm andsolver implementation. Mixed-dimensional hyperelasticity.
To enable a unified simulationframework for coupled volumetric bodies, shells, rods and particles,we compute mass and volume for all codimensional elements bytreating them as continuum regions with respect to standard dis-cretizations and so construct the total elasticity energy ฮจ ( ๐ฅ ) for theIP as ฮจ ( ๐ฅ ) = ฮจ vol ( ๐ฅ ) + ฮจ shell ( ๐ฅ ) + ฮจ rod ( ๐ฅ ) . (3)Here, without loss of generality, as representative examples, weapply fixed Corotated elasticity [Stomakhin et al. 2012] for volumes( ฮจ vol ); the Discrete Shells hinge bending energy [Grinspun et al.2003; Tamstorf and Grinspun 2013] combined with either isotropicor orthotropic StVK [Chen et al. 2018; Clyde et al. 2017] or neo-Hookean membrane models for shells ( ฮจ shell ); and the DiscreteRods stretch and bending model [Bergou et al. 2008] for rods ( ฮจ rod ).We select these models as best for comparison and evaluation giventhe models standard in existing, available codes. While, more gener-ally, the C-IPC framework is agnostic to a broad range of elasticityand time stepper choices as demonstrated in Section 6 and [Li et al. odimensional Incremental Potential Contact โข 7 Here we begin by constructing a new constitutive barrier modelthat directly enforces strain limiting while maintaining rest-stateconsistency. We start with a formulation that augments existingmembrane energies with an added strain-limiting potential for thegeneral, isotropic case (Section 4.1). We then demonstrate its exten-sion to directly augment orthotropic STVK membranes with strainlimits in Section 4.2.Because of its potential-based model, our isotropic (respectivelyanisotropic) constitutive strain limiting is then easily and directlyapplied in C-IPC, as just a simple addition to (respectively modifica-tion of) the elasticity potential, ฮจ , in (2). This provides a mesh-basedsimulation framework that produces trajectories simultaneously en-suring intersection-free and strain-limit satisfying steps. In turn,the resulting dynamics satisfy accurate momentum balance, fullycoupling frictional contact forces, strain-limiting and elasticity. We define isotropic strain-limiting constraints per element ๐ก as ๐ ๐ก๐ < ๐ , โ ๐, (4)where the singular value decomposition of each triangle ๐ก โs deforma-tion gradient, ๐น ๐ก = ๐ ๐ก ฮฃ ๐ก ๐ ๐ก๐ , gives ฮฃ ๐ก = diag ( ๐ ๐ก , ๐ ๐ก ) . The imposedconstraint bound ๐ is then the requested strain limit with practicalbounds generally selected for cloth with ๐ โ [ . , . ] [Bridsonet al. 2002; Provot et al. 1995].Inspired by IPCโs contact barrier formulation, we then constructa local energy that enforces this unilateral strain-limit. We couldbegin with a typical log-barrier strategy. This would give energies ๐ ๐ โ๏ธ ๐ โ ln ( ๐ โ ๐ ๐ก๐ ) , (5)with ๐ ๐ defining the barrier stiffness for limit ๐ . However, whileensuring strain-limits this barrier is problematic. First, from themodeling perspective, it applies forcing everywhere. Specifically,constitutive behavior is unnecessarily changed far from strain lim-its where the underlying membrane model should be untouchedand, worse yet, applies ghost forces at rest. Second, from the com-putational perspective, because this barrier is globally non-zero, Fig. 7.
Funnel.
Three cloth panels ( ๐ = . , 26K nodes each, strain upperbound at . ) are dropped on a funnel. Under friction they rest stably ontop until a sharp, scripted collision object rapidly pulls the panels down. evaluating it (and its derivatives) in optimization is prohibitivelyexpensive.From these issues we then observe that strain-limit forcing fortriangles far from the strain limit is then both unacceptable andunnecessary. Instead, we want local support for nonzero forces justclose to where these limits are achieved. One simple and seeminglyplausible remedy, to achieve this localized support region, is toclamp the barrier in (5) to zero beyond a small strain threshold ห ๐ (e.g. ห ๐ =
1) giving barrier energies ๐ ๐ โ๏ธ ๐ ๐ ๐ ๐ก๐ ( ๐ฅ ) , ๐ ๐ ๐ก๐ ( ๐ฅ ) = (cid:40) โ ln ( ๐ โ ๐ ๐ก๐ ๐ โ ห ๐ ) ๐ ๐ก๐ > ห ๐ ๐ ๐ก๐ โค ห ๐ . (6)Unfortunately, while providing local support these energies areonly ๐ถ continuous and so do not allow for practical gradient-basedoptimization necessary for the efficient solution of the IP.To smooth our strain-limiting while providing local support wethen propose a ๐ถ smoothly clamped barrier for strain-limiting ๐ ๐ ๐ก๐ ( ๐ฅ ) = (cid:40) โ( ห ๐ โ ๐ ๐ก๐ ๐ โ ห ๐ ) ln ( ๐ โ ๐ ๐ก๐ ๐ โ ห ๐ ) ๐ ๐ก๐ > ห ๐ ๐ ๐ก๐ โค ห ๐ . (7)Now, with barriers in hand, we could potentially next considertreating these barriers directly as constraints (as in primal barrierand interior point methods). However, doing so would then simplysum these barriers over elements and so would obtain inconsistentbehavior as we change meshes. Instead, to provide consistent be-havior we impose strain limiting constitutively as an energy densityintegrated over the volume of cloth to obtain the potential ฮจ SL ( ๐ฅ ) = โ๏ธ ๐ ๐ ๐ โซ ฮฉ ๐ ๐ ๐ก๐ ( ๐ฅ ) ๐๐ โ ๐ ๐ โ๏ธ ๐ก,๐ ๐ ๐ก ๐ ๐ ๐ก๐ ( ๐ฅ ) . (8)For our approximation, we use ๐ ๐ก = ๐ด ๐ก ๐ ๐ก , with ๐ด ๐ก and ๐ ๐ก respec-tively the area and thickness of triangle ๐ก . Our final isotropic strain-limiting potential, ฮจ SL , is then ๐ถ with local support and can simplybe added to our total potential in Equation 3. In turn, strain limitingis then directly handled at each time step by optimization of the IPCEquation 2, while ensuring that strain-limit barriers are not violatedin each line-search step. During steps when most triangles remainunder the strain limit threshold, little additional computation is thenrequired. While, when stretch increases, we can adjust the necessarynonzero terms from newly activated barriers and so, as we show inSection 6, strictly enforce strain limits while balancing all appliedforces. โข Li et al. Strain-Limit Stiffness.
With our barrier potential defined we nowspecify setting its stiffness. The mollified clamping strategy we applyfor our strain-limit model is inspired by IPCโs smoothing of contactbarrier energies [Li et al. 2020a]. Here notice an extra ingredient forstrain-limiting is the applied 1 /( ๐ โ ห ๐ ) factor. This scaling enables usto apply our barrier directly to a limit-normalized strain, measuredby ๐ฆ = ( ห ๐ โ ๐ ๐ก๐ )/( ๐ โ ห ๐ ) with ๐ ๐ ๐ก๐ ( ๐ฆ ) = (cid:40) โ ๐ฆ ln ( + ๐ฆ ) ๐ฆ < , ๐ฆ โฅ . (9)In turn, as our applied strain-limit ( ๐ ) and/or clamping threshold(ห ๐ ) is varied with application, the barrier function w.r.t. ๐ฆ remainsunchanged. Only the gradient of the linear map from ๐ ๐ก๐ to ๐ฆ varies.This allows us to apply a single, consistent initial barrier stiffness ๐ ๐ (we use 1 ๐พ๐๐ for all examples) across all choices of differing strainlimits and clamping thresholds. Here the barrier potential curve isthen simply linearly rescaled each time, to a different strain range,and so provides consistent conditioning to the system. Finally, toavoid numerical issues introduced by tiny gaps between a currentstrain and the imposed strain limit (e.g., when extreme boundaryconditions are imposed as in Figures 9 and 1) we adapt barrierstiffness when needed. Starting with our initial ๐ ๐ , we increase itby 2 ร (up to max bound of ๐ ๐ = . ๐๐๐ ) whenever the strain gap ๐ โ ๐ ๐ก๐ of a triangles ๐ก is less than 10 โ ( ๐ โ ห ๐ ) over two consecutiveiterations. Above we have constructed strain-limiting as an isotropic constitu-tive model. We formed a barrier energy that can be integrated andso directly added to potential energy in order to augment existingmembrane models with hard strain limits. The key to making thiswork is the application of our ๐ถ continuous clamping so that theapplication of strain-limits does not alter rest-shape gradients norrest-shape Hessians.Alternately, we can apply a comparable constitutive strain-limitingstrategy to directly modify membrane elasticity models to includestrain limits. To do so we simply substitute an anisotropic mem-brane model with a barrier energy that prevents violation of strainlimits, while matching the original membrane energy gradient andHessian at rest.This second strategy is particularly motivated by the need toincorporate strain-limiting into available data-driven models thathave been constructed to carefully fit measured cloth data. Herewe demonstrate specific application to the anisotropic, data-drivenmodel of Clyde and colleagues[Clyde 2017; Clyde et al. 2017].Their constitutive model energy isห ๐ ( ห ๐ธ , ห ๐ธ , ห ๐ธ ) = ๐ ๐ ( ห ๐ธ ) + ๐ ๐ ( ห ๐ธ )+ ๐ ๐ ( ห ๐ธ ห ๐ธ ) + ๐บ ๐ ( ห ๐ธ ) , (10)where ห ๐ธ = ๐ท ๐ ๐ธ๐ท is the reduced and aligned Green-Lagrangianstrain, ห ๐ธ ๐ = ห ๐ธ ๐ =
0, and column vectors of ๐ท are the tangent andnormal bases of the shell in material space. Functions ๐ are then asum of polynomial functions with real-valued exponents ๐ผ ๐๐ that satisfy ๐ ( ) = ๐ โฒ ( ) =
1. Here the first constraint enforces zero-energy, zero-stress rest configurations, while the second constraintallows a natural correspondence between the parameters ๐ ๐ผ๐ฝ and ๐บ and linear elasticity at infinitesimal strain. Clyde et al. use ๐ ๐ ( ๐ฅ ) = ๐ ๐ โ๏ธ ๐ = ๐ ๐๐ ๐ผ ๐๐ (( ๐ฅ + ) ๐ผ ๐๐ โ ) . Their measured data is then restricted to deformations within thefracture limit or elasticity range of the cloth. Beyond this rangemeaningful extrapolation is then unlikely, due to possible overfit-ting inside the range (exponents ๐ผ ๐๐ from the fitting can range upto 10 ). To address this limitation Clyde et al. propose a quadraticextrapolation for simulation. However such extrapolation is notphysically meaningful. To usefully apply data-driven modeling ei-ther fracture should be captured beyond this regime or else a stableand controllable strain limit imposed to stay in bounds. Here, wefocus on controllably imposing strain limit while respecting theunderlying model fitting. Barrier Formulation.
Starting from the same constitutive modelwe redefine the ๐ basis with barriers ๐ ( ๐ธ ) = โ ๐ธ ๐๐๐ฅ log (( ๐ธ ๐๐๐ฅ โ ๐ธ )/ ๐ธ ๐๐๐ฅ ) . For consistency note this basis again satisfies both ๐ ( ) = ๐ โฒ ( ) = ๐ธ ๐๐๐ฅ and so ensures that ๐ธ never exceeds measured strain limitbounds. It thus captures data-driven strain limits and at the sametime avoids extrapolation. Here, as the underlying model is anisotropic,warp, weft, and shearing directions can all impose different, mea-sured ๐ธ ๐๐๐ฅ bounds, enabling preservation of measured anisotropyfor all strain limits.Free model parameters ๐ ๐ผ๐ฝ and ๐บ , are then matched at โ ๐ ( , , ) with the underlying modelโs fit; see Appendix for details. Note thisformulation frees us from computing the sensitive (and potentiallyexpensive) polynomials with large real-valued exponents. See Sec-tion 6.2 for experiments on the proposed model. Finally we notethat, in contrast to the strategy we demonstrated in our isotropiccase from the last section (where we worked directly on an upperbound w.r.t. the deformation gradient), here this second energy issymmetric for stretch and compression. As discussed above, CCD-based filtering is critical to ensure line-search is consistent with our contact barriers. Now that we addstrain-limiting barriers we must additionally ensure that every linesearch will also not violate imposed strain limits. For isotropic limits,although 2 ร odimensional Incremental Potential Contact โข 9 For our anisotropic model we also currently apply the same bisectionstrategy. However, while currently not needed for efficiency, we donote that here we can formulate quadratic equations amenable todirectly computing largest feasible strain-limit satisfying step sizes.
In its original form IPC maintains intersection-free paths for volu-metric models by enforcing positivity of unsigned distances ๐ ๐ , as aninvariant between all non-adjacent and non-incident surface primi-tive pairs ๐ [Li et al. 2020a]. This is suitable for volumetric contactwhere this constraint permits arbitrarily close, but never intersect-ing surfaces. For codimensional models, however, this constraint isno longer sufficient. When the 3D deformation of thin materials isreduced to the deformation of a 2D surface or 1D curve, elasticitycan be well-resolved on surfaces and curves, but contact can not.Neglecting to account for finite thickness in codimensional contactgenerates unacceptable artifacts (see e.g. Figure 2) and clearly failsto capture geometries formed by thin-structure interactions (see e.g.Figures 18 and 8).We begin by observing that applying a larger ห ๐ , to codimensionalgeometries in IPC, models a thin responsive layer that resists com-pression in the normal direction and so forms an elastic thickness.In concert with this elastic layer we also generally require a corethickness beyond which further compression is not possible. Thisis needed to guarantee a minimum finite (and so visible) thicknesseven when deformation is extreme (see e.g. Figure 2c).To model thickness in contact for codimensional objects we buildan inelastic thickness model that combines the purely elastic layerwith a hard, and so inelastic offset to the contact barrier. For anonzero offset, ๐ , geometries are then guaranteed to be separatedfrom each other by ๐ . Elastic contact forces exert when distancesare below ๐ + ห ๐ and diverge when distances are at ๐ . Here largerห ๐ generates greater elastic response, while nonzero ๐ guarantees aminimum thickness. This holds even under extreme compression(see e.g. Figure 2e).We equip each surface element ๐ with a finite thickness ๐ ๐ . Thedistance constraint for primitive pairs ๐ on the surface, formedbetween element primitives ๐ and ๐ , are then ๐ ๐ ( ๐ฅ ) > ๐ ๐ = ๐ ๐ + ๐ ๐ . Boundaries of codimensional materials are thenapproximated with a rounded cross-section, whilefor interaction between zero-thickness materialsour distance constraints reduce to that of the origi-nal IPC constraint (e.g. for volume-to-volume con-tact). Finally for volume-to-codimensional con-tact, volumes can continue to maintain zero-thickness boundarieswhile interacting with finite thickness codimensional boundaries.For scenarios where modeling thickness matters (see e.g. Figures2, 3-23), we then set ๐ ๐ to the true thicknesses of the material(s)or slightly smaller, and then set ห ๐ near ๐ ๐ to compensate for theremaining thickness depending on how much compression, and soelastic response, is reasonable or required per application. For a numerically robust and efficient implementation, IPC appliessquared distances (with appropriate rescalings) to compute withequivalent barrier functions. This change swaps the modelโs derivedcontact barrier ๐ ( ๐ ๐ ( ๐ฅ ) , ห ๐ ) to an equivalent rescaled implmentationusing ๐ ( ๐ ๐ ( ๐ฅ ) , ห ๐ ) [Li et al. 2020a]. Here we directly apply ourthickened contact barrier via squared distances ๐ ๐ ( ๐ ๐ ( ๐ฅ ) , ห ๐ ) = ๐ ( ๐ ๐ ( ๐ฅ ) โ ๐ ๐ , ๐ ๐ ห ๐ + ห ๐ ) , so that contact forces correctly diverge at ๐ ๐ ( ๐ฅ ) = ๐ ๐ , and nonzerocontact forces are only applied between mid-surface pairs closerthan ห ๐ + ๐ ๐ .Here then the barrier function, gradient, and Hessian only needto be evaluated with the modified input distance, while distancegradient and Hessian remain unchanged as ๐ ( ๐ ๐ ( ๐ฅ ) โ ๐ ๐ ) ๐๐ฅ = ๐๐ ๐ ( ๐ฅ ) ๐๐ฅ . Evaluation of contact force magnitudes for computing friction forcesas well as adaptive stiffness ( ๐ ) updates for contact barriers followsimilarly.Most constraint set computations then require simple and com-parable modifications. For spatial hashing, bounding boxes of allmid-surface primitive ๐ are extended by ๐ ๐ / ๐ to remain unchanged. Next, for broad-phasecontact pair detection, query distances ห ๐ are now updated to ห ๐ + ๐ ๐ to check for bounding box overlap. Or else, equivalently, boundingbox primitives are extended as in the spatial hash. Finally, to ac-celerate continuous collision detection (CCD) queries, spatial hashconstruction also requires similar bounding box extension, whilefor its broad phase, applied gaps are ๐ ๐ (rather than 0). While the above initial modifications for thickness are straight-forward, finite thickness and codimensional DOF introduce new (b)(c) (d)(a)
Fig. 8. โPrecisionโ card shuffle.
We divide fifty-four playing cards intotwo separate piles and and bend them in preparation for a โprecisionโ bridgeshuffle. Unlike a human-performed bridge finish, where cards begin inter-leaved, here we increase challenge by keeping the two bridged piles apartand then (a) precisely, one-by-one, bottom up they are rapidly released inturn to shuffle the deck (b). We partially square the deck in (c), obtaining afully shuffled stack with non-intersecting, interleaved cards that is 15.4 mmheigh and so well-matching the height of a standard deck (d). (a) (b)(c) (d)
Fig. 9.
Cloth over codimensional needles.
To test contact across sharpboundaries we drop and then drag a cloth panel across a needle bed formedby line segments. With a large time step size of โ = . ๐ and strain limitof . , the cloth comes safely to a static rest on the segments withoutjittering (a). We then pull the left side of the cloth at ๐ / ๐ across the segmenttips (b) and observe rapid sliding along (c), and then off (d) the top of thebed without snagging nor stretching artifacts (strain limits are satisfiedthroughout). computational challenges for CCD queries. Here we analyze thesechallenges and develop a new CCD method to address them.Standard-form IPC applies position updates in each iteration toobtain minimal separation distances of ๐ ร the current separationdistances, ๐ ๐๐ข๐๐ , at time of CCD evaluations. Here ๐ โ ( , ) is aconservative rescaling factor (generally set to ๐ = . .
1) thatallows CCD queries to avoid intersection, even when surfaces arevery close [Li et al. 2020a]. To similarly resolve finite thicknesses,our conservative distance bound is now ๐ ( ๐ ๐๐ข๐๐ โ ๐ ๐ ) + ๐ ๐ .Concretely, the barrier evaluated parameters, ๐ ๐ ( ๐ฅ ) โ ๐ ๐ , mustalways remain positive. In turn we require all CCD evaluations, foreach displacement ๐ , to provide us with sufficiently accurate times ๐ก โ ( , ] to satisfy positivity. If impact would occur for a pair along ๐ , we require a time ๐ก so that the scaled displacement, ๐ก๐ , ensuresdistance remains greater than ๐ ๐ and is as close as possible to thetarget separation distance ๐ ๐ + ๐ ( ๐ ๐๐ข๐๐ โ ๐ ๐ ) .Doing this with finite thickness is then much more challengingthan without thickness. In contact ๐ ๐๐ข๐ โ ๐ ๐ can be as small as10 โ ๐ (e.g. under large compression) so that for ๐ โฅ . โ ๐ are acceptable. While, in practice, ๐ ๐ isthen regularly at the scale of 10 โ ๐ (e.g. thickness values of โผ ร โ ๐ for cloth). Without thickness ( ๐ ๐ =
0) updated distances aretargeting ๐ ๐ ๐๐ข๐๐ (and are only required to be strictly greater than 0).In other words any step size that prevents primitive-pair intersectionis valid and so relative errors less than 100% โ โ ๐ /( ๐ โ ๐ ) areacceptable. On the other hand, with thickness, updated distancesare targeting ๐ ๐ + ๐ ( ๐ ๐๐ข๐๐ โ ๐ ๐ ) and so, for standard values of ๐ ๐ ,would require relative errors from the CCD evaluations near 0 . โ โ ๐ /( โ ๐ + ๐ ( โ ๐ )) in order to avoid interpenetration.Obtaining CCD evaluations to this accuracy is then extremelychallenging for available methods. As a starting example we testedthe floating point CCD solver from Li et al. [2020a], requesting a ๐ ( ๐ ๐๐ข๐ โ ๐ ๐ ) + ๐ ๐ minimal separation for two challenging codimen-sional examples with thickness (see Figures 10 and 21) with ๐ ๐ > ๐ก = ๐ = Next we derive a useful lower bound value for CCD queries thatis numerically robust and can be efficiently evaluated in floatingpoint. This lower bound provides a conservative, guaranteed safestep size and also a clear measure to test a CCD queryโs validity:any CCD-type evaluation with smaller TOI has clearly failed. InSection 5.4 we then apply this bound to derive a simple, effectiveand numerically accurate explicit CCD solver that is robust andefficient for progress even in the challenging CCD evaluations werequire for thin material simulations.Here, without loss of generality, we will focus on the edge-edgecase between edge pairs ( ๐ฅ , ๐ฅ ) and ( ๐ฅ , ๐ฅ ) with correspondingdisplacements ๐ , ๐ , ๐ and ๐ . The distance function between ar-bitrary points parameterized respectively by ๐พ and ๐ฝ on each edgeat any time ๐ก is then ๐ ( ๐ก, ๐พ, ๐ฝ ) = || ๐ ( ๐ก, ๐พ, ๐ฝ )|| , with ๐ ( ๐ก, ๐พ, ๐ฝ ) = ( โ ๐พ ) ๐ฅ ( ๐ก ) + ๐พ๐ฅ ( ๐ก ) โ (( โ ๐ฝ ) ๐ฅ ( ๐ก ) + ๐ฝ๐ฅ ( ๐ก )) ,๐ฅ ๐ ( ๐ก ) = ๐ฅ ๐ ( ) + ๐ก๐ ๐ , and ๐ก, ๐พ, ๐ฝ โ [ , ] . (11)A CCD evaluation then seeks the smallest positive real ๐ก satisfying ๐ ( ๐ก ) = min ๐พ,๐ฝ ๐ ( ๐ก, ๐พ, ๐ฝ ) = . If such ๐ก exists we call it ๐ก ๐ . Parameters ( ๐พ ๐ , ๐ฝ ๐ ) = argmin ๐พ,๐ฝ ๐ ( ๐ก ๐ , ๐พ, ๐ฝ ) in turn give the corresponding points colliding at time ๐ก ๐ on eachboth edges.We can then express ๐ก ๐ as ๐ก ๐ = ๐ ( , ๐พ ๐ , ๐ฝ ๐ )||( โ ๐พ ๐ ) ๐ + ๐พ ๐ ๐ โ (( โ ๐ฝ ๐ ) ๐ + ๐ฝ ๐ ๐ )|| . Challenges to CCD evaluations then occur in determining ๐พ ๐ and ๐ฝ ๐ , when degeneracies and numerical error make floating-pointmethods both prone to false positives and false negatives.We can however, directly and accurately perform a distance queryw.r.t. the primitives at start (t=0) of evaluation, to find the distance, ๐ ( ) , that certainly satisfies ๐ ( ) โค ๐ ( , ๐พ ๐ , ๐ฝ ๐ ) . Then triangleinequality gives ||( โ ๐พ ๐ ) ๐ + ๐พ ๐ ๐ โ (( โ ๐ฝ ๐ ) ๐ + ๐ฝ ๐ ๐ )|| โค ||( โ ๐พ ๐ ) ๐ + ๐พ ๐ ๐ || + ||( โ ๐ฝ ๐ ) ๐ + ๐ฝ ๐ ๐ || , and since ๐พ ๐ , ๐ฝ ๐ โ [ , ] weget ||( โ ๐พ ๐ ) ๐ + ๐พ ๐ ๐ || โค max (|| ๐ || , || ๐ ||) and ||( โ ๐ฝ ๐ ) ๐ + ๐ฝ ๐ ๐ || โค max (|| ๐ || , || ๐ ||) . Put together we then have the practicaland directly computable bound on TOI ๐ก ๐ โฅ ๐ ( ) max (|| ๐ || , || ๐ ||) + max (|| ๐ || , || ๐ ||) . (12)More generally, for any query between any two primitives types, ourbound on ๐ก ๐ is correspondingly given by the appropriate distancefunction in the numerator (e.g., ๐ ( ) = ๐๐๐ ๐พ ๐ ,๐ฝ ๐ ,๐ผ ๐ ๐ ( , ๐พ ๐ , ๐ฝ ๐ , ๐ผ ๐ ) for point-triangle and ๐ ( ) = ๐๐๐ ๐พ ๐ ๐ ( , ๐พ ๐ ) for point-edge) and thesum of the max displacement from each of the two primitives inthe denominator. We then note that even when there is no smallestpositive time satisfying ๐ ( ๐ก ) = odimensional Incremental Potential Contact โข 11 (a) initial (b) static frame (d) static frame (polyline view)(c) static frame (lower view) Fig. 10.
Noodles. (a) We drop 625 ๐๐ -long noodles, modeled by discrete rods with C-IPC offset, into a bowl. (b) We set offset and ห ๐ to model noodlethickness enabling the bowl to fill when they reach equilibrium. (c) We then remove the bowl to show the final intertwined, intersection-free configuration,noting that thickness offsets are satisfied throughout. In (d) we then remove the volume render to highlight the geometry is simulated solely with discreterods polyline geometry. We next observe that, perhaps surprisingly, state-of-the-art floating-point CCD solvers can and will return TOI results smaller than ourlower bound, and so are clearly in error (See Figure 19).Finally, to improve our bound, we observe that it holds indepen-dently of the choice of the reference frame. Thus we can furthertighten the bound on each CCD query independently by pickingframes that reduce the norm of displacement vectors ๐ ๐ . For exampleby subtracting each ๐ ๐ by the average (cid:205) ๐ ๐ ๐ . We now apply our lower bound to build a new CCD algorithmthat iteratively updates and adds our lower bound over successiveconservative steps towards TOI. The resulting additive
CCD (ACCD)method then robustly solves for a bounded TOI, monotonicallyapproaching each CCD solution, while only requiring explicit callsto evaluate distances between updated primitive positions.At the start of each CCD query, to initialize the ACCD algo-rithm (see Algorithm 1), we first center the collision stencilโs dis-placement at origin to reduce our boundโs denominator, e.g., ๐ ๐ = max (|| ๐ || , || ๐ ||) + max (|| ๐ || , || ๐ ||) for edge-edge pairs, and so in-crease the step size we can safely take. If there is no relative motion( ๐ ๐ = ๐ = ๐ ( โ๏ธ ๐ ๐ ๐๐ โ ๐ ) to the offset surface based on the current squareddistance ๐ ๐ ๐๐ and the scaling factor ๐ . For this we use a formula thatis more robust to cancellation error; see lines 8-9 in Algorithm 1.Starting with a most conservative time of impact ๐ก = local scratch-pad copy of nodal positions, ๐ฅ ๐ , and initializethe current lower bound step, ๐ก ๐ , with Equation 12 (line 11).We then enter our iterative refinement loop (lines 12-21) to mono-tonically improve our TOI estimate ๐ก . At each iteration, we updateour local copy of nodal positions ๐ฅ ๐ with the current step ๐ก ๐ (lines13-14). If this new position achieves our target distance to the offsetsurface (becomes smaller than ๐ ) we have converged and the pre-vious ๐ก is the time of impact that brings distance just up to ๐ (line17). If not, we update our TOI estimate by adding the current ๐ก ๐ to ๐ก (line 18). Note that we always add our first lower bound step to ๐ก (line 16) as it is guaranteed to not bring distance closer than ๐ . Ifour TOI is now larger than ๐ก ๐ , the current minimum first time of Algorithm 1
Additive CCD procedure AdditiveCCD( ๐ฅ ๐ , ๐ ๐ , ๐ , ๐ , ๐ก ๐ , ๐ก ) ยฏ ๐ โ (cid:205) ๐ ๐ ๐ / for ๐ in { , , , } do ๐ ๐ โ ๐ ๐ โ ยฏ ๐ ๐ ๐ โ max ๐ โ (|| ๐ ๐ ||) + max ๐ โ (|| ๐ ๐ ||) if ๐ ๐ = then return false ๐ ๐ ๐๐ โ computeSquaredDistance ( ๐ฅ , ๐ฅ , ๐ฅ , ๐ฅ ) ๐ โ ๐ ( ๐ ๐ ๐๐ โ ๐ )/( โ๏ธ ๐ ๐ ๐๐ + ๐ ) ๐ก โ ๐ก ๐ โ ( โ ๐ )( ๐ ๐ ๐๐ โ ๐ )/(( โ๏ธ ๐ ๐ ๐๐ + ๐ ) ๐ ๐ ) while true do for ๐ in { , , , } do ๐ฅ ๐ โ ๐ฅ ๐ + ๐ก ๐ ๐ ๐ ๐ ๐ ๐๐ โ computeSquaredDistance ( ๐ฅ , ๐ฅ , ๐ฅ , ๐ฅ ) if ๐ก > ( ๐ ๐ ๐๐ โ ๐ )/( โ๏ธ ๐ ๐ ๐๐ + ๐ ) < ๐ then break ๐ก โ ๐ก + ๐ก ๐ if ๐ก > ๐ก ๐ then return false ๐ก ๐ โ . ( ๐ ๐ ๐๐ โ ๐ )/(( โ๏ธ ๐ ๐ ๐๐ + ๐ ) ๐ ๐ ) return trueimpact (or can be simply set to 1), we can return no collision (lines19-20). Otherwise we compute a new local, lower bound, ๐ก ๐ , fromthe updated configuration (line 21) and begin the next iteration.ACCD thus provides an exceedingly simple-to-implement, nu-merically robust CCD evaluation. It requires only explicit calls todistance evaluations and so no numerically challenging root-findingoperations. In turn, ACCD is able to support thickness offsetting andalso controllable accuracy and so flexible tuning for performance vs.accuracy trade-offs in CCD applications. In Section 6.6 we compareACCD with state-of-the-art CCD solvers, there we see that theyall fail severely, resulting in intersection or tiny TOI that stalls the optimization in challenging examples. In turn, in easier cases wherealternate CCD methods do succeed, we then see that ACCD achievessimilar and often improved timing performance. Finally, we pointout that ACCDโs stopping criteria requires ๐ > ๐ = . Worst-case performance.
As discussed above, we find in practiceand particularly throughout all our challenging testing that ACCDconverges with solid performance (comparable or faster times whenother methods succeed, efficient times when all other methods fail).It is worth analyzing, however, that it should be possible to haveunbounded, worst-case performance. Recall we pick our referenceframe so that displacements sum to 0. However, if a primitive hasa diverging displacement field which cannot be cancelled out, ourdenominator ๐ ๐ can remain large. Then, if a primitive also has asmall starting distance ๐ and so correspondingly small numerator,iteration count for ACCD on this case can be large. While we seethat this scenario does not occur in practice for elastodynamics, ex-tensions of ACCD to accelerate convergence for these possible cases,and so obtain bounded performance guarantees is an interestingfuture step. We implement our methods in C++, applying CHOLMOD [Chenet al. 2008], compiled with Intel MKL LAPACK and BLAS for linearsolves and Eigen for remaining linear algebra routines [Guennebaudet al. 2010]. Necessary derivatives for C-IPC and our algebraic simpli-fications applied for efficiency and numerical robustness are detailedin our appendices. To enable future applications, development andtesting we will release our implementation of C-IPC as an opensource project. All our experiments and evaluations are executed oneither an Intel 16-Core i9-9900K CPU 3.6GHz ร ร ร Experiments.
Below we begin with a study illustrating and an-alyzing membrane locking behaviors for standard cloth materialsand meshes (Section 6.1). We then demonstrate C-IPCโs ability tostrictly satisfy strain-limits while fully coupling all physical forcesin Section 6.2. To our knowledge C-IPC is the first method to bothenable strict satisfaction of strain-limits and to support full couplingof strain-limiting, elasticity and contact forces. To analyze the im-pact of strain-limit satisfaction and full coupling we next considercomparison against prior methods. As discussed in Section 2, exist-ing methods in strain-limiting generate artifacts including locking,jittering and interpenetration. These issues result from two sources:1) inaccuracies from applied splitting models and 2) inability tosatisfy strain-limits in computation. Here, for the first time, we firstseparately analyze the problems that are created by splitting models(Sec. 6.3) independent of solver accuracy errors. Then in Section 6.4we consider the artifacts and inaccuracies (generated by state-of-theart cloth code) that jointly result from both splitting errors and (a) 1x sti๏ฌness (b) 0.1x sti๏ฌness(c) 0.01x sti๏ฌness (d) 0.01x sti๏ฌness + SL
Fig. 11.
Membrane Locking.
A cloth on sphere simulation (8K-node) with(a) ร stiff membrane of cotton material ( . ๐๐๐ ) suffers from membranelocking where bending behavior is artificially stiffened and creased. Reduc-ing membrane stiffness to (b) . ร that of measured cotton values removesartificial bending stiffness but we still observe overly sharp edge creasing.Next an even softer membrane at (c), . ร cottonโs, then avoids observablemembrane locking artifacts altogether, but results in nonphysical stretchinginappropriate for cotton. However, (d) applying a softer . ร membranetogether with our constitutive strain limiting model to enforce measuredcotton strain limits prevents membrane locking while obtaining naturalstretch limits for the material. inability of the method to enforce requested strain-limits. In Section6.5 and 6.6 we analyze C-IPCโs thickness modeling, comparing withboth existing cloth simulation codes and prior methods for CCDevaluation. Finally, with all components of C-IPC evaluated, weassess C-IPCโs application to garment simulation (Section 6.7) andits resolution of previous stress-test cloth benchmarks (Section 6.8).We then demonstrate C-IPC on simulations of fully coupled systemsof arbitrary codimension with consistent thickness modeling (Sec-tion 6.9) and finally consider its performance on a new set of clothsimulation stress tests designed to exercise robustness and accuracy(Section 6.10). Here we first illustrate and examine the impact of membrane lock-ing on real-world materials, and strain-limitingโs ability to mitigatethem. In Penava et al. [2014], density, thickness, and directionally de-pendent membrane stress-strain curves are measured and validatedfor a range of cloth materials. However, directly applying thesereal-world cloth parameters in simulation with standard-resolutionmeshes can produce severe membrane locking. This is unavoidableunless exceedingly high and so often impractically large mesh sizesare used. Here we examine this locking behavior, first using ourIPC model without strain-limiting as a demonstration. We also notethat such locking behaviors are independent of algorithm and so,for example, are also easy to demonstrate in other codes like ARC-Sim [Narain et al. 2012] with their real-world captured materialparameters [Wang et al. 2011]; see Figure 15d.We start by considering the behavior of a simple 1 ๐ ร ๐ un-structured mesh dropped on a fixed sphere and ground plane (withfriction of ๐ = . odimensional Incremental Potential Contact โข 13 density (472 . ๐๐ / ๐ ) and thickness (0 . ๐๐ ) from Penava et al.[2014] and then consider varying membrane stiffness while keepingbending Youngโs modulus at 0 . ๐๐๐ and Poisson ratios for bothmembrane and bending to the measured value for warp directionas 0 . . ๐๐๐ to 32 . ๐๐๐ for varying in-plane di-rections.For our example, we then find that even when applying an isotropicmembrane model, with the smallest determined membrane stiffness(0 . ๐๐๐ ), we observe severe locking artifacts (see e.g. Figure 11a)where bending is artificially stiffened. If we next lower this smallestmeasured membrane stiffness by 0 . ร , we then see that artificialbending artifacts are largely eliminated but we still obtain sharpcreasing artifacts forced by the dominating membrane energy (seee.g. Figure 11b). Next, if we try an even smaller 0 . ร scaling ofmembrane stiffness, observable membrane locking effects are nowgone, but as expected the resulting material is much too stretchedand so not even close to the desired material behavior (Figure 11c).This simple example nicely demonstrates the challenges of mem-brane locking when simulating stiff real-world cloth materials. Wethen note that these artifacts are only exacerbated in more challeng-ing simulations with, for example, moving boundary conditions andtight contact.Next we apply a strict strain-limiting with C-IPCโs isotropic modelto constrain strain within the elasticity range measured by Penavaet al. [2014]. The broadest range in all directions for cotton allowsa stretch factor up to 6 . . ร scaling of membrane stiffness, we now regain a simulationfree from membrane locking and, since restricted to measured strainlimits, also free of unnatural stretching artifacts (Figure 11d). Above we have demonstrated the well-known importance and im-pact of applying strain-limiting for simulating cloth materials. Herewe investigate C-IPCโs ability to accurately enforce tight strainlimits across both our isotropic and anisotropic models.
Accuracy across decreasing strain limits.
To confirm accuracy, con-trollability and robustness of C-IPCโs constitutive strain limitingwe test increasingly severe strain limits of 1% and 0 .
1% (well belowmeasured limits in standard cloth materials [Clyde et al. 2017; Pe-nava et al. 2014]). We apply them to both the sphere drape exampleconsidered in the last section as well as a simple membrane lockingtest with a two-corner pinned cloth [Chen et al. 2019; Jin et al. 2017].Here, applying C-IPC, see Figure 12, we observe stable simulationoutput with stretches on all triangles satisfying their prescribedlimit constraints. Closer examination of the draped sphere testsin Figure 12 also confirms that both are free from artificial bend-ing stiffness. For contrast, compare these results with the artificialstiffening in Figure 11a, where the maximum stretch exceeds 4 . unilaterally , this agreeswith Jin et al.โs [2017] argument that unilateral strain-limits canavoid membrane locking better than bilateral enforcement โ at leastto a certain degree. Even unilateral enforcement is not a perfectpanacea if limits are too tight. For an extreme example we can push (c) 1% (d) 0.1%(a) 1% (b) 0.1% Fig. 12.
Extreme Strain Limits.
We confirm C-IPCโs constitutive strainlimiting model provides robust and guaranteed satisfaction down to tightestmeasured and practically applied limits for cloth (and well below). Here wetest C-IPC with strain limits of % (already well below most measured limitsfor cloth) in (a) and (c). To push the challenge further in (b) and (d) we thentest C-IPC with an excessively small, largely nonphysical strain limit of . %.Here the applied simulation meshes are low resolution, containing respec-tively only 8K (for a and b) and 2K (c and d) nodes, and so are challengingto avoid locking artifacts with real-world materials. We confirm that C-IPCpreserves requested limits across all elements in all simulation time steps.We also note that, as expected, with the excessively small, nonphysicalstrain-limits applied in (d) more constraints become active w.r.t. to the DOFcount resulting in a small amount of visible membrane locking artifacts; seeour discussion below. strain-limits to a very small and largely unrealistic (0 . Anisotropic strain limiting.
For anisotropy, the story is similar.In Figure 13 we demonstrate C-IPCโs anisotropic model with boththe sphere drape and two-corner hang tests. Again we observethat applying measured stiffness values, here from Clyde [2017],generates clear membrane locking artifacts (see Figure 13a). Nextin Figures 13b and c we scale down membrane stiffness by 0 . ร and 0 . ร respectively. Throughout we confirm C-IPC enforces Clydeโs measured strain limits and again obtaining the expectedimproved results, both without locking artifacts and without non-physical stretching, at 0 . ร the measured stiffness. Similarly, inFigure 13d (again at 0 . ร reported stiffness) we apply C-IPC topreserve measured strain limits, obtaining artifact-free wrinklingwith the anisotropic model. Existing strain-limiting methods introduce errors from two primarysources. First, at the modeling level, they split strain limiting fromthe resolution of elasticity. Second, irrespective of the splitting modelemployed, algorithms applied to solve them are inaccurate andoften unable, as we will see in the following section, to satisfy evenmoderate requested strain-limits. All prior methods then introduceerrors from both of these sources and so it has remained unclearwhat problems originate from each source. Here, we first separatelyanalyze the problems that are created by splitting models, by solvingeach step of the split to tight accuracy. This allows us to show thatthe splitting itself introduces errors that are unavoidable irrespectiveof how accurately limit constraints could be enforced. Then, next inSection 6.4, we will examine solution accuracy and see that errors inthe strain-limit solve itself then produce inconsistent and so oftenuncontrollable results for simulation.To examine splitting errors we begin by applying the standardstrain-limiting, splitting strategy and so divide each time step solveinto two sequential steps. The first solves a predictor step thatincludes all forces for the whole system, with the exception of thestrain-limit , to obtain an intermediate configuration ห ๐ฅ . Next, thesecond step projects the predictor to satisfy the strain-limit. To do sowe minimize our constitutive strain limiting energy summed withthe mass-weighted ๐ฟ distance from the final timestep solution to ห ๐ฅ . (a) 1x sti๏ฌness with SL (b) 0.1x sti๏ฌness with SL(c) 0.01x sti๏ฌness with SL (d) 0.01x sti๏ฌness with SL Fig. 13.
Anisotropic strain limiting.
We apply the C-IPC anisotropicstrain-limiting model to the sphere drape test in (a) through (c), varying onlyscaling of the measured membrane stiffnesses, and to the two-corner clothhang test in (d). All are time stepped at โ = . ๐ with their final equilibriashown here. In (c) and (d), with a . ร scaling of the measured stiffness,we see that C-IPCโs tight enforcement of the measured strain limits [Clyde2017] obtains artifact-free drapes with differing wrinkle patterns resultingfrom the anisotropic model. In the simplest case we can consider the effect of this splittingwhen there are no contact forces. We start with a pinned and pulledsquare cloth. We fix its two top corners and apply heavy weights topull its two bottom corners downwards; see Figure 14a and b. Forcomparison, without strain limiting the cloth is stretched well over10 ร (not shown) while with our constitutive strain limits, strainsare restricted to the measured elasticity range and, in this fullycoupled solution, we see the expected vertically aligned wrinklesfrom stretching at the two bottom pulled corners; see Figure 14a.On the other hand, in Figure 14b we see that splitting strain lim-iting from the elasticity solve produces obvious biasing artifactsat the two bottom corners where the wrinkles are aligned in non-physical directions. These errors, resulting from the decoupling ofstrain-limiting forces, are then increasingly severe as time step sizesincrease (here we show with โ = . .
7) then we indeed find the resulting splitting solution is nowfree from severe compression artifacts. Now however, due to theerrors in membrane and bending, the cloth in the split solution stillremains unnaturally flat against the floor; see Figure 14e.Of course as with all time-step splitting methods, splitting errorscan be reduced by applying increasingly smaller time step sizes.Here, for example, we find that visually noticeable errors betweenthe split model and our fully coupled solve disappear at โ = . ๐ .However, as expected, the resultingly large and often unaccept-able increase in compute time wipes out any expected gains fromsplitting (we will see this theme again in the next sectionsโ compar-isons with existing cloth codes). Likewise, the necessary decreasein time step size then varies with example and scene so that robust,controllable time-stepping with splitting, does not appear possible. Above we analyze the modeling errors generated by time-step split-ting with strain limiting. To do so we compare models for whicheach stage is consistently solved to tight accuracy in our commonbenchmark code. Here we next consider the joint impact that split-ting errors combined with inaccurate constraint solves then have inexisting cloth simulation codes. To do so we compare with ARCSim[Narain et al. 2012] which is, to our knowledge, the most robust odimensional Incremental Potential Contact โข 15 (a) C-IPC (b) split (c) C-IPC (d) split (e) split plus lower bound
Fig. 14.
Comparison of C-IPCโs fully coupled strain limiting model with standard splitting models highlights the differences between the fullycoupled solutions in (a) and (c) and the artifacts generated by splitting strain-limiting solves from elastodynamics in (b), (d) and (e). Here we see that withsplitting resulting large errors in elasticity (membrane and bending) can generate incorrect wrinkling directions (b), severe compression artifacts (d), or evennumerical softening as in (d) and (e). (a) C-IPC (b) ARCSim (c) C-IPC (d) ARCSim 1x sti๏ฌ (e) ARCSim 0.1x sti๏ฌ (๎ ) ARCSim 0.01x sti๏ฌ
Fig. 15.
Comparison of C-IPC and ARCSim results illustrate the different behaviors of C-IPCโs fully coupled constitutive strain limiting solutions in (a)and (c) and the output of ARCSimโs split model and augmented Lagrangian strain-limiting solver in (b) and (d-f). Here frames from ARCSim near equilibriumexhibit excessive stretch due to lack of limit enforcement in (b), (e) and (f). At default material settings in (d) ARCSim exhibits the expected locking artifacts.Then, as we decrease ARCSimโs membrane stiffness to avoid membrane locking in (e) and (f) its strain-limiting can not preserve bounds, while the split modelintroduces additional artifacts so that the over-compression and stretching artifacts are now most obvious. For reference compare with C-IPCโs strain-limited,fully coupled solutions in (a) and (c) at comparable stiffness ( . ร ) and same times and see Figure 16 for summary of timings and limits achieved. shell simulator currently available to support strain limits with real-world cloth material parameters. See Figure 16 for statistics on thestrain limits and timings achieved by ARCSim and C-IPC here.We again begin by considering the simpler, contact-free case. Todo so we consider an even simpler cloth hang example: by droppingthe pinned cloth without weighted bottom corners. We apply ARC-Simโs default cotton material and default strain limit settings whichrestrict to the range [ . , . ] . Here, even with this wide range ofpermissible strains and at full membrane stiffness, ARCSim resultsviolate the limit bounds significantly at each time step. For visualreference in Figure 15a we illustrate the fully coupled C-IPC so-lution at equilibrium satisfying the strain-limit bound. Comparewith the ARCSim output at equilibrium in Figure 15b where weeasily see the artifacts resulting from over-stretching near the fixedcorners forming a larger arc. Here both frames are taken at 4 ๐ andare time-stepped with โ = . ๐ .As we then decrease ARCSimโs time-step size to โ = . ๐ theseerrors still remain significant, and it is not until we reach at stepsize of โ = . ๐ that strain limits are mostly (although still notentirely) satisfied. However, to do so ARCSim then requires 87 min-utes to simulate this simple 4 ๐ simulation sequence, resulting in awell-over 10 ร slowdown compared to C-IPCโs fully coupled, strain-limit-satisfying solution stepped at โ = . ๐ . However, beyondthe requisite slowdown, the required decrease in step size changesper scene and material and so it is always unclear, without manytime-consuming simulation tests per scene, how to determine thenecessary time step decrease to ensure ARCSim results satisfy a pre-scribed strain limit. In contrast C-IPC again maintains strain-limitconstraint satisfaction across time step, material and scene settings. To investigate ARSimโs constraint errors further we also observethat the augmented Lagrangian strain-limiting solver [Narain et al.2012] employed in ARCSim applies a hard-coded max-iteration capset to 100. Experimentally increasing this limit (and otherwise leav-ing the ARCSim code unchanged) we confirm that even increasingthis cap by an order of magnitude, strain limit satisfaction is stillnot achieved.Next we add contact for our ARCSim comparison and considerthe 8K-node sphere drape test. Here we again observe that applyingstrain-limiting to reduce membrane locking does not work and in-stead see that strain-limit satisfaction worsens and artifacts actuallyincrease as membrane stiffness decreases.Specifically we start by applying ARCSimโs default cotton ma-terial and the same default strain-limit settings as above. For thisstiff material we observe the expected membrane locking issues;see Figure 15d. However, we do note that with this stiff membrane,strain limit constraints are certainly reasonably well satisfied exceptfor a few time steps.Next, as standard, we reduce membrane stiffness in ARCSimwith hope of mitigating the locking with the expectation that strain-limiting will compensate. If we reduce stretching stiffness by 10 ร , thelocking artifacts are indeed less; however, ARCSimโs strain-limitingdoes not provide the intended compensation for the reduced stiff-ness. Instead, there are significantly more triangles violating thestrain limit throughout the simulation (see Figure 15e) resulting instretching artifacts. If we then further reduce stiffness to 0 . ร thedefault cotton value (generally necessary in this example to removeall visible locking artifacts; e.g. compare with our comparable stiff-ness, C-IPC result in Figure 11d and 15c) ARCSimโs results in Figure Test Method h (s) Max s max
Avg s max
Min s min
Avg s min
Time (m)
Hang X C-IPC
Hang X ARCSim
Hang X ARCSim
Hang X ARCSim
Sphere X C-IPC
Sphere X ARCSim
Sphere X ARCSim
Sphere X ARCSim
Fig. 16.
Strain limiting comparison with ARCSim.
We compare ARC-Sim and C-IPC strain limiting on two cloth simulation test cases: a pinnedhanging cloth (without contact) and fixed sphere drape (with contact). Forboth methods we apply ARCSimโs default strain limit settings targetingstrains to [ . , . ] . Above we summarize statistics of the achieved max, ๐ max , (respectively, min ๐ min ) strain-limit over all elements in a timestepas the averaged and max (respectively min) over all time steps. We applyC-IPC stepped at โ = . with stiffness rescaled by . ร (necessary toavoid membrane locking artifacts and so requiring accurate strain-limitingto compensate) for both examples. We correspondingly consider ARCSimacross a range of time-step sizes and stiffness re-scalings ( ร , . ร and . ร ). See Figure 15 for visualization of results. In the hang test, even forthe full membrane stiffness, we see ARCSimโs strain limits are far fromsatisfied. Then, as step size decreases we observe ARCSim results bettersatisfy strain limits โ although never fully. Correspondingly, in the spheredrape, as we reduce membrane stiffness towards a value required to mitigatelocking, we ARCSim results increasingly violate strain limits more signifi-cantly. In comparison we observe C-IPC satisfies strain limit satisfaction atthe necessary soft membrane energy, across all elements and time steps forboth tests. Next we examine modeling small but finite thicknesses in contact.To illustrate differences between the C-IPC model and existing state-of-the-art shell simulators we examine a simple, yet challengingcloth stacking benchmark. We construct a scene with a stack often 8K-node square cloth panels aligned with vertical spacing androtated orientation, dropped simultaneously onto a fixed squareboard. See Figure 17a for initial set up. This should form a cloth pilewhere as we vary material thickness the overall pile height and bulkproperties should also correspondingly change. (a) initial frame (b) ARCSim 5mm (c) ARCSim 10mm(d) Argus 1mm (e) Argus 5mm (f) Argus 10mm
Fig. 17.
Cloth stack comparisons.
Ten 8K-node square cloth panels aredropped in stack upon onto a square board with friction ๐ = . (a). HereARCSim [Narain et al. 2012] and ARGUS [Li et al. 2018a] are tested forvarying cloth thicknesses. Time stepping at โ = . ๐ ARCSim fails toresolve contact for a thickness of ๐๐ (b), and at ๐๐ there is no longerinterpenetration but the dynamics and geometry have noticeable artifacts(c). ARGUS can successfully simulate cloth with ๐๐ and ๐๐ thicknesses(e,f) at โ = . ๐ with some smaller artifacts; but at ๐๐ thickness it failsto resolve contact altogether. See the top row Figure 18 for our correspondingIPC results for this example. ARCSim [Narain et al. 2012] computes collision response apply-ing Harmon et alโs [2008] non-rigid impact zones and a repulsionthickness applied to both model cloth thickness and to stabilizecollision-processing constraints. ARCSimโs thickness parameter de-fault is set to 5mm and is exposed as a configuration that can bechanged by users. This parameter is then associated in code withadditional thickness parameters directly applied for collision pro-jection, analysis and contact force computation. Starting with thedefault 5 ๐๐ thickness setting, and testing with increasingly smallertime steps down to โ = . ๐ , ARCSim consistently reports col-lision handling failures and demonstrates severe artifacts in thesimulation results; see Figure 17b. Next applying a 10 ๐๐ thicknessARCSim reports success in resolving contacts at โ = . ๐ , butstill generates large artifacts in both geometry and dynamics; seeFigure 17c. We did not push ARCSim further to time steps smallerthan โ = . ๐ , as this already makes ARCSim excessively slow forcomputation; e.g. as compared to C-IPC. ARGUS [Li et al. 2018a] updates ARCSim with improved con-tact and friction processing for shell simulation. As in ARCSim italso applies and exposes to users the same thickness parameters.At โ = . ๐ ARGUS is able to simulate the cloth stack with thedefault 5 ๐๐ material thickness (Figure 17e) and also with a thicker10 ๐๐ material (Figure 17f). Here we also observe that some heightdifferences in the piles are captured. However, as we decrease thick-ness towards thinner, more standard garment material thicknesses,e.g. setting thickness to 1 ๐๐ , ARGUS produces severe artifactsas shown in Figure 17d. As with ARCSim, for all three thicknesssettings ARGUSโ timing remains much slower than C-IPC.In summary we observe that state-of-the-art methods requiresmall and often impractical time step sizes to avoid failures as mod-eled geometric thickness decreases. With smaller thickness and/orincreasing collision speeds, the required time step size to avoidfailure decreases, correspondingly increasing simulation cost and We use the most recent ARCSim v0.3.1 in our testing. odimensional Incremental Potential Contact โข 17 making it unclear what time step size is required for success perscene and thickness.
C-IPC in contrast, for the same thickness benchmark consistentlymodels the varying thickness by simply setting ห ๐ (recall the distancewe start exerting contact forces at) to effective thickness values 1 ๐๐ ,5 ๐๐ , and 10 ๐๐ . In the first row of Figure 18 C-IPC successivelymodels all varying material thicknesses and corresponding heightsof the cloth stacks. C-IPCโs model then further provides constitutivebehavior for thickness in the normal direction. This is demonstratedin middle and bottom rows of Figure 18 where we then drop anelastic ball on the stack. In the middle row we show maximumcompression for each stackโs collision, observing the increasingindentation with thickness. Also note the wrinkles forming solely inthe thickest, 10 ๐๐ case. In the bottom row we show the final frameat equilibrium for each simulation after collision, illustrating thedifferent rest heights of the shell piles weighted by the volumetricFE ball model. This example also illustrates C-IPCโs natural couplingof differing elastic models which we explore in detail in Section 6.9. For the last sectionโs comparisons it was sufficient to model thick-ness effects in C-IPC solely via ห ๐ and so purely as an elastic behavior.More generally, as discussed in Section 5, as thicknesses becomesmaller and contacts tighter, C-IPCโs inelastic thickness model uti-lizing both ห ๐ and our offset ๐ is required.For simulating codimensional objects these two mechanisms arecombined together. Here larger ห ๐ introduces increased elastic re-sponse while a nonzero ๐ provides guarantees for minimum thick-ness even under extreme compression; see e.g. Figure 2. For sceneswhere consistently modeling thickness matters; see e.g. Figures 2, 3,10, 21, 22, and 23, we then set ๐ to the reported material thicknessvalue of the object or slightly smaller and then set ห ๐ near ๐ to controlthe degree of elastic response required.As discussed in Section 5.3 simulating these examples producesmuch more degenerate contact pairs and requires higher precision Fig. 18.
Sphere on cloth stack.
Ten 8K-node square cloth are droppedonto a square board with friction ( ๐ = . ). Top: Varying IPC ห ๐ per column,C-IPC is able to provide controllable thickness modeling for cloth. Middle:A soft elastic volumetric sphere ( ๐ธ = ๐พ๐๐ ) is dropped onto to the clothstacks and is illustrated here at maximum compression; here with ๐๐ thickness (far right) we observe the resulting wrinkling enabled by thethicker depression. Bottom: at rest equilibrium heights continue to varywith thickness under the weight of the ball. Example ACCD IPC Root Parity MinSep BSC T-BSC
Sprinkles โ tiny toi โ โ Noodles โ tiny toi โ โ Hair clusters โ tiny toi โ โ Braids โ tiny toi โ โ Card Shu ๏ฌ le โ tiny toi โ โ Needle bed โ Kick
Stack, 1mm โ Stack, 5mm โ Stack, 10mm โ Fig. 19.
CCD comparison statistics.
Here we summarize the per-Newtoniteration, total CCD query time (excluding spatial hash) in milliseconds forten simulation benchmark examples across CCD methods. We compareour ACCD method, IPCโs conservative floating point CCD [Li et al. 2020a],MinSep [Harmon et al. 2011; Lu et al. 2019], Root Parity [Brochu et al. 2012],BSC [Tang et al. 2014], and T-BSC [Wang et al. 2015]. First five examples inthe benchmark include finite thickness offset, which can not be supportedby Root Parity, BSC, and T-BSC methods. When simulation cannot proceedfrom CCD errors we report type. Here โ โ โ indicates unable to find non-intersecting step (endless line-search) and โtiny toiโโ too small a step size(no progress). to capture thickness offsets and so makes great demands on CCDaccuracy. Here the conservative CCD strategy from IPC [Li et al.2020a] is no longer sufficient. Employing standard floating pointCCD routines can and will simply return 0 TOI for many queries,and so incorrectly stop the IPC optimization progress altogether.Likewise alternative CCD methods with higher reported accuraciessimilarly fail in most such cases.This motivates our new ACCD method which, as demonstratedbelow, remains robust and accurate โ obtaining stable, forwardprogress in IPC-optimization for even our most complex exampleswhere all available CCD alternatives fail. At the same time ACCDcontinues to provide comparable and generally faster timings oneasier examples where alternate CCD methods are able to succeed.For comparison with ACCD we test a range of CCD methods:Minimal Separation (MinSep) CCD [Harmon et al. 2011; Lu et al.2019], Root Parity [Brochu et al. 2012], BSC [Tang et al. 2014], T-BSC[Wang et al. 2015] and IPCโs conservative floating point CCD[Liet al. 2020a] on a set of ten benchmark examples with challengingcodimensional collisions. The first five examples in the benchmarkhave thickness offsets and the remaining five do not; see Figure 19.We are able to directly test MinSep by applying it as a base solverwithin the same conservative CCD strategy as IPC [Li et al. 2020a].Here minimal separation is set to ๐ ร ๐ cur with the current scalingfactor ๐ and distance ๐ cur . If the returned TOI is less than 10 โ thequery is performed again without scaling ( ๐ =
0) in order to attemptto make the query easier. To test the Root Parity, BSC and T-BSCmethods, which only decide whether an interval "has collision" or"has no collision" without computing TOI, we apply them within amonotonic bisection for each query until finding a step size withoutcollision. When found, we then return a conservative scaling of the
Example Method Nodes h (s) Time (m) SL Contacts (avg., max.)
Cores ยต
Arabesque C-IPC
C-IPC
Argus
Clubbing C-IPC
C-IPC
Argus
Rhumba Layer
C-IPC 43.3K 0.04 1.2 1.134 31.7K, 43.7K 6 0.2
Kick Layer C-IPC
Rhumba Knife
C-IPC 12.8K 0.04 0.4 1.134 15.4K, 21.3K 6 0.2
Fig. 20.
Garment simulation tests.
A summary of garment simulationstatistics a pair of test dance sequences, Arabesque and Clubbing, fromARGUS and three additional tests with garment patterns not handled byARGUS. Here SL denotes the strain limit applied. step by 1 โ ๐ . As bisection cannot directly extend these methods tocompute times or step sizes that first bring a distance of a contactpair to ๐ , we can not test these three methods on our examples withthickness offsets; and so restrict our testing of them to just the lastfive test examples in Figure 19).We summarize our comparisons in Figure 19. IPCโs conservativestrategy with floating-point CCD works well on simpler scenes,albeit with slightly slower timings than our ACCD method, but thenfails on the even more complex ones. We also note that while con-servative floating point CCD handles degenerate cases like paralleledge-edge by querying additional point-edge and point-point pairs,ACCD only needs to query general pairs (e.g. point-triangle andedge-edge pairs for surface-only scenes). On the other hand, MinSepalways returns tiny values, sometimes even 0, for all ten benchmarkexamples and so makes no progress. Root Parity performs well withnice efficiency on most benchmark examples it can be applied to(those without thickness offsets) but unfortunately misses somecollisions altogether leading to unacceptable interpenetrations inexamples like Needle Bed and Garment. For BSC we find it reportsruntime errors for all benchmark examples it can be applied to withan internal message "Inflection points are not handled exactly inBSC", while T-BSC returns tiny values for the Garment example, andis unable to determine step sizes without collision (we confirm bothBSC variants work well on simple cases) and so becomes trapped inline search in the other five benchmark examples it can be appliedto. In contrast we see that ACCD enables all examples to completeto convergence while, for examples where alternate CCD routinesare able to complete, it offers fastest times. Finally, as a preliminaryinvestigation of ACCD application outside of the C-IPC frameworkwe apply ACCD in place of ARGUSโs default CCD routines and testwith both 1mm and 5mm cloth stack examples. While, as expected,the simulation results are nearly identical (with no visual difference)CCD costs decrease by โผ We next apply C-IPC to simulate garment drapes and dynamics.
ARGUS Comparisons.
We start by considering C-IPC on a pairof test dance sequences, Arabesque and Clubbing, proposed forARGUS [Li et al. 2018a]. Both examples use the highest resolutioninput meshes tested in Li et al. [2018a]. Here Clubbing exercisescontact processing more extensively, with faster motions and tran-sitions than Arabesque. For comparable results we run C-IPC andre-run ARGUS on the same machine, with fixed mesh resolutionfor ARGUS. Here we observe robust, smooth simulation resultsacross the entire sequences for C-IPC both with and without strainlimits. Notably, even with the addition of strict non-intersectionenforcement, convergent, fully implicit time steps and resolving allcontact stencils, not solely node-node as in ARGUSโs treatment, westill observe 3.7X and 1.5X speedup for C-IPC on Arabesque andClubbing over ARGUS results. Likewise, if we also enable C-IPCโsconstitutive strain-limiting, speedups remain similarly at 2.7X and1.4X. See Figure 20 for comparative timings and our supplementalvideo for results.
Garment Simulation Challenges.
Two primary challenges in gar-ment simulation are resolving complex, multilayered designs andsimulating clothing worn by rapidly moving characters. To test therobustness of C-IPC for these cases we construct two new examples.Due to their complex stitching patterns these garments can not besimulated in the ARGUS and ARCSim codes.We simulate the yellow dress (15.7K nodes) produced by FoldS-ketch [Li et al. 2018b] with challenging-to-resolve knife pleats cre-ated where the cloth is folded upon itself and then stitched tightin contact. We start C-IPC with the flat pattern staged around the12.8K-node mannequin (Figure 5a). We then drape by stepping C-IPC forward with the stitching forces (springs) seaming the designtogether at large time step. C-IPC obtains the static equilibriumof the draped garment after just a few large time steps obtainingsharp knife pleats while preserving the correct layering order, freeof intersection (Figure 5b). We then move the knife pleat throughthe steps of the Rhumba, preserving the flowing motion of the dresswith knife pleats remain crisp and properly colliding throughout(Figure 5c).We then apply C-IPC to similarly initialize a multilayer skirtpattern (30.5K nodes) from Sensitive Couture [Umetani et al. 2011].This obtains the initial pose static drape in Figure 5d. Next we scriptthe mannequin mesh dressed in the layered garment through a fast-moving martial arts sequence with large motions. In the bottom rowof Figure 5 we show three frames of the resulting C-IPC simulationfeaturing a portion of a rapid kicking sequence. Please also see oursupplemental video. C-IPC stably resolves these rapid movements,continuing at a large time step of โ = . ๐ , while capturing intricategarment details, enforcing strain limits and non-intersection (andso preserving the layering order) throughout. Finally we take thesame garment through the slower Rhumba sequence and observea correspondingly reduced cost as less iterations are required toresolve the time steps with the gentler motion and correspondinglylower contact count. odimensional Incremental Potential Contact โข 19 (a) initial (b) static frame (c) static frame (lower view) (d) static frame (segment view) Fig. 21.
Sprinkles.
We drop 25K ๐๐ -long and ๐๐ -thick โsprinklesโ into a bowl (a). Setting IPCโs offset and ห ๐ to 1.5mm and 0.5mm respectively to modeleach sprinkleโs thickness C-IPC fills the bowlโs volume at rest (b). In (c) we remove the bowl to highlight the non-interpenetrating rest configuration while in(d) we remove volume rendering as well to show the underlying edge segments simulated. In this section we confirm C-IPCโs performance on challengingbenchmark tests from prior work.
Cloth on rotating sphere.
In Figure 6 we test spinning sphere-typeexamples from Bridson et al. [2002] designed to exercise friction,contact processing and tight wrinkling. We start by dropping an86 ๐พ node square cloth (unstructured mesh with ๐ = .
01 and 1 . ๐ = . ร ) higher frequency wrinklesresolve in the folds. Funnel.
To test layering, friction and tight contact with strainlimiting we extend the funnel examples from Tang et al. [2018] andHarmon et al. [2008]. In Figure 7 we drop three cloth panels ( ๐ = . . Reef knot.
To combine extreme compression, contact and frictionwith large stresses we match the initial set up of the reef knot exam-ple from Harmon et al. [2009]. Two ribbons are initially intertwinedand then stretched to form a knot in Figure 4. Here we extend thechallenge of the original test by simulating with 100K nodes (10 ร that of the original), applying a strain limit of 1 . ๐ = .
02, and pulling much tighter โ stretching the ribbons to nearlytheir strain limit. This captures the frictional dynamics of the tra-jectory (without jitter) and results in extreme stress at the center,forming a much smaller knot than originally demonstrated.
In this section we demonstrate C-IPC simulation across all codomains.This covers volumetric FE materials, shells, rods and even, whendesired, particles. Because these domains interact via IPC barrierswe are able to seamlessly and directly couple all such mesh-based models in the same simulation. These tests also highlight the advan-tages and suitability of our consistent and controllable thicknessmodeling for these thin material interactions.
Card shuffle.
We test accuracy and precision in very thin shellfrictional contact during a complex, interleaving shuffle. We beginwith fifty-four stiff shell playing cards (63 . ๐๐ ร . ๐๐ ) withIPC offset and ห ๐ set to 0 . ๐๐ and 0 . ๐๐ respectively to resolvecard thickness. In Figure 8a we break the cards into two piles andapply moving boundary conditions on their sides to bend themin preparation for a โprecisionโ bridge shuffle. Unlike a human-performed bridge finish, where cards begin interleaved, here weincrease challenge by keeping the two bridged piles apart. We thenrelease the held boundary conditions card-by-card from the bottomup, alternating between left and right side piles (Figure 8b). Releasedcards quickly regain flat rest shapes and fall on the ground to forma new, shuffled, intersection-free stacked pile. We then partiallysquare the deck with two scripted collision boundaries (Figure 8c)obtaining a final pile 15.4 mm high and so well matching that of areal deck (Figure 8d). Hairs.
Modeling hair with self-collision is a severe challenge forrod contact processing โ in part due to their extremely small cross-sections. We examine C-IPCโs behavior on two challenging hairsimulation tests. Hair thickness is small, but its bulk volumetriceffects are important and have remained difficult to capture. HereC-IPC captures hair thickness by enforcing an offset of 0 . ๐๐ for all rods with ห ๐ = . ๐๐ . In Figure 3(a)-(e) we first test C-IPCwith braiding. We form a tight braid by twisting the bottom DOFsof two clusters of hair about each other. Despite undergoing largestress the resulting braid remains intersection free and preserve thethickness offset throughout. We also confirm this in dynamics byreleasing the bottom held boundary conditions and letting the braidunwind without snagging artifacts nor instabilities. Next, we testhair self-collision and friction further by dropping a cluster of hairswith one free end across a similar perpendicularly hanging cluster[McAdams et al. 2009] and over a fixed sphere in Figure 3(f)-(k).C-IPC again simulates this scene with convergent, large time stepsolves throughout. (d)(a) (b) (c) (e) ( f ) (g) Fig. 22.
Granules on cloth.
We drop a column of 50K particles onto a cloth with fixed corners. Each particle is represented by a single node, with C-IPCoffset and ห ๐ set to ๐๐ and ๐๐ respectively in order to capture the granuleโs volume at rest while the C-IPC barrier ensures particles do not tunnel througheach other nor the cloth shell. Particles at rest capture tight spherical packings (see zoom in) confirming offsets are preserved throughout. Finally we releaseone corner of the cloth to let the granules fall onto the ground. Noodles.
A large bowl of noodles is delicious but also challengingto simulate due to the large numbers of DOF and high contact counts.In Figure 10a we drop a stack of 625 40 ๐๐ -long cooked noodles intoa 13 . ๐๐ -diameter bowl (fixed geometry). Each noodle is a 200-segment discrete rod with ๐ = . ๐ set to 1 ๐๐ and0 . ๐๐ respectively to model the noodlesโ 1 . ๐๐ thickness. WithC-IPCโs thickness handling the dropped rods come to a rest (theirtotal volume should be 4 . ร โ ๐ ) largely filling the bowl (Figure10b) whose volume is 6 . ร โ ๐ . At equilibrium In Figure 10cwe remove the bowl to show the final intertwined, interpenetration-free configuration. Finally in Figure 10d we remove the volumerender to highlight that this complex non-intersecting volume issimulated solely with discrete rod polylines. Sprinkles.
Similarly we can do the same for just edge segments.We drop a grid of 25K โsprinklesโ each formed by 6 ๐๐ -long, singleedge segment into a bowl (Figure 21a). The sprinklesโ 2 ๐๐ thicknessare modeled by a C-IPC offset of 1 . ๐๐ and ห ๐ of 0 . ๐๐ . Making theoffset 3 ร larger than ห ๐ models a strongly inelastic response for theircollisions. As with the noodles we can estimate the total expectedvolume of the sprinkles as 5 . ร โ ๐ which at equilibrium fillthe bowlโs 6 . ร โ ๐ volume without interpenetration (Figure21b and c). Again in Figure 21d we remove the volume rendering tohighlight the geometry is simulated solely with codimensional linesegments. Spherical Granules on cloth.
In the extreme C-IPC can even sim-ulate small spherical granules by treating individual vertices asparticles endowed with a strict thickness offset. We test resolu-tion of non-intersecting particle trajectories by dropping a columnof 50K particles onto a cloth with fixed corners (Figure 22). Col-lisions remain non-intersecting even at high-speed and with thecloth. Particles at rest capture tight spherical packings (see zoomin) confirming offsets are preserved. Once particles pile in the cen-ter we release a cloth corner letting them flow to the ground. Tomodel each particle we compute mass with ๐ = ๐๐ / ๐ and ๐ = / ๐๐ โ . ร โ ๐ , and set thickness with a C-IPC offset of2 ๐๐ and ห ๐ of 1 ๐๐ . Simulating particles this way then resembles theDiscrete Element Method (DEM) [Cundall and Strack 1979; Yue et al.2018] where particles are โcoatedโ with an interaction zone (analo-gously ๐ < ๐ < ๐ + ห ๐ in our setting) for contact forces determinationbased on the depth, area, or volume of the zonesโ intersections [Jianget al. 2020]. DEM is a critical tool in geomechanics where it pro-vides a reasonable simplification by ignoring particle deformationwhile accurately capturing momentum and inter-particle contact forces. Here C-IPC relates contact forces directly to distance via ourbarrier function and so enables seamless coupling with deformablematerials as in the cloth here. Unlike DEM, C-IPC then also addition-ally guarantees intersection-free trajectories for particle volumesdefined by offsets. All-in.
Finally we test coupling all codimension together. We dropan FE armadillo volume, a shell cloth and a particle sand columnonto a net formed by interleaved rods (Figure 23a), applying a C-IPCoffset of 0 . ๐๐ and ห ๐ of 1 ๐๐ for all domains. We successively dropthese codomains one at a time. First the dropped armadillo is tangledand caught between the netโs rods (Figure 23b). Next the cloth fallsupon the armadillo (Figure 23c) and finally the sand column poursonto cloth and flows into the two resulting folds of cloth supportedby the armadillo and the rod net (Figure 23d, e, and f (top view)). Finally we propose and evaluate C-IPC on a set of three new, chal-lenging cloth simulation examples designed to require extreme ro-bustness and accuracy in modeling elastodynamics with frictionalcontact, strain limiting and thickness. To our knowledge C-IPCis the first method to achieve simulations of these scenarios andbehaviors.
Pulling cloth over codimensional needles.
Cloth contacts againstsharp geometries notoriously stress contact simulations with snag-ging and sedropping and then dragging a cloth panel across a needlebed formed by line segments (Figure 9). With a large time step sizeof โ = . ๐ and strain limit of 1 . ๐ / ๐ across the segment tips and observe rapidsliding along and then off the top of the bed of needles withoutsnagging nor stretching artifacts (strain limits are satisfied through-out).vere stretching artifacts. Likewise, if these contact are slidingacross asperities with large displacement (i.e., large time step and/orlarge speed) these challenges only grow. Here we test C-IPCโs han-dling of extreme cases with both challenges, by In cases like these,were we have extreme stretch from sharp contacts, C-IPCโs adaptivestrain-limiting stiffness is especially important to ensure that timestep solves remain numerically feasible and efficient. Tablecloth trick.
The classic tablecloth trick attempts to pull acloth out from under a table setting with the aim of keeping thedinnerware upright and ideally staying mostly in place. Key tomaking this trick work is a pull with sufficient speed so that slidingacceleration of the cloth can overcome friction forces with little or odimensional Incremental Potential Contact โข 21 (b) (c) (d) (e)(a) initial ( f ) Fig. 23.
All-in.
We simulate objects of all codimensions coupled together via C-IPC contact in a scene with an armadillo volume, a shell cloth, and a sandcolumn formed by individual particles released in sequence to fall onto a net formed by interleaved rods. no pull on the settings . Simulation of this effect then requires tightcoupling of elasticity and strain limiting with accurate resolutionof friction and contact as well as robust handling of sliding contactwith sharp edges. This is needed to correctly model the balancesbetween friction, the large stresses and forces from pulling, andthe normal forces from heavy dinnerware. In Figure 1 we set atable for C-IPC with both heavy and stiff dinnerware (volumetricFE models). We simulate pulling the strain-limited (1 . โ = . ๐ , with increasing speeds.As we start with lower pulling speeds, as expected, objects do notfollow the cloth. As we increase speed we then get less messy resultsand finally, with a pull at 4 ๐ / ๐ , the full setting remains upright andon the table with just a bit of movement (we also note detailedfolds begin to be formed by the strain-limited cloth pulled at highspeed). Finally, with a high-speed pull of 8 ๐ / ๐ the tablecloth pullsout smoothly with almost no change to the setting at all. Here,for the latter two successful pulls we observe extreme forces fromthe moving boundary conditions pulling the cloth at a high speedand the large friction forces exerted by contact with the heavierdinnerware. Twisted cylinder.
Finally, we test thickness modeling, contact res-olution and buckling under extreme and increasing stress. We startwith a 1 ๐ -wide cloth cylinder (0 . ๐ radius) modeled by a shellwith 88K nodes, with a thickness offset of 1 . ๐๐ and ห ๐ of 1 ๐๐ , time-stepped at โ = . ๐ . We then simultaneously twist the cylinder andslowly move the two sides together at 72 โฆ / ๐ and 5 ๐๐ / ๐ respectivelyto introduce wrinkling, folding, and eventual tight buckling. In orderto allow this scripted torture to continue for the entire 38 secondslength of this sequence we do not apply strain limiting โ so C-IPCmust resolve exceedingly high stretching forces here. Likewise, toclarify buckling behavior we do not apply gravity (avoiding droop).In the first second of simulation we immediately obtain interestingglobal folds (Figure 2a). Soon after a thick central cylinder of woundcloth forms. The volume of this cylinderโs structure is supported byC-IPCโs finite thickness offset (Figure 2b). For contrast consider inFigure 2c a frame captured at the same time step but now simulated without C-IPCโs offset. Here we clearly see that correct geometry cannot be formed without C-IPCโs thickness offset (nor can it capturethe materialโs later buckling behavior โ see below) further clarifyingthe importance of consistent thickness modeling for codimensionalmodels. Next, in Figure 2d as we continue our C-IPC simulation with offset for 32 . ๐ of further twisting, our offset-thickened cloth continues to support complex behaviors including the final buckledgeometry in Figure 2e. C-IPC extends the incremental potential contact framework fromvolumetric deformables [Li et al. 2020a] to codimensional and mixed-dimensional structures with controllable thickness, fully coupledstrain limiting, accurate frictional behavior and non-intersectionmaintained throughout all time steps. Of course as in the originalIPC method, in order to provide these non-intersection guaranteesIPC requires an initial non-interpenetrating geometry at start ofsimulation. Currently, e.g. for garments, this can often be easilyachieved by staging. However, loosening these constraints to achieveguarantees from a starting tangled state is an important and excitingdirection of future exploration.Along with enabling new and improved applications in cloth,hair and many other codimensional simulation tasks we are alsoexcited to extend C-IPC further for geomechanics applications inthe simulation of granular flows โ especially for complex granuleshapes which have so far proven challenging to resolve.Here we have first focused on providing a reliable simulator forall materials and conditions. As noted in previous works, in manyapplications the biggest overhead is the many laborious simulationpasses required to hand-tune collision parameters to make a scenework at all, and then often many more time-consuming steps tostill correct output. As such, reliability leads to overall increasedspeed of output. However, although performance for C-IPC, as an-alyzed above, is more than competitive with state of the art clothcodes when they can offer comparable accuracies, much more canbe done to improve C-IPC performance. We look forward to bothbetter leveraging parallel architectures and improving underlyingIPC solver methods. Similarly, here we are investigating and exten-sively testing with implicit methods. Another fruitful direction ofexploration is explicit time stepping. It would be exciting to con-sider relative behavior of C-IPC and ACM [Harmon et al. 2009]for side-by-side comparison between state of the art implicit andexplicit methods providing non-intersection guarantees. Likewise,we look forward to leveraging the differentiability of our approach.Here we expect that the combined benefits of C-IPCโs reliable per-formance across sweeps of input and its differentiability, shouldtogether enable many design and training tasks for garments andcomplex structures. Finally, here we have focused on providing Towards this end we have been communicating with the ACM authors for over yeartowards obtaining a representative code for comparison. After hard work on both endswe are hopeful, but licensing concerns required for access to the code have hamperedus so far. We hope to have the comparison in future work.
Examples ๏ฌ set (m) mu, eps_v (m/s) contacts avg. (max.) per timestep Cloth on sphere (aniso SL)
Sphere on cloth stack (1mm)
Sphere on cloth stack (5mm)
Sphere on cloth stack (10mm)
Cloth on rotating sphere (85K)
Cloth on rotating sphere (246K)
Funnel
Ribbon knot
Pulling cloth on needles
Table cloth pull (4m/s)
Twisting cylinder (with o ๏ฌ set) Cards shu ๏ฌ ling Braids
Hair clusters
Noodles
Sprinkles
Sand on cloth
All in
Fig. 24.
Simulation Statistics for IPC on a subset of our benchmark examples. Remaining garment examples are summarized in Figure 20. For each simulationwe report geometry, time step, strain-limit enforced (SL), thickness settings ( ห ๐ and offset), friction parameters ( ๐ and accuracy ๐ ๐ฃ ), number of contactsprocessed per time step and machine, as well as average timing and number of Newton iterations per time step solve. (*The tested anisotropic model hasstrain limits at . for both diagonal entries and . for the off-diagonal entry of the Green-Lagrangian strain [Clyde 2017].) strain-limiting solely with upper bounds and while we have not yetseen cases requiring lower bound limits these are easy to apply andso should be interesting for future exploration.While proven critical for C-IPC we have also shown that ACCDis an exceedingly simple replacement for pre-existing complex, andoften sensitive CCD modules. Here in preliminary testing we havefound that in application outside of C-IPC framework ACCD canalso provide significant speed-up. We look forward to further appli-cations of ACCD as an extremely easy-to-implement and efficientalternative CCD algorithm for tasks in both simulation and geome-try processing.In summary, independent from elasticity model and time stepsize, C-IPC guarantees out-of-the-box intersection-free simulationwith a strict satisfaction of strain limits (confirmed down to 0 . REFERENCES
David Baraff and Andrew Witkin. 1998. Large steps in cloth simulation. In
Proceedingsof the 25th annual conference on Computer graphics and interactive techniques . 43โ54.Jan Bender, Dan Koschier, Patrick Charrier, and Daniel Weber. 2014. Position-basedsimulation of continuous materials.
Computers & Graphics
44 (2014), 1โ10.Jan Bender, Matthias Mรผller, and Miles Macklin. 2015. Position-Based SimulationMethods in Computer Graphics.. In
Eurographics (tutorials) . 8.Jan Bender, Daniel Weber, and Raphael Diziol. 2013. Fast and stable cloth simulationbased on multi-resolution shape matching.
Computers & Graphics
37, 8 (2013),945โ954. Miklรณs Bergou, Max Wardetzky, Stephen Robinson, Basile Audoly, and Eitan Grinspun.2008. Discrete elastic rods. In
ACM SIGGRAPH 2008 papers . 1โ12.Kiran S Bhat, Christopher D Twigg, Jessica K Hodgins, Pradeep Khosla, Zoran Popovic,and Steven M Seitz. 2003. Estimating cloth simulation parameters from video. (2003).Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014.Projective Dynamics: Fusing Constraint Projections for Fast Simulation.
ACMTransactions on Graphics (TOG)
33, 4 (2014).Robert Bridson, Ronald Fedkiw, and John Anderson. 2002. Robust treatment of collisions,contact and friction for cloth animation. In
Proceedings of the 29th annual conferenceon Computer graphics and interactive techniques . 594โ603.Tyson Brochu, Essex Edwards, and Robert Bridson. 2012. Efficient geometrically exactcontinuous collision detection.
ACM Transactions on Graphics (TOG)
31, 4 (2012),1โ7.Rui PR Cardoso, Jeong Whan Yoon, Made Mahardika, S Choudhry, RJ Alves de Sousa,and RA Fontes Valente. 2008. Enhanced assumed strain (EAS) and assumed naturalstrain (ANS) methods for one-point quadrature solid-shell elements.
Int. J. Numer.Methods Eng.
75, 2 (2008), 156โ187.Juan J Casafranca, Gabriel Cirio, Alejandro Rodrรญguez, Eder Miguel, and Miguel AOtaduy. 2020. Mixing Yarns and Triangles in Cloth Simulation. In
Computer GraphicsForum , Vol. 39. Wiley Online Library, 101โ110.Jumyung Chang, Fang Da, Eitan Grinspun, and Christopher Batty. 2019. A UnifiedSimplicial Model for Mixed-Dimensional and Non-Manifold Deformable ElasticObjects.
Proceedings of the ACM on Computer Graphics and Interactive Techniques arXiv preprint arXiv:1911.05204 (2019).Hsiao-Yu Chen, Arnav Sastry, Wim M van Rees, and Etienne Vouga. 2018. Physicalsimulation of environmentally induced thin shell deformation.
ACM Transactionson Graphics (TOG)
37, 4 (2018), 1โ13.Ming Chen and Kai Tang. 2010. A fully geometric approach for developable clothdeformation simulation.
The visual computer
26, 6-8 (2010), 853โ863.Yanqing Chen, Timothy A Davis, William W Hager, and Sivasankaran Rajamanickam.2008. Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization andupdate/downdate.
ACM Transactions on Mathematical Software (TOMS)
35, 3 (2008),1โ14. odimensional Incremental Potential Contact โข 23
Byoungwon Choe, Min Gyu Choi, and Hyeong-Seok Ko. 2005. Simulating com-plex hair with robust collision handling. In
Proceedings of the 2005 ACM SIG-GRAPH/Eurographics symposium on Computer animation . 153โ160.David Clyde. 2017.
Numerical Subdivision Surfaces for Simulation and Data DrivenModeling of Woven Cloth . Ph.D. Dissertation. UCLA.David Clyde, Joseph Teran, and Rasmus Tamstorf. 2017. Modeling and data-driven parameter estimation for woven fabrics. In
Proceedings of the ACM SIG-GRAPH/Eurographics Symposium on Computer Animation . 1โ11.Peter A Cundall and Otto DL Strack. 1979. A discrete numerical model for granularassemblies. geotechnique
29, 1 (1979), 47โ65.Gilles Daviet. 2020. Simple and scalable frictional contacts for thin nodal objects.
ACMTransactions on Graphics (TOG)
39, 4 (2020), 61โ1.Gilles Daviet, Florence Bertails-Descoubes, and Laurence Boissieux. 2011. A hybrid iter-ative solver for robustly capturing coulomb friction in hair dynamics. In
Proceedingsof the 2011 SIGGRAPH Asia Conference . 1โ12.Crispin Deul, Tassilo Kugelstadt, Marcel Weiler, and Jan Bender. 2018. Direct position-based solver for stiff rods. In
Computer Graphics Forum , Vol. 37. Wiley Online Library,313โ324.Ellis Harold Dill. 1992. Kirchhoffโs theory of rods.
Archive for History of Exact Sciences
44, 1 (1992), 1โ23.Elliot English and Robert Bridson. 2008. Animating developable surfaces using non-conforming elements. In
ACM SIGGRAPH 2008 papers . 1โ5.Theodore F Gast, Craig Schroeder, Alexey Stomakhin, Chenfanfu Jiang, and Joseph MTeran. 2015. Optimization integrator for large time steps.
IEEE transactions onvisualization and computer graphics
21, 10 (2015), 1103โ1115.Rony Goldenthal, David Harmon, Raanan Fattal, Michel Bercovier, and Eitan Grinspun.2007. Efficient simulation of inextensible cloth. In
ACM SIGGRAPH 2007 papers .49โes.Eitan Grinspun, Anil N Hirani, Mathieu Desbrun, and Peter Schrรถder. 2003. Discreteshells. In
Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Com-puter animation . Citeseer, 62โ67.Gaรซl Guennebaud, Benoรฎt Jacob, et al. 2010. Eigen v3. http://eigen.tuxfamily.org.Qi Guo, Xuchen Han, Chuyuan Fu, Theodore Gast, Rasmus Tamstorf, and JosephTeran. 2018. A material point method for thin shells with frictional contact.
ACMTransactions on Graphics (TOG)
37, 4 (2018), 1โ15.David Harmon, Daniele Panozzo, Olga Sorkine, and Denis Zorin. 2011. Interference-aware geometric modeling.
ACM Transactions on Graphics (TOG)
30, 6 (2011), 1โ10.David Harmon, Etienne Vouga, Breannan Smith, Rasmus Tamstorf, and Eitan Grinspun.2009. Asynchronous contact mechanics. In
ACM SIGGRAPH 2009 papers . 1โ12.David Harmon, Etienne Vouga, Rasmus Tamstorf, and Eitan Grinspun. 2008. Robusttreatment of simultaneous collisions. In
ACM SIGGRAPH 2008 papers . 1โ4.R Hauptmann and K Schweizerhof. 1998. A systematic development of โsolid-shellโelement formulations for linear and non-linear analyses employing only dis-placement degrees of freedom.
Int. J. Numer. Methods Eng.
42, 1 (1998), 49โ69.F Hernandez, G Cirio, AG Perez, and MA Otaduy. 2013. Anisotropic strain limiting. In
Proc. of Congreso Espaรฑol de Informรกtica Grรกfica , Vol. 2.Chenfanfu Jiang, Theodore Gast, and Joseph Teran. 2017. Anisotropic elastoplasticityfor cloth, knit and hair frictional contact.
ACM Transactions on Graphics (TOG)
36, 4(2017), 1โ14.Yupeng Jiang, Minchen Li, Chenfanfu Jiang, and Fernando Alonso-Marroquin. 2020.A hybrid material-point spheropolygon-element method for solid and granularmaterial interaction.
Int. J. Numer. Methods Eng. (2020).Ning Jin, Wenlong Lu, Zhenglin Geng, and Ronald P Fedkiw. 2017. Inequality cloth. In
Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation .1โ10.Couro Kane, Jerrold E Marsden, Michael Ortiz, and Matthew West. 2000. Variationalintegrators and the Newmark algorithm for conservative and dissipative mechanicalsystems.
Int. J. Numer. Methods Eng.
49, 10 (2000), 1295โ1325.Danny M Kaufman and Dinesh K Pai. 2012. Geometric numerical integration ofinequality constrained, nonsmooth Hamiltonian systems.
SIAM J. on Sci. Comp.
ACM Trans. on Graphics (SIGGRAPH 2014) (2014).Ted Kim. 2020. A finite element formulation of baraff-witkin cloth. In
Symposium onComputer Animation .Tassilo Kugelstadt and Elmar Schรถmer. 2016. Position and orientation based Cosseratrods.. In
Symposium on Computer Animation . 169โ178.Cheng Li, Min Tang, Ruofeng Tong, Ming Cai, Jieyi Zhao, and Dinesh Manocha. 2020b.P-cloth: interactive complex cloth simulation on multi-GPU systems using dynamicmatrix assembly and pipelined implicit integrators.
ACM Transactions on Graphics(TOG)
39, 6 (2020), 1โ15.Jie Li, Gilles Daviet, Rahul Narain, Florence Bertails-Descoubes, Matthew Overby,George E Brown, and Laurence Boissieux. 2018a. An implicit frictional contactsolver for adaptive cloth simulation.
ACM Transactions on Graphics (TOG)
37, 4(2018), 1โ15. Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, DanielePanozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020a. Incremental PotentialContact: Intersection- and Inversion-free Large Deformation Dynamics.
ACMTransactions on Graphics
39, 4 (2020).Minchen Li, Ming Gao, Timothy Langlois, Chenfanfu Jiang, and Danny M Kaufman.2019. Decomposed optimization time integrator for large-step elastodynamics.
ACMTransactions on Graphics (TOG)
38, 4 (2019), 1โ10.Minchen Li, Alla Sheffer, Eitan Grinspun, and Nicholas Vining. 2018b. FoldSketch:Enriching Garments with Physically Reproducible Folds.
ACM Transaction onGraphics
37, 4 (2018). https://doi.org/10.1145/3197517.3201310Junbang Liang, Ming Lin, and Vladlen Koltun. 2019. Differentiable cloth simulation forinverse problems. In
Advances in Neural Information Processing Systems . 772โ781.Tiantian Liu, Sofien Bouaziz, and Ladislav Kavan. 2017. Quasi-newton methods forreal-time simulation of hyperelastic materials.
ACM Transactions on Graphics (TOG)
36, 3 (2017), 1โ16.Libin Lu, Matthew J Morse, Abtin Rahimian, Georg Stadler, and Denis Zorin. 2019.Scalable simulation of realistic volume fraction red blood cell flows through vascu-lar networks. In
Proceedings of the International Conference for High PerformanceComputing, Networking, Storage and Analysis . 1โ30.Mickaรซl Ly, Jean Jouve, Laurence Boissieux, and Florence Bertails-Descoubes. 2020.Projective Dynamics with Dry Frictional Contact.
ACM Transactions on Graphics(TOG)
39, 4 (2020).Miles Macklin, Matthias Mรผller, and Nuttapong Chentanez. 2016. XPBD: position-basedsimulation of compliant constrained dynamics. In
Proceedings of the 9th InternationalConference on Motion in Games . 49โ54.Sebastian Martin, Peter Kaufmann, Mario Botsch, Eitan Grinspun, and Markus Gross.2010. Unified simulation of elastic rods, shells, and solids.
ACM Transactions onGraphics (TOG)
29, 4 (2010), 1โ10.Aleka McAdams, Andrew Selle, Kelly Ward, Eftychios Sifakis, and Joseph Teran. 2009.Detail preserving continuum simulation of straight hair.
ACM Transactions onGraphics (TOG)
28, 3 (2009), 1โ6.Eder Miguel, Derek Bradley, Bernhard Thomaszewski, Bernd Bickel, Wojciech Matusik,Miguel A Otaduy, and Steve Marschner. 2012. Data-driven estimation of clothsimulation models. In
Computer Graphics Forum , Vol. 31. Wiley Online Library,519โ528.Matthias Mรผller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007. Positionbased dynamics.
Journal of Visual Communication and Image Representation
18, 2(2007), 109โ118.Matthias Mรผller, Tae-Yong Kim, and Nuttapong Chentanez. 2012. Fast simulation ofinextensible hair and fur.
VRIPHYS
12 (2012), 39โ44.Rahul Narain, Armin Samii, and James F Oโbrien. 2012. Adaptive anisotropic remeshingfor cloth simulation.
ACM transactions on graphics (TOG)
31, 6 (2012), 1โ10.Miguel A Otaduy, Rasmus Tamstorf, Denis Steinemann, and Markus Gross. 2009. Im-plicit contact handling for deformable objects. In
Computer Graphics Forum , Vol. 28.Wiley Online Library, 559โ568.Matthew Overby, George E Brown, Jie Li, and Rahul Narain. 2017. ADMM โ ProjectiveDynamics: Fast Simulation of Hyperelastic Models with Dynamic Constraints.
IEEETransactions on Visualization and Computer Graphics
23, 10 (2017), 2222โ2234.ลฝeljko Penava, Diana ล imiฤ-Penava, and ลฝ Knezic. 2014. Determination of the elasticconstants of plain woven fabrics by a tensile test in various directions.
Fibres &Textiles in Eastern Europe (2014).Xavier Provot. 1997. Collision and self-collision handling in cloth model dedicated todesign garments. In
Computer Animation and Simulationโ97 . Springer, 177โ189.Xavier Provot et al. 1995. Deformation constraints in a mass-spring model to describerigid cloth behaviour. In
Graphics interface . Canadian Information Processing Society,147โ147.Alessio Quaglino. 2012.
Membrane locking in discrete shell theories . Ph.D. Dissertation.Georg-August-Universitรคt Gรถttingen.Stefanie Reese. 2007. A large deformation solid-shell concept based on reduced in-tegration with hourglass stabilization.
Int. J. Numer. Methods Eng.
69, 8 (2007),1671โ1716.Nikolas Schmitt, Martin Knuth, Jan Bender, and Arjan Kuijper. 2013. Multilevel ClothSimulation using GPU Surface Sampling.
VRIPHYS
13 (2013), 1โ10.Andrew Selle, Michael Lentine, and Ronald Fedkiw. 2008. A mass spring model for hairsimulation. In
ACM SIGGRAPH 2008 papers . 1โ11.Georg Sperl, Rahul Narain, and Chris Wojtan. 2020. Homogenized yarn-level cloth.
ACM Transactions on Graphics (TOG)
39, 4 (2020), 48โ1.Alexey Stomakhin, Russell Howes, Craig Schroeder, and Joseph M Teran. 2012. En-ergetically consistent invertible elasticity. In
Proceedings of the 11th ACM SIG-GRAPH/Eurographics conference on Computer Animation . 25โ32.Rasmus Tamstorf and Eitan Grinspun. 2013. Discrete bending forces and their Jacobians.
Graphical models
75, 6 (2013), 362โ370.Rasmus Tamstorf, Toby Jones, and Stephen F McCormick. 2015. Smoothed aggregationmultigrid for cloth simulation.
ACM Transactions on Graphics (TOG)
34, 6 (2015),1โ13.
Min Tang, Dinesh Manocha, Sung-Eui Yoon, Peng Du, Jae-Pil Heo, and Ruo-Feng Tong.2011. VolCCD: Fast continuous collision culling between deforming volume meshes.
ACM Transactions on Graphics (TOG)
30, 5 (2011), 1โ15.Min Tang, Ruofeng Tong, Rahul Narain, Chang Meng, and Dinesh Manocha. 2013. AGPU-based streaming algorithm for high-resolution cloth simulation. In
ComputerGraphics Forum , Vol. 32. Wiley Online Library, 21โ30.Min Tang, Ruofeng Tong, Zhendong Wang, and Dinesh Manocha. 2014. Fast and exactcontinuous collision detection with bernstein sign classification.
ACM Transactionson Graphics (TOG)
33, 6 (2014), 1โ8.Min Tang, Huamin Wang, Le Tang, Ruofeng Tong, and Dinesh Manocha. 2016. CAMA:Contact-aware matrix assembly with unified collision handling for GPU-based clothsimulation. In
Computer Graphics Forum , Vol. 35. Wiley Online Library, 511โ521.Min Tang, Tongtong Wang, Zhongyuan Liu, Ruofeng Tong, and Dinesh Manocha. 2018.I-cloth: incremental collision handling for GPU-based interactive cloth simulation.
ACM Transactions on Graphics (TOG)
37, 6 (2018), 1โ10.Demetri Terzopoulos, John Platt, Alan Barr, and Kurt Fleischer. 1987. Elastically de-formable models. In
Proceedings of the 14th annual conference on Computer graphicsand interactive techniques . 205โ214.Bernhard Thomaszewski, Simon Pabst, and Wolfgang Strasser. 2009. Continuum-basedstrain limiting. In
Computer Graphics Forum , Vol. 28. Wiley Online Library, 569โ576.Nobuyuki Umetani, Danny M Kaufman, Takeo Igarashi, and Eitan Grinspun. 2011.Sensitive couture for interactive garment modeling and editing.
ACM Trans. Graph.
30, 4 (2011), 90.Pascal Volino and N Magnenat Thalmann. 2000. Implementing fast cloth simulationwith collision response. In
Proceedings Computer Graphics International 2000 . IEEE,257โ266.Huamin Wang, James OโBrien, and Ravi Ramamoorthi. 2010. Multi-resolution isotropicstrain limiting.
ACM Transactions on Graphics (TOG)
29, 6 (2010), 1โ10.Huamin Wang, James F OโBrien, and Ravi Ramamoorthi. 2011. Data-driven elasticmodels for cloth: modeling and measurement.
ACM transactions on graphics (TOG)
30, 4 (2011), 1โ12.Xinlei Wang, Minchen Li, Yu Fang, Xinxin Zhang, Ming Gao, Min Tang, Danny MKaufman, and Chenfanfu Jiang. 2020. Hierarchical optimization time integrationfor cfl-rate mpm stepping.
ACM Transactions on Graphics (TOG)
39, 3 (2020), 1โ16.Zhendong Wang, Min Tang, Ruofeng Tong, and Dinesh Manocha. 2015. TightCCD:Efficient and robust continuous collision detection using tight error bounds. In
Computer Graphics Forum , Vol. 34. Wiley Online Library, 289โ298.Zhendong Wang, Tongtong Wang, Min Tang, and Ruofeng Tong. 2016. Efficient androbust strain limiting and treatment of simultaneous collisions with semidefiniteprogramming.
Computational Visual Media
2, 2 (2016), 119โ130.Zhendong Wang, Longhua Wu, Marco Fratarcangeli, Min Tang, and Huamin Wang.2018. Parallel Multigrid for Nonlinear Cloth Simulation.
Computer Graphics Forum (2018).Kelly Ward and Ming C Lin. 2003. Adaptive grouping and subdivision for simulatinghair dynamics. In
IEEE, 234โ243.Nicholas J Weidner, Kyle Piddington, David IW Levin, and Shinjiro Sueda. 2018. Eulerian-on-lagrangian cloth simulation.
ACM Transactions on Graphics (TOG)
37, 4 (2018),1โ11.Zangyueyang Xian, Xin Tong, and Tiantian Liu. 2019. A Scalable Galerkin MultigridMethod for Real-time Simulation of Deformable Objects.
ACM Transactions onGraphics (TOG)
38, 6 (2019).Hongyi Xu, Funshing Sin, Yufeng Zhu, and Jernej Barbiฤ. 2015. Nonlinear materialdesign using principal stretches.
ACM Transactions on Graphics (TOG)
34, 4 (2015),1โ11.H.T.Y. Yang, S. Saigal, A. Masud, and R.K. Kapania. 2000. A survey of recent shell finiteelements.
Int. J. Numer. Methods Eng.
47 (2000).Yonghao Yue, Breannan Smith, Peter Yichen Chen, Maytee Chantharayukhonthorn,Ken Kamrin, and Eitan Grinspun. 2018. Hybrid grains: adaptive coupling of discreteand continuum simulations of granular media.
ACM Transactions on Graphics (TOG)
37, 6 (2018), 1โ19.
APPENDIX A - ISOTROPIC STRAIN LIMITING DERIVATIVEDERIVATIONS
Derivatives of singular value decomposition (SVD) will be needed: ๐๐ ๐๐ ๐๐น = ๐ข ๐ ๐ฃ ๐๐ , ๐ ๐ ๐๐ ๐๐น = ๐ ( ๐ข ๐ ๐ฃ ๐๐ ) ๐๐น = ๐๐ข ๐ ๐๐น ๐ฃ ๐๐ + ๐ข ๐ ๐๐ฃ ๐ ๐๐น ๐ (13)where ๐๐ข ๐ ๐๐น and ๐๐ฃ ๐ ๐๐น can be found in Xu et al. [2015]. Then if we denote the barrier ๐ ๐ก,๐ for ๐ -th singular value of triangle ๐ก , we have ๐๐ ๐ก,๐ ๐๐ฅ = ๐๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐๐ ๐ก,๐๐ ๐๐น ๐ก ๐๐น ๐ก ๐๐ฅ๐ ๐ ๐ก,๐ ๐๐ฅ = ๐๐ ๐ก,๐๐ ๐๐ฅ ๐ ๐ ๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐๐ ๐ก,๐๐ ๐๐ฅ + ๐๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐ ๐ ๐ก,๐๐ ๐๐ฅ (14)where ๐๐น ๐ก ๐๐ฅ can be found in anisotropic strain limiting implementa-tion detail section.For the barrier derivatives, we have ๐๐ ๐ก,๐ ๐๐ ๐ก,๐๐ = ๐ ๐ ( ๐ โ ห ๐ ) ( ( ห ๐ โ ๐ ๐ก,๐๐ ) ln ( ๐ โ ๐ ๐ก,๐๐ ๐ โ ห ๐ ) + ( ห ๐ โ ๐ ๐ก,๐๐ ) ๐ โ ๐ ๐ก,๐๐ ) ๐ ๐ ๐ก,๐ ๐๐ ๐ก,๐๐ = ๐ ๐ ( ๐ โ ห ๐ ) (โ ( ๐ โ ๐ ๐ก,๐๐ ๐ โ ห ๐ ) โ ๐ โ ๐ ๐ก,๐๐ ๐ โ ๐ ๐ก,๐๐ + ( ห ๐ โ ๐ ๐ก,๐๐ ๐ โ ๐ ๐ก,๐๐ ) ) , (15)then together with ๐๐ ๐ก,๐๐ ๐๐ฅ = ๐๐ ๐ก,๐๐ ๐๐น ๐ก ๐๐น ๐ก ๐๐ฅ and ๐ ๐ ๐ก,๐๐ ๐๐ฅ = ๐๐น ๐ก ๐๐ฅ ๐ ๐ ๐ ๐ก,๐๐ ๐๐น ๐ก ๐๐น ๐ก ๐๐ฅ (16)we can efficiently compute ๐๐ ๐ก ๐๐ฅ = โ๏ธ ๐ ( ๐๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐๐ ๐ก,๐๐ ๐๐น ๐ก ) ๐๐น ๐ก ๐๐ฅ๐ ๐ ๐ก,๐ ๐๐ฅ = ๐๐น ๐ก ๐๐ฅ ๐ โ๏ธ ๐ ( ๐ ๐ ๐ก,๐ ๐๐น ๐ก ) ๐๐น ๐ก ๐๐ฅ (17)where ๐ ๐ ๐ก,๐ ๐๐น ๐ก = ๐๐ ๐ก,๐๐ ๐๐น ๐ก ๐ ๐ ๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐๐ ๐ก,๐๐ ๐๐น ๐ก + ๐๐ ๐ก,๐ ๐๐ ๐ก,๐๐ ๐ ๐ ๐ก,๐๐ ๐๐น ๐ก .For SPD projection we can process the 6x6 matrix (cid:205) ๐ ( ๐ ๐ ๐ก,๐ ๐๐น ๐ก ) . APPENDIX B - ANISOTROPIC STRAIN LIMITINGDERIVATIONSB.1 Quadratically Approximating Data-Driven Model
The 1st-order derivatives of ๐ can be computed as ๐๐๐ ห ๐ธ = ๐ ๐ โฒ ( ห ๐ธ ) ห ๐ธ + ๐ ๐ โฒ ( ห ๐ธ ห ๐ธ ) ห ๐ธ ๐๐๐ ห ๐ธ = ๐ ๐ โฒ ( ห ๐ธ ) ห ๐ธ + ๐ ๐ โฒ ( ห ๐ธ ห ๐ธ ) ห ๐ธ ๐๐๐ ห ๐ธ = ๐บ ๐ โฒ ( ห ๐ธ ) ห ๐ธ . (18)When ห ๐ธ =
0, element is at rest shape, and all gradients above areequal to 0. odimensional Incremental Potential Contact โข 25
For the 2nd-order derivatives, we have ๐ ๐๐ ห ๐ธ = ๐ ๐ โฒโฒ ( ห ๐ธ ) ห ๐ธ + ๐ ๐ โฒ ( ห ๐ธ ) + ๐ ๐ โฒโฒ ( ห ๐ธ ห ๐ธ ) ห ๐ธ ๐ ๐๐ ห ๐ธ ๐ ห ๐ธ = ๐ ๐๐ ห ๐ธ ๐ ห ๐ธ = ๐ ๐ โฒโฒ ( ห ๐ธ ห ๐ธ ) ห ๐ธ ห ๐ธ + ๐ ๐ โฒ ( ห ๐ธ ห ๐ธ ) ๐ ๐๐ ห ๐ธ = ๐ ๐ โฒโฒ ( ห ๐ธ ) ห ๐ธ + ๐ ๐ โฒ ( ห ๐ธ ) + ๐ ๐ โฒโฒ ( ห ๐ธ ห ๐ธ ) ห ๐ธ ๐ ๐๐ ห ๐ธ = ๐บ ๐ โฒโฒ ( ห ๐ธ ) ห ๐ธ + ๐บ ๐ โฒ ( ห ๐ธ ) , (19)and other terms are all 0. When ห ๐ธ =
0, we have ๐ ๐๐ ห ๐ธ = ๐ , ๐ ๐๐ ห ๐ธ = ๐ , ๐ ๐๐ ห ๐ธ = ๐บ ๐ ๐๐ ห ๐ธ ๐ ห ๐ธ = ๐ ๐๐ ห ๐ธ ๐ ห ๐ธ = ๐ , (20)all constants. Thus to quadratically approximate the data-drivenmodel at ห ๐ธ = B.2 Implementation Details
For simplicity, if we align the cloth weft and warp directions to thecoordinate axes, we can easily compute ห ๐ธ from the 3x2 deformationgradient ๐น = [ ๐ฅ โ ๐ฅ , ๐ฅ โ ๐ฅ ] [ ๐ โฒ โ ๐ โฒ , ๐ โฒ โ ๐ โฒ ] โ asห ๐ธ ๐ ๐ = ( ๐น ๐๐ ๐น ๐๐ โ ๐ฟ ๐ ๐ ) , so we have ๐ ห ๐ธ ๐ ๐ ๐๐น ๐๐ = ( ๐ฟ ๐๐ ๐น ๐ ๐ + ๐น ๐๐ ๐ฟ ๐๐ ) ๐ ห ๐ธ ๐ ๐ ๐๐น ๐๐ ๐๐น ๐๐ = ๐ฟ ๐๐ ( ๐ฟ ๐๐ ๐ฟ ๐๐ + ๐ฟ ๐๐ ๐ฟ ๐๐ ) . (21)Then together with ๐๐ above and ๐๐น / ๐๐ฅ = { (cid:20) โ ๐ผ โ ๐ผ (cid:21) (cid:20) ๐ผ (cid:21) (cid:20) ๐ผ (cid:21) } : ๐ต where ๐ต = [ ๐ โฒ โ ๐ โฒ , ๐ โฒ โ ๐ โฒ ] โ so that ๐๐น๐๐ฅ = ๏ฃฎ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฐ โ( ๐ต + ๐ต ) โ โโ( ๐ต + ๐ต ) โ โ ๏ฃน๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃป , ๏ฃฎ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฐ ๐ต โ โ ๐ต โ โ ๏ฃน๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃป , ๏ฃฎ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฏ๏ฃฐ ๐ต โ โ ๐ต โ โ ๏ฃน๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃบ๏ฃป (22)where โ means identical to the top-left entry, we are able to imple-ment the constitutive model via ๐๐๐๐ฅ = ๐๐๐ ห ๐ธ ๐ ห ๐ธ๐๐น ๐๐น๐๐ฅ and ๐ ๐๐๐ฅ = ( ๐๐น๐๐ฅ ) ๐ ๐ ๐๐๐น ๐๐น๐๐ฅ (23)where ๐ ๐๐๐น = ( ๐ ห ๐ธ๐๐น ) ๐ ๐ ๐๐ ห ๐ธ ๐ ห ๐ธ๐๐น + ๐๐๐ ห ๐ธ ๐ ห ๐ธ๐๐น APPENDIX C - ROBUST AND EFFICIENTIMPLEMENTATION OF IPC FRICTION
To compute friction potential, ๐ is well defined without any divisionby 0.For friction gradient/force, a robust implementation requires al-gebraically deriving ๐ ( | | ๐ข ๐ | |)| | ๐ข ๐ | | , which does not contain division by 0by construction of ๐ .For friction Hessian, the inner 2 ร ๐ โฒ (|| ๐ข ๐ ||)|| ๐ข ๐ || โ ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ข ๐ ๐ข ๐๐ + ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ผ , one needs to first algebraically derive ๐ โฒ ( | | ๐ข ๐ | |) | | ๐ข ๐ | |โ ๐ ( | | ๐ข ๐ | |)| | ๐ข ๐ | | whichdoes not contain division by 0 by construction, and note that ๐ข ๐ ๐ข ๐๐ | | ๐ข ๐ | | isalways bounded and it equals 0 when || ๐ข ๐ || =
0, in which case since ๐ ( | | ๐ข ๐ | |)| | ๐ข ๐ | | > ๐ข ๐ = ( ๐ฅ, ๐ฆ ) , ๐ โฒ (|| ๐ข ๐ ||)|| ๐ข ๐ || โ ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ข ๐ ๐ข ๐๐ + ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ผ = ๐ โฒ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ข ๐ ๐ข ๐๐ + โ ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || (cid:20) ๐ฅ ๐ฅ๐ฆ๐ฅ๐ฆ ๐ฆ (cid:21) + ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || (cid:20) ๐ฅ + ๐ฆ ๐ฅ + ๐ฆ (cid:21) = ๐ โฒ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ข ๐ ๐ข ๐๐ + ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || (cid:20) ๐ฆ โ ๐ฅ๐ฆ โ ๐ฅ๐ฆ ๐ฅ (cid:21) = ๐ โฒ (|| ๐ข ๐ ||)|| ๐ข ๐ || ๐ข ๐ ๐ข ๐๐ + ๐ (|| ๐ข ๐ ||)|| ๐ข ๐ || ยฏ ๐ข ๐ ยฏ ๐ข ๐๐ (24)where ยฏ ๐ข ๐ = (โ ๐ฆ, ๐ฅ ) . So for || ๐ข ๐ || โฅ ๐ ๐ฃ โ when we have ๐ โฒ (ยท) = ๐ (ยท) =
1, the matrix is always SPD (must compute using ยฏ ๐ข ๐ ,the original formula can still be not SPD due to numerical error ofsubtraction!); otherwise we project the 2x2 matrix. APPENDIX D - ROD BENDING MODULUS
In discrete rods model [Bergou et al. 2008], the bending energy isdirectly integrated over the length of the rod as ๐ธ bend ( ๐ฅ ) = ๐ โ๏ธ ๐ = ๐ผ ( ๐ ๐ ๐ ) ยฏ ๐ ๐ (25)Here ๐ ๐ ๐ is the curvature binormal at a vertex, ยฏ ๐ ๐ = | ยฏ ๐ ๐ โ | + | ยฏ ๐ ๐ | where | ยฏ ๐ ๐ | is the rest length of edge ๐ , and ๐ผ is an adjustable pa-rameter to control the stiffness of bending. As we know when arod becomes thicker, it is harder to bend, which means we need tomanually increase ๐ผ accordingly and this makes seting up examplesinconvenient.According to Kirchoff rod theory [Dill 1992], the bending energyof rod at a unit length is defined as12 ๐ธ๐ ๐ ๐ where ๐ธ is the bending modulus in ๐๐ unit, ๐ is the radius of therod, and ๐ = ๐ ๐ ยฏ ๐ / is the point-wise curvature [Bergou et al. 2008]. Integrating this quantity over rod length then gives us ๐ธ bend ( ๐ฅ ) = โซ ฮฉ ๐ธ๐ ๐ ๐ ๐๐ โ โ๏ธ ๐ ยฏ ๐ ๐ ๐ธ๐ ๐ ( ๐ ๐ ๐ ยฏ ๐ ๐ / ) = โ๏ธ ๐ ๐ธ๐ ๐ ( ๐ ๐ ๐ ) ยฏ ๐ ๐ (26) Combining with Equation 25 we finally know that the ๐ผ in Bergouet al. [2008] is associated with the rodโs radius as ๐ผ = ๐ธ๐ ๐ ๐ธ can be automatically reflected in the bending stiffness. Thisenables us to easily setup rod examples starting from using real ma-terial thickness and Youngโs modulus of e.g. hairs, iron, etc, withouthaving to tune bending stiffness ๐ผ๐ผ