The bunny and the spheres fell into the water

Disclaimer: This method for “The bunny and the spheres fell into the water” is not my creative idea. I was inspired by the article Solving General Shallow Wave Equations on Surfaces by Huamin Wang, Gavin Miller and Greg Turk. And all the contents below comes from my understanding of this paper. The sphere is also shown in the original article, while I use the method for stanford bunny.

Bunny	Sphere	Both

The source cuda code can be found here, which you should noted is that the stb_image.h, stb_image_write.h and the stanford-bunny.obj must be placed at the same file level.

1. Top-View Shallow Water Method Overview

Primary target:

Reproduce SPHERES-style coupling behavior (waves/splashes from floating objects).
Keep identical mass but different volumes for spheres.
Replace sphere geometry with a Stanford Bunny option.
Render/export only top-view frames.

2. Connection to the Original GSWE Formulation

The paper writes (in continuous form):

u_{t} = - (u \cdot \nabla) u - \frac{1}{ρ} \nabla P_{e x t} + a_{e x t}

h_{t} + \nabla \cdot ((h - b) u) = 0

and derives a second-order height update with implicit force terms, finally solved in matrix form:

(A_{g} + A_{s} + I - I_{c}) h^{t} = b - b_{c}

In this implementation:

Surface is planar regular grid (not curved mesh GSWE).
Surface tension is disabled ( $γ = 0$ ), matching SPHERES table style.
Coupling with rigid bodies is kept through per-cell constraints (an $I_{c}, b_{c}$ -like term).

3. State Variables (Planar Grid + Rigid Bodies)

3.1 Fluid grid

On a 2D grid $(i, j)$ :

$h_{i, j}$ : water height
$u_{i, j}, v_{i, j}$ : horizontal velocity

3.2 Rigid bodies

Each body stores:

position $(x, y, z)$
velocity $(v_{x}, v_{y}, v_{z})$
mass $m$
equivalent radius $R_{e q}$
shape type: sphere or bunny

4. Time Integration

Each simulation step uses operator splitting:

Update rigid bodies and build splash source.
Build coupling maps from body-water contacts.
Advect height (semi-Lagrangian).
Build right-hand side with artificial viscosity term.
Solve implicit height system by Jacobi iterations.
Update velocity from height gradient + friction damping.
Apply boundary damping.

5. Height Advection and RHS

5.1 Semi-Lagrangian advection

For cell center $x_{ij}$ :

h_{a d v} (x_{ij}) = h^{n} (x_{ij} - Δ t u^{n} (x_{ij}))

5.2 Prospective RHS (paper-style $τ$ term + external source)

Code form:

b_{ij} = h_{a d v, ij} + (1 - τ) (h_{a d v, ij} - h_{ij}^{n - 1}) + Δ t S_{ij}

where $S$ is splash/body forcing written into the height equation.

In current settings: $τ = 0$ .

6. Implicit Height Solve with Coupling

The solver uses a local nonlinear Jacobi update:

σ_{ij} = \frac{Δ t ^{2} g}{Δ x ^{2}} max (h_{ij}, h_{min})

h_{ij}^{k + 1} = \frac{b _{ij} + c _{ij}^{w} c _{ij}^{t} + σ _{ij} ( h _{L} + h _{R} + h _{D} + h _{U} )}{1 + 4 σ _{ij} + c _{ij}^{w}}

$c^{w}$ : coupling weight map (constraint strength)
$c^{t}$ : coupling target (object bottom height proxy)

Interpretation:

$1 + 4 σ$ is the implicit gravity stiffness.
$c^{w}, c^{t}$ is the rigid/fluid coupling term (similar role to $I_{c}, b_{c}$ ).

7. Velocity Update and Friction

After solving $h^{n + 1}$ :

u \leftarrow u - Δ t g \partial_{x} h, v \leftarrow v - Δ t g \partial_{y} h

Then apply paper-style depth-dependent friction damping:

d_{f r i c} (h) = \frac{h}{h + d _{0}}, (u, v) \leftarrow d_{f r i c} (h) (u, v)

with constant $d_{0}$ in code.

8. Rigid-Body Dynamics Used for SPHERES/BUNNY

8.1 Buoyancy proxy

A sphere-cap proxy is used for submerged volume:

V_{c a p} (R, h) = π h^{2} (R - \frac{h}{3}), h \in [0, 2 R]

For bunny, the cap estimate is scaled by volume ratio:

V_{s u b}^{b u nn y} \approx V_{c a p} (R_{p ro x y}, \cdot) \cdot \frac{V _{b u nn y}}{V _{s p h ere} ( R _{p ro x y} )}

Buoyancy force:

F_{b} = ρ g V_{s u b}

8.2 Motion update

Vertical:

m \overset{z}{¨} = F_{b} - m g - c_{z} \overset{z}{˙}

Horizontal (wave-slope-driven drift):

\overset{x}{¨} = - k_{s} g \partial_{x} h - c_{x y} \overset{x}{˙}, \overset{y}{¨} = - k_{s} g \partial_{y} h - c_{x y} \overset{y}{˙}

plus floor and wall collision response with damping/restitution.

9. Coupling Map Construction

For each body, for each covered cell:

Compute an approximate object bottom elevation $z_{b o t} (x, y)$ .
Write:

c_{ij}^{w} = max (c_{ij}^{w}, k_{co u pl e}), c_{ij}^{t} = min (c_{ij}^{t}, z_{b o t})

For spheres:

z_{b o t} = z_{c} - R^{2} - r^{2}

For bunny:

Transform world $(x, y)$ into bunny local coordinates.
Sample precomputed profile $z_{b o t}^{l oc a l} (u, v)$ from OBJ-derived map.
Convert back to world height.

10. Splash Source Term

A local Gaussian source is added to height RHS:

S (x) + = A exp (- \frac{∥ x - x _{b} ∥ ^{2}}{2 σ ^{2}})

Amplitude $A$ is driven by:

downward impact speed
lateral speed
body size

This mimics the SPHERES visual effect (waves + splash rings) without introducing particles.