arXiv CSopen access

MoSA: Motion-constrained Stress Adaptation for Mitigating Real-to-Sim Gap in Continuum Dynamics via Learning Residual Anisotropy

MoSA: Motion-constrained Stress Adaptation for Mitigating Real-to-Sim Gap in Continuum Dynamics via Learning Residual Anisotropy

Jiaxu Wang† * 1 2 Junhao He * 1 Jingkai Sun 3 Yi Gu 1 Yunyang Mo 1 Jiahang Cao 3 Qiang Zhang 1 Renjing Xu† 1

arXiv:2605.22597v1 [cs.LG] 21 May 2026

Abstract

transfer. Project Page is available at https: //mercerai.github.io/MoSA/.

Learning real-world dynamics from visual observations is crucial for various domains. A common strategy is to calibrate simulators by estimating physical parameters, yet accuracy is ultimately bounded by the underlying physical models, which often assume materials are homogeneous and isotropic. Even if reasonable, realworld objects typically exhibit mild anisotropy and heterogeneity. After the near-isotropic backbone is well calibrated, these residual effects become the key bottleneck for further closing the real-to-sim gap. Although neural networks can fit dynamics end-to-end, such black-box modeling discards strong physical priors, leading to poor data efficiency and overfitting. Therefore, we propose MoSA, a motion-constrained stress adaptation framework that targets these residual effects to further improve real-to-sim dynamics learning. MoSA uses an isotropic model as a physics prior and learns residual stress operators to capture mild anisotropy and heterogeneity. It progressively adapts stresses via microplane-constrained redistribution in a physics-informed cascaded network. We further impose motion constraints by supervising temporal and spatial derivatives of the deformation field. Experimentally, our learned dynamics achieves superior accuracy, generalization, and robustness, while learning physically meaningful residual anisotropy. Finally, we validate MoSA in a robot manipulation setting, showing that better real-to-sim dynamics modeling translates into more reliable sim-to-real

1. Introduction Simulating realistic dynamics of diverse materials is essential for many applications in graphics, robotics, and embodied agents, and it is a key ingredient for building interactive digital twins of the real world. Visual observations are widely available and naturally capture rich deformation and motion signals. As a result, learning real-world dynamics from visual data has become an active research direction. A prevailing paradigm is to represent dynamics with classical physical models and identify objects by estimating a small set of physical parameters (Cai et al., 2024; Li et al., 2023). In practice, most pipelines adopt homogeneous and isotropic material models, since they provide a simple and often effective approximation for many everyday objects and enable stable differentiable simulation. However, this approximation is rarely exact: even when an object can be reasonably treated as homogeneous and isotropic, it typically exhibits small but systematic anisotropy and heterogeneity in the real world (Li et al., 2025; Fu et al., 2025). These deviations, which we call residual effects in continuum, can arise from internal structure, gradual wear, etc. As a result, estimating parameters within an isotropic model can capture the dominant backbone, but it often fails to account for these residual effects, which become a key bottleneck for further closing the real-to-sim gap. An alternative is to learn dynamics end-to-end with neural networks (Mittal et al., 2025; Cao et al., 2025; Zhobro et al., 2025). While this can model missing physics, directly fitting full dynamics discards strong physical priors and does not disentangle the dominant isotropic backbone from weak residual effects. The network must re-learn the backbone and capture subtle residuals at the same time, which often leads to poor data efficiency, unstable training, and overfitting.

* Equal contribution 1 Hong Kong University of Science and Technology, Guangzhou, China 2 MMLab, Chinese University of Hong Kong, Hong Kong SAR 3 The University of Hong Kong, Hong Kong SAR. Correspondence to: Jiaxu Wang <[email protected]>, Renjing Xu <[email protected]>.

Motivated by this bottleneck, we take a targeted modeling approach: we keep a calibrated near-isotropic backbone to explain the dominant behavior, and learn only the residual effects that arise from mild anisotropy and heterogene-

rd

Proceedings of the 43 International Conference on Machine Learning, Seoul, South Korea. PMLR 306, 2026. Copyright 2026 by the author(s).

1

Title Suppressed Due to Excessive Size

ity. Concretely, we treat an isotropic constitutive law as a physics prior and introduce a structured residual stress adaptation operator that progressively corrects the prior stress response and encourage physically consistent stress redistribution. This design preserves strong physical inductive bias while allocating learning capacity to the missing residual dynamics.

ture states from current geometry (Jing et al., 2026; Xue et al., 2024; Mittal et al., 2025; Zhobro et al., 2025; Chen et al., 2025c). Graph-based simulators (Wei & Freris, 2025; Zhang et al., 2024) employ GNNs to achieve plausible results across materials, but can overfit due to their high parameter count. To mitigate this, some methods inject physical structure by embedding neural components into PDEs; for instance, (Mittal et al., 2025) learns a unified latent constitutive law to represent material types. NCLaw (Ma et al., 2023) replaces constitutive laws with neural networks. DEL (Wang et al., 2024a) extends the Discrete Element Method with neural operators to model diverse materials. MASIV (Zhao et al., 2025) supervises particle trajectories with isotropic neural dynamics, while our method uses higher-order motion constraints and learns anisotropic residual stress corrections. OmniPhysGS (Lin et al., 2025) adopts discrete isotropic experts for plausible dynamics, whereas our method models continuous heterogeneity and anisotropic residuals for accurate real-to-sim dynamics. Recent efforts such as (Cao et al., 2025; Chen et al., 2025b) further emphasize incorporating physical priors to improve generalization. Despite this progress, many neural dynamics models effectively learn full dynamics or full operators. This often underuses strong explicit priors and does not explicitly separate the dominant near-isotropic behavior from the weaker residual effects that matter most in real-to-sim settings, which can hurt data efficiency and generalization. In contrast, we keep an isotropic constitutive prior as the backbone. We then learn residual stress corrections that target mild anisotropy and heterogeneity to breakup the bottleneck.

Another challenge in learning dynamics from vision is that supervision is often indirect and underconstrained. Earlier attempts assume access to full 3D geometry trajectories (Ma et al., 2023; Jaques et al., 2019), which is hard to obtain in real-world settings. More recent methods (Chen et al., 2025a) connect videos to 3D simulation through differentiable simulation (Bolliger et al., 2025) and neural rendering (Kerbl et al., 2023; Wang et al., 2024b). However, optimization driven only by image reconstruction remains ill-conditioned: different 3D deformation trajectories and material responses can produce similar image appearances. To mitigate this, we add motion-aware constraints that provide more direct supervision than pixel reconstruction alone. We first obtain a dynamic 3D reconstruction from the input videos, which provides reliable cues about how the object moves and deforms over time. We then use these cues to supervise the temporal and spatial derivatives of the deformation field. Together with our targeted residual stress adaptation operator, this yields our framework, MoSA. • We propose a physics-informed residual stress adaptation module, which augments an isotropic constitutive prior with bounded residual operators to capture the mild anisotropy and heterogeneity commonly present in realworld objects. By explicitly correcting the prior stress response in a structured and progressive manner, it improves simulation fidelity while preserving physical inductive bias and partial interpretability. • We introduce a motion-constrained optimization strategy that leverages motion cues from dynamic reconstruction to supervise temporal and spatial derivatives of the deformation field. This higher-order supervision provides more direct constraints than pixel reconstruction alone, improving data efficiency and reducing overfitting in video-based dynamics learning. • Experiments on both synthetic and real data demonstrate superior accuracy, generalization, and robustness, and analyses confirm that our model learns physically meaningful residual anisotropy. We further validate the practical impact of improved real-to-sim dynamics in a robot manipulation setting.

System Identification. Earlier works have explored estimating physical parameters directly from visual input (Liu et al., 2025; Chen et al., 2024; Xu et al., 2025b), aided by advances in differentiable physics simulation . These methods (Jiang et al., 2025; Vasile et al., 2025; Xu et al., 2025a) often compare rendered outputs with 2D ground truth, achieving promising results but are generally limited to elastic or rigid materials and can suffer from reconstruction artifacts. PAC-NeRF (Li et al., 2023) jointly optimizes geometry and physical parameters from multi-view videos but is hindered by rendering and geometry errors. SpringGS (Zhong et al., 2024) uses a spring-mass model for simulatable elastic objects, while GIC (Cai et al., 2024) improves physical estimation using shape constraints from dynamic 3DGS. Vid2Sim (Chen et al., 2025a) adopts a ViT for parameter estimation with a Simplicits (Modi et al., 2024) for refinement. These methods estimate parameters with a fixed physical model family, so accuracy is bounded by model misspecification. Our approach retains the isotropic backbone and learns residual dynamics beyond it.

2. Related Works Neural Physical Dynamics. Many prior works use neural networks to model 3D dynamics, often predicting fu2

Title Suppressed Due to Excessive Size

Neural-Adapted

Adaptation Progressive Adaptation

Figure 1. Overview of the pipeline. (a) Two-stage dynamic reconstruction. (b)Simulation with progressive anisotropic stress adaptation (c) Motion-constrained optimization strategy

3. Methodology

and a set of motion primitives in GIC (Cai et al., 2024). For fair comparisons, we adopt the method from GIC as our reconstruction tool.

Problem Definition. This work aims to learn real-to-sim dynamics from multi-view videos by preserving an isotropic constitutive prior as the backbone and modeling the residual effects induced by mild anisotropy and heterogeneity. Given multi-view video sequences Vi | i = 1, 2, . . . , n of a deforming object with known camera parameters, our goal is to estimate the global physical parameters of the isotropic prior and learn residual stress correction operators that compensate for the mismatch between simplified constitutive assumptions and real materials. The overall pipeline is illustrated in Fig. 1. Notably, we do not aim to model strongly anisotropic materials from scratch; instead, we target mild continuum residuals that often account for the remaining real-to-sim error when the common isotropic pipelines are adopted to model real-world dynamics.

Constitutive Model. Similar to prior works, we represent object geometry as particles and simulate using the Material Point Method (MPM) (Jiang et al., 2016). A core component of MPM is the constitutive model, which defines how objects deform under external forces by relating strain (ϵ) to stress (σ), typically as σ = σ(ϵ), where ϵ is derived from the deformation gradient F (e.g., Green strain: ϵ = 21 (FT F − I)). In elastoplastic mechanics, constitutive models include an elastic response and a plastic projection that maps excessive strain back to admissible configurations, expressed as ϵ̂ = ϵ(ϵ) or F̂ = ϵ(F ). More details are provided in the Appendix P. 3.1. A Structured, Progressive Stress Adaptation for Modeling Anisotropic Residual Effect

Dynamic Gaussian Splatting. 3DGS represents scenes explicitly by rendering a set of Gaussians using efficient differentiable rasterization. Each Gaussian is defined by its mean x0 , covariance Σ, opacity o, and color via spherical T −1 1 harmonics, with the form G(x) = e− 2 (x−x0 ) Σ (x−x0 ) . To ensure Σ is positive semi-definite, it is decomposed as Σ = RSS T RT , where R is a rotation matrix and S a scale matrix. During rendering, the color and mask are computed by blending N ordered Gaussians overlapping each pixel. X X I(u) = Ti αi ci , A(u) = Ti αi , (1) i∈N

Problem Analysis about Anisotropic Modeling. Classic constitutive models often assume materials are homogeneous and isotropic. Under these assumptions, mechanical behavior can be described by only a few parameters (e.g., Young’s modulus and Poisson’s ratio), which greatly simplifies modeling and enables stable calibration. However, this approximation is rarely exact in the real world. Even when an object is reasonably treated as near-isotropic at a coarse level, it typically exhibits mild anisotropy and spatial heterogeneity due to internal structure, fabrication, or wear. After the near-isotropic backbone is calibrated, these residual effects can dominate the remaining real-to-sim error. Among them, anisotropy is generally harder to express and identify than heterogeneity, since it depends on directional responses rather than only spatial variation. Our method therefore preserves an isotropic constitutive prior as the backbone and learns structured residual corrections for both effects, with particular focus on residual anisotropy as the more challenging component.

i∈N

Qi−1

where Ti = j=1 (1 − αj ). αi (x) is the projected 2D occupancy αi (x) = oi exp(− 12 (x − µi )T ΣTi (x − µi )). The common paradigm of dynamic 3DGS is to use a deformation model to map static Gaussian splats to other timestamps, as in: (δx, δs, δq) = Fθ (γ(x), γ(t)).

(2)

Different methods define Fθ differently, e.g., an MLP in (Yang et al., 2024), a neural triplane in (Wu et al., 2024), 3

Title Suppressed Due to Excessive Size y plane

Here we take the simple linear anisotropy as an example. Let the strain at a material point be ϵkl , a symmetric second-order tensor for k, l = 1, 2, 3. In conventional linear isotropic models, the stress is computed as σkl = Ckl ϵkl , where Ckl is a 3 × 3 stiffness matrix defined by just two independent parameters, Young’s modulus and Poisson’s ratio, yielding 2 degrees of freedom. In contrast, linear anisotropic models use tensor contraction: σij = Cijkl ϵkl , where Cijkl is a 3 × 3 × 3 × 3 fourth-order stiffness tensor. Accounting for symmetry properties (Cijkl = Cjikl and Cijkl = Cijlk ), this tensor still has 36 independent parameters, making it far more expressive—but also much more complex to model. In other words, the anisotropic stiffness tensor has 36 degrees of freedom, much higher than the isotropic case, making manual specification impractical. And this is only for linear anisotropy; the complexity increases further in nonlinear settings.

𝜎𝑦𝑥 𝜎𝑥𝑦 𝜎𝑥𝑥

3

𝑟𝑑 ෍ (𝐼𝑦𝑗 + 𝑐𝑦𝑗 )𝜎𝑗𝑦 j

𝜎𝑦𝑧 𝜎𝑧𝑦 𝜎𝑥𝑧 𝜎𝑧𝑥

Redistribute

3

𝑟𝑑 ෍ (𝐼𝑧𝑗 + 𝑐𝑧𝑗 )𝜎𝑗𝑧

𝜎𝑧𝑧 Rescale

j

3

z plane

x plane

𝑟𝑑 ෍ (𝐼𝑥𝑗 + 𝑐𝑥𝑗 )𝜎𝑗𝑥 j

z plane

x plane

Figure 2. Rescaling and redistributing of stress tensor

complementary matrix since we assume that the prior constitutive equation, even when describing isotropy, still has some deficiencies. Therefore, it is expected to make slight adjustments to the isotropic stress. We produce this term by: Lpre = W pre Φ(U, σp ) + bpre ,

(4)

where Φ is an MLP with two hidden layers. σp refers to the prior stress σp = σ(ϵ). U is the Right Polar Tensor of the deformation gradient obtained by Singular Value Decomposition (F = RU). W pre and bpre transform latent vectors into the stress domain.

Anisotropic-informed Architecture Design.To address the challenges outlined above, we do not formulate a physical equation entirely from scratch. Instead, we correct and refine a prior isotropic constitutive model to capture anisotropy and heterogeneity. Specifically, we design a physics-informed network as a projection function that maps isotropic stress to anisotropic configurations based on the current material state. The projection is defined as σij = fθ (Ckl ϵkl |s) in which s denotes the current state (including deformation gradient, prior stress, position, etc.). Moreover, rather than modeling the function f directly with a full black-box network, we define it in a specific form: σ̂ij = Tσ−1

y plane

𝜎𝑦𝑦

Second, C is the fourth-order tensor for contracting with the isotropic stress, which plays the most important role in this pipeline. Unlike the previously stated fixed elastic stiffness tensor, C is an adaptive correction term dynamically predicted by the network based on the current state of each material point. Directly predicting a full fourthorder tensor is difficult and can lead to numerical instability. To address this, we decompose C into four second-order tensors Cx , Cy , Cz , and Cxyz with a microplane-based stress redistribution constraint (Bažant Z P, 1996), and apply them sequentially to stress correction. We provide the detailed derivation of the microplane-based constrains in Appendix O. This design improves learning stability and efficiency, while providing better physical grounding. Specifically, we rescale and redistribute each component of the stress element, as shown in Fig. 2, considering only those lying on the same microplane during adjustment. That is, each corrected stress is a linear combination of all stress components on its corresponding microplane. For clarity, we use the x-plane as an example in the following.

X Φ,pre (I +CΨ )+LΘ,post , (3) ijkl )Tσ ·(σ(ϵ)kl +Lkl ij kl

where σ(ϵ)kl is the prior isotropic constitutive model, I is an identity matrix, σ̂ij is the corrected stress incorporating anisotropy. The Tσ is defined to map stress tensors from the global coordinate system to the material (or object) coordinate system. This is essential because the constitutive relations for anisotropic materials are inherently defined along their principal directions. The correction consists of three learnable terms: the redistribution term CΨ ijkl , the preΦ,pre linear term Lkl , and the post-linear term LΘ,post , each ij predicted by a cascaded neural network with parameters Ψ, Φ, and Θ, respectively. These corrections are applied progressively to adjust the prior stress. This formulation can be interpreted as first correcting the isotropic prior stress by LΦ,pre and then rescaling and redistributing it along differkl Θ,post ent directions by CΨ . Details are illustrated ijkl and Lij in the lower panel of Fig. 1 and the following section.

T

[σ̂xx , σ̂xy , σ̂xz ] = (I + Cx ) [σxx , σxy , σxz ]

T

(5)

where I is an identity matrix, Cx is a redistribution coefficient matrix associated with the x-plane. When C is set to zero, the stress degrades to the prior stress. This allows the network to focus only on learning the residual values. If we expand the σ̂xx in Eq. 5, we can obtain: x,y,z X

First, CΨ ijkl is a fourth-order tensor, and we decompose the tensor contraction into four matrix multiplications for learning efficiency. Further details will be provided later. The architecture of the cascaded network can be seen in the lower panel of Fig. 1. LΦ,pre is defined as an isotropic

(Ixj +cxj )σxj = (1+cxx )σxx +cxy σxy +cxz σxz (6)

j

The first term serves as the rescaling while the subsequent terms handle the redistribution. 4

Title Suppressed Due to Excessive Size

Similarly, we apply Cy and Cz to the y- and z-planes. Next, the mutual equivalence of shear stresses theorem (σyx = σxy , σzx = σxz ) is considered to merge paired shear stress components, i.e. σ̂xy = 0.5(σ̂xy + σ̂yx ). We then reuse Eq. 5 with another redistribution tensor Cxyz to further redistribute the updated normal stress components T σ̂xx , σ̂yy , σ̂zz , that is (I + Cxyz ) [σxx , σyy , σzz ] . This step enhances the modeling of anisotropic effects along the three principal directions. Since the refinement at each microplane occurs at the same time, we use a network to jointly predict Cx , Cy and Cz (Eq. 7).

This approach yields a smoother and more physically consistent representation, avoiding the discontinuities and abrupt jumps that may arise from particle-wise parameterization. To further stabilize the learning process, we apply a 0.2 · tanh(·) activation on the network output, constraining it to the range [−0.2, 0.2], and encourage the network predictions to have zero mean and minimal variance. This design can not only capture the global material parameter while allowing controlled local variations. 3.2. Motion-constrained Optimization Strategy Most prior works supervise dynamics learning using rendering losses that align rendered images and masks with the input videos. Such losses provide only indirect and coarse constraints on the deformation field, mainly through silhouettes and boundaries. Since image formation is manyto-one, different 3D deformations and material responses can lead to similar rendered appearances, making the optimization ill-conditioned. This motivates introducing additional motion-aware constraints that provide more direct supervision on the underlying deformation.

l1 = Ψ1 (U, σp + Lpre ), Cx , Cy , Cz = α1 · tanh(W c l1 + bc ),

(7)

Cxyz , l2 = α2 · tanh(Ψ2 ([l1 , σ̂xx , σ̂yy , σ̂zz ])). Ψ1 and Ψ2 are two MLPs. α1 and α2 are two empirical coefficients that govern the effect of the prior stress; the larger α, the more deviation from the prior stress is allowed. We constantly set the two αs as 0.1; the reason is presented in the Appendix D. Applying these Cs separately is equivalent to performing a tensor contraction with C, which we derive in the Appendix O. Similar to Lpre , LΘ,post is another linear adjustment applied after the C is contracted. Lpost = W post Θ(l2 , U, σ̂) + bpost .

Dynamic 3D sequences reconstructed from multi-view videos provide informative motion cues for learning object dynamics. We leverage these cues to better constrain dynamics optimization. Following (Cai et al., 2024), we reconstruct dynamic scenes from multi-view videos with dynamic 3DGS. To ensure accurate initial geometry and improve the capture of local deformations, we adopt a twostage reconstruction procedure that decouples static geometry from dynamic deformation.

(8)

in which l is the latent code from Eq. 7. σ̂ is the corrected stress from the previous step. Since the stress tensor should be symmetric, both Lpre and Lpost are symmetrized as L = 12 (L + LT ). To ensure all C matrices are invertible, we predict an upper triangular matrix Uc and a lower triangular matrix Lc for each C, and construct C by multiplying them. To guarantee invertibility, we add a small diagonal offset to both Lc and Uc .

We first reconstruct a static 3DGS at the initial timestamp and apply the filling strategy in (Cai et al., 2024) to obtain the initial particle representation at t = 0, denoted as P0 . At this stage, we model Gaussian primitives as isotropic, with equal scales and identity rotations, which provides a clean and stable initialization of the geometry. We then fix P0 as the point scaffold and train a dynamic 3DGS on top of it. During dynamic training, we allow the Gaussian primitives to have anisotropic covariances, so that they can deform and rotate freely over time through Fθ (Eq. 2). This decoupling prevents the dynamic model from absorbing static geometry errors and encourages Fθ to focus on learning motion and deformation. Therefore, we further introduce higher-order supervision that directly constrains the temporal and spatial derivatives of the deformation field using motion cues from the dynamic reconstruction.

By the same logic, this paradigm can also be applied to correct the plastic yield projection. Λ̂ij = Tσ−1

X Φ,pre (I + CΨ ) + LΘ,post , (9) ij )Tσ (Σ(ϵ̃) + Lj i kl

where Σ(ϵ̃) represents the singular value vector of the deviatoric principle strain tensor. The corrected version can be interpreted as an equivalent deviatoric principal strain tensor that conforms to the new anisotropic plasticity criterion. Heterogeneity Modeling. To model material heterogeneity, instead of assigning an independent set of material parameters to each particle (as was done in previous works (Dagli et al., 2025; Xu et al., 2025b)), we estimate a single global material parameter that characterizes the overall behavior of the object, together with a continuous implicit field that provides local variations of the global parameter at each spatial location within the object. The physical parameters for position x are defined as P(x) = Pglobal ·(1+0.2·tanh(η(x))).

We begin with the analysis in (Xie et al., 2024) to bridge the deformation of Gaussian splats to the time-varying deformation gradients in simulation. if a Taylor series expansion of the deformation field ϕ̃p is performed, ϕ̃p (X, t) = xp + Fp (X −Xp ) can be obtained, in which Fp is the deformation gradient at particle p. Substituting the formula into Gaus-

5

Title Suppressed Due to Excessive Size T

−1 −1 −1 1 sians: Gp (x) = e− 2 (ϕ (x,t)−Xp ) Ap (ϕ (x,t)−Xp ) , a F -dependent covariance matrix Ap,t = Fp,t Ap,0 (Fp,t )T is obtained. Combining all into the tensor format can yield:

At = Ft A0 FTt .

rotational components can be regularized as follows: Lr = Et (||δRt − Ut VtT ||2F ),

(10)

where δRt is derived from δq in Eq. 2. Upt and VptT are extracted from SVD of F on particle p at time t. ||.||F represents the Frobenius norm. Further discussion can be seen in Appendix J.

This formulation allows 3DGS rendering to be expressed as a function of both simulated particle positions and their corresponding deformation gradients. This forms the foundation of our derivative-based regularization strategy. Therefore, the rendering loss can be written as: n 1X Limg/mask = L1 (Rimg/mask (Pt , Ft ), It /Mt ), n t (11) where R∗ refers to rendering functions, Pt , Ft are simulation particle states, It , Mt are image and mask groundtruth.

As stated in Heterogeneous modeling of Sec. 3.1, we encourage the output of the local material implicit field to have zero mean and minimal variance via: 2

Lhet = λµ E[η(x)] 2   + λvar max 0, Var[η(x)] − 0.52 .

Ψ, Θ, Φ, P, η = arg

min

Ψ,Θ,Φ,P,η

Ltotal .

(16)

4. Experiments Dataset. We thoroughly evaluate our method on both synthetic and real-world datasets. For synthetic evaluation, we use the dataset released by PAC-NeRF (Li et al., 2023), which provides multi-view RGB sequences together with 3D ground-truth geometry trajectories across diverse materials. This synthetic benchmark allows us to quantitatively evaluate 3D trajectory accuracy, prove the dynamics grounding, and to verify that our residual stress corrections do not introduce spurious anisotropy when an isotropic backbone is sufficient. Current research still relies heavily on synthetic data due to the scarcity of real-world multi-view dynamic datasets. To fill this gap, we collect a real-world multi-view dataset using an advanced light-field capture system with synchronized industrial cameras. The dataset contains 7 objects, including four elastic objects (Mandarin, Chick1, Chick2, Peanut), two elastoplastic objects (Rainbowball, Rabbit), and one plastic object (Gorilla). Objects are released from random heights to fall vertically onto a platform, and 12 cameras capture synchronized surround-view RGB sequences. More details are provided in Appendix F.

where Rgf (·, ·) represents the Gaussian flow maps rendered from the Gaussian splats updated by the simulation results. xt = x0 + δxt , At is the covariance matrix derived from δst and δqt from Eq. 2. Then we introduce how we guide the spatial derivative of deformations. As stated above, the updated At can be obtained from the dynamic reconstruction results. A0 is the initial covariance and fixed for all timestamps. Intuitively, if we can solve F from this, we can build point-wise supervision for the deformation gradient. This equation can be interpreted as a similarity transformation from A0 to At . Due to the high degrees of freedom, deriving an explicit expression for F is challenging. Additionally, since the solution to a similarity transformation is not unique, the scale of the resulting F may differ. To solve this, we supervise the rotational and scaling components of F separately. For the scaling component, we focus only on the relative scale, disregarding the absolute scale. To sum up, we decompose F into Λ and U, V via SVD, and define the scaling and rotation regularization respectively. δSt ) − norm(Λ − I)). S0 + ϵ

(15)

This loss enforces an unbiased overall distribution. The global parameter governs the average response, and the local field contributes only subtle refinements. The total loss is: Ltotal = Lf low + β1 (Limg + Lmsk ) + β2 (Lr + Ls ) + β3 Lhet .The final optimization target is:

The derivatives of the deformation field to time and spatial coordinates are essentially the particle velocities v = ∂X ∂t and particle deformation gradients F = ∂X ∂x . As for the former, we utilize Gaussian Flow (Gao et al., 2024) to obtain the flow map for each camera view, and use it to supervise the flow obtained by simulation, as in Eq. 12. In this context, the Gaussian flow acts as an implicit supervision for particle velocity. This not only results in a sharper silhouette, serving as an implicit shape constraint, but more importantly, allows for precise supervision within the silhouette. n
1X Lf low = L1 Rgf (Pt , Ft ), (Rgf (xt , At )) , (12) n t

Ls = Et (norm(

(14)

Baselines and Metrics. We compare our method with PACNeRF (Li et al., 2023), DEL (Wang et al., 2024a), GIC (Cai et al., 2024), NeuMA (Cao et al., 2025), and Vid2Sim (Chen et al., 2025a). PAC-NeRF, GIC, and Vid2Sim estimate global physical parameters under classic physical models, while DEL and NeuMA learn dynamics operators with neural components. We evaluate on two settings: object dynamics grounding on the PAC-NeRF dataset and initial-state generalization on our real-world dataset. We report Chamfer

(13)

This regularization constrains the relative magnitude of scaling around the vicinity of particle p in the simulation. The 6

Title Suppressed Due to Excessive Size Table 1. Dynamic Grounding on PAC-NeRF dataset. All metrics are scaled by 100 for clarification. torus

cat

playdoh

droplet

Cream

Bird

Letter

Mean

Methods CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ CD↓ EMD↓ 21.8 21.7 20.2 20.1

11.6 10.7 9.9 9.8

9.8 7.9 7.6 7.3

14.4 12.8 12.6 12.4

18.6 12.2 12.3 11.4

5.6 2.5 2.5 2.3

10.4 9.8 10.2 9.6

3.2 1.7 1.9 1.5

20.5 19.8 19.5 19.4

12.7 9.8 10.1 9.5

19.3 17.8 16.5 16.3

21.1 20.2 19.5 19.2

12.8 12.6 10.3 10.1

8.5 7.2 7.5 6.6

16.2 14.5 13.8 13.5

11.0 9.3 9.1 8.8

Vid2Sim

GIC

NeuMA

DEL

PAC DEL GIC Ours

Ours

Figure 4. Comparisons of long-term dynamics (a) and sensitivity of physical laws (b).

GT

4.2. Real-world Generalization Real-world evaluation is essential because, even when an isotropic backbone is a reasonable approximation, common objects still exhibit mild anisotropy and heterogeneity that can dominate the remaining real-to-sim error. In this experiment, we learn dynamics from video recordings and evaluate generalization on newly captured sequences with different initial orientations and drop heights. We render the simulated results into each camera view and compare them with the held-out footage to measure how well each method transfers beyond its training data.

Figure 3. Qualitative Comparisons between our methods and baselines. More Visualization can be seen in the Appendix.

Distance (CD) and Earth Mover Distance (EMD) for 3D geometry error, and PSNR and SSIM for rendering quality. For PSNR and SSIM, we crop images around the object to avoid background dominance. On the real-world dataset, we report only PSNR and SSIM since ground-truth 3D particle trajectories are unavailable.

Quantitative results are shown in Table 2, and qualitative comparisons appear in Fig. 3. Additional visualizations and videos are available in the Appendix and Sup. Mat. Our method consistently outperforms all baselines across material types in real-world settings. Fully neural approaches such as DEL perform worst, likely because modeling full state transitions with large networks overfits under limited real data. NeuMA improves a pretrained NCLaw prior via LoRA, but degrades when the pretraining distribution differs from our real materials; moreover, both NCLaw and NeuMA assume isotropy, limiting their ability to capture anisotropic behaviors. Among system identification methods, GIC performs best, likely due to its use of geometry and motion cues from dynamic 3DGS. Vid2Sim performs poorly, likely because its forward predictor is trained only on synthetic data and does not generalize to our real-world scenes. Overall, our residual stress adaptation remains strong on real data, indicating a high performance ceiling for real-to-sim dynamics learning.

Implementation Details. We use the same loss weights for all experiments: β1 = 1.0, β2 = 0.1, β3 = 0.1, λµ = 1.0, and λvar = 1.0. Here, the image and flow losses provide the primary supervision, while the deformation regularizers and the heterogeneity constraint serve as auxiliary guidance. 4.1. Results on Synthetic Data Object Dynamics Grounding. We evaluate all methods on the object dynamics grounding task using the PAC-NeRF dataset. The goal is to assess how accurately each approach captures object dynamics under multi-view video supervision. Following common practice, all methods first reconstruct the initial object shape at the first timestamp using 3DGS, and then simulate the deformation with their corresponding dynamics modules. Quantitative results are reported in Table 1, where our method achieves the best overall performance. 7

Title Suppressed Due to Excessive Size Table 2. Quantitative comparisons of the initial state generalizations on the real-world dataset. All SSIM is scaled by 100. Chick1 Gorilla Mandarin Chick2 Peanut Rabbit RBball Mean Methods PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ DEL 28.32 91.6 29.11 90.4 31.92 92.4 28.59 90.9 28.38 91.3 28.57 91.4 30.95 91.7 29.41 91.4 NeuMA 30.73 92.4 29.78 91.1 31.85 92.4 28.92 91.0 30.70 91.8 28.30 92.3 30.74 91.1 30.00 91.7 GIC 30.93 92.5 29.75 91.0 29.54 91.6 28.17 90.9 32.88 91.3 28.08 91.3 30.78 91.7 30.02 91.5 Vid2Sim 26.71 90.9 29.02 90.0 25.85 89.8 25.70 89.5 30.69 94.6 26.85 90.3 31.71 91.7 28.08 91.0 Ours 32.05 92.7 30.19 91.9 32.83 92.7 30.17 91.6 33.01 92.0 30.35 92.1 32.06 92.7 31.35 92.3 Table 3. Ablation studies on three real-world scenes

Ours

Prior

Ours

Prior

Rabbit Gorilla Rainbowball PSNR↑ SSIM↑ PSNR↑ SSIM↑ PSNR↑ SSIM↑ no C 28.42 91.3 29.66 91.0 30.67 91.8 no Het 29.67 91.5 29.98 91.2 31.15 91.8 no Lpre 29.97 91.7 30.12 91.3 31.18 91.9 no Lpost 29.88 91.7 30.18 91.6 31.21 91.6 no Lf low 29.13 91.6 29.83 90.9 30.29 91.0 no Lscale 30.05 91.9 30.12 91.1 30.84 92.1 no Lrot 30.14 91.3 29.99 91.2 31.98 92.0 Full model 30.35 92.1 30.19 91.9 32.06 92.7

4.3. Additional Experiments and Analysis Ablation Studies. We conduct ablations on three real-world examples (rabbit, gorilla, and rainbowball) to assess the contribution of each component. Results are summarized in Table 3. We remove key terms in the stress adapter, including C, Lpre , Lpost , and the heterogeneity modeling (no Het). The redistribution term C has the largest impact, confirming the importance of anisotropic stress redistribution, while the linear correction terms also provide consistent gains. The heterogeneity module is also important. Removing it leads to consistent performance drops across all three scenes: PSNR decreases from 30.35 to 29.67 on rabbit, from 30.19 to 29.98 on gorilla, and from 32.06 to 31.15 on rainbowball. This shows that controlled spatial variation complements residual anisotropic stress adaptation and improves real-world dynamics modeling.

Figure 5. Applications of zero-shot robot manipulation transfer.

parameters to compensate for model errors. In contrast, our method refines the constitutive prior itself, making it less dependent on model selection. GIC exhibits large performance variation across different priors, whereas our method remains stable. Notably, even starting from a simple linear stress–strain model, our approach can progressively correct the prior and achieve higher PSNR, demonstrating strong robustness and flexibility.

Predictions of Long-term Dynamics. We present comparisons of long-term simulation in Fig. 4 (a), which include the SSIM at each timestep for the ”rabbit” scenario. The comparison demonstrates that our method consistently outperforms the others throughout the entire simulation. Additionally, our method exhibits a much slower increase in loss over time compared to the other approaches, maintaining better stability in long-term predictions.

Application on Downstream Manipulation Tasks. We further evaluate whether learning more accurate real-world dynamics benefits downstream robot manipulation. Specifically, we define two manipulation tasks and first learn their object dynamics from video interactions. We then deploy the learned dynamics in simulation to train a policy, and directly transfer the policy to the real robot for evaluation. As a baseline, we use a purely isotropic prior physical model in the simulator. As shown in Fig. 5, our learneddynamics policies achieve substantially higher real-world success rates: on the elastic-rabbit placement task (placing a deformable rabbit onto a white box), our policy succeeds in 68 trials per 100, compared to 42% for the isotropic prior; on the tower-hanging task, our policy succeeds in 82 trials, compared to 55% for the isotropic baseline. These results suggest that accounting for continuum residual effects can be a key factor for further improving real-to-sim dynamics and sim-to-real transfer.

Robustness to Prior Physical Models. We study how the choice of prior physical law affects performance on Chick1. As shown in Fig. 4(b), parameter-calibration baselines are sensitive to the selected model, since they can only adjust

Sensitivity to Dynamic Reconstruction Quality. Our motion constraints rely on dynamic reconstruction cues, so we further evaluate their sensitivity to reconstruction quality. Specifically, we train MoSA using dynamic 3DGS recon-

We further ablate the motion constraints, including the flow loss and the scale and rotation losses. All these terms improve accuracy and generalization, confirming the importance of derivative-level supervision for stable learning from videos.

8

Title Suppressed Due to Excessive Size

Figure 7. Spatial visualization of the normalized learned heterogeneity field η(x) on test objects. Reddish regions correspond to larger learned local stiffness, while bluish regions indicate softer regions. Figure 6. Analysis of learned anisotropy with directional Jacobian

visualize the normalized learned η(x) on three held-out test objects in Fig. 7. Since η(x) modulates the global material parameter locally, its spatial distribution reflects learned material heterogeneity. The resulting fields show smooth and object-dependent patterns, with reddish regions corresponding to higher local stiffness. This confirms that our continuous field learns physically meaningful spatial variations, rather than merely fitting unstructured residual noise.

Table 4. Sensitivity to dynamic reconstruction quality on the Rabbit scene. R.Steps denotes the progress of reconstruction.

R.Steps

PSNR↑

SSIM↑

∆PSNR

5k 7k 10k 15k 20k 30k

29.30 29.65 29.95 30.15 30.28 30.35

91.6 91.7 91.8 91.9 92.0 92.1

+0.17 +0.52 +0.82 +1.02 +1.15 +1.22

More experimental analysis including system identification, future prediction, etc., are reported in our Appendix.

5. Conclusion

structions obtained at different optimization stages on the Rabbit scene, and report the downstream simulation performance rather than the reconstruction quality itself. As shown in Table 4, motion constraints consistently improve performance over the model without motion constraints, whose PSNR is 29.13. Even early-stage reconstructions provide positive gains, indicating that the proposed motion constraints benefit downstream physics learning instead of simply amplifying reconstruction errors.

In conclusion, our results indicate that further reducing the real-to-sim gap for everyday deformable objects is less about re-learning the dominant isotropic backbone model, and more about capturing the mild but systematic continuum residual effects caused by the anisotropy and heterogeneity that standard simulators ignore. By explicitly constraining learning to residual stress corrections and strengthening supervision with motion cues, our approach achieves consistently better stability and transfer in both dynamics metrics and downstream robot manipulation. We hope this work encourages future dynamics learning methods to combine strong physical backbones with targeted residual modeling for scalable real-world deployment.

Further Analysis of What the Model Learns. We analyze whether our residual stress adaption module learns physically meaningful effects rather than overfitting a single sequence. Since PAC-NeRF is generated by an isotropic simulator, it contains no true anisotropy; thus our model should not introduce spurious directionality. We define the directional Jacobian response of Eq. 3: w(θ) = ∂ σ̂ ∂σ iso · vec(n ⊗ n) 2 , which characterizes the stress redistribution trend induced by the Progressive Stress Adaptation module. Fig. 6 shows that on PAC-NeRF, w(θ) stays near zero with no clear directional pattern, indicating only mild adjustments to the isotropic prior. In contrast, on real data (Mandarin), w(θ) exhibits pronounced directionality aligned with the object axis, suggesting that stress adaptation module adaptively captures missing anisotropic effects when present. More analysis is provided in Appendix A.

Impact Statement MoSA improves the fidelity of real-to-sim continuum dynamics by modeling mild residual anisotropy and heterogeneity beyond standard isotropic assumptions. This can benefit a broad range of applications that rely on reliable physical simulation, including graphics, embodied interaction, and robotics, by enabling more accurate and data-efficient dynamics identification and reducing failures caused by model mismatch. In the long term, stronger simulation reliability may support safer design, testing, and deployment of systems that interact with the physical world,

Spatial Variation of the Learned Heterogeneity. We 9

Title Suppressed Due to Excessive Size

and may facilitate higher-fidelity physics-consistent digital twins.

Fu, Z., Mao, H., and Yin, B. Inverse design of lattice metamaterials for fully anisotropic elastic constants: A datadriven and gradient-based method. Composite Structures, 359:118975, 2025.

As with many advances in physical modeling, improved simulation can be used in both beneficial and potentially harmful ways, including in autonomous systems. Our work focuses on foundational algorithmic improvements for dynamics learning and does not involve sensitive personal data or human subjects. We encourage responsible use and evaluation in safety-critical settings.

Gao, Q., Xu, Q., Cao, Z., Mildenhall, B., Ma, W., Chen, L., Tang, D., and Neumann, U. Gaussianflow: Splatting gaussian dynamics for 4d content creation. arXiv preprint arXiv:2403.12365, 2024. Jaques, M., Burke, M., and Hospedales, T. Physics-asinverse-graphics: Unsupervised physical parameter estimation from video. arXiv preprint arXiv:1905.11169, 2019.

References Bažant Z P, Xiang Y, P. P. C. Microplane model for concrete. i: Stress-strain boundaries and finite strain. Journal of Engineering Mechanics, 122:245–254, 1996.

Jiang, C., Schroeder, C., Teran, J., Stomakhin, A., and Selle, A. The material point method for simulating continuum materials. In Acm siggraph 2016 courses, pp. 1–52. ACM, 2016.

Bolliger, D., Fadini, G., Bambach, M., and Rupenyan, A. Differentiable material point method for the control of deformable objects. arXiv preprint arXiv:2512.13214, 2025.

Jiang, H., Hsu, H.-Y., Zhang, K., Yu, H.-N., Wang, S., and Li, Y. Phystwin: Physics-informed reconstruction and simulation of deformable objects from videos. arXiv preprint arXiv:2503.17973, 2025.

Cai, J., Yang, Y., Yuan, W., He, Y., Dong, Z., Bo, L., Cheng, H., and Chen, Q. Gic: Gaussian-informed continuum for physical property identification and simulation. arXiv preprint arXiv:2406.14927, 2024.

Jing, C., Bandi, J. K., Ye, J., Duan, Y., Abbeel, P., Wang, X., and Yi, S. Contact-aware neural dynamics. arXiv preprint arXiv:2601.12796, 2026.

Cao, J., Guan, S., Ge, Y., Li, W., Yang, X., and Ma, C. Neuma: Neural material adaptor for visual grounding of intrinsic dynamics. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2025.

Kerbl, B., Kopanas, G., Leimkühler, T., and Drettakis, G. 3d gaussian splatting for real-time radiance field rendering. ACM Transactions on Graphics, 42(4), July 2023. URL https://repo-sam.inria.fr/ fungraph/3d-gaussian-splatting/.

Chen, C., Dou, Z., Wang, C., Huang, Y., Chen, A., Feng, Q., Gu, J., and Liu, L. Vid2sim: Generalizable, video-based reconstruction of appearance, geometry and physics for mesh-free simulation. In Proceedings of the Computer Vision and Pattern Recognition Conference, pp. 26545– 26555, 2025a.

Li, X., Qiao, Y.-L., Chen, P. Y., Jatavallabhula, K. M., Lin, M., Jiang, C., and Gan, C. Pac-nerf: Physics augmented continuum neural radiance fields for geometry-agnostic system identification. arXiv preprint arXiv:2303.05512, 2023.

Chen, G., Suri, S., Wu, Y., Voulga, E., Levin, D. I., and Pai, D. K. Learning simulatable models of cloth with spatially-varying constitutive properties. arXiv preprint arXiv:2507.21288, 2025b.

Li, X., Qiao, Y., Chen, Y., Du, C., Li, Z., and Zhang, C. A machine learning-based data-driven approach for modelling anisotropic and tension-compression asymmetry behavior of elastoplastic materials using limited anisotropic experiment data. European Journal of Mechanics-A/Solids, pp. 105733, 2025.

Chen, Q., Gao, T., Jiang, C., and Wu, T. Eqcollide: Equivariant and collision-aware deformable objects neural simulator. arXiv preprint arXiv:2506.05797, 2025c.

Lin, Y., Lin, C., Xu, J., and Mu, Y. Omniphysgs: 3d constitutive gaussians for general physics-based dynamics generation. arXiv preprint arXiv:2501.18982, 2025.

Chen, S., Xu, Y., Yu, C., Li, L., and Hsu, D. Differentiable particles for general-purpose deformable object manipulation. arXiv preprint arXiv:2405.01044, 2024.

Liu, Z., Ye, W., Luximon, Y., Wan, P., and Zhang, D. Unleashing the potential of multi-modal foundation models and video diffusion for 4d dynamic physical scene simulation. In Proceedings of the Computer Vision and Pattern Recognition Conference, pp. 11016–11025, 2025.

Dagli, R., Xiang, D., Modi, V., Loop, C., Tsang, C. F., Chen, A. H., Hu, A., State, G., Levin, D. I., and Shugrina, M. Vomp: Predicting volumetric mechanical property fields. arXiv preprint arXiv:2510.22975, 2025. 10

Title Suppressed Due to Excessive Size

Ma, P., Chen, P. Y., Deng, B., Tenenbaum, J. B., Du, T., Gan, C., and Matusik, W. Learning neural constitutive laws from motion observations for generalizable pde dynamics. In International Conference on Machine Learning, pp. 23279–23300. PMLR, 2023.

with lmms. In Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV), pp. 7231– 7240, October 2025b. Xue, H., Torralba, A., Tenenbaum, J., Yamins, D., Li, Y., and Tung, H.-Y. 3d-intphys: towards more generalized 3d-grounded visual intuitive physics under challenging scenes. Advances in Neural Information Processing Systems, 36, 2024.

Mittal, H., Zhuang, P., Lee, H.-Y., and Tulsiani, S. Uniphy: Learning a unified constitutive model for inverse physics simulation. In Proceedings of the Computer Vision and Pattern Recognition Conference, pp. 16208–16218, 2025. Modi, V., Sharp, N., Perel, O., Sueda, S., and Levin, D. I. Simplicits: Mesh-free, geometry-agnostic elastic simulation. ACM Transactions on Graphics (TOG), 43(4):1–11, 2024.

Yang, Z., Gao, X., Zhou, W., Jiao, S., Zhang, Y., and Jin, X. Deformable 3d gaussians for high-fidelity monocular dynamic scene reconstruction. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 20331–20341, 2024.

Vasile, F., Qiu, R.-Z., Natale, L., and Wang, X. Gaussianaugmented physics simulation and system identification with complex colliders. arXiv preprint arXiv:2511.06846, 2025.

Zhang, K., Li, B., Hauser, K., and Li, Y. Adaptigraph: Material-adaptive graph-based neural dynamics for robotic manipulation. arXiv preprint arXiv:2407.07889, 2024.

Wang, J., Jingkai, S., Zhang, Z., He, J., Zhang, Q., Sun, M., and Xu, R. Del: Discrete element learner for learning 3d particle dynamics with neural rendering. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024a.

Zhao, Y., Chen, H., Liu, C., Li, Z., Herrmann, C., Hur, J., Li, Y., Yang, M.-H., Raj, B., and Xu, M. Toward material-agnostic system identification from videos. In Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 5944–5956, 2025.

Wang, J., Zhang, Z., He, J., and Xu, R. PFGS: High fidelity point cloud rendering via feature splatting. In Proceedings of the European Conference on Computer Vision (ECCV), pp. 193–209, 2024b.

Zhobro, M., Geist, A. R., and Martius, G. Learning 3dgaussian simulators from rgb videos. arXiv preprint arXiv:2503.24009, 2025. Zhong, L., Yu, H.-X., Wu, J., and Li, Y. Reconstruction and simulation of elastic objects with spring-mass 3d gaussians. In European Conference on Computer Vision, pp. 407–423. Springer, 2024.

Wei, L. and Freris, N. M. Multi-scale graph neural network for physics-informed fluid simulation. The Visual Computer, 41(2):1171–1181, 2025. Wu, G., Yi, T., Fang, J., Xie, L., Zhang, X., Wei, W., Liu, W., Tian, Q., and Wang, X. 4d gaussian splatting for real-time dynamic scene rendering. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 20310–20320, 2024. Xie, T., Zong, Z., Qiu, Y., Li, X., Feng, Y., Yang, Y., and Jiang, C. Physgaussian: Physics-integrated 3d gaussians for generative dynamics. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4389–4398, 2024. Xu, Q., Liu, J., Yu, S., Wang, Y., Zhou, Y., Zhou, J., Cui, J., Ong, Y.-S., and Zhang, H. Neuspring: Neural spring fields for reconstruction and simulation of deformable objects from videos. arXiv preprint arXiv:2511.08310, 2025a. Xu, X., Ge, W., Qiu, D., Chen, Z., Yan, D., Liu, Z., Zhao, H., Zhao, H., Zhang, S., Liang, J., and Chen, Y.-C. Gaussianproperty: Integrating physical properties to 3d gaussians 11

Title Suppressed Due to Excessive Size

Figure 8. Physically grounded analysis of learned anisotropy

A. Verification of Learned Anisotropy Knowledge To verify that the proposed physics-informed Progressive Stress Adaptation module truly learns the underlying physical anisotropy rather than merely benefiting from increased parameters, we design two validation experiments based on the same simulation. A uniaxial compression test is conducted on a cylindrical specimen with a prescribed strong XY-plane anisotropy in its elastic modulus. During loading, deformation videos are recorded to evaluate the learned mechanical behavior. First, we analyze the Jacobian matrix in Eq. 3. Physically, this Jacobian characterizes the internal stress redistribution trend learned by the stress adaptation module. By scanning directions on the XY plane, we compute the directional ∂ σ̂ , and compare its pattern with the GT anisotropic modulus EGT (θ), as Jacobian response w(θ) = ∂σ iso · vec(n ⊗ n) 2 shown in Fig. 8(a). The two curves exhibit that the learned redistribution trend closely aligns with the physical stiffness anisotropy. Next, we derive the effective directional modulus from the model-predicted stress–strain relations and compare it with two baselines: (i) an isotropic model and (ii) a neural model that directly learns the anisotropic mapping, i.e., σ̂ = N N (σ iso , U, θ). As shown in Fig. 8 (b), our model (orange line) accurately reproduces the GT anisotropic modulus (black line), while the baselines (grey and green lines) fail to capture the directional dependence. These results confirm that the superior performance of the stress adapter primarily arises from its physics-aware design rather than from a mere increase in network parameters.

B. Evaluation on Spring-GS dataset To further indicate the superiority of our method, we follow the setting in Spring-GS to use 20 frames as training data to predict the rest 10 frames. We first train these models the training data, then use the trained models or estimated parameters to simulate the full trajectory. We report the results in Table 5. It is observed that our method stably performs better than other counterparts. The performance of NeuMA may not be as expected because the data used for its pretraining has a significant gap compared to the Spring-GS dataset. Table 5. Quantitative Comparisons of the Future State Prediction on Spring-GS dataset. CD↓ Method SpGS NeuMA GIC Ours torus 2.38 1.68 0.75 0.68 cross 1.57 2.12 1.09 1.02 cream 2.22 1.01 0.94 0.93 apple 1.87 0.16 0.22 0.15 paste 7.03 6.22 2.79 2.61 chess 2.59 2.12 0.77 0.83 banana 18.50 0.58 0.12 0.11

EMD↓ SpGS NeuMA GIC Ours 0.087 0.043 0.034 0.032 0.051 0.059 0.058 0.049 0.094 0.050 0.050 0.046 0.076 0.029 0.030 0.028 0.126 0.097 0.096 0.094 0.095 0.113 0.059 0.056 0.135 0.093 0.017 0.015

12

PSNR↑ SpGS NeuMA GIC Ours 16.83 17.83 20.24 20.75 16.93 23.52 30.51 30.76 15.42 18.96 19.15 19.22 21.55 26.35 26.89 27.16 14.71 16.81 16.31 16.95 16.08 16.37 18.44 19.15 17.89 22.08 29.29 29.39

SSIM↑ SpGS NeuMA GIC Ours 0.919 0.913 0.942 0.949 0.940 0.928 0.939 0.941 0.862 0.877 0.909 0.912 0.902 0.927 0.948 0.950 0.872 0.896 0.894 0.903 0.881 0.838 0.912 0.916 0.904 0.913 0.964 0.968

Title Suppressed Due to Excessive Size

C. Network Architecture As we stated in the main text, we adopt four lightweight neural networks to separately model the three neural corrected modules. Φ and Θ networks aim to produce the two linear adjustment term. Ψ1 and Ψ2 networks generate the four redistribution matrix which are decomposed from the fourth-order tensor contraction C. All subnetworks include two hidden layers, each with 64 hidden dimensions. The U in each input is the Right Polar matrix of the deformation gradient F, which is a symmetric matrix. Hence we only input the upper triangle of the U (6 independent components) to these networks. We adopt the singular value decomposition to compute the U . Similarly, we include the prior stress into the input tensor, which is also symmetric. The output of each network is a 12-dimensional vector, where the first six elements and the last six elements form an upper triangular matrix and a lower triangular matrix, respectively. These two matrices are then each added to an identity matrix and multiplied together to obtain the final output. According to LU decomposition, this ensures that the output remains an invertible matrix with desirable properties.

D. Parameter Search for the α In this section, we discuss the selection of hyperparameters α1 and α2 from the Equation 7 in our main paper. To guide our selection, we performed a coarse parameter search with a step size of 0.1 in the Chick1 scenario, as shown in the Fig. 9. The heatmap illustrates the variation in PSNR values based on different combinations of α1 and α2 , where darker colors correspond to higher PSNR values. From this exploration, we observed that the optimal settings for both α1 and α2 lie around 0.15, which is where we chose to set them. Even with these relatively coarse settings, the PSNR remains high. Notably, even when the values of both α1 and α2 are set to 0.5, the PSNR still reaches 31.29, which is higher than the baseline method. This suggests that our method is not overly sensitive to small changes in the values of α1 and α2 , and reasonably chosen settings still yield competitive results.

Figure 9. Grid Search for the hyperparameter α1,2

E. Effect of the Prior Stress We evaluate the effect of removing prior stress, specifically σ(ϵ)kl in Eq. 3 of the main paper, and relying only on C and L (as shown in Table 3). While the ”no prior” model fits the training data well, its test performance decreases, indicating overfitting. In contrast, incorporating prior stress improves generalization by facilitating optimization, as evidenced by the higher PSNR and SSIM values in both training and testing for our method. Notably, even when prior stress is omitted, our method still outperforms isotropic models like GIC and DEL (shown in Table 3 of the main paper). This highlights the importance of anisotropic modeling in capturing more accurate and generalizable dynamics. The results in Table 6 demonstrate that, while prior stress enhances model performance, our method’s ability to capture complex dynamics without it still surpasses simpler isotropic models. 13

Title Suppressed Due to Excessive Size Table 6. Comparisons between with and without prior stress Episode Method Train Test

Chick1 PSNR↑ SSIM↑ NoPriorStress 33.15 0.937 Ours 33.81 0.938 NoPriorStress 31.03 0.924 Ours 32.05 0.927

Chick2 PSNR↑ SSIM↑ 30.82 0.921 30.85 0.922 28.99 0.911 30.17 0.916

F. Real-World Dataset F.1. Data Collection

(a) Light field reconstruction system

(b) Release objects from predetermined height

Figure 10. Visual illustration of our real-world data collection pipeline

(a) Chick

(b) Rabbit

(c) Mandarin

(d) Chick2

(e) Peanut

(f) Gorilla

(g) RainbowBall

Figure 11. The visualizations of all objects in our real-world dataset

To achieve high-quality, real-world dynamic data, we utilize an advanced light field reconstruction system, as shown in Fig.10, to gather multi-view RGB sequences. This system is equipped with ultra-high-precision industrial cameras arranged in a spherical configuration, allowing for synchronized, surround-view data collection with a delay of less than 1 ms. During the data collection phase, we release the objects to fall vertically from a predetermined height onto a platform, with the light field reconstruction system capturing synchronized multi-view data. After capturing the footage, we conduct a meticulous frame-by-frame inspection to ensure there are no flips or occlusions during the fall. We employ 12 to 15 industrial cameras, each capturing a sequence of 15 RGB frames from different fixed viewpoints. This approach ensures the acquisition of high-quality, realistic, and complete free-fall dynamics. After collecting the initial data, we proceeded with post-processing the raw images. The images are cropped to remove redundant static elements from the background. This ensured that the focus remained solely on the dynamic aspects of the captured data. After that, we employed GroundingSAM techniques to generate precise masks for each frame, facilitating the accurate separation of the foreground and background. F.2. Objects Description Given that the dataset is specifically designed to capture real-world dynamics during free fall, we carefully select seven distinct objects, each with unique physical properties, to represent a wide range of materials, as illustrated in Fig. 11. These objects are: 14

Title Suppressed Due to Excessive Size

Gorilla: This object contains deformable materials such as space sand, which undergo irreversible plastic deformations during free fall. Unlike elastic objects will return to their original shape after deformation, Gorilla exhibits permanent, non-recoverable changes. We select this object to demonstrate the superior performance of our model in learning the physical dynamics of plastic material only from visual data. Chick1, Chick2, Peanut, and Mandarin: These objects are elastic but differ in their levels of softness and elasticity and vary in their textures and shapes. By including a variety of elastic materials with distinct properties, we ensure that our dataset captures a broad spectrum of deformation behaviors, which is essential for understanding diverse physical dynamics. Rainbowball, Rabbit: These objects are made of elastoplastic materials, and exhibit both elastic and plastic behaviors. When subjected to forces during free fall, these objects initially deform elastically, returning to their original shape when the stress is small. However, they undergo plastic deformations that are irreversible when the force is beyond a certain threshold, similar to purely plastic materials. Rainbowball and Rabbit represent a more complex material behavior that bridges the gap between purely elastic and purely plastic materials.

G. Effect and Ablation of the Prior Stress We evaluate the effect of removing prior stress, specifically σ(ϵ)kl in Eq. 3 of the main paper, and relying only on C and L (as shown in Table 3). While the ”no prior” model fits the training data well, its test performance decreases, indicating overfitting. In contrast, incorporating prior stress improves generalization by facilitating optimization, as evidenced by the higher PSNR and SSIM values in both training and testing for our method. Notably, even when prior stress is omitted, our method still outperforms isotropic models like GIC and DEL (shown in Table 3 of the main paper). This highlights the importance of anisotropic modeling in capturing more accurate and generalizable dynamics. The results in the following Table demonstrate that, while prior stress enhances model performance, our method’s ability to capture complex dynamics without it still surpasses simpler isotropic models. Table 7. Comparisons between with and without prior stress Episode Method Train Test

Chick1 PSNR↑ SSIM↑ NoPriorStress 33.15 0.937 Ours 33.81 0.938 NoPriorStress 31.03 0.924 Ours 32.05 0.927

Chick2 PSNR↑ SSIM↑ 30.82 0.921 30.85 0.922 28.99 0.911 30.17 0.916

H. Complexity Analysis with Other Learned Simulators We conducted a complex analysis comparing the training and inference times of our method with two other learning-based methods, NeuMA and DEL. As shown in Table 8, our method outperforms the others in terms of both training and inference efficiency. Specifically, our method reduces training time to 3.5 hours, significantly faster than NeuMA (6.5 hours) and DEL (9.0 hours). In terms of inference, our method also leads with only 0.76 seconds per inference, compared to 1.05 seconds for NeuMA and 1.09 seconds for DEL. These results demonstrate that our approach not only achieves superior performance but also operates with greater efficiency, making it a more scalable solution for real-time applications. Table 8. Comparisons of Training and Inference Times.

NeuMA DEL NoMotion NoStress NoFlow Ours Train (h) 6.5h 9.0h 4.5h 5.8h 4.3h 3.5h Infer (s) 1.05s 1.09s – – – 0.76s

I. Results of System Identification Although our method is designed to both estimate global physical parameters and refine local constitutive behavior through stress correction, it still achieves superior performance even when evaluating only the accuracy of global parameter estimation. This demonstrates that our approach is not just capable of correcting local dynamics but also highly effective at identifying global material properties. Out of 23 parameter estimations across different scenes, our method achieved the 15

Title Suppressed Due to Excessive Size

highest accuracy in 21 cases. This improvement can be attributed to two key components in our design. First, the stress correction module simplifies the learning objective by offloading local modeling errors from the global parameter estimation. By refining the prior constitutive model at the local level, it allows the global parameters—such as Young’s modulus and Poisson’s ratio—to focus on capturing the overall material behavior without being forced to compensate for local discrepancies. Second, our motion-constrained optimization strategy plays a crucial role. By incorporating motion cues from dynamic 3DGS, this strategy not only implicitly supervises particle positions (first-order dynamics), but also regularizes the spatial and temporal derivatives of the deformation field. These richer constraints help guide the optimization toward more physically plausible solutions, leading to more accurate parameter estimation. Overall, the combination of these two components allows our method to outperform baselines, even in standard system identification settings. Table 9. System Identification Performance on PAC-NeRF Dataset. All the values and quantities in the table are scaled based on the magnitude of the ground truth values.

Droplet Letter Cream Toothpaste Torus Bird Playdoh Cat Trophy

PAC-NeRF µ, κ = 2.09, 1.085 µ, κ = 83.85, 1.35 µ, κ = 1.21, 1.57 τY , η = 3.16, 5.6 µ, κ = 6.51, 2.22 τY , η = 228, 9.77 E, ν = 1.04, 0.322 E, ν = 2.78, 0.273 E, ν, τY = 3.84, 0.272, 1.69 E, ν, τY = 1.61, 0.293, 3.57 θf0 ric = 36.1◦

GIC µ, κ = 2.01, 0.18 µ, κ = 95.05, 1.00 µ, κ = 1.03, 1.48 τY , η = 2.98, 6.6 µ, κ = 4.19, 9.24 τY , η = 226, 9.1 E, ν = 0.99, 0.295 E, ν = 3.08, 0.284 E, ν, τY = 1.58, 0.322, 1.56 E, ν, τY = 0.98, 0.296, 3.76 θf0 ric = 38.0◦

Ours µ, κ = 2.01, 1.07 µ, κ = 97.00, 1.00 µ, κ = 1.01, 1.39 τY , η = 2.97, 7.75 µ, κ = 4.52, 9.37 τY , η = 212, 9.63 E, ν = 0.99, 0.298 E, ν = 3.02, 0.29 E, ν, τY = 1.72, 0.289, 1.55 E, ν, τY = 0.98, 0.297, 3.77 θf0 ric = 38.5◦

Ground Truth µ, κ = 2, 1 µ, κ = 100, 105 µ, κ = 1, 1 τY , η = 3, 10 µ, κ = 5, 10 τY , η = 200, 10 E, ν = 1, 0.3 E, ν = 3, 0.3 E, ν, τY = 2, 0.3, 1.54 E, ν, τY = 1, 0.3, 3.85 θf0 ric = 40◦

Figure 12. Comparisons of simulation results on PAC-NeRF dataset with or without metric regularization (e.g. Lscale , Lrot )

J. Explanation of Scale Regularization In this section, we establish the relationship between the scale changes of Gaussian splats and the simulated deformation gradient. Since the static Gaussians are assumed to be isotropic, all anisotropic deformations during dynamic reconstruction are captured by the deformation network Fθ in dynamic 3DGS paradigm introduced in the main paper, Eq. 2. Therefore, our goal is to derive a way to regularize the deformation gradient using the output of this network. Let Λ denote the singular values of the deformation gradient F, and let δS be the scale change predicted by the deformation network. Based on this, we derive the following: ΛS0 := S0 + δS δS (Λ − I) := S0

16

(17)

Title Suppressed Due to Excessive Size Table 10. Quantitative comparison on the public real-world subset of Spring-GS. SSIM values are scaled by 100. PSNR↑

Method Spring-GS NeuMA GIC Ours

SSIM↑

Bun

Burger

Dog

Pig

Bun

Burger

Dog

Pig

30.69 31.27 34.68 41.12

34.01 23.78 40.45 43.25

32.10 25.61 37.17 39.21

34.97 25.40 38.32 42.16

99.2 99.4 99.5 99.7

99.4 99.3 99.7 99.8

99.4 99.5 99.7 99.8

99.6 99.5 99.7 99.9

Since we only supervise the relative deformation magnitudes along three directions, we incorporate a normalization function to ensure the scale of the deformation is consistent across these directions.

norm(Λ − I) := norm(

δS ) S0 + ϵ

(18)

K. Qualitative Ablations of the two Deformation Regularization Losses We present the results of the ablation about the effect of the two deformation regularizations, e.g., the scale and rotation loss, in Fig. 12, comparing two cases: one with the two losses and the other without, observed from two different viewpoints. We can observe that the ring-like shape with metric regularization looks more stable, with a uniform and smooth surface. In contrast, the shape without metric regularization shows significant deformation. The inclusion of metric regularization results in a smoother and more consistent geometric shape, reducing distortion and irregularity. This leads to better structural stability and consistency, which is particularly beneficial for modeling or optimization tasks. Therefore, metric regularization proves to be effective in enhancing model performance, especially when dealing with complex geometric forms.

L. Limitation Our current pipeline relies on accurate multi-view capture, calibration, and dynamic reconstruction; failures in these steps can weaken the motion cues and reduce training stability. In addition, MoSA is designed to model mild residual anisotropy and heterogeneity beyond an isotropic backbone; cases with extremely strong anisotropy or rapidly varying material properties may require richer priors or additional sensing signals. As future work, we plan to improve robustness to imperfect observations, extend the framework to broader material regimes, and explore incorporating complementary sensors or learned uncertainty to better handle challenging real-world captures.

M. Additional Qualitative Comparisons We show additional qualitative comparisons for the other scenarios in the real-world dataset in Fig. 13, Fig. 14, Fig. 15, Fig. 16, Fig. 17, Fig. 18. The images illustrate the effectiveness of our method, as it consistently produces clear, accurate dynamics over time for different objects. Our approach captures the object behavior in a detailed and realistic manner. The fine-grained motion details highlight the robustness and high performance of the proposed method in dynamic learning from real-world data. Vid2Sim performs worse than reported in its original paper, mainly due to several limitations. It is trained on synthetic data and heavily relies on the accuracy of the initial LBS weight predictions. Additionally, it requires strict input formats, such as exactly 16 video frames. In contrast, GIC, as a general-purpose optimization-based method, outperforms Vid2Sim on our real-world dataset. DEL and NeuMA, trained with only a single sequence, show limited generalization during testing. DEL even occasionally fails to maintain basic geometric structure.

N. Evaluation on the Spring-GS Real-World Subset To further validate the real-world generalization ability of MoSA, we conduct additional experiments on the public real-world subset of Spring-GS. As shown in Table 10, our method achieves the best performance on all scenes in terms of both PSNR and SSIM, demonstrating that the advantage of MoSA also transfers to public real-world benchmarks. 17

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Title Suppressed Due to Excessive Size

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Figure 13. Qualitative comparison of different baselines and our methods on the chick2 scenario.

Figure 14. Qualitative comparison of different baselines and our methods on the chick1 scenario.

18

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Title Suppressed Due to Excessive Size

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Figure 15. Qualitative comparison of different baselines and our methods on the rabbit scenario.

Figure 16. Qualitative comparison of different baselines and our methods on the peanut scenario.

O. Derivation of Microplane Parametrization of the Fourth-Order Correction Operator In this appendix, we show, in a self-contained way, that the proposed microplane parametrization using (Cx , Cy , Cz , Cxyz ) is mathematically equivalent to using a single fourth-order tensor Cijkl acting on the stress tensor in Eq. (3). The key idea is that all three stages of our Constructions are linear maps on the stress components, and any composition of linear maps is itself a linear map. In Voigt notatio,n this means that the The whole pipeline can always be written as a single 6 × 6 matrix, which in turn is equivalent to a fourth-order tensor. We show that our microplane parametrization with (Cx , Cy , Cz , Cxyz ) is equivalent to a single fourth-order correction tensor as in Eq. 3. 19

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Title Suppressed Due to Excessive Size

Figure 17. Qualitative comparison of different baselines and our methods on the gorilla scenario.

Fourth-order tensor contraction in Voigt form. The core correction in Eq. 3 can be written as σ̂ij = σij + Cijkl σkl ,

(19)

where Cijkl is a fourth-order tensor. Because σij is symmetric, we can collect the six independent components into a Voigt vector     σxx σ̂xx σyy  σ̂yy       σzz   σ̂zz     . σ= σ̂ =  (20) ,   σyz   σ̂yz  σzx  σ̂zx  σxy σ̂xy Then Cijkl is in one-to-one correspondence with a 6 × 6 matrix C, and Eq. 19 is exactly the matrix–vector multiplication σ̂ = (I6 + C) σ,

(21)

where each row of (I6 + C) lists the coefficients of the linear combination that produces one corrected stress component from all six original components. In other words, Eq. 21 is just a linear map R6 → R6 in Voigt space. Redistribution within each microplane. In Eqs. 5 and 6, we first redistribute stresses within each microplane. We group the components belonging to the same plane as simple 3-vectors, for example       σyy σzz σxx σ x = σxy  , σ y = σyz  , σ z = σzx  . (22) σyx σzy σxz On each plane, we apply a learned 3 × 3 linear map: σ̃ x = (I3 + Cx ) σ x ; σ̃ y = (I3 + Cy ) σ y ; σ̃ z = (I3 + Cz ) σ z ,

20

(23)

GT

Ours

Vid2Sim

GIC

NeuMA

DEL

Title Suppressed Due to Excessive Size

Figure 18. Qualitative comparison of different baselines and our methods on the Rainbowball scenario.

which is exactly what Eq. 5 does. Each entry of σ̃ x , for example, is just a linear combination of (σxx , σxy , σxz ). Stacking all intermediate components back into a Voigt vector σ (1) (and using only index reordering), this whole “within-plane” redistribution is again a 6 × 6 matrix acting on σ: σ (1) = Bplane σ,

Bplane ∈ R6×6 .

(24)

Thus Eqs. 5-6 already define a valid linear map in Voigt form. Redistribution between microplanes. Eq. 7 then mixes only the normal stresses from different planes, while keeping the shear components unchanged. Writing  (1)    σxx σ̂xx  (1)  n(1) = σyy n̂ = σ̂yy  , (25) , (1) σ̂ zz σzz Eq. 7 can be written as n̂ = (I3 + Cxyz ) n(1) ,

(26) (1)

(1)

(1)

which is again a simple 3 × 3 linear map. The shear components are just copied over: σ̂yz = σyz , σ̂zx = σzx , σ̂xy = σxy . We deliberately
apply cross-plane redistribution only to the normal stresses but not to the shear components. Normal stresses σxx , σyy , σzz control volumetric response and opening/closure of microplanes, and their interaction across directions (e.g., through Poisson-like effects) is physically meaningful. In contrast, each shear component (e.g., σxy , σyz , σzx ) represents 21

Title Suppressed Due to Excessive Size

sliding on a specific microplane and is already fully mixed within that plane by the per-plane maps Cx , Cy , Cz and the subsequent shear symmetrization. Introducing an additional cross-plane mixing layer for the shears would couple fundamentally different sliding modes (e.g., mixing σxy with σyz ), which is hard to justify physically and mainly increases the number of parameters without a clear benefit. Therefore, we restrict Cxyz to act only on the normal block in Voigt space, leading to a clean block structure of the global 6 × 6 redistribution matrix. Putting this together in Voigt form gives another 6 × 6 matrix Bnorm such that σ̂ = Bnorm σ (1) ,

Bnorm ∈ R6×6 ,

(27)

where the upper 3 × 3 block of Bnorm is I3 + Cxyz and the lower block is the identity. Putting everything together and the global redistribution matrix. Combining the “within-plane” redistribution Eq. 24, the fixed shear symmetrization (a matrix S with entries 0, 1, 12 ), and the “between-plane” redistribution Eq. 27, we obtain the overall Voigt-form multiplication σ̂ = Bnorm S Bplane σ = A σ ∈ R6×6 . Equivalently, we can write this global redistribution explicitly as     σ̂xx σxx σyy  σ̂yy        σ̂zz    = A6×6  σzz  .   σyz   σ̂yz      σzx  σ̂zx  σ̂xy σxy

(28)

(29)

Each entry aαβ is the final coefficient of σβ in the linear combination that produces σ̂α ; It is obtained by “filling” contributions from the within-plane maps Cx , Cy , Cz and the between-plane map Cxyz through the product Bnorm SBplane . The key point is simple: Bplane , S, and Bnorm are all linear operators, and their product is still a linear operator A. Therefore, there exists a C that satisfies A = I6 + C,

⇒

σ̂ = (I6 + C) σ,

(30)

which has exactly the same form as (21) and thus corresponds to an equivalent fourth-order tensor Cijkl . Intuitively, we do not change the fact that each new stress component is a linear combination of all old stress components. We simply factor the large matrix (I6×6 + C) into three steps: first redistributing stresses within each microplane (controlled by Cx , Cy , Cz ), then symmetrizing shear pairs, and finally redistributing the normal components across microplanes (controlled by Cxyz ). This keeps the mapping linear and physically reasonable, while turning the otherwise opaque fourth-order tensor Cijkl into a sequence of more interpretable operations.

P. Constitutive Models as Physical Priors A constitutive model explains how a material reacts to stress, strain, or other external forces. It defines the material’s behavior by connecting stress and strain through equations, which can describe complex behaviors like elasticity, plasticity, and fracture. The MPM simulation can model a wide range of materials by using different constitutive models. Our approach also relies on a constitutive model as a prior, which is then used to perform adaptive redistribution and correction of the stress based on this model. For any method requiring a constitutive model, we select the appropriate model from those listed in this chapter using a large language model, and it remains the same across different models. Generally, Cauchy stress σ is represented as 1 ∂Ψ E E T (F )F (31) σ= J ∂F in which J = det(F ), Ψ is the energy density function to describe the behavior of a certain material, which can be seen as a function of the elastic component of F (F E ). The F E is projected by F E = ψ(F ), which is the so-called return mapping function. Both ψ and Ψ are referred to as constitutive models. We define the Cauchy stress σ for these models as follows: 22

Title Suppressed Due to Excessive Size

Linear Cauchy Stress. The simplest practical constitutive model is linear elasticity, defined in terms of the small strain tensor, or the infinitesimal strain tensor. This strain tensor will give rise to a computationally lightweight constitutive model. Jσ = µ(F + F T − 2I)F T + λ tr(F − I)F T

(32)

Neo-Hookean Cauchy Stress. NeoHookean elastic model is a widely used nonlinear hyperelastic framework, for simulating material elasticity and predicting deformations, which can be expressed as:
Jσ = µ F F T − I + λ log(J)I (33) Fixed Corotated Cauchy Stress. Another simple and widely used model that is defined from the Singular Value Decomposition (SVD) is the so-called fixed corotated model. Jσ = 2µ (F − R) F T + λ(J − 1)JI

(34)

where R = U V T and F = U ΛV T , which represent the singular value decomposition of the elastic deformation gradient. J is the determinant of F . St. VK Cauchy Stress. St. Venant-Kirchhoff model offers significant benefits over a linear elastic model, which is a rotationally invariant model. Jσ = U (2µϵ + λ sum(ϵ)1) V T F T (35) in which ϵ = log(Λ). Other symbols remain the same meaning with previously defined. In our study, the deformation gradient is multiplicatively decomposed as F = F E F p , based on a yield stress condition in all plasticity models. We also can obtain the plastic projection process as F E = ψ(F ). Then a hyperelastic constitutive model is applied to compute the Cauchy stress σ. For a purely elastic continuum, we set F E = F = ψ(F ). Drucker-Prager Mapping. The return mapping of Drucker-Prager plasticity can be defined as : F E = U Z(Σ)V T  P  if (ϵ) > 0,  1, P Z(Σ) = Σ,   if δγ ≤ 0, (ϵ) ≤ 0,   exp ϵ − δγ ϵ̂ , otherwise. ∥ϵ̂∥ q   P 2 sin φ (ϵ) , α = 23 · 3−sin φff , and ϕf is the friction angle. where δγ = ∥ϵ̂∥ + α (dλ+2µ) 2µ Von Mises Mapping. Von Mises shares the same mapping function as the Drucker-Prager (36). ( Σ,   δγ ≤ 0, Z(Σ) = ϵ̂ exp ϵ − δγ ∥ϵ̂∥ , otherwise.

(36)

(37)

(38)

Y Here δγ = ∥ϵ̂∥F − τ2µ , τY is the yield stress.

Herschel-Bulkley Mapping. We implement Continuum foam and we follow the simple version of that described in trial PhysGaussian, s can be recovered as s = s · ∥sstrial ∥ , !  r  η 2 trial trial s=s − s − σY / 1 + (39) 3 2µ∆t The corresponding Kirchhoff stress can be: τ= in which bE = FE FE

i h
1 κ 2 J − 1 I + µ dev det(bE )− 3 bE 2

T

23

(40)

Title Suppressed Due to Excessive Size

Q. Overview of Material Point Method The Material Point Method (MPM) is a numerical simulation method for solving continuum mechanics problems. MPM uses particles and a background grid as discrete elements of the simulation domain, enabling the modeling of a wide range of materials. It solves problems by transferring mass and momentum between the particles and the background grids during simulation and performing computations on the grid. It is gonvered by the momentum and mass conservation in the Eulerian form: dρ = −ρ∇ · v, ρa = ∇ · σ + ρb. (41) dt where ρ denotes the density, v refers to the velocity, a is the acceleration. σ represents the Cauchy stress tensor and b corresponds to body forces. Mass conservation is naturally preserved in MPM by transporting the Lagrangian particles. In practical MPM implementations, the integrals are computed using particle-based quadrature, where each material point contributes to the integrals based on its mass and position. The test functions w are typically chosen from a set of shape functions defined on the computational background mesh. Spatial Discretization To numerically solve the weak form derived earlier, the Material Point Method (MPM) performs spatial discretization by leveraging shape functions defined over the computational background grid. These shape functions serve to transfer information between the Lagrangian material points and the Eulerian background mesh, thus enabling a hybrid discretization approach. Specifically, the test functions w and the acceleration field a are approximated using a set of basis functions {Na } defined on the grid nodes. This results in the following semi-discrete equation for each grid node a: NG X

Mab a(b) = −

NP X

V0 (i)τ (i)∇Na (x(i))+

i=1

b=1

NP X

(42) M (i)Na (x(i))b(i),

i=1

where NG is the total number of grid nodes and NP is the number of material points. The left-hand side of the equation involves the mass matrix Mab , which accounts for the contribution of particle mass to the grid system. It is defined as: Mab =

NP X

M (i) · Na (x(i)) · Nb (x(i)),

(43)

i=1

where M (i) is the mass of particle i, and Na , Nb are the grid basis functions evaluated at the current position x(i) of the material point. The discretization scheme effectively transfers the continuous balance of momentum equation into a system of equations posed on the grid. MPM uses this system to compute the nodal accelerations, which are then interpolated back to the particles during the update step. Temporal Discretization To advance the solution in time, the MPM typically employs an explicit time integration scheme due to its simplicity and ease of implementation. In this formulation, we adopt the explicit forward Euler method, which approximates time derivatives using finite differences over a discrete time step ∆t. Starting from the semi-discrete system obtained from spatial discretization, the temporal discretization updates the nodal velocity v(b) at each grid node b from time t to t + 1 as follows: NG X b=1

Mab

NP X vt+1 (b) − vt (b) =− V0 (i)τt (i)∇Na (xt (i)) ∆t i=1

+

NP X

(44)

M (i)Na (xt (i))bt (i),

i=1

In practice, the mass matrix Mab is often diagonalized (lumped) to decouple the equations for each node, allowing for a simple explicit velocity update: ∆t vt+1 (a) = vt (a) + fnet (a), Maa 24

Title Suppressed Due to Excessive Size

where fnet (a) includes both internal and external force contributions. The explicit Euler scheme is conditionally stable, requiring small enough ∆t to satisfy a CFL-like condition that depends on wave speed and grid spacing.

25