arXiv CSopen access

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration Shih-Yu Lai1,2∗ , Sung-Han Tien1∗ , Jui-I Huang1 , Yen-Chen Tseng1 , Yi-Ting Chiu1 , Siyuan Luo3 , Ziqiu Zeng3 , Fan Shi3 , Peter Yichen Chen4 , Tiantian Liu5 , Yu-Lun Liu6 , Bing-Yu Chen1†∗ 1 National Taiwan University

arXiv:2605.14526v1 [cs.GR] 14 May 2026

4 The University of British Columbia

2 MoonShine Animation Studio 5 Independent Researcher

Homogeneous Model

Heterogeneous Model

DiffPD

DiffPhD

3 National University of Singapore

6 National Yang Ming Chiao Tung University

Unified Contact-Rich GPU-Acceleration Hard

Differentiable Simulation of Soft Body

Soft

Hard

center of mass target location

Hyperelasticity

Figure 1: We present DiffPhD, a GPU-accelerated differentiable Projective Dynamics (PD) solver for heterogeneous hyperelastic materials: 1. Heterogeneous Materials—stiffness-aware projective weights baked directly into the global matrix handle extreme stiffness contrasts (crab); 2. Differentiable Simulation—trust-region eigenvalue filtering, adaptively switching between absolutevalue and clamping projections, is lifted onto the proximal-map Hessian inside the backward pass to yield accurate Neo-Hookean gradients at high Poisson’s ratios and large deformations (bunny, armadillo); 3. Unified Contact-Rich GPU Acceleration—a single sparse inverse factor reused across forward, backward, and frictional-contact computations preserves gradient accuracy (robot manipulator). Gray arrows: applied forces; green vertices: loaded DoFs; dark purple regions are 100× stiffer than cyan.

Abstract Differentiable simulation of soft bodies is a foundation for system identification, trajectory optimization, and Real2Sim transfer. Yet, existing methods such as the differentiable Projective Dynamics (DiffPD) struggle when faced with heterogeneous materials with extreme stiffness contrasts, hyperelasticity under large deformations, and contact-rich interactions, which are common scenarios in the real world. We present DiffPhD, a unified GPU-accelerated differentiable Projective Dynamics framework for heterogeneous materials that tackles these intertwined challenges simultaneously. Our key insight is a careful integration of: (i) stiffness-aware projective weights to embed heterogeneity into the global system; (ii) trustregion eigenvalue filtering lifted to the backward pass for stable hyperelastic gradients and a type-II Anderson Acceleration scheme

with dual-gate convergence to stabilize forward iteration under large stiffness contrasts; and (iii) a unified GPU pipeline that reuses a single sparse factor across forward, backward, and contact computations, with stiffness-amplified Rayleigh damping folded into the same factor for heterogeneity-aware dissipation at zero recurring cost. DiffPhD achieves strict gradient accuracy while delivering up to an order-of-magnitude speedup over prior differentiable solvers on heterogeneous, hyperelastic, contact-rich benchmarks. Crucially, this speedup does not come at the cost of stability: DiffPhD remains convergent on stiffness contrasts up to 100× where prior PD solvers degrade. This unlocks end-to-end gradient-based optimization on regimes previously bottlenecked by either solver fragility or per-iteration cost—shell–joint composite creatures, soft characters wielding stiff weapons, and soft-gripper robotic manipulation—all handled within a single forward–backward pass.

† Corresponding author. Emails: [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected] ∗∗ Equal contribution.

Lai, et al.

INTRODUCTION

Differentiable simulation of soft bodies has emerged as a foundation for inverse problems in computer graphics and robotics, enabling system identification, trajectory optimization, and sim-to-real transfer through end-to-end gradient flow [Du et al. 2021; Li et al. 2022b; Ma et al. 2023; Nava et al. 2022; Zeng et al. 2025]. As applications push toward realistic scenarios—characters with bones embedded in soft tissue, contacting with rigid bodies, or robotic graspers manipulating heterogeneous objects—two ingredients become indispensable: (i) material heterogeneity, where Young’s modulus varies by orders of magnitude across a single mesh and the parameters are themselves often the design target; and (ii) hyperelasticity with frictional contact, capturing the large, non-linear, volume-preserving deformations of real-world flesh, rubber, and tissue. Heterogeneity is not a way to make a body stiffer—it is the structural property that makes spatially varying response, region-specific identification, and per-zone inverse design possible at all. Yet their combination is precisely where existing differentiable solvers break down. Why heterogeneity and hyperelasticity are hard for differentiable PD? Projective Dynamics (PD) [Bouaziz et al. 2014] and its differentiable extension DiffPD [Du et al. 2021] are attractive precisely because the global stiffness matrix 𝑨 is constant across iterations, admitting a one-time factorization that amortizes every forward and backward solve—an advantage that depends on 𝑨 being well-conditioned. When stiffness contrasts span 10× or more—routine in shell–joint composites and layered biomechanical models [Modi et al. 2021; Trusty et al. 2023]—the condition number of 𝑨 explodes, iterative solves stall, and gradients become unreliable. Compounding this, Neo-Hookean energies are non-convex: the local step requires eigenvalue projection, and the choice between absolute-value filtering [Chen et al. 2024b] and clamping [Teran et al. 2005] is itself problem-dependent under high Poisson’s ratios and large volume changes [Chen et al. 2024a]—a wrong choice destabilizes forward iteration, a fixed choice forfeits convergence speed across the state diversity a differentiable rollout traverses. Quasi-Newton variants [Cai et al. 2025] inherit the same pathology, and contact-rich rollouts force the backward to repeatedly apply 𝑨−1 against contact Jacobians—turning the factorization win into a per-iteration bottleneck. Existing tools each cover only a slice: FBA [Zeng et al. 2025] delivers a nested-dissection sparse inverse 𝑨−1 = 𝑺𝑇 𝑺 with unified Signorini–Coulomb NCP contact but is forward-only; the IPC-based differentiable solver [Huang et al. 2024] delivers end-to-end gradients at an order-of-magnitude wall-clock cost above PD; DiffPD itself retains gradient flow only in the homogeneous, corotated regime, falling back to a non-PD direct solver under Neo-Hookean. Why is integration the real difficulty? The components we draw on—trust-region filtering [Chen et al. 2024a], Type-II Anderson Acceleration [Peng et al. 2018], FBA’s sparse factor [Zeng et al. 2025]—were each designed under assumptions the others violate: trust-region filtering targets the forward element elastic Hessian and must be re-targeted onto the prox-map Hessian inverted by the backward IFT, with an adaptive rule avoiding per-step Hessian re-assembly and confined to the backward pass alone, since any forward prox-map modification corrupts the AA fixed-point; Anderson Acceleration assumes an autonomous fixed-point map

and uses a wide history window, both of which heterogeneous curvature breaks across regions, demanding bounded history and a convergence test immune to AA-induced plateaus; and FBA’s forward-only factor must serve both passes under a single sparsity pattern, accommodating contact Jacobians and a per-element prox-map differential absent from the forward NCP without refactorization that would destroy PD’s amortization. Resolving these tensions without sacrificing convergence speed, gradient accuracy, or factor reuse is the core issue. Our approach. We present DiffPhD, a unified GPU-accelerated differentiable solver for heterogeneous, hyperelastic, contact-rich elastodynamics within a single gradient-consistent pipeline (Fig. 1). Our design rests on one observation: PD’s matrix 𝑨 is the shared backbone of every expensive operation—forward global step, Delassus contact compliance, backward adjoint—so we organize the solver around making 𝑨 cheap to invert and robust under contrast. (i) We encode heterogeneity structurally at assembly through stiffness-aware projective weights, so contrast lives inside 𝑨 rather than in any per-element local solve. (ii) Trust-region filtering is lifted onto the prox-map Hessian in the backward pass with a stateadaptive rule whose quadratic model reuses 𝑨 at zero extra cost; the forward fixed-point is stabilized by bounded-history Type-II AA and a dual-gate convergence rule. (iii) FBA’s nested-dissection pipeline [Zeng et al. 2025] is lifted into a single GPU-resident loop sharing one persistent factor across forward, backward, and contact stages. Our work makes three core contributions:

• Heterogeneity via Stiffness-Aware Projective Assembly. Stable, gradient-consistent simulation across ordersof-magnitude stiffness contrasts—surpassing FEM, Newton, and Mixed-FEM [Trusty et al. 2023] baselines in forward speed and backward conditioning. • Differentiable Hyperelasticity via Proximal-Map TrustRegion Filtering. The first lift of trust-region filtering onto PD’s prox-map Hessian in the backward pass, with closedform differential and a state-adaptive rule reusing the forward PD matrix as its quadratic model—complemented by bounded-history Type-II Anderson Acceleration and a dualgate convergence criterion—yielding stable Neo-Hookean gradients beyond regimes prior differentiable PD covers. • Unified Contact-Rich GPU Forward/Backward Loop. A single persistent sparse-inverse factor shared across forward, backward, and contact stages—an order-of-magnitude wall-clock improvement over prior end-to-end-differentiable contact solvers on contact-rich heterogeneous Neo-Hookean rollouts.

We validate DiffPhD on heterogeneous, hyperelastic, contactrich benchmarks—spanning system identification, trajectory and material-design optimization, and Real2Sim manipulation—where prior differentiable solvers cannot stably handle, demonstrating strict gradient accuracy and substantial speedups over state-of-theart baselines.

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

2 RELATED WORK 2.1 Differentiable Simulation of Deformable Solids Differentiable physics splits along the time-integration scheme (explicit vs. implicit) and the spatial discretization (particle, FEM, or learned). Explicit schemes pair MPM with autodiff (ChainQueen [Hu et al. 2019], DiffTaichi [Hu et al. 2020]) but require tiny timesteps and aggressive checkpointing for long rollouts; meshless alternatives include projective peridynamics [Lu et al. 2024] and Simplicits [Modi et al. 2024]. Implicit solvers obtain gradients via the adjoint method on a linearised dynamics—used for Real2Sim viscoelastic identification [Hahn et al. 2019], multi-body frictional contact (ADD [Geilinger et al. 2020]), and differentiable cloth [Li et al. 2022a; Liang et al. 2019; Qiao et al. 2020]—where the linear solve dominates per-step cost. Neural surrogates [Li et al. 2019b; Ma et al. 2023; Sanchez-Gonzalez et al. 2020] bypass the solver, trading accuracy for differentiability by construction. DiffPD [Du et al. 2021] occupies a distinct middle ground: implicitFEM accuracy with the prefactorised linear backbone of Projective Dynamics [Bouaziz et al. 2014; Liu et al. 2017] reused in the backward pass, extended to aquatic locomotion with neural hydrodynamics [Nava et al. 2022] and contact-point discovery for soft-body manipulation [Li et al. 2022b]. DiffQN [Cai et al. 2025] recasts the same elastodynamics in a quasi-Newton framework with lowrank Hessian updates, cheaper per step but with no remedy for Neo-Hookean indefiniteness. Adjacent soft-robotics efforts span gradient-free VoxCAD [Hiller and Lipson 2014], modal SOFA [Allard et al. 2007], learning-in-the-loop co-design [Spielberg et al. 2019], biomimetic swimmers [Min et al. 2019], fabricated foam quadrupeds [Bern et al. 2019], reduced-order soft-character control [Barbič and James 2005; Thieffry et al. 2019], dexterous volumetric manipulation [Lee et al. 2018], evolutionary design [Cheney et al. 2013; Corucci et al. 2016a,b], and differentiable IPC [Huang et al. 2024]. A complementary line compresses the deformation space rather than the solver: subspace integration [Barbič and James 2005], modal actuation [Benchekroun et al. 2024], and adaptive subspaces for contact-rich elastodynamics [Trusty et al. 2024] with stochastic force-dual bases [Benchekroun et al. 2025]. Also, Subspace MixedFEM [Trusty et al. 2023] is closest in spirit to our heterogeneity contribution but is forward-only; gradient flow under simultaneous heterogeneity, Neo-Hookean stiffening, and frictional contact has not been established by any of the above, which our pipeline addresses end-to-end.

2.2

Hyperelastic Energies and Hessian Stabilization

Hyperelastic models [Kumar and Rao 2016; Mooney 1940; Ogden 1984; Rivlin 1947] achieve visual realism through non-convex potentials whose element-wise Hessians are routinely indefinite under large strain. Stable Neo-Hookean [Smith et al. 2018a] addresses rest-state stability and inversion robustness, with adaptations to dynamic-production variants [Kim and Eberle 2020], Cauchy–Green ARAP [Lin et al. 2022], edge-based SVK [Kikuuwe et al. 2009], and fibre-aligned anisotropic muscles [Modi et al. 2021].

Projected-Newton solvers restore convergence by replacing the indefinite Hessian with a PSD surrogate. Per-element eigenvalue clamping [Teran et al. 2005] sets non-positive eigenvalues to a floor; absolute-value filtering [Chen et al. 2024b; Gill et al. 2019; Nocedal and Wright 2006] flips their sign instead, dramatically more efficient on stiff Neo-Hookean materials but prone to over-damping near the optimum. Diagonal additions [Fu and Liu 2016] and Tikhonov schemes [Paternain et al. 2019] share this family but introduce problem-specific thresholds; projection-on-demand [Longva et al. 2023] sidesteps eigenanalysis by inflating the mass term until Cholesky succeeds. Beyond Newton-type schemes, second-order stencil descent [Lan et al. 2023] and preconditioned NCG [Shen et al. 2024] exploit locality and matrix-free structure for real-time interiorpoint hyperelasticity. Our adaptive eigenvalue filter follows trustregion and regularized-Newton theory [Chen et al. 2024a; Kanzow and Steck 2023; MORE 1993; Moré and Sorensen 1983; Pong and Wolkowicz 2014; Sorensen 1982; Ueda and Yamashita 2014], conceptually related to saddle-free Newton [Dauphin et al. 2014]. ADMM [Narain et al. 2016; Overby et al. 2017] and position-based [Macklin et al. 2016; Müller et al. 2007] reformulations offer complementary optimization views but do not address gradient flow via the per-element proximal map—the property our hyperelasticity contributes.

2.3

Linear Solvers, Preconditioning, and Frictional Contact

For large-scale elastodynamics, the global linear solve dominates wall-clock cost. Within PD, acceleration strategies include Chebyshev semi-iteration [Wang 2015], GPU descent [Wang and Yang 2016], parallel Gauss–Seidel with randomized graph coloring [Fratarcangeli et al. 2016], reduced-subspace formulations [Brandt et al. 2018], and a domain-decomposed CPU cloth solver [Lu et al. 2025]. For the SPD system, multigrid—geometric [Dick et al. 2011; Wang et al. 2018; Xian et al. 2019; Zhu et al. 2010] and algebraic [Naumov et al. 2015; Tamstorf et al. 2015]—is the standard tool, joined by multilevel preconditioners for Laplacians [Krishnan et al. 2013] and illconditioned matrices via multiscale Cholesky [Chen et al. 2021]. Direct factorizations stay competitive when topology is fixed [Herholz and Sorkine-Hornung 2020], with adaptive algebraic reuse [Zarebavani et al. 2025] for dynamic sparsity; both scale poorly with problem size and re-meshing. Domain decomposition offers an orthogonal axis. In Additive Schwarz [Cai and Saad 1996; Dryja and Widlund 1989; Frommer and Szyld 1999; Schwarz 1870] and its multilevel extensions [Dolean et al. 2015; Dryja and Widlund 1990, 1992; Zhang 1992] add coarsespace corrections for inter-domain coupling; restricted variants [Cai and Sarkis 1999], balancing decomposition [Mandel 1993], FETI [Farhat et al. 2001; Farhat and Roux 1991], and Schur-complement preconditioners [Haase et al. 1991; Li and Saad 2017] occupy related niches. In graphics, decomposition has driven reduced-space character simulation [Barbič et al. 2012; Kim and James 2011; Wu et al. 2015; Yang et al. 2013], decomposed implicit elastodynamics [Li et al. 2019a], discontinuous-Galerkin coarse grids [Edwards and Bridson 2015], and high-resolution cloth wrinkles [Wang 2021]. The GPU MAS preconditioner [Wu et al. 2022] shows that small non-overlapping subdomains with Nicolaides-style coarse spaces [Nicolaides 1987] outperform multigrid for cloth and deformables; spectral coarse

Lai, et al.

spaces [Efendiev et al. 2012; Spillane and Rixen 2012] could improve heterogeneity handling but are hard to assemble on-the-fly. Frictional contact divides into two camps. Penalty formulations [Bridson et al. 2002; Macklin et al. 2020; Wu et al. 2020] introduce a fictitious repulsive energy with per-scene stiffness tuning. Complementarity formulations enforce non-penetration and Coulomb friction directly: IPC [Li et al. 2020] and its non-distance-barrier GPU extension [Lan et al. 2024] provide intersection- and inversionfree dynamics through a smoothed barrier, forming the basis for the differentiable extension [Huang et al. 2024] and the unified geometric-contact potential [Huang et al. 2025]; [Ly et al. 2020] integrate Signorini–Coulomb into PD’s local-global loop for cloth, and [Zeng et al. 2025] extend this to GPU-resident hyperelastic solids via a unified NCP. In-timestep remeshing [Ferguson et al. 2023] reformulates contact-rich elastodynamics around adaptive resolution rather than solver design; dynamics-aware coarsening [Chen et al. 2017] and parallel spatial-hashing collision detection [Tang et al. 2018] address adjacent concerns. Yet none co-design contact with material heterogeneity and end-to-end differentiability—the gap our unified pipeline (Sec. 4) closes.

3 BACKGROUND 3.1 Problem Formulation A deformable solid with 𝑛 𝑣 vertices, 𝑛𝑒 finite elements, and spatial dimension 𝑑 ∈ {2, 3} has generalised position 𝒒 ∈ R𝑑𝑛 𝑣 and velocity ¤ Each element carries an independent Young’s modulus 𝐸𝑒 𝒗 = 𝒒. and a shared Poisson ratio 𝜈: 𝜇𝑒 =

𝐸𝑒 , 2(1 + 𝜈)

𝜆𝑒 =

𝐸𝑒 𝜈 . (1 + 𝜈)(1 − 2𝜈)

(1)

The solid interacts with 𝐾𝑛 unilateral contacts, 𝐾𝑏 bilateral attachments, and 𝐾 𝑓 frictional pairs. Implicit Euler gives 𝑴 ˜ = 𝒇ela (𝒒𝑡 +ℎ ) − (𝛼𝑴 + 𝑩 𝜷 )𝒗𝑡 +ℎ + 𝑱𝑛𝑇 𝝀𝑛 + 𝑱𝑏𝑇 𝝀𝑏 + 𝑱 𝑓𝑇 𝝀 𝑓 , (𝒒𝑡 +ℎ − 𝒒) ℎ2 (2) P with 𝑩 𝜷 = 𝑒 𝛽𝑒 𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 the heterogeneity-weighted stiffness-damping operator (Sec. 4.5) and inertial free-fall target 𝒒˜ = 𝒒𝑡 + ℎ𝒗𝑡 + ℎ 2 𝑴 −1 (𝒇ext + 𝒇state (𝒒𝑡 , 𝒗𝑡 )).

(3)

The contact-free case is the stationarity condition of the primal objective [Bouaziz et al. 2014; Martin et al. 2011] Φ(𝒒) =

1 ˜ 2𝑴 + Ψ(𝒒), ∥𝒒 − 𝒒∥ 2ℎ 2

Ψ(𝒒) =

∑︁

𝑉𝑒𝜓𝑒 (𝑭𝑒 (𝒒)),

(4)

𝑒

analysed for its forward/backward implications in Sec. 3.3.

3.2

Heterogeneous Hyperelastic Energy

We adopt a compressible Neo-Hookean energy density at element 𝑒, 𝜇𝑒 𝜆𝑒 ∥𝑭𝑒 ∥ 2𝐹 − 𝑑 − 𝜇𝑒 ln 𝐽𝑒 + (ln 𝐽𝑒 )2, (5) 𝜓𝑒 (𝑭𝑒 ) = 2 2 chosen because (i) the ln 𝐽𝑒 barrier resists inversion (𝐽𝑒 → 0+ ⇒ 𝜓𝑒 → +∞); (ii) it reduces to linear elasticity in the small-strain limit, preserving the shear/bulk interpretation of (𝜇𝑒 , 𝜆𝑒 ); and (iii) it is 𝐶 2 on 𝐽𝑒 > 0, supporting closed-form derivatives required by the

Table 1: Principal notation. Symbol

Meaning

𝑴 ∈ R𝑑𝑛 𝑣 ×𝑑𝑛 𝑣 ℎ 𝑨 ∈ R𝑑𝑛 𝑣 ×𝑑𝑛 𝑣 2 𝑮𝑒 ∈ R𝑑 ×𝑑𝑛 𝑣 𝑤𝑒 , 𝑉𝑒 2 𝒑𝑒 ∈ R𝑑 M𝑒 𝜇𝑒 , 𝜆𝑒 ¯ 𝑘¯ ¯ 𝜆, 𝜇, 𝜇 ref = max𝑒 𝜇𝑒 𝛼 𝛽𝑒 𝑭𝑒 ∈ R𝑑 ×𝑑 𝐽𝑒 = det 𝑭𝑒 𝝈𝐹 , 𝝈 ∗ 𝑯𝜓 (𝝈 ) 𝑛𝑒𝑣 𝑯𝑒 prox 𝑯𝑒 𝜏 ∈ [0, 1] 𝜌 𝜀 TR ∆𝒒∗ 𝑚 ∆𝑮, ∆𝑸 𝜸∗ 𝑷 𝑺 𝑱𝑛 , 𝑱𝑏 , 𝑱 𝑓 𝝀𝑛 , 𝝀𝑏 , 𝝀 𝑓 𝑾𝑎𝑏 = 𝑱𝑎 𝑨−1 𝑱𝑏𝑇 Ω𝑛 , Ω 𝑓 , 𝑬𝑛 , 𝑬 𝑓 𝛿𝑛,𝑐 , 𝛿 𝑓 ,𝑐 , 𝜇 𝑗 Φ(𝒒)

lumped mass matrix time-step size PD global stiffness matrix (SPD, het+damp) element deformation-gradient operator per-element PD surrogate weight (𝑤𝑒 ∝ 𝜇𝑒 ), rest volume PD local auxiliary (projection) variable constraint manifold of element 𝑒 per-element Lamé parameters mesh-wide scalar means used in local prox damping normalisation reference scalar Rayleigh mass-damping coefficient per-element stiffness-damping coefficient (∝ 𝜇𝑒 ) deformation gradient of element 𝑒 volume-change ratio singular values of 𝑭𝑒 , converged prox stretches NH stretch-space Hessian, Eq. (15) nodes per element (𝑛𝑒𝑣 = 4 for tetrahedra, 𝑑 = 3) raw element elastic Hessian, ∇𝒒2 (𝑉𝑒 𝜓𝑒 ) ¯ PD prox-map Hessian, 𝑯𝜓 + 𝑘𝑰 backward TR blend parameter TR ratio, Eq. (28) TR tolerance (= 0.1) last converged fixed-point increment, Eq. (26) AA window size (default 𝑚 = 1 on het, 𝑚 = 5 on homo) AA residual / iterate difference matrices AA mixing coefficients METIS nested-dissection permutation FBA sparse inverse factor: 𝑨−1 = 𝑺𝑇 𝑺 contact / bilateral / friction Jacobians corresponding multipliers Delassus compliance sub-block NCP weight / regulariser matrices normal gap, tangential slip, Coulomb coefficient primal objective, Eq. (4)

forward Newton solve and the backward IFT. With elastic force 𝒇ela (𝒒) = −∇𝒒 Ψ, the first Piola–Kirchhoff stress is 𝑷𝑒 (𝑭𝑒 ) = 𝜇𝑒 𝑭𝑒 − 𝑭𝑒−𝑇 + 𝜆𝑒 ln(𝐽𝑒 ) 𝑭𝑒−𝑇 . (6) Heterogeneity sets the spectral conditioning. In a single-material ¯ the tangent-stiffness spectrum is governed ¯ 𝜆), body (𝜇𝑒 , 𝜆𝑒 ) ≡ (𝜇, by mesh quality alone; in our setting 𝜇𝑒 varies over orders of magnitude—composite shells with embedded rigid plates [Trusty et al. 2023], layered foams under impact, and shape-changing soft robots with stiff actuators [Modi et al. 2021] all exhibit 𝜇 max /𝜇min ≥ 100× (Fig. 2a). Two consequences follow: (C1) Tangent-stiffness conditioning of ∇2 Ψ is determined by 𝜇 max /𝜇min rather than mesh resolution (Fig. 2b), which motivates routing 𝜇𝑒 into PD’s projective weights 𝑤𝑒 so that 𝑨 inherits stiffness-aware preconditioning at zero per-step cost (Sec. 4.2); and (C2) Local non-convexity is elementlocalised (Fig. 2c): soft elements near 𝜎𝑖 → 0 produce indefinite per-element Hessians (Eq. (15)) while stiff neighbours remain convex, so uniform Hessian modification [Chen et al. 2024b; Teran et al. 2005] under-serves both regimes—motivating a state-adaptive trust-region filter [Chen et al. 2024a] on the per-element prox-map

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

(a) Heterogeneous Neo-Hookean Energy

rest shape

per-element

barrier

The reassembled 𝑯𝑒+ = 𝑼𝑒 diag(𝜅𝑘+ ) 𝑼𝑒𝑇 enters

deformed shape

𝑴 ∑︁ + 𝑯¯ (𝒒) = 2 + 𝑯𝑒 . ℎ 𝑒

Stretch

soft

stiff

inversion-resistant

twist

(b) C1: Stiffness Contrast Sets Global Conditioning route

into

PD weights stiffness-aware

state-adaptive trust-region filter

Figure 2: Two consequences of material heterogeneity in NeoHookean PD. (a) Heterogeneous Neo-Hookean energy. A bar with stiff (red) and soft (blue) regions deforms qualitatively differently under stretch and twist, and the energy’s logarithmic barrier prevents element inversion. (b) Consequence C1 (global). The spectrum of the global tangent stiffness widens directly with the stiffness contrast 𝜇 max /𝜇 min ; absorbing this contrast into the assembled operator 𝑨 keeps its inversion robust (Sec. 4.2). (c) Consequence C2 (local). Non-convexity is element-localised: a soft element near collapse develops an indefinite local Hessian, while its stiff neighbour stays positive definite. A uniform per-element filter under-serves both regimes, motivating our state-adaptive trust-region filter applied only in the backward pass (Sec. 4.4).

Hessian in the differentiable backward pass, adapting to local nonconvexity without breaking PD’s global fixed-point contraction (Sec. 4.4).

Variational Form and Indefinite Hessian Challenge

Per-element Hessian filtering for SPD restoration. Projected Newton replaces each 𝑯𝑒 by an SPD surrogate built from its eigendecomposition 𝑯𝑒 = 𝑼𝑒 Λ𝑒 𝑼𝑒𝑇 with Λ𝑒 = diag(𝜅 1, . . . , 𝜅𝑑𝑛𝑒𝑣 ). Two element-wise filters dominate, eigenvalue clamping [Teran et al. 2005] and absolute-value filtering [Chen et al. 2024b]:

absolute-value: 𝜅𝑘+ = |𝜅𝑘 |.

Projective Dynamics (PD) [Bouaziz et al. 2014; Liu et al. 2017] replaces the non-convex potential in Eq. (4) by a separable quadratic 2 surrogate. For each element 𝑒, an auxiliary variable 𝒑𝑒 ∈ R𝑑 is constrained to a manifold M𝑒 encoding the rotational (corotated) or stretch (Neo-Hookean) component of the deformation; the surrogate energy ∑︁ 1 ˜ 2𝑴 + ˜ {𝒑𝑒 }) = 2 ∥𝒒 − 𝒒∥ 𝑔(𝒒, 𝑤𝑒 𝑉𝑒 ∥𝑮𝑒 𝒒 − 𝒑𝑒 ∥ 2 (10) 2ℎ 𝑒 is minimised by alternating a local step—per-element projection 𝒑𝑒∗ = arg min𝒑 ∈ M𝑒 ∥𝑮𝑒 𝒒 − 𝒑∥ 2 , embarrassingly parallel—with a global step that fixes {𝒑𝑒∗ } and minimises in 𝒒, yielding the SPD linear system 𝑨𝒒 = 𝒃 with ∑︁ 𝑴 𝑴 ∑︁ 𝑤𝑒 𝑉𝑒 𝑮𝑒𝑇 𝒑𝑒∗ . (11) 𝑤𝑒 𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 , 𝒃 = 2 𝒒˜ + 𝑨 = 2 + ℎ ℎ 𝑒 𝑒 Because 𝑨 is independent of 𝒒, its factorisation is amortised across all forward iterations of a timestep and—as DiffPD [Du et al. 2021] exploits and we extend—across the differentiable backward pass (Secs. 4.8–4.6). Per-element PD surrogate weight. The weight 𝑤𝑒 in Eq. (11) matches the tangent modulus of 𝜓𝑒 at the rest configuration,

External and state-dependent forces 𝒇ext, 𝒇state are absorbed into 𝒒˜ P via Eq. (3), so Φ in Eq. (4) depends only on Ψ(𝒒) = 𝑒 𝑉𝑒𝜓𝑒 (𝑭𝑒 (𝒒)), which is generally non-convex under Neo-Hookean hyperelasticity. Newton-type solvers therefore confront an indefinite Hessian ∇2 Φ = 𝑴/ℎ 2 + ∇2 Ψ, whose per-element block 𝑯𝑒 = ∇𝒒2 (𝑉𝑒𝜓𝑒 ) becomes indefinite under high Poisson’s ratio or large volume change— the canonical non-convex regime of Eq. (5).

clamping: 𝜅𝑘+ = max(𝜅𝑘 , 0),

Clamping matches the local quadratic well near convexity but stalls under stiff non-convexity; absolute filtering preserves negativeeigen curvature information and is orders of magnitude faster on stiff Neo-Hookean instances, yet mildly damps convergence in benign regimes [Chen et al. 2024a,b]. The optimal filter is thus problem- and state-dependent. The same indefiniteness governs the adjoint linear system in the backward pass (Sec. 4.8); we resolve the filter ambiguity through a state-adaptive trust-region rule [Chen et al. 2024a] lifted into the backward pass (Sec. 4.4).

4 METHOD 4.1 Projective Dynamics: A Local-Global Surrogate

Per-Elelement Backward Pass

3.3

(9)

(7) (8)

𝑤𝑒 = 𝑘¯ (𝜇𝑒 , 𝜆𝑒 , 𝜀𝜎 ) ,

(12)

where 𝑘¯ fits a linear spring to 𝜕𝑷𝑒 /𝜕𝑭𝑒 over a singular-value range 𝜀𝜎 around the identity ([Du et al. 2021], Eq. 7). Since 𝑤𝑒 ∝ 𝜇𝑒 , stiffer elements receive proportionally larger surrogate stiffness, so 𝑨 naturally admits per-element material parameters at assembly time—a routing capability we extend in Sec. 4.2. Local step — corotated. M𝑒 = SO(𝑑) and 𝒑𝑒∗ = 𝑹𝑒 from the polar decomposition 𝑭𝑒 = 𝑹𝑒 Σ𝑒 . The prox-map Jacobian 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 is the polar-decomposition differential, intrinsically SPD on SO(𝑑). Local step — Neo-Hookean proximal operator. M𝑒 = R𝑑 ×𝑑 , and ¯ the local step is the proximal map of 𝜓𝑒 with penalty 𝑘, 𝒑𝑒∗ = 𝑎𝑟𝑔 𝑚𝑖𝑛 𝒑 ∈𝑅𝑑 ×𝑑

𝑘¯ ∥𝑭𝑒 − 𝒑∥ 2𝐹 + 𝑉𝑒 𝜓𝑒 (𝒑). 2

(13)

Lai, et al.

The SVD 𝑭𝑒 = 𝑼 diag(𝝈𝐹 )𝑽 𝑇 decouples the minimisation into a 𝑑-dimensional system over the principal stretches 𝝈 ∗ ∈ R𝑑>0 : ¯ )𝝈 ∗ = 𝑘¯ 𝝈𝐹 , (𝑯𝜓 (𝝈 ∗ ) + 𝑘𝑰

(14)

with stretch-space Hessian 𝜆¯ (𝑖 ̸= 𝑗). 𝜎𝑖 𝜎 𝑗 (15) Eq. (14) is solved by Newton iteration with backtracking line search and positivity enforcement; the converged projection is 𝒑𝑒∗ = 𝑼 diag(𝝈 ∗ ) 𝑽 𝑇 . ¯ is also the implicit operator that the The prox-map Hessian 𝑯𝜓 + 𝑘𝑰 differentiable backward pass must invert through the IFT; its potential indefiniteness near element inversion motivates the trust-region treatment of Sec. 4.4. Q 1 − ln 𝑗 𝜎 𝑗 [𝑯𝜓 ] 𝑖𝑖 = 𝜇¯ 1 + 𝜎𝑖−4 + 𝜆¯ , 𝜎𝑖2

[𝑯𝜓 ] 𝑖 𝑗 =

Volume-barrier surrogate. The default heterogeneous-NH PD assembly decomposes the projection as 𝒑𝑒∗ = (𝑹𝑒 , 𝒑𝑒vol,∗ ) with a corotated rotation step and a quadratic (𝐽𝑒 −1)2 volume step. For workloads where inversion (𝐽𝑒 → 0+ ) is a concern, we additionally provide a logarithmic volume barrier !2 ∑︁ 𝜆𝑒 ∑︁ log ln 𝜎 𝑗 , 𝜙𝑒 (𝝈) = −𝜇𝑒 ln 𝜎𝑖 + (16) 2 𝑗 𝑖 solved at the principal stretches of 𝑭𝑒 via the same Newton-onstretches step (Eq. (14) with 𝑯𝜓 replaced by the analytical Hessian log of 𝜙𝑒 ) and registered with weight 𝑤𝑒 = 𝜆𝑒 . The barrier diverges as 𝐽 → 0+ , preventing inversion; dropping the deviatoric 𝐼 1 term (present in the full Neo-Hookean 𝜓𝑒 ) lets it compose cleanly with the corotated rotation step without double-counting. This variant is opt-in via a material flag and does not alter the default assembly.

4.2

AA extrapolates oscillation

Figure 4: Stiffness-aware projective assembly. (a) Route heterogeneity through 𝑨, not the local prox. Per-element material parameters enter the global PD operator 𝑨 through projective weights 𝑤𝑒 ∝ 𝜇𝑒 (Eq. (12)); the local prox sees only mesh-wide scalar means. (b) Why per-element local prox fails. Mixing soft and stiff elements in the local Newton makes its prox-map Hessian element-varying: 𝒑𝑒∗ oscillates across PD iterations, and Anderson Acceleration amplifies the drift into meshlevel divergence. (c) 𝑨-only routing preserves contraction. Stiff weights enter 𝑨 and widen its spectrum with contrast, but every element still runs the same Newton-on-stretches; the PD fixed-point therefore stays contractive, leaving only a nonlinear residual for iteration to resolve.

Heterogeneity via Stiffness-Aware Projective Assembly

Where heterogeneity enters: in 𝑨 only, not in the local prox. For a heterogeneous mesh (𝜇𝑒 varying over orders of magnitude), two routings are conceivable: place per-element (𝜇𝑒 , 𝜆𝑒 ) inside Eq. (14) ¯ 𝑘¯ with 𝜇𝑒 , 𝜆𝑒 , 𝑘¯𝑒 ), or inject 𝜇𝑒 only through the pro¯ 𝜆, (replacing 𝜇, jective weights 𝑤𝑒 in Eq. (11). We adopt the second exclusively (Fig. 4a): 𝑤𝑒 = 𝑘¯ (𝜇𝑒 , 𝜆𝑒 , 𝜀𝜎 ) ∝ 𝜇𝑒 , (17) ¯ as in Eq. (12); the local prox Eq. (14) uses the mesh-wide with 𝑘(·) ¯ 𝑘¯ identically across elements. ¯ 𝜆, means 𝜇, Why scalar means in the local prox. PD’s local-global iteration is contractive iff 𝜌(𝑰 − 𝑨−1 𝑲pd ) < 1, with 𝑲pd the local-step stiffness from the prox-map Hessian. Per-element (𝜇𝑒 , 𝜆𝑒 ) in the local step ¯ element-varying, and any inversion-prone element makes 𝑯𝜓 + 𝑘𝑰 (small 𝜎𝑖 ) drives the local Newton into an indefinite regime that desynchronises with 𝑨’s constant operator: 𝒑𝑒∗ ping-pongs across PD iterations, which Type-II Anderson Acceleration (Sec. 4.3) extrapolates into mesh-level oscillation (Fig. 4b). Routing 𝜇𝑒 through 𝑤𝑒 alone preserves the cross-element contraction (the local map is the same stretch-space Newton in every element), while 𝑨 automatically encodes heterogeneity at assembly time since 𝑤𝑒 ∝ 𝜇𝑒 (Fig. 4c).

Fixed-point contraction without the elastic force. The PD right-hand side 𝒃 in Eq. (11) intentionally excludes 𝒇ela (𝒒𝑘 ). If included, the fixed-point map 𝑇 (𝒒) = 𝑨−1 (𝒃(𝒒) + 𝒇ela (𝒒)) has Lipschitz constant ∥𝑨−1 (∇𝒒 𝒃+∇𝒒𝒇ela )∥ 2 ; for stiff Neo-Hookean materials, the spectrum of ∇𝒒𝒇ela dominates 𝜆min (𝑨) and drives it above one. The choice 𝑤𝑒 ∝ 𝜇𝑒 ensures 𝑨 absorbs the dominant linear stiffness; the nonlinear residual is resolved by iteration. This design follows Algorithm 4 of [Zeng et al. 2025] and is essential when stiffness varies over orders of magnitude.

Figure 3: Coupled failure of heterogeneity and hyperelasticity. A bi-material Armadillo (soft upper body, stiff lower body) under twist: wrist drift originates from per-element prox-map Hessians turning indefinite at soft/stiff boundaries, and the resulting 𝝈 ∗ ping-pong is amplified by Anderson Acceleration into mesh-level oscillation. Addressed by bounded-window AA (Sec. 4.3) and backward trust-region filtering (Sec. 4.4).

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

4.3 Forward: Bounded Anderson Acceleration + Dual-Gate The PD fixed-point 𝒒𝑘+1 = 𝑇 (𝒒𝑘 ) converges linearly at rate 𝜌(𝑰 − 𝑨−1 𝑲pd ), which is close to one for stiff Neo-Hookean materials. We accelerate it with Type-II Anderson Acceleration [Peng et al. 2018; Walker and Ni 2011] with bounded window 𝑚, paired with a dual-gate convergence criterion on position and projection-residual stationarity. Bounded-window AA.. With iterate step 𝒈𝑘 = 𝒒𝑘 − 𝒒𝑘 −1 , the last 𝑚 successive differences form ∆𝑮 = [𝒈 1 − 𝒈 0, . . . , 𝒈𝑘 −1 − 𝒈𝑘 −2 ] ∈ R𝑑𝑛 𝑣 ×𝑚 ,

(18)

𝑘 −1

(19)

∆𝑸 = [𝒒 − 𝒒 , . . . , 𝒒

𝑘 −2

−𝒒

𝑑𝑛 𝑣 ×𝑚

]∈R

with Tikhonov-regularised mixing coefficients 𝜸 ∗ = (∆𝑮𝑇 ∆𝑮 +𝜌 aa 𝑰 ) −1 ∆𝑮𝑇 𝒈𝑘 ,

𝜌 aa = max(

10−6 ∥∆𝑮 ∥ 2𝐹 , 10−12 ), 𝑚 (20)

and AA-extrapolated iterate 𝒒AA = (𝒒𝑘 −1 + 𝒈𝑘 ) − (∆𝑸 + ∆𝑮)𝜸 ∗ .

(21)

Two safeguards complete the scheme. A reliability guard discards the history and falls back to 𝒒𝑘 whenever ∥𝜸 ∗ ∥ 2 > 10, since large mixing coefficients indicate that AA is extrapolating beyond the trust region of the underlying linear model. A window bound adapts 𝑚 to material heterogeneity: under stiffness contrast 𝜇 max /𝜇 min ≥ 100×, a wide history extrapolates per-element drift into mesh-level oscillation (Fig. 3). Specifically, per-element 𝒑𝑒∗ trajectories have iterate-space curvature differing by orders of magnitude across regions, so a wide history (e.g., 𝑚 = 5) extrapolates a mean trend biased toward the soft-region high-curvature direction and excites cross-region oscillation. We therefore default to 𝑚 = 1 on heterogeneous meshes (“one-step momentum”), reverting to 𝑚 = 5 on homogeneous meshes where uniform curvature makes longer history safe. Dual-gate convergence criterion. A single step-norm test ∥𝒒𝑘+1 − 𝒒𝑘 ∥ ≤𝜀 fires falsely at 𝑘 = 0 (when the first PD step from 𝒒0 is small) and during contact-state transitions where AA satisfies the step norm at a degenerate non-fixed-point energy state. We therefore require both relative criteria simultaneously: ∥𝒒𝑘+1 − 𝒒𝑘 ∥ ≤1 𝜀 rel ∥𝒒𝑘 ∥ + 𝜀 abs

and

∥𝒃(𝒒𝑘+1 ) − 𝒃(𝒒𝑘 )∥ ≤ 1, 𝜀 rel ∥𝒃(𝒒𝑘 )∥ + 𝜀 abs

𝑘 ≥ 1. (22)

The PD-residual gate detects fixed-point convergence independently of iterate-step magnitude, eliminating both failure modes.

4.4

Hyperelasticity: Backward Trust-Region on the Proximal-Map Hessian

Forward Projected Newton (Sec. 3.3) must choose between eigenvalue clamping (Eq. (7), [Teran et al. 2005]) and absolute-value filtering (Eq. (8), [Chen et al. 2024b]), with no single filter uniformly optimal across the deformation regime. We lift this choice into the differentiable backward pass via a state-adaptive trust-region rule (Fig. 5), applied not to the raw element elastic Hessian 𝑯𝑒 prox of Eq. (9) but to the PD prox-map Hessian 𝑯𝑒 that drives the implicit-function-theorem (IFT) operator of the backward pass.

¯ from The backward pass applies the IFT operator 𝑯𝜓 (𝝈 ∗ ) + 𝑘𝑰 ∗ Eq. (14) when assembling 𝜕𝒑𝑒 /𝜕𝑭𝑒 , which drives 𝜕𝒃/𝜕𝒒 in Eq. (42) and ultimately the adjoint system Eq. (38). Corotated PD is automatically SPD (SO(𝑑) keeps the prox-map Jacobian well-conditioned). For Neo-Hookean PD, the per-element prox-map Hessian prox

𝑯𝑒

¯ ∈ R𝑑 ×𝑑 (𝝈 ∗ ) ≡ 𝑯𝜓 (𝝈 ∗ ) + 𝑘𝑰

(23)

becomes indefinite whenever 𝑯𝜓 has eigenvalues more negative ¯ than −𝑘—the canonical regime under high Poisson’s ratio or large volume change (Fig. 5a), and the visible cause of the wrist drift (Fig. 3). Vanilla inversion then yields a non-descent direction for the adjoint and numerically unreliable gradients. This motivates lifting trust-region eigenvalue filtering, originally proposed for forward Projected Newton [Chen et al. 2024a], into the backward pass. The projection is applied in the backward pass only (Fig. 5c): prox forward modifications of 𝑯𝑒 (PSD clamping, absolute filtering, blended filters with eigenvalue floors, per-DoF under-relaxation) all violate the fixed-point assumption of Type-II AA in Eqs. (18)–(21)— perturbing 𝒑𝑒∗ (𝒒) as a function of 𝒒 corrupts 𝜸 ∗ and destabilises the global iteration. Trust-region blend on the prox-map Hessian. For each element 𝑒, the family of SPD prox-map surrogates is prox prox prox 𝑯˜ 𝑒 (𝜏) = (1 − 𝜏) 𝑯𝑒 + 𝜏 |𝑯𝑒 |,

𝜏 ∈ [0, 1],

(24)

where |·| is per-element absolute-value projection in stretch space (eigendecompose, take |𝜅𝑖 |, recompose). Following Chen et al. [2024a], 1 the discrete sweep 𝜏 ∈ {0, , 1} recovers canonical filters (Fig. 5b): 2 𝜏 = 0 is the unprojected IFT, 𝜏 = 1 is absolute-value filtering, and 1 𝜏 = coincides with eigenvalue clamping ([Chen et al. 2024a], 2 Lemma 4.2). Closed-form differential of the filtered prox-map. With 𝜏 ∗ fixed (rule below), Eq. (14) is replaced by prox 𝑯˜ 𝑒 (𝜏 ∗ ) d𝝈 ∗ = 𝑘¯ d𝝈𝐹 ,

(25)

from which the filtered 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 assembles via the SVD differential as in Sec. 4.1. The downstream 𝜕𝒃/𝜕𝒒 in Eq. (42) inherits SPD structure element-wise, restoring well-posedness of Eq. (38) in regimes where DiffPD’s pure-PD adjoint degrades. State-adaptive selection on the last fixed-point increment. We drive 𝜌 with the last converged PD increment, not the whole-frame displacement: ∆𝒒∗ = 𝒒∗ − 𝒒∗−1 (26) is the difference between the last two PD iterates (restricted to free DoFs, after AA mixing). Whole-frame 𝒒𝑡 +ℎ − 𝒒𝑡 contains inertial drift unrelated to elastic stationarity, whereas the last PD increment is the pure Newton-direction signal the backward IFT actually inverts—making its quadratic-model fidelity the right adaptation criterion. With actual and PD-model decreases 1 ∆Φact = Φ(𝒒∗−1 ) − Φ(𝒒∗ ), ∆Φmod = |(∆𝒒∗ )𝑇 𝑨 ∆𝒒∗ |, (27) 2 trust-region ratio 𝜌 =

∆Φact , ∆Φmod

𝜌 ← 1 if ∆Φmod < 10−12,

(28)

Lai, et al.

Figure 5: Trust-region filter on the prox-map Hessian, applied in the backward pass only. (a) The IFT operator can be indefinite. Under high deformation or near inversion, the per-element prox-map Hessian picks up negative eigenvalues (Eq. (23)), making vanilla inversion through the backward IFT unstable. (b) Trust-region blend. A scalar 𝜏 ∈ [0, 1] interpolates between three canonical filters (Eq. (24)): unpro1 jected (𝜏 = 0), eigenvalue clamping (𝜏 = ), and absolute-value 2 filtering (𝜏 = 1). A state-adaptive rule selects 𝜏 ∗ from the trustregion ratio evaluated on the last PD increment (Eq. (29)). (c) Backward-only, 𝑨 reused. The forward fixed-point is left untouched—any forward modification would perturb the AA mixing and destabilise iteration. The filtered Hessian feeds the adjoint chain (Eq. (42)), and the trust-region ratio’s quadratic model reuses the persistent factor of 𝑨 at zero extra cost (Sec. 4.6).

Figure 6: One persistent factor, shared across forward, backward, and contact. (a) One persistent factor. The SPD operator 𝑨 is factorised once via METIS nested dissection into 𝑨−1 = 𝑺𝑇 𝑺 (Eq. (32)); the factor pair sits in persistent GPU buffers. (b) Two-SpMV inverse application. Any 𝑨−1 𝒗 is evaluated as two sparse multiplications 𝑺𝑇 (𝑺𝒗) (Eq. (33)); the dense inverse is never materialised. (c) Shared across all stages. The same factor pair serves the forward global solve, the contact Delassus operator (Sec. 4.7), and the backward adjoint (Sec. 4.8)—refactorisation is triggered only on topology, material, or damping change.

we extend with per-element stiffness coefficients for heterogeneity: ∑︁ 𝒇𝑑 = −(𝛼𝑴 + 𝑩 𝜷 )𝒗, 𝑩𝜷 = 𝛽𝑒 𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 , (30) 𝑒

and rule  1    ∗ 𝜏 = 2  1 

|𝜌 − 1| ≤ 𝜀 TR

(quadratic model accurate),

otherwise

(model inadequate; prefer absolute filter), (29) we take 𝜀 TR = 0.1 matching [Chen et al. 2024a]. Reusing 𝑨 as the quadratic model. ∆Φmod in Eq. (27) uses PD’s constant 𝑨, not a fresh per-step Newton Hessian. Since the PD surrogate is constant in 𝒒, 𝑨 provides a state-independent reference quadratic across the entire fixed-point trajectory at zero extra cost: the SpMV 𝑨∆𝒒∗ reuses the persistent factor pairs {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )} already populated by the forward pass (Sec. 4.6). TR filtering is invoked once per backward solve; the per-element 𝑑 ×𝑑 symmetric eigendecomposition in Eq. (24) is a closed-form Jacobi sweep on 3×3 matrices for tetrahedra, fully parallelised across 𝑛𝑒 .

4.5

Stiffness-Amplified Rayleigh Damping

High-frequency ringing in soft regions during force-peak transitions (observed as cross-boundary oscillation when stiff and soft regions couple at low Poisson’s ratio) is dissipated by a Rayleightype damping force—originally a linear combination of mass- and stiffness-proportional terms [Hughes 2012; Rayleigh 1894], which

with scalar mass coefficient 𝛼 and per-element stiffness coefficient 𝛽𝑒 . Implicit Euler absorbs 𝒇𝑑 into the linearly augmented PD operator: 𝛽𝑒 (1 + 𝛼ℎ)𝑴 ∑︁ (𝑤𝑒 + )𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 . 𝑨 = + (31) 2 ℎ ℎ 𝑒 Stiffness-amplified routing: 𝛽𝑒 ∝ 𝜇𝑒 . Standard heterogeneity-aware tuning suggests 𝛽𝑒 ∝ 1/𝜇𝑒 , so soft regions (where compliant highfrequency modes live) receive more dissipation. We adopt the opposite: 𝛽𝑒 = 𝛽 0 𝜇𝑒 /𝜇 ref with 𝜇 ref = max𝑒 𝜇𝑒 , so stiff regions are damped more aggressively. The rationale is the cross-boundary transfer pathway: in heterogeneous twist scenarios the high-frequency ringing contaminating the soft region is seeded by stiff-region ringing transmitted across the material interface; damping at the source preserves soft-region force-peak dynamics (its arms keep their elastic snap) while suppressing the seed. Empirically, 𝛽𝑒 ∝ 1/𝜇𝑒 overdamped soft-region deformation under force peak, while 𝛽𝑒 ∝ 𝜇𝑒 preserved amplitude and removed visible cross-boundary oscillation. Damping folds into 𝑨 at zero recurring cost. Both 𝛼𝑴 and 𝛽𝑒 𝑮𝑒𝑇 𝑮𝑒 are constant in 𝒒 and time-step independent. Eq. (31) is assembled once when {𝑤𝑒 } are computed from {𝜇𝑒 }, and the same nesteddissection 𝑺𝑺𝑇 factor pairs of Sec. 4.6 represent its inverse—no

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

refactorisation, no per-step damping solve, and the backward pass inherits the same factored representation automatically.

4.6

Unified GPU Pipeline: Persistent 𝑺𝑇 𝑺 across Forward, Backward, and Contact

FBA nested-dissection sparse inverse. The SPD 𝑨 from Eq. (31) is reordered by METIS nested dissection [Karypis and Kumar 1997] via permutation 𝑷 , then factorised [Zeng et al. 2025] (Fig. 6a) 𝑷 𝑨𝑷 𝑇 = 𝑳𝑫𝑳𝑇 ,

𝑺 = diag(𝑑𝑖−1/2 )𝑳 −1 𝑷,

𝑨−1 = 𝑺𝑇 𝑺.

(32)

The dense 𝑺𝑇 𝑺 is never materialised; it is applied via two SpMVs. For FEM meshes with METIS ordering, nnz(𝑺)/𝑛 2𝑣 ≈ 2–3%, and the 𝑗-th column of 𝑺 is supported on { 𝑗 }∪{ancestors of 𝑗 } in the elimination tree, computed column-independently in parallel. Two-SpMV evaluation. Once 𝑺 and 𝑺𝑇 are stored in CSR on the GPU, any 𝑨−1 𝒗 is two SpMVs (Fig. 6b), 𝑨−1 𝒗 = 𝑺𝑇 (𝑺𝒗).

(33)

Since 𝑨 decouples across 𝑑 coordinate axes, we maintain 𝑑 independent factor pairs {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )}𝑑𝑘=1 for per-axis SpMV without memory contention. These pairs sit in persistent file-scope GPU buffers populated once by the forward pass and accessed by the backward pass without re-factorisation (Fig. 6c). Why no algebraic preconditioner. Eq. (33) is an exact application of 𝑨−1 within floating-point precision; layering an algebraic preconditioner (Jacobi, Schwarz, or MAS with Nicolaides aggregation) adds per-level subdomain solves and interlevel transfers without improving accuracy, and requires outer PCG to recover fidelity. Empirically, MAS preconditioning was strictly slower than bare two-SpMV on our heterogeneous Neo-Hookean benchmarks at equal forward accuracy, because 𝑨 is constant across PD iterations and 𝑺𝑺𝑇 already represents 𝑨−1 exactly. We therefore reuse the same persistent factor for the forward global step and the backward adjoint solve. Batch Delassus computation. With contact Jacobian 𝑱 = [𝑱𝑛𝑇 , 𝑱𝑏𝑇 , 𝑱 𝑓𝑇 ]𝑇 ∈

computed via two GPU SpMMs: 𝒀 = 𝑺 𝑱 𝑇 , 𝑾 = 𝒀 𝑇 𝒀 . The columns 𝒂𝑐 = 𝑨−1 𝒋𝑐 extracted from 𝑺𝑇 𝒀 are cached for reuse in the position update (Sec. 4.7). Sparsity-preserving SpMM for basis-vector batches. The columns of 𝑱 𝑇 in the Delassus batch SpMM 𝒀 = 𝑺 𝑱 𝑇 are elementary basis vectors with a single non-zero at the contact node, so 𝒀 and the 𝒂𝑐 extracted from 𝑺𝑇 𝒀 inherit their sparsity from 𝑺 rather than from the RHS—each column carries on the order of nnz(𝑺)/𝑛 𝑣 nonzeros (2–3 % of 𝑛 𝑣 ) under METIS ordering. Collapsing this fill onto a single diagonal would degenerate 𝑾 = (𝑺 𝑱 𝑇 )𝑇 (𝑺 𝑱 𝑇 ) to near-diagonal and the downstream 𝑴 sys = Ω𝑾 Ω𝑇 +𝑬 would lose the off-diagonal coupling that propagates contact impulses, with 𝜅(𝑴 sys ) rising by orders of magnitude. We therefore require an SpMM preserving the sparsity of 𝑺 across all 𝐾 basis-vector RHSs rather than treating them as dense with output thresholding: in our implementation, a CSR-aware batched SpMM with an explicit pattern-analysis pass via cuSPARSE analyses the CSR pattern of 𝑺 (𝑘) once per refactorisation (Eq. (32)) and amortises it across all 𝐾 columns and 𝑑 axes. The same persistent {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )} then serves forward, backward, and contact stages with no extra data structure, matching the columnwise SpMV evaluation Eq. (33) to machine precision. Refactorisation policy. 𝑨 is refactorised only on (i) topology change (remesh), (ii) material change ({𝐸𝑒 } updated, e.g. inverse design), or (iii) damping coefficient change (Fig. 6c). Within a single forward rollout or backward adjoint solve at fixed {𝐸𝑒 , 𝛼, 𝛽 0 }, the factor is computed once and reused for the entire trajectory: this is the unified-factor amortisation DiffPD [Du et al. 2021] applies to forward, FBA [Zeng et al. 2025] extends to contact, and we extend further to the differentiable backward adjoint (Sec. 4.8).

4.7

Figure 7: Mesh-to-mesh deformable contact: a dominant failure mode. Penalty-based and pure-PD models leave residual error on non-convex two-body impacts—surface separation (top) or interpenetration (middle)—contaminating the backward gradients. We adopt the Signorini–Coulomb non-smooth NCP of [Zeng et al. 2025], integrated into the persistent factor pipeline (Sec. 4.6).

R𝐾 ×𝑑𝑛 𝑣 , 𝐾 = 𝐾𝑛 + 𝐾𝑏 + 𝐾 𝑓 , the Delassus operator 𝑾 = 𝑱 𝑨−1 𝑱 𝑇 factors as 𝑾 = (𝑺 𝑱 𝑇 )𝑇 (𝑺 𝑱 𝑇 ),

Signorini–Coulomb Complementarity Contact

We adopt the non-smooth NCP formulation of [Zeng et al. 2025] (Fig. 7), encoding unilateral contact, bilateral attachment, and Coulomb friction in a unified Fischer–Burmeister framework. For each contact 𝑐, 𝛿𝑛,𝑐 = 𝒏𝑇𝑐 𝒒 𝑣𝑐 − 𝑔𝑐 is the signed normal gap and 𝜹 𝑓 ,𝑐 = 𝑱 𝑓 ,𝑐 𝒗 the tangential slip. The KKT structure 𝛿𝑛,𝑐 ≥ 0, 𝜆𝑛,𝑐 ≥ 0, 𝛿𝑛,𝑐 𝜆𝑛,𝑐 = 0 and ∥𝝀 𝑓 ,𝑐 ∥ ≤ 𝜇 𝑗 𝜆𝑛,𝑐 is encoded [Fischer 1992; Zeng et al. 2025] as √︃ 2 + 𝑟 2 𝜆 2 = 0, 𝜑𝑛 (𝛿𝑛,𝑐 , 𝑟𝑛,𝑐 , 𝜆𝑛,𝑐 ) = 𝛿𝑛,𝑐 + 𝑟𝑛,𝑐 𝜆𝑛,𝑐 − 𝛿𝑛,𝑐 (34) 𝑛,𝑐 𝑛,𝑐 with 𝑟𝑛,𝑐 =ℎ 2𝑊𝑛,𝑐𝑐 , and analogously 𝜑 𝑓 (𝑠𝑐 , 𝑟 𝑓 ,𝑐 , 𝜉𝑐 ) for friction with 𝑠𝑐 = ∥𝜹 𝑓 ,𝑐 ∥ and 𝜉𝑐 = 𝜇 𝑗 𝜆𝑛,𝑐 − ∥𝝀 𝑓 ,𝑐 ∥. Differentiating Eq. (34) yields the per-contact NCP weights 𝛿𝑛,𝑐 𝑟𝑛,𝑐 𝜆𝑛,𝑐 𝜔𝑛,𝑐 = 1 − √︁ , 𝐸𝑛,𝑐 = 1 − √︁ 𝑟𝑛,𝑐 , (35) 2 2 2 2 2 2 𝛿𝑛,𝑐 + 𝑟𝑛,𝑐 𝜆𝑛,𝑐 𝛿𝑛,𝑐 + 𝑟𝑛,𝑐 𝜆𝑛,𝑐 and analogously 𝜔 𝑓 ,𝑐 , 𝐸 𝑓 ,𝑐 ; stacking gives diagonal matrices Ω, 𝑬. Reduced contact system. With block-weighted Jacobian 𝑱¯ = Ω𝑱 , the per-iteration linearised NCP [Zeng et al. 2025] reads 𝑴 sys ∆𝝀 = 𝒉 vec − 𝑱¯𝑨−1 𝒓 con,

(36)

𝑴 sys = Ω𝑾 Ω𝑇 + 𝑬, 𝒓 con = 𝒃 + 𝑱 𝑇 (Ω𝝀𝑘 ). The 𝐾 ×𝐾 system Eq. (36) is solved by dense LDLT since 𝐾 ≪𝑛 𝑣 . Multipliers update via block

Lai, et al.

projection Π C = diag(Π ≥0, 𝑰 , Π C𝜇 ); the position is recovered without an additional global solve via the cached contact columns: 𝒒𝑘+1 = 𝑺𝑇 (𝑺𝒃) +

𝐾 ∑︁

𝜔𝑐 𝜆𝑐𝑘+1 𝒂𝑐 .

(37)

𝑐=1

4.8

Differentiable Backward Pass

Given a scalar loss L(𝒒𝑡 +ℎ , 𝒗𝑡 +ℎ ) and incoming gradients 𝒒¯𝑡 +ℎ , 𝒗¯𝑡 +ℎ , we propagate sensitivities via the adjoint method on the converged (𝒒∗, 𝝀 ∗ ) with the active contact set fixed. The residuals 𝑹 𝒒 = 𝒒∗ − 𝑨−1 (𝒃(𝒒∗ ) + 𝑱 𝑇 (Ω𝝀 ∗ )) and 𝑹𝝀 = 𝑴 sys 𝝀 ∗ − (𝒉 vec − 𝑱¯𝑨−1 𝒓 con ) vanish identically; differentiating yields the block KKT adjoint system (𝜕𝑹 𝒒 /𝜕𝒒)𝑇 (𝜕𝑹𝝀 /𝜕𝒒)𝑇 𝒚𝒒 𝒒¯ = , (38) 𝒚𝝀 0 (𝜕𝑹 𝒒 /𝜕𝝀)𝑇 𝑴 sys with 𝒒¯ = 𝒒¯𝑡 +ℎ + 𝒗¯𝑡 +ℎ /ℎ and 𝜕𝑹 𝒒 𝜕𝒃 [𝛿𝒒] = 𝛿𝒒 − 𝑨−1( [𝛿𝒒]), 𝜕𝒒 𝜕𝒒 𝜕𝑹 𝒒 [𝛿𝝀] = −𝑨−1 𝑱 𝑇 (Ω 𝛿𝝀), 𝜕𝝀 𝜕𝑹𝝀 𝜕𝒃 [𝛿𝒒] = 𝑱¯ 𝛿𝒒 − 𝑱¯𝑨−1( [𝛿𝒒]). 𝜕𝒒 𝜕𝒒 The local-step Jacobian through per-element projection is ∑︁ 𝜕𝒑𝑒∗ 𝜕𝒃 [𝛿𝒒] = 𝑤𝑒 𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 𝛿𝒒, 𝜕𝒒 𝜕𝑭𝑒 𝑒

(39) (40) (41)

(42)

where 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 for Neo-Hookean is assembled from Eq. (25) with prox prox 𝑯𝑒 replaced by 𝑯˜ 𝑒 (𝜏 ∗ ) (Sec. 4.4); for corotated PD it is the polar-decomposition differential [Du et al. 2021], intrinsically SPD. The 𝑑𝑛 𝑣×𝑑𝑛 𝑣 𝒒-block is solved by direct two-SpMV via the persistent factor pairs (Eq. (33)); the 𝐾×𝐾 𝝀-block by dense LDLT. When 𝐾 = 0, the entire backward solve reduces to a single GPU two-SpMV call. Backbone adjoint and gradient routing. The backbone adjoint 𝝁 = 𝑨−1𝒚𝒒 = 𝑺𝑇 (𝑺𝒚𝒒 ) is evaluated via Eq. (33). Gradients flow as 𝜕𝒇state 𝑇 𝜕L 1 𝑴 ) 𝑴𝝁 − 𝒗¯𝑡 +ℎ , (43) += 2 𝝁 + ( 𝜕𝒒𝑡 ℎ 𝜕𝒒𝑡 ℎ 𝜕L 𝜕L ¯ += ℎ 𝒒, += ℎ 2 𝑴𝝁, (44) 𝜕𝒗𝑡 𝜕𝒇ext 𝜕L += 𝑉𝑒 (𝝁𝑇 𝑮𝑒𝑇 𝒑𝑒∗ − 𝝁𝑇 𝑮𝑒𝑇 𝑮𝑒 𝒒∗ ), (45) 𝜕𝑤𝑒 and per-element Young’s modulus 𝜕L/𝜕𝐸𝑒 via the chain rule on Eqs. (1), (12).

4.9

Algorithm Summary

Algorithms 1 and 2 together integrate one forward step (Sec. 4.3– 4.7) and the corresponding backward step (Sec. 4.8). The factor pairs {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )} are populated once per material/topology change and persist across both passes through the forward cache F . Computational complexity. Per forward iteration: local step O(𝑛𝑒 ); global step O(nnz(𝑺)) via two GPU SpMVs (nnz(𝑺)/𝑛 2𝑣 ≈ 2–3%); contact O(𝐾 2 ) dense LDLT; AA mixing O(𝑚 2𝑑𝑛 𝑣 ) (𝑚 ≤ 5). Delassus pre-computation costs O(𝐾 · nnz(𝑺)) via GPU SpMM, amortised over all iterations of a timestep. Backward: one GPU two-SpMV reusing forward buffers for the 𝑑𝑛 𝑣 𝒒-block, plus per-element 𝑑 ×𝑑

Algorithm 1 DiffPhD — Forward Step Require: 𝒒𝑡 , 𝒗𝑡 , 𝒇ext, ℎ, {(𝐸𝑒 , 𝑉𝑒 , 𝜈)}, 𝛼, 𝛽 0, 𝑚, 𝜀 rel, 𝜀 abs Ensure: 𝒒𝑡 +ℎ , 𝒗𝑡 +ℎ , 𝝀 ∗ ; forward cache F = {𝒒∗−1, 𝒒∗, 𝝀 ∗, {𝒑𝑒∗, 𝝈𝑒∗, 𝑼𝑒 , 𝑽𝑒 }, 𝑾, {𝒂𝑐 }, A} 1: if 𝑨 stale then 2: {material/topology/damping change} ¯ 𝑒 , 𝜆𝑒 ) ∝ 𝜇𝑒 ; 𝛽𝑒 ← 𝛽 0 𝜇𝑒 /𝜇 ref 3: 𝑤𝑒 ← 𝑘(𝜇 (12),(30) P 4: 𝑨← (1+𝛼ℎ)𝑴/ℎ 2 + 𝑒 (𝑤𝑒 + 𝛽𝑒 /ℎ)𝑉𝑒 𝑮𝑒𝑇 𝑮𝑒 (31) Factorise 𝑷 𝑨𝑷 𝑇 = 𝑳𝑫𝑳𝑇 ; 𝑺 ← diag(𝑑𝑖−1/2 )𝑳 −1 𝑷 ; upload 5: {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )} to GPU (32) 6: end if 7: 𝒒˜ ← 𝒒𝑡 + ℎ𝒗𝑡 + ℎ 2 𝑴 −1 (𝒇ext + 𝒇state ) (3) 8: if 𝐾 > 0 then 9: 𝒀 ← 𝑺 𝑱 𝑇 ; 𝑾 ← 𝒀 𝑇 𝒀 ; 𝒂𝑐 ← (𝑺𝑇 𝒀 )·𝑐 ; 𝑟𝑐 ←ℎ 2𝑊𝑐𝑐 10: end if 11: 𝒒 0 ← 𝒒𝑡 ; 𝝀 0 ← 0; H ← ∅ 12: for 𝑘 = 0, . . . , 𝐾max −1 do Local step (parallel over 𝑒): solve Eq. (14) with means 13: ¯ 𝑘¯ for 𝝈 ∗ ; 𝒑 ∗ ← 𝑼𝑒 diag(𝝈 ∗ )𝑽 𝑇 ¯ 𝜆, 𝜇, (13) 𝑒 𝑒 𝑒 𝑒 P 14: 𝒃 ← (𝑴/ℎ 2 )𝒒˜ + 𝑒 𝑤𝑒 𝑉𝑒 𝑮𝑒𝑇 𝒑𝑒∗ (11) 15: if 𝐾 > 0 then Ω, 𝑬 ← NCP-weights(𝜹, 𝒓, 𝝀𝑘 ); 𝒓 con ←𝒃 +𝑱 𝑇 Ω𝝀𝑘 ; 𝑴 sys ← 16: Ω𝑾 Ω𝑇 +𝑬 (36) −1 (𝒉 𝑇 (𝑺𝒓 𝑘+1 ← Π (𝝀𝑘 + ∆𝝀) ¯ 17: ∆𝝀 ← 𝑴 sys − 𝑱 𝑺 )); 𝝀 vec con C 18: end if P 19: 𝒒ˆ ← 𝑺𝑇 (𝑺𝒃) + 𝑐 𝜔𝑐 𝜆𝑐𝑘+1 𝒂𝑐 (37) 20: {Bounded-window AA (window 𝑚, default 𝑚 = 1 on het)} 21: 𝒈𝑘 ← 𝒒ˆ − 𝒒𝑘 ; update H (keep last 𝑚 pairs) 22: if |H |≥ 2 then 23: 𝜸 ∗ ← (∆𝑮𝑇 ∆𝑮 + 𝜌 aa 𝑰 )−1 ∆𝑮𝑇 𝒈𝑘 (20) ˆ H ← ∅ (21) 24: 𝒒𝑘+1 ← 𝒒ˆ − (∆𝑸 + ∆𝑮)𝜸 ∗ if ∥𝜸 ∗ ∥ ≤ 10, else 𝒒, 25: else 26: 𝒒𝑘+1 ← 𝒒ˆ 27: end if 28: Enforce Dirichlet BCs 29: if 𝑘 ≥ 1 and dual gate (22) holds then 30: break 31: end if 32: end for 33: 𝒒𝑡 +ℎ ← 𝒒𝑘+1 ; 𝒗𝑡 +ℎ ← (𝒒𝑡 +ℎ − 𝒒𝑡 )/ℎ; 𝝀 ∗ ←𝝀𝑘+1 34: Cache F ← {𝒒 ∗−1 , 𝒒 ∗ , 𝝀 ∗ , {𝒑𝑒∗ , 𝝈𝑒∗ , 𝑼𝑒 , 𝑽𝑒 }, 𝑾, {𝒂𝑐 }, A} for backward

symmetric eigendecomposition for the trust-region projection. The 𝐾 = 0 fast path bypasses all contact adjoint computation. Damping coefficients 𝛼, 𝛽 0 enter Eq. (31) at assembly only and incur no additional per-iteration cost in either pass.

RESULTS & FINDINGS

We evaluate DiffPhD along the three axes that motivate the framework: (i) forward stability and speed on heterogeneous, contact-rich elastodynamics where prior PD solvers diverge; (ii) gradient accuracy and end-to-end speedup on differentiable inverse problems –system identification, initial-state optimization, and trajectory optimization; and (iii) transfer to robotics via heterogeneous contact-rich

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Algorithm 2 DiffPhD — Backward Step (uses cache F and persistent {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )} from Alg. 1) Require: Forward cache F ; persistent GPU factor pairs {(𝑺 (𝑘), (𝑺𝑇 )(𝑘) )}; incoming gradients 𝒒¯𝑡 +ℎ , 𝒗¯𝑡 +ℎ ; tolerance 𝜀 TR Ensure: 𝜕L/𝜕𝒒𝑡 , 𝜕L/𝜕𝒗𝑡 , 𝜕L/𝜕𝒇ext , 𝜕L/𝜕𝑤𝑒 , 𝜕L/𝜕𝐸𝑒 1: 𝜕L/𝜕𝒒𝑡 += −¯ 𝒗𝑡 +ℎ /ℎ; 𝒒¯ ← 𝒒¯𝑡 +ℎ + 𝒗¯𝑡 +ℎ /ℎ 2: if Neo-Hookean then 3: ∆𝒒∗ ← 𝒒∗ − 𝒒∗−1 ; ∆Φact ← Φ(𝒒∗−1 ) − Φ(𝒒∗ ) (26) 1 ∗ 𝑇 ∗ (27) 4: ∆Φmod ← |(∆𝒒 ) 𝑨∆𝒒 | 2 5: 𝜌 ← ∆Φact /∆Φmod (28) 1 6: 𝜏 ∗ ← if |𝜌 − 1|≤ 𝜀 TR else 1 (29) 2 7: end if 8: for all element 𝑒 in parallel do 9: if corotated then 10: 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 ←polar-decomposition differential 11: else prox ¯ 12: 𝑯𝑒 ← 𝑯𝜓 (𝝈𝑒∗ ) + 𝑘𝑰 (23) prox prox prox ∗ ˜ 13: 𝑯𝑒 ← (1 − 𝜏 )𝑯𝑒 + 𝜏 ∗ |𝑯𝑒 | (24) prox 14: 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 ←SVD diff. from 𝑯˜ 𝑒 and (𝑼𝑒 , 𝑽𝑒 ) (25) 15: end if 16: 𝑲𝑒 ←𝑤𝑒 𝑉𝑒 𝑮𝑒𝑇 (𝜕𝒑𝑒∗ /𝜕𝑭𝑒 )𝑮𝑒 (42) 17: end for P 18: Assemble 𝜕𝒃/𝜕𝒒 ← 𝑒 𝑲𝑒 via FEM scatter 19: if |A|= 0 then 20: 𝒚𝒒 ←apply (𝜕𝑹 𝒒 /𝜕𝒒)−𝑇 𝒒¯ via two-SpMV (39),(33) 21: else 22: Restore 𝑾, Ω, 𝑬, 𝑴 sys, 𝑱¯ from F 23: Solve Eq. (38) (𝒒-block: two-SpMV; 𝝀-block: dense LDLT) 24: end if 25: 𝝁 ← 𝑺𝑇 (𝑺𝒚 𝒒 ) (33) 26: Route gradients via Eqs. (43)–(45) and chain rule on Eqs. (1), (12)

and the adjoint pass. Eigenvalue filtering on the per-element proxmap Hessian is fused into the same element-parallel kernel that assembles 𝜕𝒑𝑒∗ /𝜕𝑭𝑒 . Hardware. All experiments were conducted on a desktop system equipped with a single NVIDIA RTX 4090 GPU (24 GB VRAM), and an Intel(R) Core(TM) i7-12700 CPU (12 cores @ 4.9 GHz) with 16 GB of RAM. Wall-clock times are reported end-to-end and include host– device transfer. Baselines. We compare against representative prior work spanning the key method families our contribution touches. NewtonCholesky provides a sparse direct implicit-Newton accuracy reference, complemented by its Neo-Hookean variant with stable energy [Smith et al. 2018b]. DiffPD [Du et al. 2021] serves as our most direct differentiable PD baseline. To evaluate performance under different preconditioning strategies, we also include a variant of DiffPD equipped with MAS [Wu et al. 2022], which replaces the standard preconditioner with a domain decomposition based approach. While our proposed DiffPhD adopts FBA [Zeng et al. 2025] as its underlying forward solver to leverage the state-of-the-art GPU PD solver with NCP frictional contact, we evaluate its efficiency against the aforementioned differentiable baselines. For the differentiable baselines we use the authors’ reference implementations and a uniform L-BFGS optimizer with identical convergence thresholds; following [Du et al. 2021], all wall-clock comparisons use a forward tolerance of 10−4 and a backward tolerance of 10−6 unless stated otherwise. Notation. homo denotes a single-stiffness model and hetero a model with at least one ≥ 10× stiffness contrast across sub-domains. n/c denotes failure to converge within 104 iterations or 24 wallclock hours; n/d denotes a non-differentiable solver (used only as a forward reference).

5.2 manipulation and a Real2Sim study. We close with a unified ablation that disentangles the contribution of each architectural component. Table 2 summarizes every benchmark; rows in bold mark heterogeneous variants where prior differentiable PD solvers either diverge or produce numerically unreliable gradients.

5.1

Experimental Setup

Implementation. DiffPhD is implemented in C++/CUDA on top of the DiffPD [Du et al. 2021] codebase. Sparse matrix–vector and matrix–matrix products run on the GPU through cuSPARSE: the per-axis application of 𝑨−1 is two SpMVs against the persistent (𝑺, 𝑺𝑇 ) factor pair (Sec. 4.5), and the batched Delassus assembly 𝑾 = 𝑱 𝑨−1 𝑱 𝑇 uses cusparseSpMM with the CSR_ALG2 algorithm and an explicit preprocess pass for sparsity-pattern analysis (the default algorithm selector returns rank-deficient output on the unit-spike RHS used here; Sec. 4.5). The dense 𝑴 sys contact block is solved by Eigen::LDLT on the CPU. Sparse factorizations of 𝑨 = 𝑺𝑇 𝑺 are computed once per topology, material, or damping-coefficient change with METIS nested dissection [Karypis and Kumar 1998] and uploaded to persistent GPU buffers reused by both the forward

Heterogeneous Forward Simulation

We first verify that DiffPhD resolves heterogeneous PD dynamics that destabilize prior solvers. The three scenes (Cantilever, Armadillo, Crab) all carry stiffness contrasts of 10×–100×. PD energy under heterogeneity. The Cantilever test partitions a beam into three segments, with the middle third set as a stiff region (red, 10×) between two soft ends (blue, 1×). The beam is initially twisted and then released to return to its equilibrium state. Under this motion, the central red section restores to its rest profile rapidly with minimal deformation, while the blue ends exhibit significant lagging and larger oscillations, matching physical intuition. Under axial loading, elongation is similarly concentrated in the blue regions while the red section resists extension (Fig. 2(a)). To capture this, DiffPhD employs a stiffness-aware weight 𝑤𝑒 ∝ 𝜇𝑒 absorbed directly into 𝑨 at assembly (Sec. 4.1). This ensures the energy landscape faithfully represents the material distribution and maintains a smooth, artifact-free interface. Neo-Hookean heterogeneity on Armadillo. We evaluate the performance of DiffPhD using a twisted Armadillo scene with varying stiffness ratios. The mesh is partitioned by height into three equal segments to introduce material heterogeneity. Across our tests, we vary the stiffness contrast between these regions ranging from 10×

Lai, et al.

Table 2: Benchmark scenes used throughout the evaluation. #DoF is the number of degrees of freedom (3× vertices in 3D) and #E is the number of finite elements; ℎ is the timestep. Bold rows indicate heterogeneous variants where prior differentiable PD solvers fail to converge, and which therefore form the core of our evaluation. Section

Scene

#DoF

ℎ (ms)

Het.

Contact

Hyper.

Cantilever Armadillo (twist; homo., het.) Crab (homo.)

2,268 54,855 2,235,655

500 44,337 172,587

10 30 30

✓ ✓ –

– ✓ ✓

Crab (het.)

2,235,655

172,587

✓

Heterogeneous PD energy Stiff/soft body partitioning Shell–joint single-material reference Shell–joint composite

Gatorman Napkin 25×25 (homo.) Napkin 25×25 (het.) Napkin 50×50 (homo.) Napkin 50×50 (het.)

20,814 4,056 4,056 15,606 15,606

24,605 625 625 2,500 2,500

2 2 2 2 2

✓ – ✓ – ✓

✓ ✓ ✓ ✓ ✓

– – – – –

Complex mesh-to-mesh contact Codimensional cloth on sphere Codimensional cloth, two-stiffness High-resolution cloth on sphere High-resolution two-stiffness cloth

Bouncing Ball (homo.) Bouncing Ball (het.)

9,132 9,132

1,288 1,288

4 4

– ✓

✓ ✓

Plant

29,763

3,863

–

System ID, single material System ID, three sub-region materials System ID, oscillatory dynamics

§5.4.2

Bunny (het.) Routing Tendon

7,062 2,475

1,601 512

1 10

✓ ✓

✓ –

– –

§5.4.3

Torus

3,204

568

–

✓

–

Locomotion (rolling)

Oreo

153,666

176,932

0.25

✓

–

Dice

3,123

4,588

✓

–

Heterogeneous cookie-stack contact Real2Sim die manipulation

§5.2

§5.3

§5.4.1

§5.5

to 100×. Through our experiments, we empirically observe that the baseline DiffPD is generally limited to a stiffness contrast of approximately 50×, whereas DiffPhD continues to provide stable results across the entire tested range up to 100× (Fig. 8). Crab: shell–joint composite. The Crab follows the canonical scenario from [Trusty et al. 2023]: featuring a stiff carapace (10×) coupled to soft joints. We evaluate the forward computation times of DiffPhD against three baselines: Newton-Cholesky, DiffPD, and MAS, with detailed performance metrics summarized in Table 3. In our experiment, the crab is pulled upward to observe its deformation and contact response (Fig. 9). While our method maintains efficiency across both forward and backward passes, we further compare DiffPhD with a specialized heterogeneous solver, Subspace MFEM [Trusty et al. 2023]. We observe that Subspace MFEM fails to handle scenarios involving contact, resulting in either a complete simulation stall or significant mesh artifacts (Fig. 10), whereas DiffPhD robustly handles these interactions while maintaining material contrast.

5.3

Task / role in the evaluation

Contact-Rich Forward Simulation

We next isolate the unified contact pipeline. Throughout this work, contact interactions are resolved using a complementarity-based contact model within our unified NCP formulation. To evaluate its efficacy, we present two benchmarks of escalating difficulty: Gatorman (a complex character mesh from [Trusty et al. 2023] striking a sphere, focusing on the stability of transient contact and

Initial-state optimization Muscle-energy backward (energyrouting)

Napkin (Fig. 12, featuring both homogeneous and heterogeneous variants). Gatorman. This benchmark evaluates our unified contact pipeline’s capability in handling complex mesh-to-mesh interactions with significant material heterogeneity. In this scenario, a high-resolution Gatorman character mesh wields a sword—with a Young’s modulus of 6 × 106 Pa (60× the reference stiffness of the Gatorman’s body) to strike a target sphere set to 1×105 Pa, matching the body’s reference stiffness (Fig. 11). Our model inherently supports the Signorini– Coulomb complementarity condition, allowing the solver to resolve both static and dynamic friction within the unified NCP framework as the sword interacts with the sphere’s surface. Despite the transient nature of the impact and the rapid re-establishment of contact sets, our stiffness-aware projective assembly maintains tight convergence without the non-physical oscillations or penetration often associated with penalty-based methods. Napkin: heterogeneous codimensional contact. The heterogeneous Napkin benchmark evaluates contact stability under significant material contrast. As shown in Fig. 12, a soft region (0.1× stiffness) drapes around an obstacle while a stiff region (1× stiffness) resists deformation, creating asymmetric contact patches that grow from 6% to 50% of the mesh—mapping to the contact ratios in Table 4. While DiffPD exhibits severe computational overhead as contact density increases—with backward costs escalating to 4198.28,ms at 50 × 50 resolution—DiffPhD maintains high efficiency via our unified GPU pipeline. At the highest contact density, DiffPhD

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Frame 10

100x stiff 5.00e6 50x stiff 3.536e6 20x stiff 2.236e6 1x & all mid 5.00e5

Frame 30

20x soft 1.118e5 50x soft 7.071e4

10x

20x

50x

100x

100x soft 5.00e4 Young’s modulus E (log scale)

Figure 8: Armadillo under varying stiffness contrast. The mesh is partitioned into three height-stacked sub-regions whose relative Young’s moduli sweep from 10× to 100× (columns); within each column, the upper torso is stiffened relative to the legs by the indicated ratio, with intermediate regions interpolated on a log scale (right colour bar, ranging from 5.00×104 Pa for the softest sub-region to 5.00×106 Pa for the stiffest). Two snapshots are shown per contrast: frame 10 (top row) captures the onset of twist where the stiff torso resists rotation while the soft legs deform; frame 30 (bottom row) the post-impact pose after the release, where the residual deformation concentrates in the soft regions. DiffPhD remains stable across the entire range with the torso–leg interface preserved and no element inversion, while DiffPD fails to converge beyond ∼ 10× as its constant PD matrix loses spectral conditioning under high stiffness contrast (Sec. 4.2). achieves up to a 23.53× speed-up in the backward pass. As visualized in Fig. 12 (bottom row), our solver captures these distinct physical behaviors with zero penetration and stable transitions. These results demonstrate that DiffPhD efficiently resolves largescale, high-density contact problems where traditional differentiable solvers become computationally prohibitive.

5.4

Differentiable Inverse Problems

We now demonstrate that the speedups of Sec. 5.2–5.3 translate into faster, more reliable optimization on real downstream tasks. We organize the experiments along the three application axes of DiffPD [Du et al. 2021]—system identification, initial-state optimization, trajectory optimization—and we add heterogeneous variants to the task classes to expose where prior solvers break. 5.4.1

System Identification.

Homogeneous Bouncing Ball. We replicate the DiffPD [Du et al. 2021] Bouncing Ball benchmark: recover the Young’s modulus and Poisson’s ratio of a falling soft sphere from a single 100-frame trajectory of floor-impact dynamics (Fig. 13, left three columns). In this scenario, DiffPhD matches the recovery accuracy of prior solvers while maintaining high computational efficiency. On our hardware, the solver achieves a per-evaluation cost of 14.65s for the forward pass and 70.62s for the backward pass. These results

confirm that our unified framework is both accurate and performant in standard, uniform-material scenarios, reaching a sharper loss minimum compared to baselines within ∼30 L-BFGS evaluations. Heterogeneous Bouncing Ball. The more demanding variant partitions the ball into three coloured sectors (red, purple, blue in Fig. 13, right three columns) with distinct (𝐸𝑖 , 𝜈𝑖 ) per sector at a ground-truth contrast of {3×, 0.1×, 1×} in Young’s modulus, and the optimiser must jointly recover all six material parameters from a single rolling/bouncing trajectory. The asymmetric mass distribution breaks the rotational symmetry of the homogeneous variant, so the per-sector deformation under impact becomes the disambiguating signal. While the material heterogeneity typically causes baseline gradients to desynchronize across sectors—leading to a premature plateau in the loss—DiffPhD successfully recovers all six material parameters. This is enabled by our adaptive trust-region filter (Sec. 4.4), which fires on contact-transition frames to ensure the adjoint gradients remain well-conditioned. The recovered trajectory tracks the ground truth pose-by-pose, demonstrating that complex material heterogeneity is now a tractable target for differentiable system identification. Plant: system identification on an articulated geometry. We evaluate system identification on a potted plant mesh whose slender

Lai, et al.

Heterogeneous

Homogeneous

Table 3: Crab benchmark: per-timestep wall-clock (s) and convergence quality. Each cell reports homogeneous / heterogeneous values; the heterogeneous Crab has a 10×-stiff red carapace coupled to 1×-stiff blue joints. The bottom block reports DiffPhD’s forward-time speed-up over each baseline (homo / hetero). DiffPD converges to identical iterates on both variants because it does not route per-element stiffness into 𝑨 (Sec. 4.2). Forward tolerance 10−1 . Lower is better; best per row in bold.

frame 1

frame 15

Figure 9: Crab pull-up: homogeneous vs. heterogeneous response. The crab is pulled upward at the top of its carapace (arrow). Top row, heterogeneous: a stiff carapace (10×, red) lifts the soft legs (1×, blue) cleanly off the ground between frame 1 and frame 15. Bottom row, homogeneous: the single-stiffness baseline cannot transmit the lift force without material contrast and produces excessive global deformation.

branches and leaves induce strong geometric non-uniformity in the dynamic response: branches at different heights deform with distinct curvature signatures under gravity-driven oscillation, giving a 100-frame trajectory rich enough to disambiguate the material parameters from a single rollout. The optimiser recovers the plant’s Young’s modulus and Poisson’s ratio from this trajectory (Fig. 14). All three methods converge close to the ground-truth modulus (∼ 106 Pa) and Poisson’s ratio (∼ 0.40) on the parameter axes (Fig. 15, left and centre panels), but they differ markedly in residual loss: DiffPhD reaches 0.029 in 23 L-BFGS evaluations at 32.56 s forward and 28.74 s backward per evaluation (Tab. 5), whereas DiffPD and MAS plateau at 0.692 after 15 evaluations—a 23.9× larger residual. The plateau (visible as the horizontal tail of the blue and red loss curves in Fig. 15, right panel) indicates L-BFGS line search failure: with the prox-map Hessian developing indefiniteness in slender branch elements at peak swing curvature, the baselines’ adjoint gradients lose descent direction, while DiffPhD’s trust-region filter (Sec. 4.4) keeps the gradient on a descent direction throughout the trajectory. 5.4.2

Initial-State Optimisation.

Bunny: inverse design. On the heterogeneous Bunny benchmark, we optimise the bunny’s initial position, orientation (Euler angles), and velocity so that its centre of mass reaches a target after 100

Method

Metric

Homo

Hetero

Newton-Cholesky

forward backward loss grad

638.07 187.78 22.20 2.60

832.37 186.49 19.47 5.86

DiffPD [Du et al. 2021]

forward backward loss grad

366.92 674.20 19.15 1.01

MAS [Wu et al. 2022]

forward backward loss grad

250.76 99.22 22.05 2063.62

547.04 342.46 19.15 274.72

DiffPhD (Ours)

forward backward loss grad

22.14 48.78 −2.04 612.12

42.22 336.63 18.88 0.00

Speed-up (forward)

vs. Newton-Cholesky vs. DiffPD [Du et al. 2021] vs. MAS [Wu et al. 2022]

28.83× 16.58× 11.33×

19.71× 8.69× 12.96×

frames of free-fall and bouncing contact (Fig. 16). The trajectory traverses several contact make/break transitions on the floor; at each transition, the active contact set changes discontinuously and the per-step Delassus operator 𝑾 = 𝑱 𝑨−1 𝑱 𝑇 rebuilds with a different contact Jacobian 𝑱 . Solvers whose backward adjoint propagates gradients through these transitions without contact-aware factor reuse (Sec. 4.6) accumulate residual error each time the active set flips, and L-BFGS line search loses descent direction. Newton-Cholesky fails to converge on this scene under the heterogeneous variant; DiffPD and MAS converge to suboptimal residuals (0.403 and 0.374 respectively), with MAS stalling because its multilevel preconditioner must refactorise on every active-set change while the heterogeneity in the bunny disrupts its subdomain conditioning (Fig. 17, red curve in loss panel). DiffPhD converges to a final loss of 0.132 in 30 L-BFGS evaluations at 6.16 s forward and 4.55 s backward per evaluation, a 1.74× per-evaluation speed-up over DiffPD and 2.78× over MAS (Tab. 5). The recovered final loss is 3.05× lower than DiffPD’s and 2.83× lower than MAS’s, showing that the speed-up is not purchased through degraded convergence. The convergence curves (Fig. 17) confirm that all three methods track similar parameter trajectories in the first few L-BFGS evaluations, but DiffPhD’s backward gradient remains well-conditioned through the bouncing transitions while MAS’s loss plateaus. The persistent 𝑺𝑇 𝑺 factor (Sec. 4.6) is reused verbatim across forward, contact-Delassus, and backward

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

FEM

Subspace MFEM

DiffPHD (Ours)

Heterogeneous

Homogeneous

Initial state

Figure 10: Crab pull-up: differentiable solver comparison. From left to right: initial state, FEM (Newton-NH), Subspace MFEM [Trusty et al. 2023], and DiffPhD. Top row, homogeneous (blue): a single stiffness across the entire body. Bottom row, heterogeneous (red shell, 10× stiffer than the blue legs): the shell–joint composite of [Trusty et al. 2023]. Under upward pulling, Subspace MFEM exhibits unbounded local stretching at the claw and shell–joint boundary (“Large Deformation”), with the homogeneous case collapsing more severely than the heterogeneous one. DiffPhD preserves shape in both regimes and is the only differentiable solver among the four to remain stable throughout the contact-rich trajectory. Table 4: Napkin benchmark: per-step wall-clock (s) and loss across contact density and mesh resolution. Each cell reports homogeneous / heterogeneous values; the heterogeneous variant has light regions 1× stiff and dark regions 0.1×. Three metrics are reported per scene: forward and backward time (s), and final loss. Contact set size denotes the fraction of surface vertices in active contact. Entries marked n/a fall into two cases: Cholesky failed to converge on the entire 50×50 mesh, while DiffPD at 50×50 with 50% contact set size was skipped due to prohibitive estimated runtime (projected ∼10–12h for a single forward-backward pass). Best per column (within a resolution block) in bold. Forward tolerance 10−4 throughout. Resolution: 25×25

Resolution: 50×50

Contact set size

Method

Metric

24%

38%

50%

24%

38%

50%

Cholesky

forward backward loss

3.78/5.30 31.77/57.68 25.46/8.16

21.70/163.01 132.58/190.23 13.19/17.10

52.08/657.49 324.97/592.59 17.87/15.58

65.28/1081.04 626.39/1205.74 20.08/22.03

n/a / n/a n/a / n/a n/a / n/a

DiffPD

forward backward loss

3.76/5.27 6.09/9.31 25.46/8.16

21.80/163.35 23.33/201.25 13.19/17.10

51.95/651.98 45.74/794.18 17.87/15.58

65.40/1085.17 68.61/1724.21 20.08/22.03

65.81/98.45 89.75/352.38 −49.97/22.28

1061.19/2586.27 1018.70/2441.01 −31.52/0.93

4598.97/9417.71 2743.15/8543.97 −0.97/12.51

4528.02/ n/a 4198.28/ n/a −29.17/ n/a

MAS

forward backward loss

3.65/3.28 30.94/55.19 25.46/4.53

21.79/21.77 30.71/56.11 13.19/16.00

52.02/47.69 30.59/55.22 17.87/16.91

65.31/46.23 28.93/53.84 20.08/22.37

65.75/47.47 114.28/211.49 −49.97/18.81

1056.17/290.04 120.02/219.07 −31.52/0.83

4636.76/445.86 122.14/214.86 −0.97/−13.15

4512.30/428.38 122.86/213.81 −29.17/−30.46

DiffPhD

forward backward loss

3.22/3.16 4.53/5.50 33.66/4.53

8.63/21.50 5.00/24.80 19.59/16.00

14.40/47.73 6.23/52.28 22.19/16.91

17.81/46.39 7.04/73.00 21.83/22.37

40.59/46.71 35.11/42.32 −32.99/−57.95

179.88/184.75 56.18/73.66 −25.78/3.84

356.48/381.22 122.35/157.95 −30.53/29.29

624.86/640.73 167.01/296.78 −30.80/37.04

Speed-up

vs. DiffPD vs. MAS

1.27/1.69× 4.46/6.76×

3.31/7.87× 3.85/1.68×

4.74/14.46× 4.01/1.03×

5.39/23.53× 3.79/0.84×

2.06/5.06× 2.38/2.91×

8.81/19.46× 4.98/19.70×

15.33/33.31× 9.94/1.23×

11.02/ n/a 5.85/0.69×

adjoint stages; when the active contact set changes between successive L-BFGS evaluations, only the per-iteration 𝑱 -dependent SpMM is recomputed, while the underlying 𝑨−1 representation persists—this is the operative mechanism that keeps the adjoint direction stable through transitions where the baselines stall.

Routing Tendon: muscle-driven articulation. A vertical beam is actuated through 16 muscle groups whose time-invariant activations are optimised so that the tip (red dot) reaches a target (blue dot) at frame 100. Two distinct fibre layouts, visible as pink/green colour stripes (Fig. 18), define different muscle-energy routings

Lai, et al.

457× worse than DiffPhD. Per-evaluation, DiffPhD is 3.11× faster than DiffPD and 3.32× faster than MAS on the forward+backward combined wall-clock, despite the baselines stalling at much higher losses. The muscle-energy backward path is the entire bottleneck on this contact-free scene—no contact set to update, no Delassus to rebuild—so the speed-up isolates the unified GPU pipeline (Sec. 4.6) applied to the muscle-actuation Jacobian 𝜕𝒃/𝜕𝒂, which dominates the backward cost through the prox-map differential.

Figure 11: Gatorman: robust mesh-to-mesh contact. A composite visualization showing four temporal snapshots of the Gatorman character executing a strike. The Gatorman wields a stiff sword (6 × 106 Pa) to hit a target sphere of reference stiffness (105 Pa), involving complex mesh-to-mesh interactions and significant material contrast. Contact set size: 24%

Contact set size: 38%

Contact set size: 50%

Homogeneous

Obstacle

Contact set size: 6%

Heterogeneous

Light: 1x stiffness Dark: 0.1x stiffness

Figure 12: Heterogeneous codimensional contact. We simulate a one-layer napkin (50 × 50 × 1 voxels) draped over four spherical obstacles of increasing sizes (top row; contact set sizes from 6% to 50%). Middle row: A homogeneous napkin exhibiting uniform and symmetric deformation. Bottom row: A heterogeneous napkin where the light half retains the baseline stiffness (1×) and the dark half is significantly softer (0.1×). The softer material exhibits more pronounced curling and wraps more tightly around the obstacle, resulting in a highly asymmetric drape profile. Detailed performance metrics and running times are reported in Table 4.

through the same volume. Starting from random activations (Fig. 18, left two panels), DiffPhD converges in 42 L-BFGS iterations to a controller that bends the beam onto the target (Fig. 18, right two panels) at a final loss of 0.441 (Tab. 5). Both DiffPD and MAS fail to converge on this scene: their L-BFGS line search diverges within the first 20 evaluations and the loss plateaus at 201.674—roughly

5.4.3 Trajectory Optimisation. We evaluate trajectory optimisation on two scenes that together stress continuous-contact gait control and muscle-routed articulation: the Torus locomotion benchmark of DiffPD [Du et al. 2021] (homogeneous, contact-rich) and the Routing Tendon muscle-driven beam (homogeneous, contact-free). Both serve as homogeneous controls: where prior differentiable solvers already converge, DiffPhD should match their solutions while being faster, isolating the contribution of the unified GPU pipeline (Sec. 4.6) on the contact-bound and contact-free regimes. Torus: continuous-contact locomotion. The Torus is driven forward by per-element segmental actuation over a 200-frame rollout against the ground, and L-BFGS optimises the activation profile to maximise forward displacement. Starting from a random initialguess actuation that only produces uncoordinated rolling (Fig. 19, top row), DiffPhD recovers a peristaltic gait in 30 L-BFGS iterations whose displacement at frame 200 is 4.1× that of the initial guess (Fig. 19, bottom row), matching the final loss of DiffPD’s reference solution within 1% and running almost the same time as DiffPD’s. The speed-up is smaller than on the static heterogeneity benchmarks of Sec. 5.2 because per-frame contact resolution dominates the wall-clock—a regime where the persistent 𝑺𝑇 𝑺 factor (Sec. 4.6) amortises the global solve but does not eliminate the new-contact-pair work, a finding consistent with [Du et al. 2021].

5.5

Real2Sim: Robot Manipulator

To demonstrate transfer to real-world robotics, we evaluate DiffPhD on two soft-contact manipulation scenarios: (i) Oreo — a Google Robot arm closing on a five-body stack with 1013 Pa rigid claws, a 108 Pa intermediate ball, and a 4.55×105 Pa soft body — and (ii) Dice — a UR5 arm poking a soft die against a tabletop (Fig. 20). Both scenes are simulated inside the SimplerEnv framework [Li et al. 2024]; Oreo acts as a capability check for our contact and material model under extreme stiffness contrast (≥ 108 × across five bodies), while Dice serves as the quantitative Real2Sim benchmark with an explicit optimisation target. Oreo: extreme-contrast capability check. Oreo is reported as a forward-only demonstration: the goal is to show that DiffPhD’s contact-aware persistent factor (Sec. 4.6) remains numerically stable across the five-body contact configuration, not to optimise a downstream loss. A single forward rollout completes in 932 s on our hardware; the recovered deformation is shown in Fig. 20 (left). We verified two adjacent soft-body moduli (4.55×105 and 5.00×105 Pa) and both converge to clean contact without claw–soft interpenetration. Neither DiffPD nor MAS completes a forward step on this scene because the PD weight assembly 𝑤𝑒 ∝ 𝜇𝑒 (Sec. 4.2) develops a condition number beyond their PCG tolerance at the 1013 Pa rigid-claw weights.

Groundtruth Optimized Initial guess

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Homogeneous

Heterogeneous

Figure 13: Bouncing Ball system identification: homogeneous vs. heterogeneous variants. Three rows show three time instants of the rolling/bouncing trajectory. Top row: initial guess. Middle row: DiffPhD recovery. Bottom row: ground truth. Left three columns—homogeneous variant: a single-material green ball whose Young’s modulus and Poisson’s ratio are recovered from a 100-frame trajectory; the optimized trajectory tracks the ground-truth pose closely. Right three columns—heterogeneous variant: the ball is partitioned into three sectors (red, blue, pink) with distinct (𝐸𝑖 , 𝜈𝑖 ) at a ground-truth contrast of {1×, 5×, 10×} in Young’s modulus, and DiffPhD jointly recovers all six parameters; per-sector deformation under impact reveals the heterogeneous material distribution. Table 5: Inverse-problem benchmark across system identification, initial-state optimisation, and trajectory optimisation. For each scene we report forward and backward wall-clock per L-BFGS evaluation (s), number of evaluations (Eval) to convergence, and final loss. The two right-most columns summarise DiffPhD’s per-evaluation (forward+backward) speed-up over each baseline. Forward tolerance 10−3 . Best per scene in bold. DiffPD [Du et al. 2021] Scene Plant Bouncing Ball Bunny (het.) Routing Tendon (het.) Torus

MAS [Wu et al. 2022]

DiffPhD (Ours)

Speed-up

forward

backward

Eval

loss

forward

backward

Eval

loss

forward

backward

Eval

loss

vs. DiffPD

vs. MAS

35.76 45.15 6.71 6.78 7.19

152.73 20.71 11.87 38.65 58.47

15 38 23 20 16

0.692 0.136 0.403 201.674 −0.146

35.78 46.30 6.66 6.65 7.08

156.78 21.36 23.13 41.86 1026.26

15 38 31 20 6

0.692 0.136 0.374 201.674 −0.145

32.56 14.65 6.16 8.77 9.43

28.74 70.62 4.55 5.85 62.77

23 30 30 42 35

0.029 17.149 0.132 0.441 −0.166

3.07× 0.77× 1.74× 3.11× 0.91×

3.13× 0.79× 2.78× 3.32× 14.41×

Dice: Real2Sim contact-location recovery. For the Real2Sim recovery we drive the dataset — PokeFlex [Obrist et al. 2025] recording of a UR5 arm pressing a soft die — and use DiffPhD to recover the in-plane contact location s = (𝑠𝑥 , 𝑠 𝑦 ) between the probe and the die. Starting from a deliberately misaligned probe pose, the gradient-based search drives the simulated die so that the actuator lands on the true contact point and the resulting deformation matches the captured one (Fig. 20, bottom row); the alignment loss collapses smoothly over a handful of iterations (Fig. 20, top right). End-to-end, DiffPhD converges in 13 L-BFGS evaluations at 6.51 s

forward + 56.02 s backward per evaluation, an average of 62.5 s total per evaluation and roughly 13.5 minutes wall-clock for the full recovery. The same pipeline supports per-zone heterogeneous materials on a single body, so the recovered scene can directly host downstream manipulation tasks with spatially varying compliance. This demonstrates that DiffPhD closes the real-to-sim gap via end-to-end gradient-based optimisation, producing digital twins whose contact behaviour is calibrated against real video rather than hand-tuned.

Optimized

Initial guess

Lai, et al.

Figure 16: Bunny inverse design. The initial guess (left) launches the bunny well past the target (blue dot); the optimised initial pose and velocity (right), recovered through DiffPhD’s differentiable Neo-Hookean backward pass, brings the bunny’s centre of mass to the target after 100 frames. Trust-region eigenvalue filtering is engaged on ∼31% of backward solves and stabilises the gradient through the bouncing contact transitions.

0.21

1.30

|Euler angles|

Figure 14: Plant: spatially-localised system identification. A potted plant with non-uniform branch stiffness is observed across four time instants (columns). Rows: initial guess (top), DiffPhD recovery (middle), ground truth (bottom). Highlighted regions (yellow: lower leaves; blue: upper branches) show the two sub-regions whose Young’s moduli are jointly recovered. The optimised trajectory tracks the ground truth in both regions simultaneously.

Optimized

|initial position|

Groundtruth

Initial guess

1.25 1.20

Poisson's ratio estimate

loss

0.45

103 101

loss

10−1

0.35

function evaluations

Ground truth

function evaluations DiffPD

MAS

function evaluations

4.25

function evaluations Ours

function evaluations

100

4.00 3.75 3.50 0

Figure 15: Plant system identification: convergence curves. Per-region Young’s modulus and loss versus L-BFGS evaluations on the Plant benchmark.

5.6

4.50

10−3

0.30

loss

0.40

106

1.15

|initial velocity|

Young's modulus estimate

Poisson's ratio

Young's modulus (Pa)

0.20

Ablation Study

We isolate two architectural axes of DiffPhD on the heterogeneous Armadillo (twist scene, forward-simulation stability) and the heterogeneous Neo-Hookean Bouncing Ball (system identification): (a) the choice of eigenvalue projection on the per-element prox-map Hessian (none as in DiffPD [Du et al. 2021], fixed clamping [Teran et al. 2005], fixed absolute-value [Chen et al. 2024b], or our stateadaptive trust-region rule of Sec. 4.4); and (b) GPU residency via the unified local–global pipeline of Sec. 4.6. The first three rows of Table 6 isolate projection alone on the CPU pipeline (DiffPD’s setting); the final row enables both GPU residency and the adaptive filter, giving the full DiffPhD configuration. Findings. Three observations emerge from Table 6. (1) Projection alone is not enough; GPU residency is what unlocks the speed-up. On the CPU pipeline, swapping the projection from none (DiffPD) to fixed clamp or abs cuts backward time

function evaluations

30 DiffPD

0 MAS

function evaluations

Ours

Figure 17: Bunny inverse design: parameter and loss convergence. Recovered initial position, initial Euler angles, and initial velocity norms (first three panels) and loss (right) versus L-BFGS evaluations. DiffPhD reaches a 3.05× lower final loss than DiffPD and 2.83× lower than MAS; MAS (red) stalls earliest because its multilevel preconditioner must refactorize on every active-set change and loses subdomain conditioning under heterogeneity. roughly in half on Armadillo (140 s → 66–73 s) but adds backward cost on Bouncing Ball (4.3 s → 11 s), and leaves forward time unchanged in both. Projection on its own is a gradient-quality choice, not a speed knob. The full DiffPhD row collapses forward time by a further 2.3× on Armadillo and 2.5× on Bouncing Ball, almost entirely from the unified GPU local–global pipeline (Sec. 4.6); without it, every CPU projection variant remains within ±5% of the DiffPD forward time. (2) The state-adaptive trust-region filter unlocks distinct optima. The three CPU rows converge to identical losses on both

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Initial guess

frame 1

Optimized

frame 100

frame 1

frame 100

Optimized

Initial guess

Figure 18: Routing Tendon: muscle-driven initial-state optimisation. A vertical beam with two fibre layouts (visible as pink/green colour stripes) is actuated through 16 muscle groups to bring its tip (red dot) to a target (blue dot) at frame 100. Left two panels: initial guess with random activations. Right two panels: the optimised controller after 30 L-BFGS iterations through DiffPhD’s prox-map differential.

frame 1

frame 100

frame 150

frame 200

Figure 19: Torus locomotion: trajectory optimisation. The Torus rolls forward under per-element segmental actuation. Top row: initial-guess actuation produces uncoordinated oscillation; the torus only travels a short distance over 200 frames. Bottom row: after 30 L-BFGS iterations through DiffPhD’s differentiable backward pass, the optimised actuation yields a peristaltic gait with 4.1× longer forward displacement at frame 200.

scenes (Armadillo −3.544; Bouncing Ball 142.967): fixed projection schemes preserve DiffPD’s optimisation landscape and reach the same fixed point. The adaptive filter changes the landscape that L-BFGS traverses—reaching a Bouncing Ball loss of 36.6 (a 3.9× reduction) while moving Armadillo to a different basin (+58.84, a saddle reachable only through the differentiable backward path that the adaptive filter stabilises). On scenes where the prox-map Hessian alternates between near-convex (clamping wins) and stiffly non-convex (abs wins) across the trajectory, the state-adaptive rule

of Eq. (29) selects per-element 𝜏 ∗ at zero additional cost: ∆Φmod is one SpMV reusing the persistent factor. (3) Each axis is independently necessary. Removing the adaptive filter (rows 2–3) preserves DiffPD’s fixed point but does not reach the lower-loss basin DiffPhD finds. Removing GPU residency (rows 1–3) leaves forward time at the DiffPD baseline. Stacking both axes yields the 2.3–2.5× forward and 6.3× backward end-to-end speed-up of the full DiffPhD over the strongest CPU-only ablation on Armadillo, with comparable improvements on Bouncing Ball.

6 DISCUSSION 6.1 Failure Modes Live at the Intersection of Properties The most consequential numerical failures we encountered while building DiffPhD were silent at the level of any single property and only became observable when three properties were exercised simultaneously: heterogeneous Neo-Hookean elasticity, Anderson acceleration on the outer fixed-point loop, and frictional contact through the Fischer–Burmeister NCP. Per-element projection passes were stable in homogeneous meshes; Anderson mixing was stable in convex-energy regimes; the NCP contact solver was stable at zero restitution. Yet the heterogeneous NH forward solver developed a mesh-scale instability concentrated at the 1–2% of elements straddling the stiff/soft interface, visible only when the indefinite stretch-space Hessian (Eq. (15)) and Anderson Type-II mixing co-occurred at the same iterate. The implication is methodological rather than mechanical: property-by-property unit tests cannot certify a unified differentiable solver, because each component’s correctness assumption can be silently violated by another component’s regime. We argue this is a generic characteristic of unified differentiable physics—the contribution-by-contribution validation discipline that has served the forward-simulation literature does not transfer to solvers that share state across heterogeneity, hyperelasticity, and contact.

6.2

Trust-Region as a Two-Sided Safety Net prox

¯ (Eq. (23)) The Neo-Hookean prox-map Hessian 𝑯𝑒 = 𝑯𝜓 (𝝈 ∗ )+𝑘𝑰 ¯ < 𝜆¯ ln Q 𝑗 𝜎 𝑗 —geometrically, goes indefinite whenever 𝜎𝑖2 (𝜇¯ + 𝜆) near element inversion or under high-Poisson compression. This indefiniteness needs filtering on both passes of the solver, but for asymmetric reasons. On the forward pass, the trust-region cap bounds the inner Newton step in stretch space so that the proxmap output 𝝈 ∗ (𝒒) remains a Lipschitz function of 𝒒—the fixed-point contraction condition that the Anderson Type-II mixing formula requires (Sec. 4.3). On the backward pass, the 𝜌-adaptive blend 𝑯˜ = (1−𝜏)𝑯 +𝜏 |𝑯 | (Sec. 4.4) ensures the implicit-function-theorem inverse used in the adjoint linear system stays in a descent cone for the loss gradient. The unifying insight is that eigenvalue-aware regularisation is needed only where indefiniteness is actually exercised: per-element, per-iterate, and only on the offending block. Applying the same treatment globally—PSD-projecting the assembled 𝑨, the most common defensive move in differentiable PD—perturbs the entire fixed-point map and either desynchronises Anderson on the

Lai, et al.

Heterogeneous Soft Contact

Hard

Real World

Real-to-Sim Alignment Loss

Ground Truth

Optimized Result

Hard Soft

Initial state

Optimized

Hard

Soft

Figure 20: Real-to-Sim alignment for a soft-contact manipulation pipeline. Left — Heterogeneous Soft Contact. A simulated Google Robot gripper closes on a foam ball with an A–B–A material layout (hard / soft / hard along the closing axis). DiffPhD’s bidirectional mesh contact resolves the per-zone deformation without introducing a visible gap between the gripper fingertips and the ball surface, illustrating that the same framework supports spatially varying compliance on a single body. Right — Real-to-Sim Alignment. From a PokeFlex [Obrist et al. 2025] recording of a UR5 arm pressing a soft die, DiffPhD recovers the in-plane contact location s = (𝑠𝑥 , 𝑠 𝑦 ) between the probe and the die. The top row compares the captured deformation field (Ground Truth) with our reconstruction (Optimized Result); the convergence plot on the far right shows the alignment loss decreasing smoothly over a handful of iterations. The bottom row visualises the search: starting from a deliberately misaligned probe pose (Initial state), the gradient-based optimisation moves the simulated die so that the probe lands on the true contact point and the deformation matches the real capture (Optimized). Table 6: Ablation study. Columns indicate which components are enabled: GPU = unified GPU local-global pipeline; Proj. = per-element prox-map eigenvalue projection (none/abs/clmp/adapt); Right columns report end-to-end optimisation time (s) and final loss on two heterogeneous benchmarks.

Configuration

Armadillo

Method

GPU

Proj.

forward (s)

DiffPD [Du et al. 2021] DiffPD w/ clamping DiffPD w/ abs DiffPhD

× × × ✓

none clmp abs adapt

191.845 195.192 194.974 84

backward (s)

loss

forward (s)

backward (s)

loss

140.203 66.272 72.875 22.263

-3.544 -3.544 -3.544 58.84

25.616 24.816 25.616 10.092

4.265 11.106 11.142 49.204

142.965 142.967 142.967 36.565

forward pass or biases the gradient on the backward pass. The locality of the trust-region decision is what lets the same machinery serve both passes without these side-effects.

6.3

Bouncing Ball (Neo-Hookean)

Aggregate Metrics Are Insufficient for Gradient Correctness

A class of numerical bugs we surfaced repeatedly during development was invisible to every aggregate metric—per-frame PD

residual, contact-set convergence, visual mesh fidelity—and only became observable through per-iteration invariants designed around a specific component’s assumption. The batched Delassus assembly 𝑾 = 𝑱 𝑨−1 𝑱 𝑇 on the GPU degenerated to a near-diagonal matrix under certain SpMM-algorithm selections, leaving forward simulation visually intact but silently destroying the off-diagonal contact coupling the backward gradient depends on; the bug was only diagnosable through the per-frame count of active contacts whose Delassus diagonal magnitude exceeded a fixed threshold, paired

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

with a deliberately slow CPU fallback used as an A/B oracle. The lesson generalises: a differentiable solver can produce a plausible loss trajectory whose gradient is wrong by a constant scale—the worstcase failure for inverse design, because the optimisation will still descend monotonically into the wrong basin. Future differentiablephysics work needs counter-based component-local invariants— active-contact count, prox-map ping-pong amplitude, per-element Hessian eigenvalue spread—as first-class correctness signals during development, not just end-to-end loss curves at evaluation time.

7 LIMITATIONS AND FUTURE WORK 7.1 Velocity-Level NCP Drifts in the Quasi-Static Limit The Fischer–Burmeister formulation we adopt (Sec. 4.7) enforces complementarity 𝝀𝑛 ⊥ 𝜹𝑛 on the velocity-level contact reaction— the natural choice for dynamic impact—but does not explicitly conserve position-level non-penetration across the entire rollout. At restitution 𝑒 → 0, the post-correction is triggered only while 𝒗 ·𝑛 < 0 ∧ 𝝀 > 0; once a body settles into static support, neither condition fires robustly, and configuration drifts ∼1–2 cm below the nominal contact plane over 50 frames at our standard ℎ = 4 ms timestep. Our bouncing-ball benchmarks use 𝑒 = 0.3, which is physically appropriate for soft rubber and keeps the solver inside the regime where velocity-level enforcement is sufficient. Extending DiffPhD to manipulation scenarios that include sustained static contact—grasping under gravity, multi-finger pinch closure, persistent surface support—requires augmenting the NCP with a positionlevel projection coupled to the PD global step, or replacing it with an impulse-integrated ground reaction that accumulates the contact force budget. The persistent 𝑺𝑇 𝑺 factor of Sec. 4.6 is orthogonal to this choice and would carry over directly, so we view the gap as engineering rather than fundamental.

7.2

Material-Parameter Backward Is Energy-Specific

The chain rule from projective weights {𝑤𝑒 } to physical material parameters (𝐸𝑒 , 𝜈𝑒 ) is currently derived analytically per energy: a corotated+volume composite uses 𝑤 0 = 2𝜇, 𝑤 1 = 𝜆 (Eq. (1) chain rule), while the single-energy Neo-Hookean form uses 𝑤𝑒 = 2𝜇 + 𝜆. The geometry- and initial-state inverse problems in our evaluation (Bunny, Routing Tendon, Torus) do not exercise this Jacobian— they optimise quantities that flow through 𝜕𝒃/𝜕𝒂, not through 𝜕𝑨/𝜕𝑬—so the material backward path is exercised only by the system-identification benchmarks (Bouncing Ball, Plant). Extending DiffPhD to material inverse design across the full suite of supported energies—anisotropic fibre, viscoelastic, rate-dependent plastic, or any future hyperelastic energy that fits the PD prox-map formulation— requires re-deriving 𝜕𝑤𝑒 /𝜕(𝐸𝑒 , 𝜈𝑒 ) for each. The natural future direction is automatic symbolic differentiation of the prox-map energy itself, so that adding a new energy to the framework automatically extends the material inverse design path; this is a refactor we sketch but do not implement.

CONCLUSION

We presented DiffPhD, a unified differentiable solver that brings heterogeneous, hyperelastic, contact-rich elastodynamics within reach of end-to-end gradient-based optimisation. Our three contributions compose around a single representational invariant—a persistent 𝑺𝑇 𝑺 factor of the PD operator that absorbs material contrast at assembly: stiffness-aware projective assembly that routes per-element 𝜇𝑒 into 𝑨 (Sec. 4.2), a state-adaptive trust-region filter on the per-element prox-map Hessian (Sec. 4.4), and a unified GPU local–global pipeline that shares the same factor across forward, contact, and backward stages (Sec. 4.6). Each contribution addresses an independent failure mode of prior differentiable PD: the forward stiffness–backward gradient coupling that destabilises DiffPD beyond ∼ 50× contrast, the indefinite prox-map Hessian that derails L-BFGS line search under hyperelastic contact, and the per-evaluation refactorisation cost that limits mesh size and contact density. Across seventeen benchmarks, DiffPhD remains stable up to 100× stiffness contrast on the Armadillo where DiffPD diverges, delivers 8.69× forward speed-up on a 2.2M-DoF heterogeneous Crab and 33.31× on dense codimensional contact, and reaches a 23.9× smaller residual loss than the strongest baseline on slendergeometry system identification. On the inverse-problem suite, DiffPhD converges to lower-loss optima than DiffPD and MAS in every scene where the baselines also converge, and to a finite optimum where they diverge (Routing Tendon, 457×). The Real2Sim manipulation pipeline completes in 13.5 minutes wall-clock, taking the framework from forward simulator into the regime where differentiable physics becomes a practical building block for robotic manipulation under heterogeneous frictional contact. We view DiffPhD as evidence for a broader claim: in differentiable elastodynamics, the design target should be operator amortisation across forward and backward stages, not arithmetic throughput on either pass alone. The contributions we present are specific to projective dynamics, but the principle of routing material structure into the operator at assembly, filtering indefiniteness adaptively at the per-element block, and sharing a single factor across all stages of a differentiable pipeline is general—and we believe it applies equally to Newton-based differentiable solvers, to position-based dynamics, and to differentiable contact formulations beyond the NCP we adopt here.

References Jeremie Allard, St ’ephane Cotin, Francois Faure, P-J Bensoussan, François Poyer, Christian Duriez, Herv ’e Delingette, and Laurent Grisoni. 2007. SOFA - an Open Source Framework for Medical Simulation. Studies in health technology and informatics 125 (02 2007), 13–8. Jernej Barbič and Doug L. James. 2005. Real-Time subspace integration for St. VenantKirchhoff deformable models. ACM Trans. Graph. 24, 3 (July 2005), 982–990. https: //doi.org/10.1145/1073204.1073300 Jernej Barbič, Funshing Sin, and Eitan Grinspun. 2012. Interactive editing of deformable simulations. ACM Trans. Graph. 31, 4, Article 70 (July 2012), 8 pages. https: //doi.org/10.1145/2185520.2185566 Otman Benchekroun, Eitan Grinspun, Maurizio Chiaramonte, and Philip Allen Etter. 2025. Force-Dual Modes: Subspace Design from Stochastic Forces. ACM Trans. Graph. 44, 6, Article 200 (Dec. 2025), 14 pages. https://doi.org/10.1145/3763310 Otman Benchekroun, Kaixiang Xie, Hsueh-Ti Derek Liu, Eitan Grinspun, Sheldon Andrews, and Victor Zordan. 2024. Actuators A La Mode: Modal Actuations for Soft Body Locomotion. In SIGGRAPH Asia 2024 Conference Papers (Tokyo, Japan)

Lai, et al.

(SA ’24). Association for Computing Machinery, New York, NY, USA, Article 53, 10 pages. https://doi.org/10.1145/3680528.3687638 James M Bern, Pol Banzet, Roi Poranne, and Stelian Coros. 2019. Trajectory Optimization for Cable-Driven Soft Robot Locomotion.. In Robotics: Science and Systems, Vol. 1. Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, and Mark Pauly. 2014. Projective dynamics: fusing constraint projections for fast simulation. ACM Trans. Graph. 33, 4, Article 154 (July 2014), 11 pages. https://doi.org/10.1145/2601097. 2601116 Christopher Brandt, Elmar Eisemann, and Klaus Hildebrandt. 2018. Hyper-reduced projective dynamics. ACM Trans. Graph. 37, 4, Article 80 (July 2018), 13 pages. https://doi.org/10.1145/3197517.3201387 Robert Bridson, Ronald Fedkiw, and John Anderson. 2002. Robust treatment of collisions, contact and friction for cloth animation. ACM Trans. Graph. 21, 3 (July 2002), 594–603. https://doi.org/10.1145/566654.566623 Xiao-Chuan Cai and Yousef Saad. 1996. Overlapping Domain Decomposition Algorithms for General Sparse Matrices. Numer. Linear Algebra Appl. 3 (1996), 221– 237. https://doi.org/10.1002/(SICI)1099-1506(199605/06)3:3<221::AID-NLA80>3.0. CO;2-7 Xiao-Chuan Cai and Marcus Sarkis. 1999. A Restricted Additive Schwarz Preconditioner for General Sparse Linear Systems. SIAM J. Sci. Comput. 21, 2 (Jan. 1999), 792–797. https://doi.org/10.1137/S106482759732678X Youshuai Cai, Chen Li, Haichuan Song, Youchen Xie, and ChangBo Wang. 2025. DiffQN: Differentiable Quasi-Newton Method for Elastodynamics. In Pacific Graphics Conference Papers, Posters, and Demos, Marc Christie, Ping-Hsuan Han, Shih-Syun Lin, Nico Pietroni, Teseo Schneider, Hsin-Ruey Tsai, Yu-Shuen Wang, and Eugene Zhang (Eds.). The Eurographics Association. https://doi.org/10.2312/pg.20251265 Desai Chen, David I. W. Levin, Wojciech Matusik, and Danny M. Kaufman. 2017. Dynamics-aware numerical coarsening for fabrication design. ACM Trans. Graph. 36, 4, Article 84 (July 2017), 15 pages. https://doi.org/10.1145/3072959.3073669 Honglin Chen, Hsueh-Ti Derek Liu, Alec Jacobson, David I.W. Levin, and Changxi Zheng. 2024a. Trust-Region Eigenvalue Filtering for Projected Newton. In SIGGRAPH Asia 2024 Conference Papers (Tokyo, Japan) (SA ’24). Association for Computing Machinery, New York, NY, USA, Article 120, 10 pages. https://doi.org/10. 1145/3680528.3687650 Honglin Chen, Hsueh-Ti Derek Liu, David IW Levin, Changxi Zheng, and Alec Jacobson. 2024b. Stabler Neo-Hookean Simulation: Absolute Eigenvalue Filtering for Projected Newton. In ACM SIGGRAPH 2024 Conference Papers. 1–10. Jiong Chen, Florian Schäfer, Jin Huang, and Mathieu Desbrun. 2021. Multiscale Cholesky preconditioning for ill-conditioned problems. ACM Trans. Graph. 40, 4, Article 81 (July 2021), 13 pages. https://doi.org/10.1145/3450626.3459851 Nick Cheney, Robert MacCurdy, Jeff Clune, and Hod Lipson. 2013. Unshackling evolution: evolving soft robots with multiple materials and a powerful generative encoding. In Proceedings of the 15th Annual Conference on Genetic and Evolutionary Computation (Amsterdam, The Netherlands) (GECCO ’13). Association for Computing Machinery, New York, NY, USA, 167–174. https://doi.org/10.1145/2463372.2463404 Francesco Corucci, Nick Cheney, Hod Lipson, Cecilia Laschi, and Josh Bongard. 2016a. Evolving swimming soft-bodied creatures. In ALIFE XV, the fifteenth international conference on the synthesis and simulation of living systems, late breaking proceedings, Vol. 6. Francesco Corucci, Nick Cheney, Hod Lipson, Cecilia Laschi, and Josh Bongard. 2016b. Material properties affect evolution’s ability to exploit morphological computation in growing soft-bodied creatures. https://doi.org/10.7551/978-0-262-33936-0-ch043 Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio. 2014. Identifying and attacking the saddle point problem in high-dimensional non-convex optimization. In Advances in Neural Information Processing Systems, Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger (Eds.), Vol. 27. Curran Associates, Inc. https://proceedings.neurips.cc/ paper_files/paper/2014/file/04192426585542c54b96ba14445be996-Paper.pdf Christian Dick, Joachim Georgii, and Rüdiger Westermann. 2011. A real-time multigrid finite hexahedra method for elasticity simulation using CUDA. Simulation Modelling Practice and Theory 19, 2 (2011), 801–816. https://doi.org/10.1016/j.simpat.2010.11. 005 Victorita Dolean, Pierre Jolivet, and Fr ’ed ’eric Nataf. 2015. An Introduction to Domain Decomposition Methods. Society for Industrial and Applied Mathematics, Philadelphia, PA. https://doi.org/10.1137/1. 9781611974065 arXiv:https://epubs.siam.org/doi/pdf/10.1137/1.9781611974065 Maksymilian Dryja and Olof B. Widlund. 1989. Towards a unified theory of domain decomposition algorithms for elliptic problems. Technical Report 486. Department of Computer Science, Courant Institute of Mathematical Sciences, New York University, New York, NY. Maksymilian Dryja and Olof B. Widlund. 1990. Multilevel additive methods for elliptic finite element problems. Technical Report 507. Department of Computer Science, Courant Institute of Mathematical Sciences, New York University, New York, NY. 58–69 pages. Maksymilian Dryja and Olof B. Widlund. 1992. Additive Schwarz methods for elliptic finite element problems in three dimensions. In Fifth International Symposium on

Domain Decomposition Methods for Partial Differential Equations, Tony F. Chan, David E. Keyes, Gérard A. Meurant, Jeffrey S. Scroggs, and Robert G. Voigt (Eds.). SIAM, Philadelphia, PA, 3–18. Tao Du, Kui Wu, Pingchuan Ma, Sebastien Wah, Andrew Spielberg, Daniela Rus, and Wojciech Matusik. 2021. DiffPD: Differentiable Projective Dynamics. ACM Trans. Graph. 41, 2, Article 13 (Nov. 2021), 21 pages. https://doi.org/10.1145/3490168 Essex Edwards and Robert Bridson. 2015. The Discretely-Discontinuous Galerkin Coarse Grid for Domain Decomposition. arXiv:1504.00907 [cs.NA] https://arxiv. org/abs/1504.00907 Yalchin Efendiev, Juan Galvis, Raytcho Lazarov, and Joerg Willems. 2012. Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms. ESAIM: Mathematical Modelling and Numerical Analysis 46, 5 (2012), 1175–1199. https://doi.org/10.1051/m2an/2011073 Charbel Farhat, Michel Lesoinne, Patrick LeTallec, Kendall Pierson, and Daniel Rixen. 2001. FETI-DP: a dual–primal unified FETI method—part I: A faster alternative to the two-level FETI method. Internat. J. Numer. Methods Engrg. 50, 7 (2001), 1523–1544. https://doi.org/10.1002/nme.76 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.76 Charbel Farhat and Francois-Xavier Roux. 1991. A method of finite element tearing and interconnecting and its parallel solution algorithm. Internat. J. Numer. Methods Engrg. 32, 6 (1991), 1205–1227. https://doi.org/10.1002/nme.1620320604 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620320604 Zachary Ferguson, Teseo Schneider, Danny Kaufman, and Daniele Panozzo. 2023. In-Timestep Remeshing for Contacting Elastodynamics. ACM Trans. Graph. 42, 4, Article 145 (July 2023), 15 pages. https://doi.org/10.1145/3592428 A. Fischer. 1992. A special newton-type optimization method. Optimization 24, 3-4 (1992), 269–284. https://doi.org/10.1080/02331939208843795 arXiv:https://doi.org/10.1080/02331939208843795 Marco Fratarcangeli, Valentina Tibaldo, and Fabio Pellacini. 2016. Vivace: a practical gauss-seidel method for stable soft body dynamics. ACM Trans. Graph. 35, 6, Article 214 (Dec. 2016), 9 pages. https://doi.org/10.1145/2980179.2982437 Andreas Frommer and Daniel B. Szyld. 1999. Weighted max norms, splittings, and overlapping additive Schwarz iterations. Numer. Math. 83, 2 (1999), 259–278. https: //doi.org/10.1007/s002110050449 Xiao-Ming Fu and Yang Liu. 2016. Computing inversion-free mappings by simplex assembly. ACM Trans. Graph. 35, 6, Article 216 (Dec. 2016), 12 pages. https: //doi.org/10.1145/2980179.2980231 Moritz Geilinger, David Hahn, Jonas Zehnder, Moritz Bächer, Bernhard Thomaszewski, and Stelian Coros. 2020. ADD: analytically differentiable dynamics for multi-body systems with frictional contact. ACM Trans. Graph. 39, 6, Article 190 (Nov. 2020), 15 pages. https://doi.org/10.1145/3414685.3417766 Philip E Gill, Walter Murray, and Margaret H Wright. 2019. Practical optimization. SIAM. G. Haase, U. Langer, and A. Meyer. 1991. The approximate Dirichlet domain decomposition method. Part I: an algebraic approach. Computing 47, 2 (July 1991), 137–151. https://doi.org/10.1007/BF02253431 David Hahn, Pol Banzet, James M. Bern, and Stelian Coros. 2019. Real2Sim: visco-elastic parameter estimation from dynamic motion. ACM Trans. Graph. 38, 6, Article 236 (Nov. 2019), 13 pages. https://doi.org/10.1145/3355089.3356548 Philipp Herholz and Olga Sorkine-Hornung. 2020. Sparse cholesky updates for interactive mesh parameterization. ACM Trans. Graph. 39, 6, Article 202 (Nov. 2020), 14 pages. https://doi.org/10.1145/3414685.3417828 Jonathan Hiller and Hod Lipson. 2014. Dynamic simulation of soft multimaterial 3d-printed objects. Soft robotics 1, 1 (2014), 88–101. Yuanming Hu, Luke Anderson, Tzu-Mao Li, Qi Sun, Nathan Carr, Jonathan RaganKelley, and Fredo Durand. 2020. DiffTaichi: Differentiable Programming for Physical Simulation. In International Conference on Learning Representations. https: //openreview.net/forum?id=B1eB5xSFvr Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua B Tenenbaum, William T Freeman, Jiajun Wu, Daniela Rus, and Wojciech Matusik. 2019. ChainQueen: A Real-Time Differentiable Physical Simulator for Soft Robotics. Proceedings of IEEE International Conference on Robotics and Automation (ICRA) (2019). Zizhou Huang, Maxwell Paik, Zachary Ferguson, Daniele Panozzo, and Denis Zorin. 2025. Geometric Contact Potential. ACM Trans. Graph. 44, 4, Article 158 (July 2025), 24 pages. https://doi.org/10.1145/3731142 Zizhou Huang, Davi Colli Tozoni, Arvi Gjoka, Zachary Ferguson, Teseo Schneider, Daniele Panozzo, and Denis Zorin. 2024. Differentiable solver for time-dependent deformation problems with contact. ACM Trans. Graph. 43, 3, Article 31 (May 2024), 30 pages. https://doi.org/10.1145/3657648 T.J.R. Hughes. 2012. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover Publications. https://books.google.com.tw/books?id= cHH2n_qBK0IC Christian Kanzow and Daniel Steck. 2023. Regularization of limited memory quasiNewton methods for large-scale nonconvex minimization. Mathematical Programming Computation 15, 3 (2023), 417–444. George Karypis and Vipin Kumar. 1997. METIS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices. Technical Report. University of Minnesota, Department of

DiffPhD: A Unified Differentiable Solver for Projective Heterogeneous Materials in Elastodynamics with Contact-Rich GPU-Acceleration

Computer Science and Engineering. https://hdl.handle.net/11299/215346 George Karypis and Vipin Kumar. 1998. A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs. SIAM J. Sci. Comput. 20, 1 (Dec. 1998), 359–392. Ryo Kikuuwe, Hiroaki Tabuchi, and Motoji Yamamoto. 2009. An edge-based computationally efficient formulation of Saint Venant-Kirchhoff tetrahedral finite elements. ACM Trans. Graph. 28, 1, Article 8 (Feb. 2009), 13 pages. https: //doi.org/10.1145/1477926.1477934 Theodore Kim and David Eberle. 2020. Dynamic deformables: implementation and production practicalities. In ACM SIGGRAPH 2020 Courses (Virtual Event, USA) (SIGGRAPH ’20). Association for Computing Machinery, New York, NY, USA, Article 23, 182 pages. https://doi.org/10.1145/3388769.3407490 Theodore Kim and Doug L. James. 2011. Physics-based character skinning using multi-domain subspace deformations. In Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation (Vancouver, British Columbia, Canada) (SCA ’11). Association for Computing Machinery, New York, NY, USA, 63–72. https://doi.org/10.1145/2019406.2019415 Dilip Krishnan, Raanan Fattal, and Richard Szeliski. 2013. Efficient preconditioning of laplacian matrices for computer graphics. ACM Trans. Graph. 32, 4, Article 142 (July 2013), 15 pages. https://doi.org/10.1145/2461912.2461992 Nomesh Kumar and V. Venkateswara Rao. 2016. Hyperelastic Mooney-Rivlin Model: Determination and Physical Interpretation of Material Constants. MIT International Journal of Mechanical Engineering 6, 1 (2016), 43–46. Lei Lan, Minchen Li, Chenfanfu Jiang, Huamin Wang, and Yin Yang. 2023. Secondorder Stencil Descent for Interior-point Hyperelasticity. ACM Trans. Graph. 42, 4, Article 108 (July 2023), 16 pages. https://doi.org/10.1145/3592104 Lei Lan, Zixuan Lu, Jingyi Long, Chun Yuan, Xuan Li, Xiaowei He, Huamin Wang, Chenfanfu Jiang, and Yin Yang. 2024. Efficient GPU Cloth Simulation with Nondistance Barriers and Subspace Reuse. ACM Trans. Graph. 43, 6, Article 226 (Nov. 2024), 16 pages. https://doi.org/10.1145/3687760 Seunghwan Lee, Ri Yu, Jungnam Park, Mridul Aanjaneya, Eftychios Sifakis, and Jehee Lee. 2018. Dexterous manipulation and control with volumetric muscles. ACM Trans. Graph. 37, 4, Article 57 (July 2018), 13 pages. https://doi.org/10.1145/3197517. 3201330 Minchen Li, Zachary Ferguson, Teseo Schneider, Timothy Langlois, Denis Zorin, Daniele Panozzo, Chenfanfu Jiang, and Danny M. Kaufman. 2020. Incremental potential contact: intersection-and inversion-free, large-deformation dynamics. ACM Trans. Graph. 39, 4, Article 49 (Aug. 2020), 20 pages. https://doi.org/10.1145/ 3386569.3392425 Minchen Li, Ming Gao, Timothy Langlois, Chenfanfu Jiang, and Danny M. Kaufman. 2019a. Decomposed optimization time integrator for large-step elastodynamics. ACM Trans. Graph. 38, 4, Article 70 (July 2019), 10 pages. https://doi.org/10.1145/ 3306346.3322951 Ruipeng Li and Yousef Saad. 2017. Low-Rank Correction Methods for Algebraic Domain Decomposition Preconditioners. SIAM J. Matrix Anal. Appl. 38, 3 (2017), 807–828. https://doi.org/10.1137/16M110486X arXiv:https://doi.org/10.1137/16M110486X Sizhe Li, Zhiao Huang, Tao Du, Hao Su, Joshua Tenenbaum, and Chuang Gan. 2022b. Contact Points Discovery for Soft-Body Manipulations with Differentiable Physics. In International Conference on Learning Representations (ICLR). Xuanlin Li, Kyle Hsu, Jiayuan Gu, Karl Pertsch, Oier Mees, Homer Rich Walke, Chuyuan Fu, Ishikaa Lunawat, Isabel Sieh, Sean Kirmani, Sergey Levine, Jiajun Wu, Chelsea Finn, Hao Su, Quan Vuong, and Ted Xiao. 2024. Evaluating Real-World Robot Manipulation Policies in Simulation. arXiv preprint arXiv:2405.05941 (2024). Yifei Li, Tao Du, Kui Wu, Jie Xu, and Wojciech Matusik. 2022a. DiffCloth: Differentiable Cloth Simulation with Dry Frictional Contact. ACM Trans. Graph. 42, 1, Article 2 (oct 2022), 20 pages. https://doi.org/10.1145/3527660 Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B. Tenenbaum, and Antonio Torralba. 2019b. Learning Particle Dynamics for Manipulating Rigid Bodies, Deformable Objects, and Fluids. In International Conference on Learning Representations. https: //openreview.net/forum?id=rJgbSn09Ym Junbang Liang, Ming Lin, and Vladlen Koltun. 2019. Differentiable Cloth Simulation for Inverse Problems. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc. https://proceedings.neurips.cc/paper_files/paper/2019/file/ 28f0b864598a1291557bed248a998d4e-Paper.pdf Huancheng Lin, Floyd M Chitalu, and Taku Komura. 2022. Isotropic ARAP Energy Using Cauchy-Green Invariants. ACM Trans. Graph. 41, 6 (2022), 275–1. Tiantian Liu, Sofien Bouaziz, and Ladislav Kavan. 2017. Quasi-Newton Methods for Real-Time Simulation of Hyperelastic Materials. ACM Trans. Graph. 36, 4, Article 116a (May 2017), 16 pages. https://doi.org/10.1145/3072959.2990496 Andreas Longva, Fabian Löschner, Jos ’e Antonio Fern ’andez-Fern ’andez, Egor Larionov, Uri M. Ascher, and Jan Bender. 2023. Pitfalls of Projection: A study of Newton-type solvers for incremental potentials. arXiv:2311.14526 [cs.GR] https://arxiv.org/abs/2311.14526 Zixuan Lu, Xiaowei He, Yuzhong Guo, Xuehui Liu, and Huamin Wang. 2024. Projective Peridynamic Modeling of Hyperelastic Membranes With Contact. IEEE Transactions on Visualization and Computer Graphics 30, 8 (Aug. 2024), 4601–4614. https:

//doi.org/10.1109/TVCG.2023.3271511 Zixuan Lu, Ziheng Liu, Lei Lan, Huamin Wang, Yuko Ishiwaka, Chenfanfu Jiang, Kui Wu, and Yin Yang. 2025. High-performance CPU Cloth Simulation Using Domaindecomposed Projective Dynamics. ACM Trans. Graph. 44, 4, Article 51 (July 2025), 17 pages. https://doi.org/10.1145/3731182 Mickaël Ly, Jean Jouve, Laurence Boissieux, and Florence Bertails-Descoubes. 2020. Projective dynamics with dry frictional contact. ACM Trans. Graph. 39, 4, Article 57 (Aug. 2020), 8 pages. https://doi.org/10.1145/3386569.3392396 Pingchuan Ma, Peter Yichen Chen, Bolei Deng, Joshua B Tenenbaum, Tao Du, Chuang Gan, and Wojciech Matusik. 2023. Learning Neural Constitutive Laws from Motion Observations for Generalizable PDE Dynamics. In International Conference on Machine Learning. PMLR. Miles Macklin, Kenny Erleben, Matthias Müller, Nuttapong Chentanez, Stefan Jeschke, and Tae-Yong Kim. 2020. Primal/Dual Descent Methods for Dynamics. Computer Graphics Forum (2020). https://doi.org/10.1111/cgf.14104 Miles Macklin, Matthias Müller, and Nuttapong Chentanez. 2016. XPBD: positionbased simulation of compliant constrained dynamics. In Proceedings of the 9th International Conference on Motion in Games (Burlingame, California) (MIG ’16). Association for Computing Machinery, New York, NY, USA, 49–54. https://doi. org/10.1145/2994258.2994272 Jan Mandel. 1993. Balancing domain decomposition. Communications in Numerical Methods in Engineering 9, 3 (1993), 233–241. https://doi.org/10.1002/cnm.1640090307 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.1640090307 Sebastian Martin, Bernhard Thomaszewski, Eitan Grinspun, and Markus Gross. 2011. Example-based elastic materials. In ACM SIGGRAPH 2011 Papers (Vancouver, British Columbia, Canada) (SIGGRAPH ’11). Association for Computing Machinery, New York, NY, USA, Article 72, 8 pages. https://doi.org/10.1145/1964921.1964967 Sehee Min, Jungdam Won, Seunghwan Lee, Jungnam Park, and Jehee Lee. 2019. SoftCon: simulation and control of soft-bodied animals with biomimetic actuators. ACM Trans. Graph. 38, 6, Article 208 (Nov. 2019), 12 pages. https: //doi.org/10.1145/3355089.3356497 V. Modi, L. Fulton, A. Jacobson, S. Sueda, and D.I.W. Levin. 2021. EMU: Efficient Muscle Simulation in Deformation Space. Computer Graphics Forum 40, 1 (2021), 234–248. https://doi.org/10.1111/cgf.14185 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.14185 Vismay Modi, Nicholas Sharp, Or Perel, Shinjiro Sueda, and David I. W. Levin. 2024. Simplicits: Mesh-Free, Geometry-Agnostic Elastic Simulation. ACM Trans. Graph. 43, 4, Article 117 (July 2024), 11 pages. https://doi.org/10.1145/3658184 M. Mooney. 1940. A Theory of Large Elastic Deformation. Journal of Applied Physics 11, 9 (09 1940), 582–592. https://doi. org/10.1063/1.1712836 arXiv:https://pubs.aip.org/aip/jap/articlepdf/11/9/582/18304743/582_1_online.pdf JORGE J. MORE. 1993. Generalizations of the trust region problem. Optimization Methods and Software 2, 3-4 (1993), 189–209. https://doi.org/10.1080/10556789308805542 arXiv:https://doi.org/10.1080/10556789308805542 Jorge J. Moré and D. C. Sorensen. 1983. Computing a Trust Region Step. SIAM J. Sci. Statist. Comput. 4, 3 (1983), 553–572. https://doi.org/10.1137/0904038 arXiv:https://doi.org/10.1137/0904038 Matthias Müller, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. 2007. Position based dynamics. J. Vis. Comun. Image Represent. 18, 2 (April 2007), 109–118. https: //doi.org/10.1016/j.jvcir.2007.01.005 Rahul Narain, Matthew Overby, and George E. Brown. 2016. ADMM ⊇ projective dynamics: fast simulation of general constitutive models. In Proceedings of the ACM SIGGRAPH/Eurographics Symposium on Computer Animation (Zurich, Switzerland) (SCA ’16). Eurographics Association, Goslar, DEU, 21–28. Maxim Naumov, M. Arsaev, Patrice Castonguay, J. Cohen, J. Demouth, Joe Eaton, S. Layton, N. Markovskiy, I.Z. Reguly, N. Sakharnykh, V. Sellappan, and R. Strzodka. 2015. AmgX: A Library for GPU Accelerated Algebraic Multigrid and Preconditioned Iterative Methods. SIAM Journal on Scientific Computing 37 (10 2015), S602–S626. https://doi.org/10.1137/140980260 Elvis Nava, John Z Zhang, Mike Yan Michelis, Tao Du, Pingchuan Ma, Benjamin F. Grewe, Wojciech Matusik, and Robert Kevin Katzschmann. 2022. Fast Aquatic Swimmer Optimization with Differentiable Projective Dynamics and Neural Network Hydrodynamic Models. In Proceedings of the 39th International Conference on Machine Learning (Proceedings of Machine Learning Research, Vol. 162), Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (Eds.). PMLR, 16413–16427. https://proceedings.mlr.press/v162/nava22a.html R. A. Nicolaides. 1987. Deflation of Conjugate Gradients with Applications to Boundary Value Problems. SIAM J. Numer. Anal. 24, 2 (1987), 355–365. http://www.jstor.org/ stable/2157562 Jorge Nocedal and Stephen J Wright. 2006. Numerical optimization. Springer. Jan Obrist, Miguel Zamora, Hehui Zheng, Ronan Hinchet, Firat Ozdemir, Juan Zarate, Robert K. Katzschmann, and Stelian Coros. 2025. PokeFlex: A Real-World Dataset of Volumetirc Deformable Objects for Robotics. Under review (2025). R.W. Ogden. 1984. Non-linear Elastic Deformations. E. Horwood. https://books.google. com.hk/books?id=-pEeAQAAIAAJ Matthew Overby, George E. Brown, Jie Li, and Rahul Narain. 2017. ADMM ⊇ Projective Dynamics: Fast Simulation of Hyperelastic Models with Dynamic Constraints. IEEE

Lai, et al.

Transactions on Visualization and Computer Graphics 23, 10 (Oct 2017), 2222–2234. https://doi.org/10.1109/TVCG.2017.2730875 Santiago Paternain, Aryan Mokhtari, and Alejandro Ribeiro. 2019. A Newton-Based Method for Nonconvex Optimization with Fast Evasion of Saddle Points. SIAM Journal on Optimization 29, 1 (2019), 343–368. https://doi.org/10.1137/17M1150116 arXiv:https://doi.org/10.1137/17M1150116 Yue Peng, Bailin Deng, Juyong Zhang, Fanyu Geng, Wenjie Qin, and Ligang Liu. 2018. Anderson acceleration for geometry optimization and physics simulation. ACM Trans. Graph. 37, 4, Article 42 (July 2018), 14 pages. https://doi.org/10.1145/3197517. 3201290 Ting Kei Pong and Henry Wolkowicz. 2014. The generalized trust region subproblem. Comput. Optim. Appl. 58, 2 (June 2014), 273–322. https://doi.org/10.1007/s10589013-9635-7 Yi-Ling Qiao, Junbang Liang, Vladlen Koltun, and Ming C. Lin. 2020. Scalable differentiable physics for learning and control. In Proceedings of the 37th International Conference on Machine Learning (ICML’20). JMLR.org, Article 727, 10 pages. J.W.S.B. Rayleigh. 1894. The Theory of Sound. Number Volume One in The Theory of Sound. Macmillan. https://books.google.com.tw/books?id=hd8EAAAAYAAJ R. S. Rivlin. 1947. Torsion of a Rubber Cylinder. Journal of Applied Physics 18, 5 (05 1947), 444–449. https://doi. org/10.1063/1.1697674 arXiv:https://pubs.aip.org/aip/jap/articlepdf/18/5/444/18307337/444_1_online.pdf Alvaro Sanchez-Gonzalez, Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter W. Battaglia. 2020. Learning to simulate complex physics with graph networks. In Proceedings of the 37th International Conference on Machine Learning (ICML’20). JMLR.org, Article 784, 10 pages. H.A. Schwarz. 1870. Ueber einen Grenzübergang durch alternirendes Verfahren. Druck von Zürcher & Furrer. https://books.google.de/books?id=tMcxYAAACAAJ Xing Shen, Runyuan Cai, Mengxiao Bi, and Tangjie Lv. 2024. Preconditioned Nonlinear Conjugate Gradient Method for Real-time Interior-point Hyperelasticity. In ACM SIGGRAPH 2024 Conference Papers (Denver, CO, USA) (SIGGRAPH ’24). Association for Computing Machinery, New York, NY, USA, Article 96, 11 pages. https://doi. org/10.1145/3641519.3657490 Breannan Smith, Fernando De Goes, and Theodore Kim. 2018a. Stable Neo-Hookean Flesh Simulation. ACM Trans. Graph. 37, 2, Article 12 (March 2018), 15 pages. https://doi.org/10.1145/3180491 Breannan Smith, Fernando De Goes, and Theodore Kim. 2018b. Stable neo-hookean flesh simulation. ACM Transactions on Graphics (TOG) 37, 2 (2018), 1–15. D. C. Sorensen. 1982. Newton’s Method with a Model Trust Region Modification. SIAM J. Numer. Anal. 19, 2 (April 1982), 409–426. https://doi.org/10.1137/0719026 Andrew Spielberg, Allan Zhao, Yuanming Hu, Tao Du, Wojciech Matusik, and Daniela Rus. 2019. Learning-In-The-Loop Optimization: End-To-End Control And Co-Design Of Soft Robots Through Learned Deep Latent Representations. In Advances in Neural Information Processing Systems, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Eds.), Vol. 32. Curran Associates, Inc. https://proceedings.neurips.cc/paper_files/paper/2019/file/ 438124b4c06f3a5caffab2c07863b617-Paper.pdf Nicole Spillane and Daniel Rixen. 2012. Automatic spectral coarse spaces for robust FETI and BDD algorithms. (11 2012). Rasmus Tamstorf, Toby Jones, and Stephen F. McCormick. 2015. Smoothed aggregation multigrid for cloth simulation. ACM Trans. Graph. 34, 6, Article 245 (Nov. 2015), 13 pages. https://doi.org/10.1145/2816795.2818081 Min Tang, Zhongyuan Liu, Ruofeng Tong, and Dinesh Manocha. 2018. PSCC: Parallel Self-Collision Culling with Spatial Hashing on GPUs. Proceedings of the ACM on Computer Graphics and Interactive Techniques 1, 1 (2018), 18:1–18. Joseph Teran, Eftychios Sifakis, Geoffrey Irving, and Ronald Fedkiw. 2005. Robust Quasistatic Finite Elements and Flesh Simulation. In Symposium on Computer Animation, D. Terzopoulos, V. Zordan, K. Anjyo, and P. Faloutsos (Eds.). The Eurographics Association. https://doi.org/10.2312/SCA/SCA05/181-190 Maxime Thieffry, Alexandre Kruszewski, Christian Duriez, and Thierry-Marie Guerra. 2019. Control Design for Soft Robots based on Reduced Order Model. IEEE Robotics and Automation Letters 4, 1 (Jan. 2019), 25–32. https://doi.org/10.1109/LRA.2018. 2876734 Ty Trusty, Otman Benchekroun, Eitan Grinspun, Danny M. Kaufman, and David I.W. Levin. 2023. Subspace Mixed Finite Elements for Real-Time Heterogeneous Elastodynamics. In SIGGRAPH Asia 2023 Conference Papers (Sydney, NSW, Australia) (SA ’23). Association for Computing Machinery, New York, NY, USA, Article 112, 10 pages. https://doi.org/10.1145/3610548.3618220 Ty Trusty, Yun (Raymond) Fei, David Levin, and Danny Kaufman. 2024. Trading Spaces: Adaptive Subspace Time Integration for Contacting Elastodynamics. ACM Trans. Graph. 43, 6, Article 227 (Nov. 2024), 16 pages. https://doi.org/10.1145/3687946 Kenji Ueda and Nobuo Yamashita. 2014. A regularized Newton method without line search for unconstrained optimization. Comput. Optim. Appl. 59, 1–2 (Oct. 2014), 321–351. https://doi.org/10.1007/s10589-014-9656-x Homer F. Walker and Peng Ni. 2011. Anderson Acceleration for Fixed-Point Iterations. SIAM J. Numer. Anal. 49, 4 (Aug. 2011), 1715–1735. https://doi.org/10.1137/ 10078356X

Huamin Wang. 2015. A chebyshev semi-iterative approach for accelerating projective and position-based dynamics. ACM Trans. Graph. 34, 6, Article 246 (Nov. 2015), 9 pages. https://doi.org/10.1145/2816795.2818063 Huamin Wang. 2021. GPU-based simulation of cloth wrinkles at submillimeter levels. ACM Trans. Graph. 40, 4, Article 169 (July 2021), 14 pages. https://doi.org/10.1145/ 3450626.3459787 Huamin Wang and Yin Yang. 2016. Descent methods for elastic body simulation on the GPU. ACM Trans. Graph. 35, 6, Article 212 (Dec. 2016), 10 pages. https: //doi.org/10.1145/2980179.2980236 Zhendong Wang, Longhua Wu, Marco Fratarcangeli, Min Tang, and Huamin Wang. 2018. Parallel Multigrid for Nonlinear Cloth Simulation. Computer Graphics Forum 37, 7 (2018), 131–141. https://doi.org/10.1111/cgf.13554 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/cgf.13554 Botao Wu, Zhendong Wang, and Huamin Wang. 2022. A GPU-based multilevel additive schwarz preconditioner for cloth and deformable body simulation. ACM Trans. Graph. 41, 4, Article 63 (July 2022), 14 pages. https://doi.org/10.1145/3528223. 3530085 Longhua Wu, Botao Wu, Yin Yang, and Huamin Wang. 2020. A Safe and Fast Repulsion Method for GPU-based Cloth Self Collisions. ACM Trans. Graph. 40, 1, Article 5 (Dec. 2020), 18 pages. https://doi.org/10.1145/3430025 Xiaofeng Wu, Rajaditya Mukherjee, and Huamin Wang. 2015. A unified approach for subspace simulation of deformable bodies in multiple domains. ACM Trans. Graph. 34, 6, Article 241 (Nov. 2015), 9 pages. https://doi.org/10.1145/2816795.2818065 Zangyueyang Xian, Xin Tong, and Tiantian Liu. 2019. A scalable galerkin multigrid method for real-time simulation of deformable objects. ACM Trans. Graph. 38, 6, Article 162 (Nov. 2019), 13 pages. https://doi.org/10.1145/3355089.3356486 Yin Yang, Weiwei Xu, Xiaohu Guo, Kun Zhou, and Baining Guo. 2013. BoundaryAware Multidomain Subspace Deformation. IEEE Transactions on Visualization and Computer Graphics 19, 10 (2013), 1633–1645. https://doi.org/10.1109/TVCG.2013.12 Behrooz Zarebavani, Danny M. Kaufman, David I. W. Levin, and Maryam Mehri Dehnavi. 2025. Adaptive Algebraic Reuse of Reordering in Cholesky Factorizations with Dynamic Sparsity Patterns. ACM Trans. Graph. 44, 4, Article 119 (July 2025), 17 pages. https://doi.org/10.1145/3731179 Ziqiu Zeng, Siyuan Luo, Fan Shi, and Zhongkai Zhang. 2025. Fast But Accurate: A Real-Time Hyperelastic Simulator with Robust Frictional Contact. ACM Trans. Graph. 44, 4, Article 88 (July 2025), 19 pages. https://doi.org/10.1145/3730834 Xuejun Zhang. 1992. Multilevel Schwarz methods. Numer. Math. 63, 1 (Dec. 1992), 521–539. https://doi.org/10.1007/BF01385873 Yongning Zhu, Eftychios Sifakis, Joseph Teran, and Achi Brandt. 2010. An efficient multigrid method for the simulation of high-resolution elastic solids. ACM Trans. Graph. 29, 2, Article 16 (April 2010), 18 pages. https://doi.org/10.1145/1731047. 1731054

Record · ID 187285 · SHA-256 61c42e926d58fa40
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.