The actual buoyancy math

Last time the DLL just returned a hardcoded 100N upward force. Time to replace that with the actual math, Hirae's closed-form integrator from the SIGGRAPH Asia 2025 paper. Now it loops through every triangle, figures out if it's above water, below water, or crossing the surface, calculates the pressure force and torque on each one, and adds them all up.

The algorithm

For each triangle, compute the area vector (one half-cross-product gives you both the area and the outward normal), then plug the vertex coordinates into closed-form equations for force and torque. No numerical integration, no quadrature error, exact up to floating point precision.

The force equation (Hirae Eq. 5) is way simpler than you'd expect:

    glm::vec3 areaVec = 0.5f * glm::cross(b - a, c - a);
    float sumY = a.y + b.y + c.y;

    float part1 = (rho * gravity / 3.0f);
    float scalar = part1 * (-3.0f * h + sumY);
    totalForce += scalar * areaVec;

h is the water height at the triangle's centroid. For flat water that's 0, so the force depends on how deep the vertices are, deeper stuff gets pushed harder.

The torque equation is where it gets ugly (Eqs. 6-9). You need area and normal separately, so split the area vector with glm::length(), to get "area"/length and then we can just divide with that to get the normal. Makes it just one sqrt per triangle. Then compute an auxiliary vector A with three components (sum-of-pairwise-products of vertex coordinates), and cross it with the normal. Looks intimidating on paper but it's just arithmetic, no trig, no loops, one pass and done.

    float S = glm::length(areaVec);
    if (S < 1e-12f) return;
    glm::vec3 n = areaVec / S;

    float sumX = a.x + b.x + c.x;
    float sumZ = a.z + b.z + c.z;

    float deltaXZ = sumY - 4.0f * h; // This was wrong in the paper, we'll get to this
    float deltaY  = 2.0f * sumY - 4.0f * h;

    glm::vec3 A;
    A.x = sumX * deltaXZ + (a.x * a.y + b.x * b.y + c.x * c.y);
    A.y = sumY * deltaY  - 2.0f * (a.y * b.y + b.y * c.y + c.y * a.y);
    A.z = sumZ * deltaXZ + (a.z * a.y + b.z * b.y + c.z * c.y);

    glm::vec3 innerPart = ((12.0f * h - 4.0f * sumY) * com) + A;
    float firstPart = part1 * (S / 4.0f);
    totalTorque += firstPart * glm::cross(innerPart, n);

The payoff, exact torque regardless of mesh resolution. The naive approach (sample pressure at the centroid) creates ghost torques on coarse meshes, according to Hirae's, this as told previously leads to having to have very polygon heavy meshes for big objects. Hirae's method matches the analytical solution from Igarashi and Nakamura with as few as 24 triangles.

Clipping

Triangles that cross the water surface can't be ignored, they contribute most of the restoring torque when an object tilts. Without clipping, a tilted cube gets force only from fully submerged faces and spins endlessly because there's no counter-torque from the partially submerged faces. We have two veresions of this that we care about:

• One vertex below, two above, submerged region is a triangle (the bottom vertex plus two intersection points).
• Two vertices below, one above, submerged region is a quad, split into two triangles.

We find intersection points by linear interpolation along the edges using signed distance to the water surface:

float dAbove = above.y - hAbove;
float dBelow = below.y - hBelow;
float t = dAbove / (dAbove - dBelow);
glm::vec3 p = glm::mix(above, below, t);

Both clipped cases feed the sub-triangles into the same accumulateBuoyancy function. No code duplication.

Bitmask classification

Classify a triangle by packing both the count and which vertices are submerged into a single integer. Each vertex gets a bit: a=bit 2, b=bit 1, c=bit 0. Result is a 3-bit mask from 0 (dry) to 7 (fully submerged). Small epsilon (0.1mm, 10^-4) prevents vertices exactly on the waterline from creating degenerate zero-area triangles.

int belowcount = ((a.y < -EPS) << 2) | ((b.y < -EPS) << 1) | (c.y < -EPS);

Switch on the mask to dispatch to the right clipping function. No if-else chains the bit pattern tells you which vertex is the odd one out.

switch (belowcount)
{
    case 0b000: continue;
    case 0b111: accumulateBuoyancy(a, b, c, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b100: oneBelowSplit(a, b, c, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b010: oneBelowSplit(b, c, a, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b001: oneBelowSplit(c, a, b, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b011: oneAboveSplit(a, b, c, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b101: oneAboveSplit(b, c, a, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
    case 0b110: oneAboveSplit(c, a, b, 0.0f, 0.0f, 0.0f, comWorld, obj->rho, obj->gravity, totalForce, totalTorque); break;
}

Things that went wrong

Triangle loop stepping bug, I had i++ instead of i += 3. Created overlapping triangles from consecutive triplets processed [0,1,2], then [1,2,3], then [2,3,4], out of bounds on the last two. Garbage triangles, wrong forces.

Winding order bug in clip dispatch, a subtle one. When vertex b was the lone vertex below water (mask 0b010), I passed (b, a, c) to the clip function. Same when b was alone above (0b101). Single swap from (a, b, c) flips the area vector's direction (cross product is antisymmetric). Clipped sub-triangles produced forces pointing backward. Cube got horizontal drift and unwanted spin whenever it tilted, also stacked which led to a spinning mess. To fix this I used cyclic permutation (b, c, a) instead of swap (b, a, c). Cyclic rotations preserve winding order, vertices stay in the same sequence around the triangle, just starting at a different point. Swaps flip handedness. Both look like "reorder vertices" but only one keeps the normal pointing outward.

Winding-order bug from the $(a,b,c)$ clip dispatch.

COM in wrong coordinate space, vertices get transformed to world space through the model matrix every frame. The center of mass... didn't. I stored rb.centerOfMass (local space) at creation and passed it directly into the torque formula alongside world-space vertices. Moment arm for every triangle was computed relative to the wrong point, offset by however far the cube was from world origin. Constant torque bias that no amount of angularDampning could beat, reapplied every FixedUpdate. So instead I transformed the COM through the model matrix the same way as vertices, right after the vertex loop, then pass world-space COM to the torque formula. Sneaky bug because "cube that never stops rotating" looks identical to "no rotational damping yet" (which I also didn't have). Only caught it after cranking angularDampning to 50 and watching the spin refuse to die. Real damping scales with angular velocity, constant bias torque doesn't, so drag can't win.

THE BIGGEST PROBLEM

Paper equation as printed. Force balance is still right, but the extra factor of $2$ in $A_x$ and $A_z$ leaves a residual torque, so the cube never settles and can also settle toward the wrong orientation in some settings.

Typo in the paper's torque formula. Even after the COM fix, the cube drifted slightly in rotation that didn't match the paper's equilibrium. Spent a long time looking at the code until I went back through the equations.

The paper writes auxiliary vector $A$ in Eqs. 7-9 with $\{-4h + 2(y_1+y_2+y_3)\}$ applied uniformly to all three components. I'd implemented it exactly. But deriving the surface integral from scratch gives a different result for $A_x$ and $A_z$ than for $A_y$.

The integral of two linear functions $f, g$ over a triangle is

Plug $f = y - h$, $g = x - x_{\mathrm{com},x}$ and expand: the bracketed factor multiplying $\sum x_i$ comes out as $\left(\sum y_i - 4h\right)$, not $\left(2\sum y_i - 4h\right)$. Same for $z$. Only $A_y$ gets the doubled term, and that's because of the algebraic identity

that lets the paper rewrite

That trick works for $Y$ because $f$ and $g$ are the same variable; for $X\cdot Y$ and $Z\cdot Y$ there's no such identity. The "2" is just wrong there. Same sanity check as in the appendix: triangle $(0,0,0)$, $(1,0,0)$, $(0,1,0)$ with $\sum x_i = 1$, $\sum y_i = 1$, and $h = 0$. Then

while the paper's printed version gives

That's exactly the extra spurious $\sum x_i\sum y_i$ term showing up in the wrong place. Same error for $A_z$.

Why doesn't this blow up? The spurious term per triangle is $\left(\sum x_i\right)\left(\sum y_i\right)$ (or $\left(\sum z_i\right)\left(\sum y_i\right)$). For a perfectly axis-aligned, origin-centered mesh it almost sums to zero across the closed surface. Tilt the cube or move it off origin and the cancellation breaks. Constant residual torque every frame. Fix: one line, split $\Delta$ into

use the right one per component. After this and the COM fix, torque at equilibrium drops to floating-point noise instead of persistent drift.

Pretty sure this is a typo in the paper. The algorithm is right, the physics is right, just the printed formula has an extra "2" in two places. Full derivation in the appendix.

How it tests

Rebuilt the DLL, hit Play. The cube falls, hits the water, starts oscillating around an equilibrium that's not face-down. One edge sits lower, it rocks around that pose. Freaked me out at first. Spent an hour convinced there was still a bug. Force magnitude was right, torque was near zero at rest, math checked out. Why wouldn't it sit flat?

Turns out this is correct physics and I just didn't know about ship stability. A floating object is only stable if its metacenter (geometric property of the hull at that waterline) is above the center of mass. For a uniform cube:

• ρ/ρ_water around 0.1 or 0.9: cube sits very high or very deep, face-down is stable, floats parallel.
• ρ/ρ_water in middle range (~0.21 to ~0.79): face-down is unstable. Metacenter sits below COM, so any tilt creates torque that increases the tilt instead of restoring it. Cube rotates until it finds stable orientation, edge-down or corner-down for mid-range densities.

Hirae's paper uses density ratio 0.75 and quotes analytical equilibrium tilt as 26.565° from Igarashi & Nakamura (2007). My cube didn't settle at a specific angle because there's no rotational damping yet, just kept oscillating around equilibrium instead of settling. But the direction was right, clear non-zero tilt that the cube kept returning to.

Intuitive way to think about it: very light flat board floats face-down, very heavy block of the same shape floats face-down too, just deeper. Something in between? Water can't decide which face to push on more, equilibrium becomes an edge pointing down where buoyancy distribution is more forgiving. People have videos of wooden cubes in sinks doing this. So the integrator isn't just producing forces, it's correctly modeling the metacenter stability physics of the hull shape. Way cooler than expected from what was supposed to just be "make boat float."

Verified: buoyant force equals mg at rest, torque oscillates around zero (no constant bias), setting cube density to 100 kg/m³ makes it settle face-flat as predicted.

Tilting the cube from its settled position produces restoring torque from clipped triangles. Debug viz shows color-coded triangles (green dry, blue wet, orange clipped), yellow normal arrows, red/blue force/torque vectors.

Fully fixed buoyancy math, correct winding, center of mass in world space, and corrected $A_x/A_z$ terms. Everything works well.

What this validates

Hirae's closed-form integrator works. Force and torque computed per-triangle with no quadrature error. Clipping handles water surface crossing correctly. Bitmask dispatch is clean, two clip functions handle all six mixed cases symmetrically. Angular momentum damping comes next right now the cube oscillates because there's no rotational damping.

Next I'll be implementing Hirae's angular momentum damping to stabilize rotation.

Appendix: deriving the auxiliary vector A

Keeping the full thing here so future-me (or anyone else hitting the same typo) doesn't have to redo it.

Side-by-side: paper vs. corrected

To make the disagreement totally unambiguous, here are the paper's printed equations next to what the derivation actually gives. $x_i, y_i, z_i$ are the three triangle vertex coordinates and $h_i$ is the water height at the triangle's centroid.

As printed in Hirae et al. (2025), Eqs. 7–9:

$A_1 = (x_1+x_2+x_3)\{-4h_i + \mathbf{2}(y_1+y_2+y_3)\} + x_1 y_1 + x_2 y_2 + x_3 y_3$

$A_2 = (y_1+y_2+y_3)\{-4h_i + 2(y_1+y_2+y_3)\} - 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

$A_3 = (z_1+z_2+z_3)\{-4h_i + \mathbf{2}(y_1+y_2+y_3)\} + z_1 y_1 + z_2 y_2 + z_3 y_3$

Corrected (what the integral actually evaluates to):

$\boxed{A_1 = (x_1+x_2+x_3)\{-4h_i + (y_1+y_2+y_3)\} + x_1 y_1 + x_2 y_2 + x_3 y_3}$

$A_2 = (y_1+y_2+y_3)\{-4h_i + 2(y_1+y_2+y_3)\} - 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

$\boxed{A_3 = (z_1+z_2+z_3)\{-4h_i + (y_1+y_2+y_3)\} + z_1 y_1 + z_2 y_2 + z_3 y_3}$

The only change: the 2 in front of $(y_1+y_2+y_3)$ inside $A_1$ and $A_3$ should be a 1. $A_2$ is unchanged. Eq. 6 of the paper (the outer torque expression that uses $A$) is also fine, only $A_1$ and $A_3$ are off.

The rest of this appendix shows where the corrected forms come from and why $A_2$ is the only one that gets to keep the "2".

Setting up the integral

Per submerged triangle $T$ with vertices $\mathbf{v}_1, \mathbf{v}_2, \mathbf{v}_3$, area $S$, outward unit normal $\mathbf{n}$, the torque about the center of mass $\mathbf{x}_{\text{com}}$ is:

$\boldsymbol{\tau} = \iint_T (\mathbf{x} - \mathbf{x}_{\text{com}}) \times (-p\,\mathbf{n})\, dS$

with hydrostatic pressure $p = \rho g(h - y)$. Since $\mathbf{n}$ is constant on the triangle, pull it out:

$\boldsymbol{\tau} = \rho g \left[\,\iint_T (y - h)(\mathbf{x} - \mathbf{x}_{\text{com}})\, dS\,\right] \times \mathbf{n}$

So the whole derivation reduces to evaluating one component-wise integral:

$I_k = \iint_T (y - h)(x_{k,i} - x_{\text{com},k})\, dS \quad \text{for } k \in \{x, y, z\}$

The integration identity

For any two functions linear over a triangle (i.e. linearly interpolated from vertex values), barycentric integration gives a clean closed form:

$\iint_T f\,g\, dS = \frac{S}{12}\left[\left(\sum_{i=1}^{3} f_i\right)\!\left(\sum_{i=1}^{3} g_i\right) + \sum_{i=1}^{3} f_i\,g_i\right]$

This is a standard finite-element result, falls out of integrating products of linear shape functions over the reference triangle.

Applying it

Set $f_i = y_i - h$ and $g_i = x_{k,i} - x_{\text{com},k}$. Let $\text{sumY} \equiv y_1+y_2+y_3$ and $\text{sumX}_k \equiv x_{k,1}+x_{k,2}+x_{k,3}$. Then:

$\sum f_i = \text{sumY} - 3h, \qquad \sum g_i = \text{sumX}_k - 3\,x_{\text{com},k}$

$\sum f_i\,g_i = \sum_i x_{k,i}\,y_i \;-\; x_{\text{com},k}\,\text{sumY} \;-\; h\,\text{sumX}_k \;+\; 3h\,x_{\text{com},k}$

Plug both into the identity, expand $(\text{sumY} - 3h)(\text{sumX}_k - 3\,x_{\text{com},k})$, and collect:

$I_k = \frac{S}{12}\Big[\;\underbrace{\text{sumX}_k\,(\text{sumY} - 4h) + \sum_i x_{k,i}\,y_i}_{\text{call this } A_k} \;+\; (12h - 4\,\text{sumY})\,x_{\text{com},k}\;\Big]$

So defining $A_k = \text{sumX}_k(\text{sumY} - 4h) + \sum_i x_{k,i} y_i$ gives the per-triangle torque contribution:

$\boldsymbol{\tau}_T = \frac{\rho g S}{12}\Big[\,\mathbf{A} - 4(\text{sumY} - 3h)\,\mathbf{x}_{\text{com}}\,\Big] \times \mathbf{n}$

This matches the paper's Eq. 6 exactly. The disagreement is only in what $\mathbf{A}$ looks like.

Component by component

For $k = x$:

$A_x = \text{sumX}\,(\text{sumY} - 4h) + (x_1 y_1 + x_2 y_2 + x_3 y_3)$

For $k = z$, by symmetry:

$A_z = \text{sumZ}\,(\text{sumY} - 4h) + (z_1 y_1 + z_2 y_2 + z_3 y_3)$

For $k = y$:

$A_y = \text{sumY}\,(\text{sumY} - 4h) + (y_1^2 + y_2^2 + y_3^2)$

These are the canonical forms straight from the integration identity. Note all three have a single $\text{sumY}$ inside the bracket, no factor of 2 anywhere.

Where the "2" in $A_2$ legitimately comes from

The Y component can be rewritten using the algebraic identity:

$(y_1+y_2+y_3)^2 = (y_1^2+y_2^2+y_3^2) + 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

which rearranges to:

$y_1^2 + y_2^2 + y_3^2 = \text{sumY}^2 - 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

Substitute this into $A_y$:

$A_y = \text{sumY}\,(\text{sumY} - 4h) + \text{sumY}^2 - 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

Combine the two $\text{sumY}^2$ terms:

$A_y = \text{sumY}\,(2\,\text{sumY} - 4h) - 2(y_1 y_2 + y_2 y_3 + y_3 y_1)$

That's exactly the paper's Eq. 8. The "2" multiplying $\text{sumY}$ is real, it pays for the $-2 \cdot \Sigma\text{(pairwise)}$ correction that comes with it.

Why $A_x$ and $A_z$ can't do the same trick

The cross-product sums $\sum_i x_i y_i$ and $\sum_i z_i y_i$ have no perfect-square identity that would let you absorb an extra $\text{sumX} \cdot \text{sumY}$ (resp. $\text{sumZ} \cdot \text{sumY}$) into them. There simply isn't an algebraic equivalent of $(y_1+y_2+y_3)^2 = \Sigma y_i^2 + 2\Sigma_{i<j} y_i y_j$ that mixes two different coordinates.

So when the paper writes $A_1$ and $A_3$ with the same $\{-4h + 2\,\text{sumY}\}$ bracket as $A_2$, those 2s have nothing compensating them. Each one adds a spurious $\text{sumX}_k \cdot \text{sumY}$ term per triangle.

The most plausible explanation is that the authors derived $A_2$ in the rewritten form first (where $f$ and $g$ are the same variable, so the perfect-square identity applies naturally), then mechanically pattern-matched the same template to $A_1$ and $A_3$ during write-up. Easy mistake to make; nasty to catch unless you re-derive from the integral.

Numerical sanity check

Triangle $(0,0,0)$, $(1,0,0)$, $(0,1,0)$. Area $S = \tfrac{1}{2}$, $h = 0$, $\mathbf{x}_{\text{com}} = \mathbf{0}$. Direct integration:

$\iint_T x\,y\, dS = \tfrac{1}{24}$

Working backwards from $I_x = \tfrac{S}{12}\,A_x$, the correct value of $A_x$ must be:

$A_x = \frac{12}{S} \cdot \frac{1}{24} = 1$

Corrected formula: $\text{sumX} = 1$, $\text{sumY} = 1$, so

$A_x = 1 \cdot (1 - 0) + (0\cdot 0 + 1\cdot 0 + 0\cdot 1) = 1 \quad (Right)$

Paper's printed formula:

$A_x = 1 \cdot (2 \cdot 1 - 0) + 0 = 2 \quad (Wrong)$

Off by exactly the spurious $1 \cdot \text{sumX} \cdot \text{sumY} = 1$, as predicted. Same factor-of-2 error appears for $A_z$ on a corresponding triangle in the YZ plane.

Why it doesn't visually explode

The spurious term per triangle is proportional to $\text{sumX}_k \cdot \text{sumY} \cdot n_k$, summed over all triangles. For a perfectly axis-aligned cube centered at the origin and floating flat at $y = 0$, opposite-face triangles contribute equal-and-opposite errors and the bias mostly cancels. Tilt the cube or move it off origin and the cancellation breaks, leaving a constant residual torque every frame.

This is also why the bug is hard to spot from a quick smoke test, a stationary, level cube looks fine. It only becomes obvious when you validate at the analytical equilibrium tilt (like Hirae's own Tables 1 and 2 do at $26.565°$): the residual torque should be on the order of $10^{-7}$ but with the printed formula it sits orders of magnitude higher.

Empirical verification

Talk is cheap, so I implemented both formulas and ran Hirae's Section 5 static-equilibrium test on its own terms: 1m³ cube, mass 750 kg (density ratio 0.75), flat water at $y=0$, no waves, angular momentum damping $\alpha = 0.97$, rotation projected onto a single horizontal axis so the cube behaves like the 2D rectangular bar that Hirae's cited analytical reference (Igarashi & Nakamura 2007) actually solves. Settle for ~10 seconds, then log the residual force, torque, and tilt every 50 frames.

The 2D-bar setup matters here. An unconstrained 3D cube at $\rho_\text{ratio}=0.75$ tilts further and settles corner-down at $54.74° = \arccos(1/\sqrt{3})$, which is a perfectly valid equilibrium but not the one Hirae's tables report. Hirae's $26.565° = \arctan(1/2)$ is the 2D answer for a square cross-section, and to compare against it I have to take the cube's extra rotational degree of freedom away. With one horizontal axis pinned, the cross-section is exactly Igarashi & Nakamura's setup.

Three things to call out.

The corrected formula matches Hirae's own analytical reference to 0.01°. $26.57°$ vs $\arctan(1/2) = 26.5651°$. That's the exact same number Hirae reports in their Tables 1 and 2, reached on the exact same scenario, with residual torque at floating-point noise.

Force is unaffected, as predicted by the derivation. Both versions give $|F_b| = mg$ to the displayed precision. The typo is only in Eqs. 7 and 9 (the torque auxiliary vector); Eq. 5 (the force) doesn't depend on $A$ at all, so this is exactly what should happen.

Torque differs by ~5 orders of magnitude, and the buggy version finds a false equilibrium. Corrected: $|\tau| \approx 2 \times 10^{-4}$ N·m, implied lever arm $\sim 3 \times 10^{-8}$ m, i.e. nanometers on a 1m cube. Paper's printed formula: a rock-steady $31.5$ N·m bias, lever arm ~4 mm. Damping doesn't kill it, it just balances the bias against the damping rate, so the cube settles at a stable-looking pose 34° away from the analytical answer. This is the failure mode that lets a typo survive peer review and casual reimplementation: the simulation doesn't crash, the cube floats stably, the force balance is right. It's only when you check the angle against an analytical reference that the bug surfaces.