Allow separate stretching and bending yield strains for flex plasticity.

The flex `<elasticity yield="...">` attribute (and the same on
`<flexcomp><elasticity>`) now takes one or two values, `yield="eps_stretch
eps_bend"`. For non-interpolated 2D flexes a single value applies to both. For
3D and interpolated flexes a single value is only the stretching yield strain.
Either way, existing models compile and simulate as before.

API:

-   `mjsFlex::stressyield` is now `double[2]` (stretch, bend), and
    `mjModel::flex_yield` is `nflex x 2`. Both are unreleased.
-   New `mjModel::flex_plasticmode` (`nflex x 1`) records which modes have
    plastic state: bit 1 is bending, bit 2 is stretching (same encoding as
    `elastic2d`).
-   Whether a mode is plastic is decided at compile time, as for joints and
    tendons. A mode gets plastic state, and its bit set, only if its yield
    strain is positive. Stretching state comes first, then bending state.
    For example, `yield="0.01 0"` keeps bending purely elastic, and its yield
    strain is never read at runtime.
-   Setting the yield strain of a plastic mode to 0 at runtime still means
    full plastic flow, as before.
-   A positive bending yield strain on a 3D or interpolated flex is a compile
    error ("bending yield requires a non-interpolated 2d flex"). 1D flexes
    still reject any positive value.

Why a new field: whether bending state exists could already be inferred (the
2D shell thickness is stored in `stiffness[6]` only when bending yields).
Nothing in the model indicated whether stretching state exists. A per-flex
bitmask was the smallest addition that makes both blocks and their offsets
explicit. It also replaces the per-mode layout formulas that were duplicated
in `mj_resetPlastic`, model validation, and `mj_printModel`
(which would otherwise read past a flex's state).

XML:

-   `yield` is `double[1..2]` with `reading=custom` in mjcf.schema.
-   For `<flex>` the reader resolves a single value from `dim` and `node`. For
    `<flexcomp>` it is resolved in `Make()`, once dim and dof are final (dim
    can come from a mesh or gmsh file).
-   The writer saves a single value when it reproduces the bending yield
    strain, and two values otherwise.

Engine, C and MJWarp: the return mapping and the passive forces locate
stretching and bending state by `flex_plasticmode`, not by the yield value.
In MJWarp, `flex_yield` is `wp.vec2` and `flex_plasticmode` is threaded into
the plastic and flex passive kernels. MJX warp shims were regenerated.

Compiler: `mjCFlex::PlasticStretchNum()` gives the size of the stretching block
of a compiled or uncompiled flex. `PlasticNum()` and keyframe scaling in a
scaled frame (`ScaleKeyframe`) both use it, so a bend-only shell scales all of
its entries as bending entries.

Generated files regenerated: XML tables/xsd/dmcontrol schema, references.h,
introspect structs.py, wasm and Unity bindings, mjx/warp, and golden
flex_model.data (only `flex_yield` gains a column and `flex_plasticmode` is
added). dm_control's MjModel indexer gains the new field.

Tests:

-   XML reader: 1 vs 2 values on 2D and 3D, on flexcomp and deformable flex;
    single values on 3D and interpolated flexes resolve to stretching only.
    Errors for 3 values, a negative bending yield strain, a bending yield
    strain on 3D, interpolated, and 1D flexes.
-   XML writer: round trip with 1 and 2 values, 2D and 3D.
-   Engine:
    -   stretch-only and bend-only yielding are independent;
    -   a mode with a zero yield strain has no state and stays elastic even
        if its yield strain is changed at runtime;
    -   a runtime yield strain of 0 on a plastic mode gives full plastic flow;
    -   absent modes allocate no state.
-   Keyframe scaling: a bend-only shell in a scaled frame, followed by another
    plastic flex or last.
-   Model validation of `flex_plasticmode`.
-   MJWarp: C-vs-MJWarp parity for separate yields, per-mode allocation, and
    runtime-zero full flow.

PiperOrigin-RevId: 1013971100
Change-Id: I26055b2b73b4e673119f0f28cfe4e757a128a38c

Author: quagla

doc/XMLreference.rst

@@ -3995,7 +3995,9 @@ element is used to adjust the properties of all edges in the flex.
 
 |body/flexcomp/elasticity attrib list|
    Same as in :ref:`flex/elasticity<flex-elasticity>`.
-   All attributes are passed through to the automatically-generated flex.
+   All attributes are passed through to the automatically-generated flex. In particular, :at:`yield` takes one or two
+   values, for stretching and bending; a single value is also the bending yield strain only if the generated flex is a
+   non-interpolated 2D flex.
 
 .. _flexcomp-pin:
 
@@ -4664,22 +4666,31 @@ except for 1D flexes with zero :ref:`damping<flex-elasticity-damping>`.
 
 .. _flex-elasticity-yield:
 
-:at:`yield`: :at-val:`real(1), "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
+:at:`yield`: :at-val:`real(1 or 2), "0"`
+   Dimensionless plastic yield strain thresholds :math:`\varepsilon_y` for stretching and bending of 2D and 3D flexes.
+   For non-interpolated 2D flexes, a single value applies to both, e.g., ``yield="0.01"`` is equivalent to
+   ``yield="0.01 0.01"``; for other flexes, a single value is the stretching threshold. A mode (stretching or bending)
+   is plastic if its threshold is positive at compile time; its plastic state is then allocated in ``mjData.plastic``
+   and the bit of the mode is set in ``mjModel.flex_plasticmode[f]`` (1: bending, 2: stretching). A mode whose threshold
+   is zero has no plastic state and remains purely elastic, e.g., ``yield="0.01 0"`` makes only stretching plastic. As
+   for joints and tendons, setting the threshold of a plastic mode to zero at runtime in ``mjModel.flex_yield`` (which
+   has two columns) corresponds to zero yield strain (full plastic flow). If the stretching (first) threshold is
+   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. If the bending (second) threshold is positive, which requires a non-interpolated 2D flex (a
+   compiler error otherwise), for 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). For example,
+   ``yield="0 0.001"`` makes a shell yield in bending while its stretching remains elastic. 1D flexes do not support
+   plasticity (both values must be zero). 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 plastic), followed by ``edgenum`` bending plastic entries in 2D when bending is plastic. 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.
+   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

@@ -35,7 +35,8 @@ Significant new features
   - Ideal elasto-plasticity for joint (:ref:`springyield<body-joint-springyield>`) and tendon
     (:ref:`spatial<tendon-spatial-springyield>`, :ref:`fixed<tendon-fixed-springyield>`) springs (left).
   - 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 (right).
+    with separate stretching and bending yield strains for 2D shells, volumetric-deviatoric (von Mises :math:`J_2`)
+    stretching and 2D shell bending return mapping (right).
 
   .. image:: https://www.gstatic.com/mujoco/doc/images/changelog/wire.gif
      :align: left

doc/includes/references.h

@@ -944,6 +944,7 @@ typedef struct mjModel_ {
   int*      flex_elemedgeadr;     // first element edge id address            (nflex x 1)
   int*      flex_bendingadr;      // first bending data address               (nflex x 1)
   int*      flex_plasticadr;      // plastic state address; -1: elastic       (nflex x 1)
+  int*      flex_plasticmode;     // plastic state bits; 1: bend, 2: stretch  (nflex x 1)
   int*      flex_shellnum;        // number of shells                         (nflex x 1)
   int*      flex_shelldataadr;    // first shell data address                 (nflex x 1)
   int*      flex_texcoordadr;     // address in flex_texcoord; -1: none       (nflex x 1)
@@ -974,7 +975,7 @@ typedef struct mjModel_ {
   mjtNum*   efm0_L;               // factor of M + (dt^2+dt*d)*K_bend         (nefm0L x 1)
   int*      efmC_perm;            // sparse factor order: vertex dof or -1    (nefmCvert x 1)
   mjtNum*   flex_damping;         // Rayleigh's damping coefficient           (nflex x 1)
-  mjtNum*   flex_yield;           // yield strain threshold; 0: elastic       (nflex x 1)
+  mjtNum*   flex_yield;           // yield strain (stretch, bend); 0: elastic (nflex x 2)
   mjtNum*   flex_edgestiffness;   // edge stiffness                           (nflex x 1)
   mjtNum*   flex_edgedamping;     // edge damping                             (nflex x 1)
   int*      flex_edgeequality;    // 0:none, 1:edges, 3:strain                (nflex x 1)
@@ -2147,7 +2148,7 @@ typedef struct mjsFlex_ {          // flex specification
   double shear;                    // 1D torsion shear modulus; 0: young / (2 (1 + poisson))
   double damping;                  // Rayleigh's damping
   double thickness;                // thickness (2D only)
-  double stressyield;              // yield strain threshold; 0: elastic
+  double stressyield[2];           // yield strain thresholds (stretch, bend); 0: elastic
   int elastic2d;                   // 2D passive forces; 0: none, 1: bending, 2: stretching, 3: both
   int elastic3d;  // experimental 3D material (mjSpec only); 0: Saint
                   // Venant-Kirchhoff, 1: Stable Neo-Hookean

include/mujoco/mjmodel.h

@@ -574,6 +574,7 @@ typedef struct mjModel_ {
   int*      flex_elemedgeadr;     // first element edge id address            (nflex x 1)
   int*      flex_bendingadr;      // first bending data address               (nflex x 1)
   int*      flex_plasticadr;      // plastic state address; -1: elastic       (nflex x 1)
+  int*      flex_plasticmode;     // plastic state bits; 1: bend, 2: stretch  (nflex x 1)
   int*      flex_shellnum;        // number of shells                         (nflex x 1)
   int*      flex_shelldataadr;    // first shell data address                 (nflex x 1)
   int*      flex_texcoordadr;     // address in flex_texcoord; -1: none       (nflex x 1)
@@ -604,7 +605,7 @@ typedef struct mjModel_ {
   mjtNum*   efm0_L;               // factor of M + (dt^2+dt*d)*K_bend         (nefm0L x 1)
   int*      efmC_perm;            // sparse factor order: vertex dof or -1    (nefmCvert x 1)
   mjtNum*   flex_damping;         // Rayleigh's damping coefficient           (nflex x 1)
-  mjtNum*   flex_yield;           // yield strain threshold; 0: elastic       (nflex x 1)
+  mjtNum*   flex_yield;           // yield strain (stretch, bend); 0: elastic (nflex x 2)
   mjtNum*   flex_edgestiffness;   // edge stiffness                           (nflex x 1)
   mjtNum*   flex_edgedamping;     // edge damping                             (nflex x 1)
   int*      flex_edgeequality;    // 0:none, 1:edges, 3:strain                (nflex x 1)

include/mujoco/mjspec.h

@@ -554,7 +554,7 @@ typedef struct mjsFlex_ {          // flex specification
   double shear;                    // 1D torsion shear modulus; 0: young / (2 (1 + poisson))
   double damping;                  // Rayleigh's damping
   double thickness;                // thickness (2D only)
-  double stressyield;              // yield strain threshold; 0: elastic
+  double stressyield[2];           // yield strain thresholds (stretch, bend); 0: elastic
   int elastic2d;                   // 2D passive forces; 0: none, 1: bending, 2: stretching, 3: both
   int elastic3d;  // experimental 3D material (mjSpec only); 0: Saint
                   // Venant-Kirchhoff, 1: Stable Neo-Hookean

include/mujoco/mjspecmacro.h

@@ -306,7 +306,7 @@
     X   ( double,         shear,         1      ) \
     X   ( double,         damping,       1      ) \
     X   ( double,         thickness,     1      ) \
-    X   ( double,         stressyield,   1      ) \
+    XVEC( double,         stressyield,   2      ) \
     X   ( int,            elastic2d,     1      ) \
     X   ( int,            elastic3d,     1      ) \
     X   ( mjtFlexSection, section,       1      ) \

include/mujoco/mjxmacro.h

@@ -482,6 +482,7 @@
     X   ( int,     flex_elemedgeadr,      nflex,         1                    ) \
     X   ( int,     flex_bendingadr,       nflex,         1                    ) \
     X   ( int,     flex_plasticadr,       nflex,         1                    ) \
+    X   ( int,     flex_plasticmode,      nflex,         1                    ) \
     X   ( int,     flex_shellnum,         nflex,         1                    ) \
     X   ( int,     flex_shelldataadr,     nflex,         1                    ) \
     X   ( int,     flex_texcoordadr,      nflex,         1                    ) \
@@ -512,7 +513,7 @@
     X   ( mjtNum,  efm0_L,                nefm0L,        1                    ) \
     X   ( int,     efmC_perm,             nefmCvert,     1                    ) \
     X   ( mjtNum,  flex_damping,          nflex,         1                    ) \
-    X   ( mjtNum,  flex_yield,            nflex,         1                    ) \
+    X   ( mjtNum,  flex_yield,            nflex,         2                    ) \
     X   ( mjtNum,  flex_edgestiffness,    nflex,         1                    ) \
     X   ( mjtNum,  flex_edgedamping,      nflex,         1                    ) \
     X   ( int,     flex_edgeequality,     nflex,         1                    ) \

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

@@ -741,12 +741,13 @@ def _next_plastic_flex_stretch(
   flex_stiffnessadr: wp.array[int],
   flex_elemedgeadr: wp.array[int],
   flex_plasticadr: wp.array[int],
+  flex_plasticmode: wp.array[int],
   flex_vertbodyid: wp.array[int],
   flex_elem: wp.array[int],
   flex_elemedge: wp.array[int],
   flexedge_length0: wp.array[float],
   flex_stiffness: wp.array[float],
-  flex_yield: wp.array[float],
+  flex_yield: wp.array[wp.vec2],
   flex_elemflexid: wp.array[int],
   # Data in:
   qvel_in: wp.array2d[float],
@@ -764,9 +765,9 @@ def _next_plastic_flex_stretch(
   padr_base = flex_plasticadr[f]
   if padr_base < 0:
     return
-  yield_ = wp.max(0.0, flex_yield[f])
+  yield_ = wp.max(0.0, flex_yield[f][0])
   dim = flex_dim[f]
-  if (dim != 2 and dim != 3) or flex_interp[f] != 0:
+  if (dim != 2 and dim != 3) or flex_interp[f] != 0 or (flex_plasticmode[f] & 2) == 0:
     return
   stiffness_adr_base = flex_stiffnessadr[f]
   if stiffness_adr_base < 0 or flex_stiffness[stiffness_adr_base] == 0.0:
@@ -912,13 +913,14 @@ def _next_plastic_flex_bend(
   flex_stiffnessadr: wp.array[int],
   flex_bendingadr: wp.array[int],
   flex_plasticadr: wp.array[int],
+  flex_plasticmode: wp.array[int],
   flex_vertbodyid: wp.array[int],
   flex_edge: wp.array[wp.vec2i],
   flex_edgeflap: wp.array[wp.vec2i],
   flexedge_length0: wp.array[float],
   flex_stiffness: wp.array[float],
   flex_bending: wp.array[float],
-  flex_yield: wp.array[float],
+  flex_yield: wp.array[wp.vec2],
   flex_edgeflexid: wp.array[int],
   # Data in:
   qvel_in: wp.array2d[float],
@@ -935,7 +937,7 @@ def _next_plastic_flex_bend(
   padr_base = flex_plasticadr[f]
   if padr_base < 0:
     return
-  if flex_dim[f] != 2 or flex_interp[f] != 0:
+  if flex_dim[f] != 2 or flex_interp[f] != 0 or (flex_plasticmode[f] & 1) == 0:
     return
   bending_adr_base = flex_bendingadr[f]
   stiffness_adr_base = flex_stiffnessadr[f]
@@ -1000,13 +1002,13 @@ def _next_plastic_flex_bend(
   if denom <= MJ_MINVAL or l0 <= MJ_MINVAL:
     return
 
-  yield_ = wp.max(0.0, flex_yield[f])
+  yield_ = wp.max(0.0, flex_yield[f][1])
   kappa_yield = 2.0 * yield_ / thickness
   omega = wp.dot(ed2, frc3)
   theta = wp.atan2(omega * l0, wp.dot(frc2, frc3))
   theta_yield = 0.5 * (two_a1 + two_a2) / l0 * kappa_yield
 
-  stretch_pnum = 3 * flex_elemnum[f] if flex_stiffness[stiffness_adr_base] != 0.0 else 0
+  stretch_pnum = 3 * flex_elemnum[f] if (flex_plasticmode[f] & 2) != 0 else 0
   padr = padr_base + stretch_pnum + eid
   b16_0 = flex_bending[bbase + 16]
   b_p = plastic_in[worldid, padr]
@@ -1037,7 +1039,7 @@ def _next_plastic_flex_interp(
   flex_nodebodyid: wp.array[int],
   flex_node0: wp.array[wp.vec3],
   flex_stiffness: wp.array[float],
-  flex_yield: wp.array[float],
+  flex_yield: wp.array[wp.vec2],
   flex_faceadr: wp.array[int],
   flex_cell_map: wp.array[wp.vec4i],
   flex_face: wp.array2d[int],
@@ -1152,7 +1154,7 @@ def _next_plastic_flex_interp(
       f1 += xpos_j * (phi_x * dphi_y * phi_z)
       f2 += xpos_j * (phi_x * phi_y * dphi_z)
 
-  yield_ = wp.max(0.0, flex_yield[f])
+  yield_ = wp.max(0.0, flex_yield[f][0])
   if shell_mode:
     normal_axis = 0 if ci < 2 * cy * cz else (1 if ci < 2 * (cy * cz + cx * cz) else 2)
     axis0 = (normal_axis + 1) % 3
@@ -1417,6 +1419,7 @@ def _advance(m: Model, d: Data, qacc: wp.array, qvel: Optional[wp.array] = None)
         m.flex_stiffnessadr,
         m.flex_elemedgeadr,
         m.flex_plasticadr,
+        m.flex_plasticmode,
         m.flex_vertbodyid,
         m.flex_elem,
         m.flex_elemedge,
@@ -1453,6 +1456,7 @@ def _advance(m: Model, d: Data, qacc: wp.array, qvel: Optional[wp.array] = None)
         m.flex_stiffnessadr,
         m.flex_bendingadr,
         m.flex_plasticadr,
+        m.flex_plasticmode,
         m.flex_vertbodyid,
         m.flex_edge,
         m.flex_edgeflap,

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

@@ -596,6 +596,8 @@ def _check_friction(name: str, id_: int, condim: int, friction, checks):
       m.tendon_plasticadr = np.full(mjm.ntendon, -1, dtype=np.int32)
     if not hasattr(mjm, "flex_plasticadr"):
       m.flex_plasticadr = np.full(mjm.nflex, -1, dtype=np.int32)
+    if not hasattr(mjm, "flex_plasticmode"):
+      m.flex_plasticmode = np.zeros(mjm.nflex, dtype=np.int32)
 
   m.opt = opt
   m.callback = types.Callback()

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

@@ -2652,6 +2652,7 @@ def kernel(
     flex_stiffnessadr: wp.array[int],
     flex_elemedgeadr: wp.array[int],
     flex_plasticadr: wp.array[int],
+    flex_plasticmode: wp.array[int],
     flex_edge: wp.array[wp.vec2i],
     flex_elem: wp.array[int],
     flex_elemedge: wp.array[int],
@@ -2693,7 +2694,7 @@ def kernel(
     ee_base = flex_elemedgeadr[f] + local_elemid * nedge
     stiffness_adr = flex_stiffnessadr[f] + local_elemid * (FLEX_STIFFNESS_3D if dim == 3 else 21)
     snh = wp.static(FLEX_STIFFNESS_3D == 24) and dim == 3 and flex_stiffness[stiffness_adr + 21] != 0.0
-    padr = flex_plasticadr[f] + local_elemid * nedge if flex_plasticadr[f] >= 0 else -1
+    padr = flex_plasticadr[f] + local_elemid * nedge if (flex_plasticmode[f] & 2) != 0 else -1
 
     elem_verts = wp.vec4i(-1, -1, -1, -1)
     vert_xpos = mat43()
@@ -2866,6 +2867,7 @@ def flex_hessian(m: Model, d: Data):
         m.flex_stiffnessadr,
         m.flex_elemedgeadr,
         m.flex_plasticadr,
+        m.flex_plasticmode,
         m.flex_edge,
         m.flex_elem,
         m.flex_elemedge,
@@ -3594,6 +3596,7 @@ def kernel(
     flex_stiffnessadr: wp.array[int],
     flex_elemedgeadr: wp.array[int],
     flex_plasticadr: wp.array[int],
+    flex_plasticmode: wp.array[int],
     flex_vertbodyid: wp.array[int],
     flex_elem: wp.array[int],
     flex_elemedge: wp.array[int],
@@ -3668,7 +3671,7 @@ def kernel(
 
     elemedge_adr = flex_elemedgeadr[f] + local_elemid * nedge
     edge_adr = flex_edgeadr[f]
-    padr = flex_plasticadr[f] + local_elemid * nedge if flex_plasticadr[f] >= 0 else -1
+    padr = flex_plasticadr[f] + local_elemid * nedge if (flex_plasticmode[f] & 2) != 0 else -1
 
     elongation_spring = vec6(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
     if not dsbl_spring:
@@ -4045,6 +4048,7 @@ def kernel(
     flex_stiffnessadr: wp.array[int],
     flex_bendingadr: wp.array[int],
     flex_plasticadr: wp.array[int],
+    flex_plasticmode: wp.array[int],
     flex_vertbodyid: wp.array[int],
     flex_edge: wp.array[wp.vec2i],
     flex_edgeflap: wp.array[wp.vec2i],
@@ -4110,9 +4114,8 @@ def kernel(
     bbase = bendingadr + 17 * eid
     c16 = flex_bending[bbase + 16]
     padr_base = flex_plasticadr[f]
-    if padr_base >= 0:
-      stiffness_adr = flex_stiffnessadr[f]
-      stretch_pnum = 3 * flex_elemnum[f] if (stiffness_adr >= 0 and flex_stiffness[stiffness_adr] != 0.0) else 0
+    if (flex_plasticmode[f] & 1) != 0:
+      stretch_pnum = 3 * flex_elemnum[f] if (flex_plasticmode[f] & 2) != 0 else 0
       c16 += plastic_in[worldid, padr_base + stretch_pnum + eid]
     frc = mat43()
     if has_spring and c16 != 0.0:
@@ -5318,6 +5321,7 @@ def passive(m: Model, d: Data):
         m.flex_stiffnessadr,
         m.flex_elemedgeadr,
         m.flex_plasticadr,
+        m.flex_plasticmode,
         m.flex_vertbodyid,
         m.flex_elem,
         m.flex_elemedge,
@@ -5354,6 +5358,7 @@ def passive(m: Model, d: Data):
         m.flex_stiffnessadr,
         m.flex_bendingadr,
         m.flex_plasticadr,
+        m.flex_plasticmode,
         m.flex_vertbodyid,
         m.flex_edge,
         m.flex_edgeflap,