Add the geometric term to the flex stretch stiffness - #3451
Merged
copybara-service[bot] merged 1 commit intoAug 4, 2026
Merged
Conversation
K_stretch was the Gauss-Newton Hessian of the stretch force, not its Jacobian. With elongation e_a = L_a^2 - L0_a^2 and force f = -sum_ab M_ab e_a grad(e_b)/2, K = 2 sum_ab M_ab (s_a d_a)(s_b d_b)^T + sum_a Me_a (Laplacian_a (x) I3) and only the first term was there. The second is proportional to the edge tension Me_a = sum_b M_ab e_b, so it vanishes at rest and grows with strain: the operator was first-order correct and no more. Finite-differencing it against -d(qfrc_passive)/dq on a mesh dilated by 5% gives 7.8% of the force scale; with the term it is exact to roundoff. Add only the tensile part. The geometric block is Me_a*[[I,-I],[-I,I]] over the edge's two vertices, which is positive semi-definite exactly when Me_a >= 0; a compressed edge would make K indefinite, and both consumers -- the CG constraint solver and the PCG in mjd_effSolve -- require SPD. The clamp is structural, so no eigendecomposition is needed, and it is confined to the operator: mj_flexPassiveStretch keeps the full Me_a, so no force changes. Both the matrix-free operator and the CSR assembly the effective metric builds from are updated, since they must agree. This changes how flexes with elastic2d="stretch" integrate under the implicit integrators and the effective metric -- bag.xml moves, poncho.xml is bit-identical because bending energy is quadratic and has no geometric term. The interpolated-flex path still uses its Gauss-Newton approximation, which FlexInterpDerivativesDeformed asserts.
smallquail
added a commit
to smallquail/mujoco
that referenced
this pull request
Aug 3, 2026
K_stretch was the Gauss-Newton Hessian of the stretch force, not its Jacobian. With elongation e_a = L_a^2 - L0_a^2 and force f = -sum_ab M_ab e_a grad(e_b)/2, K = 2 sum_ab M_ab (s_a d_a)(s_b d_b)^T + sum_a Me_a (Laplacian_a (x) I3) and only the first term was there. The second is proportional to the edge tension Me_a = sum_b M_ab e_b, so it vanishes at rest and grows with strain. Only its tensile part is added: the geometric block is Me_a*[[I,-I],[-I,I]], positive semi-definite exactly when Me_a >= 0, and its consumers need SPD. The clamp is structural, so no eigendecomposition is needed, and the force keeps the full Me_a. Carried here as its own commit because the FEM path below needs it: without it the humanoid-in-bag scene diverges within 0.3 s. Duplicated from PR google-deepmind#3451 against upstream, so drop this commit when that lands.
yuvaltassa
approved these changes
Aug 4, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
K_stretchwas the Gauss-Newton Hessian of the stretch force rather than its Jacobian. With elongatione_a = L_a^2 - L0_a^2and forcef = -sum_ab M_ab e_a grad(e_b)/2, the Hessian isand only the first term was present. The second is proportional to the edge tension
Me_a = sum_b M_ab e_b, so it vanishes at rest and grows with strain — the operator was first-order correct and no more.Finite-differencing against
-d(qfrc_passive)/dqon a mesh dilated by 5% gives an error of 7.8% of the force scale; with the geometric term it is exact to roundoff. The existingFlexStretchDerivativesdid not catch this because it deforms at 5e-4 strain, where the two agree — its comments documented the omission as intended behaviour.Only the tensile part is added. The geometric block is
Me_a*[[I,-I],[-I,I]]over the edge's two vertices, which is positive semi-definite exactly whenMe_a >= 0. A compressed edge would makeKindefinite, and both consumers — the CG constraint solver and the PCG inmjd_effSolve— require an SPD operator. Because the sign condition is structural, the projection is amax(Me, 0)and needs no eigendecomposition. It is confined to the operator:mj_flexPassiveStretchkeeps the fullMe_a, so no force changes.Both the matrix-free operator and the CSR assembly that the effective metric builds from are updated, since they must agree.
Impact
This changes how flexes with
elastic2d="stretch"integrate under the implicit integrators and the effective metric.model/flex/bag.xmlmoves;model/flex/poncho.xmlis bit-identical, because bending energy is quadratic in position and has no geometric term at all. Models with no flex elasticity are untouched.The interpolated-flex path still uses its own Gauss-Newton approximation, which
FlexInterpDerivativesDeformedasserts explicitly; that is left alone.Testing
DerivativeTest.FlexStretchDerivativesTensiledilates the mesh so every edge is in tension (the clamp inactive, so the operator is exact) and finite-differences the operator against the force. It fails without this change at 7.02 vs a 0.009 tolerance.