Contact-rich inverse problems require efficient gradients and reliable enforcement of non-penetration. Existing differentiable methods provide fast simulation and rigid-deformable coupling, but contact treatments based on penalties or velocity-level constraints can permit penetration or positional drift. Incremental Potential Contact (IPC) maintains non-penetration through geometric barriers and {continuous collision detection (CCD)}. However, repeated forward and backward evaluations remain costly in full-space finite-element (FEM) implementations, particularly when nearly rigid components retain many vertex {degrees of freedom (DoFs)}. We present DiffGIPC, a GPU-accelerated differentiable IPC framework for optimizing interacting deformable and nearly rigid bodies. Our affine-body derivatives enable {rest-shape} design, initial-state optimization, and articulated control while retaining only twelve motion DoFs per affine body. These gradients account for how changes in geometry, placement, and actuation affect other objects through collision and friction. We preserve these interactions in one coupled FEM-affine adjoint system per frame, reducing the backward-system size without discarding the dependencies that drive the inverse problem. GPU-optimized backward computation further reduces the cost of repeated gradient evaluations. A reusable PyTorch interface connects supported physical parameters to user-defined parameter maps and losses across design, identification, motion, placement, and articulated control. Experiments demonstrate peak speedups over DiffIPC of approximately 150.0x with full-space FEM and 4133.3x with affine bodies for combined forward simulation and backward differentiation.
Figures & tables
Figure 1. Optimizing the placement of seven affine bodies. Keeping their shapes fixed, DiffGIPC adjusts only the horizontal initial positions to form a stable stack. The panels show the unoptimized configuration and the simulated results after 10, 50, and 78 iterations. Four panels show the initial stone placement and optimization updates 10, 50, and 78 at selected rendered states.
Figure 2. DiffGIPC’s backward pipeline. Each reverse frame assembles one coupled FEM-ABD adjoint system, retaining collision, friction, and joint interactions with twelve DoFs per affine body. The GPU solve and local accumulation propagate gradients to earlier states and physical parameters for optimization through PyTorch. Saved FEM and affine states feed local derivative evaluation, matrix assembly, one coupled adjoint solve, and gradient accumulation. The PyTorch loss supplies gradients, and the parameter map returns physical-parameter gradients to the optimizer.
Figure 3. Reverse-time accumulation in the frame-wise adjoint. Dashed arrows indicate forward state dependencies. Blue arrows carry contributions from the solved λt to the earlier right-hand sides bt−1 and bt−2 . For the dynamic predictor shown, inertia and friction contribute to bt−1 , while inertia contributes −Mλt to bt−2 . Each right-hand side also includes its direct loss gradient. Three time levels are stacked vertically in a single-column figure. Saved states occupy the left lane, and the coupled solve and earlier adjoint right-hand sides occupy the right lane. Dashed arrows point toward newer states. Blue arrows carry the current adjoint's contributions to the two earlier frames.
Figure 4. Garment Hanger shape design. Two shared Bézier control scales modify both shoulders of the affine rest surface under frictional cloth contact. The initial and optimized configurations are shown at frame 33. The optimized hanger retains both garment straps. Initial and optimized garment hangers shown side by side with their supported cloth at frame 33.
Figure 5. PCG computation and execution control. (a) Five stages of one iteration (Algorithm 2 ). (b) A host-controlled loop retrieves a convergence value after each GPU iteration. (c) A conditional CUDA Graph keeps iteration and stopping on the GPU. Frame restoration, Hessian assembly, and gradient accumulation remain outside the graph. Both execution modes use the same PCG recurrence and stopping criterion. The top row shows matrix multiplication, solution and residual update, preconditioning, reduction and search-direction update, and the stopping test. Two CPU and GPU lane diagrams compare a host-controlled loop with a device-controlled conditional graph. Only the host-controlled loop returns a convergence scalar to the CPU at every iteration.
Figure 6. Collision-mediated initial-velocity optimization. The first bunny’s initial velocity is adjusted so that the second bunny’s horizontal centroid reaches its target. The panels compare frame 40 before optimization and after 1, 3, and 10 iterations. Four panels show the initial placement and optimization updates 1, 3, and 10 at simulation frame 40.
Figure 7. Trajectory matching through contact between two deformable bodies. The first body’s initial velocity is optimized to reproduce the target motion. Opaque surfaces show the simulation and translucent surfaces the target at frame 100, before optimization and after 1, 3, and 10 iterations. Trajectory matching for two contacting deformable bodies. The shared initial velocity of the first body is optimized to reproduce a prescribed motion. Opaque surfaces show the simulated bodies and translucent surfaces show the target configuration. The panels show frame 100.
Example
System
Forward
Backward
Forward + Backward
Linear
CCD
Hessian
Search
Other
Total
Hessian
Linear
Other
Total
Bunny collision (Fig. 6 )
DiffIPC
4.5
10
4.6
1.1
1.9
22
1.4
1.1
0.44
2.9
25
DiffGIPC
0.65
0.044
0.037
0.017
0.002
0.75
0.0069
0.22
0.0057
0.24
0.98
Speedup
6.9×
232.0×
123.9×
67.7×
963.8×
30.0×
204.7×
4.9×
76.3×
12.5×
25.8×
Tentacle matching (Fig. 7 )
DiffIPC
140
14
120
66
11
350
7.4
11
2.8
21
370
DiffGIPC
8.1
2.4
2.2
1.4
0.011
14
0.04
1
0.024
1.1
15
Table 1. Mean time per forward and backward evaluation (seconds) on the shared FEM benchmarks. Forward + Backward gives the sum of the two per-evaluation means.
Figure 8. Rest-shape optimization of a loaded FEM support. Under fixed constraints and downward loading, the objective minimizes squared first-Piola stress. The panels show the deformed support at frame 20 before optimization and after 20, 50, and 200 iterations. Four front-view stress renderings compare the initial support and the designs after 20, 50, and 200 updates.
Figure 9. Fitting the material response of a constrained cube to marker positions. Each column pairs side (top) and oblique (bottom) views at frame 19, before optimization and after 2, 5, and 20 iterations. Opaque surfaces denote the simulation, translucent surfaces the target, and arrows the Dirichlet constraints. A two-row, four-column image sequence for material identification. Side views form the top row and oblique views the bottom row. Each column pairs the two views of the same optimization iteration, ordered 0, 2, 5, and 20 from left to right. Solid surfaces show the simulated deformation, translucent surfaces show the target, and arrows show the constraint directions.
Figure 10. Friction optimization on a curved slide. A bunny-slide friction coefficient is adjusted so that the bunny reaches its target position at frame 100. The panels show the unoptimized simulation and results after 1, 3, and 5 iterations. A deformable bunny slides along a curved surface. The bunny-slide friction coefficient is optimized so that its position at the specified time matches a target position. The panels show final frame 100 for the initial coefficient and after one, three, and five optimization updates.
Figure 11. Optimization histories for the shared FEM benchmarks. Blue and red denote DiffGIPC and DiffIPC, respectively. The horizontal axes count accepted parameter updates from the initial parameters (iteration 0). Loss scales are task-specific, and Loaded Support is remeshed between updates. Five side-by-side plots show loss against optimization iteration for bunny collision, tentacle matching, loaded-support shape design, material fitting, and friction optimization. Blue curves denote DiffGIPC, and red curves denote DiffIPC.
Figure 12. Cloth rest-scale fitting on a fixed mannequin. The blue panels show the initial drape and results after 1, 3, and 5 updates. The orange panel shows the target. Five garment-focused views with one common crop defined from the target. The four blue panels show the initial, intermediate, and optimized drapes. The orange panel shows the target. Sleeves and hems remain visible so changes in garment size can be compared at a common scale. Orange identifies the target, not stress.
Figure 13. Grasp-angle control. Initial (left) and optimized (right) grasps at frame 60. Joint-target optimization moves the grasp toward the 0.45 -rad target. Coral guides follow the fingers. Dashed blue guides indicate visual reference directions. Initial grasp beside optimized grasp, with finger and reference guides.
Figure 14. Articulated cup catching through plate-trajectory optimization. With the arm held at rest (left), the cup falls from the plate. After optimization (right), the plate catches and retains it. The middle column shows the initial IK trajectory, which also fails to retain the cup. Four rows show the same saved frames (94, 111, 120, and 239). Rest-arm controls are on the left, the failing initial IK trajectory is in the middle, and optimized controls are on the right.
Figure 15. Dual-arm towel grasp by optimizing one vertical offset applied to both grippers’ target positions. Before optimization (left), the arms fail to grasp the towel. After optimization (right), they grasp the towel and then lift and fold it. Rows show the motion at 0 , 2.00 , 2.40 , and 3.30 s. Four rows pair the same saved frames of the dual-arm towel grasp. Before optimization is on the left, and after optimization is on the right. The optimized arms grasp, lift, and fold the towel.
Figure 16. Optimization histories for the additional tasks. Stone placement changes its centroid tolerance after 50 updates. Garment Hanger, cup catching, and Dual-arm grasp use constant offsets that set their final plotted values to zero. Six loss plots for stone placement, cloth scale, Garment Hanger shape design, grasp angle, cup catching, and dual-arm grasp.
Example
System
Dim
Relative FD error
Stone placement (Fig. 1 )
DiffGIPC
14
7.53×10−11
Garment Hanger (Fig. 4 )
DiffGIPC
2
5.45×10−3
Bunny collision (Fig. 6 )
DiffGIPC
3
1.21×10−6
DiffIPC
3
1.01×10−4
Tentacle matching (Fig. 7 )
DiffGIPC
3
3.33×10−5
DiffIPC
3
3.83×10−11
Table 2. Finite-difference checks at the initial optimization parameters.
Resolution / case
System
Forward
Backward
Forward + Backward
Linear
CCD
Hessian
Search
Other
Total
Hessian
Linear
Other
Total
Bunny collision (Fig. 6 ): resolution study
933 vertices
DiffGIPC ABD+ABD
0.025
0.039
0.041
0.021
0.0023
0.13
0.0081
0.0042
0.0032
0.016
0.14
DiffGIPC ABD+FEM
0.58
0.049
0.058
0.03
0.0025
0.72
0.011
0.22
0.0085
0.24
0.96
DiffGIPC FEM+FEM
0.65
0.044
0.037
0.017
0.002
0.75
0.0069
0.22
0.0057
0.24
0.98
DiffIPC FEM+FEM
4.5
10
4.6
1.1
1.9
22
1.4
1.1
0.44
2.9
25
Table 3. Mean forward and backward evaluation times (seconds) for the bunny-collision resolution study and Garment Hanger rest-shape design. Forward + Backward gives the sum of the two per-evaluation means.
Example
Nfwd
Nbwd
Forward mean
Backward mean
Teval
Linear
CCD
Hessian
Search
Other
Total
Hessian
Linear
Other
Total
PCG
Others
PCG
Others
Stone placement (Fig. 1 )
80
78
0.097
0.0093
0.079
0.12
0.024
0.0088
0.34
0.051
0.067
0.003
0.023
0.14
39
Bunny collision (Fig. 6 )
15
12
0.62
0.027
0.044
0.037
0.017
0.002
0.75
0.0069
0.22
0.005
0.0057
0.24
14
Tentacle matching (Fig. 7 )
11
10
6.8
1.3
2.4
2.2
1.4
0.011
14
0.04
0.99
0.024
0.024
1.1
170
Garment Hanger (Fig. 4 )
32
30
2.8
0.87
0.88
1.5
0.8
0.0088
6.9
0.033
0.2
0.01
0.022
0.27
230
Table 4. Complete optimization workloads and mean phase times (seconds) of the eleven examples. Nfwd counts forward evaluations, including line-search trials and verification runs. Nbwd counts complete gradient evaluations used for optimization. Linear comprises PCG and Others, and Hessian combines local and global assembly.
Example
Graph
Forward (s)
Backward (s)
Total (s)
Linear
Other
Total
Linear
Other
Total
PCG
Others
PCG
Others
Stone placement (Fig. 1 )
On
0.097
0.0093
0.24
0.34
0.067
0.003
0.074
0.14
0.49
Off
0.18
(1.8×)
0.0094
0.24
0.42
(1.2×)
0.12
(1.8×)
0.003
0.073
0.2
(1.4×)
0.62
(1.3×)
Bunny collision (Fig. 6 )
On
0.62
0.027
0.1
0.75
0.22
0.005
0.013
0.24
0.98
Off
0.97
(1.6×)
0.026
0.1
1.1
(1.5×)
0.35
(1.6×)
0.0051
0.011
0.36
(1.5×)
1.5
(1.5×)
Table 5. CUDA Graph enabled (On) versus disabled (Off). Times are mean seconds per evaluation; parentheses in the Off rows give speedups from enabling the graph. Linear comprises PCG and Others; the separate Other is Total minus Linear. The rightmost Total combines forward and backward time.
Appendix figures & tables2 assets
Supplementary material from the paper’s appendix.
Appendix
Example
h
Horizon
drel
κ
df
#DOF
Newton ϵN
Friction
E
κA
Stone placement (Fig. 1 )
0.03
1000
10−3
105
10−3
84
10−3
0.3
–
107
Bunny collision (Fig. 6 )
0.05
40
10−3
1010
–
5598
3×10−4
0
4×107
–
Tentacle matching (Fig. 7 )
0.005
100
10−3
107
10−3
20688
10−2
1
2×104
–
Garment Hanger (Fig. 4 )
1/120
600
10−3
105
10−3
51357
10−2
0.05
5×104
108
Loaded Support (Fig. 8 )
0.02
20
3×10−4
106
–
22545
10−3
0
108
–
Material fitting (Fig. 9 )
0.1
10
10−3
107
–
12567
10−2
0
106
–
Appendix
Table 6. Simulation settings for the design, identification, and control tasks in Section 7 . The columns report time step h , relative collision activation distance drel , barrier coefficient κ , friction-smoothing velocity threshold df , motion DoFs, relative Newton tolerance ϵN , friction coefficient, Young’s modulus E for volumetric FEM or cloth stretching, and affine stiffness κA . For Dual-arm grasp, 0.3 is the arm–cloth friction coefficient and 0.2 is the cloth–table coefficient.
Physical group
Differentiated quantities
Role
Tetrahedral material
μ , λ
Stable Neo-Hookean material response
Initial FEM motion
Selected initial vertex velocities
Initial dynamic state
Deformable rest shape
Selected rest coordinates
Volumetric and cloth rest geometry
Cloth material
Stretch, shear, and bending coefficients
Cloth constitutive response
Friction
Mesh self and ground friction coefficients
Frictional contact
Initial affine motion
q0,b , q˙0,b
Initial affine state
Appendix
Table 7. Physical parameter families exposed to task-specific design variables through the interface in Section 6.3 . Each family contributes local derivatives to the shared coupled adjoint.
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) trust-region 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 100x 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.
Shih-Yu Lai, Sung-Han Tien, Jui-I Huang +9
National Taiwan University · MoonShine Animation Studio · National University of Singapore +3
Generating intelligent robot behavior in contact-rich settings is a research problem where zeroth-order methods currently prevail. Developing methods that make use of first/second order information about rigid-body dynamics in the presence of contact holds great promise in terms of increasing the solution speed and computational efficiency. The main bottleneck in this research direction is the difficulty in obtaining gradients and Hessians that are actually useful for numerical optimization, due to pathologies in all three steps of a common simulation pipeline: i) collision detection, ii) contact dynamics, iii) time integration. This abstract proposes a method that aims to address the collision detection part of the puzzle, via a novel pipeline designed from scratch with smooth (i.e. twice) differentiability and massive vectorizability on GPUs as the main priorities. This is in contrast to standard collision detection routines that are instead optimized for runtime on CPUs and minimal memory footprint, but do employ logic and control flow that hinder differentiability and vectorization. The proposed pipeline consists of the following contributions: i) highly expressive and compute efficient SDF representations, ii) differentiable broad-phase and narrow-phase routines that use these representations to generate vertex-SDF and edge-SDF contacts, iii) a differentiable routine for convex decomposition based contact blending.
We present IsaacIPC, a robotic simulation framework that couples GPU accelerated incremental potential contact (IPC) with IsaacSim/Lab. IsaacIPC maps simulated deformation between simulation and visual meshes, enabling real-time realistic rendering with applications to data collection and policy evaluation. For tactile sensing, we introduce the geometric mortar contact potential (GMCP), which defines a barrier potential over contact samples on tactile surfaces to better resolve contact-pressure distributions. We evaluate GMCP on contact benchmarks and demonstrate IsaacIPC on rigid-deformable robotic simulations including a quadruped robot, a dexterous hand, and a universal manipulation interface (UMI) gripper.
Qixin Liang, Zhongqing Han
Anker Humanoid Lab · Shenzhen, China · The University of Hong Kong +1