Declare the environment as spatial fields over convex cells, to lift the single-medium restriction

[teerthsharma]

Declare the environment as spatial fields over convex cells, to lift the single-medium restriction

The feature, motivation and pitch

The entire specification of the physical environment in MuJoCo is five global constants: opt.gravity (one 3-vector), opt.density (one scalar), opt.viscosity (one scalar), opt.wind (one 3-vector), opt.magnetic (one 3-vector). There is no spatial field of any kind and no per-region override; body_gravcomp is a per-body scale on the single global vector, not a field. A simulator whose environment is a constant represents exactly one medium, and cannot represent a boundary between two media, because a boundary is the one thing a constant cannot have.

This proposal replaces the constants with fields declared over a decomposition of space into convex cells, each cell carrying its own gravity, density, viscosity and (optionally) pressure. A geom is clipped against the cells it meets and the force is the sum of per-cell contributions computed on the clipped pieces. One scene then contains vacuum, air and water simultaneously; a body may span all three; and buoyancy, partial submersion, ballast response, and gravity-gradient attitude dynamics all follow from the declaration rather than from per-body attributes. A regularity theorem (below) guarantees the resulting force is C¹ in pose except at a finite, compile-time-computable set of configurations, so the boundary does not introduce chatter. All measurements below are against MuJoCo 3.11.0 (Python bindings) with model sources from checkout 747d36da, on Windows 11 x86-64; the analytic reference where one is needed is scipy 1.17.1 DOP853 at rtol 1e-13, atol 1e-15.

The six states a body can be in

A single sphere of radius 0.15 m at exactly the fluid density — mass 14.137167 kg, displacing 14.137167 kg — passed through every configuration a robot encounters. The buoyant force column is ρ_f·g·V_sub evaluated in closed form; the MuJoCo column is measured.

state	body centre relative to surface	V_sub [m³]	physics	MuJoCo 3.11.0
normal, in air	+0.40 m	0.0000000	0.000 N	0 N
entering, breaking the surface	+0.10 m	0.0010472	10.273 N	0 N
semi-submerged	0.00 m	0.0070686	69.343 N	0 N
submerged	−0.40 m	0.0141372	138.686 N	0 N
escaping, surfacing	−0.10 m	0.0130900	128.413 N	0 N
space-faring, vacuum	—	—	gravity-gradient torque up to 1.593 N·m	0 N·m

The four interior rows — the entire transit through the interface, which is where every amphibious, diving, surfacing and breaching manoeuvre happens — are not approximated badly. They do not exist. A neutrally buoyant body, which should hover, falls at −9.810000 m/s² at every one of these heights, and reaches −8.6849 m after 2 s from a start at −5.0 m.

The lower panel is the whole proposal in one curve: the physical force falls smoothly from full displacement to zero as the body rises through the interface, MuJoCo is the flat line at zero, and the two dashed verticals — the heights at which the body is exactly fully submerged and exactly fully emerged — are the only places the slope can jump. That last claim is proved below and holds for every convex shape.

Problem statement

There is no interface anywhere in the engine. A sphere of radius 0.1 m with fluidshape="ellipsoid" in a model with density="1000" returns a vertical fluid force of +17.592919 N with its centre at z = +0.5, at +0.05, at 0.0, at −0.05 and at −0.5 — the same number above, at, and below any surface one might draw. Searching src/ for buoyan, hydrostat, archimed, submerg, waterline, free surface returns two hits, both in experimental USD schema files, neither in physics code. The fluid model computes blunt and slender drag, angular drag, Kutta lift and Magnus lift; it does not compute the force that makes things float.

A perfectly neutrally buoyant body falls at g. A sphere of radius 0.15 m at exactly the fluid density (mass 14.137167 kg, displaced mass 14.137167 kg) at rest, fully submerged: qfrc_fluid = [0,0,0,0,0,0], qfrc_passive = [0,0,0,0,0,0], vertical acceleration −9.810000 m/s² against a physical 0. After 2 s it has fallen from −5.0 m to −8.6849 m.

The sanctioned proxy cancels its own moment. The changelog describes gravcomp as implementing "gravity compensation and buoyancy" and balloons.xml ships as the demonstration. What the mechanism computes is decided by two lines of engine_passive.c: the force is −body_mass · gravcomp · gravity, applied at xipos — the centre of mass — with a zero torque argument. A force through the centre of mass produces no moment about the centre of mass, so the restoring couple ρgV·BG·sin θ that stabilizes every real vessel is structurally zero. Measured on a deliberately self-righting hull (ellipsoid with a dense keel 0.30 m below centre, gravcomp="1"):

roll	gravcomp generalized torque	angular acceleration	physical angular acceleration
5°	0.171908 N·m	0	0.121878 rad/s²
15°	0.510499 N·m	0	0.361837 rad/s²
30°	0.986209 N·m	0	0.699093 rad/s²
45°	1.394710 N·m	0	1.118721 rad/s²

The torque column is nonzero — which makes the failure easy to miss — and the response column is zero because gravity's torque cancels it identically at every attitude. Across a 0–90° sweep in 2.5° steps the largest angular acceleration returned anywhere is 1.56e-15 rad/s². Released at 30° roll and integrated 20 s: still 30.0000°, |ω| = 0.000e+00. The hull is in neutral equilibrium at every attitude and can never right itself.

Ballast does nothing. Buoyancy is proportional to displaced volume; the proxy is proportional to mass; so the net force is (γ−1)·m·g and mass cancels out of the trim. A 0.30 m float trimmed neutral and then ballasted:

ballast taken on	net vertical force, MuJoCo	net vertical force, physics
0.5 kg	0.000000 N	−4.905000 N
1.0 kg	0.000000 N	−9.810000 N
3.0 kg	0.000000 N	−29.430000 N

Variable-buoyancy depth control — how essentially every AUV regulates depth — has no plant to act on. The workaround, rewriting gravcomp = ρV/m on every mass change, lives in mjModel rather than mjData, so it is outside the control path, invisible to inverse dynamics, and shared across a batch where one model backs many environments.

The added-mass acceleration terms are documented, implemented, and never executed. doc/computation/fluid.rst:338-390 derives them in full from Lamb (1932); the compiler computes and stores the virtual masses and inertias per geom (geom_fluid[6..11]; measured for a 5:1 prolate hull: virtual_mass = [0.001238, 0.018729, 0.018729]); and the sole call site, engine_passive.c:1249, passes NULL for the acceleration argument under the in-tree comment "disabled due to dependency on qacc but included for completeness". The branch has been dead since March 2023. Measured consequences, impulsive force on a body at rest, fluid density 1000:

case	m [kg]	m_added [kg]	a MuJoCo	a physical	error
sphere, neutrally buoyant	2.0944	2.0944	4.77465	2.38732	100.0%
prolate 5:1, axial	20.9440	1.2382	0.47746	0.45081	5.9%
prolate 5:1, lateral	20.9440	18.7294	0.47746	0.25206	89.4%
prolate 5:1, lateral, ρ_body = 500	10.4720	18.7294	0.95493	0.34245	178.9%

The error is worst exactly where underwater robotics lives: slender, neutrally buoyant hulls in sway. It is a model error, not an integration error: with all drag and lift coefficients zeroed, the body's own kinetic energy is conserved to machine precision while the physical invariant ½v′(M+M_A)v drifts by 43.060728%, 43.060765% and 43.060765% at timesteps 1e-3, 1e-4 and 1e-5 — identical to six significant figures across three decades, so no timestep refinement touches it. The mass matrix is bit-identical with fluid on and off (max |ΔM| = 0). ρM_A is diagonal and positive semidefinite on all six geom types checked (min eigenvalue ≥ 0), so it is safe to add to the effective metric.

The orbital case: a constant vector is not a weaker model, it is no model. A constant gravity field has zero gradient, and the gradient carries the attitude physics. The gravity-gradient torque τ = 3(μ/R³) R̂×(IR̂) — the passive mechanism that has held elongated satellites pointed at Earth since the 1960s — is identically zero for constant g, and measured as exactly [0,0,0] bias torque on an elongated body at tilts 0/15/30/45/60/90°. At ISS altitude (μ/R³ = 1.279502e-06 s⁻²):

body	principal moments [kg·m²]	max gravity-gradient torque	MuJoCo
6U CubeSat	0.020, 0.035, 0.040	3.839e-08 N·m	0
ESPA-class satellite	80, 120, 150	1.343e-04 N·m	0
servicing arm + payload	150, 900, 950	1.535e-03 N·m	0
ISS truss segment	1.2e5, 9.0e5, 9.5e5	1.593e+00 N·m	0

Tidal (differential) acceleration across an extended body, 2(μ/R³)L: 2.559e-05 m/s² across a 10 m servicing arm, 2.789e-04 across a 109 m truss, 2.559e-03 across a 1 km tether — all exactly zero under a constant vector. These are small numbers acting continuously in a fixed direction: the regime where a simulator that reports zero is misleading rather than imprecise, because the effect accumulates instead of averaging out.

neutral body sinking at g; keeled-hull restoring response, MuJoCo vs ρgV·BG·sinθ/I; ballast response

gravity-gradient torque physics vs MuJoCo; tidal acceleration vs body extent

Proposed design

Fields, declared once, over convex cells. Space is decomposed by a finite set of planes into convex cells; each cell carries (ρ, μ, ψ) freely and g under the admissibility restriction below. Nothing is written onto bodies: the water is declared once and every geom intersecting it receives the correct force from its own geometry. Adding a body adds no environment data; changing a mass changes no environment data; ballast becomes a mass change, which is what it physically is. The per-body count of hand-maintained environment numbers goes from n to zero.

Per-geom clipping. A geom fully inside one cell uses its compile-time volume and centroid — no clipping at all. A geom straddling a boundary is clipped, and the hydrostatic force is the cell's ρg times the clipped volume, applied at the clipped centroid:

• Sphere: spherical cap, V = πh²(3r−h)/3, centroid in closed form, C¹ in h throughout.
• Ellipsoid: the affine map to the unit ball reduces it to the sphere; volume scales by the determinant, centroid maps forward.
• Box, cylinder, capsule: closed-form case enumerations (the box by vertex count below the plane; the cylinder via circular/elliptical segments; the capsule adds two spherical caps).
• Mesh: clip the convex hull the compiler already builds; faces below the plane contribute directly, crossing faces are clipped, the cut edges assemble into the single closed loop guaranteed by convexity, and the cap closes the cycle. Volume and centroid then follow from the divergence-theorem sums the compiler already uses for mesh inertia (V = ⅓∮x·n dA, V·c_i = ½∮x_i²n_i dA). The closed-cycle condition ∂(∂Ω_sub) = 0 is the correctness condition of the algorithm: an implementation that clips triangles independently and forgets the cap returns plausible, silently wrong volumes.

Added inertia enters the metric, not the right-hand side. The second consumer of the same clipped-volume machinery is the dead added-mass branch: the correct treatment of ρM_A is M̃ = M + Σ_g J_gᵀ R_g (ρ_f M_A^(g)) R_gᵀ J_g, which is exactly the discretization shipped for flex in 3.11.0 (mjd_effBuild assembling M̃ = M + (h²+hd)K) and the shape of the two open proposals #3443 (actuator ∂τ/∂q) and #2433 (custom passive forces). Fluid added mass is a further PSD-safe force class from the follow-up list in #3386, and the only one with no proposal behind it. Applying it instead as an explicit force from the previous step's acceleration is the classical added-mass instability of partitioned fluid–structure schemes and diverges precisely when m_a/m approaches 1 — the measured ratio for a light slender hull is 1.789 — so the metric route is the only correct resurrection of the branch.

the effective-metric family: flex shipped, #3443 and #2433 open, fluid added mass unclaimed with a dead branch

Working implementation

A reference implementation of the single-boundary case exists and is measured below. It adds src/engine/engine_hydrostatic.c, one enable bit mjENBL_HYDROSTATIC (mjtype.h, mjNENABLE 6 → 7), and six lines in mj_passive that call it after mj_fluid so the force accumulates into qfrc_fluid. Built with MSVC 19.44.35224 x64, Release, Ninja, from 747d36da.

The core is one dispatch, mju_geomSubmerged(cen, geomtype, size, n, off), returning the volume below the plane dot(n,x) = off in the geom frame and writing the volume centroid. It has three tiers:

1. Support-value test, O(1). For a convex body the support value in direction n decides containment outright: sup <= off means fully submerged, -sup >= off means fully emerged. Both are closed form per primitive — size[0] for a sphere, |size·n| for an ellipsoid, Σ size_i |n_i| for a box. A fully submerged geom therefore needs no clipping at all and returns the compile-time volume with the centroid at the geom origin.
2. Closed form for smooth primitives. The sphere cap integrates directly, V = π(r²·off − off³/3 + 2r³/3) with first moment M = π(r²off²/2 − off⁴/4 − r⁴/4), so the centroid is M/V — no textbook special cases and no division by a vanishing height. The ellipsoid reduces to the sphere by the affine map that takes it to the unit ball: the plane pulls back to dot(A n, y) <= off, the volume scales by det A, and the centroid pushes forward. No new integration.
3. Polyhedral clip for the box and, by extension, hull meshes. Each face is clipped by the half-space, the cut points are collected, ordered into a loop about their centroid in the plane, and appended as the cap. The boundary is then integrated by the divergence theorem, V = ⅙ Σ dot(p₀, (p₁−p₀)×(p₂−p₀)) over the fan triangles, with the first moment from the same cross products. The cap is what makes the boundary a closed cycle; omitting it returns a plausible, wrong number.

Correctness, 13 of 13 checks at machine precision.

check	got	want	rel. error
sphere r=1 fully submerged	4.1887902	4π/3	0
sphere r=1 exactly half	2.0943951	2π/3	0
sphere half, centroid depth	−0.375	−3/8	0
sphere cap h=0.5	0.654498469	πh²(3r−h)/3	1.7e-16
box 2×2×2 fully submerged	8	8	0
box 2×2×2 exactly half	4	4	0
box half, centroid depth	−0.5	−0.5	0
box 2×2×2 quarter	2	2	3.33e-16
box, plane through centre along (1,1,1)/√3	4	4	0
ellipsoid 1×2×3 fully submerged	25.1327412	4π·6/3	0
ellipsoid 1×2×3 exactly half	12.5663706	2π·6/3	0

A 2,000-point sweep of a box cut by a plane along (1,1,1)/√3 from fully emerged to fully submerged returns 0 non-monotone steps, consistent with the regularity result below.

The physics that was previously zero. Keeled hull, ellipsoid 0.5×0.1×0.1 with a 20000 kg/m³ keel sphere 0.30 m below centre, released at roll θ, mj_forward:

roll	angular acceleration, flag on	angular acceleration, MuJoCo today
0°	0 (equilibrium)	0
5°	−1.11183 rad/s²	0
15°	−3.30171 rad/s²	0
30°	−6.37841 rad/s²	0
45°	−9.02043 rad/s²	0

Restoring, monotone in roll, and zero at zero — the qualitative behaviour a hull must have and did not.

Ballast response, 0.30 m float, net vertical force:

ballast	flag on	closed form
1.0 kg	−9.768908 N	−9.810000 N
2.0 kg	−19.578908 N	−19.620000 N
3.0 kg	−29.388908 N	−29.430000 N

Every row differs from the naive closed form by +0.041092 N exactly, which is ρ g · (4/3)π(0.01)³ — the buoyancy of the small marker geom carrying the ballast mass, which the implementation correctly counts and the hand-written expectation omitted. Corrected for that, the agreement is exact.

Cost. Best of 5 × 2000 steps, single thread, same binary, flag toggled at runtime. On models with no fluid the path early-exits on opt.density <= 0:

model	nv	flag off	flag on	ratio
humanoid.xml	27	0.01735	0.01774	1.0225
humanoid100.xml	627	0.36930	0.37842	1.0247
cube_3x3x3.xml	66	0.19983	0.19858	0.9937
hammock.xml	312	0.24159	0.24174	1.0006
balloons.xml	30	0.04137	0.04136	0.9998
car.xml	8	0.00459	0.00473	1.0313

Two of six measure faster with the flag on, so the spread is run-to-run noise and the early exit is free.

When the path does run — fully submerged scenes with the existing fluid model active, best of 7 × 1500 steps:

geom	n	nv	flag off	flag on	ratio	µs per geom
sphere	50	300	0.01596	0.01934	1.212	0.0677
sphere	200	1200	0.07537	0.09007	1.195	0.0735
sphere + ellipsoid fluid	200	1200	0.08327	0.09654	1.159	0.0664
box	50	300	0.01608	0.02255	1.402	0.1294
box	200	1200	0.08282	0.10402	1.256	0.1060
ellipsoid	200	1200	0.08973	0.10436	1.163	0.0731

The marginal cost agrees between n=50 and n=200 for every geom type, so it is linear in geom count as the construction requires. The support-value fast accept is worth 48% on the box path — 0.2049 → 0.1060 µs per geom at n=200 — because a fully submerged box then skips the clip entirely. The residual 0.066–0.073 µs per geom for the closed-form primitives is dominated by mj_applyFT, not by the volume computation: applying a force at a point costs what it costs, and the existing fluid model pays the same per geom.

Regularity: the boundary does not chatter

A field means bodies straddle boundaries; straddling means clipping; and a restoring force whose slope jumps is the standard ingredient of waterline limit cycles. For convex geoms the question closes completely.

Theorem (slice regularity). Let K be a compact convex body in R³ with nonempty interior and n a unit vector. Put t_min = min over K of x·n and t_max = max; let A(t) be the area of the section K∩{x·n = t}, V(t) the volume of the cap K∩{x·n ≤ t}, and F = K∩{x·n = t_max} the extremal face. Then (a) V is the integral of A and V′ = A at every continuity point of A; (b) A^1/2 is concave on [t_min, t_max] — for t = λt₁+(1−λ)t₂, convexity gives λK_{t₁} + (1−λ)K_{t₂} ⊆ K_t and planar Brunn–Minkowski gives A(t)^1/2 ≥ λA(t₁)^1/2 + (1−λ)A(t₂)^1/2 — hence A is continuous on the open interval and V is C¹ there; (c) the one-sided limit of A at t_max is area(F), so V′ jumps there by area(F), zero iff F is a point or a segment; (d) consequently V′ can fail to be continuous only at t_min and t_max.

The count is topological: for smooth strictly convex K the height function on the boundary 2-sphere is Morse with exactly two critical points, whose critical values are t_min and t_max; for a polytope the two critical values persist and the jump magnitude is the extremal face area, which is what (c) computes. The buoyant force on any convex geom is therefore C¹ through the entire partially-submerged regime except at exactly two heights per (geom, plane) pair — the support values in the plane normal direction, computable at compile time from the support function, with no tolerance and no configuration-dependent case analysis. The exceptional set in configuration space is a finite union of codimension-one sets sized by the number of (geom, plane) incidences.

Verified by direct section-area measurement on eight convex bodies (cube along face normal / body diagonal / edge direction, tetrahedron along face normal / generic, random hulls of 12 and 40 points, 400-point sphere): zero interior discontinuities in every case; the cube along a face normal gives interior variation exactly 0.00000 with endpoint limits 4.00000 = the face area; the tetrahedron along a face normal gives 3.46410 = √12 = the equilateral face area at one end and 0 at the vertex end, as (c) predicts.

A(t) and V(t) for four bodies: continuous inside, jumps only at ends with extremal faces

dV/dz is the waterplane area; it jumps exactly at the two critical heights for the box and upright cylinder, and nowhere for the sphere and horizontal cylinder

Algorithm and cost

Clipping is a flood fill on the hull graph, and it is provably correct. For a convex polytope and any threshold t, the cap vertex set {v : v·n > t} induces a connected subgraph of the edge graph — the monotone-path property (from any vertex, an edge path exists along which a linear functional strictly increases to the maximiser; such a path never descends, so it never leaves the cap). Three consequences: the search starts at the support vertex the existing hill-climb (mjc_hillclimbSupport) already computes; the cost is output-sensitive (a barely-submerged geom costs almost nothing); and the result is independent of the starting vertex, so warm starting cannot change the answer. Tested directly: 0 disconnected caps and 0 hill-climb failures in 4,200 cuts over 420 directions across seven bodies including a flat slab and a needle, the aspect ratios where such properties fail first.

The machinery already exists, from two recent contributions. The disjoint-set island rewrite (#3396, merged) exports mj_dsuMerge / mj_dsuRoot / mj_dsuAssign / mj_floodFill through the private MJAPI, and its canonical minimum-index roots give the (geom, cell) incidence decomposition a deterministic accumulation order — which is what makes the per-cell force sum bitwise reproducible rather than merely correct in exact arithmetic. The hull-graph rewrite (#3450, open) builds the graph the flood fill walks in linear time (probes 3V²−6V → 8V−12, 4.21x end-to-end at V = 40,962), which is the precondition for walking it every step rather than once at load. The dependency is real in both directions: the clipping cannot be written without the exported component helpers and cannot be afforded without the linear construction.

the dependency chain; the same primitive, connected components, at two scales: kinematic trees and hull vertices

The per-step budget is exact, not asymptotic. The planes that matter for a geom are the k that cross it, and k planes cut a convex body into at most R(k) = C(k,0)+C(k,1)+C(k,2)+C(k,3) = (k³+5k+6)/6 pieces — verified by sign-vector enumeration over 400,000 ball samples: 1, 2, 4, 8, 15, 26, 42, attained in general position, strictly fewer (14, 22, 32) for concurrent planes. The whole-scene budget is Σ_g R(k_g)·O(3V_g−6), the second factor by Euler's formula on the hull — the identity #3450 already uses as its cost model. Realistic declarations are tiny: one free surface is k = 1 (2 pieces); surface plus a wall is k = 2 (4); the environment complex for an air–water scene has two cells and one shared face. It is a declaration, not a discretization, and does not grow with scene size or resolution.

Admissibility: not every field assignment is physical

Density, viscosity and pressure may be piecewise constant freely — they are not forces derived from a potential, and a discontinuous ρ at an air–water interface is exactly the point. Gravity may not. For piecewise-constant g the potential in cell i is U(x) = −g_i·x + c_i, and continuity of U across a shared face with unit normal n_ij requires

g_j − g_i ∥ n_ij for every adjacent pair — only the normal component of gravity may jump.

Measured by integrating the work around a closed rectangular loop crossing a boundary: Earth-g and Moon-g side by side across a vertical plane (the configuration a user writes first) yields +8.190000 J/kg per circuit — free energy, unbounded, a robot cycling the boundary accelerates forever; the same gravity difference across a horizontal face yields −2.8e-11 J/kg (quadrature error); a difference parallel to the vertical face's normal yields exactly 0. The recommended resolution is not the face-by-face compile check but the stronger restriction: gravity fields limited to smooth analytic forms — uniform, central −μx/|x|³, or linear in position — which are conservative by construction, cover both the planetary-surface and orbital cases, and cannot be misconfigured. A central field also supplies the gravity-gradient torque and tidal terms of the tables above as consequences rather than features.

Measured summary

check	result
repeat determinism, shipped corpus	82/82 models bitwise identical, 3 × 200 steps; 0 failed to compile
hydrostatic term in engine sources	0 occurrences in physics code
neutrally buoyant body, vertical acceleration	−9.810000 m/s² (physical: 0)
keeled hull, max angular accel, 0–90° sweep	1.56e-15 rad/s² (physical at 45°: 1.118721)
keeled hull released at 30°, after 20 s	30.0000°, |ω| = 0.000e+00
net force from 3 kg ballast	0.000000 N (physical: −29.430000 N)
depth dependence of fluid force	identical at −1000 m and +100 m
added-mass lateral acceleration error	89.4% to 178.9%
added-mass energy-invariant violation	43.060728 / 43.060765 / 43.060765 % at dt = 1e-3/1e-4/1e-5
ρM_A positive semidefinite, 6 geom types	min eigenvalue ≥ 0 in all cases
gravity-gradient bias torque, elongated body	exactly [0,0,0] at 6 tilts (physical: up to 1.593 N·m)
slice-regularity: interior discontinuities	0 across 8 bodies; endpoint jumps = extremal face areas to 5 digits
cap connectivity	0 disconnected caps, 0 hill-climb failures, 4,200 cuts, 420 directions
pieces from k cutting planes	1, 2, 4, 8, 15, 26, 42 = R(k) exactly, general position
non-conservative side-by-side gravity	+8.190000 J/kg per circuit, matching (Δg_z)·h·w/2 in closed form
baseline cost for the benchmark comparison	humanoid.xml 0.0177 ms/step (nv=27); humanoid100.xml 0.3086 (nv=627); cube_3x3x3.xml 0.1776 (nv=66); hammock.xml 0.2120 (nv=312), single thread, mean of 3 × 1000 steps

Questions

1. Framing. The changelog describes gravcomp as implementing buoyancy and balloons.xml ships as the demonstration — the one case where the proxy is nearly right, since a balloon's centres of buoyancy and mass nearly coincide. Is the position that gravcomp is the intended buoyancy mechanism with the centre-of-mass application point a known trade, or is a real hydrostatic force welcome? The answer decides whether this is a feature proposal or a documentation correction.
2. The dead branch. Was the NULL at engine_passive.c:1249 intended as permanent, or as a placeholder pending exactly the effective-metric machinery that has since shipped for flex (Implicit flex elasticity in the CG constraint solver via an effective metric #3386) and is proposed for actuators (Make actuator kp implicit (∂τ/∂q) to lift the stiffness cap #3443) and passive forces (Mujoco (implicit) integration does not offer a mechanism to take into account custom passive forces in qDeriv calculation #2433)? If the latter, the term is implementable as originally written, as one more PSD-safe force class behind the general-effective-metric flag.
3. Environment specification. Is a first-class environment field object (cells + per-cell parameters in mjModel, XML-declared) acceptable, or should the first landing be the degenerate single-boundary case — one opt-level free-surface plane with (ρ, μ) above and below — which already covers every underwater and amphibious use case in the tracker (Buoyancy force problem with fluid dynamic #3066, Modelling Shallow Water Correctly #1984, How to implement fin/thruster dynamics #2655) with no new complex machinery?
4. Gravity forms. Restricting gravity to uniform / central / linear (conservative by construction) versus enforcing the tangential-continuity check on arbitrary per-cell constants: is ruling out the piecewise-constant form acceptable? The measured +8.19 J/kg free-energy loop is the argument that it should be.
5. Staging. Five separable stages, each useful alone: (1) the added-mass metric term; (2) compile-time geom_volume/geom_centroid, no dynamics change; (3) fully-submerged hydrostatics (force at the true centre of buoyancy — most of the value, none of the clipping risk); (4) the free surface and partial submersion (the clipping and the two-critical-heights machinery); (5) the general cell complex. Is this ordering reviewable?

Alternatives

Document gravcomp's limits and stop. Fixes discoverability; fixes nothing else. No documentation makes a force through the centre of mass produce a righting moment, or a mass-proportional force respond to ballast.

User callbacks via mj_applyFT. What the unanswered subsea user (#2655) eventually did. Puts submerged-volume computation on every user independently, outside inverse dynamics, invisible to sensors, absent from MJX/Warp. The geometry lives in the engine; the computation belongs beside it.

A plugin. Viable for the hydrostatic force alone; fails for added mass, which must enter the mass matrix and factorization; and splits the shared volume/centroid machinery across a boundary.

Depth-spring proxies. Reproduce neither the attitude dependence nor the volume dependence, and introduce a stiffness that must be tuned against the timestep — trading a missing force for a stability constraint.

Full CFD. Out of scope, and the fluid documentation says so. Nothing here solves a flow field: hydrostatics needs the displaced volume; added mass is a compile-time constant of the geometry.

Explicit added-mass force from the previous step's qacc. The classical partitioned-scheme instability; diverges as m_a/m → 1, the exact regime of a neutrally buoyant vehicle (measured ratio 1.789). The metric or nothing.

Additional context

A full-history search (2,551 issues and PRs plus 830 discussions; 3,381 of 3,450 numbers recovered; every empty result re-verified against the comment-inclusive GitHub index) returns zero hits for: hydrostatic, archimedes, centre/center of buoyancy, metacentric, righting moment, submerged, waterline, free surface, neutrally buoyant, displaced volume, virtual mass, Munk moment, potential flow, acceleration reaction, spatially varying gravity, gravity field, density field, submerged volume, half-space clipping. The only issue with buoyancy in its title (#3066) reports a version regression inside the gravcomp path and does not question the proxy; #1984 (amphibious, "only half the body should be in water") and #2655 (subsea, buoyancy/thrusters/fins) describe the missing capability and received no maintainer reply. No maintainer has declined buoyancy as out of scope in any thread found.

Estimated implementation, conservatively: ~2,100 lines including ~670 of tests and benchmarks — closed-form clipped volumes and centroids per primitive (sphere/ellipsoid ~160, box ~150, cylinder/capsule ~130), hull-mesh clip with cut-loop assembly (~170), compile-time volume/centroid (~85), force assembly and free-surface plumbing (~195), added inertia in the metric with PSD guard (~180), inverse mirror (~60), model fields and XML I/O (~160), bindings tables (~55), tests (closed-form vs quadrature, righting moment vs ρgV·BG·sinθ, C¹ continuity away from the two critical heights, flag-off bit-identity across the 82-model corpus, PSD and energy-invariant checks, ~670), benchmark entry (~50), documentation (~120). The hard parts are the box case enumeration, the mesh cut-loop assembly (where the closed-cycle condition must be enforced rather than assumed), and the near-tangent numerics.

Validation on the reference implementation: the full CTest suite passes 1,339 of 1,339 entries, 0 failed, built with MSVC 19.44.35224 x64 Release from 747d36da with the new file and enable bit in place. One entry, doc_test, failed on the first run because doc/includes/references.h had not been regenerated for the new enable bit — the suite catches exactly the documentation-plumbing omission listed in the size estimate above, and passes once that file is updated. Two entries are skipped by the harness on this platform (RecompileCompareTest, WriteReadCompareTest), the same two as on an unmodified checkout. Still outstanding before review: ASan and UBSan over the new tests — UBSan matters here because the clipping divides by face normals that can be near-parallel to the plane — and per-model trajectory hashes against the parent commit to confirm flag-off bit-identity across the 82-model corpus.

Limits

Scope is hydrostatics and added inertia in piecewise environments with a static free surface: waves, currents, sloshing, surface tension, wave-making resistance, compressibility and depth-varying density are all out of scope. Concave meshes are clipped via their convex hulls, which overestimates displaced volume for concave shapes — exact for convex meshes, an explicit documented approximation otherwise. The buoyant force on a box or an upright cylinder is C⁰ but not C¹ at the two heights where a flat face lies in the surface (measured waterplane-area steps 0.1714 m² and 0.0898 m² on the tested geometries); the theorem guarantees these are the only such heights and locates them at compile time, but they cannot be removed without changing the geometry. The added-mass half changes the trajectories of existing models that use the ellipsoid fluid shape and cannot claim no-regression; it belongs behind its own flag, presented as making the implementation match its own documentation. The added-mass model itself is the potential-flow ellipsoid approximation, exact only for an ellipsoid in an unbounded inviscid fluid, and MuJoCo's existing inscribed-ellipsoid approximation (−47.6% volume for a cube) becomes more visible once forces depend on volume. The corpus argument is necessarily weaker than for an accuracy fix — no shipped model uses a capability that does not exist — so the case rests on the application class (#3066, #1984, #2655) rather than on existing models. MJX and Warp would initially reject the new fields, a parity gap of a kind the project already carries. The per-step cost claims are bounds and baselines, not flag-on measurements; the benchmark comparison is committed to, not asserted.

Relates to #3386, #3443, #2433, #3066, #1984, #2655, #3396, #3450.

[teerthsharma]

Update: #3450 merged on 3 August, so two of the three pieces this proposal builds on are now in main — the disjoint-set component helpers from #3396, and the linear hull-graph construction from #3450. The cost model in the section above is therefore stated against merged code rather than against pending work, and the figure showing #3450 as open is out of date as of today updating diagram shortly.

Thanks @yuvaltassa — the turnaround was faster than anticipated.

[teerthsharma]

Draft implementation of the field half, split into two stacked branches so the data and the reader can be reviewed apart.

• Add environment layers to mjModel, MJCF and the spec API #3460 — nlayer and the layer_* arrays in mjModel, the <option><layer> element, and mjs_addLayer/mjs_getLayer/mjs_deleteLayers. Nothing reads them yet; a model with no <layer> compiles to nlayer == 0 and is unchanged. 1977/1977 on Linux.
• Sample the environment layer field in the fluid models #3461 — mj_envLayer and mj_envSample, with the two fluid models and their derivatives sampling per geom instead of reading mjOption. 1991/1991, and 14/14 for the field tests under ASan, UBSan and LeakSanitizer.

Both leave gravity global. It does not enter as a force at a point: mj_rne sets the world frame's fictitious acceleration to -g and propagates it down the tree, so each body picks it up through cinert * cacc, an integral over its mass distribution with no sample point to name. mj_gravcomp must cancel that term exactly and mj_energyPos must match it, so sampling any one of them at xipos would desynchronise the others for every body spanning a boundary. A gravity field needs the integral, not a sample.

Neither branch carries the buoyant force from the issue body. That depends on clipped volumes, which is where the harder geometry sits, and it is left for a later change once the field itself is settled.

Two measurements from the branches that bear on the argument above. Step time with nlayer == 0 has a geometric mean of 0.9999 against main across 13 models, inside the envelope of main measured against a second build of itself (0.9981), with bit-identical qpos+qvel checksums on all 13 — so the fallback is the old path rather than an approximation of it. mj_sizeModel grows by exactly 8 bytes, the nlayer counter, and nbuffer is byte-identical, so the arrays cost nothing unused. When layers are declared the lookup costs +0.88 ns per binary-search step.

Both branches need #3458 before they can be built under a sanitizer: two constructs in mjsan.h stop any ASan or UBSan build of MuJoCo compiling under GCC 15.2.0 or clang 21.1.8. The ASan figures above were produced with that fix applied locally.