Add 2D and 3D stretching elasto-plasticity (yield) for interpolated flexes (dof="trilinear" and dof="quadratic").

### Physical Model
- **Cell / Face Plastic Engineering Strain (`d->plastic`):**
  - **3D Volume Interpolation (`interp > 0`, `elastic2d="none"`):** Each grid cell $e$ stores a 6-component symmetric plastic engineering strain tensor $\varepsilon^p\_e = (\varepsilon^p\_{xx}, \varepsilon^p\_{yy}, \varepsilon^p\_{zz}, \gamma^p\_{xy}, \gamma^p\_{xz}, \gamma^p\_{yz})$ (`6 * cx * cy * cz` entries per flex, allocated for the full `cellcount` grid).
  - **2D Shell Interpolation (`interp < 0`, `dim="2"` with `elastic2d="stretch"` or `"both"`):** Each outer grid face $f$ stores a 3-component in-plane plastic engineering strain tensor $\varepsilon^p\_f = (\varepsilon^p\_{11}, \varepsilon^p\_{22}, \gamma^p\_{12})$ (`6 * (cy*cz + cx*cz + cx*cy)` entries per flex, i.e. 3 per outer face). With `elastic2d="bend"`, `yield` is ignored (bending plasticity is not supported for interpolated flexes) and the flex remains purely elastic.
- **Corotated Polar Decomposition on the Physical Jacobian:**
  - In `mj_flexPassiveInterp` (and `mujoco_warp`'s `compute_interp_cell_quat` / `compute_interp_face_quat`), the cell/face rotation $R\_e$ is the standard polar decomposition of the physical Jacobian $J\_e = \partial x / \partial X$ at the cell/face center, obtained by scaling the columns of the parametric gradient $\partial x / \partial \xi$ by the reference cell extents $(1/\Delta x, 1/\Delta y, 1/\Delta z)$ (and the face tangents $t\_1 / \Delta\_0$, $t\_2 / \Delta\_1$ before the normal cross product in 2D). It does not read `d->plastic`: since the plastic strain is additive in the corotated frame, $R\_e$ is the exact rigid rotation at any rotated plastic rest configuration, including non-cubic cells.
  - For large rotations (non-positive trace of the Jacobian), the iterative polar decomposition (`mju_flexPolarRotation` / `mat33_to_quat_polar`) is warm-started from the Gram-Schmidt rotation of the Jacobian columns, so it does not stall at the identity saddle point (e.g. 180 degree rotations, which are exact for cells with consistent physical scaling).
  - Using the precomputed cell/face stiffness matrix $K\_e$, the nodal displacement vector in the local corotated frame is offset by the plastic deformation $\Delta d\_i = R\_e^T x\_i - X\_i - \varepsilon^p\_e (X\_i - X\_0)$, leaving $K\_e$ invariant while shifting the zero-force configuration.
- **Return Mapping (`advanceFinish` / `_next_plastic_flex_interp`):**
  - At the end of each step, the cell-center corotated total strain $\varepsilon^{\text{tot}} = \mathrm{sym}(R\_e^T J\_e - I)$ is decomposed into elastic and plastic parts $\varepsilon^{\text{el}}\_e = \varepsilon^{\text{tot}} - \varepsilon^p\_e$, and radial return projects $\varepsilon^{\text{el}}\_e$ onto the von Mises ($J\_2$) yield surface (3D isochoric in 3D, plane-stress with out-of-plane thinning in tension in 2D) with threshold $\varepsilon\_y =$ `flex_yield[f]`.
  - To prevent degenerate cell Jacobians under extreme repeated loading, $\varepsilon^p\_e$ is bounded in Frobenius norm ($\|\varepsilon^p\_e\|\_F \le 0.35$, guaranteeing $\lambda\_{\min}(I + \varepsilon^p\_e) > 0.71$).
- **Dimensionless Keyframe Scaling:** Because $\varepsilon^p\_e$ is a dimensionless engineering strain, `ScaleKeyframe` leaves interpolated flex entries invariant under spatial scaling (`flex_interp[f] != 0`).

PiperOrigin-RevId: 1010658301
Change-Id: I41a16150e9dfa56610794d27669e43d32599b6ca

Author: quagla

doc/XMLreference.rst

@@ -4876,15 +4876,21 @@ and **quadratic**.
 .. _flex-elasticity-yield:
 
 :at:`yield`: :at-val:`real(1), "0"`
-   Dimensionless plastic yield strain threshold :math:`\varepsilon_y` for non-interpolated 2D and 3D flexes. If zero,
-   the flex is purely elastic. If positive, each element's elastic stretching strain tensor is bounded by the von Mises
-   (:math:`J_2`) yield criterion (3D isochoric in 3D, plane-stress with out-of-plane thinning in tension in 2D) via
-   radial return mapping, and for 2D shells with bending enabled (:ref:`elastic2d<flex-elasticity-elastic2d>` **bend**
-   or **both**), each initially flat interior hinge's outer-fiber bending strain
-   :math:`\frac{t}{2} |\Delta\kappa_e| \le \varepsilon_y` is bounded via return mapping in dihedral angle space. The
-   plastic state in ``mjData.plastic`` starting at ``flex_plasticadr[f]`` stores ``3 * elemnum`` (2D) or
-   ``6 * elemnum`` (3D) stretching entries when stretching is enabled, followed by ``edgenum`` bending entries when 2D
-   bending is enabled (with boundary edges remaining ``0``).
+   Dimensionless plastic yield strain threshold :math:`\varepsilon_y` for 2D and 3D flexes. If zero, the flex is purely
+   elastic. If positive, each element's (or interpolated cell's/face's) elastic stretching strain tensor is bounded by
+   the von Mises (:math:`J_2`) yield criterion (3D isochoric in 3D, plane-stress with out-of-plane thinning in tension
+   in 2D) via radial return mapping, and for non-interpolated 2D shells with bending enabled
+   (:ref:`elastic2d<flex-elasticity-elastic2d>` **bend** or **both**), each interior hinge's outer-fiber bending strain
+   :math:`\frac{t}{2} |\Delta\kappa_e| \le \varepsilon_y` is bounded via return mapping in dihedral angle space
+   (assuming an initially flat rest configuration). In ``mjData.plastic`` at ``flex_plasticadr[f]``, a non-interpolated
+   flex stores its stretching plastic metric offsets first (``3 * elemnum`` entries in 2D or ``6 * elemnum`` entries in
+   3D when stretching is active), followed by ``edgenum`` bending plastic entries in 2D when bending is active. An
+   interpolated flex with :ref:`cellcount<body-flexcomp-cellcount>` :math:`(c_x, c_y, c_z)` stores one dimensionless
+   plastic stretching strain per grid cell (``6 * cx * cy * cz`` entries in 3D) or per outer grid face in 2D shell mode
+   (``3 * 2 * (cy*cz + cx*cz + cx*cy)`` entries); the elastic strain is measured in the corotated frame of the polar
+   decomposition of the physical deformation gradient :math:`J_e = \partial x / \partial X`. Bending plasticity is not
+   supported for interpolated flexes: with :ref:`elastic2d<flex-elasticity-elastic2d>` **bend**, :at:`yield` is
+   ignored and the flex remains purely elastic.
 
 .. _flex-elasticity-elastic2d:
 

doc/changelog.rst

@@ -25,7 +25,7 @@ Significant new features
    :width: 40%
    :target: https://github.com/google-deepmind/mujoco/blob/main/model/flex/clay.xml
 
-- Added continuum elasto-plasticity (:ref:`yield<flex-elasticity-yield>`) for non-interpolated 2D and 3D
+- Added continuum elasto-plasticity (:ref:`yield<flex-elasticity-yield>`) for 2D and 3D
   :ref:`flexes<deformable-flex>`, with volumetric-deviatoric (von Mises :math:`J_2`) stretching and 2D shell bending
   return mapping in ``mjData.plastic``.
 

mjx/mujoco/mjx/third_party/mujoco_warp/_src/constraint.py

@@ -1643,6 +1643,7 @@ def kernel(
       cell_quat = face_quat_in[worldid, cell_idx]
     else:
       cell_quat = support.compute_interp_cell_quat(
+        flex_node0,
         flexnode_xpos_in,
         1,
         ci,

mjx/mujoco/mjx/third_party/mujoco_warp/_src/derivative.py

@@ -2563,6 +2563,7 @@ def kernel(
     flex_nodeadr: wp.array[int],
     flex_stiffnessadr: wp.array[int],
     flex_nodebodyid: wp.array[int],
+    flex_node0: wp.array[wp.vec3],
     flex_stiffness: wp.array[float],
     flex_damping: wp.array[float],
     flex_centered: wp.array[bool],
@@ -2621,7 +2622,7 @@ def kernel(
     if shell_mode:
       cell_quat = face_quat_in[worldid, flex_faceadr[f] + ci]
     else:
-      cell_quat = support.compute_interp_cell_quat(flexnode_xpos_in, order_abs, ci, cj, ck, cy, cz, ny_g, nz_g, nstart, worldid)
+      cell_quat = support.compute_interp_cell_quat(flex_node0, flexnode_xpos_in, order_abs, ci, cj, ck, cy, cz, ny_g, nz_g, nstart, worldid)
     R = wp.quat_to_matrix(cell_quat)
     RT = wp.transpose(R)
 
@@ -2900,6 +2901,7 @@ def kernel(
     flex_nodeadr: wp.array[int],
     flex_stiffnessadr: wp.array[int],
     flex_nodebodyid: wp.array[int],
+    flex_node0: wp.array[wp.vec3],
     flex_stiffness: wp.array[float],
     flex_damping: wp.array[float],
     flex_faceadr: wp.array[int],
@@ -2969,7 +2971,7 @@ def kernel(
     if shell_mode:
       cell_quat = face_quat_in[worldid, flex_faceadr[f] + ci]
     else:
-      cell_quat = support.compute_interp_cell_quat(flexnode_xpos_in, order_abs, ci, cj, ck, cy, cz, ny_g, nz_g, nstart, worldid)
+      cell_quat = support.compute_interp_cell_quat(flex_node0, flexnode_xpos_in, order_abs, ci, cj, ck, cy, cz, ny_g, nz_g, nstart, worldid)
     cell_quat_inv = wp.quat(-cell_quat[0], -cell_quat[1], -cell_quat[2], cell_quat[3])
 
     op1 = order_abs + 1
@@ -3475,6 +3477,7 @@ def eff_build(m: Model, d: Data):
           m.flex_nodeadr,
           m.flex_stiffnessadr,
           m.flex_nodebodyid,
+          m.flex_node0,
           m.flex_stiffness,
           m.flex_damping,
           m.flex_centered,
@@ -3608,6 +3611,7 @@ def eff_shift(m: Model, d: Data):
         m.flex_nodeadr,
         m.flex_stiffnessadr,
         m.flex_nodebodyid,
+        m.flex_node0,
         m.flex_stiffness,
         m.flex_damping,
         m.flex_faceadr,
@@ -3963,6 +3967,7 @@ def eff_mul_m(
           m.flex_nodeadr,
           m.flex_stiffnessadr,
           m.flex_nodebodyid,
+          m.flex_node0,
           m.flex_stiffness,
           m.flex_damping,
           m.flex_faceadr,