physics.c file

Functions

void vm_euler_semi(vm_float_t* x, vm_float_t* v, const vm_float_t* a, const int n, const vm_float_t dt): Semi-implicit Euler: v += a dt, then x += v dt.
void vm_verlet(const vm_acc_fn acc, vm_float_t* x, vm_float_t* v, vm_float_t* a, const int n, const vm_float_t dt, void* ctx): Velocity Verlet with an acceleration callback.
void vm_rk2(const vm_ode_fn f, vm_float_t* y, const int n, const vm_float_t dt, void* ctx): Explicit midpoint RK2 for y' = f(y).
void vm_rk4(const vm_ode_fn f, vm_float_t* y, const int n, const vm_float_t dt, void* ctx): Classic RK4 for y' = f(y).
auto vm_cfl_dt(const vm_float_t cfl, const vm_float_t dx, const vm_float_t speed) -> vm_float_t: CFL timestep dt = cfl * dx / (|u| + ε).
auto mat3_chol(const matrix3 a, matrix3* L) -> bool: 3×3 Cholesky A = L Lᵀ.
static auto mat3_chol_solve(const matrix3 L, const vector3 b, vector3* x) -> bool: Solve L Lᵀ x = b given a lower-triangular Cholesky factor L.
auto mat3_spd_solve(const matrix3 a, const vector3 b, vector3* x) -> bool: Solve the 3×3 SPD system A x = b via Cholesky.
auto vm_inertia_world(const matrix3 ib, const quaternion q) -> matrix3: World-frame inertia I_w = R I_b Rᵀ from a body tensor and orientation.
auto vm_omega_from_angmom(const matrix3 I, const vector3 L) -> vector3: Recover ω from angular momentum L = I ω.
auto vm_rigid_energy(const vm_float_t mass, const vector3 v, const matrix3 I, const vector3 w) -> vm_float_t: Rigid kinetic energy ½ m |v|² + ½ ω · (I ω).
void vm_rigid_step(vector3* x, vector3* v, quaternion* q, vector3* w, const vector3 F, const vector3 tau, const vm_float_t mass, const matrix3 I_body, const vm_float_t dt): One symplectic-Euler rigid step.
void vm_baumgarte_correct(vector3* x, vector3* v, const vector3 n, const vm_float_t C, const vm_float_t beta, const vm_float_t gamma, const vm_float_t dt): Single-constraint Baumgarte correction along a unit normal.
auto vm_grid3_make(const int nx, const int ny, const int nz, const vm_float_t dx, const vm_float_t dy, const vm_float_t dz) -> vm_grid3: Build a 3-D grid descriptor.
auto vm_grid_ncells(const vm_grid3 g) -> int: Returns the total number of cells in the grid.
auto vm_grid_cell(const vm_grid3 g, const int i, const int j, const int k) -> int: Compute linear index of a cell in a 3-D grid stored in row-major order.
auto vm_mac_nu(const vm_grid3 g) -> int: Returns the number of u-velocity MAC face values for the grid.
auto vm_mac_nv(const vm_grid3 g) -> int: Returns the number of MAC grid v-velocity components.
auto vm_mac_nw(const vm_grid3 g) -> int: Returns the number of MAC grid faces in the z (vertical) direction.
auto vm_mac_u(const vm_grid3 g, const int i, const int j, const int k) -> int: Compute linear index of u-velocity on MAC grid.
auto vm_mac_v(const vm_grid3 g, const int i, const int j, const int k) -> int: Computes the linear index for the v-component of a MAC grid velocity.
auto vm_mac_w(const vm_grid3 g, const int i, const int j, const int k) -> int: Compute linear index of MAC grid w-component at cell (i,j,k).
static auto vm_grid_on_boundary(const vm_grid3 g, const int i, const int j, const int k) -> bool: Test if a grid cell is on the domain boundary.
auto vm_grid_laplacian(vm_spmat* A, const vm_grid3 g, const vm_bc_t bc) -> bool: Assemble the SPD operator -∇² on a cell-centered grid.
void vm_mac_div(vm_float_t* div, const vm_float_t* u, const vm_float_t* v, const vm_float_t* w, const vm_grid3 g): Cell-centered divergence of a MAC velocity field.
void vm_mac_grad(vm_float_t* gu, vm_float_t* gv, vm_float_t* gw, const vm_float_t* p, const vm_grid3 g): MAC face gradient of a cell-centered scalar (pressure).
void vm_mac_curl_z(vm_float_t* cz, const vm_float_t* u, const vm_float_t* v, const vm_grid3 g): Cell-centered z-vorticity (∂v/∂x − ∂u/∂y) from MAC u, v.

Function documentation

void vm_euler_semi(vm_float_t* x, vm_float_t* v, const vm_float_t* a, const int n, const vm_float_t dt)

Semi-implicit Euler: v += a dt, then x += v dt.

Parameters
x	Position vector (in/out); length `n`.
v	Velocity vector (in/out); length `n`.
a	Acceleration vector; length `n`.
n	State dimension.
dt	Timestep.

void vm_verlet(const vm_acc_fn acc, vm_float_t* x, vm_float_t* v, vm_float_t* a, const int n, const vm_float_t dt, void* ctx)

Velocity Verlet with an acceleration callback.

Parameters
acc	Acceleration callback `acc(x, a, ctx)`.
x	Position vector (in/out); length `n`.
v	Velocity vector (in/out); length `n`.
a	Acceleration vector (in/out); length `n`.
n	State dimension.
dt	Timestep.
ctx	User context passed to `acc`.

Uses the incoming a at x, advances x, re-evaluates acc, then completes the velocity half-kick.

void vm_rk2(const vm_ode_fn f, vm_float_t* y, const int n, const vm_float_t dt, void* ctx)

Explicit midpoint RK2 for y' = f(y).

Parameters
f	ODE right-hand side `f(y, dy, ctx)`.
y	State vector (in/out); length `n`.
n	State dimension.
dt	Timestep.
ctx	User context passed to `f`.

void vm_rk4(const vm_ode_fn f, vm_float_t* y, const int n, const vm_float_t dt, void* ctx)

Classic RK4 for y' = f(y).

Parameters
f	ODE right-hand side `f(y, dy, ctx)`.
y	State vector (in/out); length `n`.
n	State dimension.
dt	Timestep.
ctx	User context passed to `f`.

vm_float_t vm_cfl_dt(const vm_float_t cfl, const vm_float_t dx, const vm_float_t speed)

CFL timestep dt = cfl * dx / (|u| + ε).

Parameters
cfl	CFL number (typically in `(0, 1]`).
dx	Characteristic cell size.
speed	Characteristic speed (e.g. `\|u\|`).
Returns	Stable timestep estimate.

bool mat3_chol(const matrix3 a, matrix3* L)

3×3 Cholesky A = L Lᵀ.

Parameters
a	SPD coefficient matrix.
L	Lower-triangular Cholesky factor (out); upper triangle set to 0.
Returns	`true` on success, `false` if `L` is NULL or `a` is not SPD.

L is lower; the upper triangle is zeroed.

static bool mat3_chol_solve(const matrix3 L, const vector3 b, vector3* x)

Solve L Lᵀ x = b given a lower-triangular Cholesky factor L.

Parameters
L	Lower-triangular Cholesky factor of an SPD 3×3 matrix.
b	Right-hand side vector.
x	Solution vector (out).
Returns	`true` on success.

Performs forward substitution for L y = b, then back substitution for Lᵀ x = y.

bool mat3_spd_solve(const matrix3 a, const vector3 b, vector3* x)

Solve the 3×3 SPD system A x = b via Cholesky.

Parameters
a	SPD coefficient matrix.
b	Right-hand side vector.
x	Solution vector (out).
Returns	`true` on success, `false` if `x` is NULL or factorization fails.

matrix3 vm_inertia_world(const matrix3 ib, const quaternion q)

World-frame inertia I_w = R I_b Rᵀ from a body tensor and orientation.

Parameters
ib	Body-frame inertia tensor.
q	Orientation quaternion.
Returns	World-frame inertia tensor.

vector3 vm_omega_from_angmom(const matrix3 I, const vector3 L)

Recover ω from angular momentum L = I ω.

Parameters
I	Inertia tensor (same frame as `L`).
L	Angular momentum.
Returns	Angular velocity, or the zero vector if the solve fails.

vm_float_t vm_rigid_energy(const vm_float_t mass, const vector3 v, const matrix3 I, const vector3 w)

Rigid kinetic energy ½ m |v|² + ½ ω · (I ω).

Parameters
mass	Mass.
v	Linear velocity.
I	Inertia tensor (same frame as `w`).
w	Angular velocity (same frame as `I`).
Returns	Kinetic energy.

I and ω must share a frame.

void vm_rigid_step(vector3* x, vector3* v, quaternion* q, vector3* w, const vector3 F, const vector3 tau, const vm_float_t mass, const matrix3 I_body, const vm_float_t dt)

One symplectic-Euler rigid step.

Parameters
x	World-frame position (in/out).
v	World-frame linear velocity (in/out).
q	Orientation quaternion (in/out).
w	Body-frame angular velocity (in/out).
F	World-frame force.
tau	Body-frame torque.
mass	Mass.
I_body	Body-frame inertia tensor (SPD 3×3).
dt	Timestep.

x, v, F are world-frame. w and tau are body-frame. I_body is the body inertia (any SPD 3×3). Orientation is advanced with quat_integrate.

void vm_baumgarte_correct(vector3* x, vector3* v, const vector3 n, const vm_float_t C, const vm_float_t beta, const vm_float_t gamma, const vm_float_t dt)

Single-constraint Baumgarte correction along a unit normal.

Parameters
x	Position to correct (in/out).
v	Velocity to correct (in/out).
n	Unit constraint normal.
C	Signed constraint value (`0` at contact).
beta	Position Baumgarte coefficient.
gamma	Velocity Baumgarte coefficient.
dt	Timestep used for the velocity correction scale.

C is the signed constraint value (0 at contact). Position is moved by -beta C n; velocity by -gamma C / dt n.

vm_grid3 vm_grid3_make(const int nx, const int ny, const int nz, const vm_float_t dx, const vm_float_t dy, const vm_float_t dz)

Build a 3-D grid descriptor.

Parameters
nx	Number of cells in x.
ny	Number of cells in y.
nz	Number of cells in z; values `<= 0` are treated as `1`.
dx	Cell spacing in x.
dy	Cell spacing in y.
dz	Cell spacing in z.
Returns	Grid descriptor with the given dimensions and spacing.

int vm_grid_ncells(const vm_grid3 g)

Returns the total number of cells in the grid.

Parameters
g	3D grid descriptor.
Returns	Total cell count (`nx * ny * nz`), or `0` if `nx` or `ny` is non-positive.

int vm_grid_cell(const vm_grid3 g, const int i, const int j, const int k)

Compute linear index of a cell in a 3-D grid stored in row-major order.

Parameters
g	grid dimensions and spacing
i	cell index along x
j	cell index along y
k	cell index along z
Returns	flattened 1-D index of the cell

int vm_mac_nu(const vm_grid3 g)

Returns the number of u-velocity MAC face values for the grid.

Parameters
g	Grid definition.
Returns	Number of u-faces, or zero if the grid is invalid.

int vm_mac_nv(const vm_grid3 g)

Returns the number of MAC grid v-velocity components.

Parameters
g	Grid dimensions and spacing.
Returns	Number of v-velocity samples.

int vm_mac_nw(const vm_grid3 g)

Returns the number of MAC grid faces in the z (vertical) direction.

Parameters
g	Grid dimensions and spacing.
Returns	Number of vertical MAC faces, or 0 if the grid is invalid.

int vm_mac_u(const vm_grid3 g, const int i, const int j, const int k)

Compute linear index of u-velocity on MAC grid.

Parameters
g	Grid descriptor.
i	Cell index in x.
j	Cell index in y.
k	Cell index in z.
Returns	Linear array index for the u-face value.

int vm_mac_v(const vm_grid3 g, const int i, const int j, const int k)

Computes the linear index for the v-component of a MAC grid velocity.

Parameters
g	Grid descriptor.
i	Cell index in x.
j	Cell index in y.
k	Cell index in z.
Returns	Linear array index for v at the given staggered location.

int vm_mac_w(const vm_grid3 g, const int i, const int j, const int k)

Compute linear index of MAC grid w-component at cell (i,j,k).

Parameters
g	Grid descriptor.
i	Cell index in x.
j	Cell index in y.
k	Cell index in z.
Returns	Linear array index for the w velocity component.

static bool vm_grid_on_boundary(const vm_grid3 g, const int i, const int j, const int k)

Test if a grid cell is on the domain boundary.

Parameters
g	3-D grid descriptor.
i	cell index in x direction.
j	cell index in y direction.
k	cell index in z direction.
Returns	true if the cell lies on any face of the grid, false otherwise.

bool vm_grid_laplacian(vm_spmat* A, const vm_grid3 g, const vm_bc_t bc)

Assemble the SPD operator -∇² on a cell-centered grid.

Parameters
A	Output sparse matrix; size `vm_grid_ncells(g)`.
g	Grid dimensions and spacing.
bc	Boundary condition (`VM_BC_DIRICHLET` or `VM_BC_NEUMANN`).
Returns	`true` on success, `false` on invalid input or allocation failure.

5-point in 2-D (nz == 1), 7-point in 3-D. Dirichlet boundary cells become identity rows. Homogeneous Neumann drops the missing neighbour (singular constant nullspace).

void vm_mac_div(vm_float_t* div, const vm_float_t* u, const vm_float_t* v, const vm_float_t* w, const vm_grid3 g)

Cell-centered divergence of a MAC velocity field.

Parameters
div	Output cell-centered divergence; length `vm_grid_ncells(g)`.
u	MAC face-centered x-velocity; length `vm_mac_nu(g)`.
v	MAC face-centered y-velocity; length `vm_mac_nv(g)`.
w	MAC face-centered z-velocity, or NULL when `g.nz <= 1`; length `vm_mac_nw(g)`.
g	Grid dimensions and spacing.

w may be NULL when g.nz <= 1. div has length vm_grid_ncells(g).

void vm_mac_grad(vm_float_t* gu, vm_float_t* gv, vm_float_t* gw, const vm_float_t* p, const vm_grid3 g)

MAC face gradient of a cell-centered scalar (pressure).

Parameters
gu	Output MAC face gradient in x; length `vm_mac_nu(g)`.
gv	Output MAC face gradient in y; length `vm_mac_nv(g)`.
gw	Output MAC face gradient in z, or NULL when `g.nz <= 1`; length `vm_mac_nw(g)`.
p	Cell-centered scalar field; length `vm_grid_ncells(g)`.
g	Grid dimensions and spacing.

Boundary faces are left at 0. gw may be NULL when g.nz <= 1.

void vm_mac_curl_z(vm_float_t* cz, const vm_float_t* u, const vm_float_t* v, const vm_grid3 g)

Cell-centered z-vorticity (∂v/∂x − ∂u/∂y) from MAC u, v.

Parameters
cz	Output cell-centered vorticity; length `vm_grid_ncells(g)`.
u	MAC face-centered x-velocity.
v	MAC face-centered y-velocity.
g	Grid dimensions and spacing.