Contact of a flex with `passive` collisions enabled was applied as an
explicit spring of fixed stiffness 1e4, which the timestep bounds: any
stiffness worth having oscillates faster than the step can resolve, so
the force was too soft to keep sheets apart and interpenetration was
routine.
Carry its curvature in the effective metric M + K instead, alongside the
flex's own stretch and bending stiffness. The contact block k*J^T*J is
appended to the per-vertex candidate list already assembled for the flex
stencils, so it costs additional entries in an existing matrix rather
than a new one, and the accompanying shift -h*K*v is what damps the
stiff modes. At a 2 ms timestep this holds roughly 50x the stiffness an
explicit force of the same step could.
With the timestep no longer setting the bound, the stiffness is chosen
as a natural frequency scaled by the participating vertex mass rather
than left at a fixed 1e4, so one value suits models of any scale.
Passive handling is scoped to contacts whose every dof is a flex vertex
carried by the metric: flex against flex, flex against itself, and flex
against static geometry, which contributes no dofs of its own. For those
the Hessian is assembled in full. Contact with a body that can move
would have that body's dofs dropped from it, and is left on the
constraint solver.
The feature now requires an integrator whose constraint solve runs in
that metric, and is rejected with an error otherwise.
Add model/flex/drape.xml as the example model, replacing sphere_passive,
whose contacts no longer demonstrated the feature.
The new spotlight attribute softness (real in [0, 1], default 0) is the
fraction of the cone, measured inward from the cutoff, over which
intensity falls to zero. It is used by physically-based lighting models;
the Phong model's corresponding knob remains exponent.
The filament renderer previously hardcoded the inner cone angle to 0,
making the entire beam penumbra: the shader attenuates by the squared
smoothstep ((cos(theta) - cos(outer)) / (cos(inner) - cos(outer)))^2, so
a cutoff-25 spot delivered its rated candela only exactly on-axis and
about a third of it averaged over the light pool, with the deficit
shrinking as the cutoff widens. The inner angle is now
(1 - softness) * cutoff, so at the default the light delivers its full
intensity everywhere inside the cone and illuminance follows E = I/d^2
independent of the cutoff. Setting softness to 1 reproduces the previous
appearance exactly (verified bit-identical), which is the migration path
for models tuned against the old behavior.
The filament light type also changes from FOCUSED_SPOT to SPOT. With
intensity given in candela and the cone set at build time the two types
produce identical output (FOCUSED_SPOT's power-conserving rescale only
applies when the cone changes after the intensity is set), but SPOT
guarantees that candela never rescales with cone angle should the cone
ever become runtime-editable.
Verified with headless renders under a linear tone mapper against an
equal-candela point light at cutoffs 25/45/80: softness 0 gives
spot/point luminance ratio 1.000 at all sampled angles inside the cone;
softness 0.2 is flat over the inner 80% of the cone; softness 1 matches
the previous renderer with zero linear-pixel difference. XML round-trip
and the [0, 1] compile-time check verified. Introspect and wasm bindings
regenerated.
PiperOrigin-RevId: 959334706
Change-Id: I0f0729781899880de1729ea9b3d8c055d715a025
Every step, the flex block of the implicit effective metric M + K was
factorized by sparse Cholesky, because K depends on the configuration. On
model/flex/bag.xml, added here, that is roughly half the step, against a
comparable share for the constraint solve it exists to accelerate.
Keep only the metric's per-vertex 3x3 diagonal blocks, prefactored. Neither
consumer needs the exact inverse: the CG constraint solver only wants a
preconditioner, and qacc_smooth can come from an iterative solve using those
blocks. They are O(n) to build and to apply, but weaker, so CG runs about twice
the iterations and qacc_smooth becomes an iteration rather than a direct solve.
Net, the bag model steps roughly twice as fast.
The preconditioner, by metric state. Inactive, meaning no flex elasticity or an
explicit integrator: M^-1, unchanged. Bending only (nefmK == 0): M^-1 plus the
exact constant bending factor from mj_setConst on the dofs it covers,
unchanged; that factor is built at model compile time and costs nothing per
step. Per-step stiffness: M^-1 plus the 3x3 blocks, where before it was a
per-step sparse Cholesky, or, when M couples across the flex block, an inner
PCG of up to 50 iterations run once per outer CG iteration.
Only models carrying per-step stretch stiffness change in wall-clock. Both
ponchos hold their timing and take slightly fewer CG iterations than before,
because the preconditioner is now symmetric: it applies M^-1 and the covered
blocks to disjoint sets of dofs, where previously the two overlapped and the
operator was not symmetric, which PCG requires.
mjd_effSolve is the accurate solve of (M + K)x = b; what used to carry that
name only preconditions and is now mjd_effPrec. Its CG guarded the division by
pAp with mjMINVAL, an absolute floor on a quantity that scales with the square
of the right-hand side, so a small b aborted the solve while the curvature was
healthy: four flex models were quietly left short of tolerance. For an SPD
metric the guard is positivity, and with that the same solves converge. The qacc_smooth call site in
mj_fwdAcceleration is textually unchanged but now reaches the iterative solve,
which converges on opt.tolerance rather than a hardcoded threshold, floored in
mjUSESINGLE builds where the squared target is unreachable in float. Reaching
the iteration cap names the ill-conditioned flex stiffness and then reports it
through mjWARN_INERTIA, rather than returning an under-converged result.
Covered dofs are located by walking the covered rows of the stiffness matrix,
as they need not be 3-aligned from dof 0: any joint declared before a flexcomp
shifts them.
mjData.efm_L_rownnz, efm_L_rowadr and efm_L_colind described the sparsity of
the deleted factorization and are removed: left NULL with nonzero mjxmacro
extents they made the Python bindings hand back uninitialized arrays.
efm_active loses the value 2 for the same reason, nothing selects a solve path
on preconditioner exactness any more. Both are recorded under breaking changes.
model/flex/bag.xml is added because no shipped model carried per-step stretch
stiffness. The ponchos are bending-only and trampoline.xml uses an explicit
integrator, so the metric never activates there. It is excluded from
WriteReadCompareTest: stretch stiffness amplifies rest geometry that XML rounds
on save.
https://youtu.be/GioWwB36XHI
The new geom attribute adhesion (units of force, signed; pair-level
override) translates the contact friction cone along its normal so
that the force origin lies strictly inside it. Consequences: each
contact can pull with up to the given force before breaking, and the
tangential friction budget becomes mu*(f_N + adhesion) -- the
Mohr-Coulomb yield condition with cohesion c = mu*adhesion -- so
lightly-squeezed grasps retain a guaranteed friction floor.
A translated cone factors exactly into {constant attractive force}
+ {original cone}, so no solver kernels change. The implementation is
this factorization: a constant attraction along contact normals
accumulated into the new mjData.qfrc_adhesion (summed into
qfrc_passive), plus a bias of adhesive contact rows' reference
acceleration (aref += R*adhesion), which makes resting penetration
exactly independent of adhesion. Contacts of adhesive pairs remain
active throughout the gap zone, producing rows with positive violation
whose reference acceleration pulls: a tether that resists pull-off
smoothly, captures objects released within the band into steady
contact, and detaches at the specified force. Adhesion values of the
two geoms combine by sum; explicit pairs override.
mj_contactForce reports the net interface force (cone force minus the
adhesive pull), whose normal component can now be negative. Negative
adhesion is allowed and produces a repulsive offset (air hockey).
PiperOrigin-RevId: 950858148
Change-Id: I879c08eba7ae501e5c0f8c2f807167344da4c2bc
https://www.youtube.com/watch?v=PdSdrqhSiZA
The new geom attribute surfacevel (6 numbers: linear and angular velocity in the geom's local frame, angular about the geom frame origin) specifies the velocity of the geom's surface material relative to the geom frame. The relative surface velocity of the two geoms is added to the tangential contact rows of efc_vel in mj_referenceConstraint, so friction drives touching bodies toward the motion of the surface: objects on a conveyor are transported at belt speed, turntables impart omega x r with torsional spin-up for condim >= 4, and surface velocities compose with each other and with body motion. The component along the contact normal is projected out: probe experiments showed that velocity-space emission chatters mass-independently and ingestion merely deepens penetration; normal-direction effects belong to force-space features.
surfacevel is interpreted in the geom frame as authored: for mesh geoms, whose compiled frame absorbs the mesh centering and principal-axes transform, the compiler re-expresses the authored value in the compiled frame.
No special interaction with sleeping: objects being transported do not fall asleep because they are moving; objects at rest on an active surface may sleep like any other resting object.
Includes showcase models (model/surfacevel/): a luggage carousel whose ring is a spinning square-profile supertorus fed by a cascade of belts with matched spinning end rollers, bags dropping in and circulating indefinitely; and a treadmill with a passive humanoid.
PiperOrigin-RevId: 948647785
Change-Id: I0c6559a91cc7ece1237eb8ac2e51986e7342d962
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
Previously shell mode required cellcount=1 along at least one axis. This CL
adds support for cellcount > 1 in all three axes by pinning interior grid nodes
to the parent body and reconstructing their positions from boundary nodes via
Transfinite Interpolation (TFI).
PiperOrigin-RevId: 924314800
Change-Id: I8c2438f4866dd4133feed65f535a1ab69f0c9188
Warp is not a true stretch mode but a parasitic shear mode so it should be scaled like bending (~thickness^3) rather than stretching (~thickness). The thickness parameter can therefore be used to set the stretch to warp stiffness ratio of the strain constraints.
PiperOrigin-RevId: 916900311
Change-Id: I52e54210fbb2a6d309abc1f084934830adce3423
Standard flex (flex_interp=0) with thin-plate bending treated bending forces purely explicitly. This caused contact-induced vertex vibrations and non-physical energy injection for flat resting sheets, because the solver treated each vertex as an independent mass during contact and contact normals are orthogonal to stretch constraints.
Fix: extend the existing preconditioned CG solver to include the constant bending stiffness K_bend in the implicit operator via matrix-free mat-vec.
PiperOrigin-RevId: 914774020
Change-Id: I45e0d6749abb6f873566203bccae956514b2576b
with a preconditioned Conjugate Gradient (CG) solver that operates
on the full system matrix.
The previous approach extracted flex DOFs into a reduced banded system,
factored it separately, and overwrote the global solve. This required
precomputed bandwidth (makeFlexBandwidth), parent-joint detection,
coupling corrections, and a FlexInterpContext struct — and only worked
for standalone flex trees without parent joints.
The new CG solver uses the already-factored global system (M - h*qDeriv)
as a preconditioner and adds the flex stiffness contribution via
matrix-free products (mjd_flexInterp_mulKD/mulK). This handles any
kinematic configuration — including flexes attached to articulated
chains or with parent joints — without sparsity pattern restrictions.
Before (`bunny_multicell`):
```
Simulation time : 50.80 s
Steps per second : 197
Realtime factor : 0.20 x
Time per step : 5080.3 µs
CG iters / step : 3.16
Contacts / step : 31.04
Constraints / step : 124.15
Degrees of freedom : 178
Dynamic memory usage : 0.4% of 100M
```
After:
```
Simulation time : 9.52 s
Steps per second : 1051
Realtime factor : 1.05 x
Time per step : 951.7 µs
CG iters / step : 3.21
Contacts / step : 30.90
Constraints / step : 123.61
Degrees of freedom : 178
Dynamic memory usage : 0.3% of 100M
```
PiperOrigin-RevId: 913758038
Change-Id: If5aa617b2d535c86aec9bd71c9e0003a2b38bdd7
This change adds a new passive force computation for flexes with elastic2d="bend" and dof="trilinear". The bending energy is based on the squared difference of normals between adjacent face elements at their shared edge midpoint. The edge data is precomputed during model compilation and stored in flex_bending.
PiperOrigin-RevId: 910772638
Change-Id: I3b12c7b7f1ba6ac1875df495d89e8cfec921ca80
Modify EigendecomposeStiffness to project out rigid body translations and rotations from the stiffness matrix eigenvectors. This prevents "ghost damping" from rigid body modes being incorrectly constrained. Update hollow_vs_solid.xml with solimp parameters for edge constraints. Adjust engine_core_constraint_test.cc to reflect the reduced number of equality constraints due to the projection.
PiperOrigin-RevId: 910631207
Change-Id: Ibaaa59bc030cfc6ae657d8f0d1b2a51002cc8d4e
When elastic2d="stretch" is set on an interpolated flexcomp, treat the bounding box boundary as membrane elements rather than volumetric cells. This computes plane-stress stiffness over the boundary faces and updates the runtime force/derivative kernels accordingly.
Interior vertex tracking (moving vertices that follow the deforming shell) is not yet implemented so all mesh vertices need to be on the bounding box surface or the background grid should have no interior nodes (i.e. cellcount should be 1 on at least one axis).
PiperOrigin-RevId: 907654080
Change-Id: I51b90e2f6a1d1b036f9604e42de20e377dc5d3f9
Do not allow a mix of `elastic2d != none` with `dof = trilinear` since the latter assumes 3d elasticity.
Also, do not assume that `flex_interp > 0` in the engine. This will enable to use, e.g., `flex_interp = -1` to mean a linear surface finite element instead of a 3d finite element which is currently identified with `flex_interp = 1`.
PiperOrigin-RevId: 903852035
Change-Id: Ia6290b4a05e9e510ffb7f36d141cd525b40d3110
This change introduces a `flex_cellcount` field to `mjModel` to specify the number of cells in each dimension for interpolated flexes. The stiffness computation, passive force calculation, and Jacobian derivatives are updated to operate on a per-cell basis, significantly improving performance by localizing computations to the nodes within each cell.
PiperOrigin-RevId: 901216393
Change-Id: Ic23132e609de11e71bb7fef8d1f139daad2ec264
The `vertcollide` attribute in flex and flexcomp contact specifications has been removed since the next mjWarp release will support flex collisions. This attribute was previously used to enable vertex-based collisions by adding sphere geoms at flex vertices. The functionality associated with this attribute has been deprecated and removed from the XML parser, documentation, and internal data structures. Affected example XML files have also been updated.
PiperOrigin-RevId: 885595503
Change-Id: I6a78975f1bc540c09ea8b5306ff3be0962b7ab24
The flex-SDF algorithm uses the same logic as mesh-SDF. Previous to this change, the flex was instead colliding with the convex hull of the mesh generated automatically with marching cubes from the SDF.
PiperOrigin-RevId: 876420302
Change-Id: I3bf6552b6b8e108ba4c4209bbcf10fe5f2b2fd90
This CL rolls back the change from changelist 794576123 since the previous version was more robust and stable. Compared to the original algorithm, this changes handles degenerate and very large triangles as well.
BEGIN_PUBLIC
Restore mesh-SDF collisions.
END_PUBIC
PiperOrigin-RevId: 875716983
Change-Id: If0de41dc7bbd60460779821d2338cf2d661f40d5
Temporary changes to the octree:
- Changed frame from mesh to geom
Before change (tolerance 1e-3):
```
Simulation time : 2.34 s
Steps per second : 4275
Realtime factor : 8.55 x
Time per step : 233.9 µs
Newton iters / step : 2.26
Contacts / step : 3.69
Constraints / step : 14.76
Degrees of freedom : 12
Dynamic memory usage : 0.7% of 14M
```
After change (6 octree levels):
```
Simulation time : 2.24 s
Steps per second : 4472
Realtime factor : 8.94 x
Time per step : 223.6 µs
Newton iters / step : 2.47
Contacts / step : 3.37
Constraints / step : 13.49
Degrees of freedom : 12
Dynamic memory usage : 0.2% of 14M
```
PiperOrigin-RevId: 781019754
Change-Id: Ib15581244dfe9e571e1c6c2cac4a101cf3d8ba3d
chenlin
Copybara-Service