# DiffCloth: Differentiable Cloth Simulation with Dry Frictional Contact

Yifei Li  MIT CSAILCambridgeMAUSA<EMAIL_ADDRESS>,  Tao Du  MIT
CSAILCambridgeMAUSA<EMAIL_ADDRESS>,  Kui Wu  Tencent Lightspeed &
Quantum StudiosUSA<EMAIL_ADDRESS>,  Jie Xu  MIT CSAILCambridgeMAUSA
<EMAIL_ADDRESS>and  Wojciech Matusik  MIT CSAILCambridgeMAUSA
<EMAIL_ADDRESS>
###### Abstract.

Cloth simulation has wide applications including computer animation, garment
design, and robot-assisted dressing. In this work, we present a differentiable
cloth simulator whose additional gradient information facilitates cloth-
related applications. Our differentiable simulator extends the state-of-the-
art cloth simulator (Ly et al., 2020) based on Projective Dynamics and with
dry frictional contact governed by the Signorini-Coulomb law. We derive
gradients with contact in this forward simulation framework and speed up the
computation with Jacobi iteration inspired by previous differentiable
simulation work (Du et al., 2021). To our best knowledge, we present the first
differentiable cloth simulator with the Coulomb law of friction. We
demonstrate the efficacy of our simulator in various applications, including
system identification, manipulation, inverse design, and a real-to-sim task.
Many of our applications have not been demonstrated in previous differentiable
cloth simulators. The gradient information from our simulator enables
efficient gradient-based task solvers from which we observe a substantial
speedup over standard gradient-free methods.

Projective Dynamics, differentiable simulation, cloth simulation

††submissionid: 241††ccs: Computing methodologies Physical simulation Figure
1. We present a differentiable cloth simulator with Coulomb frictional contact
and demonstrate its efficacy in various applications: estimating cloth
material parameters and the frictional coefficient (left two), putting on a
hat and wearing a sock with robotic manipulators (middle two), designing and
controlling the dynamical motion of a dress (second to right), and
reconstructing a digital twin of real-world cloth from its motion-captured
data (rightmost).

##  1\. Introduction

Clothing is ubiquitous in our daily lives. With the widespread appearance of
clothing in the fashion industry, film industry, computer animation, and video
games, simulating cloth has been an active research topic for more than two
decades. Today, research advancement in cloth simulation has unlocked various
applications such as virtual try-on (Guan et al., 2012), garment design
(Umetani et al., 2011; Wang, 2018; Bartle et al., 2016; Montes et al., 2020),
fold design (Li et al., 2018b), garment grading (Brouet et al., 2012),
sagging-free inversion (Ly et al., 2018), and robot-assisted dressing (Clegg
et al., 2020; Clegg et al., 2018).

Inspired by the recent development of differentiable physics simulation and
its success in rigid-body systems (de Avila Belbute-Peres et al., 2018a;
Geilinger et al., 2020), fluidic systems (Hu et al., 2020; Du et al., 2020),
and deformable-body systems (Hahn et al., 2019; Geilinger et al., 2020), we
argue in this paper that a large number of cloth-related applications would
also benefit from a high-quality differentiable cloth simulator. The critical
ingredient in previous differentiable simulators is their ability to compute
gradients by backpropagating any differentiable performance metrics through
simulation. Such additional gradient information unlocks gradient-based
continuous optimization methods, which often bring a substantial speedup
compared with their gradient-free counterparts in downstream applications.

Compared with the recent active development of differentiable simulators for
rigid-body and soft-body dynamics, research work about differentiable cloth
simulation is still relatively sparse. Indeed, cloth simulation introduces
unique challenges from its co-dimensional dynamics and, in particular, rich
contact events. While many differentiable simulators have provided solutions
to deriving gradients with contact models of varying complexity, their
techniques typically do not expect contacts to be as frequent as in cloth
simulation. Deriving gradients with frequent contacts and self-collisions in
cloth simulation is not fully resolved even in the state-of-the-art
differentiable cloth simulators (Liang et al., 2019), and our work attempts to
fill this gap.

In this work, we present a differentiable cloth simulator with extra care of
its contact model. We choose to base our simulator on the state-of-the-art
cloth simulator described by Ly et al. (2020) and employs both its Projective
Dynamics (PD) simulation method and its dry frictional contacts described by
the Signorini-Coulomb law. Therefore, our differentiable cloth simulator
inherits both the speedup and robustness from Projective Dynamics and the
physical accuracy from the dry frictional contact model. Unlike previous
papers relying on automatic differentiation tools to derive gradients, we
modify the Jacobi solver in Du et al. (2021) to accommodate the dry frictional
contact model. We show that the modified Jacobi solver leads to a substantial
speedup in the gradient computation.

We demonstrate the efficacy of our simulator in various applications,
including system identification of material parameters from the cloth data
captured in the real world, inverse design of garments for computer animation,
and planning of robotic manipulators in robot-assisted dressing. Many of these
applications would either be impossible or require a much longer time with
existing methods. With the extra gradient information from our differentiable
simulator, we unlock gradient-based optimizers to solve these problems with a
much higher sample efficiency than traditional gradient-free methods such as
evolutionary strategies (ES) or reinforcement learning (RL).

To summarize, our paper contributes the following:

  1. (1)

We develop a differentiable cloth simulator that supports gradients with dry
frictional contacts based on a state-of-the-art cloth simulation method (Ly et
al., 2020);

  2. (2)

We present a Jacobi solver modified from existing work (Du et al., 2021) and
proves empirically its speedup in gradient computation with rich contacts in
cloth simulation;

  3. (3)

we evaluate the efficacy of our simulator in various applications, including
system identification, trajectory optimization, assisted dressing, and real-
to-sim transfer.

##  2\. Related Work

Our work is closely related to cloth simulation and its applications in
computer graphics. It is also relevant to the more recent differentiable
physics methods developed in the graphics and machine learning communities.

#### Cloth Simulation

Physics-based cloth simulation has been a popular topic in the graphic
community for decades since clothing has been widely used in our daily life.
The implicit Euler integration is used to simulate cloth robustly with large
timesteps (Terzopoulos et al., 1987; Baraff and Witkin, 1998) while
introducing excessive numerical damping. To alleviate this problem, implicit
and explicit methods (IMEX) (Bridson et al., 2003; Stern and Grinspun, 2009)
explicitly integrate elastic forces and implicitly integrate damping forces.
Researchers have also introduced other variational integrators, e.g., BDF2
(Choi and Ko, 2002) and Symplectic integrator (Stern and Desbrun, 2006), to
conserve the total energy of the system.

For cloth modeled as a mass-spring system, Liu et al. (2013) treat the
implicit Euler integration as an energy minimization problem. The global
linear system then remains constant on run-time, and each spring constraint
can be solved separately in a local step. That global/local solver idea is
generalized as Projective Dynamics (PD) (Bouaziz et al., 2014), which supports
material models whose elastic energy has a specific quadratic form. Projective
Dynamic is further extended to support more general materials (Overby et al.,
2017; Liu et al., 2017). Though a local step in PD can be processed in
parallel, the global step still needs to maintain a large pre-factorized
sparse matrix and do back substitutions in each step. Another relevant idea is
Position-Based Dynamics (PBD) (Müller et al., 2007; Macklin et al., 2016),
which iteratively projects each constraint in a non-linear Gauss-Seidel-like
fashion, leading to highly parallelizable computation on GPUs.

Both PD and PBD have led to a few follow-up works on more advanced speedup
techniques applicable to cloth simulation. Wang (2015) and a follow-up work
(Wang and Yang, 2016) propose a Chebyshev acceleration technique that can be
applied to both PD and PBD to speed up convergence. Fratarcangeli et al.
(2016) also introduce a parallel randomized Gauss-Seidel method that re-
organizes the unknowns of the sparse linear system into a few independent
blocks, which can be solved in parallel with a single Gauss-Seidel step.
Recently, the computation of integration is further accelerated by the
multigrid method, e.g., geometric multigrid scheme (Wang et al., 2018),
algebraic multigrid scheme (Tamstorf et al., 2015), and Galerkin multigrid
scheme (Xian et al., 2019).

Lastly, our work is also relevant to previous attempts to matching cloth
simulation with real fabric behaviors (Miguel et al., 2012; Wang et al., 2011;
Miguel et al., 2013; Clyde et al., 2017), which is typically done by fitting
constitutive material models with real material properties and can benefit
from extra gradient information from a differentiable simulator (Hahn et al.,
2019).

#### Cloth Contact and Friction

Contact and friction are key ingredients in modern cloth simulation. Provot
(1997) proposes a method called “impact zones” to collect the nodes involved
in multiple collisions into impact zones, which are treated as rigid bodies.
Harmon et al. (2008) improve the failsafe of impact zones by allowing some
sliding motion of the incriminated vertices. Bridson et al. (2002) introduce a
hybrid method to combine the idea of applying repulsion impulses and friction
and impact zones to handle cloth collision robustly, which is still widely
used in modern cloth simulators (Narain et al., 2012; Li et al., 2020b). To
model friction fully implicitly, Brown et al. (2018) treats friction as an
additional dissipative term in optimization. Regarding contact handling,
repulsive forces, or penalty methods, have been widely used (Bouaziz et al.,
2014; Wang, 2015; Geilinger et al., 2020; Macklin et al., 2020) since it is
easy to implement. However, these methods often need high stiffness in the
penalty energy, leading to disturbing jitter artifacts that demand careful
tuning. Recently, Li et al. (2020a); Li et al. (2021) define a smooth
dissipative potential for friction using barrier methods. Their method can
guarantee interpenetration-free states after each timestep without the need
for parameter tuning.

Unlike penalty-based methods, constraint-based collision handling methods
formulate contact as constraints in the physics system to provide more
physics-based simulation results. Otaduy et al. (2009) first formulate cloth
contacts as a sparse Linear Complementarity Problem (LCP). Their key idea is
to interleave frictional contact iterations with normal contact iterations.
Instead of using a pyramidal approximation of the Coulomb friction cone
(Otaduy et al., 2009), Li et al. (2018a) rely upon the exact Coulomb friction
cone and adaptively refine nodes to ensure an accurate treatment of frictional
contact. Although constraint-based methods typically ensure the physics-based
constraints characterizing contact and friction are satisfied after time
integration, their computation cost is typically much more expensive than a
penalty-based method. Recently, Ly et al. (2020) propose an efficient
algorithm to incorporate frictional contact into Projective Dynamics so that
non-penetration and Coulomb constraints are satisfied simultaneously in a
semi-implicit way.

#### Differentiable Physics

Despite the recent advances in differentiable physics, e.g., differentiable
simulators in rigid-body dynamics (Popović et al., 2003; de Avila Belbute-
Peres et al., 2018b; Toussaint et al., 2018; Degrave et al., 2019; Qiao et
al., 2020), soft-body dynamics (Hu et al., 2019, 2020; Hahn et al., 2019;
Geilinger et al., 2020; Du et al., 2021), and fluid dynamics (Treuille et al.,
2003; McNamara et al., 2004; Wojtan et al., 2006; Schenck and Fox, 2018; Holl
et al., 2020), differentiable cloth simulation still lacks a good solution.
One natural idea is to use particle-based strategies that approximate a nodal
system of cloth with graph neural networks (Li et al., 2019; Sanchez-Gonzalez
et al., 2020). Although the neural networks are naturally differentiable, the
physical accuracy is hard to guarantee.

To accurately predict the behavior of real-world objects, recent papers Hahn
et al. (2019); Geilinger et al. (2020) present differentiable soft-body
systems and mass-spring systems with implicit Euler time integration and
penalty-based contacts. Although a few recent papers (Macklin et al., 2020; Li
et al., 2020a; Li et al., 2021) have introduced differentiable contact
handling methods, none of them show a clothing example using the
differentiability besides demonstrating the differentiability of their methods
in theory. Liang et al. (2019) are the first to introduce a fully functional
differentiable cloth simulation with contact, friction, and self-collision.
Instead of constructing a static LCP problem, they develop a Quadratic-
Programming (QP) problem to minimize the change between the collision-free
state and the original mesh state. In this paper, we build upon the
differentiable PD framework (Du et al., 2021) to simulate cloth dynamics with
physics-based frictional contact (Ly et al., 2020) and augment it with
gradient computation.

##  3\. Background

In this section, we briefly review the Projective Dynamics cloth simulation
method with dry frictional contact describe in Ly et al. (2020), on which our
differentiable cloth simulator is based. Our core contribution lies in the
development of gradient computation in this PD framework with dry frictional
contact, which we detail in the next section.

###  3.1. Projective Dynamics Method

#### Implicit time integration

We model cloth as a nodal system with $m$ 3D nodes. Let $\mathbf{x}(t)$ and
$\mathbf{v}(t)$ be two vector functions from $\mathbf{R}^{+}$ to
$\mathbf{R}^{3m}$ indicating the positions and velocities of all nodes at time
$t$. We consider the standard implicit time-stepping scheme discretizing
$\mathbf{x}(t)$ and $\mathbf{v}(t)$ as follows:

(1) $\displaystyle\mathbf{x}_{n+1}=$
$\displaystyle\mathbf{x}_{n}+h\mathbf{v}_{n+1},$ (2)
$\displaystyle\mathbf{v}_{n+1}=$
$\displaystyle\mathbf{v}_{n}+h\mathbf{M}^{-1}[\mathbf{f}_{\text{int}}(\mathbf{x}_{n+1})+\mathbf{f}_{\text{ext}}],$

where $h$ is the time step, $\mathbf{M}$ a positive diagonal mass matrix of
size $3m\times 3m$, and $\mathbf{f}_{\text{int}}$ and
$\mathbf{f}_{\text{ext}}$ the internal and external forces at each node
stacked as two $3m$-dimensional vectors, respectively. We use the cloth
material model described in Bouaziz et al. (2014) to define
$\mathbf{f}_{\text{int}}$. The implicit time-stepping scheme connects the
state of the nodal system $(\mathbf{x}_{n},\mathbf{v}_{n})$ at time $t_{n}$ to
$(\mathbf{x}_{n+1},\mathbf{v}_{n+1})$ at new time $t_{n+1}=t_{n}+h$.

#### Optimization view

As discussed by Martin et al. (2011), implicit time-stepping scheme can be
rephrased as finding a saddle point of the following energy minimization
problem:

(3)
$\displaystyle\min_{\mathbf{x}_{n+1}}\underbrace{\frac{1}{2h^{2}}(\mathbf{x}_{n+1}-\mathbf{y})^{\top}\mathbf{M}(\mathbf{x}_{n+1}-\mathbf{y})+W(\mathbf{x}_{n+1})}_{g(\mathbf{x}_{n+1})},$

where
$\mathbf{y}=\mathbf{x}_{n}+h\mathbf{v}_{n}+h^{2}\mathbf{M}^{-1}\mathbf{f}_{\text{ext}}$
is a constant vector that can be precomputed at the beginning of the current
timestep. The potential energy $W$ defines the internal force
$\mathbf{f}_{\text{int}}$ by its spatial gradients:
$\mathbf{f}_{\text{int}}=-\nabla W(\mathbf{x}_{n+1})$. It is easy to verify
that setting the gradient of the objective function to zero leads to a system
of equations identical to Eqns. (1) and (2).

#### Local and global solvers in PD

The key assumptions in PD is that the internal energy $W$ can be written as a
sum of quadratic forms (Bouaziz et al., 2014; Liu et al., 2017):

(4) $\displaystyle W_{i}(\mathbf{x})=$
$\displaystyle\min_{\mathbf{p}_{i}\in\mathcal{M}_{i}}\underbrace{\frac{w_{i}}{2}\|\mathbf{A}_{i}\mathbf{x}-\mathbf{p}_{i}\|_{2}^{2}}_{\tilde{W}_{i}(\mathbf{x},\mathbf{p}_{i})},$
(5) $\displaystyle W(\mathbf{x})=$ $\displaystyle\sum_{i}W_{i}(\mathbf{x}).$

Here, each energy $W_{i}$ projects $\mathbf{A}_{i}\mathbf{x}$, a linear
transformation of $\mathbf{x}$, to its closest point in the set
$\mathcal{M}_{i}$ and scales its squared distance by a prespecified stiffness
$w_{i}$. Both $\mathbf{A}_{i}$ and $\mathcal{M}_{i}$ are predefined and
independent of $\mathbf{x}$.

With such an assumption on $W$, PD proposes a local-global solver to minimize
a surrogate objective function:

(6)
$\displaystyle\tilde{g}(\mathbf{x}_{n+1},\mathbf{p})=\frac{1}{2h^{2}}(\mathbf{x}_{n+1}-\mathbf{y})^{\top}\mathbf{M}(\mathbf{x}_{n+1}-\mathbf{y})+\sum_{i}\tilde{W}_{i}(\mathbf{x}_{n+1},\mathbf{p}_{i}),$

where $\mathbf{p}$ stacks up all $\mathbf{p}_{i}$ from each $W_{i}$. In the
local step, PD fixes $\mathbf{x}_{n+1}$ and projects each $\mathbf{p}_{i}$ to
its corresponding $\mathcal{M}_{i}$ by solving Eqn. (4). Such a local step can
be done in parallel for each $\tilde{W}_{i}$. In the global step, PD fixes
$\mathbf{p}$ and minimizes $\tilde{g}$ as a function of $\mathbf{x}_{n+1}$,
which turns out to have a closed-form solution:

(7)
$\displaystyle\underbrace{(\mathbf{M}+h^{2}\sum_{i}w_{i}\mathbf{A}_{i}^{\top}\mathbf{A}_{i})}_{\mathbf{P}}\mathbf{x}_{n+1}=\underbrace{\mathbf{M}\mathbf{y}+h^{2}\sum_{i}w_{i}\mathbf{A}_{i}^{\top}\mathbf{p}_{i}}_{\mathbf{b}(\mathbf{p})}.$

By alternating between the local and global steps, PD monotonically decreases
the surrogate energy $\tilde{g}$ until convergence, which can be shown to
agree with the saddle point of $g$ (Liu et al., 2017). The source of
efficiency in PD comes from the observation that $\mathbf{P}$ is a constant
matrix that can be prefactorized at the beginning of simulation, leading to an
efficient global step requiring back-substitution only.

###  3.2. Dry Frictional Contact in Projective Dynamics

#### Signorini-Coulomb law

Ly et al. (2020) augments the standard PD framework described above with non-
penetration collision and Coulomb friction by assuming contact applies to
nodes only. At each timestep after a collision detection algorithm returns a
set of contact nodes, the following Signorini-Coulomb law (Brogliato, 2016)
needs to be enforced at each contact node $j$:

(8a) $\displaystyle[left=\empheqlbrace\,]$ $\displaystyle\text{(take off)
}\textbf{r}_{j}=\mathbf{0},\mathbf{u}_{j|N}>0,$ (8b)
$\displaystyle\text{(stick)}\|\textbf{r}_{j|T}\|<\mu\textbf{r}_{j|N},\mathbf{u}_{j}=\mathbf{0},$
(8c) $\displaystyle\text{(slip)
}\|\textbf{r}_{j|T}\|=\mu\textbf{r}_{j|N},\mathbf{u}_{j|N}=0,\exists\alpha>0,\mathbf{u}_{j|T}=-\alpha_{j}\textbf{r}_{j|T}.$

Here, $\mu$ is the frictional coefficient and $\mathbf{u}_{j}$ and
$\textbf{r}_{j}$ the nodal velocity and contact force represented in the local
contact frame spanned by the tangential plane and contact normal on the
contact surface. The notations $j|T$ and $j|N$ represent the tangential and
normal components of a 3D vector defined in the local frame.

#### Implicit time integration with dry frictional contact

The original implicit time integration now needs to be augmented with
additional constraints describing contact conditions (Ly et al., 2020):

(9)
$\left\\{\begin{array}[]{lr}\mathbf{x}_{n+1}=\mathbf{x}_{n}+h\mathbf{v}_{n+1},\\\
\mathbf{v}_{n+1}=\mathbf{v}_{n}+h\mathbf{M}^{-1}[\mathbf{f}_{\text{int}}(\mathbf{x}_{n+1})+\mathbf{f}_{\text{ext}}+\mathbf{J}_{n}^{\top}\textbf{r}_{n+1}],\\\
\mathbf{u}_{n+1}=\mathbf{J}_{n}\mathbf{v}_{n+1},\\\
(\textbf{r}_{n+1,j},\mathbf{u}_{n+1,j})\text{ satisfies
Eqn.}(\ref{eqn:contact_all}),\forall j\in\mathcal{C}_{n}.\end{array}\right.$

Here, the Jacobian matrix $\mathbf{J}_{n}$ of size $3|\mathcal{C}_{n}|\times
3m$ selects global vectors defined on the contact nodes and computes their
coordinates in the local contact frames. The contact set $\mathcal{C}_{n}$ at
this timestep is dynamically determined by a collision detection algorithm.

#### Projective Dynamics with dry frictional contact

The core idea in Ly et al. (2020) is an additional local step in PD that
solves each contact node independently. Noting that the contact conditions are
primarily defined on nodal velocities instead of nodal positions, Ly et al.
(2020) first rewrite the global step in Eqn. (7) as an equation of velocities:

(10)
$\displaystyle\mathbf{P}\mathbf{v}_{n+1}=\underbrace{\frac{1}{h}[\mathbf{b}(\mathbf{p})-\mathbf{P}\mathbf{x}_{n}]}_{\tilde{\mathbf{b}}(\mathbf{p})}.$

The velocity-based PD global-local steps then alternate between this new
global solver and the original local step, which is unaffected by this change
of variables. By comparing Eqn. (9) and Eqn. (10), we can see that
incorporating the contact force adds an additional impulse
$h\mathbf{J}_{n}^{\top}\textbf{r}_{n+1}$ to the right-hand side of the global
step:

(11)
$\displaystyle\mathbf{P}\mathbf{v}_{n+1}=\tilde{\mathbf{b}}(\mathbf{p})+h\mathbf{J}_{n}^{\top}\textbf{r}_{n+1}.$

The goal of the global solver is now to find
$(\mathbf{v}_{n+1},\textbf{r}_{n+1})$ that satisfy contact conditions at each
$j$. Noting that in Eqn. (11), $\mathbf{P}$ is the only operator that couples
unknown variables from different $j$, Ly et al. (2020) propose the following
Jacobi solver based on the decomposition $\mathbf{P}=\mathbf{M}+\mathbf{C}$:

(12)
$\displaystyle\mathbf{M}\mathbf{v}_{n+1}^{k+1}=\tilde{\mathbf{b}}{\mathbf{p}}-\mathbf{C}\mathbf{v}_{n+1}^{k}+h\mathbf{J}_{n}^{\top}\textbf{r}_{n+1},$

where the superscript $k$ and $k+1$ indicate two consecutive iterations in the
Jacobi solver. It is easy to see the proposed Jacobi solver fully decouples
the momentum equation on each node, allowing us to enforce the Signorini-
Columb law by adjusting $(\mathbf{v}_{n+1,j},\textbf{r}_{n+1,j})$ on each
contact node $j$ independently in a straightforward manner. When the Jacobi
solver converges, Ly et al. (2020) prove that the resultant solution is indeed
a solution to Eqn. (9).

##  4\. Differentiable Cloth Simulation

We now describe in detail how we extend the forward simulator in Sec. 3 to
build a differentiable cloth simulator. We start by differentiating implicit
time integration in Projective Dynamics without contact, followed by
explaining how contact gradients can be added to this framework. Compared with
other differentiable simulators, our simulator is very unique because of its
treatment on the contact-rich nature in cloth simulation: Many existing
differentiable simulators focus on sparse contact events in rigid-body or
deformable-body dynamics (Geilinger et al., 2020; Hu et al., 2019; Du et al.,
2021), and the state-of-the-art differentiable cloth simulation method (Liang
et al., 2019) handles rich contact events but without physics-based contact
forces or friction. To our best knowledge, our work is the first to present a
differentiable cloth simulator that can handle rich contact events with
Coulomb’s law of friction.

###  4.1. Gradients without Contact

#### Differentiating implicit time integration

The core step in building a differentiable simulator based on implicit time-
stepping scheme is to backpropgate gradients through the implicit integration
described in Eqns. (1) and (2), or equivalently, to derive the Jacobian of the
output $(\mathbf{x}_{n+1},\mathbf{v}_{n+1})$ with respect to the input
$(\mathbf{x}_{n},\mathbf{v}_{n+1})$. Mature techniques such as sensitivity
analysis, adjoint method, and implicit function theorem have proven to be
successful in computing such gradients (Hahn et al., 2019; Geilinger et al.,
2020; Du et al., 2021). To sketch the idea, we plug $\mathbf{v}_{n+1}$ from
Eqn. (2) into Eqn. (1):

(13) $\displaystyle\mathbf{x}_{n+1}=$
$\displaystyle\mathbf{x}_{n}+h\mathbf{v}_{n}+h^{2}\mathbf{M}^{-1}[\mathbf{f}_{\text{int}}(\mathbf{x}_{n+1})+\mathbf{f}_{\text{ext}}]$
$\displaystyle=$
$\displaystyle\mathbf{y}+h^{2}\mathbf{M}^{-1}\mathbf{f}_{\text{int}}(\mathbf{x}_{n+1}),$

which is essentially the first-order optimality condition in Eqn. (3).
Differentiating $\mathbf{x}_{n}$ on both sides gives

(14) $\displaystyle\frac{\partial\mathbf{x}_{n+1}}{\partial\mathbf{x}_{n}}=$
$\displaystyle[\mathbf{I}+h^{2}\mathbf{M}^{-1}\nabla^{2}W(\mathbf{x}_{n+1})]^{-1}$
$\displaystyle=$
$\displaystyle[\mathbf{M}+h^{2}\nabla^{2}W(\mathbf{x}_{n+1})]^{-1}\mathbf{M}.$

In backpropagation, such a Jacobian matrix is coupled with the gradients of a
loss function $L$ which are passed to the previous state (assuming gradients
below are both column vectors):

(15) $\displaystyle\frac{\partial L}{\partial\mathbf{x}_{n}}=$
$\displaystyle(\frac{\partial\mathbf{x}_{n+1}}{\partial\mathbf{x}_{n}})^{\top}\frac{\partial
L}{\partial\mathbf{x}_{n+1}}$ $\displaystyle=$
$\displaystyle\mathbf{M}\underbrace{[\mathbf{M}+h^{2}\nabla^{2}W(\mathbf{x}_{n+1})]^{-1}\frac{\partial
L}{\partial\mathbf{x}_{n+1}}}_{\mathbf{z}_{n+1}}.$

Here, we obtain the adjoint vector $\mathbf{z}_{n+1}$ by solving the linear
system of equations with the matrix
$\mathbf{M}+h^{2}\nabla^{2}W(\mathbf{x}_{n+1})$ on the left-hand side, which
avoids an explicit inversion of the large and sparse matrix. Backpropagating
gradients of $L$ with respect to $\mathbf{v}_{n}$ and $\mathbf{v}_{n+1}$ can
be derived similarly. In fact, it solves the same linear system but with a
different vector $\frac{\partial L}{\partial\mathbf{v}_{n+1}}$ on the right-
hand side.

#### Differentiating with Projective Dynamics

With assumptions on the energy form $W$ in Projective Dynamics, Du et al.
(2021) show that we can speed up the computation in $\mathbf{z}_{n+1}$ by
exploiting the special form of $\nabla^{2}W$. The first- and second-order
derivatives in Eqn. (5) are given by the following equations (Du et al.,
2021):

(16) $\displaystyle\nabla W_{i}(\mathbf{x})=$ $\displaystyle
w_{i}\mathbf{A}_{i}^{\top}(\mathbf{A}_{i}\mathbf{x}-\mathbf{p}_{i}^{*}(\mathbf{x})),$
(17) $\displaystyle\nabla^{2}W_{i}(\mathbf{x})=$ $\displaystyle
w_{i}\mathbf{A}_{i}^{\top}\mathbf{A}_{i}-w_{i}\mathbf{A}_{i}^{\top}\frac{\partial\mathbf{p}_{i}^{*}}{\partial\mathbf{x}},$

where
$\mathbf{p}_{i}^{*}(\mathbf{x})=\arg\min_{\mathbf{p}_{i}\in\mathcal{M}_{i}}\tilde{W}_{i}(\mathbf{x},\mathbf{p}_{i})$
is the projection of $\mathbf{A}_{i}\mathbf{x}$ onto $\mathcal{M}_{i}$
obtained in the local step. We plug them to the definition of
$\mathbf{z}_{n+1}$ and rewrite the linear system as follows:

(18)
$\displaystyle(\mathbf{P}-\underbrace{h^{2}\sum_{i}w_{i}\mathbf{A}_{i}^{\top}\frac{\partial\mathbf{p}_{i}^{*}}{\partial\mathbf{x}}\Bigr{|}_{\mathbf{x}_{n+1}}}_{\Delta\mathbf{P}})\mathbf{z}_{n+1}=\frac{\partial
L}{\partial\mathbf{x}_{n+1}}.$

Noting that a prefactorization of $\mathbf{P}$ exists in Projective Dynamics,
Du et al. (2021) propose the following Jacobi method to solve
$\mathbf{z}_{n+1}$:

(19)
$\displaystyle\mathbf{P}\mathbf{z}_{n+1}^{k+1}=\Delta\mathbf{P}\mathbf{z}_{n+1}^{k}+\frac{\partial
L}{\partial\mathbf{x}_{n+1}}.$

Similarly to the local-global steps in Projective Dynamics, a local step can
evaluate $\Delta\mathbf{P}$ across each $\nabla^{2}W_{i}$ in parallel, and a
global step, which reuses $\mathbf{P}$, can solve $\mathbf{z}_{n+1}^{k+1}$
with backsubstitution only. Du et al. (2021) show that such a local-global
solver is empirically faster than directly solving Eqn. (15), which resembles
the source of efficiency in the original Projective Dynamics method.

###  4.2. Gradients with Contact

While Du et al. (2021) have discussed extensively how to fuse backpropagation
into the Projective Dynamics framework, its support of contact is limited to
non-penetration conditions without a proper treatment on friction. On the
other hand, other prior papers have provided solutions to deriving gradients
with contacts governed by Signorini-Coulomb law, but they focus on either
small-scale problems (e.g., rigid bodies in de Avila Belbute-Peres et al.
(2018a)) or sparse contact events on a static ground only (Geilinger et al.,
2020). Here, we present our differentiable cloth simulator that both inherits
the speedup from Projective Dynamics and handles gradients with complicated
contact events like self-collisions, which are largely overlooked in existing
differentiable simulators but very common in cloth simulation.

#### Differentiating implicit time integration with contact

We first rewrite the implicit time integration in Eqn. (9) by removing
$\mathbf{x}_{n+1}$ and all inactive inequality constraints in Eqn. (8):

(20)
$\left\\{\begin{array}[]{lr}\mathbf{v}_{n+1}=\mathbf{v}_{n}+h\mathbf{M}^{-1}[\mathbf{f}_{\text{int}}(\mathbf{x}_{n}+h\mathbf{v}_{n+1})+\mathbf{f}_{\text{ext}}+\mathbf{J}_{n}^{\top}\textbf{r}_{n+1}],\\\
\mathbf{C}_{n}(\mathbf{J}_{n}\mathbf{v}_{n+1},\textbf{r}_{n+1})=\mathbf{0},\end{array}\right.$

where
$\mathbf{C}(\mathbf{u},\textbf{r}):\mathbf{R}^{3|\mathcal{C}_{n}|}\times\mathbf{R}^{3|\mathcal{C}_{n}|}\rightarrow\mathbf{R}^{3|\mathcal{C}_{n}|}$
encodes all the equality constraints in Eqn. (8) when $\mathbf{x}_{n+1}$ and
$\mathbf{v}_{n+1}$ are solved in forward simulation. The output dimension of
$\mathbf{C}$ is justified by the fact that each case in Eqn. (8) has exactly
three equality conditions. By definition, $\mathbf{C}_{n}$ is separable
because the contact condition on each contact node depends on its own velocity
and contact force only. Differentiating with respect to $\mathbf{v}_{n}$ on
both sides of the first equation gives

(21)
$\displaystyle\frac{\partial\mathbf{v}_{n+1}}{\partial\mathbf{v}_{n}}=[\mathbf{M}+h^{2}\nabla^{2}W(\mathbf{x}_{n+1})-\underbrace{h\mathbf{J}_{n}^{\top}\frac{\partial\textbf{r}_{n+1}}{\partial\mathbf{v}_{n+1}}}_{\Delta\mathbf{R}^{\top}}]^{-1}\mathbf{M}.$

Comparing it to Eqn. (14), we see the matrix to be inverted now has an
additional component dependent on
$\frac{\partial\textbf{r}_{n+1}}{\partial\mathbf{v}_{n+1}}$, which we obtain
from differentiating the constraint $\mathbf{C}_{n}=\mathbf{0}$:

(22)
$\displaystyle\frac{\partial\textbf{r}_{n+1}}{\partial\mathbf{v}_{n+1}}=-[\frac{\partial\mathbf{C}_{n}}{\partial\textbf{r}}\Bigr{|}_{\textbf{r}_{n+1}}]^{+}\frac{\partial\mathbf{C}_{n}}{\partial\mathbf{u}}\Bigr{|}_{\mathbf{J}_{n}\mathbf{v}_{n+1}}\mathbf{J}_{n},$

where $[\cdot]^{+}$ indicates the Moore-Penrose inverse of a matrix. We stress
that computing both partial derivatives of $\mathbf{C}_{n}$ and the inverse
are trivial because $\mathbf{C}_{n}$ is separable. Indeed, both partial
derivatives are $3\times 3$ block diagonal matrices that can be parallelized
among contact nodes. Backpropagation through $\mathbf{v}_{n+1}$ to
$\mathbf{v}_{n}$ can be implemented with the same adjoint method before, which
we give in the equation below for completeness (the notation
$\mathbf{z}_{n+1}$ is overloaded):

(23) $\displaystyle\frac{\partial
L}{\partial\mathbf{v}_{n}}=\mathbf{M}\underbrace{[\mathbf{M}+h^{2}\nabla^{2}W(\mathbf{x}_{n+1})-\Delta\mathbf{R}]^{-1}\frac{\partial
L}{\partial\mathbf{v}_{n+1}}}_{\mathbf{z}_{n+1}}.$

#### Speedup with Projective Dynamics

With additional information about $W$ from Projective Dynamics, we can rewrite
the adjoint vector $\mathbf{z}_{n+1}$ in Eqn. (23) by comparing it to Eqn.
(18):

(24)
$\displaystyle(\mathbf{P}-\Delta\mathbf{P}-\Delta\mathbf{R})\mathbf{z}_{n+1}=\frac{\partial
L}{\partial\mathbf{v}_{n+1}},$

which naturally leads to the following Jacobi iteration:

(25)
$\displaystyle\mathbf{P}\mathbf{z}_{n+1}^{k+1}=(\Delta\mathbf{P}+\Delta\mathbf{R})\mathbf{z}_{n+1}^{k}+\frac{\partial
L}{\partial\mathbf{v}_{n+1}}.$

Comparing it to Eqn. (19), we see the role of $\Delta\mathbf{P}$ is replaced
with $\Delta\mathbf{P}+\Delta\mathbf{R}$. Since we have shown that computing
$\Delta\mathbf{R}$ can be massively parallelized among contact nodes, it is
suitable for contact-rich cloth simulation and preserves the efficiency of the
local-global solver in the Jacobi iteration.

#### Convergence

In theory, the Jacobi iteration is guaranteed to converge from any initial
guesses of $\mathbf{z}_{n+1}$ when the spectral radius
$\rho[\mathbf{P}^{-1}(\Delta\mathbf{P}+\Delta\mathbf{R})]<1$. Empirically, we
notice in our cloth simulation that divergence is uncommon especially when
high-precision simulation results are not required (Sec. 6). When the Jacobi
solver fails to converge, we switch back to a direct sparse matrix solver to
solve Eqn. (24).

#### Extensions

We deliberately skipped the definition of the Jacobian matrix $\mathbf{J}_{n}$
in all equations above in order to clearly present the main idea behind our
differentiable cloth simulation. We now elaborate on the role of
$\mathbf{J}_{n}$ and discuss its implications on different types of
collisions. In the simplest case where the contact surface is a static flat
plane, $\mathbf{J}_{n}$ is a constant matrix that can be assembled beforehand,
and no further modifications are needed to the derivation of gradients so far.
A slightly more interesting case is when the contact surface is static but
non-planar, e.g., a spherical surface. Such a contact surface has spatially
varying surface normals. Therefore, $\mathbf{J}_{n}$ is dependent on the
positions of the nodes where the contact events occur. In our implementation,
we define $\mathbf{J}_{n}$ as a function of $\mathbf{x}_{n}$, i.e., we
estimate the contact normal based on the positions of the contact nodes at the
beginning of the timestep when the contact events occur. While it is
technically possible to define $\mathbf{J}_{n}$ with the exact contact
location computed based on $\mathbf{x}_{n}$, $\mathbf{x}_{n+1}$, or
$\mathbf{v}_{n+1}$, we have found that defining $\mathbf{J}_{n}$ as a function
of $\mathbf{x}_{n}$ simplifies the gradient computation substantially without
introducing significant visual artifacts in our cloth simulation. Replacing
$\mathbf{J}_{n}$ with $J_{n}(\mathbf{x}_{n})$ does not require extra changes
to the gradients in Sec. 4.2, as it has no dependencies on $\mathbf{v}_{n}$ or
$\mathbf{v}_{n+1}$.

The last type of collisions we consider is self-collisions between nodes. In
this case, contacts occur between pairs of nodes where the contact normals are
defined on the relative position between two nodes in the pair. Therefore, we
can still define $\mathbf{J}_{n}$ as a function of $\mathbf{x}_{n}$ with each
row block now consists of two blocks of nonzero elements corresponding to the
two contact nodes in a pair. The gradient derivation in Sec. 4.2 remains
unchanged in this case. As we inherit from Ly et al. (2020) the contact model
on nodes only, we also inherit its limitation of not handling other types of
self-collisions like edge-edge collisions or vertex-face collisions. We leave
the derivation of gradients with such cases as future work.

##  5\. Evaluations

The two key components that distinguish our differentiable cloth simulator
from others are its differentiable contact model and its Jacobi iterations in
backpropagation. We devote this section to a comprehensive analysis on these
two key ingredients.

#### Contact model in differentiable cloth simulation

We compare the contact model in our differentiable simulator with that in the
state-of-the-art differentiable cloth simulator Liang et al. (2019). Both
models ensures penetration-free simulation and detects self-collisions (node-
node collisions in ours and vertex-face and edge-edge collisions in Liang et
al. (2019)). However, that contact model in Liang et al. (2019) does not take
into account complementarity conditions on either contact forces or frictional
forces, which may lead to undesirable artifacts. To visualize this issue, we
consider a test scenario in which a napkin falls freely onto the inner surface
of a bowl (Fig. 2). The differentiable simulators from both Liang et al.
(2019) and our work manage to simulate the napkin without numerical explosion.
However, the dry frictional contact model in our differentiable simulator
shows more physically realistic motion (Fig. 2 top), whereas the contact model
from Liang et al. (2019) leads to more drastic changes in the size of the
napkin with a popping artifact after its contact with the bowl (Fig. 2 middle
and bottom). These artifacts are explainable because the collision handling
algorithm in Liang et al. (2019) modifies node positions after penetration
without verifying whether such an update requires sticky contact forces. When
the napkin becomes in contact with the concave inner surface of the bowl, such
a direct modification injects extra elastic energies. This experiment confirms
supporting a more advanced contact and friction model in differentiable cloth
simulation is both viable and beneficial.

Figure 2. Napkin. We simulate the motion sequence of the napkin ($50\times
50\times 2$ triangles, $h=5$ms, $500$ timesteps) with the contact models from
two differentiable cloth simulators: ours (top) and Liang et al. (2019)
(middle). The area ratio (bottom) is the ratio between the napkin’s area at
the current timestep and in the rest shape. Please refer to our video for the
full motions.

#### Iterative vs. direct solver in backpropagation

Another key component in our differentiable cloth simulation is the usage of
the Jacobi iteration in Eqn. (24) that utilizes the prefactorized $\mathbf{P}$
in Projective Dynamics. In a standard differentiable simulator, solving Eqn.
(24) would be done with a sparse matrix solver treating
$(\mathbf{P}-\Delta\mathbf{P}-\Delta\mathbf{R})$ as a whole. To understand the
performance of the proposed Jacobi solver, we design two benchmark tests: a
“Wind” test where a hanging napkin moves under synthetic wind and a “Slope”
test where a ribbon slides on a slanted plane (Fig. 3). These two tests
represent two extreme cases in contact handling: the napkin in “Wind” barely
has any contact or self-collisions, whereas every single node of the ribbon in
“Slope” is in contact with the plane. We then compare the time cost of our
Jacobi solver with a sparse LU solver in both tests. We implement both solvers
using Eigen (Guennebaud et al., 2010) and choose LU factorization because
$(\mathbf{P}-\Delta\mathbf{P}-\Delta\mathbf{R})$ is usually not a symmetric or
positive definite matrix, preventing us from using more specialized sparse
matrix solvers like Cholesky factorization or Conjugate Gradient methods. For
the Jacobi solver, we report two results from low-precision (1e-4) and high-
precision (1e-6) convergence threshold that controls the termination condition
in the Jacobi iterations. We find these thresholds by varying it from 1e-1 to
1e-9 until the gradients computed from the Jacobi solver start to agree with
those returned by the direct solver, which we treat as the ground-truth
gradients. We repeat the experiments with three mesh resolutions to test the
scalability of our method.

We summarize the statistics from both tests in Table 1 and 2. Overall, using
Jacobi iteration in backpropagation is faster than using the direct solver,
and the speedup becomes more substantial as we increase the mesh resolution.
The speedup is also more evident with low-precision threshold, which is
consistent with Projective Dynamics. A further decomposition of time confirms
that solving the linear system takes up the majority of backpropagation time,
so any improvement in the choice of solver has a dominating positive effect.
Although the low-precision results may not match the ground-truth gradients as
closely as their high-precision counterparts, their backpropagation is
significantly faster and can still benefit gradient-based optimization. This
is because even an imperfect gradient can still guide gradient-descent
algorithms to converge as long as it is along the descending direction. This
is confirmed empirically by our experiments in Sec. 6, in most of which we use
low-precision convergence threshold and still observe successful gradient-
based optimization results. It is also worth mentioning that contact-related
operators, whose time cost is included in the “Other” and “Contact” columns in
Table 1 and 2, add very little extra computation to backpropagation. Finally,
we observe uncommon but non-negligible failure cases of the Jacobi iterations
in one experiment ($48\times 48$ with 1e-6 as the threshold in Table 1). For
the $15\%$ failed timesteps, we switch to the sparse LU solver, adding an
extra $11.5\%$ time in backpropagation.

Table 1. Comparisons between the Jacobi iteration (ours) and the sparse LU direct solver in backpropagation under various mesh resolutions in the “Wind” test. The numbers in the “Res.” column report the mesh resolution. The number in parenthesis in the “Solver” column indicates the epsilon value controlling the convergence of the Jacobi solver. The “Backprop. Time” column reports the net time in backpropagation and is further decomposed into the next four columns, from which the sum of ratios is $100\%$: “$\Delta\mathbf{P}$” reports the time spent on assembling $\Delta\mathbf{P}$, “Jacobi Solve” and “Direct Solve” shows the time cost by either solver, and operations not covered by the three categories are in the “Other” column. The “Failure” column reports the ratio between the number of timesteps seeing nonconvergence in the Jacobi iterations and the number of total timesteps. The “Speedup” column is the ratio between the direct solver’s time and our time in “Backprop. Time” in each test. Res. | Solver | Backprop. Time (s) |  $\Delta\mathbf{P}$ (%) | Jacobi Solve (%) | Direct Solve (%) | Other (%) | Failure (%) | Speedup
---|---|---|---|---|---|---|---|---
12x12 | Direct | 0.606 | 13.9 | - | 84.4 | 1.7 | 0 | -
Ours (1e-4) | 0.212 | 39.9 | 55.5 | - | 4.6 | 0 | 2.9x
Ours (1e-6) | 0.424 | 19.2 | 78.5 | - | 2.3 | 0 | 1.4x
24x24 | Direct | 3.360 | 9.9 | - | 89.1 | 1.0 | 0 | -
Ours (1e-4) | 0.591 | 58.2 | 36.3 | - | 5.5 | 0 | 5.7x
Ours (1e-6) | 2.858 | 11.7 | 87.2 | - | 1.1 | 0 | 1.2x
48x48 | Direct | 24.001 | 6.8 | - | 92.5 | 0.7 | 0 | -
Ours (1e-4) | 2.262 | 63.2 | 29.7 | - | 7.1 | 0 | 10.6x
Ours (1e-6) | 29.255 | 5.4 | 82.6 | 11.5 | 0.6 | 15 | 0.8x
Table 2. Comparisons between the Jacobi iteration (ours) and the sparse LU direct solver in backpropagation under various mesh resolutions in the “Slope” test. The numbers in the “Res.” column report the mesh resolution. The number in parenthesis in the “Solver” column indicates the epsilon value controlling the convergence of the Jacobi solver. The “Backprop. Time” column reports the net time in backpropagation and is further decomposed into the next five columns, from which the sum of ratios is $100\%$: “$\Delta\mathbf{P}$” reports the time spent on assembling $\Delta\mathbf{P}$, “Jacobi Solve” and “Direct Solve” shows the time cost by either solver, “Contact” reports the time cost by computing gradients in the contact model, and operations not covered by the three categories are in the “Other” column. The Jacobi solver converges without any failure cases in this test. The “Speedup” column is the ratio between the direct solver’s time and our time in “Backprop. Time” in each test. Res. | Solver | Backprop. Time (s) |  $\Delta\mathbf{P}$ (%) | Jacobi Solve (%) | Direct Solve (%) | Contact (%) | Other (%) | Speedup
---|---|---|---|---|---|---|---|---
12x12 | Direct | 0.978 | 13.0 | - | 82.4 | 3.4 | 1.2 | -
Ours (1e-4) | 0.312 | 70.1 | 16.1 | - | 9.4 | 4.5 | 3.1x
Ours (1e-6) | 0.643 | 18.9 | 76.8 | - | 2.7 | 1.6 | 1.5x
24x24 | Direct | 5.331 | 10.4 | - | 87.4 | 1.5 | 0.7 | -
Ours (1e-4) | 0.781 | 70.1 | 16.1 | - | 9.4 | 4.5 | 6.8x
Ours (1e-6) | 1.507 | 34.7 | 58.1 | - | 4.9 | 2.4 | 3.7x
48x48 | Direct | 37.296 | 6.5 | - | 92.5 | 0.7 | 0.4 | -
Ours (1e-4) | 3.095 | 68.1 | 19.8 | - | 8.2 | 3.8 | 12.0x
| Ours (1e-6) | 6.825 | 31.2 | 63.6 | - | 3.5 | 1.8 | 5.5x
Figure 3. Wind and Slope. Motion sequences from the two benchmark tests for
comparing iterative and direct solvers in backpropagation. Top: the “Wind”
test ($h=1/90$s, 200 timesteps) where a napkin moves under synthetic wind.
Bottom: the “Slope” test ($h=1/100$s, 300 timesteps) where a ribbon slides
along a slanted plane.

##  6\. Applications

In this section, we demonstrate a variety of cloth-related applications that
can benefit from our proposed differentiable cloth simulation with dry
frictional contact. We run all optimizations on a workstation of 80 CPU cores
and 80G memory. Depending on problem complexity, the wallclock time of running
these applications varies from less than 30 minutes to 2 hours for all
methods. We write the backbone of our simulator in C++ and use Eigen for
matrix and vector operations. Each application defines an optimization problem
which we solve with L-BFGS-B, a classic gradient-based optimizer that can
leverage the differentiability of our cloth simulator. We report the
optimization results and comparisons with other gradient-free algorithms in
Table 3. Below we describe each application in detail.

Table 3. Comparisons between the performance of gradient-based and gradient-free optimization on all examples in Sec. 6. The “Init. Loss” and “Final Loss” columns report the loss before and after optimization, which is averaged over multiple random seeds if possible. The “Forward” and “Backward” columns report the timesteps in forward simulation and backpropagation used by each optimization strategy. Note that “Backward” columns are not included in CMA-ES, (1+1)-ES, or RL methods because they do not require gradients. Name | Init. Loss | Gradient-Based (Ours) | CMA-ES | (1+1)-ES | RL
---|---|---|---|---|---
Final Loss | Forward | Backward | Final Loss | Forward | Final Loss | Forward | Final Loss | Forward
T-shirt | 26.07 | 1.39 | 4000 | 4000 | 2.05 | 8000 | 1.30 | 8000 | - | -
Sphere | 1.63 | 0.38 | 3000 | 3000 | 0.01 | 6000 | 0.06 | 6000 | - | -
Hat | 40.46 | 2.22 | 6000 | 6000 | 19.54 | 12000 | 10.51 | 12000 | 6.79 | 1.2e6
Sock | 35.34 | 3.86 | 15000 | 15000 | 11.25 | 30000 | 7.25 | 30000 | - | -
Dress | 1.49 | 0.77 | 4000 | 4000 | 0.96 | 8000 | 0.84 | 8000 | - | -
Flag | 4.42 | 0.82 | 10000 | 10000 | 0.90 | 20000 | 1.05 | 20000 | - | -

###  6.1. System Identification

We start by showing two system identification examples: “T-shirt” (Fig. 4) and
“Sphere” (Fig. 5). In the “T-shirt” example, we are given a sequence of
motions of a hanging T-shirt under synthetic wind generated from a sinusoidal
function. The goal is to estimate a material parameter in the cloth (1 DoF)
and identify the wind parameters (5 DoFs controlling the amplitude, the phase,
and the frequency of the sinusoidal function in three dimensions) from the
motion data. We define the loss function as the L2-distance between nodal
positions of the T-shirt from simulation and the given motion sequence. Unlike
the “Wind” example in Sec. 5, contacts and frictions are much more frequent in
this example due to self-collision between the front and back layer of the
T-shirt under the wind. We show the simulation of the T-shirt with parameters
before and after optimization in Fig. 4 and compare three optimizers in Table
3: L-BFGS-B, CMA-ES, and (1+1)-ES. Both CMA-ES and (1+1)-ES are standard
gradient-free evolutionary strategies (ES). We can conclude from Fig. 4 and
Table 3 that all three methods manage to optimize system parameters leading to
motion sequences visually identical to the given input, but L-BFGS-B converges
much faster due to the extra knowledge of gradients.

To highlight the effect of dry frictional contact in our simulator, we create
the “Sphere” example with the goal of matching a motion sequence of a cloth
interacting with a sphere by estimating the frictional coefficient of the
sphere. In this example, we let a piece of square cloth fall freely on a
sphere, whose motion after being in contact with the sphere is largely
controlled by the frictional coefficient. This example involves both self-
collisions between nodes on the cloth and contacts with the sphere. Similar to
the example above, we run both L-BFGS-B and two ES baselines and report their
statistics in Table 3. All methods can optimize to a frictional coefficient
that generates a motion sequence visually identical to the given input.
However, this optimization problem is special in that it only has one
parameter. In this case, evolutionary strategies can achieve a lower final
loss because the samples can freely explore the full range of frictional
contact, while L-BFGS-B can get trapped in a local optimum.

Figure 4. T-shirt (4278 DoFs, $h=1/90$s, 250 timesteps). Estimating the cloth
material and wind parameters based on a given synthetic motion sequence of the
T-shirt. From left to right we show the simulated T-shirt at 0, 80, 160, and
250 timesteps. Top: the ground-truth motion sequence; Middle: simulation with
the initial guess on the cloth and wind parameters; Bottom: simulation after
optimization with L-BFGS-B. Figure 5. Sphere (1875 DoFs, $h=1/180$s, 350
timesteps). Estimating the frictional coefficient of a piece of square cloth
based on its contact with a sphere. From top to bottom we show the simulation
at 0, 100, 200, and 300 timesteps. Top: the ground-truth motion sequence;
Middle: simulation with the initial guess on the cloth parameters; Bottom:
simulation after optimization with L-BFGS-B.

###  6.2. Robot-Assisted Dressing

Another line of research that can benefit from a differentiable cloth
simulator is robot-assisted dressing. The mainstream solution to these tasks
is typically gradient-free methods like evolutionary strategies or
reinforcement learning before a differentiable cloth simulator becomes
available. We present two examples to demonstrate the usage of gradients in
robot-assisted dressing: a “Hat” (Fig. 6) and a “Sock” (Fig. 7). In both
examples, the goal is to find trajectories for a kinematic robotic manipulator
to put on the hat or the sock. The end effector of the manipulator picks a few
prespecified vertices on the cloth meshes and pulls them along kinematic
trajectories, which are parametrized as B-splines. By optimizing the control
points of the B-splines (18 DoFs in “Hat” and 36 DoFs in “Sock”), we can
direct the manipulator to move the hat until it reaches the target location on
top of the sphere. We define the loss function in “Hat” as the L2-distance
between the hat’s final position at the last timestep and a predefined target
position, which we generate by translating the hat’s rest shape onto the top
of the sphere. The loss function for “Sock” is defined as the L2-distance
between the desired and simulated locations of a few key points manually
chosen on the sock and evaluated at the middle and the end of the simulation.

With the gradient information at hand, we run the L-BFGS-B optimizer to tune
the control points of the trajectories and compare its performance with CMA-ES
and (1+1)-ES (Table 3). We notice that within the same time budget, L-BFGS-B
converges substantially faster than the gradient-free baselines with a better
solution (Fig. 6 and Fig. 7).

Another widely used strategy in solving robot-assisted dressing is deep
reinforcement learning (RL), which typically formulates such tasks as a Markov
Decision Process (MDP) and trains a closed-loop neural network controller.
Although RL and the above optimization methods solve robot-assisted dressing
from different perspectives and comparisons are not well-grounded, we try our
best to tune PPO (Schulman et al., 2017), a state-of-the-art RL algorithm, on
both tasks and provide its statistics in Table 3. RL attempts to solve a more
difficult closed-loop control problem, and understandably, the performance of
the solution found by RL is worse than the optimization methods above within a
similar time budget.

By comparing to both ES and RL, we can safely conclude that the fast
convergence in L-BFGS-B unlocked by our differentiable cloth simulator is a
clear advantage over gradient-free methods. Additionally, it may benefit from
exploration strategies in ES or RL to avoid local minima, which we leave as
future work.

Figure 6. Hat (1737 DoFs, $h=1/100$s, 400 timesteps). Optimizing trajectories
for a manipulator so that it moves a hat onto the sphere. Top left: one
initial trajectory before optimization overlaid with intermediate hat
positions in simulation. Top right: the optimized trajectory from L-BFGS-B
(which shares a similar behavior as the ones optimized from ES algorithms).
Bottom left: RL algorithm converges to a failure control trajectory in some
seeds. Bottom right: the loss-timestep curves for all methods. Figure 7. Sock
(3165 DoFs, $h=1/160$s, 400 timesteps). Optimizing trajectories for a
manipulator so that it puts on a sock to the foot model (rigid two-capsule
object in grey). Left: the starting locations of the foot and the sock.
Middle: the control trajectory optimized by L-BFGS-B and overlaid with
intermediate sock positions. Right: the optimized trajectory successfully put
the sock to the foot. Bottom right: the loss-timestep curves for all methods.

###  6.3. Inverse Design

The next application we present is “Dress”, an inverse design example that
aims to optimize cloth material parameters in a dress so that its dynamic
motion can satisfy certain design intents. Specifically, we consider a
pirouetting ballet dancer wearing a dress with material parameters to be
optimized so that after the dress spins, the apex angle of the cone-like dress
agrees with a desired value (50 degrees in our experiment). We define the loss
function as the difference between the hemline height corresponding to the
desired apex angle and the estimated apex angle from points on the hemline of
the dress at the last frame of the simulation. We report the optimization
results from L-BFGS-B and two ES baselines in Table 3 and visualize the
simulation results before and after optimization in Fig. 8. We notice a
speedup of 15x (measured by the number of timesteps) from L-BFGS-B over ES
with better optimized results.

Figure 8. Dress (10902 DoFs, $h=1/120$s, 125 timesteps). Designing the
material parameters of a dress for a pirouetting ballet dancer so that the
apex angle of the cone formed by the spinning dress reaches a desired value
(50 degrees in our experiment). The leftmost figure shows the rest shape of
the dress under gravity, followed by the spinning dress at 125 timesteps with
material parameters before (middle) and after optimization using L-BFGS-B
(right).

###  6.4. A Real-to-Sim Example

We end this section with a real-to-sim “Flag” example (Fig. 9). In this
example, we use the real-world motion sequence captured on a flag flapping in
the wind from previous work (White et al., 2007) and aim to reconstruct a
digital twin of the scene in simulation. This includes not only estimating the
material parameters of the flag but also modeling the wind condition, which is
particularly challenging due to its intricate, stochastic model with unknown
degrees of freedom. We model the wind force at each time step as a 3D force
applied near the center of the scene and spatially decayed proportional to the
inverse distance to the center. To model the transient nature of the wind
force, we modulate magnitude of the 3D vector with a sinusoidal wave as a
function of time with trainable frequency and phase offset. Together, the
material and wind models define an eight-dimensional parameter space to be
optimized. We define the loss function as the L2-distance between the
positions of all nodes from each timestep in simulation and the ground-truth
motion sequence.

We solve this task using L-BFGS-B and the two ES methods and report their
performance in Table 3 measured when all methods are given the same number of
total timesteps (20000 in this case). The loss numbers are averaged over 10
random seeds that sample the initial guess of all parameters. All three
methods substantially reduce the loss after optimization and achieve a similar
finial loss, with L-BFGS-B performing slightly better. We visualize the motion
of the simulated flag before and after optimization in Fig. 9, in which we
also plot the trajectories of a few key nodes with reference trajectories in
blue from the motion-captured data. By comparing the left and right images in
Fig. 9, we can see that the L-BFGS-B optimizer reduces the discrepancy between
simulated and actual trajectories after optimization. The sim-to-real matching
is still imperfect as indicated by the nonzero final loss, which we suspect
could be due to the capacity of our synthetic wind model. A more sophisticated
and expressive wind model, e.g., a neural network, may serve a better role in
modeling and matching the real-world physics, which we leave as future work.

Figure 9. Flag. A real-to-sim example that reconstructs a digital flapping
flag ($540$ vertices, $1026$ triangles, $h=1/120$s, 100 timesteps) based on
motion data captured from real-world experiments. We plot the trajectories of
a few key points on the cloth from the ground-truth motion (blue curves) and
the simulation results with initial guesses on the material and wind
parameters (red curves, left) and after optimization (green curves, right).

##  7\. Conclusions

In this paper, we presented a differentiable cloth simulator built on
Projective Dynamics with Signorini-Coulomb frictional contact. Our
differentiable simulator is different from existing papers in its simultaneous
accommodation of rich and frequent (self-)contact, Signorini-Coulomb contact
law, and differentiability in cloth simulation. We additionally presented a
Jacobi solver that exploits the contact gradients to speed up the
backpropagation and observed a substantial speedup (up to an order of
magnitude with low-precision simulation) over direct solvers. Our
differentiable cloth simulator enabled gradient-based optimization methods in
a diverse set of applications, for which traditional gradient-free methods are
generally much less sampling efficient.

#### Limitations and future work

There are still quite a few limitations in our method that are worth further
investigation. First, since our method is built on Projective Dynamics, it
also limits the choice of material models. It would be useful to generalize
the current framework and support more physically accurate cloth material
models, e.g., the piecewise linear elastic model described in Wang et al.
(2011).

Second, the contact model we use from Ly et al. (2020) does not take into
account vertex-face or edge-edge self-collisions. While we empirically
observed that handling only vertex-vertex self-collisions managed to produce
plausible results with medium mesh resolutions, having vertex-face and edge-
edge collision detection and handling is still highly desirable for a more
physically realistic cloth simulator.

Third, there is no theoretical guarantee on the convergence of the Jacobi
iterations we implemented in backpropagation. Although non-convergence is
uncommon in our experiments, it introduces a costly switch to the slower
direct solver, which we hope to fully resolve in future work.

Finally, many of our applications were in simulation only. It would be more
exciting if these results could be replicated in real-world settings. We
consider connecting our differentiable cloth simulator to more real-world
applications, including real-world robot-assisted dressing, material parameter
identification for real-world fabric samples, computational design for sports
suits, and so on. Closing the sim-to-real gap is a non-trivial problem, in
which we believe our differentiable simulator could play a beneficial role.

## References

  * (1)
  * Baraff and Witkin (1998) David Baraff and Andrew Witkin. 1998.  Large Steps in Cloth Simulation. In _Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques_ _(SIGGRAPH ’98)_. Association for Computing Machinery, New York, NY, USA, 43–54.
  * Bartle et al. (2016) Aric Bartle, Alla Sheffer, Vladimir G. Kim, Danny M. Kaufman, Nicholas Vining, and Floraine Berthouzoz. 2016.  Physics-Driven Pattern Adjustment for Direct 3D Garment Editing.  _ACM Trans. Graph._ 35, 4, Article 50 (July 2016), 11 pages.
  * Bouaziz et al. (2014) Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014.  Projective Dynamics: Fusing Constraint Projections for Fast Simulation.  _ACM Trans. Graph._ 33, 4, Article 154 (July 2014), 11 pages.  https://doi.org/10.1145/2601097.2601116
  * Bridson et al. (2002) Robert Bridson, Ronald Fedkiw, and John Anderson. 2002\.  Robust Treatment of Collisions, Contact and Friction for Cloth Animation.  _ACM Trans. Graph._ 21, 3 (July 2002), 594–603.
  * Bridson et al. (2003) R. Bridson, S. Marino, and R. Fedkiw. 2003.  Simulation of Clothing with Folds and Wrinkles. In _Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (San Diego, California) _(SCA ’03)_. Eurographics Association, Goslar, DEU, 28–36.
  * Brogliato (2016) Bernard Brogliato. 2016\.  _Nonsmooth Lagrangian Systems_.  Springer International Publishing, Cham, 241–370.
  * Brouet et al. (2012) Remi Brouet, Alla Sheffer, Laurence Boissieux, and Marie-Paule Cani. 2012.  Design Preserving Garment Transfer.  _ACM Trans. Graph._ 31, 4, Article 36 (July 2012), 11 pages.
  * Brown et al. (2018) George E. Brown, Matthew Overby, Zahra Forootaninia, and Rahul Narain. 2018.  Accurate Dissipative Forces in Optimization Integrators.  _ACM Trans. Graph._ 37, 6, Article 282 (Dec. 2018), 14 pages.
  * Choi and Ko (2002) Kwang-Jin Choi and Hyeong-Seok Ko. 2002.  Stable but Responsive Cloth.  _ACM Trans. Graph._ 21, 3 (July 2002), 604–611.
  * Clegg et al. (2020) Alexander Clegg, Zackory Erickson, Patrick Grady, Greg Turk, Charles C. Kemp, and C. Karen Liu. 2020\.  Learning to Collaborate From Simulation for Robot-Assisted Dressing.  _IEEE Robotics and Automation Letters_ 5, 2 (2020), 2746–2753.
  * Clegg et al. (2018) Alexander Clegg, Wenhao Yu, Jie Tan, C. Karen Liu, and Greg Turk. 2018.  Learning to Dress: Synthesizing Human Dressing Motion via Deep Reinforcement Learning.  _ACM Trans. Graph._ 37, 6, Article 179 (Dec. 2018), 10 pages.
  * Clyde et al. (2017) David Clyde, Joseph Teran, and Rasmus Tamstorf. 2017\.  Modeling and Data-Driven Parameter Estimation for Woven Fabrics. In _Proceedings of the ACM SIGGRAPH / Eurographics Symposium on Computer Animation_ (Los Angeles, California) _(SCA ’17)_. Association for Computing Machinery, New York, NY, USA, Article 17, 11 pages.
  * de Avila Belbute-Peres et al. (2018a) Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and J Zico Kolter. 2018a.  End-to-end differentiable physics for learning and control.  _Advances in neural information processing systems_ 31 (2018), 7178–7189.
  * de Avila Belbute-Peres et al. (2018b) Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and Zico Kolter. 2018b.  End-to-End Differentiable Physics for Learning and Control. In _Advances in Neural Information Processing Systems_. 7178–7189.
  * Degrave et al. (2019) Jonas Degrave, Michiel Hermans, Joni Dambre, and Francis wyffels. 2019.  A Differentiable Physics Engine for Deep Learning in Robotics.  _Frontiers in Neurorobotics_ 13 (2019), 6.
  * Du et al. (2021) Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, and Wojciech Matusik. 2021.  DiffPD: Differentiable Projective Dynamics with Contact.  arXiv:2101.05917 [cs.LG]
  * Du et al. (2020) Tao Du, Kui Wu, Andrew Spielberg, Wojciech Matusik, Bo Zhu, and Eftychios Sifakis. 2020\.  Functional Optimization of Fluidic Devices with Differentiable Stokes Flow.  _ACM Trans. Graph._ 39, 6, Article 197 (Dec. 2020), 15 pages.  https://doi.org/10.1145/3414685.3417795
  * Fratarcangeli et al. (2016) Marco Fratarcangeli, Valentina Tibaldo, and Fabio Pellacini. 2016.  Vivace: A Practical Gauss-Seidel Method for Stable Soft Body Dynamics.  _ACM Trans. Graph._ 35, 6, Article 214 (Nov. 2016), 9 pages.
  * Geilinger et al. (2020) Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz Bächer, Bernhard Thomaszewski, and Stelian Coros. 2020.  ADD: analytically differentiable dynamics for multi-body systems with frictional contact.  _ACM Transactions on Graphics (TOG)_ 39, 6 (2020), 1–15.
  * Guan et al. (2012) Peng Guan, Loretta Reiss, David A. Hirshberg, Alexander Weiss, and Michael J. Black. 2012.  DRAPE: DRessing Any PErson.  _ACM Trans. Graph._ 31, 4, Article 35 (July 2012), 10 pages.
  * Guennebaud et al. (2010) Gaël Guennebaud, Benoît Jacob, et al. 2010\.  Eigen v3.  http://eigen.tuxfamily.org.
  * Hahn et al. (2019) David Hahn, Pol Banzet, James M Bern, and Stelian Coros. 2019\.  Real2sim: Visco-elastic parameter estimation from dynamic motion.  _ACM Transactions on Graphics (TOG)_ 38, 6 (2019), 1–13.
  * Harmon et al. (2008) David Harmon, Etienne Vouga, Rasmus Tamstorf, and Eitan Grinspun. 2008.  Robust Treatment of Simultaneous Collisions. In _ACM SIGGRAPH 2008 Papers_ (Los Angeles, California) _(SIGGRAPH ’08)_. Association for Computing Machinery, New York, NY, USA, Article 23, 4 pages.
  * Holl et al. (2020) Philipp Holl, Nils Thuerey, and Vladlen Koltun. 2020\.  Learning to Control PDEs with Differentiable Physics. In _International Conference on Learning Representations_.
  * Hu et al. (2020) Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan Ragan-Kelley, and Frédo Durand. 2020.  DiffTaichi: Differentiable Programming for Physical Simulation.  _International Conference on Learning Representations_ (2020).
  * Hu et al. (2019) Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua B Tenenbaum, William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik. 2019\.  Chainqueen: A real-time differentiable physical simulator for soft robotics. In _2019 International conference on robotics and automation (ICRA)_. IEEE, 6265–6271.
  * Li et al. (2020b) 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 Dynamic Matrix Assembly and Pipelined Implicit Integrators.  _ACM Trans. Graph._ 39, 6, Article 180 (Nov. 2020), 15 pages.
  * Li et al. (2018a) Jie Li, Gilles Daviet, Rahul Narain, Florence Bertails-Descoubes, Matthew Overby, George E. Brown, and Laurence Boissieux. 2018a.  An Implicit Frictional Contact Solver for Adaptive Cloth Simulation.  _ACM Trans. Graph._ 37, 4, Article 52 (July 2018), 15 pages.
  * Li et al. (2020a) Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020a.  Incremental Potential Contact: Intersection-and Inversion-Free, Large-Deformation Dynamics.  _ACM Trans. Graph._ 39, 4, Article 49 (July 2020), 20 pages.
  * Li et al. (2021) Minchen Li, Danny M. Kaufman, and Chenfanfu Jiang. 2021\.  Codimensional Incremental Potential Contact.  _ACM Trans. Graph. (SIGGRAPH)_ 40, 4, Article 170 (2021).
  * Li et al. (2018b) Minchen Li, Alla Sheffer, Eitan Grinspun, and Nicholas Vining. 2018b.  Foldsketch: Enriching Garments with Physically Reproducible Folds.  _ACM Trans. Graph._ 37, 4, Article 133 (July 2018), 13 pages.
  * Li et al. (2019) Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua Tenenbaum, and Antonio Torralba. 2019.  Learning Particle Dynamics for Manipulating Rigid Bodies, Deformable Objects, and Fluids. In _International Conference on Learning Representations_.
  * Liang et al. (2019) Junbang Liang, Ming Lin, and Vladlen Koltun. 2019.  Differentiable Cloth Simulation for Inverse Problems. In _Advances in Neural Information Processing Systems_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc.  https://proceedings.neurips.cc/paper/2019/file/28f0b864598a1291557bed248a998d4e-Paper.pdf
  * Liu et al. (2013) Tiantian Liu, Adam W. Bargteil, James F. O’Brien, and Ladislav Kavan. 2013.  Fast Simulation of Mass-Spring Systems.  _ACM Trans. Graph._ 32, 6, Article 214 (Nov. 2013), 7 pages.
  * Liu et al. (2017) Tiantian Liu, Sofien Bouaziz, and Ladislav Kavan. 2017\.  Quasi-newton methods for real-time simulation of hyperelastic materials.  _Acm Transactions on Graphics (TOG)_ 36, 3 (2017), 1–16.
  * Ly et al. (2018) Mickaël Ly, Romain Casati, Florence Bertails-Descoubes, Mélina Skouras, and Laurence Boissieux. 2018.  Inverse Elastic Shell Design with Contact and Friction.  _ACM Trans. Graph._ 37, 6, Article 201 (Dec. 2018), 16 pages.
  * Ly et al. (2020) Mickaël Ly, Jean Jouve, Laurence Boissieux, and Florence Bertails-Descoubes. 2020.  Projective Dynamics with Dry Frictional Contact.  _ACM Trans. Graph._ 39, 4, Article 57 (July 2020), 8 pages.
  * Macklin et al. (2020) M. Macklin, K. Erleben, M. Müller, N. Chentanez, S. Jeschke, and T. Y. Kim. 2020\.  Primal/Dual Descent Methods for Dynamics. In _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Virtual Event, Canada) _(SCA ’20)_. Eurographics Association, Goslar, DEU, Article 9, 12 pages.
  * Macklin et al. (2016) Miles Macklin, Matthias Müller, and Nuttapong Chentanez. 2016\.  XPBD: Position-Based Simulation of Compliant Constrained Dynamics. In _Proceedings of the 9th International Conference on Motion in Games_ (Burlingame, California) _(MIG ’16)_. Association for Computing Machinery, New York, NY, USA, 49–54.
  * Martin et al. (2011) Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011.  Example-Based Elastic Materials.  _ACM Trans. Graph._ 30, 4, Article 72 (July 2011), 8 pages.  https://doi.org/10.1145/2010324.1964967
  * McNamara et al. (2004) Antoine McNamara, Adrien Treuille, Zoran Popović, and Jos Stam. 2004.  Fluid Control Using the Adjoint Method.  _ACM Trans. Graph._ 23, 3 (Aug. 2004), 449–456.
  * Miguel et al. (2012) E. Miguel, D. Bradley, B. Thomaszewski, B. Bickel, W. Matusik, M. A. Otaduy, and S. Marschner. 2012.  Data-Driven Estimation of Cloth Simulation Models.  _Comput. Graph. Forum_ 31, 2pt2 (May 2012), 519–528.
  * Miguel et al. (2013) Eder Miguel, Rasmus Tamstorf, Derek Bradley, Sara C. Schvartzman, Bernhard Thomaszewski, Bernd Bickel, Wojciech Matusik, Steve Marschner, and Miguel A. Otaduy. 2013\.  Modeling and Estimation of Internal Friction in Cloth.  _ACM Trans. Graph._ 32, 6, Article 212 (Nov. 2013), 10 pages.
  * Montes et al. (2020) Juan Montes, Bernhard Thomaszewski, Sudhir Mudur, and Tiberiu Popa. 2020.  Computational Design of Skintight Clothing.  _ACM Trans. Graph._ 39, 4, Article 105 (July 2020), 12 pages.
  * Müller et al. (2007) Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007.  Position Based Dynamics.  _J. Vis. Comun. Image Represent._ 18, 2 (April 2007), 109–118.
  * Narain et al. (2012) Rahul Narain, Armin Samii, and James F. O’Brien. 2012\.  Adaptive Anisotropic Remeshing for Cloth Simulation.  _ACM Trans. Graph._ 31, 6, Article 152 (Nov. 2012), 10 pages.
  * Otaduy et al. (2009) Miguel A. Otaduy, Rasmus Tamstorf, Denis Steinemann, and Markus Gross. 2009.  Implicit Contact Handling for Deformable Objects.  _Computer Graphics Forum_ 28, 2 (2009), 559–568.
  * Overby et al. (2017) Matthew Overby, George E. Brown, Jie Li, and Rahul Narain. 2017\.  ADMM $\supseteq$ Projective Dynamics: Fast Simulation of Hyperelastic Models with Dynamic Constraints.  _IEEE Transactions on Visualization and Computer Graphics_ 23, 10 (Oct 2017), 2222–2234.
  * Popović et al. (2003) Jovan Popović, Steven Seitz, and Michael Erdmann. 2003\.  Motion Sketching for Control of Rigid-Body Simulations.  _ACM Trans. Graph._ 22, 4 (Oct. 2003), 1034–1054.
  * Provot (1997) Xavier Provot. 1997\.  Collision and self-collision handling in cloth model dedicated to design garments. In _Computer Animation and Simulation ’97_ , Daniel Thalmann and Michiel van de Panne (Eds.). Springer Vienna, Vienna, 177–189.
  * Qiao et al. (2020) Yi-Ling Qiao, Junbang Liang, Vladlen Koltun, and Ming Lin. 2020\.  Scalable Differentiable Physics for Learning and Control. In _International Conference on Machine Learning_.
  * Sanchez-Gonzalez et al. (2020) Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. 2020.  Learning to Simulate Complex Physics with Graph Networks. In _International Conference on Machine Learning_.
  * Schenck and Fox (2018) Connor Schenck and Dieter Fox. 2018.  SPNets: Differentiable Fluid Dynamics for Deep Neural Networks.  _Conference on Robot Learning (CoRL)_ (2018).
  * Schulman et al. (2017) John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov. 2017.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ (2017).
  * Stern and Desbrun (2006) Ari Stern and Mathieu Desbrun. 2006.  Discrete Geometric Mechanics for Variational Time Integrators. In _ACM SIGGRAPH 2006 Courses_ (Boston, Massachusetts) _(SIGGRAPH ’06)_. Association for Computing Machinery, New York, NY, USA, 75–80.
  * Stern and Grinspun (2009) Ari Stern and Eitan Grinspun. 2009.  Implicit-Explicit Variational Integration of Highly Oscillatory Problems.  _Multiscale Modeling & Simulation_ 7, 4 (2009), 1779–1794.
  * Tamstorf et al. (2015) Rasmus Tamstorf, Toby Jones, and Stephen F. McCormick. 2015\.  Smoothed Aggregation Multigrid for Cloth Simulation.  _ACM Trans. Graph._ 34, 6, Article 245 (Oct. 2015), 13 pages.
  * Terzopoulos et al. (1987) Demetri Terzopoulos, John Platt, Alan Barr, and Kurt Fleischer. 1987\.  Elastically Deformable Models.  _SIGGRAPH Comput. Graph._ 21, 4 (Aug. 1987), 205–214.
  * Toussaint et al. (2018) Marc Toussaint, Kelsey Allen, Kevin Smith, and Joshua Tenenbaum. 2018.  Differentiable Physics and Stable Modes for Tool-Use and Manipulation Planning. In _Robotics: Science and Systems_ , Vol. 2.
  * Treuille et al. (2003) Adrien Treuille, Antoine McNamara, Zoran Popović, and Jos Stam. 2003.  Keyframe Control of Smoke Simulations.  _ACM Trans. Graph._ 22, 3 (July 2003), 716–723.
  * Umetani et al. (2011) Nobuyuki Umetani, Danny M. Kaufman, Takeo Igarashi, and Eitan Grinspun. 2011.  Sensitive Couture for Interactive Garment Modeling and Editing.  _ACM Trans. Graph._ 30, 4, Article 90 (July 2011), 12 pages.
  * Wang (2015) Huamin Wang. 2015\.  A Chebyshev Semi-Iterative Approach for Accelerating Projective and Position-Based Dynamics.  _ACM Trans. Graph._ 34, 6, Article 246 (Oct. 2015), 9 pages.
  * Wang (2018) Huamin Wang. 2018\.  Rule-Free Sewing Pattern Adjustment with Precision and Efficiency.  _ACM Trans. Graph._ 37, 4, Article 53 (July 2018), 13 pages.
  * Wang et al. (2011) Huamin Wang, James F. O’Brien, and Ravi Ramamoorthi. 2011\.  Data-Driven Elastic Models for Cloth: Modeling and Measurement. In _ACM SIGGRAPH 2011 Papers_ (Vancouver, British Columbia, Canada) _(SIGGRAPH ’11)_. Association for Computing Machinery, New York, NY, USA, Article 71, 12 pages.
  * Wang and Yang (2016) Huamin Wang and Yin Yang. 2016.  Descent Methods for Elastic Body Simulation on the GPU.  _ACM Trans. Graph._ 35, 6, Article 212 (Nov. 2016), 10 pages.
  * Wang et al. (2018) Zhendong Wang, Longhua Wu, Marco Fratarcangeli, Min Tang, and Huamin Wang. 2018.  Parallel Multigrid for Nonlinear Cloth Simulation.  _Computer Graphics Forum_ 37, 7 (2018), 131–141.
  * White et al. (2007) Ryan White, Keenan Crane, and David Forsyth. 2007.  Capturing and Animating Occluded Cloth. In _ACM Transactions on Graphics (SIGGRAPH)_.
  * Wojtan et al. (2006) Chris Wojtan, Peter Mucha, and Greg Turk. 2006.  Keyframe Control of Complex Particle Systems Using the Adjoint Method. In _Proceedings of the 2006 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Vienna, Austria) _(SCA ’06)_. Eurographics Association, Goslar, DEU, 15–23.
  * Xian et al. (2019) Zangyueyang Xian, Xin Tong, and Tiantian Liu. 2019\.  A Scalable Galerkin Multigrid Method for Real-Time Simulation of Deformable Objects.  _ACM Trans. Graph._ 38, 6, Article 162 (Nov. 2019), 13 pages.