Affine Body Dynamics:
Fast, Stable and Intersection-free Simulation of Stiff Materials

Given a mass density distribution, $\rho$, over the body domain, $\Omega$, the kinetic energy of each affine body is then

with the generalized mass matrix for each affine body defined as:

Here we note an additional convenience of ABD: for flat affine coordinates we obtain a constant mass matrix. In turn this ensures that the equations of motion for affine bodies, unlike rigid bodies, do not add nonlinear Coriolis-type forces. Instead, they are embodied as generalized internal forces.

With $V$ the total potential energy the free ABD is then simply the equations of motion:

where external forces $f_{k}\in\mathbb{R}^{3}$, applied at material points $k$, are included as $f=\sum_{k}\mathsf{J}(\bar{x}_{k})^{T}f_{k}$.

In place of SE(3) coordinates we rigidify each affine body with a stiff orthogonality potential

scaled by the stiffness $\kappa$ and the body’s volume $\nu$. We apply a large stiffness ($\kappa>100$GPa) to ensure the deformation on the body is sufficiently suppressed and negligible.

As we relax the rigidity constraint, we instead apply a highly stiff penalty term. Generally engineering rule of thumb suggests that stiff penalties should be avoided in simulation, as they tend to exacerbate the numerical difficulties of computing with nonlinear potentials. However, we observe that contact and collision forces (irrespective of whether they are treated via constraints, penalties springs, or barriers) introduce a much more dominant stiffness to the system. In our case this already requires handling system solves with a robust Newton-type algorithm so that the overhead of a single, additional stiff orthogonality potential per body is indifferent. In practice, as we show in Section 5, directly stiffened affine systems significantly improve in performance over comparable systems that resolve rigid motion explicitly.

The stiff potential itself provides the effective constitutive model for the affine bodies – modulating their collision response upon impact (and so intrinsically handling restitution). While $V_{\perp}$ is a natural choice for close-to-rigid motion (Moser and Veselov, 1991), it is never exclusively suited. Indeed affine bodies could alternately be equipped with the ARAP energy (Igarashi et al., 2005; Alexa et al., 2000; Sorkine and Alexa, 2007), i.e., $\|\mathsf{A}-\mathsf{R}(\mathsf{A})\|^{2}_{F}$ for $\mathsf{A}=\mathsf{R}\mathsf{S}$ being the polar decomposition, to enforce rigidity, or else neo-Hookean, or any number of other rotation-invariant hyperelastic energies (Bonet and Wood, 1997). However, we find that the orthogonality potential is both effective and efficient. $V_{\perp}$ requires no expensive decompositions (as opposite to ARAP) – it can be computed as a polynomial

leading to more efficient evaluations of energy gradient and Hessian:

In comparison, energy operations for affine bodies with $V_{\perp}$ are over $43\%$ and $178\%$ faster than applying ARAP and neo-Hookean models respectively.

To ensure that every search fulfils the nonintersection guarantee, a significant (and for rigid bodies dominant) cost of each time step’s Newton solve is the repeated evaluation of CCD for primitive pairs in $\mathscr{C}$. Edge-edge or vertex-face pairs are defined on vertices $x_{j}\in\mathbb{R}^{3}$ for $j=1$ to $4$ belonging to one of two bodies in $\mathscr{B}$. At any iteration $\ell$ of a Newton solve, new positions of these vertices $x_{j}$ are proposed, per participating body $b$, as $\mathsf{A}_{b}^{\ell}\bar{x}_{j}+p_{b}^{\ell}$, and so we search along the direction formed by $\Delta\mathsf{A}_{b}^{\ell}=\mathsf{A}_{b}^{\ell}-\mathsf{A}_{b}^{\ell-1}$ and $\Delta p_{b}^{\ell}=p_{b}^{\ell}-p_{b}^{\ell-1}$. The corresponding trajectories are then just $(\mathsf{A}_{b}^{\ell-1}+\alpha\Delta\mathsf{A}_{b}^{\ell})\bar{x}_{j}+p_{b}^{\ell}+\alpha\Delta p_{b}^{\ell-1}$, with $\alpha\in[0,1]$. In other words, we are testing displacements (from $0$ to $1$) along $\Delta\mathsf{A}_{b}^{\ell}\bar{x}_{j}+\Delta p_{b}^{\ell}$ starting from $x_{j}^{\ell-1}=\mathsf{A}_{b}^{\ell-1}\bar{x}_{j}+p_{b}^{\ell-1}$ i.e., the previous iteration footstep and so can apply standard linear CCD rather than the expensive, curved CCD required for rigid-body trajectories. Throughout we employ Additive CCD (ACCD) method (Li et al., 2021b) for all affine body CCD evaluations.

The IPC framework inverts traditional contact-processing practices by integrating collision detection within every iterate, inside each nonlinear time step solve (rather than once per time step, after completing a solve – as has been a standard practice in rigid body pipelines). Above we have already seen the implications of this choice on CCD and the advantages for ABD. We next address culling of contact pair evaluations and then below, in the following section, an integrated approach to efficiently compute the local energy Hessian evaluations and their assembly.

With the affine bodies’ stiff potentials applied to model close-to-rigid motion, we eliminate, as in rigid body models, the need for self-collision processing on intra-body surface primitives pairs. To further cull inter-body surface pairs from downstream collision processing, Rigid body methods commonly employ precomputed per-body (or per unique mesh instance) BVH which only requires a rigid (or, in our setting, affine) transform to evaluate at each configuration update. To further improve collision detection queries, most popular rigid body libraries (e.g., Bullet, Mujoco, PhysX, Flex) additionally require that all nonconvex surface meshes be replaced by approximations – convex decompositions (Mamou et al., 2016) so that evaluations can utilize convexity assumptions.

In the context of stiff-body models we propose a simple strategy that significantly enhances contact-pair culling. We start with the observation that contact-dense regions between close-to-rigid meshes are most often local. Unlike deforming meshes (e.g., consider a cloth model), in common, contact-rich configurations we see that a large portion of each body’s surface remains free of contacts. A BVH evaluation over each body’s full surface is then a bit too aggressive (and a potentially unnecessarily expensive precomputation cost) – especially as we consider higher resolution models and/or increasing numbers of bodies.

We instead begin with a coarse AABB per body, offset to account for the small $\hat{d}$ parameter of the contact accuracy. While clearly not providing a tight bound we only need to evaluate potential contacts in each pairwise overlap between these initial, intersecting AABB (i-AABB) volume. With stiff affine transformations the overlap volume of each i-AABB is most often sufficiently small – the total number of surface primitives within our i-AABB is often one order smaller than that obtained from full-mesh BVH queries. Then, utilizing the guarantee that i-AABBs will not intersect one another, this already provides effective culling for us to directly proceed with parallel (multi-threaded) primitive-pair collision-check evaluations within each i-AABB.

In the most extreme, heavily entangled configurations, e.g., see Fig. 2, we find i-AABB volumes can be too conservative for CPU implementations. Here we simply build a shallow, three-level AABB hierarchy (a binary AABB tree) for each i-AABB. Our observation is that for GPU implementations and most (even contact-intense) CPU examples the single-level i-AABB is highly effective (and exceedingly simple to implement). We utilize the i-AABB hierarchies solely for CPU examples where dense, entangled contacts are consistently encountered, such as the geared system in Fig. 1 and the interlocked collision of complex models in Fig. 2. As a representative improvement we note that in the latter example the ship and the helicopter models have $227$K and $116$K surface triangles respectively; here the three-level hierarchical i-AABB is over $90\%$ more efficient in culling surface pairs in comparison to per-body BVH.

It is, by far, most efficient to simulate a set of noncontacting affine bodies. Along with reduced cost for collision detection, the global (system-wide) Hessian for $E$ is then simply a $12|\mathscr{B}|\times 12|\mathscr{B}|$ block-diagonal matrix with just a separable $12\times 12$ block for each affine body’s mass and orthogonality energy ($V_{\perp}$) contributions. When the system includes contacts, however, off-diagonal terms from active, inter-body contact and friction potentials necessarily pollute the Hessian to account for contact coupling.

Standard FEM-type evaluation and assembly would suggest iterating across all active contact potentials. However, here the surface mesh resolution, and so the corresponding number of surface pairs forming contact potentials, is generally much larger than the number of affine bodies we simulate. Traditional assembly, in this contact-oblivious way, is neither necessary nor efficient. We instead integrate our Hessian evaluation and assembly with the i-AABB hierarchy for contact-aware parallelization in both multi-core CPU and GPGPU implementations.

We start with the easiest part. The global Hessian’s default non-zero diagonal blocks are given by a constant mass term and the orthogonality potentials’ Hessian (Eq. (8)). This can be computed trivially in parallel. Next comes the expensive part. The Hessian of the contact potential, $V_{C}$, then has an unpredictable pattern which varies widely with contact states. Recall that, when bodies $i$ and $j$ are sufficiently close ($\leq\hat{d}$), barrier and friction forces between them activate, resulting in non-zero contributions to both their respective $i$-th and $j$-th $12\times 12$ diagonal blocks and to the off-diagonal blocks linking the corresponding body coordinates in the Hessian.

Here we observe that the configuration-varying sparsity of the global Hessian is effectively determined by our i-AABB culling. If the leaves of the i-AABB tree between bodies $i$ and $j$ are empty, they certainly do not contact, and the barrier and friction potentials make no contribution to the Hessian. Otherwise, we can conservatively allocate space for the corresponding $12\times 12$ off-diagonal blocks to store all potential (active contacts are not yet certain) non-zero contact Hessian contributions between bodies $i$ and $j$. Thus, utilizing the i-AABB structures we apply a two-pass strategy to compute and assemble barrier and friction terms for the global Hessian. The first pass iterates across all surface primitives pairs within each i-AABB; their local Hessians are computed in parallel and cached. Here we also account for the Hessian’s symmetry to reduce total memory consumption. Our second pass then accumulates the local Hessians from the surface pairs. Here, as each culled i-AABB is independent and non-intersecting, accumulation is parallelized at each corresponding non-zero element of the global Hessian.

We find that i-AABB Hessian construction is significantly faster than default sequential parallelization of Hessian computation (e.g., as applied in rigid-IPC). We observe speedups of up to two orders, especially when the contacting system is composed of large numbers of bodies. Here, for example, the simulation in Fig. 3 provides a representative example, with a $188\times$ speedup over contact-oblivious, sequential Hessian construction.

Both ABD and rigid-IPC (Ferguson et al., 2021) exploit the IPC (Li et al., 2020) model for contact processing and friction modeling. Here we carefully compare ABD with rigid-IPC across several representative simulation scenarios of varying complexities, as illustrated in Figs. 3 through 7.

The snapshots of the first comparison are given in Fig. 3. In this test, a heavy ball linked by a chain of metal rings dashes into a stack of wooden blocks, which are scattered by the collision. There are in total $575$ bodies in this example. Both rigid-IPC and ABD produce high-quality simulation results without any inter-penetration. However, ABD is $124\times$ faster. The performance of ABD is not sensitive to the time step size. Doubling or even quadrupling the time step size (from $1/100$ sec to $1/50$ sec and $1/25$ sec) lead to similar results and performance with ABD. Detailed timing statistics can be found in Section 5.5. Similar observations are received in chain net examples. For the small chain net case (Fig. 5), we have few bodies, and the relative velocities among bodies are small. In this “entry-level” test, ABD offers a $34\times$ speedup.

An “upgraded” experiment is reported in Fig. 4, where the resolution of the chain net is set to $16\times 16$. We also drop a ball, which bounces back and forth on the net triggering interesting dynamics. The collisions and contacts also become more complicated than the plain simple chain net in Fig. 5. The difference between our method and rigid-IPC becomes more significant. At some instances when sharp collisions occur under high relative velocities, rigid-IPC could take multiple hours to simulate one frame ($\Delta t=1/100$ sec). On the other hand, our method needs seconds at most. On average, ABD is over $1,200\times$ faster than rigid-IPC. This speedup exceeds $4,000\times$ from time to time during the simulation. In ABD, affine CCD, i-AABB, as well as the integrated Hessian assembly are all parallelization-ready. Therefore, ABD typically receives one more orders performance gain on CUDA. That makes ABD over $10,000\times$ faster than rigid-IPC in the example of Fig. 4.

The friction handling in ABD and rigid-IPC is similar as both follow the variational friction model originally proposed in (Li et al., 2020). Nevertheless, ABD still exhibits superior performance in simulations dominated by frictions. We hereby report three more experiments in Figs. 6 and 7. In those comparisons, both ABD and rigid-IPC use the same simulation settings with the time step size of $1/100$ sec. The screw example (Fig. 6 left) is relatively simple – the surface geometry only has $7$K triangles and $5$K edges. We rotate the bolt into the nut. On average, rigid-IPC needs $2.59$ sec to simulate one step, and our method only needs $87$ ms. The arch (Fig. 6 right) consists of $100$ brick blocks. In this example, we mainly have static friction to deal with. Rigid-IPC runs faster than the screw case and needs $0.67$ sec for simulating one frame. Our method however, only uses $8.7$ ms. The card of house example starts with a 10-level stack of $155$ cards. The stack is first held up under the frictions between cards and gets crashed by two falling boxes. In this example, it takes about $8.9$ sec for rigid-IPC to simulate one frame, and ABD needs $86$ ms.

Bullet is a widely used rigid body engine known for the simplicity and efficiency. Bullet is also based on the classic rigid body model formulation. However, Bullet uses the so-called convex collision resolution, which approximates the geometry of the body via convex decomposition proxies. The collision is resolved using impulse-based method at the convex proxies. Therefore, it is not surprising to see Bullet fails in simulations under lasting and intense collisions and contacts. We have shown that our method is orders of magnitude faster than rigid-IPC with further improved robustness. For smaller simulation problems, ABD is nearly as efficient as Bullet. It quickly outperforms Bullet from almost all perspectives as the complexity of the simulation escalates.

We compare ABD and Bullet using the wrecking ball setup (i.e., see Fig. 3) but under different block counts and time step sizes. In this set of comparisons, the ball is no longer attached to the chain as in Fig. 3. This is because under high-velocity movements, collision resolution in Bullet frequently fails, and the rings on the chain would “break out” so that the ball may not hit the stack. The results are reported in Fig. 8, where we have three configurations: light, medium, and full. The light test only has $16$ blocks on the stack; the medium test has $142$ blocks; and the full test is the same as in Fig. 3 with $560$ blocks.

The detailed timing comparison is listed in Tab. 1. In order to ensure the comparison is fair, we turn off multi-threading in ABD and Bullet. This is because our method processes collision on the surface triangles, which will receive more acceleration under parallelization. It is noted that under the light test ($\Delta t=1/100$ sec), Bullet is rather efficient and only needs $2$ ms to simulate one frame. Our method is slower than Bullet and takes $3$ ms to simulate one frame. This difference diminishes under a smaller time step. For instance, when the time step size is set as $1/240$ sec, which is the default setting in Bullet, the difference of per-frame simulation time is less than one millisecond ($2.2$ ms for ABD and $1.5$ ms for Bullet). If the time step size is more conservatively set to $1/1000$ sec, ABD becomes faster than Bullet by a small margin ($2$ ms for ABD and $3$ ms for Bullet). However, simulations using Bullet in all time step sizes have intersections (even with $\Delta t=1/1000$ sec). Our method on the other hand is guaranteed to be free of inter-penetration. Another interesting observation is per-frame simulation using ABD consistently becomes more efficient under smaller time steps, which is not the case for Bullet. This may be because a smaller time step could expose more collisions to Bullet solver, which are otherwise missed under a bigger step. In addition, Bullet fails all the friction experiments including the screwing bolt, arch, and house of cards (Figs. 6 and 7) even with highly conservative time step size ($\Delta t=1/10000$ sec).

The story is similar for the medium/full test: ABD ($92$ ms in medium and $657$ ms in full) is slightly slower than Bullet ($68$ ms in medium and $629$ ms in full) when $\Delta t=1/100$ sec²²2Basically, ABD and Bullet have the same FPS in the full test with $\Delta t=1/100$ sec. Yet ABD takes the lead under smaller time steps of $1/240$ sec and $1/1000$ sec. In addition, Bullet has significantly more inter-penetration instances in the medium/full tests than in the light test, which could end up with undesired artifacts in many simulation tasks.

As discussed, ABD’s collision processing directly handles each body’s input surface mesh boundary with all guarantees, including non-intersection, applying directly to those meshes. This means that control of resolution and so quality is available via standard, well optimized geometry pipelines, e.g. mesh decimation. Bullet, however, requires pre-processing to convert all input meshes body into convex decompositions. In turn Bullet will then simulate these decomposed, simplified models as proxies, resolving collision handling on them, for each rigid body’s surface rather than the original mesh. Bullet users often use V-HACD (Mamou et al., 2016), an automated library for computing convex decomposition, for this process. Utilizing libraries like V-HACD, can be a slow and often time-consuming iterative process to hand-tune quality parameters in order to obtain a reasonable a reasonable geometric approximation and resolution. Even so important features can be lost while details and symmetries are often unnecessarily broken. As such, many advanced users often resort to laborious hand-crafting of proxies when capturing important surface features is critical in an application (e.g., for design or robotics). Here, we demonstrate a comparative example of utilizing detailed, bone dragons in Fig 9 dropped on a highly featured geometry. After careful hand-tuning of V-HACD parameters, we present in Fig 9 right the resulting simulation using our best resulting V-HACD proxies. Here we see the convex decomposition geometry is still insufficiently detailed to capture local collision behavior between sharp asperities and convexities. In contrast, in Fig 9 left, ABD directly simulates the detailed geometry with tight tolerance so that all the original input meshes’ detailed affordances are kept. In turn we see that this allows the dragons’ concavities to tightly entangle and also slide directly onto and be caught by the sharp points – all while remaining intersection free.

An additional benefit of ABD is convenient constraint handling – especially for the those prescribing rotational DOFs. For instance, it is common to require multiple objects to obey a given kinematic relation following the joint connecting them. Such constraints are often nonlinear as they are formulated from rigid body rotational DOF. Enforcing them would then require computing the derivatives of the rotational DOFs (either via constrained Lagrangian methods or generalized coordinates). Here relaxation from rigidity to affinity also eases the processing of such constraints.

It is known that a non-degenerate tetrahedron uniquely defines an affine transform. This suggests our generalized coordinate $q$ can be mapped to any linear tetrahedron. From this perspective, an affine body simulation can also be viewed as a single-element FEM (with a simplified strain energy). The geometry of this element however, can be setup flexibly even its position deviates away significantly from the object. Let $\phi$ denote the map between $q$ and this virtual tetrahedron such that $q=\phi(\mathsf{P})$. $\mathsf{P}\in\mathbb{R}^{3\times 4}$ stores the deformed vertex positions (i.e., $p_{1}$, $p_{2}$, $p_{3}$, and $p_{4}$) of the element. Let the rest shape position of the element be $\bar{\mathsf{P}}$, and one can easily verify that:

Here, $\text{vec}(\cdot)$ denotes the vectorization of a matrix. Since $\bar{\mathsf{P}}$ is given, we could use $\text{vec}(\mathsf{P})$ as the new generalized coordinate of the system. Clearly $\phi$ is a linear function of $\mathsf{P}$. $\partial\phi/\partial p_{i}$ only depends on $\bar{\mathsf{P}}$. Therefore, the Jacobi of the system remains constant:

In practice, we exploit the fact that $\bar{\mathsf{P}}$ could be any tetrahedron to manipulate such constraints intuitively. For instance as shown in Fig. 12, the gear is constrained to rotate around a prescribed axis (blue dash line). In ABD, we map $q$ to a tetrahedron. Since the tetrahedron can be chosen freely, we make sure that one of its edges is parallel to the rotation axis. After that, the axis constraint can be simply posed as a linear equality constraint fixing two vertices of the edge ($p_{2}$ and $p_{4}$ in the figure). Indeed, such formulation is over-constraining as it eliminates six freedoms instead of two (we should also allow $p_{2}$ and $p_{4}$ to move laterally along the edge). However, such “over-constrainment” still gives correct simulation because ABD uses more DOFs to simulate a rigid object. Other types constraints can be dealt with in a similar way (e.g., the constraint over the rotation angle can be converted to a linear shearing constraint).

Two of such examples are reported in Fig. 10. The pendulum is a basic mechanism structure with two rigid links constrained by a hinge. The octopus has eight tentacles, and each of it consists of five joints. We also compared ABD with rigid-IPC with these models. ABD does not only enjoy a more handy formulation but also runs significantly faster ($371\times$ and $28\times$ speedups for the pendulum and octopus respectively). A more challenging experiment is reported in Fig. 1, where we simulate a more complex mechanism device of a set of 28 gears. Those gears are coupled by tooth contacts and shafts (see Fig. 14). There are totally $2,450$K triangles and $3,090$K edges on the surface (i.e., those could participate in culling and CCD). A component-wise breakdown showing the inter-connectivity of gears is given in the right of the figure. As the driving torque is applied at one of the gear (the one with red arrow), the entire device moves forward. This simulation is very demanding for the robustness of the simulation algorithm. Because gear tooth are close to each other, highly localized and detailed CCD is massively used. All rigid body simulation algorithms we have tested fail in this case including rigid-IPC, MuJuCo, and Bullet. As mentioned, Bullet requires building a convex proxy for collision processing, which does not capture the zigzag geometry at the gear tooth. MuJuCo fails the test even under highly conservative settings (e.g., very small time step size). Rigid-IPC also fails this experiment: after first few steps, the Newton iteration loops infinitely because the cured CCD is unable to find a usable TOI (time of impact). ABD manages to simulate this scene without any issues. Each frame takes about $12$ sec on the CPU. As the majority of the computation is at CCD processing, we port our i-AABB algorithm to CUDA, and the simulation of the gear set reaches an interactive rate (at 5 – 10 FPS).

Another benefit of ABD, when viewed as a reduced single-tetrahedron FEM, is to simulate models with both rigid and soft components. Such objects are vastly available in real world. While existing rigid body methods can also be augmented to incorporate such hybrid simulations, it is particularly effortless in the framework of ABD. We show a simulation of such scene. The barbarian ship has a rigid ship body and two deformable canvas. Each ship has $225$K surface triangles and $341$K edges. The helicopter is also hybrid with a rigid body, two soft rotors, and four soft wheels. There are $116$K triangles and $176$K edges on each helicopter. After two ships fall into the glass tank, five helicopters follow, producing interesting animation effects of both rigid and deformable dynamics. Note that the collision between rigid and soft bodies can be handled uniformly using barrier-based penalties. In such hybrid simulation, the Hessian of elastic potential on the deformable components are much smaller (by several orders) than the Hessian of the orthogonality potential ($V_{\perp}$), and the positive definiteness of the global system matrix is numerically compromised. Therefore, LLT factorization (SimplicialLLT) may fail occasionally. In this case, we switch to LDLT Cholesky for each Newton solve. A possible remedy for increase numerically stability is to use Schur complement like formulation (Peiret et al., 2019) to somehow decouple DOFs from rigid (affine) and deformable components.

Detailed time statistics of other the experiments is reported in Tab. 2. In all the experiments, we uniformly scale the scene to a $1\times 1\times 1$ box and set $\hat{d}$ as $1/1000$. In other words, if the size of the model is around one meter, the contact accuracy is guaranteed to be less than one millimeter, and all the models do not intersect with each other at any time during the simulation. We report comparative timing benchmark of ABD and rigid-IPC under $\Delta t=1/100$ in most experiments. Normally, rigid-IPC is numerically robust under larger time steps. However, its performance is highly sensitive to a bigger $\Delta t$. This is because the underlying curved CCD quickly becomes prohibitive when examining over a wider trajectory gap. Therefore, we do not report our speedup over rigid-IPC for any time steps bigger than $\Delta t=1/100$ sec. For chain net example shown in Fig. 4, the literal speedup exceeds four orders in general if we set $\Delta t=1/25$. It is not our intention to oversell just by tweaking the time step size.

The majority portion of our performance improvement is brought by relaxing the rigidity constraint. This can be reflected by two important metrics from Tab. 2: the iteration count ($\#$ Iter.) and the time needed for CCD processing (CCD). It is clearly shown that as long as the contact frequency in the simulation is intense, ABD always requires much fewer iterations to converge than rigid-IPC does. This difference is “extremized” in the gear set example (Fig. 1), where ABD needs $17$ iterations for each time step, and rigid-IPC needs infinite iterations. CCD processing is another major game changer. As mentioned, rigid-IPC strictly follows the rigidity constraint making per-step trajectory curved. This curved trajectory is split and converted back to piece-wise line segments again during the CCD. This conversion is fully avoided in ABD, and one can use any existing CCD algorithm to detect potential intersection between primitives and compute TOI. We would like to remind that CCD needs to be carried out at each Newton iteration in order to make sure the line search does not introduce inter-penetrations. Therefore, the performance gap between A-IPC and rigid-IPC is further scaled by the iteration count. Together, those two factors contribute to $90\%$ of our speedup in complex simulations such as the big chain net (Fig. 4), house of cards (Fig. 7) and wrecking ball (Fig. 3). In other relatively lighter simulation experiments, our improved culling and Hessian assembly become more profitable. Our culling is up to to $90\%$ more effective than conventional BVH-based strategy for models with complex geometry. It is on average $30\%$ more effective for simpler shapes (e.g., the chain net or the octopus). In terms of system solve, ABD is typically slower than existing rigid body methods due to inflated DOFs. However, this disadvantage is invisible because ABD always enjoys a much fewer iterations and much faster CCD processing. This is our core inspiration of designing ABD.

We suspect increasing the body count will grant the lead to rigid-IPC at a certain point and simulate a huge chain net model with $27,645$ bodies. The result is opposite – we estimate ABD has a $\sim 5,000\times$ speedup (on CPU). This is just a rough assessment as we are never able to finish rigid-IPC in this stress test, which will need several years while ABD finishes the simulation in five days.