# Vertex Block Descent

Anka He Chen  0000-0002-5819-3453<EMAIL_ADDRESS>University of UtahSalt
Lake CityUTUSA ,  Ziheng Liu  ****-****-****-**** NA University of UtahSalt
Lake CityUTUSA ,  Yin Yang  ****-****-****-**** NA University of UtahSalt Lake
CityUTUSA and  Cem Yuksel  0000-0002-0122-4159<EMAIL_ADDRESS>University of
Utah & RobloxSalt Lake CityUTUSA

###### Abstract.

We introduce vertex block descent, a block coordinate descent solution for the
variational form of implicit Euler through vertex-level Gauss-Seidel
iterations. It operates with local vertex position updates that achieve
reductions in global variational energy with maximized parallelism. This forms
a physics solver that can achieve numerical convergence with unconditional
stability and exceptional computation performance. It can also fit in a given
computation budget by simply limiting the iteration count while maintaining
its stability and superior convergence rate.

We present and evaluate our method in the context of elastic body dynamics,
providing details of all essential components. Then, we discuss how it can be
used for other simulation problems, including particle-based simulations and
rigid bodies.

physics-based simulation, elastic body, rigid body, time integration

††ccs: Computing methodologies Physical simulation††ccs: Computing
methodologies Collision detection

Figure 1. Example simulation results using our solver, both of those methods
involve more than 100 million DoFs and 1 million active collisions.

##  1\. Introduction

Physics-based simulation is the cornerstone of most graphics applications and
the demands from simulation systems to deliver improved stability, accelerated
computational performance, and enhanced visual realism are ever-growing.
Particularly in real-time graphics applications, the stability and performance
requirements are so strict that realism can sometimes be begrudgingly
considered of secondary importance.

Notwithstanding the substantial amount of research and groundbreaking
discoveries made on physics solvers over the past decades, existing methods
still leave some things to be desired. They either deliver high-quality
results, but fail to meet the computational demands of many applications or
fit in a given computation time by sacrificing quality. Stability, on the
other hand, is always a challenge, particularly with strict computation
budgets.

In this paper, we introduce _vertex block descent (VBD)_ , a physics solver
that offers unconditional stability, superior computational performance than
prior methods, and the ability to achieve numerical convergence. Though our
method is a general solution that can be used for a variety of simulation
problems, we present and evaluate it in the context of elastic body dynamics.
Then, we briefly discuss how our method can be applied to some other
simulation systems, including particle-based simulations and rigid bodies.

Our VBD method is based on block coordinate descent that performs vertex-based
Gauss-Seidel iterations to solve the variational form of implicit Euler. For
elastic body dynamics, each iteration runs a loop over the mesh vertices,
adjusting the position of a single vertex at a time, temporarily fixing all
others. This offers maximized parallelism when coupled with vertex-based mesh
coloring, which can achieve an order of magnitude fewer colors (i.e.
serialized workloads) as compared to element-based parallelization. Our local
position-based updates can ensure that the variational energy is always
reduced. Therefore, our method maintains its stability even with a single
iteration per time step and large time steps, operating with unconverged
solutions containing a large amount of residual. With more iterations, it
converges faster than its alternatives. Thus, it can more easily fit in a
given computation budget, while maintaining stability with improved
convergence.

We present all essential components of using VBD for elastic body dynamics,
including formulations for damping, constraints, collisions, and friction. We
also introduce a simple initialization scheme to warm-start the optimization
and improve convergence. In addition, we discuss momentum-based acceleration
techniques and parallelization in the presence of dynamic collisions. Our
evaluation includes large simulations (Figure 1) and stress tests (e.g. Figure
2 and 3) that demonstrate VBD’s performance, scalability, and stability.

Figure 2. Twisting two beams together, totaling 97 thousand vertices and 266
thousand tetrahedra, demonstrating complex frictional contact and buckling.

##  2\. Related Work

There is a large body of work on physics-based simulation in computer
graphics. Here we only discuss the recent and the most relevant work to our
method.

Implicit time integrators are widely accepted as the primary methods for
simulating elastic bodies in computer graphics due to their exceptional
stability, especially when addressing stiff problems. Among these options,
backward Euler (Baraff and Witkin, 1998; Hirota et al., 2001; Volino and
Magnenat-Thalmann, 2001; Martin et al., 2011) is the most commonly utilized
method, though other approaches like implicit Newmark (Kane et al., 2000;
Bridson et al., 2002, 2005), BDF2 (Choi and Ko, 2005; Hauth and Etzmuss,
2001), and implicit-explicit (Eberhardt et al., 2000; Stern and Grinspun,
2009) have also been explored. Backward Euler is often approximated as a
single Newton step, solving a linear system of equations (Baraff and Witkin,
1998). Line search can be applied to improve robustness (Hirota et al., 2001).
Preconditioning (Ascher and Boxerman, 2003) or positive-definite projection
(Teran et al., 2005) can be used to improve convergence. To circumvent a full
linear solve for every Newton step, Cholesky factorization (Hecht et al.,
2012) emerges as a viable strategy. Techniques like multi-resolution (Grinspun
et al., 2002; Capell et al., 2002) or multigrid (Bolz et al., 2003; Tamstorf
et al., 2015; Xian et al., 2019; Wang et al., 2020; Zhu et al., 2010) solvers
project finer details onto a coarser grid with fewer degrees of freedom,
effectively reducing the computational cost of the linear system solver. Our
solution linearizes the forces locally, avoiding the global linear system, and
it converges to the same result as backward Euler with multiple Newton steps.

Additionally, stiffness warping (Müller et al., 2002) repurposes the stiffness
matrix from the rest shape to handle significant rotational deformations.
Using a quasi-Newton method (Liu et al., 2017) with an approximate Hessian can
be far more cost-effective for computing its inverse with prefactoring than
the actual Hessian. Example-based dynamics has also been explored (Chao et
al., 2010; Müller et al., 2005; Martin et al., 2011), where the deformation
energy is defined based on the nearest point in the example space. Projective
Dynamics (Bouaziz et al., 2023) represents deformation energy via a series of
constraints that can be solved independently, then synchronized either through
a prefactorized global step to accelerate convergence.

The relation between dynamics, energy, and minimization has been leveraged in
variational integrators (Kane et al., 2000, 1999; Simo et al., 1992; Lew et
al., 2004; Kharevych et al., 2006). Reformulating backward Euler to an
optimization problem combined with optimization techniques to allows the usage
of large steps (Martin et al., 2011; Gast et al., 2015). Domain-decomposed
optimization for solving the nonlinear problems of implicit numerical time
integration can accelerate convergence Li et al. (2019). Yet, the optimization
formulation has its drawbacks, notably the varying problem formulation and
initial positions in each step. Consequently, achieving consistent convergence
within a fixed time budget becomes challenging. Real-time simulators
compromise by accepting partially converged results, which visually resemble
the final solution. Unfortunately, this compromise can lead to significant
visual artifacts mainly due to retained residuals from earlier steps, which
can accumulate across frames and can jeopardize the solver’s stability. Our
method, as illustrated in Figure 16, demonstrates exceptional stability even
when retaining a substantial amount of residual with a limited number of
iterations.

Figure 3. Stress tests that begin simulations under extreme deformations:
(top) a perfectly flattened armadillo model with 15 thousand vertices and 50
thousand tetrahedra, and (bottom) a Utah teapot model with 2 thousand vertices
and 8.5 thousand tetrahedra, deformed by randomly placing its vertices onto
the surface of a sphere. Both models quickly recover to their original shape
shortly after the simulation starts. Both simulations use accelerated
iterations with $\rho=0.95$.

Position-based dynamics methods (PBD (Müller et al., 2007) and XPBD (Macklin
et al., 2016)) convert forces into (soft) constraints and directly update the
positions with Gauss-Seidel iterations operating on one constraint at a time.
These position-based updates result in exceptional stability, which is often
exploited by limiting the number of iterations to fix the computation cost, an
effective strategy for real-time simulations. Akin to PBD, our method also
works with position updates, but operates directly using the force
formulations without converting them to constraints. Parallelization with XPBD
is achieved by graph coloring the constraints (i.e. the dual graph)
(Fratarcangeli and Pellacini, 2015; Fratarcangeli et al., 2016; Ton-That et
al., 2023). However, this dual graph contains multiple times more connections
(depending on the types of constraints) than the original graph of vertices,
severely limiting the level of parallelism. In comparison, our method is
parallelized by coloring the original graph, which leads to much fewer colors
(i.e. computation groups that must be processed sequentially) and thereby
better parallelism. More importantly, the approximations in XPBD’s formulation
introduce errors that make it diverge from the solution of implicit Euler and
can degrade realism, particularly with large time steps and limited iteration
count, which are common in practice. In addition, XPBD particularly struggles
with high mass ratios. Our method, on the other hand, has none of these
problems.

In recent years, a growing effort has been placed on accelerating simulations
using GPUs (Wang, 2015; Huthwaite, 2014; Wang and Yang, 2016; Lan et al.,
2022; Macklin et al., 2020; Li et al., 2020b, 2023). Among them, the first-
order descent methods (Wang and Yang, 2016; Macklin et al., 2020) have gained
popularity due to their excellent parallelism. These methods employ a Jacobi-
style preconditioned first-order descent on the backward Euler minimization
formulation, enabling full vertex-level parallelization. However, Jacobi-style
iterations typically converge substantially slower than Gauss-Seidel
iterations. Also, such methods necessitate a line search to avoid overshooting
and ensure stability.

Our method can be categorized as a coordinate descent method for optimization
(Wright, 2015). In graphics, this technique has been used for geometric
processing (Naitsat et al., 2020) and simulation with a barrier function (Lan
et al., 2023). Recently, Lan et al. (2023) proposed a hybrid scheme where
Gauss-Sediel and Jacobi iterations are combined at each parallel call. In
comparison, our method uses blocks of coordinates based on vertices instead of
blocks based on elements, which results in much better parallelism and smaller
local linear systems to solve, leading to faster convergence.

##  3\. Vertex Block Descent for Elastic Bodies

Our vertex block descent method is essentially a solver for optimization
problems. Therefore, it can be applied to various simulation problems,
particularly if they can formulated as an optimization problem. In this
section, we explain our method in the context of elastic body dynamics for
objects represented by a set of vertices that carry mass and a set of force
elements and constraints that act on them. Generalization of our method to
other example simulation problems is discussed later in Section 5.

We begin with deriving our global optimization method that splits the
simulated system into vertex-level local systems (Section 3.1). Then, we
discuss the methods we use for solving the local systems (Section 3.2) and
describe how we incorporate damping (Section 3.3) and constraints (Section
3.4). We present our collision formulation (Section 3.5) and friction
formulation (Section 3.6). After explaining our methods for warm-starting our
optimization to improve convergence (Section 3.7), we describe how to
incorporate momentum-based acceleration to improve convergence (Section 3.8).
Finally, we discuss the improved parallelization that our method provides
along with methods for efficiently parallelizing dynamically-introduced force
elements due to varying collision events (Section 3.9).

###  3.1. Global Optimization

Our vertex block descent method is formulated based on the optimization form
of implicit Euler time integration (Martin et al., 2011; Gast et al., 2015).
Given a system with $N$ vertices, we represent the state of our simulation at
step $t$ as ($\vec{x}^{t},\vec{v}^{t})$, where $\vec{x}^{t}\in\mathbb{R}^{3N}$
and $\vec{v}^{t}\in\mathbb{R}^{3N}$ are the concatenated position and velocity
vectors, respectively. The resulting optimization problem can be written as

(1) $\displaystyle\vec{x}^{t+1}$
$\displaystyle=\underset{\vec{x}}{\operatorname{argmin}}\;G(\vec{x})\;,$

for evaluating the positions at the end of the time step $\vec{x}^{t+1}$ that
minimizes the _variational energy_ $G(\vec{x})$ in the form of

(2) $\displaystyle G(\vec{x})$
$\displaystyle=\frac{1}{2h^{2}}\lVert\vec{x}-\vec{y}\rVert_{M}^{2}+E(\vec{x})\;.$

The first term here is the _inertia potential_ , which contains the time step
size $h$, the mass-weighted norm $\lVert\cdot\rVert_{M}$, and

(3) $\displaystyle\vec{y}$
$\displaystyle=\vec{x}^{t}+h\vec{v}^{t}+h^{2}\vec{a}_{\text{ext}}$

with a fixed external acceleration $\vec{a}_{\text{ext}}$, such as gravity.
The second term $E(\vec{x})$ is the _total potential energy_ evaluated at
$\vec{x}$.

We propose an optimization technique that falls under the category of
coordinate descent methods to efficiently minimize this energy $G$. If we only
modify a single vertex at a time, fixing all other vertices, the part of the
energy term $E(\mathbf{x})$ that is affected only includes the set of force
elements $\mathcal{F}_{i}$ that are acting on (or using the position of)
vertex $i$. Thus, we define the _local variational energy_ $G_{i}$ around
vertex $i$ as

(4) $\displaystyle G_{i}(\vec{x})$
$\displaystyle=\frac{m_{i}}{2h^{2}}\lVert\vec{x}_{i}-\vec{y}_{i}\rVert^{2}+\sum_{j\in\mathcal{F}_{i}}E_{j}(\vec{x})\;,$

where $m_{i}$ is the mass of the vertex and $E_{j}$ is the energy of force
element $j$ in $\mathcal{F}_{i}$.

Note that $G$ is not equal to the sum of local variational energies, i.e.
$G(\vec{x})\neq\sum_{i}G_{i}(\vec{x})$, simply because the force elements
appear multiple times in this sum (once for each of its vertices). However,
when we modify the position of a single vertex only, the reduction in $G_{i}$
is equal to the resulting reduction in $G$.

Based on this observation, our method operates on a single vertex at a time
and updates its position by minimizing the local energy

(5) $\displaystyle\vec{x}_{i}$
$\displaystyle\leftarrow\underset{\vec{x}_{i}}{\operatorname{argmin}}\;G_{i}(\vec{x})$

and solves the global system using Gauss-Seidel iterations. Each local
minimization for a vertex effectively finds a descent step for $G$ using the
degrees of freedom (DoF) for the vertex as a block of coordinates, hence the
name _vertex block descent_ (VBD). After each iteration, the total reduction
in $G$ is equal to the sum of all reductions in $G_{i}$. In other words, the
energy change of each vertex position adjustment is accumulated to the system
energy. Consequently, if we can make sure that each local energy drops $G_{i}$
when we are adjusting vertex $i$, we can guarantee that we are descending the
system energy $G$.

Thus, our system directly operates on vertex positions and the resulting
velocities are calculated following the implicit Euler formulation

(6)
$\displaystyle\vec{v}^{t+1}=\frac{1}{h}\big{(}\vec{x}^{t+1}-\vec{x}^{t}\big{)}\;.$

###  3.2. Local System Solver

The position updates per vertex (in Equation 5) involve solving a local system
that only depends on the position change of a single vertex, represented by
$G_{i}$. Note that Equation 4 only has 3 DoF, so the cost of evaluating and
inverting its Hessian is much cheaper compared to the global problem in
Equation 2. Therefore, we can fully utilize the second-order information and
use Newton’s method to minimize the localized energy $G_{i}$, which involves
solving the 3D linear system

(7) $\displaystyle\vec{H}_{i}\,\Delta\vec{x}_{i}$
$\displaystyle=\mathbf{f}_{i}\;,$

where $\Delta\mathbf{x}_{i}$ is the change in position, $\mathbf{f}_{i}$ is
the total force acting on the vertex, calculated using

(8) $\displaystyle\mathbf{f}_{i}$ $\displaystyle=-\frac{\partial
G_{i}(\vec{x})}{\partial\vec{x}_{i}}=-\frac{m_{i}}{h^{2}}\left(\mathbf{x}_{i}-\mathbf{y}_{i}\right)\;-\sum_{j\in\mathcal{F}_{i}}\frac{\partial
E_{j}(\vec{x})}{\partial\vec{x}_{i}}\;,$

and $\vec{H}_{i}\in\mathbb{R}^{3\times 3}$ is the Hessian of $G_{i}$ with
respect to the DoF of vertex $i$, such that

(9) $\displaystyle\vec{H}_{i}$
$\displaystyle=\frac{\partial^{2}G_{i}(\vec{x})}{\partial\vec{x}_{i}\partial\vec{x}_{i}}=\frac{m_{i}}{h^{2}}\mathbf{I}\;+\sum_{j\in\mathcal{F}_{i}}\frac{\partial^{2}E_{j}(\vec{x})}{\partial\vec{x}_{i}\partial\vec{x}_{i}}\;.$

Here, the first term is a diagonal matrix and the second one is the sum of
Hessians of the force elements with respect to vertex $i$. Intuitively, the
solution of this linear system is the extreme point of the quadratic
approximation for the localized energy $G_{i}$. By reducing $G_{i}$ with each
iteration, we can guarantee a reduction in $G$, thereby iteratively solving
the global system in Equation 2.

We can analytically solve this system using
$\Delta\mathbf{x}_{i}=\mathbf{H}_{i}^{-1}\mathbf{f}_{i}$. For such a small
system, the analytical solver is very efficient and stable. We found it to be
faster than solvers based on Conjugate Gradient or LU/QR decomposition.
Another advantage of the analytical solver is that it does not require the
Hessian to be positive-definite. Of course, when the Hessian is not positive-
definite, the direction given by Equation 7 may not be a descent direction.
Nevertheless, we opt for this direction regardless, recognizing that even when
$\vec{H}_{i}$ is not positive-definite, solving Equation 7 still brings us
towards the extremum of the quadratic approximation for the localized energy
$G_{i}$. This solution is close to where the gradient of the inertia and the
potential terms balance out and it is usually a stable state of the system.
Also, the motivation of the variational form of implicit Euler (Equation 2),
is to find a point where ${dG(\vec{x})/d\vec{x}=0}$. Therefore, any extreme
point is a valid solution of implicit Euler, and it does not have to be a
local minimum. In all our experiments, including those specifically designed
stress tests (see Figure 2 and Figure 3), we have consistently observed that
this scheme does not pose any issues concerning system stability or
convergence.

An alternative solution to this is the PSD Hessian projection (Teran et al.,
2005). However, it is exceptionally rare for the Hessian to not be positive-
definite, and the PSD projection process is notably costly due to multiple SVD
decompositions. Engaging in this costly operation to prevent such rare events
seems unjustified, especially considering that these occurrences do not
jeopardize system stability or convergence significantly.

Another challenge with the analytical solver arises when encountering a
(nearly) low-rank Hessian. To address this, we propose a simple solution: if
$|\det(\vec{H}_{i})|\leq\epsilon$ for some small threshold $\epsilon$, we opt
to bypass adjusting this particular vertex for that iteration. Given that its
neighboring vertices are likely to undergo adjustments before the next
iteration, it is improbable that its Hessian will remain low-rank in
subsequent iterations. With this simple solution, in the extreme scenario
where all vertices possess a low-rank Hessian, the system could potentially
become frozen. Yet, it is crucial to note that such a case is highly
improbable, since the Hessian of the inertia potential is always full-rank.
One potential remedy for this would be switching to the modified Conjugate
Gradient solver (Lan et al., 2023) when such a case happens. However, doing
this will add additional branching to the code and can slow down the solver.
Thus, we have not included it in our implementation, but, depending on
specific use cases, there is always the flexibility to opt for the Conjugate
Gradient solver as a backup solution.

The linear system in Equation 7 corresponds to a single Newton step, so it
does not necessarily provide the optimal solution for Equation 5. In fact,
since it is just a second-order approximation of $G_{i}$, it does not even
guarantee a reduction in $G_{i}$. To ensure the descent of energy with this
single step, we can incorporate a backtracking line search along the descent
direction. Note that, unlike global line searches in descent-based simulation
methods (e.g. Wang and Yang (2016)), this line search operates locally. It
specifically verifies the descent of $G_{i}$, which in turn guarantees a
descent in $G$ without the need for evaluating the global system.

Line search avoids over-shooting and, thereby, ensures stability. In practice,
however, the additional computation cost of line search may not be justified.
In our experiments, we found that line search costs an extra 40% computation
time, while not providing any measurable benefits. This is because VBD can
maintain stability even without line search. Therefore, the results we present
in this paper do not include line search, though it is an option available.

###  3.3. Damping

Damping plays a crucial role in simulations. It prevents excessive
oscillations and also enhances system stability. Despite the inherent
numerical damping introduced by the implicit Euler method, providing users
with manual control over damping is highly desirable. To address this, we have
integrated a simplified Rayleigh damping model into our solver (Sifakis and
Barbic, 2012). This process is both straightforward and efficient, as it also
operates locally within the $3\times 3$ system and utilizes the precomputed
stiffness matrix.

Since we are relying on implicit Euler, we can represent the velocity as
position change, using $\vec{v}_{i}=(\vec{x}_{i}-\vec{x}_{i}^{t})/h$. Then, we
can add the damping term to the Hessian in Equation 9, resulting

(10) $\displaystyle\vec{H}_{i}$
$\displaystyle=\frac{m_{i}}{h^{2}}\mathbf{I}\;+\sum_{j\in\mathcal{F}_{i}}\frac{\partial^{2}E_{j}(\vec{x})}{\partial\vec{x}_{i}\partial\vec{x}_{i}}+\left(\sum_{j\in\mathcal{F}_{i}}\frac{k_{d}}{h}\frac{\partial^{2}E_{j}(\vec{x})}{\partial\vec{x}_{i}\partial\vec{x}_{i}}\right)\;,$

where $k_{d}$ is the damping coefficient. Finally, we add the damping force to
Equation 8 using the same damping term, such that

(11) $\displaystyle\mathbf{f}_{i}$
$\displaystyle=-\frac{m_{i}}{h^{2}}\left(\mathbf{x}_{i}-\mathbf{y}_{i}\right)-\\!\\!\sum_{j\in\mathcal{F}_{i}}\frac{\partial
E_{j}(\vec{x})}{\partial\vec{x}_{i}}-\left(\sum_{j\in\mathcal{F}_{i}}\frac{k_{d}}{h}\frac{\partial^{2}E_{j}(\vec{x})}{\partial\vec{x}_{i}\partial\vec{x}_{i}}\right)\left(\vec{x}_{i}-\vec{x}_{i}^{t}\right).$

###  3.4. Constraints

Since our method directly manipulates the position of each vertex, managing a
constraint on a vertex becomes straightforward. Constraints generally fall
into two categories: unilateral ($C(\mathbf{x})\leq 0$) or bilateral
($C(\mathbf{x})=0$). With biliteral constraints, if a vertex position is
directly set to a specific value, we simply skip updating its position.
Otherwise, it is constrained to a (usually linear) subspace. Our approach
involves representing the constrained vertex position using the subspace
basis. This transforms both the vertex position and gradient into an
$L$-dimensional vector, where $L$ is the subspace dimension. Consequently,
handling local steps for constrained vertices involves solving an $L\times L$
system. Regarding unilateral constraints, we allow compromises and define
potential energy to be solved alongside other potentials. This method handles
world box constraints in our simulations.

(a) Edge-edge

(b) Vertex-triangle

Figure 4. Two collisions types: (a) edge-edge can have at most two pairs and
(b) vertex-face can have at most one pair with the same color, since vertices
on the same side of a collision must have different colors.

###  3.5. Collisions

Collisions can be handled by simply introducing a quadratic collision energy
per vertex, based on the penetration depth $d$, such that

(12) $\displaystyle E_{c}(\vec{x})$ $\displaystyle=\frac{1}{2}\,k_{c}\,d^{2}$
with $\displaystyle d$
$\displaystyle=\max\Big{(}0,(\vec{x}_{b}-\vec{x}_{a})\cdot\vec{\hat{n}}\Big{)}\;,$

where $k_{c}$ is the collision stiffness parameter, $\mathbf{x}_{a}$ and
$\mathbf{x}_{b}$ are the two _contact points_ on either side of the collision,
and $\mathbf{\hat{n}}$ is the contact normal.

There are two collision types for triangle meshes (Figure 4):

  * •

Edge-edge collisions use continuous collision detection (CCD).
$\mathbf{x}_{a}$ and $\mathbf{x}_{b}$ correspond to the intersection points on
either edge and the contact normal is the direction between them, i.e.
${\mathbf{\hat{n}}=\mathbf{n}/\|\mathbf{n}\|}$, where
${\mathbf{n}=\mathbf{x}_{b}-\mathbf{x}_{a}}$.

  * •

Vertex-triangle collisions are detected either by CCD or discrete collision
detection (DCD). In this case, $\mathbf{x}_{a}$ is the colliding vertex and
$\mathbf{x}_{b}$ is the corresponding point on either the collision point (for
CCD) or the closest point (for DCD) on the triangle (Chen et al., 2023).
$\mathbf{\hat{n}}$ is the surface normal at $\mathbf{x}_{b}$.

In our implementation, we perform a DCD at the beginning of the time step
using $\mathbf{x}^{t}$ to identify vertices that have already penetrated, and
the rest of the collisions use CCD. We simplify the computation of the
gradient and the Hessian of the collision energy by not differentiating
through $\mathbf{\hat{n}}$, i.e. assuming that $\mathbf{\hat{n}}$ is constant.

Performing collision detection at every iteration using CCD can be expensive
and can easily become the bottleneck. Therefore, in our implementation, we
perform CCD once every $n_{\text{col}}$ iterations. While this has the risk of
missing some collision events, they are likely to be captured via DCD in the
next time step. All detected collisions remain as force elements until the
next collision detection. For vertex-triangle collisions detected with DCD, it
is important to recompute the closest point ($\mathbf{x}_{b}$) before
computing the gradient and Hessian of $E_{c}$.

###  3.6. Friction

To compute friction for collision $c$, we must consider the relative motion of
the contact points defined as

(13)
$\delta\mathbf{x}_{c}=\big{(}\mathbf{x}_{a}-\mathbf{x}_{a}^{t}\big{)}-\big{(}\mathbf{x}_{b}-\mathbf{x}_{b}^{t}\big{)}\;,$

where $\mathbf{x}_{a}^{t}$ and $\mathbf{x}_{b}^{t}$ are the positions of
$\mathbf{x}_{a}$ and $\mathbf{x}_{b}$ at the beginning of the time step.

With this $\delta\mathbf{x}_{c}$, we can use the friction model of
_incremental potential contact_ (IPC) (Li et al., 2020a). First, we project
$\delta\mathbf{x}_{c}$ to the 2D contact tangential space, using a
transformation matrix $\mathbf{T}_{c}\in\mathbb{R}^{3\times 2}$, to evaluate
the tangential relative translation
$\mathbf{u}_{c}=\vec{T}_{c}^{T}\delta\vec{x}_{c}$, where
$\mathbf{T}_{c}=[\mathbf{\hat{t}}~{}\mathbf{\hat{b}}]$ is formed by a tangent
$\mathbf{\hat{t}}$ and binormal $\mathbf{\hat{b}}$ vectors orthogonal to
$\mathbf{\hat{n}}$. The signed magnitude of the collision force applied on
vertex $i$ is $\lambda_{c,i}=\frac{\partial
E_{c}}{\partial\vec{x}_{i}}\cdot\mathbf{\hat{n}}$. Note that the sign of
$\lambda_{c,i}$ is different for vertices on different sides of the collision.
Let $\mu_{c}$ be the friction coefficient. We can then calculate the friction
force using

(14) $\displaystyle\mathbf{f}_{c,i}$
$\displaystyle=-\mu_{c}\,\lambda_{c,i}\,\frac{\partial\delta\vec{x}_{c}}{\partial\vec{x}_{i}}\,\vec{T}_{c}\,f_{1}(\|\vec{u}_{c}\|)\frac{\vec{u}_{c}}{\|\vec{u}_{c}\|}\;\text{,
where}$ (15) $\displaystyle f_{1}(u)$
$\displaystyle=\begin{cases}2\left(\frac{u}{\epsilon_{v}h}\right)-\left(\frac{u}{\epsilon_{v}h}\right)^{2},&u\in(0,h\epsilon_{v})\\\
1,&u\geq h\epsilon_{v}\end{cases}.$

Here, $f_{1}$ serves as a smooth transition function between static and
dynamic friction. When the relative velocity exceeds a small threshold
$\epsilon_{v}$, dynamic friction is applied. Conversely, if the relative
velocity is below this threshold, static friction is applied, scaling between
the range of $[0,1]$.

In our formulation, we need the Hessian of the friction term, which is the
derivative of this function. We approximate the derivative by not
differentiating through $\|\vec{u}_{c}\|$ for a more stable friction force
formulation (Macklin et al., 2020), such that

(16) $\displaystyle\frac{\partial\mathbf{f}_{c,i}}{\partial\mathbf{x}_{i}}$
$\displaystyle\approx-\mu_{c}\,\lambda_{c,i}\,\frac{\partial\delta\mathbf{x}_{c}}{\partial\mathbf{x}_{i}}\,\vec{T}_{c}\,\frac{f_{1}(\|\vec{u}_{c}\|)}{\|\vec{u}_{c}\|}\vec{T}_{c}^{T}\,\left(\frac{\partial\delta\mathbf{x}_{c}}{\partial\mathbf{x}_{i}}\right)^{T}\;.$

Without this approximation, PSD projection and line search might be needed to
ensure stability. Finally, all friction forces $\mathbf{f}_{c,i}$ and their
derivatives $\partial\mathbf{f}_{c,i}/\partial\mathbf{x}_{i}$ are added to
$\mathbf{f}_{i}$ and $\mathbf{H}_{i}$, respectively.

###  3.7. Initialization

Since our method is an iterative solver, we begin with an _initial guess_ for
$\vec{x}$. The closer this initial guess is to the resulting $\vec{x}^{t+1}$,
the fewer number of iterations we would need to converge. Regardless of the
initialization strategy for determining this initial guess, because we are
solving an optimization problem, the system is expected to converge to the
same result, though it may require a different number of iterations.

Providing a good initial guess (one that is close to $\vec{x}^{t+1}$) is
particularly important for applications with a limited computation budget,
e.g. using a fixed number of iterations. In such applications, the initial
guess can strongly impact the remaining residual at the end of the final
iteration.

(a) Previous
position

(b) Inertia

(c) Inertia & acceleration

(d) Adaptive

(e) Reference

Figure 5. Different initialization options for a swinging elastic pendulum
dropped from the same height (blue line) simulated with our method using only
20 iterations per frame, showing the same frame of the simulation. Notice that
initializing using (a) previous position and (b) inertia fail to properly move
under gravity, while (c) inertia and acceleration leads to accessive
stretching (red line) when VBD does not run to convergence. (d) Our adaptive
solution closely matches (e) the reference generated by fully converged
Newton’s method (green line).

We begin with considering three simple options for initialization:

(a) Previous position: $\displaystyle\vec{x}$ $\displaystyle=\vec{x}^{t}$ (b)
Inertia: $\displaystyle\vec{x}$ $\displaystyle=\vec{x}^{t}+h\mathbf{v}^{t}$
(c) Inertia and acceleration: $\displaystyle\vec{x}$
$\displaystyle=\mathbf{x}^{t}+h\mathbf{v}^{t}+h^{2}\mathbf{a}_{\text{ext}}=\vec{y}$

Option (a) struggles with substantially stiff materials. As we adjust each
vertex separately, assuming the others remain fixed, our method must lean on
the inertia potential’s gradient (i.e.
${m_{i}(\mathbf{x}_{i}-\mathbf{y}_{i})/2h^{2}}$) to march toward the local
minimum. With stiff materials, this method results in notably slower
convergence rates due to the inertia potential being considerably less stiff.
Consequently, it encounters challenges in simulating scenarios that resemble
free fall, like a swinging elastic pendulum at its maximum height, as shown in
5(a).

Option (b) allows the system to start with the inertia of the previous step,
which helps, but terminating the iterations prior to convergence can again
result in local material stiffness overpowering external acceleration, as
shown in 5(b).

Option (c) is similar to the initialization of position-based dynamics, and
performs notably better as it effectively preserves inertia and properly
reacts to external acceleration. However, with a limited number of iterations,
materials behave softer than they should, often resulting in overstretching or
collapsing under gravity. An example of this can be seen in 5(c), where the
pendulum extends more than it should. Most notably, it struggles with steady
objects at rest in contact, consistently initializing them into a penetrating
state, as if they are in free fall. In such cases, the contact forces must
entirely undo the position change of initialization during the iterations.
This not only creates unnecessary computation, but also places considerable
strain on the accuracy of the collision detection and handling methods
(including friction). Therefore, properly simulating objects that are stacked
on top of each other becomes a major challenge with this initialization
option.

1.0
0.8
0.6
0.4
0.2
0.0
gravity initialization factor ($\tilde{a}$)

Figure 6. The ratio of gravity $\tilde{a}$ used with adaptive initialization
during the swinging of an elastic pendulum. The model is a single-piece
tetrahedral mesh. It automatically distinguishes vertices in approximate free-
fall (red) and those where elasticity counteracts gravity (blue).

We propose an adaptive initialization scheme that combines options (b) and
(c), taking advantage of the freedom that VBD provides in the choice of
initialization. This scheme uses

(d) Adaptive: $\displaystyle\vec{x}$
$\displaystyle=\mathbf{x}^{t}+h\mathbf{v}^{t}+h^{2}\mathbf{\tilde{a}}$

replacing the external acceleration $\mathbf{a}_{\text{ext}}$ in (c) with an
estimated acceleration term $\mathbf{\tilde{a}}$, determined by exploiting the
typical similarity between two consecutive time steps. We begin with the
acceleration of the previous frame
${\vec{a}^{t}=(\mathbf{v}^{t}-\mathbf{v}^{t-1})/h}$ and compute its component
$a^{t}_{\text{ext}}$ along the external acceleration direction
${\mathbf{\hat{a}}_{\text{ext}}=\mathbf{a}_{\text{ext}}/\|\mathbf{a}_{\text{ext}}\|}$,
such that
${a^{t}_{\text{ext}}=\mathbf{a}^{t}\cdot\mathbf{\hat{a}}_{\text{ext}}}$. Then,
we simply make sure that the estimated acceleration does not exceed the
external acceleration, using

(17) $\displaystyle\mathbf{\tilde{a}}$
$\displaystyle=\tilde{a}\,\mathbf{a}_{\text{ext}}\;,$ where
$\displaystyle\tilde{a}$ $\displaystyle=\begin{cases}1\;,&\text{if
$a^{t}_{\text{ext}}>\|\mathbf{a}_{\text{ext}}\|$}\\\ 0\;,&\text{if
$a^{t}_{\text{ext}}<0$}\\\
a^{t}_{\text{ext}}/\|\mathbf{a}_{\text{ext}}\|\;,&\text{otherwise.}\end{cases}$

This adaptive approach includes external acceleration in the initialization
when the motion of a vertex resembles free fall. When an object is stationary,
however, as in rest-in-contact, it maintains the previous position in the
initialization, preventing undesired penetrations before the first iteration.
It also successfully avoids excessive stretching, as can be seen in 5(d).
Visualizations of different $\tilde{a}$ values in this simulation are shown in
Figure 6. In short, our adaptive initialization is a simple but effective
strategy and it is possible because VBD does not dictate a particular
initialization (unlike XPBD, for example).

###  3.8. Accelerated Iterations

We use the Chebyshev semi-iterative approach (Wang, 2015) to improve the
convergence of our method, though other momentum-based acceleration
techniques, such as the Nesterov’s method (Golub and Van Loan, 2013) can be
applied as well. The Chebyshev method iteratively computes an acceleration
ratio based on the approximation of the system’s spectral radius. Instead of
directly using the output vertex positions of Gauss-Seidel
${\mathbf{\bar{x}}}^{(n)}$ after iteration $n$, it recomputes the positions at
the end of the iteration using

(18) $\displaystyle\mathbf{x}^{(n)}$
$\displaystyle=\omega_{n}({\mathbf{\bar{x}}}^{(n)}-\mathbf{x}^{(n-2)})+\mathbf{x}^{(n-2)}\;,$

where $\omega_{n}$ is the acceleration ratio that changes at each iteration as

(19) $\displaystyle\omega_{n}$
$\displaystyle=\frac{4}{4-\rho^{2}\omega_{n-1}}$ with
$\displaystyle\omega_{1}$ $\displaystyle=1$ and $\displaystyle\omega_{2}$
$\displaystyle=\frac{2}{2-\rho^{2}}\;,$

where $\rho\in(0,1)$ is the estimated spectral radius, which can be set
manually, or automatically tuned using the technique introduced in (Wang and
Yang, 2016). Note that this position recomputation procedure is performed
globally after each Gauss-Seidel iteration is completed, not after each local
solver step.

(a) Not accelerated

(b) Accelerated

(c) Reference

Figure 7.  Demonstrating the accelerator’s impact in a collision-intensive
scene: a squishy ball (230K vertices, 700K tetrahedra) dropping and bouncing.
Both (a) and (b) use $h=1/120\text{ seconds}$ with a constant number of $120$
iterations per time step, taking $30\text{ms}$ of average computation time per
frame. (a) Without acceleration 120 iterations appear to be insufficient. (b)
Our acceleration scheme ($\rho=0.95$), skipping colliding vertices, manages to
resolve complex collisions, notably enhancing elasticity convergence for much
stiffer outcomes, closely matching (c) the reference computed using 2000
iterations.

This accelerator was originally developed for solving linear systems, assuming
the energy to be smooth and (nearly) quadratic. Elastic energy generally
fulfills these criteria. However, collision energy tends to be discontinuous
and highly stiff, making the use of an accelerator in collision-intensive
scenes prone to overshot and compromise the system’s stability. To address
this, we propose a simple yet highly effective solution for accelerating
scenes with collisions: skipping the accelerations for actively colliding
vertices. Note that the acceleration must be a continuous process. If a vertex
is detected colliding at a certain iteration, we will skip the acceleration
for it in all the following iterations in the same step, regardless of whether
the collision has been resolved. This approach has minimal impact on the
convergence of elasticity, since typically only a small fraction of vertices
are in collision. Also, for those colliding vertices, the elasticity is
usually overpowered by the collision forces. Thus, this solution maintains the
stability of the system, while effectively accelerating the convergence of
elasticity, as shown in Figure 7.

###  3.9. Parallelization

Gauss-Seidel-type iterative methods are often parallelized using graph
coloring by determining groups (i.e. colors) that can be handled in parallel
without impacting the sequential nature of the Gauss-Seidel loop. Obviously,
the same can be applied to VBD by simply coloring vertices such that no force
element uses multiple vertices of the same color.

(a) Vertex colors: 3

(b) Dual graph & element colors: 7

(c) Vertex colors: 8

(d) Element colors: 76

Figure 8. Coloring vertices vs. elements: (a) vertex coloring needs 3 colors
for 10 vertices while (b) element coloring (i.e. coloring the vertices of the
dual graph) needs 7 colors for 9 triangles in this example. The difference is
more pronounced for tet-meshes: (c) our vertex coloring uses only 8 colors for
3,891 vertices while (d) our element coloring implementation needs 76 colors
for 14,802 tetrahedra in this example.

The advantage of VBD here is that, because it colors vertices, it typically
results in much fewer colors as compared to techniques that color
constraints/force elements, such as PBD. This is because coloring these
elements is equivalent to coloring the vertices of the dual graph, which not
only has more vertices, but, more importantly, also has a lot more connections
per vertex in general. Examples of this are shown in Figure 8. Since different
colors must be handled sequentially, fewer colors means better parallelism.

When all force elements are known ahead of the simulation, graph coloring can
be performed as a preprocess. However, dynamically generated constraints/force
elements, such ones due to collisions, cannot be known ahead of time,
requiring dynamic recoloring.

In our implementation for elastic body dynamics, we avoid the cost of
recoloring by precomputing coloring only for material forces, ignoring
dynamically generated force elements due to collisions. This means that these
collision forces may use multiple vertices of the same color. Therefore, we
cannot simply run a parallel loop over all vertices of the same color and
update them, because handling a vertex with a dynamically generated force
element may run into race conditions (with partially updated vertex positions)
when accessing other vertices of the same color.

We resolve this by having an auxiliary vertex position buffer that stores the
updated vertex position. When executing each local VBD position update, we
write the updated vertex positions to the auxiliary buffer instead of directly
overwriting the original vertex position buffer. After the position updates
for each color, we copy the updated positions from the auxiliary buffer to the
original vertex position buffer. This method prevents the race condition that
arises from simultaneous read and write operations on a vertex position.

Under this mechanism, dynamically generated force elements using multiple
vertices of the same color result in (partially) Jacobi-style iterations for
those vertices. This is because the process relies on the positions from the
previous iteration of vertices sharing the same color. For different-colored
vertices, it is equivalent to a GS approach, considering the updated positions
of other colors. Note that this differentiation is implicitly handled through
our racing-free scheme, without explicitly separating the approaches in our
code.

(a) No recoloring (ours)

(b) Recoloring

Figure 9. Handling collisions using (a) our scheme without recoloring and (b)
recoloring to achieve perfect Gauss-Seidel iterations, both simulated using
friction forces and accelerated iterations with $\rho=0.95$. Notice that the
results are highly similar, though not identical.

One might expect this solution to negatively impact the convergence of VBD.
Fortunately, however, such Jacobi-style information exchanges are relatively
rare. This is because, as shown in Figure 4, with vertex-face collisions at
most one pair and with edge-edge collisions at most two pairs of vertices can
have the same color. Also, vertices with the same colors must be located on
different sides of the collision; therefore, their elastic energies are
usually decoupled. This ensures that the majority of the information exchange
follows the Gauss-Seidel order and thereby the impact of our solution on
convergence is minimal. Figure 9 presents an example with a large number of
collisions, comparing our solution of skipping recoloring to proper Gauss-
Seidel iterations with recoloring, showing that the differences are minor.

Our solution also works with other types of topological changes, such as
tearing and fracturing. Deleting force elements does not require any changes
to vertex coloring. When an object is split by duplicating vertices, as in the
case of tearing a piece of cloth along some edges (see Figure 10), duplicated
vertices can safely inherit the colors of their original vertices.

###  3.10. GPU Implementation

In this section, we describe a GPU implementation specifically designed to
leverage the inherent parallelism of modern GPUs, which consists of two
hierarchical levels: block-level and thread-level parallelism. Block-level
parallelism facilitates large-scale parallel operations, with the stipulation
that each block operates independently. On the other hand, thread-level
parallelism provides finer, single-instruction-multiple-thread (SIMT) style
parallelism, allowing for inter-thread communication and synchronization
within the same block.

Reflecting on Vertex Buffer Data (VBD), we observed that it naturally aligns
with this hierarchical architecture. We have thousands of vertices within a
single color category that operate independently, and each vertex is
associated with multiple force elements, which can be processed in parallel.
Hence we provide the algorithm described by Algorithm 1. The cornerstone of
this strategy involves assigning a thread block to each vertex (line 7), as
opposed to using a single thread. Within these blocks, each thread
concurrently computes the force and Hessian for an adjacent force element,
storing the outcomes in the shared memory (lines 8$\sim$11). We allocate a
fixed number of threads for each thread block, denoted by $n_{t}$. When the
total number of adjacent force elements surpasses $n_{t}$, individual threads
are tasked with looping over their assigned elements. During this process,
they calculate the forces and Hessians for these elements and then
collectively sum these values in their assigned shared memory. At last, the
forces and Hessians of all the threads are combined using a local reduction
sum method (lines 12,13).

This parallelization mechanism turns out to be remarkably efficient. In our
tests, it runs nearly an order of magnitude faster than the trivial
parallelization, where a thread is assigned to each vertex to sequentially
process its adjacent force elements. The primary advantage lies in the
optimization of the memory access pattern, a common bottleneck in GPU
programs. This method notably minimizes memory divergence within warps. Given
that elements adjacent to a vertex create an n-dimensional disc (2 for
triangular mesh and 3 for tetrahedral mesh), these neighboring elements often
share multiple vertices. For instance, in a tetrahedral mesh, a vertex
$\vec{x}_{a}$ adjacent to another vertex $\vec{x}_{b}$ is shared by at least
three, but more commonly over five, neighboring tetrahedra of $\vec{x}_{b}$.
This shared access can propagate to all threads within the same warp, making a
single global memory access far more efficient. Approximately, in our
parallelism scheme, one global memory access can be shared by about 5 threads.

Furthermore, this strategy facilitates parallel processing in evaluating the
force and Hessian of adjacent elements, which are then written to the
significantly faster shared memory. This bypasses the slower global memory and
enables the parallel aggregation of force and Hessian values, enhancing the
overall efficiency of the process.

Figure 10. Tearing a piece of cloth with 2500 vertices and 4800 triangles.

Input: $\vec{x}^{t}$: the positions of the previous step; $\vec{v}^{t}$: the
velocities of the previous step; $\vec{a}_{\text{ext}}$: the external
acceleration

Output: This step’s position $\vec{x}^{t+1}$ and velocity $\vec{v}^{t+1}$.

1 $\vec{y}\leftarrow\vec{x}^{t}+h\vec{v}^{t}+h^{2}\vec{a}_{\text{ext}}$

2 Initial DCD using $\vec{x}^{t}$

3 $\vec{x}\leftarrow$ initial guess with adaptive initialization

4 for each  _iteration $n\leq n_{\max}$_ do

5 if  _$n\\!\\!\mod{n_{\text{col}}}=1$_ then CCD using $\vec{x}$

6

7 for each  _color $c$_ do

// block level parallelization

8 parallel for each  _vertex $i$ in color $c$_ do

// thread level parallelization

9 parallel for each  _$j\in\mathcal{F}_{i}$_ do

// variables in shared memory

10 $\mathbf{f}_{i,j}=-\frac{\partial E_{j}}{\partial\mathbf{x}_{i}}$

11
$\mathbf{H}_{i,j}=\frac{\partial^{2}E_{j}}{\partial\mathbf{x}_{i}\partial\mathbf{x}_{i}}$

12

13 end

// local reduction sums

14 $\mathbf{f}_{i}=\sum_{j\in\mathcal{F}_{i}}\mathbf{f}_{i,j}$

15 $\mathbf{H}_{i}=\sum_{j\in\mathcal{F}_{i}}\mathbf{H}_{i,j}$

16

17 $\Delta\vec{x}_{i}\leftarrow\mathbf{H}_{i}^{-1}\mathbf{f}_{i}$

18 $\Delta\vec{x}_{i}\leftarrow$ optional line search from
$\vec{x}_{i}+\Delta\vec{x}_{i}$ to $\mathbf{x}_{i}$

19 $\vec{x}_{i}^{\text{new}}\leftarrow\vec{x}_{i}+\Delta\vec{x}_{i}$

20

21 end

22

// Copy updated positions back to the vertex buffer to achieve GS iteration

23 parallel for each  _vertex $i$ in color $c$_ do

24 $\mathbf{x}_{i}=\vec{x}_{i}^{\text{new}}$

25 end

26

27 end

// Optional: accelerated iteration

28 parallel for each  _vertex $i$_ do

29 Update $\mathbf{x}_{i}$ using Equation 18.

30 end

31

32 end

33 $\vec{v}=(\vec{x}-\vec{x}^{t})/h$

34 return $\vec{x}$, $\vec{v}$

Algorithm 1 VBD simulation for one time step on GPU.

Figure 11. Simulation of 216 squishy balls with tentacles, a total of 48
million vertices and 151 million tetrahedra, dropped into a Utah teapot,
forming a stable pile with active frictional contacts. The average and maximum
computation time per time step is 3.6 and 3.9 seconds, respectively, using
$S=4$ substeps per frame and $n_{\max}=40$ iterations per step. The final
frame of this simulation is shown in Figure 1.

$\uparrow$ This platform
will be removed.

Figure 12. Simulation of 10,368 deformable objects, totaling over 36 million
vertices and 124 million tetrahedra, dropped onto a platform inside a box
container. Then, the platform is suddenly removed and the objects collectively
fall onto the ground, forming stable piles both before and after the platform
is removed. The average and maximum computation times per time step are 4.2
and 4.7 seconds, respectively, using $S=2$ substeps and $n_{\max}=60$
iterations per step. The final frame of this simulation is shown in Figure 1.

##  4\. Results

We evaluate our method with elastic body dynamics qualitatively with various
tests and quantitatively with direct comparisons to alternative methods. We
use Neo-Hookean(Smith et al., 2018) materials (without the logarithmic term)
for our volumetric objects, StVK (Volino et al., 2009) for clothes, and linear
springs for elastic rods.

We use a fixed frame time of $1/60$ seconds and a fixed iteration count
$n_{\max}$, instead of estimating convergence after every iteration. Each
frame is computed using $S$ substeps, such that ${h=1/(60S)}$ seconds. Using
smaller time steps increases accuracy and reduces numerical damping with any
implicit Euler method. With VBD, a smaller time step only requires fewer
iterations per step for similar visual quality, but it can also achieve a
smaller residual with the same number of total iterations per frame (i.e.,
$Sn_{\max}$).

We use no line search in our tests, as none of our tests required it for
stability, and running line search can result in a noticeable drop in
performance. In our experiments, we apply CCD only in the first iteration
(i.e. $n_{\text{col}}=n_{\max}$), unless otherwise specified.

In our implementation, we handle collisions on the CPU using Intel’s Embree
library (Wald et al., 2014). The two phases with parallel loops are
implemented on the GPU using CUDA. All timing results are generated on a
computer with an AMD Ryzen 5950X CPU, 64GB DDR3 RAM, and an NVIDIA RTX 4090
GPU.

###  4.1. Large-Scale Tests

In Figure 1 we present two large-scale test, showcasing our method’s
performance, scalability, and stability in scenarios involving a large number
of complex collisions, including stacking and rest in contact.

The first one has 216 squishy balls with tentacles, totaling 48 million
vertices and 151 million tetrahedra acting as force elements, dropped into a
Utah teapot. Intermediate frames of this simulation are shown in Figure 11.

The second one includes more than 10 thousand deformable objects, totaling
over 36 million vertices and 124 million tetrahedra, dropped into a box and
piled on a platform, which is then suddenly removed, making the pile
collectively fall onto the ground. The intermediate frames are shown in Figure
12.

As can be seen in our supplemental video, both of these simulations exhibit
stable motion, quickly forming static piles, while maintaining rest-in-contact
behavior with over 1 million active collisions. VBD’s parallelism and fast
convergence resulted in an average computation time of 16.4 and 7.2 seconds
per frame in these simulations, respectively.

Figure 13. Visual convergence with different numbers of iterations per frame
for different material stiffness (with accelerated iterations using
${\rho=0.93,0.86,0.75}$ top to bottom), simulating a beam with 463 vertices
and 1.5 thousand tetrahedra.

###  4.2. Unit Tests

The convergence rate of VBD depends on the stiffness of the simulated system.
This is demonstrated in Figure 13, comparing VBD with different numbers of
iterations per frame to the converged solution computed using Newton’s method.
As expected, VBD converges slower for stiffer materials, which is common for
descent-based solvers. In this example, 15 iterations are more than sufficient
for the softest material, while stiffer ones require more. As can be seen in
our supplemental video, even though VBD can qualitatively imitate the behavior
of stiff materials with a relatively small number of iterations, the motion
can quickly diverge from the converged solution, due to the remaining
residual, unless a sufficient number of iterations are used.

$\mu_{c}=0.0$

$\mu_{c}=0.3$

$\mu_{c}=0.6$

$\mu_{c}=0.9$

$\mu_{c}=0.0$

$\mu_{c}=0.2$

$\mu_{c}=0.4$

$\mu_{c}=0.6$

Figure 14. Testing different friction coefficients $\mu_{c}$ for (top) an
elastic cube with 400 vertices and 1.45 thousand tetrahedra, initially resting
on an incline, and (bottom) 4 elastic octopus models, totaling 15.6 thousand
vertices and 60 thousand tetrahedra, dropped into a box.

We present our tests with different friction coefficients $\mu_{c}$ for
friction forces in Figure 14. Notice that, $\mu_{c}$ impacts the motion, as
expected, and with sufficiently high $\mu_{c}$, we can properly preserve the
position on an incline and form taller piles.

In our supplemental video, we also include a comparison of different damping
stiffness $k_{d}$, showing that, despite numerical damping of implicit Euler,
without damping we can preserve kinetic energy for a long time. As we increase
$k_{d}$, the motion subsides quicker, as expected.

###  4.3. Stress Tests

We present a challenging frictional contact case in Figure 2, twisting two
thin beams together. This example includes extreme deformations, generating
strong material forces that compete with self-collisions and collisions
between the two beams. It is simulated with collision detection occurring once
every 5 iterations (i.e. $n_{\text{col}}=5$). Notice that VBD can properly
handle such strong deformations with frictional contact.

We show two simulations in Figure 3 for stress-testing the stability of VBD
under extreme deformations. The first one shows an armadillo model that is
perfectly flattened and the second one is a Utah teapot model with all
vertices randomly scrambled and placed on the surface of a sphere. Even though
the simulations begin with these extremely unstable energy configurations, VBD
quickly recovers these models without performing a line search and using 100
iterations per frame.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15. A stress test with extreme stretching: a Stanford bunny model with
1.8 thousand vertices and 5.9 thousand tetrahedra is stretched by slowly
pulling 10 vertices away, which are then suddenly released. The model recovers
its shape after going through considerable deformations and high-velocity
motion, simulated with self-collisions and using $S=5$ substeps and
$n_{\max}=10$ iterations per step.

Another stress test is shown in Figure 15. In this case, 10 vertices of a
Stanford bunny model are first slowly pulled away, generating a state with
considerable potential energy, and then suddenly released (right after Figure
15b), causing severe deformations. Once again, VBD successfully handles this
challenging simulation case, involving self-collisions with high-velocity
impacts, using only $n_{\max}=10$ iterations per step and $S=5$ per frame.

Figure 16. A stress test using only a single iteration per frame (i.e. a time
step of $h=1/60$ seconds and $n_{\max}=1$). One vertex on the armadillo
model’s nose is pulled while the finger and toe vertices are fixed. The model
has 15 thousand vertices and 50 thousand tetrahedra.

Figure 16 presents a stability test under large residuals by using only a
single iteration per frame (i.e. $S=1$ and $n_{\max}=1$). Notice that our
method produces stable deformations with extreme stretching, even though the
simulation lacks a sufficient iteration count to properly reduce the residual
at each frame.

[0.005] [0pt][0pt] [0.005] [0pt][0pt]

Figure 17. Convergence of different descent methods for simulating an
armadillo model with 15 thousand vertices and 50 thousand tetrahedra with
(top) a relatively soft material and (bottom) a 10$\times$ stiffer material.
Vertices near the top inside the glass block are fixed and the models are
initially stretched, as shown on the left, by pulling down foot vertices.
Then, the position constraints on foot vertices are suddenly removed, allowing
the model to deform for 33 ms. The deformation is computed using a single time
step of $h=33$ ms. The graphs show relative loss over iterations and
computation time. All methods are implemented on the GPU using the same
framework with single precision (32-bit) floating-point numbers, except for
Newton’s method with Cholesky factorization, which runs on the CPU using
double precision (64-bit). Accelerated versions use $\rho=0.95$.

###  4.4. Convergence Rate

To evaluate the convergence rate of VBD, we present a simple test shown in
Figure 17, where an armadillo model that was previously stretched is suddenly
released. Here, we calculate the relative loss after each iteration and
compare it to alternative solvers.

The first alternative is preconditioned gradient descent (GD) (Wang and Yang,
2016) implemented on GPU within the same framework as ours. GD requires a form
of line search for stability, which is implemented as testing the variational
energy after every 8 iterations and, when needed, reducing the optimization
step size and goes back to the state 8 iterations ago (following the
implementation of Wang and Yang (2016)). We also include a version of GD that
is accelerated using the Chebyshev semi-iterative approach (Wang, 2015), as
our method. Although GD iterations are approximately 30% faster than our
method, it necessitates a line search to prevent explosion at the given state.
With line search, it is 10% slower than our method. Its convergence rate per
iteration is considerably slower, as it corresponds to Jacobi iterations. In
this example, when combined with acceleration, GD performs similarly to our
method without acceleration but lags significantly behind our method with
acceleration.

We also compare to a version of our method that uses Jacobi iterations, called
_Block Jacobi_ , implemented by computing the position change for all vertices
in parallel and then applying the position update simultaneously to all
vertices at the end of each iteration. We incorporate the same line search
scheme as GD for Block Jacobi, as it is necessary to prevent explosion. Block
Jacobi outperforms GD, as it uses vertex blocks for computation, which
corresponds to using a diagonal Hessian block as a precondition, as opposed to
just a diagonal line that GD uses. Without the line search, it runs 20% faster
than the GS version, due to its perfect vertex-level parallelization (without
any coloring). Combined with line search, it is 17% slower than the GS
version. Because it uses Jacobi iterations, its convergence is hindered, as
compared to our Gauss-Seidel iterations.

Furthermore, we compare our method’s convergence to Newton’s method using a
direct Cholesky (LDLT) factorization solver provided by Intel’s MKL library
(Wang et al., 2014). To ensure convergence, we do a PSD projection for each
tet’s Hessian of elasticity. Newton’s method uses a line search for each
iteration. Since its iterations are slow, the computational overhead of line
search is negligible. Though the convergence of Newton’s method per iteration
is far superior to all others, because of its expensive computation time per
iteration, in these examples it lags far behind. Nonetheless, for relatively
tame experiments with less stretching and motion, and especially for highly
stiff and high-resolution simulations that are much more expensive to
simulate, we would expect Newton’s method to eventually overtake all
alternatives beyond a certain level of convergence.

Table 1. Performance results and simulation parameters.

Experiment Name | Number of | Mateiral | Contact &Fiction | Simulation Parameters | Time per step
---|---|---|---|---|---
| Vert. | Tet. | color | Type | Stiffness | Damping | $k_{c}$ | $\mu_{c},\epsilon_{v}$ |  $h$(sec.) | Iterations | avg./max
Twisting Thin Beams (Figure 2) | 97K | 266K | 8 | NeoHookean | $\mu=5e4,\lambda=1e6$ | 1e-6 | 1e6 | 0.1, 1e-2 | 1/300 | 100 | 60/78ms
Flattening initialization (Figure 3) | 15K | 50K | 8 | NeoHookean | $\mu=2e6,\lambda=1e7$ | 1e-6 | NA | NA | 1/60 | 100 | 3.3/3.8ms
Random initialization (Figure 3) | 2K | 8.5K | 8 | NeoHookean | $\mu=2e6,\lambda=1e7$ | 1e-6 | NA | NA | 1/60 | 100 | 2.6/2.8ms
Tetmesh Pendulum (Figure 6) | 304 | 755 | 6 | NeoHookean | $\mu=1e7,\lambda=1e8$ | 0 | NA | NA | 1/300 | 20 | ¡0.1ms
Squishy Ball Drops (Figure.7,9,18) | 230K | 700K | 8 | NeoHookean | $\mu=2e6,\lambda=2e7$ | 1e-7 | 1e7 | 0.1, 1e-2 | 1/120 | 120 | 15/17ms
Tearing Cloth (Figure 10) | 2500 | 4800 | 3 | StVK | $\mu=1e4,\lambda=1e4$ | 1e-5 | NA | NA | 1/300 | 20 | 11.2/11.5 ms(cpu)
Dropping 216 Squshy Balls (Figure 11) | 48M | 151M | 9 | NeoHookean | $\mu=2e6,\lambda=2e7$ | 1e-7 | 1e7 | 0.1, 1e-2 | 1/240 | 40 | 3.6/3.9s
Dropping 10368 Models (Figure 12) | 36M | 124M | 8 | NeoHookean | $\mu=1e6,\lambda=1e7$ | 1e-7 | 1e7 | 0.1, 1e-2 | 1/120 | 60 | 4.2/4.7s
Beam Sagging (Figure 13) | 463 | 1.5K | 6 | NeoHookean | $\mu=1e6/3e6/1e7$ | 1e-6 | NA | NA | 1/300 | 3/5/10 | avg.: 0.08/0.12/0.24ms
|  |  |  |  | $\lambda=1e7/3e7/1e8$ |  |  | NA |  |  | max: 0.08/0.16/0.31ms
Cude Sliding (Figure 14) | 800 | 2.9K | 6 | NeoHookean | $\mu=1e6,\lambda=1e7$ | 1e-6 | 1e7 | 0/0.3/0.6/0.9, 1e-2 | 1/300 | 10 | 0.10/0.17ms
Octopi Stacking (Figure 14) | 15.6K | 60K | 8 | NeoHookean | $\mu=1e6,\lambda=1e7$ | 1e-6 | 1e7 | 0/0.1/0.4/0.6, 1e-2 | 1/300 | 10 | 1.1/1.3ms
Extreme Stretch (Figure 16) | 1.8K | 5.9K | 8 | NeoHookean | $\mu=2e6,\lambda=1e7$ | 1e-6 | 1e7 | 0.2, 1e-2 | 1/300 | 10 | 0.36/1.02ms
2 Cube colliding (Figure 19) | 800 | 2.9K | 6 | NeoHookean | $\mu=1e6,\lambda=1e7$ | 1e-6 | 1e7 | 0.3, 1e-2 | 1/300 | 10 | 0.21/0.27ms

###  4.5. Comparisons to XPBD

Our method has an entirely different formulation than XPBD, but there are some
strong similarities, as both methods operate with position updates using
Gauss-Seidel iterations. Here we provide two direct comparisons to highlight
some important differences.

(a) XPBD
$h=1/120$ sec.
$n_{\max}=120$
$n_{\text{col}}=120$
_( $0.32$ sec./frame)_

(b) XPBD
$h=1/3000$ sec.
$n_{\max}=5$
$n_{\text{col}}=125$
_( $0.35$ sec./frame)_

(c) XPBD
$h=1/3000$ sec.
$n_{\max}=5$
$n_{\text{col}}=5$
_( $3.1$ sec./frame)_

(d) VBD (ours)
$h=1/120$ sec.
$n_{\max}=120$
$n_{\text{col}}=120$
_( $0.03$ sec./frame)_

Figure 18. A squishy ball with tentacles, comprising 230 thousand vertices and
700 thousand tetrahedra, dropped on the ground, simulated using (a) XPBD with
a large time step and 240 iterations per frame, (b) XPBD with a $25\times$
smaller time step and 250 total iterations per frame, (c) XPBD with the same
small time step and iteration count but with $25\times$ more frequent
collision detection, and (d) VBD with a large time step and 240 iterations per
frame. Comparing (a) and (d), VBD is faster than XPBD with the same settings.
XPBD’s solution approaches VBD as the time step decreases, but it also
requires more frequent collision detection to achieve a visually similar
result to VBD.

XPBD replaces the Hessian matrix and uses only the Hessian of inertia
potential. This omission is justified by using a small time step, because the
significance of the inertia potential increases quadratically as the time step
decreases. Nonetheless, with complex examples, the impact of this
approximation can be severe, even with small time step. This is demonstrated
in Figure 18 with a challenging collision-rich scenario involving a squishy
ball with tentacles dropped to the ground. Comparing XPBD with 120 iterations
per step (18(a)) to our method with the same number of iterations (18(d)), we
can see that our method not only achieves a more stable animation, it also
performs faster because of its improved parallelism, as compared to XPBD.
Reducing the time step helps XPBD even when using a similar total number of
iterations per frame (18(b)). However, simply reducing the time step is not
sufficient in this case, as XPBD also needs more frequent collision detection
(18(c)). Using collision detection with the same frequency as ours while
taking small time step (18(b)) leads to collisions that are detected too late
and cause stability issues in this case. This is not only because the
collisions that are detected too late are deeper, but also because smaller
time step lead to higher vertex velocities when resolving stiff collisions.

(a) Initial State

(b) XPBD

(c) VBD (ours)

Figure 19. Dropping a large and heavy elastic cube onto a smaller and much
lighter box with a mass ratio of 1:2000. Each cube has 400 vertices and 1.5
thousand tetrahedra.

One of the fundamental challenges of XPBD is handling high mass ratios. This
is demonstrated with a simple example in Figure 19, where a large and heavy
elastic cube is dropped onto a smaller and much lighter cube, with a mass
ratio of 1:2000. In this example, XPBD’s collision constraints, even with
infinite stiffness, cannot overcome the mass ratio and the smaller cube is
entirely crushed upon contact. This is because of the dual formulation of XPBD
(Macklin et al., 2020). Our method, on the other hand, has no such
difficulties with handling high mass ratios.

##  5\. VBD for Other Simulation Systems

We have described our method in Section 3 in the context of elastic body
dynamics. Yet, VBD is not limited to such simulations and can be used to solve
various optimization problems. Here, we consider some other example simulation
systems and briefly discuss how our method can be applied. This is not
intended as an exhaustive list but merely as examples that could guide the
reader to discern how their specific simulation problem could utilize VBD.

###  5.1. Particle-Based Simulations

Particle-based simulations can easily use VBD by simply replacing the vertices
in our description above with particles. Since VBD needs the Hessian of the
force element energies, implementations would require computing the
derivatives of all forces acting on a particle wrt. its position.

Parallelizing particle-based simulations also involves additional
considerations. Mass-spring type simulations, such as peridynamics (Levine et
al., 2015), can use our parallelization approach with vertex coloring.
However, simulations involving disjoint or loosely-joined particles, such as
particle-based fluid simulation (Müller et al., 2003; Takahashi et al., 2015;
Peer et al., 2015), would not only require recoloring at each time step but
also using a conservative neighborhood definition (including potential
position change within a time step) for coloring, since position updates can
alter the set of neighboring particles that interact with each particle.

stiff soft
springs springs

Figure 20. 20 particles attached with springs, forming a swinging chain,
simulated using VBD with a $S=1$ substep and $100$ iterations per step. The
particle on one end of the chain is fixed and the particle on the other end
has 1000$\times$ more mass than the others. (Left) using sufficiently stiff
springs, they expand no more than 0.7% of their rest lengths, despite the
substantial mass difference. (Right) using 100$\times$ less stiff springs, the
chain undergoes a visible expansion as it swings.

Figure 20 shows a simple example where 20 particles, including one that is
1000$\times$ heavier, are connected with springs of two different stiffness,
simulated using VBD.

###  5.2. Rigid Body Simulation

For handling rigid body simulations with VBD, we can replace each vertex in
our formulation with an entire rigid body, using the variational formulation
of rigid body dynamics (Ferguson et al., 2021). Unlike a vertex that has only
3 DoF, a rigid body also has rotational DoF, resulting in 6 DoF. Therefore, in
our local system, we must solve a larger problem, where
$\mathbf{x}_{i}\in\mathbb{R}^{6}$ and $\vec{H}_{i}\in\mathbb{R}^{6\times 6}$,
including Hessians of all force elements wrt. all 6 DoF of $\mathbf{x}_{i}$.
Note that, in this case, these force elements are not internal material
forces, but external forces acting on the rigid body, due to collisions or
other constraints.

Other than this additional complexity, we can follow the same procedure with
VBD. Parallelization with coloring depends on the nature of the rigid body
simulation. For example, pre-coloring, as we used in our examples for elastic
bodies might work for problems like a rigid body chain. For disjoint rigid
bodies interacting through collisions only, dynamic recoloring might be
needed.

Articulated rigid bodies can be handled by defining joint constraints with an
elastic potential. Infinitely stiff constraints are also possible, but VBD
cannot guarantee that they will be satisfied using a fixed number of
iterations. Another alternative is hard constraints can be introduced by
reducing the total DoF in the system and replacing the vertices in our
formulation with an articulated rigid body, having more than 6 DoF. Obviously,
this would lead to an even larger local system, requiring modifications to the
variational formulation.

Figure 21. 5 rigid bodies, each with 6 DoF, forming a chain through
collisions, simulated using our VBD formulation for rigid body dynamics.

Figure 22. Dropping 60 rigid bodies into a Utah teapot, showcasing collisions
and frictional contact. Remarkably, one rigid body stays on the spout due to
friction.

Example rigid body simulations are shown in Figure 21 and Figure 22, simulated
using VBD, as described above.

Note that since our collision formulation is based on penetration potential,
it corresponds to penalty forces. We leave the exploration of handling
impulse-based collisions (Mirtich, 1996) with VBD to future work.

###  5.3. Unified Simulations

Unified simulation systems are useful for handling scenarios that involve
different material types. Typical unified simulation systems use a fundamental
building block, such as a particle, to represent all supported materials
(Müller et al., 2004; Solenthaler et al., 2007; Becker et al., 2009; Martin et
al., 2010; Macklin et al., 2014). We can form a unified simulation system
using VBD without representing all materials using the same building block.
For any simulation system described above, we can combine it with another,
provided that we can define the information exchange as an energy potential.
For example, when the interactions take place as collisions, we can easily
join rigid body simulations with elastic bodies or particles via the collision
potential. Joint constraints with elastic potential would be another easy way
to combine different simulation systems. The advantage here is that a large
rigid body, for example, can be represented as a single object with just 6
DoF, as opposed to using multiple building blocks that are constrained to move
as a rigid construction. This way, a heterogeneous collection of
representations can be joined within the same integrator using VBD.

On the other hand, this form of defining a unified simulation system may be
challenging for other types of information exchange, such as evaluating
buoyancy. Exploring such problems would be another interesting direction for
future research.

##  6\. Limitations

VBD is a descent-based method that operates through local iterations.
Therefore, it may not be a good solution for problems that would benefit from
a global treatment.

The speed of information travel with VBD depends on the connections of
vertices and the number of iterations used. A perturbation applied on a vertex
can impact other vertices of a connected chain through force elements at most
as far as the number of colors within a single iteration. Therefore, VBD is
not ideal for high-resolution stiff systems, as it may require too many
iterations for a perturbation of a vertex to travel across the system. In such
cases, a global solution using Newton’s method may prove to be more effective.

Our collision formulation for VBD is based on penetration potential.
Therefore, it cannot guarantee penetration-free results. In fact, penetrations
are almost never completely resolved, as some amount of penetration is needed
to maintain some collision force. Exploring penetration-free collisions with
VBD would be an interesting direction for future research.

In addition, defining a similar collision energy for codimensional objects,
particularly for self-collisions, can be a challenge.

(a) 100 iterations

(b) Converged

(c) 100 iterations

(d) Converged

Figure 23. A chain of particles connected with soft springs (orange) and
10,000$\times$ stiffer (blue) springs. Simulations with VBD using (a,c) 100
iterations per frame fail to converge and result in excessive extensions, as
compared to (b,d) convreged results.

VBD is a primal solver (Macklin et al., 2020), so it can easily handle high
mass ratios (see Figure 6, 19, and 20), but it struggles with high stiffness
ratios. This is shown in Figure 23 using a stiffness ratio of 1:10000, where
VBD has poor convergence behavior.

##  7\. Conclusion

We have presented vertex block descent, an efficient iterative descent-based
solution for optimization problems, and described how it can be used for
physics-based simulations with implicit Euler integration defined through a
variational formulation. We have explained all essential details of elastic
body dynamics using VBD, including handling damping, constraints, collisions,
and friction. We have defined an adaptive initialization technique, enabled by
VBD’s formulation, and discussed how to use momentum-based acceleration to
improve convergence. We have also presented effective methods for
parallelization using VBD, considering dynamically introduced/removed force
elements, and explained how its vertex-level computation improves the
parallelization of its Gauss-Seidel iterations.

Our results show that VBD can handle highly complex simulation cases (Figure
1), it remains stable under extreme stress tests (Figure 2, 3, 15, and 16),
and offers fast convergence (Figure 17).

In addition, we have summarized how VBD can be used for other types of
simulation problems, such as particle systems and rigid bodies, including
unified simulations. We have also mentioned some related future research
directions and discussed VBD’s limitations.

## References

  * (1)
  * Ascher and Boxerman (2003) Uri M Ascher and Eddy Boxerman. 2003.  On the modified conjugate gradient method in cloth simulation.  _The Visual Computer_ 19 (2003), 526–531.
  * Baraff and Witkin (1998) David Baraff and Andrew Witkin. 1998.  Large Steps in Cloth Simulation.  (1998).
  * Becker et al. (2009) Markus Becker, Markus Ihmsen, and Matthias Teschner. 2009.  Corotated SPH for deformable solids. In _Proceedings of the Fifth Eurographics Conference on Natural Phenomena_ (Munich, Germany) _(NPH’09)_. Eurographics Association, Goslar, DEU, 27–34.
  * Bolz et al. (2003) Jeff Bolz, Ian Farmer, Eitan Grinspun, and Peter Schröder. 2003.  Sparse matrix solvers on the GPU: conjugate gradients and multigrid.  _ACM transactions on graphics (TOG)_ 22, 3 (2003), 917–924.
  * Bouaziz et al. (2023) Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2023.  Projective dynamics: Fusing constraint projections for fast simulation.  In _Seminal Graphics Papers: Pushing the Boundaries, Volume 2_. 787–797.
  * Bridson et al. (2002) Robert Bridson, Ronald Fedkiw, and John Anderson. 2002.  Robust treatment of collisions, contact and friction for cloth animation. In _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_. 594–603.
  * Bridson et al. (2005) Robert Bridson, Sebastian Marino, and Ronald Fedkiw. 2005.  Simulation of clothing with folds and wrinkles.  In _ACM SIGGRAPH 2005 Courses_. 3–es.
  * Capell et al. (2002) Steve Capell, Seth Green, Brian Curless, Tom Duchamp, and Zoran Popović. 2002.  A multiresolution framework for dynamic deformations. In _Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation_. ACM, San Antonio Texas, 41–47.  https://doi.org/10.1145/545261.545268
  * Chao et al. (2010) Isaac Chao, Ulrich Pinkall, Patrick Sanan, and Peter Schröder. 2010.  A simple geometric model for elastic deformations.  _ACM Transactions on Graphics_ 29, 4 (July 2010), 1–6.  https://doi.org/10.1145/1778765.1778775
  * Chen et al. (2023) He Chen, Elie Diaz, and Cem Yuksel. 2023.  Shortest Path to Boundary for Self-Intersecting Meshes.  42, 4, Article 146 (2023), 15 pages.  https://doi.org/10.1145/3592136
  * Choi and Ko (2005) Kwang-Jin Choi and Hyeong-Seok Ko. 2005.  Stable but responsive cloth. In _ACM SIGGRAPH 2005 Courses on - SIGGRAPH ’05_. ACM Press, Los Angeles, California, 1.  https://doi.org/10.1145/1198555.1198571
  * Eberhardt et al. (2000) B. Eberhardt, O. Etzmuß, and M. Hauth. 2000.  Implicit-Explicit Schemes for Fast Animation with Particle Systems.  In _Computer Animation and Simulation 2000_ , W. Hansmann, W. Purgathofer, F. Sillion, Nadia Magnenat-Thalmann, Daniel Thalmann, and Bruno Arnaldi (Eds.). Springer Vienna, Vienna, 137–151.  https://doi.org/10.1007/978-3-7091-6344-3_11 Series Title: Eurographics.
  * Ferguson et al. (2021) Zachary Ferguson, Minchen Li, Teseo Schneider, Francisca Gil-Ureta, Timothy Langlois, Chenfanfu Jiang, Denis Zorin, Danny M. Kaufman, and Daniele Panozzo. 2021.  Intersection-free Rigid Body Dynamics.  _ACM Transactions on Graphics (SIGGRAPH)_ 40, 4, Article 183 (2021).
  * Fratarcangeli and Pellacini (2015) M. Fratarcangeli and F. Pellacini. 2015.  Scalable Partitioning for Parallel Position Based Dynamics.  _Computer Graphics Forum_ 34, 2 (May 2015), 405–413.  https://doi.org/10.1111/cgf.12570
  * Fratarcangeli et al. (2016) Marco Fratarcangeli, Valentina Tibaldo, and Fabio Pellacini. 2016.  Vivace: a practical gauss-seidel method for stable soft body dynamics.  _ACM Transactions on Graphics_ 35, 6 (Nov. 2016), 1–9.  https://doi.org/10.1145/2980179.2982437
  * Gast et al. (2015) Theodore F. Gast, Craig Schroeder, Alexey Stomakhin, Chenfanfu Jiang, and Joseph M. Teran. 2015.  Optimization Integrator for Large Time Steps.  _IEEE Transactions on Visualization and Computer Graphics_ 21, 10 (Oct. 2015), 1103–1115.  https://doi.org/10.1109/TVCG.2015.2459687
  * Golub and Van Loan (2013) Gene H Golub and Charles F Van Loan. 2013.  _Matrix computations_.  JHU press.
  * Grinspun et al. (2002) Eitan Grinspun, Petr Krysl, and Peter Schröder. 2002.  CHARMS: A simple framework for adaptive simulation.  _ACM transactions on graphics (TOG)_ 21, 3 (2002), 281–290.
  * Hauth and Etzmuss (2001) Michael Hauth and Olaf Etzmuss. 2001.  A High Performance Solver for the Animation of Deformable Objects using Advanced Numerical Methods.  _Computer Graphics Forum_ 20, 3 (Sept. 2001), 319–328.  https://doi.org/10.1111/1467-8659.00524
  * Hecht et al. (2012) Florian Hecht, Yeon Jin Lee, Jonathan R. Shewchuk, and James F. O’Brien. 2012.  Updated sparse cholesky factors for corotational elastodynamics.  _ACM Transactions on Graphics_ 31, 5 (Aug. 2012), 1–13.  https://doi.org/10.1145/2231816.2231821
  * Hirota et al. (2001) G. Hirota, S. Fisher, A. State, C. Lee, and H. Fuchs. 2001.  An implicit finite element method for elastic solids in contact. In _Proceedings Computer Animation 2001. Fourteenth Conference on Computer Animation (Cat. No.01TH8596)_. IEEE Comput. Soc, Seoul, South Korea, 136–254.  https://doi.org/10.1109/CA.2001.982387
  * Huthwaite (2014) Peter Huthwaite. 2014.  Accelerated finite element elastodynamic simulations using the GPU.  _J. Comput. Phys._ 257 (2014), 687–707.
  * Kane et al. (1999) C. Kane, J. E. Marsden, and M. Ortiz. 1999.  Symplectic-energy-momentum preserving variational integrators.  _J. Math. Phys._ 40, 7 (July 1999), 3353–3371.  https://doi.org/10.1063/1.532892
  * Kane et al. (2000) Couro Kane, Jerrold E Marsden, Michael Ortiz, and Matthew West. 2000.  Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems.  _International Journal for numerical methods in engineering_ 49, 10 (2000), 1295–1325.
  * Kharevych et al. (2006) Liliya Kharevych, Weiwei Yang, Yiying Tong, Eva Kanso, Jerrold E. Marsden, Peter Schröder, and Matthieu Desbrun. 2006.  Geometric, Variational Integrators for Computer Animation. In _ACM SIGGRAPH / Eurographics Symposium on Computer Animation_ , Marie-Paule Cani and James O’Brien (Eds.). The Eurographics Association.  https://doi.org/10.2312/SCA/SCA06/043-051
  * Lan et al. (2023) Lei Lan, Minchen Li, Chenfanfu Jiang, Huamin Wang, and Yin Yang. 2023.  Second-order Stencil Descent for Interior-point Hyperelasticity.  _ACM Transactions on Graphics_ 42, 4 (Aug. 2023), 1–16.  https://doi.org/10.1145/3592104
  * Lan et al. (2022) Lei Lan, Guanqun Ma, Yin Yang, Changxi Zheng, Minchen Li, and Chenfanfu Jiang. 2022.  Penetration-free projective dynamics on the GPU.  _ACM Transactions on Graphics (TOG)_ 41, 4 (2022), 1–16.
  * Levine et al. (2015) J. A. Levine, A. W. Bargteil, C. Corsi, J. Tessendorf, and R. Geist. 2015.  A peridynamic perspective on spring-mass fracture. In _Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Copenhagen, Denmark) _(SCA ’14)_. Eurographics Association, Goslar, DEU, 47–55.
  * Lew et al. (2004) A. Lew, J. E. Marsden, M. Ortiz, and M. West. 2004.  Variational time integrators.  _Internat. J. Numer. Methods Engrg._ 60, 1 (May 2004), 153–212.  https://doi.org/10.1002/nme.958
  * 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 Transactions on Graphics (TOG)_ 39, 6 (2020), 1–15.
  * Li et al. (2020a) Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy R 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 (2020), 49.
  * Li et al. (2019) Minchen Li, Ming Gao, Timothy Langlois, Chenfanfu Jiang, and Danny M. Kaufman. 2019.  Decomposed optimization time integrator for large-step elastodynamics.  _ACM Transactions on Graphics_ 38, 4 (Aug. 2019), 1–10.  https://doi.org/10.1145/3306346.3322951
  * Li et al. (2023) Xuan Li, Yu Fang, Lei Lan, Huamin Wang, Yin Yang, Minchen Li, and Chenfanfu Jiang. 2023.  Subspace-Preconditioned GPU Projective Dynamics with Contact for Cloth Simulation. In _SIGGRAPH Asia 2023 Conference Papers_. 1–12.
  * 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_ 36, 3 (June 2017), 1–16.  https://doi.org/10.1145/2990496
  * 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.  _Computer Graphics Forum_ 39, 8 (Dec. 2020), 89–100.  https://doi.org/10.1111/cgf.14104
  * Macklin et al. (2014) Miles Macklin, Matthias Müller, Nuttapong Chentanez, and Tae-Yong Kim. 2014.  Unified particle physics for real-time applications.  _ACM Trans. Graph._ 33, 4, Article 153 (jul 2014), 12 pages.  https://doi.org/10.1145/2601097.2601152
  * 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_. ACM, Burlingame California, 49–54.  https://doi.org/10.1145/2994258.2994272
  * Martin et al. (2010) Sebastian Martin, Peter Kaufmann, Mario Botsch, Eitan Grinspun, and Markus Gross. 2010.  Unified simulation of elastic rods, shells, and solids.  _ACM Trans. Graph._ 29, 4, Article 39 (jul 2010), 10 pages.  https://doi.org/10.1145/1778765.1778776
  * Martin et al. (2011) Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011.  Example-based elastic materials.  In _ACM SIGGRAPH 2011 papers_. 1–8.
  * Mirtich (1996) Brian Vincent Mirtich. 1996.  _Impulse-based dynamic simulation of rigid body systems_.  Ph. D. Dissertation.  AAI9723116.
  * Müller et al. (2003) Matthias Müller, David Charypar, and Markus Gross. 2003.  Particle-based fluid simulation for interactive applications. In _Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (San Diego, California) _(SCA ’03)_. Eurographics Association, Goslar, DEU, 154–159.
  * Müller et al. (2002) Matthias Müller, Julie Dorsey, Leonard McMillan, Robert Jagnow, and Barbara Cutler. 2002.  Stable real-time deformations. In _Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation_. 49–54.
  * Müller et al. (2007) Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007.  Position based dynamics.  _Journal of Visual Communication and Image Representation_ 18, 2 (2007), 109–118.
  * Müller et al. (2005) Matthias Müller, Bruno Heidelberger, Matthias Teschner, and Markus Gross. 2005.  Meshless deformations based on shape matching.  _ACM transactions on graphics (TOG)_ 24, 3 (2005), 471–478.
  * Müller et al. (2004) M. Müller, R. Keiser, A. Nealen, M. Pauly, M. Gross, and M. Alexa. 2004.  Point based animation of elastic, plastic and melting objects. In _Proceedings of the 2004 ACM SIGGRAPH/Eurographics Symposium on Computer Animation_ (Grenoble, France) _(SCA ’04)_. Eurographics Association, Goslar, DEU, 141–151.  https://doi.org/10.1145/1028523.1028542
  * Naitsat et al. (2020) Alexander Naitsat, Yufeng Zhu, and Yehoshua Y Zeevi. 2020.  Adaptive block coordinate descent for distortion optimization. In _Computer Graphics Forum_ , Vol. 39. Wiley Online Library, 360–376.
  * Peer et al. (2015) Andreas Peer, Markus Ihmsen, Jens Cornelis, and Matthias Teschner. 2015.  An implicit viscosity formulation for SPH fluids.  _ACM Trans. Graph._ 34, 4, Article 114 (jul 2015), 10 pages.  https://doi.org/10.1145/2766925
  * Sifakis and Barbic (2012) Eftychios Sifakis and Jernej Barbic. 2012.  FEM simulation of 3D deformable solids: A practitioner’s guide to theory, discretization and model reduction.  _ACM SIGGRAPH 2012 Courses, SIGGRAPH’12_ (08 2012).  https://doi.org/10.1145/2343483.2343501
  * Simo et al. (1992) J.C. Simo, N. Tarnow, and K.K. Wong. 1992.  Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics.  _Computer Methods in Applied Mechanics and Engineering_ 100, 1 (Oct. 1992), 63–116.  https://doi.org/10.1016/0045-7825(92)90115-Z
  * Smith et al. (2018) Breannan Smith, Fernando De Goes, and Theodore Kim. 2018.  Stable neo-hookean flesh simulation.  _ACM Transactions on Graphics (TOG)_ 37, 2 (2018), 1–15.
  * Solenthaler et al. (2007) Barbara Solenthaler, Jürg Schläfli, and Renato Pajarola. 2007.  A unified particle model for fluid–solid interactions: Research Articles.  _Comput. Animat. Virtual Worlds_ 18, 1 (feb 2007), 69–82.
  * Stern and Grinspun (2009) Ari Stern and Eitan Grinspun. 2009.  Implicit-Explicit Variational Integration of Highly Oscillatory Problems.  _Multiscale Modeling & Simulation_ 7, 4 (Jan. 2009), 1779–1794.  https://doi.org/10.1137/080732936 arXiv:0808.2239 [math].
  * Takahashi et al. (2015) Tetsuya Takahashi, Yoshinori Dobashi, Issei Fujishiro, Tomoyuki Nishita, and Ming C. Lin. 2015.  Implicit Formulation for SPH-based Viscous Fluids.  _Comput. Graph. Forum_ 34, 2 (may 2015), 493–502.  https://doi.org/10.1111/cgf.12578
  * Tamstorf et al. (2015) Rasmus Tamstorf, Toby Jones, and Steve McCormick. 2015.  Smoothed Aggregation Multigrid for Cloth Simulation.  _ACM Transactions on Graphics_ 34 (Oct. 2015), 1–13.  https://doi.org/10.1145/2816795.2818081
  * Teran et al. (2005) Joseph Teran, Eftychios Sifakis, Geoffrey Irving, and Ronald Fedkiw. 2005.  Robust quasistatic finite elements and flesh simulation. In _Proceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation_. 181–190.
  * Ton-That et al. (2023) Quoc-Minh Ton-That, Paul G. Kry, and Sheldon Andrews. 2023.  Parallel block Neo-Hookean XPBD using graph clustering.  _Computers & Graphics_ 110 (Feb. 2023), 1–10.  https://doi.org/10.1016/j.cag.2022.10.009
  * Volino and Magnenat-Thalmann (2001) P. Volino and N. Magnenat-Thalmann. 2001.  Comparing efficiency of integration methods for cloth simulation. In _Proceedings. Computer Graphics International 2001_. IEEE Comput. Soc, Hong Kong, China, 265–272.  https://doi.org/10.1109/CGI.2001.934683
  * Volino et al. (2009) Pascal Volino, Nadia Magnenat-Thalmann, and Francois Faure. 2009.  A simple approach to nonlinear tensile stiffness for accurate cloth simulation.  _ACM Trans. Graph._ 28, 4, Article 105 (sep 2009), 16 pages.  https://doi.org/10.1145/1559755.1559762
  * Wald et al. (2014) Ingo Wald, Sven Woop, Carsten Benthin, Gregory S Johnson, and Manfred Ernst. 2014.  Embree: a kernel framework for efficient CPU ray tracing.  _ACM Transactions on Graphics (TOG)_ 33, 4 (2014), 1–8.
  * Wang et al. (2014) Endong Wang, Qing Zhang, Bo Shen, Guangyong Zhang, Xiaowei Lu, Qing Wu, Yajuan Wang, Endong Wang, Qing Zhang, Bo Shen, et al. 2014\.  Intel math kernel library.  _High-Performance Computing on the Intel® Xeon Phi™: How to Fully Exploit MIC Architectures_ (2014), 167–188.
  * Wang (2015) Huamin Wang. 2015.  A chebyshev semi-iterative approach for accelerating projective and position-based dynamics.  _ACM Transactions on Graphics_ 34, 6 (Nov. 2015), 1–9.  https://doi.org/10.1145/2816795.2818063
  * Wang and Yang (2016) Huamin Wang and Yin Yang. 2016.  Descent methods for elastic body simulation on the GPU.  _ACM Transactions on Graphics_ 35, 6 (Nov. 2016), 1–10.  https://doi.org/10.1145/2980179.2980236
  * Wang et al. (2020) Xinlei Wang, Minchen Li, Yu Fang, Xinxin Zhang, Ming Gao, Min Tang, Danny M. Kaufman, and Chenfanfu Jiang. 2020.  Hierarchical Optimization Time Integration for CFL-Rate MPM Stepping.  _ACM Transactions on Graphics_ 39, 3 (June 2020), 1–16.  https://doi.org/10.1145/3386760
  * Wright (2015) Stephen J Wright. 2015.  Coordinate descent algorithms.  _Mathematical programming_ 151, 1 (2015), 3–34.
  * 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 Transactions on Graphics_ 38, 6 (Dec. 2019), 1–13.  https://doi.org/10.1145/3355089.3356486
  * Zhu et al. (2010) Yongning Zhu, Eftychios Sifakis, Joseph Teran, and Achi Brandt. 2010.  An efficient multigrid method for the simulation of high-resolution elastic solids.  _ACM Transactions on Graphics_ 29, 2 (March 2010), 1–18.  https://doi.org/10.1145/1731047.1731054