Mastering Multi-State Dynamics: Continuum Damage and Viscoplastic MPM for RL-Based Dexterous Manipulation

Introduction

Dexterous manipulation of multi-state materials includes cutting fruit, scooping ice cream, shoveling sand, and splitting brittle solids. These tasks require a simulator that handles large deformation, topology change, and two-way contact between tools and materials while remaining fast enough for online reinforcement learning. The Material Point Method (MPM) has become the natural particle-grid substrate for graphics-scale soft, brittle, granular, and viscoplastic materials [SSC*13,SSJ*14,KGP*16]. CD-MPM-style approaches add a damage variable so that fracture is driven by stress rather than by mesh topology updates [33 WCL*20]. In parallel, robotic cutting benchmarks have begun to target heterogeneous materials with force feedback [38, 40 KSC*26, 19].

Two practical obstacles remain when this line of work is plugged into an RL loop. First, MPM cost grows with the number of active particles, but direct tool contact is confined to a small region of the material: a uniform high-resolution grid spends the rollout budget on regions the tool never reaches. UBSoft addresses this for elastoplastic media by adapting MPM resolution to grid proximity, but does not model damage or cutting [LFW*24]. Second, CD-MPM resolves fracture physically through a phase-field damage field [33], but in its current form runs on a single global grid, and is too costly to drop directly into PPO-style training that demands millions of frames.

We bridge this gap with a CD-MPM-style cutting simulator wrapped as a Gym-like RL environment. The simulator drives separation through a tool-aware damage layer, so the policy decides where cracks form rather than following a kinematic blade path over a static crack model. MGPCG solves the damage smoothing equation within a tool-centered active region, making the implicit phase-field damage solve cost scale with the patch volume rather than the full grid. MPM contact impulses are converted into forces

@page 1/16 starts_mid=0 ends_mid=1

Introduction

Dexterous manipulation of multi-state materials includes cutting fruit, scooping ice cream, shoveling sand, and splitting brittle solids. These tasks require a simulator that handles large deformation, topology change, and two-way contact between tools and materials while remaining fast enough for online reinforcement learning. The Material Point Method (MPM) has become the natural particle-grid substrate for graphics-scale soft, brittle, granular, and viscoplastic materials [SSC*13,SSJ*14,KGP*16]. CD-MPM-style approaches add a damage variable so that fracture is driven by stress rather than by mesh topology updates [33 WCL*20]. In parallel, robotic cutting benchmarks have begun to target heterogeneous materials with force feedback [38, 40 KSC*26, 19].

Two practical obstacles remain when this line of work is plugged into an RL loop. First, MPM cost grows with the number of active particles, but direct tool contact is confined to a small region of the material: a uniform high-resolution grid spends the rollout budget on regions the tool never reaches. UBSoft addresses this for elastoplastic media by adapting MPM resolution to grid proximity, but does not model damage or cutting [LFW*24]. Second, CD-MPM resolves fracture physically through a phase-field damage field [33], but in its current form runs on a single global grid, and is too costly to drop directly into PPO-style training that demands millions of frames.

We bridge this gap with a CD-MPM-style cutting simulator wrapped as a Gym-like RL environment. The simulator drives separation through a tool-aware damage layer, so the policy decides where cracks form rather than following a kinematic blade path over a static crack model. MGPCG solves the damage smoothing equation within a tool-centered active region, making the implicit phase-field damage solve cost scale with the patch volume rather than the full grid. MPM contact impulses are converted into forces and torques on the rigid tool, allowing material resistance to feed back into policy decisions.

Concretely, our contributions are:

• A unified GPU-accelerated MPM environment for soft solids, brittle materials, granular media, and viscoplastic or semi-fluid media, with a shared constitutive layer and an interface for reinforcement learning.

• A coupled phase-field damage formulation with a local solve around the tool, accelerated by conjugate gradients with geometric multigrid preconditioning. The background MPM grid remains uniform.

• Bidirectional coupling between material and tool dynamics through contact impulses, with tools held by friction in a gripper. Cutting policies are trained with PPO through rigid-body pretraining followed by fine-tuning in the MPM environment.

Related Work

MPM Simulation for Robot-Oriented Multi-State Materials

The Material Point Method (MPM) was introduced as a particle-in-cell extension for history-dependent solid mechanics [21, 29] and has become a standard graphics tool for large deformation, contact, and topology change across snow, sand, phase-change materials, and viscoplastic media [SSC∗13, SSJ∗14, KGP∗16, JST∗16]. Two threads in this body of work matter most for tool-driven cutting. MPM fracture methods keep a particle-grid formulation but drive separation through damage rather than through topology updates. CD-MPM evolves a Ginzburg-Landau phase-field damage variable. AnisoMPM replaces the phase field with non-local geometric damage and structural tensors for anisotropic materials. Ductile fracture has been formulated through elastoplasticity with Rankine or von Mises damage and a particle-cloud fracture-surface visualization. A recent system pairs an optimized CD-MPM with 3D Gaussian splatting for mixed-material brittle fracture [WFL∗19, WCL∗20, WDG∗19, HCL∗26]. A parallel thread closes the loop with the tool: MLS-MPM with CPIC supports cutting and two-way rigid-body coupling directly by allowing particle-grid transfers across displacement discontinuities, building on a longer history of coupling between rigid and deformable bodies in physics simulation [7 SSF08].

Robot-oriented MPM simulators bring this physics into specific manipulation regimes: tactile contact, fluid manipulation, soft-robot co-design, deformable-object reconstruction, surgical or food cutting, and multiphysics RL stability [XZX∗23, WMS∗23, SZB∗24, CHS∗26, KSC∗26, OT25, XLO25], and engine-level systems such as Genesis target general embodied-AI physics [5]. Closest to our scaling concern, UBSoft drops a multi-level grid into MLS-MPM and adapts cell size to robot proximity, but is restricted to elastoplastic media (sand, snow) and does not model damage or cutting [14]. Across this line, no prior system combines a CD-MPM-style damage layer for tool-driven separation with a tool-centered localized damage solve inside a single uniform-grid MLS-MPM solver. Our simulator provides this combination and exposes it through a Gym-like RL interface while preserving two-way rigid-tool coupling.

Cutting Robotics: Simulators, Real Robots, and Policy Learning

Robotic cutting combines a separation model on the simulator side with control and feedback on the policy side. Mesh-based simulators rely on topology updates, virtual nodes, or remeshing on tetrahedral and polyhedral discretizations [18 WVD15]. Among recent differentiable cutting simulators, DiSECt augments FEM with a continuous-contact SDF and a damage-spring cohesive zone along a prescribed cutting plane and makes the whole pipeline differentiable [9]. The MPM side has grown in parallel: RoboNinja learns a tolerance-modulated cutting policy for soft objects with rigid cores on top of MPM dynamics with sparse force-collision feedback [XXL∗23]. CRESSim-MPM provides a GPU MPM library with contact between soft and rigid bodies tailored to surgical cutting and suturing [19]. CulinaryCut-VLAP builds an MLS-MPM food-cutting environment in which contact forces are estimated from particle-grid impulse exchange [13]. TopoCut extends this line to multi-step cutting, pairing a particle-based elastoplastic solver with a damage-driven topology-discovery mechanism and a discrete diffusion policy guided by spectral rewards [30].

On the policy side, model-free reinforcement learning, particularly PPO, remains the practical backbone for continuous robotic control [25]. Broader deformable-object benchmarks have extended learning from cloth and fluids to plasticine, dough-like materials, dexterous deformables, and graph-based shaping [LOWH21, HHD∗21, LHC∗23, CXY∗24, SXH∗22, SXC∗23]. Real-robot cutting work emphasizes force feedback and multi-material structure: a recent chicken-deboning system pairs a custom multi-material cutting simulator with a hardware testbed and uses residual reinforcement learning on discretized force observations [40], and the simulator side cited above is consistently motivated by hardware-relevant force-aware behavior [XXL∗23, KSC∗26, OT25]. Across this line, simulators are typically tuned for fracture fidelity in a fixed-resolution scene or for a single application domain, while real-robot pipelines often replay kinematic tool trajectories or operate without a fast simulator behind them. Our framework closes this gap by exposing a CD-MPM-style multi-material cutting model and a tool-centered localized damage solve behind a Gym-like interface with two-way rigid-tool feedback, so the PPO policies of § 5 train on the same contact physics the tool will feel at deployment.

Method

We call our method MG-CD-MPM. It extends the CD-MPM update of [WFL∗19] in two ways for RL-driven tool manipulation. First, we solve the stress-driven phase field only in a tool-centered region $\mathcal{N}(T)$, as described in § 4.3. Second, we precondition PCG with a 3D geometric-multigrid V-cycle, yielding MGPCG. Crack initiation and propagation remain stress-driven; only the solve region is localized. The same multigrid approach has recently been applied to real-time graphics simulation in a mixed-FEM elasticity setting [20], complementary to the MPM substrate we target here. The MPM background grid itself stays uniform at cell size $h_0$ everywhere. Per timestep the simulator performs the standard three-stage MPM cycle of P2G, grid update, and G2P, with the localized MGPCG damage solve sitting inside Stage 2; the full step is summarized in Alg. 1. Figure 2 visualizes the corresponding momentum and phase-field streams and the tool–policy feedback loop.

Coupled Equations of Motion

Following CD-MPM [WFL∗19], we couple momentum to a phase-field damage variable $d \in[0,1]$, with $d=0$ an undamaged material point and $d=1$ a fully cracked one. We use the damage convention $d=1-c$, where $c$ is the intact phase field used in CD-MPM.

Energy split. A brittle solid should not crack under compression alone, so the elastic energy is decomposed into a tensile part $\psi^{+}$ that drives damage and a compressive part $\psi^{-}$ that does not:

with $\mathbf{P}^{\pm}=\partial\psi^{\pm}/\partial\mathbf{F}$. The factor $(1-d)^{2}$ acts only on the tensile stress $\mathbf{P}^{+}$. The compressive response $\mathbf{P}^{-}$ remains intact even at $d=1$. A fully cracked region therefore resists compression but has no tensile or shear resistance. For the brittle material class (§ 4.4) we adopt the volumetric–deviatoric split of Amor et al. [1],

with bulk modulus $K$, shear modulus $\mu$, infinitesimal strain $\boldsymbol{\epsilon}=\frac{1}{2}(\mathbf{F}+\mathbf{F}^{\mathsf{T}})-\mathbf{I}$, and its deviator $\boldsymbol{\epsilon}_{\mathrm{dev}}$. The ramp functions $\langle x\rangle_{\pm}=(x\pm|x|)/2$ separate the trace of strain into volumetric expansion, which feeds $\psi^{+}$, and volumetric compression, which feeds $\psi^{-}$. The deviator enters $\psi^{+}$ in full, since shear distortion produces shear-driven cracking regardless of the sign of the volumetric part.

Evolution equations. The momentum and damage fields then satisfy

Equation (4) is the standard balance of linear momentum with body forces $\mathbf{b}$; damage enters the dynamics only by softening $\mathbf{P}^{+}$ wherever $d>0$. Equation (5) is the Ginzburg–Landau phase-field equation. The drive $2(1-d)\mathcal{H}$ increases $d$ when the past tensile energy $\mathcal{H}$ is high, and vanishes as $d\to1$. The restoring term $(G_{c}/\ell_{c})d$ opposes growth and balances the drive near $\mathcal{H}\sim G_{c}/(2\ell_{c})$, which sets the crack-initiation threshold. The Laplacian term $G_{c}\ell_{c}\Delta d$ spreads damage over a band of width $\sim\ell_{c}$ instead of a zero-width discontinuity. The viscosity $\eta$ stabilizes the update; we keep it small so that $d$ follows the stress state within a substep.

Driving energy and irreversibility. The driver $\mathcal{H}$ is the maximum past tensile elastic energy density at each particle,

The running maximum of $\psi^{+}$ retains the loading history and continues to drive $d$ after the material relaxes. The damage field itself is held to the same constraint by an explicit projection,

applied in Stage 3, which ensures that even numerical oscillations of Eq. (5) cannot reverse an opened crack.

Tool-local activation. We restrict Eq. (5) to a neighborhood of the rigid tool,

where $\phi_{\mathbf{T}}$ is the tool SDF and $r_{d}$ a small activation radius (a few finest cells $h_{0}$ around the blade). Outside $\mathcal{N}(\mathbf{T})$, we set $\dot{d}=0$ and freeze $\mathcal{H}$. The particle state is then preserved: a particle that leaves with $d_{p}=0.6$ keeps that value until the tool returns, while an untouched particle remains at $d_{p}=0$.

The restriction in Eq. (8) controls only when new damage can grow; it does not gate force exchange. The momentum update, damaged stress, particle–grid transfers, and rigid-tool projection still run over the full domain. Far-field material can therefore push on the tool, and previously damaged material transmits the stress set by its stored damage. Section 3.5 analyzes the resulting wrench.

Three-Stage MPM Update with Multigrid Damage Solve

A timestep $t^{n}\to t^{n+1}$ advances particle and grid state in the canonical P2G / grid-update / G2P sequence on a uniform background grid of cell size $h_{0}$. The momentum part of the step is standard MLS-MPM [7]; our damage solve runs in Stage 2, where the implicit phase-field damage system is solved by MGPCG on a tool-centered patch (§ 4.3).

Stage 1: Particles to grid. Particle $p$ contributes to nodes $i$ in its quadratic-B-spline neighborhood with weights $w_{ip}^{n}=N(\mathbf{x}_{i}-\mathbf{x}_{p}^{n})$. Following MLS-MPM [7], mass and momentum are transferred with the per-particle affine matrix $\mathbf{C}_{p}^{n}$,

and for brittle particles in $\mathcal{N}(\mathbf{T})$ damage is transferred mass-weighted,

Stage 2: Grid update. Grid momentum is integrated at every node $i$ under the internal force from the damaged-split stress of Eq. (1) together with external forces $\mathbf{f}_{i}^{\mathrm{ext},n}$,

Eqs. (12)–(13) run over the entire grid: every particle, whether inside or outside $\mathcal{N}(\mathbf{T})$, contributes its current damaged-split stress to the internal force, and every grid node receives a momentum update. This is what makes far-field elastic material capable of pushing back on the tool, and what lets a previously-damaged region transmit a stress reduced by its frozen $d_p$.

Damage itself is advanced only on the active nodes $i \in\mathcal{N}(\mathcal{T})$. Treating $\mathcal{H}_i^n$ and $d_i^n$ in the source term explicitly while keeping $d_i^{n+1}$ and the Laplacian implicit, Eq. (5) discretizes as

Over the active nodes $\{i \in\mathcal{N}(\mathcal{T})\}$ this is a sparse, symmetric positive-definite linear system

with $L_\ell$ the standard 7-point (3D) discrete Laplacian on the level-$\ell$ stencil of an internal multigrid hierarchy and homogeneous-Neumann conditions at $\partial\mathcal{N}(\mathcal{T})$. The original CD-MPM paper solves this implicit phase-field damage system by Jacobi-preconditioned conjugate gradient on a single fine grid [33].

We solve Eq. (15) with PCG and a geometric multigrid V-cycle preconditioner [2]. Each cycle smooths the fine grid, restricts the residual to coarser grids by $2^d$ full-weighting, recurses to the coarsest level, and prolongs the correction back. We then post-smooth on the way up. The production MGPCG solver stops when $\lVert r\rVert_2/\lVert b\rVert_2 < 10^{-5}$ or after 20 outer iterations. Section 4.4.1 compares it with localized Jacobi-PCG and the other solver baselines. The use of multigrid as an inner solver for MPM-side implicit systems is itself well-established for implicit time-integration [34]; we transplant it to the damage-smoothing problem under the single-resolution-grid setting, restricted to the tool-tracking active region of Sec. 3.3. Outside $\mathcal{N}(\mathcal{T})$ we hold $d_i^{n+1}=d_i^n$. This is controlled localization, not an adaptive MPM discretization. Its accuracy depends on keeping the patch boundary far enough from the evolving damage band (Sec. 4.4.1).

Independently of the damage solve, rigid-tool boundary conditions on $v_i^{n+1}$ are enforced through $\phi_T$ at every grid node $i \in\partial T$ for which the SDF predicts penetration. The per-node impulse absorbed at the boundary is given by

Here $\mathbf{v}_i^{n+1,\mathrm{free}}$ is the free velocity from Eq. (13) before projection. The impulse is the momentum the tool extracts from the material at that node. It feeds the wrench accumulation of § 4.5, providing the policy with its force-feedback signal whether or not the contacting node lies inside $\mathcal{N}(\mathrm{T})$.

Stage 3: Grid to particles. Particles gather velocity, affine state, and deformation gradient. Particles satisfying $p \in\mathcal{N}(\mathrm{T})$ also gather damage:

We then enforce the irreversibility constraint of Eq. (7) as a per-particle projection,

and refresh the history variable $\mathcal{H}_p$ via Eq. (6) using the updated deformation gradient $\mathbf{F}_p^{n+1}$. When $d_p^{n+1}$ crosses the separation threshold $d^*$, particle $p$ is tagged separated, and subsequent P2G transfers across its cluster boundary use the CPIC discontinuity machinery of [HFG*18].

Tool-Tracking Active Solve Region

We exploit a structural property of tool-driven cutting: at any given moment the damage field is non-trivial only in a small neighborhood of the tool tip; far from the tool the material is intact and the damage scalar stays at its equilibrium value $c \equiv1$ (equivalently $d \equiv0$). In such an intact region, the implicit phase-field system of Eq. (15) has a zero damage state and no active source, so a full-grid V-cycle would spend most of its work in cells where it has nothing to update.

We therefore restrict the damage solve to an axis-aligned patch centered on the tool,

with $(i_{\mathrm{T}},j_{\mathrm{T}},k_{\mathrm{T}})=\left\lfloor(\mathbf{x}_{\mathrm{T}}-\mathbf{x}_0)/h_0\right\rfloor$ the integer cell index of the tool tip, and $(h_x,h_y,h_z)$ half-extents chosen so that the regularized damage band of width $\sim\ell_c$ never reaches $\partial\mathcal{N}(\mathrm{T})$. In practice, the damage field is smoothed every $k$ MPM substeps. We (i) copy the patch $\mathbf{b}\leftarrow d_{\mathrm{grid}}[\mathcal{N}(\mathrm{T})]$, (ii) run residual-controlled MGPCG on the patch with homogeneous Neumann boundary conditions, and (iii) write back $d_{\mathrm{grid}}[\mathcal{N}(\mathrm{T})]\leftarrow$ MGPCG output. Outside $\mathcal{N}(\mathrm{T})$, we keep the historical field unchanged. Nodes that enter the patch later are initialized from this stored field. Localization is accurate when the patch boundary stays far from the evolving damage band; § 5.4.1 measures the resulting error as the patch size varies. The damage-smoothing cost becomes proportional to $|\mathcal{N}(\mathrm{T})|$ rather than to the full grid, an order-of-magnitude saving when the tool occupies a small fraction of the scene, which is the operating regime of every task in § 5.

Figure 3 illustrates how the active damage-solve region follows the tool while the background MPM discretization remains uniform.

Solver localization. Eq. (23) adapts the solver, not the discretization. The MPM background grid stays uniform at cell size $h_0$ everywhere; P2G and G2P transfers of § 4.2 are unchanged. We do not subdivide cells near the tool, do not maintain a sparse hierarchical grid in the sense of UBSoft [LFW*24], and do not perform particle split / merge across resolution levels. The “hierarchy” we exploit is purely internal to the V-cycle preconditioner used by PCG: each preconditioner application recursively visits a ladder of progressively coarser copies of $\mathcal{N}(\mathrm{T})$ at cell sizes $h_\ell=2^\ell h_0$, $\ell=0,\ldots,L$, applying a Jacobi smoother on each level and pro-

Algorithm 1 One simulation step at time $t^n \to t^{n+1}$
Stage 1: Particles to grid
1: for each particle $p$ scattering to nodes $i$ on the background grid do
2: scatter mass and affine momentum
	(Eqs. 9, 10)
3: if $p \in\mathcal{N}(\mathrm{T})$ then
4: scatter damage to $d_i^n$
	(Eq. 11)
5: end if
6: end for
Stage 2: Grid update
7: for each grid node $i$ do
8: assemble $\mathbf{f}_i^{\mathrm{int},n}$ from split damaged stress
	(Eq. 12)
9: update $\mathbf{v}_i^{n+1}$
	(Eq. 13)
10: enforce rigid-tool SDF boundary, record impulse $\mathbf{j}_i$
11: end for
12: solve sparse SPD system $\mathbf{A}_\ell d^{n+1}=\mathbf{b}^n$ on $\mathcal{N}(\mathrm{T})$ via MGPCG
	(Eq. 15)
Stage 3: Grid to particles, then resize
13: for each particle $p$ do
14: gather $\mathbf{v}_p^{n+1},\mathbf{C}_p^{n+1},\mathbf{F}_p^{n+1},\mathbf{x}_p^{n+1}$
	(Eqs. 17–21)
15: if $p \in\mathcal{N}(\mathrm{T})$ then
16: gather $d_p^{n+1}$ from grid
	(Eq. 20)
17: irreversibility: $d_p^{n+1}\leftarrow\max(d_p^{n+1},d_p^n)$
	(Eq. 22)
18: update history $\mathcal{H}_p^{n+1}=\max(\mathcal{H}_p^n,\psi^+(\mathbf{F}_p^{n+1}))$
	(Eq. 6)
19: tag separated if $d_p^{n+1}>d^*$
20: end if
21: end for
Two-way tool integration
22: accumulate $(\mathbf{F}_{\mathrm{T}},\boldsymbol{\tau}_{\mathrm{T}})$ from $\{\mathbf{j}_i\}$
	(Eq. 24)
23: integrate tool 6-DoF state under wrench plus policy action $\mathbf{a}^n$

longing the correction back to the finest patch. Long-wavelength damage error that Jacobi does not efficiently remove on the finest grid becomes short-wavelength error on a coarser one where Jacobi removes it in a handful of iterations.

Material Constitutive Layer

A per-particle tag $m_p \in\{\mathrm{S},\mathrm{B},\mathrm{G},\mathrm{V}\}$ selects how $\mathbf{P}(\mathbf{F}_p)$ in Eq. (4) and the optional return mapping are evaluated:

• Soft elastic-plastic solids (S): fixed-corotated elasticity with an optional von Mises return mapping; damage is disabled ($d_p = 0$).

• Brittle solids (B): linear elasticity with the volumetric–deviatoric energy split of Eqs. (2)–(3), coupled to the phase-field damage of Eq. (5).

• Granular media (G): Drucker–Prager elastoplasticity in the formulation of [KGP* 16]; $d_p = 0$.

• Viscoplastic / semi-fluid materials (V): Herschel–Bulkley with apparent viscosity $\eta_{\mathrm{HB}}(\dot{\gamma})$; $d_p = 0$.

The tag dispatches to the right routine in a single pass over particles, so a single rollout can mix all four classes.

Two-Way Rigid-Tool Coupling

The rigid tool $\mathrm{T}$ is a 6-DoF body with mass $M$ and inertia tensor $\mathbf{I}$. Its SDF $\phi_{\mathrm{T}}$ is attached to the tool frame. Stage 2 records the boundary impulses $\mathbf{j}_i$ (Eq. (16)). We combine them with the node positions $\mathbf{x}_i$ to form a wrench at the tool’s center of mass:

The sums in Eq. (24) include every SDF-boundary node, whether or not it lies in $\mathcal{N}(\mathrm{T})$. The wrench is therefore the full mechanical contact, not a damage-gated signal. We integrate the tool state under this wrench and the policy action. Its motion thus reflects material resistance rather than a replayed kinematic trajectory.

What damage does to the wrench. Damage affects the wrench through the softened stress $\mathbf{P}_d$, which determines the impulses $\mathbf{j}_i$ in Stage 2 (Eq. (12)). The policy sees three phases during a cut. (i) Approach. Near intact material ($d_p \approx0$), the tool feels the full stiffness $\mathbf{P}^{+}+\mathbf{P}^{-}$, so $\lVert\mathbf{F}_{\mathrm{T}}\rVert$ rises with engagement depth. (ii) Damage growth. The MGPCG damage solve grows $d_p$ in $\mathcal{N}(\mathrm{T})$. The factor $(1-d_p)^2$ reduces $\mathbf{P}^{+}$, while $\mathbf{P}^{-}$ preserves compressive support. (iii) Separation. When $d_p$ reaches $d^{*}$, CPIC transfers [HFG* 18] stop momentum sharing across the cut, and the wrench drops as the blade clears. This force profile distinguishes a clean cut from stalled penetration and is available to the policy through its force observations (§ 4.6).

Gym Interface and Policy Backbones

The simulator provides the standard Gym methods step(action) and reset(). The robot is a 7-DoF FR3 arm; the cutting / scooping / shoveling tool is held by friction in the parallel gripper rather than rigidly welded to the end-effector. The policy must therefore both execute the manipulation and maintain enough grip force to keep the tool from slipping under reactive contact wrench.

Observation. The observation $\mathbf{o}_t \in\mathbb{R}^{21+n_{\mathrm{task}}}$ concatenates joint state and a coarse tool/scene summary:

The vector contains the arm position $\mathbf{q}_t$ and velocity $\dot{\mathbf{q}}_t$, the tool-tip position $\mathbf{x}_{\mathrm{T}}^{\mathrm{tip}}$ and orientation $\mathbf{q}_{\mathrm{T}}^{\mathrm{tip}}$, plus a small task summary $e_t^{\mathrm{task}}$. The task summary usually contains at most four scalars, such as alive-particle ratio, tip-to-target distance, and tip height. The wrench $(\mathbf{F}_{\mathrm{T}},\boldsymbol{\tau}_{\mathrm{T}})$ is not given to the policy. It enters the rigid-body dynamics and deflects the arm from its commanded targets; the resulting velocity appears in the next observation. PPO normalizes $\mathbf{o}_t$ with running mean and variance estimates.

Action. The action $\mathbf{a}_t \in[-1,1]^8$ concatenates a normalized 7-DoF FR3 joint-position delta around a fixed home pose with a 1-DoF gripper-opening command:

Here $\odot$ is elementwise multiplication and $\mathbf{r}=(2.5,1.7,2.5,1.7,2.5,2.0,2.5)$ rad is the FR3 joint range. $w_{\max}$ is the commanded finger separation. Joint and gripper targets drive stiff PD controllers at 50 Hz, while the MPM update uses multiple substeps for stability. Because the tool is held by friction, the policy must coordinate arm motion and grip closure; if the reactive wrench overcomes friction, the tool slips and the cut is lost.

Reward. Each task defines its own reward as a weighted sum of interpretable shaping terms computed from the MPM field and rigid-body state. Common ingredients are tip-to-target proximity, depth progress into the material, blade-axis alignment, a damage-driven cut-quality term that rewards a clean separation band, and joint-effort and tool-slip regularizers; the per-task formulas and weights are given in § 5.2.

Policy backbone: PPO with optional BC warm-start. We use PPO [25] as the primary RL backbone, optionally warm-started from a scripted action prior. PPO is on-policy: each update uses a fresh rollout, so the critic and advantage depend on the latest batch. On the tofu cutting task, a randomly initialized single-env actor usually stays near the “do-nothing” return of $-27$ to $+27$; random exploration rarely forms a clean sequence of descent, engagement, and cutting. Two remedies help. A 64-env rollout averages gradient noise and reaches mean return $\sim+328$, removing $\sim2600$ particles per episode. A BC warm start initializes the actor from a scripted action prior: the scripted action reaches $\sim+325$, and subsequent single-env PPO refinement holds a deterministic return of $\sim+278$. Hyperparameters are given in § 5.3.

Experiments

Figure 4 shows the shared robot–tool arrangement used by the two closed-loop cutting case studies.

We evaluate the simulator on six tool-manipulation demos (Figure 1) that together span the multi-state regime introduced in § 4.4: four cutting tasks across a soft-to-brittle stiffness/toughness spectrum, plus a viscoplastic scooping scene and a granular shoveling scene.

Here, coverage measures how completely the tool slices through the material along the target cut plane. It is the fraction of eligible bins on that plane that show persistent separation on both sides.

We use simulation coverage below only for the breadth of scripted material demonstrations.

Tasks

We exercise MG-CD-MPM on six demonstrations grouped into three roles: two RL-trained cutting tasks (cut_tofu, cut_peanut_brittle), two additional cutting demos that fill in the soft-to-brittle spectrum (cut_jelly, cut_wax), and two non-cutting demos that exercise the viscoplastic and granular material classes (scoop_icecream, shovel_sand). All six share the Gym-like observation/action layout of § 4.6 and the per-particle constitutive layer of § 4.1; the three tools (knife, spoon, shovel) used across the six demos are shown in Figure 5. We narrow the closed-loop RL evaluation to two representative cutting tasks. A full MLS-MPM substep with the MGPCG damage solve in $\mathcal{N}(T)$ costs an order of magnitude more wall time than a comparable rigid-body rollout, so each end-to-end training run consumes several hours of single-GPU compute even with the patch optimization of § 4.3; this scaling is discussed in § 6.

Cutting (RL-trained). The two RL-trained cutting tasks (cut_tofu and cut_peanut_brittle) share the same kinematic setup: a planar knife held in the FR3 gripper cuts through a slab on a table. They differ in slab material parameters (Table 1) and reward weighting. cut_tofu uses an extra-firm tofu slab ($E = 200\,\mathrm{kPa}$, $\nu= 0.40$, $\sigma_F/E = 0.15$) paired with a connectivity-based cut-quality reward (clean band of removed particles with left/right symmetry) and a topological split bonus (weight 60, fires only when the alive-particle field has exactly two connected components), with 200-step episodes. cut_peanut_brittle uses the candy slab ($E = 1\,\mathrm{MPa}$, $\nu= 0.25$, $\sigma_F/E = 0.015$): a stiffer, much lower-toughness material paired with a simpler reward (rate of newly-separated particles minus proximity and effort penalties), with 250-step episodes, so success on cut_peanut_brittle is attributable to the damage solver of § 4.2 rather than to topological-quality shaping. Both tasks share the tip-to-target proximity term $\exp(-\lVert\mathbf{x}^{\mathrm{tip}}-\mathbf{x}^{\star}\rVert/0.08)$, the vertical depth-into-slab term, the blade-axis alignment term $\cos^2(\hat{\mathbf{n}}_{\mathrm{T}},-\hat{\mathbf{z}})$, and the joint-velocity, joint-effort and tool-slip penalties of § 4.6. Full reward breakdowns are in § 5.2.

Demo	$E$ (Pa)	$\nu$	$\sigma_F/E$
jelly	$6 \times10^4$	0.30	0.20
tofu	$2 \times10^5$	0.40	0.15
wax	$5 \times10^5$	0.30	0.05
peanut brittle	$1 \times10^6$	0.25	0.015
ice cream	$4 \times10^4$	0.30	0.30
sand	$2 \times10^3$	0.20	0.10

Cutting (simulation coverage). We also render two scripted cutting scenes with jelly and wax to sample additional material settings (Table 1). cut_jelly uses a very soft gel ($E = 60\ \mathrm{kPa}$, roughly $3\times$ softer than cut_tofu) with an elevated failure-stress ratio ($\sigma_F/E = 0.20$ vs. 0.15 for cut_tofu) that captures the high failure strain of gelatinous solids: the blade descends through the gel without producing a sharp crack tip and the two halves wobble apart under their own elasticity, exercising the regime where damage propagates through elastic shear rather than brittle fracture. cut_wax sits between cut_tofu and cut_peanut_brittle on the same axis ($E = 500\ \mathrm{kPa}$, $\sigma_F/E = 0.05$): a stiffer ductile slab whose damage profile is intermediate between the soft tearing of jelly/tofu and the clean snap of peanut brittle, confirming that the CD-MPM damage layer interpolates smoothly across stiffness/toughness combinations rather than only matching the two trained tasks. Together with the RL-trained pair, these four cutting demos cover roughly 1.5 orders of magnitude in stiffness (60 kPa to 1 MPa) within the same per-particle damage layer.

Non-cutting simulation coverage. The remaining two demonstrations (scoop_icecream and shovel_sand) exercise the viscoplastic and granular material classes of § 4.4 under scripted tool motion, both still on the same MLS-MPM+MG-CD-MPM backbone as the cutting demos; only the per-particle material parameters in Table 1 change. scoop_icecream (class V) is a semi-fluid lift: a hemispherical scoop enters an ice-cream tub, carves out a portion, and lifts it without the chunk fragmenting. The ice cream is modelled as a very soft elastic substrate ($E = 40\ \mathrm{kPa}$, $\nu= 0.30$, $\rho= 900\ \mathrm{kg}/\mathrm{m}^3$) paired with a high damage threshold ($\sigma_F/E = 0.30$), so the material flows around the spoon under sustained shear but the lifted chunk stays intact below threshold. The same damage solver provides a practical proxy for Herschel–Bulkley rheology without a separate viscoplastic kernel. shovel_sand (class G) drives a perforated sieve disc through

a 30 cm $\times$ 30 cm $\times$ 10 cm pile and the grains build a heap in front of the blade. We model the pile as a heavy ($\rho= 1600\ \mathrm{kg}/\mathrm{m}^3$) near-cohesionless medium: a deliberately low Young’s modulus ($E = 2\ \mathrm{kPa}$, $\nu= 0.20$) lets the bulk particles flow rather than rebound, and a moderate damage threshold ($\sigma_F/E = 0.10$) lets individual grains detach from the bulk as they cross the blade rim, so that the collective behaviour reads as granular under the same per-particle damage layer used for cutting and scooping. Together, the four cutting demos and these two non-cutting scenes confirm that the same Gym-like interface, observation/action layout, FR3 rigid-body coupling, and MLS-MPM+MG-CD-MPM solver support all four material classes by adjusting per-particle material parameters alone. Reward functions and RL training for the four scripted demos (cut_jelly, cut_wax, scoop_icecream, shovel_sand) are deferred to future work (§ 6).

Reward Functions

We express the rewards for the two RL-trained cutting tasks, cut_tofu and cut_peanut_brittle, as weighted sums of interpretable shaping terms computed from the MPM particle field and the rigid-body state. We use the following notation throughout. $\mathbf{x}^{\mathrm{tip}}=(x^{\mathrm{tip}},y^{\mathrm{tip}},z^{\mathrm{tip}})$ is the tool tip position in world coordinates and $\hat{\mathbf{n}}_{\mathrm{T}}$ is the tool’s body $\hat{\mathbf{z}}$ axis for orientation. $\mathbf{x}^{c}=(c_x,c_y,c_z)$ is the geometric centre of the material block, and $z_{\mathrm{top}},z_{\mathrm{bot}}$ are the block’s top and bottom $z$-coordinates with height $h=z_{\mathrm{top}}-z_{\mathrm{bot}}$. The FR3 tool centre point (TCP) is the standard end-effector reference frame defined between the gripper fingers, as supplied by the FR3 URDF; $z^{\mathrm{TCP}}$ denotes its world $z$, and $z^{\mathrm{knife}}$ denotes the world $z$ of the knife body’s centre of mass. When the gripper closes correctly on the knife handle, the centre of mass sits 2.25 cm above the TCP origin ($z_{\mathrm{knife}}=z^{\mathrm{TCP}}+0.0225\ \mathrm{m}$), which is the equilibrium configuration we use as the slip baseline below. $\mathbf{q}\in\mathbb{R}^{7}$ is the arm joint-velocity vector, $n_t^{\mathrm{alive}}$ is the number of MPM particles still active at step $t$, $\Delta n_t^{\mathrm{dead}}=n_{t-1}^{\mathrm{alive}}-n_t^{\mathrm{alive}}$, and $\Omega_t^{\mathrm{alive}}\subset\mathbb{R}^{3}$ is the voxelised support of the alive-particle density at step $t$; we write $\#\mathrm{cc}(\Omega_t^{\mathrm{alive}})$ for the number of 6-connected components of that voxel set. $\mathbf{a}_t$ is the action emitted by the policy at step $t$ (incremental joint targets, § 4.6), and $\mathbf{a}_t^{\mathrm{demo}}$ is the scripted demonstration action at the same step (the open-loop sequence of descent, holding, and retraction used for the behavior-cloning warm-start of § 5.3). $\mathbb{1}[\cdot]$ denotes the indicator function (1 if the predicate holds, 0 otherwise). Per-task reward terms below use SI units (metres, radians); all numerical constants inside indicators and exponentials are in metres unless stated otherwise, and no global per-step clipping is applied.

Soft slab reward. The cut_tofu reward decomposes as

with components

Here $b_{\mathrm{clear}}=1-n_t^{\mathrm{band}}/n_0^{\mathrm{band}}$ measures how much of an axial band of half-width 1.5 cm centered on the cut plane $y=c_y$ has been cleared of alive particles, and $b_{\mathrm{bal}}=1-|n_L-n_R|/n_{\mathrm{far}}$ measures left/right balance of the alive particles outside the band ($n_L,n_R$ count alive particles with $y<c_y$ and $y>c_y$ respectively). The split bonus $r_{\mathrm{split}}$ counts the connected components of the alive-density volume on a voxel grid and fires only when there are exactly two, since "the two halves are balanced" does not imply "the two halves are topologically separated". An uncut bottom strip can leave a single component even when the cut quality looks high. The action-tracking term $r_{\mathrm{track}}$ is gated to the cut zone (tip within 10 cm above the block top): without the gate, an early failure mode is for the policy to mimic the scripted demonstration while waving the tool in free air, harvesting tracking reward without ever engaging the material. The two velocity regularizers act at different scales: $r_{\mathrm{vel}}$ is the main penalty for wrist shaking ($\sim\pm0.1$/step at typical joint speeds) while $r_{\mathrm{eff}}$ provides a smooth downward pressure on the action norm without dominating the gradient. The drop penalty $r_{\mathrm{drop}}$ fires when the knife body’s $z$-coordinate falls below the TCP’s $z$ (with the 2.25 cm rest offset accounted for). With the policy-controlled gripper width $a_t^{\mathrm{grip}}$ of Sec. 3.6, this term doubles as both a friction-loss penalty (reactive wrench overcoming the commanded grip) and a direct deterrent against the policy opening the gripper to drop the knife: a 5 cm slip costs 0.1/step, an order of magnitude below the cutting reward, so it nudges the policy toward keeping the tool gripped rather than dominating the learning signal.

Choice of outer weights and length scales. The relative magnitudes of the six shaped weights encode a coarse curriculum. The dense approach terms $r_{\mathrm{reach}},r_{\mathrm{depth}},r_{\mathrm{vert}}\in[0,1]$ are weighted at 1.5, 1.0, 0.5 respectively: together they provide the bulk of the early-episode gradient while the tool is still approaching, and they saturate near 1 once the blade reaches the block, at which point $r_{\mathrm{cut-q}}$ at weight 15 becomes the only dense signal that continues to grow and takes over the gradient. The terminal sparse bonus $r_{\mathrm{split}}$ is weighted at 60 so that a single successful split contributes roughly twice the integrated cut-quality reward over a $\sim200$-step episode, letting topological success dominate near-misses; $r_{\mathrm{track}}$ is held at 0.5 because it is an auxiliary descent prior, and a larger weight would make the policy reproduce the scripted demo even when material engagement has failed. The inner length scales follow the workspace geometry. The 0.08 m inside $r_{\mathrm{reach}}$ is the radius at which the exponential reach reward falls to $e^{-1}$, chosen so the gradient is informative across a 5∼10 cm reach error rather than collapsing to zero. The 0.02 m offset in $x^\star=(c_x,c_y,z_{\mathrm{top}}-0.02)$ places the reach target 2 cm below the block top, so a converged policy penetrates the surface rather than hovering above it. The decay constant 2 inside the $r_{\mathrm{track}}$ exponential gives half-reward at an action deviation of $\sim0.35$ rad/step from the demo, on the same scale as one PPO update’s policy noise. The 1.5 cm band half-width in $b_{\mathrm{clear}}$ matches the blade thickness plus a small empirical tolerance for damage spread, so the band contains exactly the particles that a successful cut should remove. Finally, the $10^{-8}$ inside the $r_{\mathrm{vel}}$ square root is a numerical smoothing constant that keeps the gradient bounded at zero joint velocity and has no effect on the magnitude of the term elsewhere.

Hard slab reward. The cut_peanut_brittle reward is built on the cut_tofu shaping above: we keep the joint-effort regularizer (reweighted $5\times$ heavier) and a tip-to-block proximity term, drop the topological-quality and alignment terms ($r_{\mathrm{cut-q}}$, $r_{\mathrm{split}}$, $r_{\mathrm{track}}$, $r_{\mathrm{drop}}$, $r_{\mathrm{vert}}$, $r_{\mathrm{depth}}$), replace the soft-block reach reward with a linear distance penalty to the block centre, and add an explicit per-step kill-count term that rewards the damage solver directly:

The reduced shaping is intentional: cut_tofu’s $r_{\mathrm{cut-q}}$ and $r_{\mathrm{split}}$ can in principle be earned by any policy that cleanly displaces the band of particles between the two halves, regardless of whether the displacement comes from CD-MPM damage propagation or from raw particle-deactivation artefacts. Stripping those terms on the stiffer slab and rewarding particle deactivation directly through $\Delta n_t^{\mathrm{dead}}$ attributes success unambiguously to the damage solver of Sec. 3.2, which is what we ablate in Sec. 4.4. The three coefficients of $r_t^{\mathrm{hard}}$ are calibrated so that the kill-count term dominates while the regularizers remain secondary: at $\Delta n_t^{\mathrm{dead}}\sim10$ particles/step typical of a working cut, $0.10\Delta n_t^{\mathrm{dead}}$ contributes $\sim+1$/step, against $\sim-0.025$/step from the $0.05\lVert x^{\mathrm{tip}}-x^c\rVert$ distance term at a typical 0.5 m workspace distance and $\sim-0.01$/step from the joint-effort term. The joint-effort coefficient 0.001 is $5\times$ heavier than cut_tofu’s $r_{\mathrm{eff}}=-2\times10^{-4}$ to account for the larger joint torques required to drive the stiffer slab, which would otherwise let brute-force motion substitute for controlled cutting.

Implementation and Hyperparameters

Simulator and robot stack. All physics runs on NVIDIA WARP [citep][Mac22] as a unified GPU backend. The CD-MPM solver of Sec. 3.2 discretizes the equations of motion with explicit MLS-MPM [citep][HFG*18]: mass and momentum are transferred to a fixed background grid with $\sim16$K active cells via APIC, the grid update is a single grid Warp pass per substep (no global linear solve), and the per-particle damage update of Sec. 3.2 is a local closed-form Cardano eigenvalue of the symmetric stress tensor. The articulated FR3 arm and the free-rigid tool are stepped by a Warp-resident MuJoCo solver (mujoco_warp) that shares the device and CUDA stream with MPM, so the two-way coupling kernels of Sec. 3.5 (compute_body_forces, subtact_body_force) exchange impulses on-device without any host round-trip. Control runs at 50 Hz, with per-control MPM and rigid substeps tuned per task to satisfy the Courant–Friedrichs–Lewy (CFL) stability condition: cut_tofu uses voxel size 0.04m with $1+1$ substeps ($\sim19$K particles per env), cut_peanut_brittle voxel 0.018m with $3+3$ substeps, scoop_icecream voxel 0.018m with $2+2$ substeps, and shovel_sand voxel 0.020m with $2+2$ substeps. The tool is held by contact friction in the FR3 fingers (gripper PD gains $k_p=5\times10^4$, $k_d=700$, friction $\mu=2.0$) rather than welded to the end-effector, which is what makes the gripper-slip penalty in Sec. 3.6 a meaningful signal. The production damage solver uses outer PCG with a full geometric V-cycle preconditioner, a relative-residual tolerance of $10^{-5}$, and a cap of 20 outer iterations. Each V-cycle applies three weighted-Jacobi pre- and post-smoothing sweeps and 24 coarse-level iterations; restriction and prolongation remain on the GPU.

Rendering. Both rendering paths share the same geometry pipeline: alive MPM particles are voxelized into a density field, Gaussian-smoothed, and converted into a triangle mesh by marching cubes. The lightweight path rasterizes the mesh through Newton’s Warp-resident viewport for training-time inspection, while the offline path imports the mesh into Mitsuba 3 [JSR*22] for Monte Carlo path-traced figures with physically based area lighting and a smooth-dielectric model for the brittle-jelly material.

Two-stage training pipeline with MPM-derived force jittering. The full MLS-MPM environment limits training throughput: a single env runs at $\sim12$ fps for raw simulation and $\sim5$–$6$ fps with PPO updates included. Training $2\times10^5$ on-policy steps at this throughput requires substantial wall-clock time. We exploit the observation that the gripper-and-arm sub-problem (track a target with the tool firmly held, reject an unknown end-effector load) is a pure rigid-body control problem and can be pre-trained on a much cheaper proxy environment. Concretely, we use a two-stage pipeline.

Stage 1: rigid-body pre-training with MPM-calibrated wrench jitter. We first run a small set of scripted MLS-MPM rollouts of the target task and collect, via mpm.collect_collider_impulses, the per-step net force $\mathbf{F}^{\mathrm{mpm}}$ and torque $\boldsymbol{\tau}^{\mathrm{mpm}}$ that the cutting material applies to the tool. From these we estimate the empirical 95th percentile of the force and torque magnitudes, written $p^F_{95}$ and $p^\tau_{95}$: that is, in 95% of simulation steps $\lvert\mathbf{F}^{\mathrm{mpm}}\rvert\le p^F_{95}$, and analogously for the torque. We use the 95th percentile rather than the maximum because the maximum is dominated by transient impact and fracture spikes that misrepresent the typical contact load, and rather than the mean because the mean is pulled down by the many steps in which the tool is not engaged with the material; on the cutting tasks the estimates are $p^F_{95}\approx0.5$ N and $p^\tau_{95}\approx5$ N$\cdot$m. We then train on a pure rigid-body environment built on mujoco_warp, with no MPM at all, in which at every control step we sample each Cartesian component of the applied force and torque independently and uniformly within these calibrated bounds:

and apply the resulting wrench to the tool body for that step. This injects a disturbance with the same characteristic magnitude that MPM contact will eventually deliver, without running the MPM simulation. Because the entire stage avoids MPM, we batch $\sim8192$ FR3+tool worlds inside a single mujoco_warp nworld model via newton.ModelBuilder.replicate() and train on one GPU; observation and action layouts match Eqs. (25)–(26) exactly so the resulting policy is checkpoint-compatible with Stage 2.

Stage 2: MG-CD-MPM fine-tuning. We load the Stage 1 checkpoint with the --init-checkpoint flag and continue PPO inside the full MG-CD-MPM environment of Sec. 3. This environment combines the tool-centered MGPCG damage solve with the MLS-MPM momentum substrate. We use a single-env, single-process Warp rollout on the same GPU. This is where the policy specializes from “track a target while rejecting an unknown end-effector load” to “engage the material, drive damage along the desired plane, release on completion.” Stage 2 uses $4\times10^4$–$10^5$ MG-CD-MPM steps to fine-tune the Stage 1 policy.

Networks. PPO’s actor-critic is a two-hidden-layer MLP of width 256 with tanh activations and orthogonal initialization (output-layer gain 0.01 on the policy head, 1.0 on the value head). The policy is a diagonal Gaussian with a learnable, state-independent log $\sigma$ initialized at $-3$ so the initial policy stays close to zero. Observations are normalized online by a running mean/variance estimate updated at every rollout step; without this normalization the joint angles ($\sim3$ rad) dominate target-offset components ($\sim5$ cm) and the policy fails to learn the millimetre-scale corrections demanded by the cutting tasks.

PPO. We use the clip variant of PPO [SWD*17] with rollout length 2048 steps per env, 10 update epochs, minibatch size 256, learning rate $3\times10^{-4}$ (Adam, default betas), clip range $\epsilon=0.2$, value coefficient 0.5, no entropy bonus, GAE($\gamma=0.99,\lambda=0.95$), gradient-norm clipping at 0.5, and a target-KL early stop at 0.02 that aborts the remaining update epochs if the empirical KL of the last epoch exceeds threshold. Advantages are recomputed per update across the $T\times N$ rollout; value-loss clipping is not used. Single batches whose gradients contain non-finite entries (occasional on the cutting tasks when a fracture event injects a spike into the reward) are detected and silently skipped to avoid invalid optimizer updates.

PPO variants. On the cutting tasks (cut_tofu, cut_peanut_brittle) a single-environment PPO actor trained from a random initialization collapses to a do-nothing return in $[-27,+27]$: random Gaussian on-policy exploration almost never produces a clean sequence of descent, engagement, and cutting, and the advantage signal is overwhelmed by noise. We address this through two variants. (i) Vectorized PPO rolls out 64 parallel environments and aggregates the gradient across them, stabilizing training to a mean return of $+328$ on cut_tofu with $\sim2600$ particles removed per episode; this is the strongest cutting policy we report. (ii) PPO with BC warm-start initializes a single-environment PPO actor from the behavior-cloning prior described next, yielding a deterministic return of $\sim+278$ that subsequent PPO refinement preserves.

Behavior cloning (BC) and critic warm-start. Stage 1 itself begins with a behavior-cloning pass that initializes the policy from a hand-scripted action sequence. The scripted trajectory is a deterministic sequence of descent, holding, and retraction for cutting tasks and analogous hand-coded primitives for shoveling and scooping; on cut_tofu it reliably produces an episode return of $\sim+325$, while a freshly-initialized random policy plateaus at $+100 \sim+150$ because random Gaussian on-policy actions almost never form a clean sequence of descent, engagement, and cutting and the critic never sees an example of success. We use the demo for two purposes. First, an actor BC pass (20 epochs at learning rate $10^{-3}$) fits the diagonal-Gaussian policy by minimizing the negative log-likelihood $-\sum_t \log\pi_\theta(a_t^{\mathrm{demo}} \mid o_t)$ on a dataset of $\sim4{,}000$–$64{,}000$ demonstration tuples collected by replaying the scripted action across 64 parallel rigid-body envs; the MPM-calibrated wrench jitter of Stage 1 gives every env a different perturbation history, providing the observation-side diversity that lets BC generalize beyond the exact demo trajectory rather than collapsing onto it. After BC, the policy’s per-dimension log$\sigma$ is floored at $-2$ ($\sigma\approx0.13$): without this floor NLL drives log$\sigma$ to very negative values on near-deterministic action dimensions, and the first PPO update’s empirical KL diverges. Second, an optional critic BC (“BC-critic”) fits the value head on the same demonstration data to the Monte-Carlo discounted return, $V_\phi(o_t) \approx\mathbb{E}\left[\sum_{k\geq0}\gamma^k r_{t+k} \mid o_t\right]$, with $\gamma$ matched to the PPO discount. Without this prefit the critic is random and the first PPO iteration produces advantages $A=R-V$ of the same order as the raw return, which once again causes a sharp increase in KL divergence and destabilizes the BC-initialized actor; with the prefit the iter-1 advantage stays bounded and PPO refines the BC policy stably. The combined procedure (BC actor + BC critic + log$\sigma$ floor) is what produces the “PPO with warm-start” rows in the comparison below.

Compute budget. All experiments in this paper were run on a single NVIDIA RTX 4090 GPU; both stages share the same device. Stage 1, despite hosting 8192 rigid-body envs in a single mujoco_warp nworld, fits on one GPU because the model contains no MPM. Stage 2 is the throughput bottleneck: a single MLS-MPM env runs at $\sim12$fps for raw simulation (scripted action, no policy) and $\sim5$–$6$fps once policy inference, GAE, and the PPO update are layered on; the Stage 2 budget is $4\times10^4$–$10^5$ MPM steps per (task, seed). We train 3 seeds per configuration and report mean$\pm$one standard deviation; checkpoints are written every 5 outer iterations and the best running mean return is preserved.

Baselines and Ablations

We evaluate four aspects: solver speed at matched residual, momentum conservation in a free collision, policy sensitivity to damage modeling, and crack quality across tools and materials. Unless noted otherwise, timings use one NVIDIA RTX 4090.

Localized MGPCG Efficiency

We separate localization from multigrid acceleration by solving the same 3D implicit phase-field system in four ways: full-grid Jacobi-PCG, localized Jacobi-PCG, localized PCG with a symmetric two-level preconditioner, and localized PCG with a full V-cycle (our MGPCG solver). All four use the same right-hand side, initial guess, operator, and relative-residual tolerance $10^{-5}$. Five timed CUDA repeats include all kernels and patch transfers. Table 2 reports the times and iteration counts, separating localization from preconditioning.

Grid	Full PCG ms	Full PCG iter.	Local PCG ms	Local PCG iter.	Two-level PCG ms	Two-level PCG iter.	Local MGPCG (Ours) ms	Local MGPCG (Ours) iter.
$192^3$	139.38	149	40.15	144	54.45	11	35.80	4
$256^3$	694.36	199	53.90	195	73.96	15	35.98	4
$320^3$	1741.15	249	115.57	246	105.61	19	38.60	4

Localization alone yields $3.47\times$, $12.88\times$, and $15.07\times$ speedups at $192^3$, $256^3$, and $320^3$; full multigrid preconditioning raises these to $3.89\times$, $19.30\times$, and $45.11\times$. At the default patch the localized field differs from the full reference by $0.65\%$. Increasing patch half-width from $3.75\ell_c$ to $7.5\ell_c$ reduces error from $1.61\%$ to $0.074\%$, while solver time increases from 35.0 to 48.4 ms. Thus localization is accurate when its boundary is separated from the active damage band, not necessarily exact in every scene.

In the production Warp solver, the same 225-frame, 806k-particle tofu trajectory reaches median residual $2.83\times10^{-6}$ (maximum $9.87\times10^{-6}$) in 8 MGPCG iterations, costing 319.56 ms for MPM and 330.44 ms total per control frame. The previous fixed three-cycle path costs 272.74 ms but stops at median residual $1.06\times10^{-2}$; it is faster only because it is under-solved. The final damage fields differ by 4.56% in relative $L_2$.

Two-Way Coupling Conservation

For conservation, we remove gravity, controls, and the floor. A translating free MPM block collides with a stationary 2-kg free rigid box using the production SDF contact and impulse-accumulation path. We use a free body rather than an FR3-held tool because the robot’s fixed base, joint constraints, and actuators exchange external impulse with the world. The free-body setup therefore isolates the equal-and-opposite exchange in the coupling kernels. The block travels at 1.5 m/s along $+x$. In the centered case, its impulse line passes through the box center of mass, up to a $10^{-4}$ m offset that avoids an ambiguous symmetric SDF face. This mainly tests linear-momentum exchange. In the eccentric case, the block center is displaced by 3.5 cm in $y$, which adds a moment arm and tests angular as well as linear momentum. Twelve runs combine these contacts with grid spacings 0.012/0.010/0.008 m and timesteps $2\times10^{-4}/10^{-4}$ s. The audit uses the CPU backend of the same kernel and has no external impulse.

At each step we compute total linear momentum $\mathbf{P}$ and angular momentum $\mathbf{L}$ about the world origin. $\mathbf{L}$ includes the MPM/APIC orbital and affine terms, plus rigid-body orbital and spin terms. We report $e_P=\|\mathbf{P}(t)-\mathbf{P}(0)\|/P_s$ and $e_L=\|\mathbf{L}(t)-\mathbf{L}(0)\|/L_s$, with $P_s=\max(\|\mathbf{P}(0)\|,\epsilon)$ and $L_s=\max(\|\mathbf{L}(0)\|,\|\mathbf{P}(0)\|\max(\|x_T(0)\|,0.1\ \mathrm{m}),\epsilon)$. The stepwise columns use the corresponding residuals $\Delta\mathbf{P}_{\mathrm{MPM}}+\mathbf{J}_T$ and $\Delta\mathbf{L}_{\mathrm{MPM}}+\mathbf{x}_T\times\mathbf{J}_T+\mathbf{K}_T$, with the same normalization.

Table 3 reports the maximum error over all grid spacings, timesteps, and simulated steps. Both the global conservation errors and per-step action–reaction residuals remain at approximately $10^{-6}$ or below. The eccentric case is the stricter angular-momentum test because of its nonzero lever arm; its maximum $e_L$ is $1.33\times10^{-6}$. These results validate equal-and-opposite linear and angular impulse exchange to numerical precision. Angular impulse is accumulated about the rigid body’s center of mass; an articulated-body bridge must shift the moment if its collider reference differs from it.

Contact	$e_P$	$e_L$	Step P	Step L
Centered	$1.06\mathrm{e}{-6}$	$8.35\mathrm{e}{-7}$	$3.36\mathrm{e}{-7}$	$2.67\mathrm{e}{-7}$
Eccentric	$1.53\mathrm{e}{-6}$	$1.33\mathrm{e}{-6}$	$6.50\mathrm{e}{-7}$	$5.69\mathrm{e}{-7}$

Does Damage Fidelity Improve the Policy?

We construct paired compact Gym ablations. They share the same observation definitions, actions, rewards, initial particles, robot and MPM configurations, PPO network, BC warm start, 10,240-step budget, and three seeds. Their observations use the same robot-state-plus-task-summary interface as the preceding tasks, with scalars for cut completion, tip-to-center distance, tip height, cut-plane error, and contact force. The internal damage model is the only change. The full model uses stress history, localized phase-field MGPCG, damage-dependent stress degradation, and thresholded CPIC separation. The simplified model assigns persistent CPIC sides directly from the blade sweep. Both models use the same randomly positioned particle region intersecting the target plane, but only the full model represents its increased fracture resistance; the simplified model is material-agnostic. The inclusion location is hidden from the policy, although contact-force feedback remains visible. For this ablation, we divide the nominal $y = 0$ cut plane into 1.5-cm $x$–$z$ bins and retain bins that initially contain particles on both sides. We count a bin as covered when persistent CPIC labels appear on both sides. Success requires at least 80% coverage and cut-plane error at most 2.5 cm.

Both policies solve the simplified test, but the full-trained policy has higher success and coverage under heterogeneous full physics. This pilot study provides a controlled comparison across three random seeds under the same training budget. The compact ablation welds the tool to isolate separation fidelity and therefore does not measure gripper slip.

Tool–Material Crack Quality

To test fracture fidelity beyond one example, we combine three controlled 3D tool-shaped driving fields (thin straight, thick straight, and wedge) with soft/tough, medium, and hard/brittle materials. For each pair, full-grid PCG and localized MGPCG solve the same phase-field system to $10^{-5}$. This is a solver validation, not nine end-to-end robot rollouts. Table 5 reports damage-field, fracture-energy, and crack-surface discrepancies.

Material	Tool	$L_2$ (%)	$E_f$ (%)	$H$ (mm)
Soft/tough	Thin	0.132	1.196	0.0
Soft/tough	Thick	0.235	1.285	2.5
Soft/tough	Wedge	0.193	0.877	2.5
Medium	Thin	0.114	0.664	2.5
Medium	Thick	0.200	0.722	2.5
Medium	Wedge	0.154	0.435	2.5
Hard/brittle	Thin	0.088	0.302	2.5
Hard/brittle	Thick	0.153	0.334	2.5
Hard/brittle	Wedge	0.108	0.171	2.5

Across all combinations, damage-field error is 0.088%–0.235%, median fracture-energy error is 0.66% and the maximum is 1.29%, and maximum Hausdorff distance is one 2.5-mm cell. Soft/tough material under the thin field produces no thresholded surface in either solver and is recorded as a consistent no-crack result. Maximum separated-volume discrepancy is $2.1 \times10^{-5}$ of the domain. Thus, across this tool–material matrix, restricting the solve to the tool-centered patch preserves the full-grid crack field, fracture energy, and extracted crack surface to within the reported discretization errors.

Two-Dimensional CD-MPM Comparison

The localized phase-field construction also applies in two dimensions. Figure 6 shows a 2D cutting example at $t = 0.50\,\mathrm{s}$ on a $63 \times41$ uniform grid. The left panel is the CD-MPM reference. The right panel applies the phase-field V-cycle only in the tool-centered $16 \times24$ region and keeps the field outside it. This produces a spatially regularized damage band. Because the example uses a prescribed number of V-cycles, it is for qualitative visualization rather than matched-accuracy timing or speedup claims; the controlled 3D comparison appears in Table 2.

Simulator baselines. The Full-PCG baseline implements the original CD-MPM phase-field system on the complete uniform grid with Jacobi-preconditioned CG. It shares our particle state, constitutive update, transfers, contact, right-hand side, and residual tolerance; only damage localization and the preconditioner differ. This is the reference used in Table 2 and the $3 \times3$ crack matrix.

Kill is the simplest legacy rule: particles swept by the blade are removed from subsequent P2G/G2P transfers. It deletes mass and stored elastic energy, produces sparse grid mass and impulse spikes, and is not a conservative or stable RL environment. Standalone Kill diverges at blade contact in our implementation, so it produces no complete rollout. The Simplified-separation baseline of Table 4 remains trainable because it preserves particles and momentum with persistent CPIC labels. Figure 7 shows the difference: ours gives a stable RL-policy rollout, whereas the Kill-based closed loop diverges at contact. For a direct visual comparison, we replay the successful trajectory open-loop with a welded tool on tofu and apply the Kill mask to that same stable MPM trajectory.

Train	Test	Success	Coverage
Full	Full	3/3	84.85±0.66%
Full	Simplified	3/3	100.0±0.0%
Simplified	Full	2/3	76.89±13.46%
Simplified	Simplified	3/3	100.0±0.0%

RL Backbone Ablations

In our main cutting tasks, successful trajectories are difficult to discover through random exploration. On cut_tofu, single-environment PPO from random initialization collapses near the do-nothing return, while a 64-environment vectorized run reaches approximately +328 mean return and a BC-warm-started single-environment policy retains approximately +278. These observations motivate the curriculum of § 5.3; they are distinct from the matched three-seed damage-model comparison in Table 4. Figure 8 summarizes these three training behaviors on the main cut_tofu task.

We also report the SAC ablation to assess whether changing the RL algorithm addresses the exploration difficulty. Across four attempted variants, the entropy-regularized actor moves the blade away from the narrow trajectory of descent, engagement, and cutting and remains at the do-nothing floor. This negative result does not establish a general disadvantage of SAC; it records the failure mode under our observation, reward, and training budget. Figure 9 shows a representative rollout exhibiting this behavior.

Conclusion

We presented MG-CD-MPM, a Gym-like MPM simulator for RL-based dexterous manipulation of multi-state materials. Its primary algorithmic contribution couples stress-driven phase-field damage to a tool-centered localized solve on a uniform MPM background grid and uses residual-controlled MGPCG to accelerate that local system. Combining it with two-way rigid-tool coupling on an FR3-arm Gym interface in which the tool is held by the gripper rather than welded to the end-effector, MG-CD-MPM covers soft solids, brittle materials, granular media, and viscoplastic media within a single environment. We exercised the simulator on six demonstrations: RL-trained PPO agents (vectorized and BC-warm-started variants) on two cutting tasks that bracket the stiffness/toughness axis (cut_tofu, cut_peanut_brittle), plus scripted simulation-coverage scenes for two intermediate cutting demos (cut_jelly, cut_wax) and the viscoplastic (scoop_icecream) and granular (shovel_sand) material classes.

Limitations and future work. The dominant limitation of MG-CD-MPM is the wall-clock cost of MPM itself. Even with MGPCG restricted to $\mathcal{N}(T)$, a single full-physics rollout costs an order of magnitude more time than a rigid-body environment, and each end-to-end RL training run on cut_tofu or cut_peanut_brittle consumes several hours of single-GPU compute. This is why we focused the RL evaluation of this paper on cutting and treated the viscoplastic (scoop_icecream) and granular (shovel_sand) scenes as simulation-coverage demonstrations rather than as RL benchmarks. Closing this gap is the natural next step: extending the two-stage pipeline of § 5.3 with broader proxy distillation, parallel MGPCG batching across patches, or amortized damage solves should make multi-task RL training tractable on a single workstation. Beyond compute, brittle damage parameters are currently tuned per material class rather than learned, friction uses a single Coulomb $\mu$ per material, the tool-tracking patch assumes a single primary gripper, and transfer to a physical robot remains to be evaluated.

Acknowledgements

We thank Libin Liu, Siyuan Luo, and Ziqiu Zeng for helpful discussions and insights.