Implicit flex elasticity in the CG constraint solver via an effective metric

This CL replaces the post-hoc implicit flex correction (`flexInterp_cgsolve`) with a **linearly-implicit effective metric** `M̃ = M + (h² + h·damping)·K` carried by the CG constraint solver itself. Contact/friction forces and implicit flex elasticity are now computed against one consistent metric, instead of the solver seeing `M` and a post-solve correction changing `qacc` behind its back.

Gate (unchanged semantics): `solver="CG"` + implicit/implicitfast integrator + pyramidal cones + flex stiffness present. Newton and PGS are untouched. `solver="CG"` remains the user-facing contract — the factorization is an implementation detail of the preconditioner.

### What's in the metric

- **mjData `efm_*`** (arena, efc-like lifetime/skip semantics; built in `mj_fwdPosition`, value-refreshed in `mj_fwdVelocity`): the per-step stiffness CSR `efm_B_*`, its reverse-Cholesky factor `efm_dofid` + `efm_L_*` (nested-dissection ordered, separators-first for the reverse factorization), and the smooth-force shift `efm_c = h·K·qvel`.
- **`mjd_flexStiff_assemble`** now assembles stretch (Gauss–Newton), standard dim-2 bending, and — via the cached corotated stiffness `d->flexelem_krot` — interp stiffness (all node bodies on simple sliders: point Jacobian is I₃, `flex_centered` not required; fixed nodes drop like pins) into one dof-level CSR. `mjd_effMulAdd`/`mjd_effSolve` apply the metric, with matrix-free operator fallbacks where assembly does not apply.
- **mjModel `efm0_*`** (`nefm0dof`/`nefm0L`): the constant part of the metric factor — currently the dim-2 bending factor, computed once in `mj_setConst` — so bending-only models pay zero per-step factorization cost. Naming mirrors mjData's `efm_*` with the standard `0`-suffix (reference/constant) idiom, and is deliberately not bending-specific: future constant contributors extend it without renames.
- The solver consumes the metric through pre-shifted `qfrc_smooth` and the metric products `Ma`/`Mv`/`Mgrad`; `qacc_smooth` becomes the unconstrained minimizer of the implicit dynamics, which makes the no-constraint shortcut and the warmstart choice consistent by construction.
- **`mj_inverse` adds `B·qacc − c`**, making inverse dynamics discrete-consistent with the gated forward dynamics — exact, since the gated path has no qDeriv term (new test `ForwardTest.GatedFlexInverseConsistency`).

### Performance

All numbers: ms/step over the same 2000-step window, models as shipped on each side (old code with the old model settings vs this CL with the new ones).

The new solver path activates on exactly two shipped models — the ponchos, the only flex models that need an implicit integrator (poncho on Euler degenerates to >200 ms/step). For them, this CL trades speed for consistency: the implicit bending solve now runs inside every solver iteration, where the contact solve can see the stiffness, instead of once after the solve. Solver iterations drop because the curvature is visible, but each iteration pays for the implicit solve:

| model | before | after | solver iters/step |
|---|---|---|---|
| poncho | 2.47 | 3.30 (1.33×) | 16.8 → 11.8 |
| poncho_edgeequality | 1.96 | 2.72 (1.39×) | 13.2 → 10.0 |

What that price buys: contact forces consistent with the implicit elasticity (previously the post-hoc correction changed `qacc` after the constraint solve), discrete-consistent inverse dynamics, and the removal of the post-hoc special case from the integration path. Raising poncho's timestep from 2 to 5 ms leaves its per-step cost nearly flat, so the consistency price can be recovered by taking fewer steps where accuracy allows.

Every other flex model was measured stable on Euler at its shipped timestep and switches to it (these models predate the post-hoc integrator; implicit was never load-bearing for them). They end up equal or faster than before: bunny_multicell 0.47 → 0.40, trampoline 0.28 → 0.25, plate 1.02 → 0.99, pancake 0.34 → 0.33.

Finally, the per-step factorization makes configurations practical that the old code could only integrate explicitly: implicit stretch elasticity (`elastic2d="stretch"`/`"both"`, dim-3 solids) and factorized interp stiffness. No before/after exists for these — stock has no implicit treatment of stretch at all.

### Behavior changes

- With the post-hoc correction deleted, interp/bending models running `solver="Newton"` (or elliptic cones, or islands) now integrate flex elasticity **explicitly** (previously: post-hoc implicit). Affects e.g. `gripper_trilinear` (stable, and faster, but different semantics). Follow-up options: Newton-side metric support, or a documented fallback.
- With the gate on, `mj_forward` outputs are timestep-dependent for gated models (they answer the linearly-implicit discrete problem); `qacc_smooth` and `mj_inverse` change accordingly. Non-gated models are bit-identical (full suite green throughout).

### Validation

- 1737/1737 tests, including new: `FlexStretchDerivatives` (FD-validated GN operator), `FlexStiffAssemble`/`FlexStiffAssembleInterp` (CSR ≡ operators), `GatedFlexInverseConsistency` (fails pre-change), equivalence tests vs the old post-hoc treatment (bending matches to 2e-11).
- Fingerprint discipline throughout: bending-only models bit-exact across every refactor; permutation/kernel changes verified iteration-identical.

### Known follow-ups (not in this CL)

3×3-block sparse Cholesky kernel (the numeric factorization is index-bound; projected ~3× on the factor); mjModel persistence of the factor's symbolic pattern (rest-pose ND makes sizes compile-time); the general effective-metric mode (all solvers, all PSD-safe force classes, behind an enable flag).

PiperOrigin-RevId: 948561856
Change-Id: I8b8e32ebd0428042af71647d0470d10773bf6daf

test/engine/engine_derivative_test.cc

@@ -30,6 +30,7 @@
 #include "src/engine/engine_forward.h"
 #include "src/engine/engine_io.h"
 #include "src/engine/engine_util_blas.h"
+#include "src/engine/engine_util_sparse.h"
 #include "test/fixture.h"
 
 namespace mujoco {
@@ -1771,5 +1772,376 @@ TEST_F(DerivativeTest, FlexInterpDerivativesDeformed) {
       << "Jacobian approximation should differ from FD when deformed";
 }
 
+// Helper: assemble the standard-flex stretch stiffness into a dense matrix,
+// column-by-column using mjd_flexStretch_mul with scale (s1 + s2*damping).
+static void stretchK_dense(mjModel* m, mjData* d, mjtNum* K, int nv,
+                           mjtNum s1, mjtNum s2) {
+  std::vector<mjtNum> e_i(nv, 0);
+  std::vector<mjtNum> col(nv, 0);
+  for (int i = 0; i < nv; i++) {
+    mju_zero(e_i.data(), nv);
+    mju_zero(col.data(), nv);
+    e_i[i] = 1.0;
+    mjd_flexStretch_mul(m, d, col.data(), e_i.data(), s1, s2);
+    for (int j = 0; j < nv; j++) {
+      K[j * nv + i] = col[j];
+    }
+  }
+}
+
+// verify mjd_flexStretch_mul (Gauss-Newton Hessian of the standard-flex
+// stretch force) against finite differences of qfrc_passive, plus symmetry,
+// positive semi-definiteness and (s1, s2) scale linearity. The model covers
+// both element edge tables (dim=2 triangles and dim=3 tets) and a pinned
+// vertex (zero-dof body guard).
+TEST_F(DerivativeTest, FlexStretchDerivatives) {
+  static const char* const kXml = R"(
+  <mujoco>
+    <option integrator="implicit"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="4 4 1" spacing="0.1 0.1 0.1"
+                radius=".01" dim="2" mass="1" pos="0 0 1">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" thickness="0.01"
+                    elastic2d="stretch" damping="50"/>
+        <pin id="0"/>
+      </flexcomp>
+      <flexcomp name="solid" type="grid" count="3 3 3" spacing="0.1 0.1 0.1"
+                radius=".01" dim="3" mass="1" pos="1 0 1">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" damping="10"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(kXml, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  int nv = model->nv;
+  ASSERT_EQ(model->nq, nv);  // all slide dofs
+
+  MjDataPtr data = MakeData(model);
+
+  // deform both flexes deterministically. Keep the strain small: the operator
+  // is the Gauss-Newton Hessian, exact to O(strain) (the geometric term is
+  // dropped, see FlexInterpDerivativesDeformed for the analogous property).
+  for (int i = 0; i < nv; i++) {
+    data->qpos[i] += 5e-4 * (mju_Halton(i, 2) - 0.5);
+  }
+  mj_forward(model.get(), data.get());
+
+  // part 1: FD verification of K*vec against qfrc_passive (qvel = 0, so the
+  // kD elongation term vanishes and qfrc_passive is the pure stretch spring)
+  {
+    std::vector<mjtNum> vec(nv), res(nv, 0);
+    for (int i = 0; i < nv; i++) {
+      vec[i] = mju_Halton(i, 2) - 0.5;
+    }
+    mjd_flexStretch_mul(model.get(), data.get(), res.data(), vec.data(), 1, 0);
+
+    mjtNum eps = MjTol(1e-7, 1e-4);
+    mjData* data_perturbed = mj_copyData(NULL, model.get(), data.get());
+    mju_addToScl(data_perturbed->qpos, vec.data(), eps, nv);
+    mj_forward(model.get(), data_perturbed);
+
+    // qfrc_passive = -dV/dq  =>  -(qfrc_new - qfrc)/eps ~= K * vec.
+    // Compare max error against the force scale: the operator omits the
+    // geometric (stress-proportional) term, so the residual is O(strain) of
+    // the overall scale and individual near-zero entries are not meaningful.
+    mjtNum max_err = 0, scale = 0;
+    for (int i = 0; i < nv; ++i) {
+      mjtNum fd =
+          -(data_perturbed->qfrc_passive[i] - data->qfrc_passive[i]) / eps;
+      max_err = mju_max(max_err, mju_abs(res[i] - fd));
+      scale = mju_max(scale, mju_abs(fd));
+    }
+    EXPECT_GT(scale, 1.0) << "test should exercise nontrivial stiffness";
+    EXPECT_LT(max_err, MjTol(5e-3, 5e-2) * scale)
+        << "stretch stiffness mismatch: max_err " << max_err
+        << " at force scale " << scale;
+    mj_deleteData(data_perturbed);
+  }
+
+  // part 2: symmetry and positive semi-definiteness of the assembled K
+  {
+    std::vector<mjtNum> K(nv * nv, 0);
+    stretchK_dense(model.get(), data.get(), K.data(), nv, 1, 0);
+
+    mjtNum max_asymmetry = 0;
+    for (int i = 0; i < nv; i++) {
+      for (int j = 0; j < i; j++) {
+        max_asymmetry =
+            mju_max(max_asymmetry, mju_abs(K[i * nv + j] - K[j * nv + i]));
+      }
+    }
+    EXPECT_THAT(max_asymmetry, MjNear(0, 1e-10, 5e-4))
+        << "K_stretch is not symmetric";
+
+    for (int trial = 0; trial < 5; trial++) {
+      std::vector<mjtNum> v(nv);
+      for (int i = 0; i < nv; i++) {
+        v[i] = mju_Halton(i + trial * nv, 3) - 0.5;
+      }
+      mjtNum vKv = 0;
+      for (int i = 0; i < nv; i++) {
+        for (int j = 0; j < nv; j++) {
+          vKv += v[i] * K[i * nv + j] * v[j];
+        }
+      }
+      EXPECT_GE(vKv, MjTol(-1e-8, -1e-5)) << "K_stretch is not PSD";
+    }
+  }
+
+  // part 3: (s1, s2) scale linearity across flexes with different damping:
+  // mul(s1, s2) == s1*mul(1, 0) + s2*mul(0, 1)
+  {
+    std::vector<mjtNum> vec(nv), a(nv, 0), b(nv, 0), c(nv, 0);
+    for (int i = 0; i < nv; i++) {
+      vec[i] = mju_Halton(i, 5) - 0.5;
+    }
+    mjtNum h = 1e-3;
+    mjd_flexStretch_mul(model.get(), data.get(), a.data(), vec.data(),
+                        h * h, h);
+    mjd_flexStretch_mul(model.get(), data.get(), b.data(), vec.data(), 1, 0);
+    mjd_flexStretch_mul(model.get(), data.get(), c.data(), vec.data(), 0, 1);
+    for (int i = 0; i < nv; i++) {
+      EXPECT_THAT(a[i], MjNear(h * h * b[i] + h * c[i], 1e-12, 1e-5))
+          << "scale linearity mismatch at DOF " << i;
+    }
+  }
+}
+
+// verify mjd_flexStiff_assemble against the matrix-free operators: the
+// assembled CSR applied to test vectors must reproduce mjd_flexBend_mul +
+// mjd_flexStretch_mul at the same state
+TEST_F(DerivativeTest, FlexStiffAssemble) {
+  static const char* const kXml = R"(
+  <mujoco>
+    <option integrator="implicit"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="4 4 1" spacing="0.1 0.1 0.1"
+                radius=".01" dim="2" mass="1" pos="0 0 1">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" thickness="0.01"
+                    elastic2d="both" damping="7"/>
+      </flexcomp>
+      <flexcomp name="solid" type="grid" count="3 3 3" spacing="0.1 0.1 0.1"
+                radius=".01" dim="3" mass="1" pos="1 0 1">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" damping="10"/>
+        <pin id="0"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(kXml, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  int nv = model->nv;
+  MjDataPtr data = MakeData(model);
+
+  // deform deterministically
+  for (int i = 0; i < nv; i++) {
+    data->qpos[i] += 2e-3 * (mju_Halton(i, 2) - 0.5);
+  }
+  mj_forward(model.get(), data.get());
+
+  // assemble both terms with a mixed (s1, s2) scale
+  mjtNum s1 = 4e-6, s2 = 2e-3;
+  std::vector<int> rownnz(nv), rowadr(nv);
+  int nnz = mjd_flexStiff_assemble(model.get(), data.get(), rownnz.data(),
+                                   rowadr.data(), NULL, NULL, s1, s2,
+                                   /*flg_bend=*/1, /*flg_stretch=*/1, NULL);
+  ASSERT_GT(nnz, 0);
+  std::vector<int> colind(nnz);
+  std::vector<mjtNum> val(nnz);
+  mjd_flexStiff_assemble(model.get(), data.get(), rownnz.data(), rowadr.data(),
+                         colind.data(), val.data(), s1, s2, /*flg_bend=*/1,
+                         /*flg_stretch=*/1, NULL);
+
+  // compare CSR apply vs operators on test vectors
+  for (int trial = 0; trial < 3; trial++) {
+    std::vector<mjtNum> vec(nv), res_op(nv, 0), res_csr(nv, 0);
+    for (int i = 0; i < nv; i++) {
+      vec[i] = mju_Halton(i + trial*nv, 3) - 0.5;
+    }
+    mjd_flexBend_mul(model.get(), data.get(), res_op.data(), vec.data(), s1,
+                     s2);
+    mjd_flexStretch_mul(model.get(), data.get(), res_op.data(), vec.data(), s1,
+                        s2);
+    for (int i = 0; i < nv; i++) {
+      mjtNum sum = 0;
+      for (int k = 0; k < rownnz[i]; k++) {
+        sum += val[rowadr[i] + k]*vec[colind[rowadr[i] + k]];
+      }
+      res_csr[i] = sum;
+    }
+    for (int i = 0; i < nv; i++) {
+      EXPECT_THAT(res_csr[i], MjNear(res_op[i], 1e-12, 2e-5))
+          << "assembly/operator mismatch at DOF " << i << " trial " << trial;
+    }
+  }
+}
+
+// verify the interp assembly mode: with the K_rot cache supplied, the assembled CSR applied
+// to test vectors must reproduce mjd_flexInterp_mul, whose sign convention is negated
+TEST_F(DerivativeTest, FlexStiffAssembleInterp) {
+  static const char* const kXml = R"(
+  <mujoco>
+    <option integrator="implicit"/>
+    <worldbody>
+      <flexcomp name="soft" type="grid" count="4 4 4" spacing="0.1 0.1 0.1"
+                radius=".01" dim="3" mass="1" pos="0 0 1" dof="trilinear">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" damping="2"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(kXml, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  int nv = model->nv;
+  MjDataPtr data = MakeData(model);
+  ASSERT_EQ(mjd_flexInterpAssemblable(model.get()), 1);
+
+  // deform deterministically, refresh kinematics, cache the corotated stiffness
+  for (int i = 0; i < nv; i++) {
+    data->qpos[i] += 2e-3 * (mju_Halton(i, 2) - 0.5);
+  }
+  mj_forward(model.get(), data.get());
+  std::vector<mjtNum> krot(model->nflexstiffness, 0);
+  mjd_flexInterp_cacheKrot(model.get(), data.get(), krot.data());
+
+  // assemble interp only
+  mjtNum s1 = 4e-6, s2 = 2e-3;
+  std::vector<int> rownnz(nv), rowadr(nv);
+  int nnz = mjd_flexStiff_assemble(model.get(), data.get(), rownnz.data(), rowadr.data(),
+                                   NULL, NULL, s1, s2, /*flg_bend=*/0, /*flg_stretch=*/0,
+                                   krot.data());
+  ASSERT_GT(nnz, 0);
+  std::vector<int> colind(nnz);
+  std::vector<mjtNum> val(nnz);
+  mjd_flexStiff_assemble(model.get(), data.get(), rownnz.data(), rowadr.data(),
+                         colind.data(), val.data(), s1, s2, /*flg_bend=*/0, /*flg_stretch=*/0,
+                         krot.data());
+
+  // compare CSR apply vs the operator called with negated scales (its convention)
+  for (int trial = 0; trial < 3; trial++) {
+    std::vector<mjtNum> vec(nv), res_op(nv, 0), res_csr(nv, 0);
+    for (int i = 0; i < nv; i++) {
+      vec[i] = mju_Halton(i + trial*nv, 3) - 0.5;
+    }
+    mjd_flexInterp_mul(model.get(), data.get(), res_op.data(), vec.data(), -s1, -s2,
+                       krot.data());
+    for (int i = 0; i < nv; i++) {
+      mjtNum sum = 0;
+      for (int k = 0; k < rownnz[i]; k++) {
+        sum += val[rowadr[i] + k]*vec[colind[rowadr[i] + k]];
+      }
+      res_csr[i] = sum;
+    }
+    for (int i = 0; i < nv; i++) {
+      EXPECT_THAT(res_csr[i], MjNear(res_op[i], 1e-12, 2e-5))
+          << "interp assembly/operator mismatch at DOF " << i << " trial " << trial;
+    }
+  }
+}
+
+// mjd_effSolve: exact-preconditioner fast path solves (M+K)x = b directly; the general
+// refinement path stays within its tolerance when exactness does not hold
+TEST_F(DerivativeTest, EffSolveExact) {
+  // relative residual of (M+K)x - b after mjd_effSolve
+  auto solve_residual = [](const mjModel* m, mjData* d) {
+    int nv = m->nv;
+    std::vector<mjtNum> b(nv), x(nv), r(nv);
+    for (int i = 0; i < nv; i++) {
+      b[i] = mju_Halton(i, 3) - 0.5;
+    }
+    mjd_effSolve(m, d, x.data(), b.data());
+    mju_mulSymVecSparse(r.data(), d->M, x.data(), nv, m->M_rownnz, m->M_rowadr, m->M_colind);
+    mjd_effMulAdd(m, d, r.data(), x.data());
+    mju_subFrom(r.data(), b.data(), nv);
+    return mju_norm(r.data(), nv) / mju_norm(b.data(), nv);
+  };
+
+  // stretch + bending cloth on world: per-step factor, exact
+  static const char* const kXmlBoth = R"(
+  <mujoco>
+    <option solver="CG" integrator="implicitfast"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="6 6 1" spacing="0.05 0.05 0.05"
+                radius=".005" dim="2" mass="0.5" pos="0 0 1" dof="full">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e3" poisson="0.2" damping="0.1" elastic2d="both" thickness="0.01"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(kXmlBoth, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  MjDataPtr data = MakeData(model);
+  mj_forward(model.get(), data.get());
+  ASSERT_GE(data->efm_active, 1);
+  EXPECT_GT(data->nefmK, 0);
+  EXPECT_GT(data->nefmdof, 0);
+  EXPECT_EQ(data->efm_active, 2);
+  EXPECT_LT(solve_residual(model.get(), data.get()), MjTol(1e-10, 1e-6));
+
+  // bending-only cloth: no CSR or per-step factor, constant factor covers, exact
+  static const char* const kXmlBend = R"(
+  <mujoco>
+    <option solver="CG" integrator="implicitfast"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="6 6 1" spacing="0.05 0.05 0.05"
+                radius=".005" dim="2" mass="0.5" pos="0 0 1" dof="full">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e3" poisson="0.2" damping="0.1" elastic2d="bend" thickness="0.01"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+  model = LoadModelFromString(kXmlBend, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  data = MakeData(model);
+  mj_forward(model.get(), data.get());
+  ASSERT_GE(data->efm_active, 1);
+  EXPECT_EQ(data->nefmK, 0);
+  EXPECT_EQ(data->nefmdof, 0);
+  EXPECT_GT(model->nefm0dof, 0);
+  EXPECT_EQ(data->efm_active, 2);
+  EXPECT_LT(solve_residual(model.get(), data.get()), MjTol(1e-10, 1e-6));
+
+  // cloth under a jointed parent: M couples across the covered block, not exact,
+  // the refinement path must still meet its tolerance
+  static const char* const kXmlMoving = R"(
+  <mujoco>
+    <option solver="CG" integrator="implicitfast"/>
+    <worldbody>
+      <body name="base" pos="0 0 1">
+        <joint type="slide" axis="0 0 1"/>
+        <geom type="sphere" size=".01" mass="1" contype="0" conaffinity="0"/>
+        <flexcomp name="cloth" type="grid" count="6 6 1" spacing="0.05 0.05 0.05"
+                  radius=".005" dim="2" mass="0.5" pos="0 0 0" dof="full">
+          <contact selfcollide="none" contype="0" conaffinity="0"/>
+          <elasticity young="1e3" poisson="0.2" damping="0.1" elastic2d="both" thickness="0.01"/>
+        </flexcomp>
+      </body>
+    </worldbody>
+  </mujoco>
+  )";
+  model = LoadModelFromString(kXmlMoving, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  data = MakeData(model);
+  mj_forward(model.get(), data.get());
+  ASSERT_GE(data->efm_active, 1);
+  EXPECT_EQ(data->efm_active, 1);
+  EXPECT_LT(solve_residual(model.get(), data.get()), 1e-4);
+}
+
 }  // namespace
 }  // namespace mujoco

test/engine/engine_forward_test.cc

@@ -2933,7 +2933,7 @@ TEST_F(ForwardTest, FlexTrilinearInstability) {
 TEST_F(ForwardTest, FlexDampingRigidMotion) {
   constexpr char xml[] = R"(
   <mujoco>
-      <option gravity="0 0 0" timestep="0.01" integrator="implicitfast"/>
+      <option gravity="0 0 0" timestep="0.01" integrator="implicitfast" solver="CG"/>
       <worldbody>
           <flexcomp name="flex" type="grid" count="3 3 3" spacing="0.1 0.1 0.1"
                     pos="0 0 0" euler="45 45 45" radius="0.01" dim="3" mass="1" dof="trilinear">
@@ -2984,7 +2984,7 @@ TEST_F(ForwardTest, FlexDampingRigidMotion) {
 TEST_F(ForwardTest, FlexParentCoupling) {
   static const char* const kXml = R"(
   <mujoco>
-    <option integrator="implicit" timestep="0.01"/>
+    <option integrator="implicit" timestep="0.01" solver="CG"/>
     <worldbody>
       <body name="parent" pos="0 0 0">
         <freejoint/>
@@ -3034,14 +3034,17 @@ TEST_F(ForwardTest, FlexParentCoupling) {
     if (diff > max_diff) max_diff = diff;
   }
 
-  EXPECT_LT(max_diff, MjTol(2e-5, 1.5e-2))
+  // tolerance rebaselined 2e-5 -> 5e-5 for the in-solver implicit flex treatment: implicit
+  // and explicit flex damping legitimately differ at O(h*damping*K/M) in this comparison, and
+  // the in-solver form lands at ~3e-5 where the old post-hoc operator landed just under 2e-5
+  EXPECT_LT(max_diff, MjTol(5e-5, 1.5e-2))
       << "Implicit integrator should match Euler at small timestep";
 }
 
 TEST_F(ForwardTest, TrilinearPinnedParentWithFreejoint) {
   static constexpr char xml[] = R"(
   <mujoco>
-  <option integrator="implicitfast"/>
+  <option integrator="implicitfast" solver="CG"/>
   <worldbody>
     <body>
       <joint type="free"/>
@@ -3601,7 +3604,7 @@ TEST_F(ImplicitIntegratorTest, FlexContactEnergy) {
 TEST_F(ImplicitIntegratorTest, BendingDampingDecaysEnergy) {
   static constexpr char xml[] = R"(
   <mujoco>
-    <option gravity="0 0 0" timestep="0.001" integrator="implicitfast">
+    <option gravity="0 0 0" timestep="0.001" integrator="implicitfast" solver="CG">
       <flag energy="enable"/>
     </option>
     <worldbody>
@@ -3661,7 +3664,7 @@ TEST_F(ImplicitIntegratorTest, BendingDampingDecaysEnergy) {
 TEST_F(ImplicitIntegratorTest, InterpStretchEnergy) {
   static constexpr char xml[] = R"(
   <mujoco>
-    <option gravity="0 0 0" timestep="0.001" integrator="implicitfast">
+    <option gravity="0 0 0" timestep="0.001" integrator="implicitfast" solver="CG">
       <flag energy="enable"/>
     </option>
     <worldbody>
@@ -3707,4 +3710,46 @@ TEST_F(ImplicitIntegratorTest, InterpStretchEnergy) {
 }
 
 }  // namespace
+// with the implicit effective metric active, inverse dynamics must recover the applied force
+// (zero here): the forward solve is (M+B)*qacc = qfrc_smooth + c + J'*f and the inverse adds
+// the same B*qacc - c terms. This is the fwd/inv consistency fence for the flex-CG dispatch.
+TEST_F(ForwardTest, GatedFlexInverseConsistency) {
+  static const char* const kXml = R"(
+  <mujoco>
+    <option solver="CG" integrator="implicitfast" tolerance="1e-14"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="6 6 1" spacing="0.1 0.1 0.1"
+                radius=".01" dim="2" mass="1" pos="0 0 1">
+        <contact selfcollide="none" contype="0" conaffinity="0"/>
+        <elasticity young="1e4" poisson="0.3" thickness="0.01"
+                    elastic2d="both" damping="0.5"/>
+      </flexcomp>
+    </worldbody>
+  </mujoco>
+  )";
+
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(kXml, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  MjDataPtr data = MakeData(model);
+  int nv = model->nv;
+
+  // deform and settle a few steps under gravity
+  for (int i=0; i < nv; i++) {
+    data->qvel[i] = 0.1 * (mju_Halton(i, 3) - 0.5);
+  }
+  for (int step=0; step < 50; step++) {
+    mj_step(model.get(), data.get());
+  }
+
+  // forward then inverse at the same state
+  mj_forward(model.get(), data.get());
+  mj_inverse(model.get(), data.get());
+
+  // no applied forces: the inverse must return ~zero, at the scale of the passive forces
+  mjtNum scale = 1 + mju_norm(data->qfrc_passive, nv);
+  EXPECT_LT(mju_norm(data->qfrc_inverse, nv), 1e-6 * scale);
+}
+
+
 }  // namespace mujoco

test/engine/engine_island_test.cc

@@ -373,6 +373,39 @@ TEST_F(IslandTest, IslandFlex) {
   mj_deleteModel(model);
 }
 
+// stiffness couples all vertices of a flex: one contact anywhere on the flex
+// must pull every vertex tree (and the contacting body) into a single island
+TEST_F(IslandTest, FlexStiffnessUnionsTrees) {
+  static const char xml[] = R"(
+  <mujoco>
+    <option solver="Newton"/>
+    <worldbody>
+      <flexcomp name="cloth" type="grid" count="4 4 1" spacing="0.1 0.1 0.1"
+                radius=".005" dim="2" mass="0.5" pos="0 0 1" dof="full">
+        <contact selfcollide="none"/>
+        <elasticity young="1e3" poisson="0.2" damping="0.1" elastic2d="both" thickness="0.01"/>
+      </flexcomp>
+      <body pos="0.1 0.1 0.96">
+        <freejoint/>
+        <geom type="sphere" size="0.05"/>
+      </body>
+    </worldbody>
+  </mujoco>
+  )";
+  char error[1024];
+  MjModelPtr model = LoadModelFromString(xml, error, sizeof(error));
+  ASSERT_THAT(model.get(), NotNull()) << error;
+  MjDataPtr data = MakeData(model);
+  mj_forward(model.get(), data.get());
+
+  // the sphere penetrates the cloth at one corner
+  ASSERT_GT(data->ncon, 0);
+
+  // one island containing every dof: 16 vertices and the free sphere
+  EXPECT_EQ(data->nisland, 1);
+  EXPECT_EQ(data->nidof, model->nv);
+}
+
 TEST_F(IslandTest, IslandEfcElliptic) {
   const std::string xml_path = GetTestDataFilePath(kIlslandEfcPath);
   char error[1024];