1b-agent/OneCAD-Simulation
1
0 Fork 0
forked from 1b-agent/fieldcraft-upstream
Code Pull Requests Activity

Lattice Boltzmann flow on the CPU: model, solver, engine (M12) (#36)

Browse Source 
Model: Newtonian fluid materials; fluid domains (box, spacing, conditions on
the box faces, obstacle bodies with conditions on their named faces, initial
velocity, body acceleration); fluid steps with caps. Deck writers refuse
fluid steps; a fluid-only model needs no mesh in model.check; one step
selector (Model::select_step) for every solver.

New L2 crate modelcraft-lbm: D3Q19 two-relaxation-time collision with Guo
forcing; halfway bounce-back walls (moving walls after Ladd), slip faces,
periodic pairs; Zou–He velocity and pressure faces with the D3Q19
transverse correction; obstacles voxelised from the face-tagged boundary
surface with exact predicates and symbolic perturbation; momentum-exchange
forces per box face and per B-rep face; unit conversion, automatic time step
and Mach/tau limits; bounded runs with frames and force histories.

Engine: model.fluid_domain.create/edit/delete, model.fluid_domain.boundary,
solve.fluid (obstacles from the open design or imported surfaces, lowerable
caps), solve.probe for fluid fields and wall force histories. Fluid Domain,
Inlet & Outlet and Lattice Boltzmann toolbar entries are ready.

Verification (analytic): plane Poiseuille exact with Λ = 3/16 and second
order (2.000) with Λ = 1/4; Couette exact with the wall shear; Taylor–Green
decay second order (2.01, 2.00); staircase pipe; Stokes drag of a sphere
array within 1.7% of Hasimoto at R = 6; Zou–He channel; mass conserved to
round-off; bit-for-bit reproducibility. Benchmark: Schäfer–Turek 2D-1.


Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm

Co-authored-by: jyviko <2501648+jyviko@users.noreply.github.com>
Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Irakli KourtisjyvikoClaude Opus 5.5
 authored and GitHub committed 2026-10-13 20:03:12 -04:00
1 parent bd55abea7f
commit a17752a49f
55 files changed
+3800 -64
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

crates/lbm/src/flow.rs
+598
View File
@@ -0,0 +1,598 @@
//! The lattice Boltzmann state and one time step: collide, stream, open boundaries.
//!
//! # Scheme
//! - **Lattice and equilibrium:** D3Q19 with the second-order equilibrium
//! (Qian, d'Humières & Lallemand 1992; Krüger et al. 2017, ch. 3).
//! - **Collision:** two relaxation times (Ginzburg, Verhaeghe & d'Humières 2008): the symmetric
//! part of each opposite pair relaxes with `1/τ⁺` (which sets the viscosity), the antisymmetric
//! part with `1/τ⁻`, related by the parameter `Λ = (τ⁺ − ½)(τ⁻ − ½)`. With `Λ = 3/16`
//! halfway bounce-back puts the wall exactly halfway between nodes and steady solutions do not
//! depend on the viscosity.
//! - **Body force:** Guo, Zheng & Shi 2002 (PRE 65:046308): the velocity is
//! `(Σ f c + F/2) / ρ` and the forcing term is split into its symmetric part, relaxed with
//! `1 − 1/(2τ⁺)`, and its antisymmetric part, with `1 − 1/(2τ⁻)` (the TRT form of Guo's
//! term, Krüger et al. 2017, ch. 10).
//! - **Walls:** halfway bounce-back at every link from a fluid cell to a solid cell or across a
//! wall face of the box; moving walls add Ladd's term `6 w_q ρ₀ c_q·u_w` (Ladd 1994, JFM
//! 271:285). Slip faces reflect populations specularly.
//! - **Inlets and outlets:** Zou & He 1997 (Phys. Fluids 9:1591) on the first layer of cells: the
//! unknown incoming populations are the opposite ones plus the equilibrium difference
//! (bounce-back of the non-equilibrium part), the density (velocity inlet) or the normal velocity
//! (pressure outlet) follow from mass and normal momentum, and the transverse momentum is
//! corrected on the diagonal unknowns as M. Hecht and J. Harting do for D3Q19 (*Implementation
//! of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations*, J. Stat.
//! Mech. 2010, P01018). The prescribed velocity is the one with the half force of Guo's scheme.
//! - **Forces:** momentum exchange on every wall link (Ladd 1994; R. Mei, D. Yu, W. Shyy, L.-S.
//! Luo, *Force evaluation in the lattice Boltzmann method involving curved geometry*, PRE
//! 65:041203, 2002): a population `f*` leaving a fluid cell toward the wall and the one `f'`
//! coming back transfer `c_out f* − c_in f'` to the wall.
//!
//! # Layout and determinism
//! Two arrays of populations, structure of arrays (`f[q n + cell]`), A–B pattern. Streaming
//! pulls: every population of a cell comes from one source (a neighbour, a periodic image, the
//! cell itself through a wall, or an open boundary). Cells whose 18 neighbours are all fluid
//! cells inside the box stream by one shifted copy per direction; the others follow per-link
//! rules. Every loop runs in a fixed order with no parallel reductions, so the result is the
//! same bit for bit on every platform.
use modelcraft_model::{BoxFace, FaceCondition, Profile};
use crate::lattice::{C, Q, W, c, equilibrium, moving_wall, opposite};
use crate::setup::{MAX_LATTICE_SPEED, Obstacle, Setup};
use crate::voxel::{Grid, TriangleTree, voxelise};
use crate::{Error, Result};
/// Where a population of a cell next to a wall or box face comes from.
#[derive(Debug, Clone, Copy, PartialEq)]
enum Rule {
/// The same direction from another cell (a neighbour or a periodic image).
Pull { from: u32 },
/// Bounce-back at a wall of force group `group`; `du` is the moving-wall term.
Bounce { group: u32, du: f64 },
/// Specular reflection at a slip face: direction `q` from cell `from`.
Reflect { from: u32, q: u8, group: u32 },
/// From outside an inlet or outlet: set by the open boundary after streaming.
Open,
}
/// A population, direction `q` of `cell`, and where it comes from.
#[derive(Debug, Clone, Copy)]
struct Link {
cell: u32,
q: u8,
rule: Rule,
}
/// An inlet or outlet cell.
#[derive(Debug, Clone, Copy)]
struct OpenCell {
cell: u32,
axis: usize,
/// +1 at a min face (the inward normal points up the axis), −1 at a max face.
inward: f64,
target: Target,
}
#[derive(Debug, Clone, Copy)]
enum Target {
/// The velocity (lattice units).
Velocity([f64; 3]),
/// The density.
Density(f64),
}
/// A force group: a box face or a face of an obstacle.
#[derive(Debug, Clone, PartialEq)]
pub struct Group {
/// `x_min` … for box faces; the face name for obstacle faces.
pub name: String,
/// The obstacle body, for obstacle faces.
pub body: Option<String>,
}
/// The flow state.
#[derive(Debug, Clone)]
pub struct Flow {
grid: Grid,
n: usize,
/// The distance between two directions' arrays in `f` and `g` (at least `n`).
stride: usize,
fluid: Vec<bool>,
f: Vec<f64>,
g: Vec<f64>,
/// The links that do not stream by the shifted copies, in cell order.
links: Vec<Link>,
open: Vec<OpenCell>,
omega_plus: f64,
omega_minus: f64,
acceleration: [f64; 3],
/// The momentum transferred to each group in the last step (lattice units).
forces: Vec<[f64; 3]>,
groups: Vec<Group>,
/// Per direction, the index offset of the neighbour a population streams from.
offset: [isize; Q],
fluid_cells: usize,
steps_done: u64,
/// Obstacles that cover no cell centre (smaller than the spacing, or outside the box).
empty_obstacles: Vec<String>,
}
const NONE: u32 = u32::MAX;
/// A vector allocated without aborting when memory runs out.
fn filled<T: Clone>(len: usize, value: T, what: &str) -> Result<Vec<T>> {
let mut v = Vec::new();
v.try_reserve_exact(len).map_err(|_| Error::Limit(format!("not enough memory for the {what} ({len} entries)")))?;
v.resize(len, value);
Ok(v)
}
impl Flow {
/// Build the grid (voxelising the obstacles), the streaming rules and the initial state (at
/// rest or the domain's initial velocity, density 1).
pub fn new(setup: &Setup, obstacles: &[Obstacle]) -> Result<Flow> {
let grid = setup.grid;
let [nx, ny, nz] = grid.cells;
let n = grid.len();
if n == 0 || n >= NONE as usize {
return Err(Error::Limit(format!("{n} cells cannot be addressed")));
}
// Each direction's array starts an odd number of cache lines (8 doubles) after the
// previous one, so the 19 arrays never fall into the same cache sets (a power-of-two grid
// would otherwise put them all in one).
let lines = n.div_ceil(8);
let stride = 8 * (lines + 1 - lines % 2);
let q_n = stride.checked_mul(Q).ok_or_else(|| Error::Limit("the grid is too large".into()))?;
// Solid cells, obstacle by obstacle (the first obstacle that covers a cell owns it).
let mut fluid = filled(n, true, "cell flags")?;
let mut owner = filled(n, NONE, "cell owners")?;
let mut trees = Vec::with_capacity(setup.obstacles.len());
let mut empty_obstacles = Vec::new();
// Force groups: the six box faces, then every face of every obstacle.
let mut groups: Vec<Group> = BoxFace::ALL.iter().map(|f| Group { name: f.name().to_string(), body: None }).collect();
let mut base = Vec::with_capacity(setup.obstacles.len());
for (k, faces) in setup.obstacles.iter().enumerate() {
let o = obstacles.get(faces.index).ok_or_else(|| Error::NotFound(format!("obstacle {}", faces.index)))?;
let mut inside = filled(n, false, "obstacle cells")?;
if voxelise(&o.surface, &grid, &mut inside)? == 0 {
empty_obstacles.push(o.body.clone());
}
for (c, &solid) in inside.iter().enumerate() {
if solid && owner[c] == NONE {
owner[c] = k as u32;
fluid[c] = false;
}
}
trees.push((TriangleTree::new(&o.surface), o.surface.tri_face()));
base.push(groups.len() as u32);
groups.extend(faces.names.iter().map(|name| Group { name: name.clone(), body: Some(o.body.clone()) }));
}
let domain = &setup.domain;
let periodic = [domain.faces.periodic(0), domain.faces.periodic(1), domain.faces.periodic(2)];
let wall_velocity = |face: BoxFace| match *domain.faces.get(face) {
FaceCondition::MovingWall { velocity } => Some(setup.to_lattice(velocity)),
FaceCondition::Wall {} => Some([0.0; 3]),
_ => None,
};
let is_slip = |face: BoxFace| matches!(domain.faces.get(face), FaceCondition::Slip {});
let dims = [nx as i64, ny as i64, nz as i64];
let index = |p: [i64; 3]| grid.index([p[0] as usize, p[1] as usize, p[2] as usize]) as u32;
let mut offset = [0isize; Q];
for (q, o) in offset.iter_mut().enumerate() {
let cq = C[q];
*o = cq[0] as isize + nx as isize * (cq[1] as isize + ny as isize * cq[2] as isize);
}
// The group and moving-wall term of a link into fluid cell `cell` (direction q) from the
// obstacle cell `solid`: the face of the triangle nearest the link's midpoint.
let obstacle_link = |cell: usize, solid: u32, q: usize| -> Rule {
let k = owner.get(solid as usize).copied().unwrap_or(NONE) as usize;
let (Some((tree, tags)), Some(faces), Some(&first)) = (trees.get(k), setup.obstacles.get(k), base.get(k)) else {
return Rule::Bounce { group: 0, du: 0.0 };
};
let (centre, cq) = (grid.centre(cell), c(q));
let mid = [0, 1, 2].map(|a| centre[a] - 0.5 * cq[a] * grid.spacing);
let tag = tree.nearest(mid).and_then(|t| tags.get(t as usize).copied()).unwrap_or(0);
let u = faces.velocities.get(tag as usize).copied().unwrap_or([0.0; 3]);
Rule::Bounce { group: first + tag, du: moving_wall(q, u) }
};
// The links that do not stream from the neighbour the shifted copies read.
let mut links: Vec<Link> = Vec::new();
let mut fluid_cells = 0;
for cell in 0..n {
if !fluid[cell] {
continue;
}
fluid_cells += 1;
let ijk = grid.coords(cell);
let p = [ijk[0] as i64, ijk[1] as i64, ijk[2] as i64];
for q in 1..Q {
let cq = C[q];
let mut s = [p[0] - i64::from(cq[0]), p[1] - i64::from(cq[1]), p[2] - i64::from(cq[2])];
let mut wrapped = false;
let mut crossed: [Option<BoxFace>; 3] = [None; 3];
for a in 0..3 {
if s[a] < 0 || s[a] >= dims[a] {
if periodic[a] {
s[a] = s[a].rem_euclid(dims[a]);
wrapped = true;
} else {
crossed[a] = Some(BoxFace::of(a, s[a] >= dims[a]));
}
}
}
let crossed_faces = || crossed.iter().flatten().copied();
let rule = match crossed_faces().next() {
None => {
let from = index(s);
match (fluid.get(from as usize).copied().unwrap_or(false), wrapped) {
(true, false) => continue,
(true, true) => Rule::Pull { from },
(false, _) => obstacle_link(cell, from, q),
}
}
Some(_) => {
if let Some(face) = crossed_faces().find(|f| wall_velocity(*f).is_some()) {
Rule::Bounce { group: face.index() as u32, du: moving_wall(q, wall_velocity(face).unwrap_or([0.0; 3])) }
} else if let Some(face) = crossed_faces().find(|f| is_slip(*f)) {
// The mirror source: the same cell layer, the tangential components shifted.
let a = face.axis();
let mut m = s;
m[a] = p[a];
let mut mirrored = C[q];
mirrored[a] = -mirrored[a];
let q2 = C.iter().position(|v| *v == mirrored).unwrap_or(opposite(q));
let single = crossed_faces().count() == 1;
let inside = (0..3).all(|b| m[b] >= 0 && m[b] < dims[b]);
match inside.then(|| index(m)) {
Some(from) if single && fluid.get(from as usize).copied().unwrap_or(false) => {
Rule::Reflect { from, q: q2 as u8, group: face.index() as u32 }
}
_ => Rule::Bounce { group: face.index() as u32, du: 0.0 },
}
} else {
Rule::Open
}
}
};
links.push(Link { cell: cell as u32, q: q as u8, rule });
}
}
// Inlets and outlets: the fluid cells of each open face's layer.
let mut open = Vec::new();
for face in BoxFace::ALL {
let a = face.axis();
let layer = if face.is_max() { grid.cells[a] - 1 } else { 0 };
let (t1, t2) = ((a + 1) % 3, (a + 2) % 3);
for j in 0..grid.cells[t2] {
for i in 0..grid.cells[t1] {
let mut ijk = [0usize; 3];
ijk[a] = layer;
ijk[t1] = i;
ijk[t2] = j;
let cell = grid.index(ijk);
if !fluid[cell] {
continue;
}
let target = match *domain.faces.get(face) {
FaceCondition::Velocity { velocity, profile } => {
let mut scale = 1.0;
if profile == Profile::Parabolic {
for t in [t1, t2].into_iter().filter(|t| !periodic[*t]) {
let s = (ijk[t] as f64 + 0.5) / grid.cells[t] as f64;
scale *= 4.0 * s * (1.0 - s);
}
}
Target::Velocity(setup.to_lattice(velocity).map(|v| v * scale))
}
FaceCondition::Pressure { pressure } => Target::Density(setup.lattice_density(pressure)),
_ => break,
};
open.push(OpenCell { cell: cell as u32, axis: a, inward: if face.is_max() { -1.0 } else { 1.0 }, target });
}
}
}
// The initial state: equilibrium at density 1 and the initial velocity.
let mut f = filled(q_n, 0.0, "populations")?;
let g = filled(q_n, 0.0, "populations")?;
for cell in 0..n {
for q in 0..Q {
f[q * stride + cell] = if fluid[cell] { equilibrium(q, 1.0, setup.initial_velocity) } else { 0.0 };
}
}
Ok(Flow {
grid,
n,
stride,
fluid,
f,
g,
links,
open,
omega_plus: 1.0 / setup.tau_plus,
omega_minus: 1.0 / setup.tau_minus,
acceleration: setup.acceleration,
forces: vec![[0.0; 3]; groups.len()],
groups,
offset,
fluid_cells,
steps_done: 0,
empty_obstacles,
})
}
pub fn grid(&self) -> &Grid {
&self.grid
}
pub fn fluid_cells(&self) -> usize {
self.fluid_cells
}
pub fn is_fluid(&self, cell: usize) -> bool {
self.fluid.get(cell).copied().unwrap_or(false)
}
pub fn groups(&self) -> &[Group] {
&self.groups
}
/// The momentum each group received in the last step (lattice units: force per step).
pub fn forces(&self) -> &[[f64; 3]] {
&self.forces
}
/// Obstacles that cover no cell centre.
pub fn empty_obstacles(&self) -> &[String] {
&self.empty_obstacles
}
pub fn steps_done(&self) -> u64 {
self.steps_done
}
/// Set every fluid cell to the equilibrium of `state(cell) = (density, velocity)` in lattice
/// units (velocity as defined with the half force).
pub fn set_equilibrium(&mut self, state: impl Fn(usize) -> (f64, [f64; 3])) {
let n = self.n;
let s = self.stride;
for cell in 0..n {
if !self.fluid[cell] {
continue;
}
let (rho, u) = state(cell);
let bare = [0, 1, 2].map(|a| u[a] - 0.5 * self.acceleration[a]);
for q in 0..Q {
self.f[q * s + cell] = equilibrium(q, rho, bare);
}
}
}
/// The density and velocity of a cell (lattice units), or `None` for a solid cell.
pub fn moments(&self, cell: usize) -> Option<(f64, [f64; 3])> {
if !self.is_fluid(cell) {
return None;
}
let s = self.stride;
let (mut rho, mut j) = (0.0, [0.0; 3]);
for q in 0..Q {
let fq = self.f.get(q * s + cell).copied().unwrap_or(0.0);
let cq = c(q);
rho += fq;
j[0] += fq * cq[0];
j[1] += fq * cq[1];
j[2] += fq * cq[2];
}
let u = [0, 1, 2].map(|a| (j[a] + 0.5 * rho * self.acceleration[a]) / rho);
Some((rho, u))
}
/// The total mass of the fluid cells (lattice units: the sum of every population).
pub fn mass(&self) -> f64 {
let n = self.n;
let s = self.stride;
let mut m = 0.0;
for cell in (0..n).filter(|c| self.fluid[*c]) {
for q in 0..Q {
m += self.f[q * s + cell];
}
}
m
}
/// One time step. Fails with [`Error::Unstable`] when a density stops being a positive
/// number or a lattice speed exceeds [`MAX_LATTICE_SPEED`].
pub fn step(&mut self) -> Result<()> {
let max_u2 = self.collide();
if !(max_u2 <= MAX_LATTICE_SPEED * MAX_LATTICE_SPEED) {
return Err(Error::Unstable(format!(
"the flow went unstable at time step {}: the lattice speed reached {:.3} (limit {MAX_LATTICE_SPEED}) or a density stopped being a positive number; lower the time step, refine the spacing or check the boundary conditions",
self.steps_done + 1,
max_u2.sqrt()
)));
}
self.stream();
self.open_boundaries();
std::mem::swap(&mut self.f, &mut self.g);
self.steps_done += 1;
Ok(())
}
/// Two-relaxation-time collision with Guo forcing, in place. Returns the largest squared
/// lattice speed, or NaN when a density is not positive.
fn collide(&mut self) -> f64 {
let n = self.n;
let s = self.stride;
let (wp, wm) = (self.omega_plus, self.omega_minus);
let (sp, sm) = (1.0 - 0.5 * wp, 1.0 - 0.5 * wm);
let a = self.acceleration;
let forced = a != [0.0; 3];
let mut max_u2: f64 = 0.0;
let mut fq = [0.0; Q];
for cell in 0..n {
if !self.fluid[cell] {
continue;
}
let (mut rho, mut j) = (0.0, [0.0; 3]);
for q in 0..Q {
let v = self.f[q * s + cell];
fq[q] = v;
let cq = C[q];
rho += v;
j[0] += v * f64::from(cq[0]);
j[1] += v * f64::from(cq[1]);
j[2] += v * f64::from(cq[2]);
}
if !(rho > 0.0 && rho.is_finite()) {
return f64::NAN;
}
let force = [rho * a[0], rho * a[1], rho * a[2]];
let u = [(j[0] + 0.5 * force[0]) / rho, (j[1] + 0.5 * force[1]) / rho, (j[2] + 0.5 * force[2]) / rho];
let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
max_u2 = max_u2.max(uu);
let uf = u[0] * force[0] + u[1] * force[1] + u[2] * force[2];
// Rest population: symmetric only.
let eq0 = W[0] * rho * (1.0 - 1.5 * uu);
let mut out0 = fq[0] - wp * (fq[0] - eq0);
if forced {
out0 += sp * W[0] * (-3.0 * uf);
}
self.f[cell] = out0;
for q in (1..Q).step_by(2) {
let o = q + 1;
let cq = C[q];
let cv = [f64::from(cq[0]), f64::from(cq[1]), f64::from(cq[2])];
let cu = cv[0] * u[0] + cv[1] * u[1] + cv[2] * u[2];
let w = W[q];
let eq_plus = w * rho * (1.0 + 4.5 * cu * cu - 1.5 * uu);
let eq_minus = w * rho * 3.0 * cu;
let plus = 0.5 * (fq[q] + fq[o]);
let minus = 0.5 * (fq[q] - fq[o]);
let relax_plus = wp * (plus - eq_plus);
let relax_minus = wm * (minus - eq_minus);
let (mut s_plus, mut s_minus) = (0.0, 0.0);
if forced {
let cf = cv[0] * force[0] + cv[1] * force[1] + cv[2] * force[2];
s_plus = sp * w * (9.0 * cu * cf - 3.0 * uf);
s_minus = sm * w * 3.0 * cf;
}
self.f[q * s + cell] = fq[q] - relax_plus - relax_minus + s_plus + s_minus;
self.f[o * s + cell] = fq[o] - relax_plus + relax_minus + s_plus - s_minus;
}
}
max_u2
}
/// Pull streaming into `g`, with the wall rules and the momentum exchange.
fn stream(&mut self) {
let n = self.n;
let s = self.stride;
// Every direction as one shifted copy: right for cells whose neighbours are all fluid
// cells in the box; the others are overwritten by their rules below.
for q in 0..Q {
let off = self.offset[q];
let (src, dst) = (&self.f[q * s..q * s + n], &mut self.g[q * s..q * s + n]);
if off >= 0 {
let o = off as usize;
if o < n {
dst[o..].copy_from_slice(&src[..n - o]);
}
} else {
let o = off.unsigned_abs();
if o < n {
dst[..n - o].copy_from_slice(&src[o..]);
}
}
}
self.forces.fill([0.0; 3]);
for link in &self.links {
let (cell, q) = (link.cell as usize, link.q as usize);
let cq = c(q);
match link.rule {
Rule::Pull { from } => self.g[q * s + cell] = self.f[q * s + from as usize],
Rule::Bounce { group, du } => {
let out = self.f[opposite(q) * s + cell];
let back = out + du;
self.g[q * s + cell] = back;
if let Some(force) = self.forces.get_mut(group as usize) {
// c_out f* − c_in f' with c_out = −c_q.
for a in 0..3 {
force[a] -= cq[a] * (out + back);
}
}
}
Rule::Reflect { from, q: q2, group } => {
let q2 = q2 as usize;
let out = self.f[q2 * s + from as usize];
self.g[q * s + cell] = out;
let c2 = c(q2);
if let Some(force) = self.forces.get_mut(group as usize) {
for a in 0..3 {
force[a] += c2[a] * out - cq[a] * out;
}
}
}
// A placeholder the open boundary replaces.
Rule::Open => self.g[q * s + cell] = self.f[opposite(q) * s + cell],
}
}
}
/// Zou–He velocity and pressure boundaries on the streamed populations.
fn open_boundaries(&mut self) {
let s = self.stride;
let a_half = self.acceleration.map(|v| 0.5 * v);
for oc in &self.open {
let cell = oc.cell as usize;
let axis = oc.axis;
let sign = oc.inward;
let ci = |q: usize| f64::from(C[q][axis]);
// Tangential (c_n = 0) and outgoing (c_n = −s, known) sums.
let (mut tangential, mut outgoing) = (0.0, 0.0);
for q in 0..Q {
let v = self.g[q * s + cell];
let cn = ci(q) * sign;
if cn == 0.0 {
tangential += v;
} else if cn < 0.0 {
outgoing += v;
}
}
// The bare velocity (without the half force) the populations must carry.
let (rho, bare) = match oc.target {
Target::Velocity(u) => {
let bare = [u[0] - a_half[0], u[1] - a_half[1], u[2] - a_half[2]];
let un = sign * bare[axis];
((tangential + 2.0 * outgoing) / (1.0 - un), bare)
}
Target::Density(rho) => {
let un = 1.0 - (tangential + 2.0 * outgoing) / rho;
let mut bare = [-a_half[0], -a_half[1], -a_half[2]];
bare[axis] = sign * un;
(rho, bare)
}
};
// Unknowns: bounce-back of the non-equilibrium part.
for q in (1..Q).filter(|q| ci(*q) * sign > 0.0) {
let cq = c(q);
let cu = cq[0] * bare[0] + cq[1] * bare[1] + cq[2] * bare[2];
self.g[q * s + cell] = self.g[opposite(q) * s + cell] + 6.0 * W[q] * rho * cu;
}
// Transverse momentum: correct the two diagonal unknowns along each tangential axis.
for t in (0..3).filter(|t| *t != axis) {
let mut m = -rho * bare[t];
for q in 0..Q {
m += self.g[q * s + cell] * f64::from(C[q][t]);
}
for q in (1..Q).filter(|q| ci(*q) * sign > 0.0 && C[*q][t] != 0) {
self.g[q * s + cell] -= 0.5 * m * f64::from(C[q][t]);
}
}
}
}
}
crates/lbm/src/lattice.rs
+141
View File
@@ -0,0 +1,141 @@
//! The D3Q19 velocity set (Qian, d'Humières & Lallemand 1992; Krüger et al. 2017, ch. 3).
//!
//! Directions come in opposite pairs: 1–2, 3–4, …, 17–18, so the opposite of an odd `q` is
//! `q + 1` and of an even `q` (not 0) is `q − 1`.
/// The number of velocities.
pub const Q: usize = 19;
/// The lattice velocities `c_q` in lattice units.
pub const C: [[i32; 3]; Q] = [
[0, 0, 0],
[1, 0, 0],
[-1, 0, 0],
[0, 1, 0],
[0, -1, 0],
[0, 0, 1],
[0, 0, -1],
[1, 1, 0],
[-1, -1, 0],
[1, -1, 0],
[-1, 1, 0],
[1, 0, 1],
[-1, 0, -1],
[1, 0, -1],
[-1, 0, 1],
[0, 1, 1],
[0, -1, -1],
[0, 1, -1],
[0, -1, 1],
];
/// The weights `w_q`: 1/3 at rest, 1/18 along the axes, 1/36 along the face diagonals.
pub const W: [f64; Q] = [
1.0 / 3.0,
1.0 / 18.0,
1.0 / 18.0,
1.0 / 18.0,
1.0 / 18.0,
1.0 / 18.0,
1.0 / 18.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
1.0 / 36.0,
];
/// The squared speed of sound in lattice units, `c_s² = 1/3`.
pub const CS2: f64 = 1.0 / 3.0;
/// The direction opposite `q`.
#[inline(always)]
pub const fn opposite(q: usize) -> usize {
if q == 0 {
0
} else if q % 2 == 1 {
q + 1
} else {
q - 1
}
}
/// `c_q` as floating-point numbers.
#[inline(always)]
pub fn c(q: usize) -> [f64; 3] {
let v = C.get(q).copied().unwrap_or([0; 3]);
[f64::from(v[0]), f64::from(v[1]), f64::from(v[2])]
}
/// The second-order equilibrium `w_q ρ (1 + 3 c·u + 9/2 (c·u)² − 3/2 u²)`.
#[inline(always)]
pub fn equilibrium(q: usize, rho: f64, u: [f64; 3]) -> f64 {
let c = c(q);
let cu = c[0] * u[0] + c[1] * u[1] + c[2] * u[2];
let uu = u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
W.get(q).copied().unwrap_or(0.0) * rho * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * uu)
}
/// Ladd's moving-wall term for a population bounced back into direction `q` off a wall moving
/// with `u` (lattice units, reference density 1): `2 w_q c_q·u / c_s²` (Ladd 1994).
#[inline(always)]
pub fn moving_wall(q: usize, u: [f64; 3]) -> f64 {
let c = c(q);
2.0 * W.get(q).copied().unwrap_or(0.0) * (c[0] * u[0] + c[1] * u[1] + c[2] * u[2]) / CS2
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn moments_of_the_velocity_set() {
// Σ w = 1, Σ w c = 0, Σ w c_a c_b = c_s² δ_ab, and opposite directions pair up.
let mut m0 = 0.0;
let mut m1 = [0.0; 3];
let mut m2 = [[0.0; 3]; 3];
for q in 0..Q {
let c = c(q);
m0 += W[q];
for a in 0..3 {
m1[a] += W[q] * c[a];
for b in 0..3 {
m2[a][b] += W[q] * c[a] * c[b];
}
}
let o = opposite(q);
assert_eq!(opposite(o), q);
assert_eq!([C[o][0], C[o][1], C[o][2]], [-C[q][0], -C[q][1], -C[q][2]]);
assert_eq!(W[o], W[q]);
}
assert!((m0 - 1.0f64).abs() < 1e-15);
for a in 0..3 {
assert!(m1[a].abs() < 1e-15);
for b in 0..3 {
let want = if a == b { CS2 } else { 0.0 };
assert!((m2[a][b] - want).abs() < 1e-15);
}
}
// The equilibrium carries the density and momentum it is built from.
let (rho, u) = (1.03, [0.02, -0.01, 0.005]);
let (mut r, mut j) = (0.0, [0.0; 3]);
for q in 0..Q {
let f = equilibrium(q, rho, u);
r += f;
for a in 0..3 {
j[a] += f * c(q)[a];
}
}
assert!((r - rho).abs() < 1e-15);
for a in 0..3 {
assert!((j[a] - rho * u[a]).abs() < 1e-15);
}
}
}
crates/lbm/src/lib.rs
+107
View File
@@ -0,0 +1,107 @@
//! ModelCraft's lattice Boltzmann flow solver (CPU).
//!
//! It solves a fluid step of the model ([`modelcraft_model::Procedure::Fluid`]) in its fluid
//! domain: weakly compressible, isothermal, Newtonian flow on a uniform grid of cubic cells that
//! fills the domain's box, less the cells inside the obstacle bodies.
//!
//! - [`voxel`]: which cells the obstacles take, decided with exact predicates from the
//! face-tagged boundary surface, and which face each wall link meets.
//! - [`setup`]: the step in lattice units, the time step, and the stability limits.
//! - [`flow`]: the D3Q19 two-relaxation-time scheme with Guo forcing, halfway bounce-back
//! walls (moving walls after Ladd), slip faces, Zou–He inlets and outlets, and the
//! momentum-exchange force on every wall.
//! - [`solve`]: the time loop, output frames and force histories in model units.
//!
//! # Sources
//! - T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. M. Viggen, *The Lattice
//! Boltzmann Method: Principles and Practice*, Springer 2017 (velocity sets, equilibria, units,
//! stream–collide, bounce-back, forcing, boundary conditions).
//! - Y. H. Qian, D. d'Humières, P. Lallemand, *Lattice BGK models for Navier–Stokes equation*,
//! EPL 17(6):479, 1992 (D3Q19).
//! - I. Ginzburg, F. Verhaeghe, D. d'Humières, *Two-relaxation-time lattice Boltzmann scheme:
//! about parametrization, velocity, pressure and mixed boundary conditions*, Commun. Comput.
//! Phys. 3(2):427–478, 2008 (TRT and the parameter Λ).
//! - Z. Guo, C. Zheng, B. Shi, *Discrete lattice effects on the forcing term in the lattice
//! Boltzmann method*, PRE 65:046308, 2002 (body force).
//! - A. J. C. Ladd, *Numerical simulations of particulate suspensions via a discretized Boltzmann
//! equation. Part 1*, JFM 271:285–309, 1994 (bounce-back, moving walls, momentum exchange).
//! - R. Mei, D. Yu, W. Shyy, L.-S. Luo, *Force evaluation in the lattice Boltzmann method
//! involving curved geometry*, PRE 65:041203, 2002 (momentum exchange).
//! - Q. Zou, X. He, *On pressure and velocity boundary conditions for the lattice Boltzmann BGK
//! model*, Phys. Fluids 9(6):1591–1598, 1997; M. Hecht, J. Harting, J. Stat. Mech. 2010,
//! P01018 (the D3Q19 transverse correction).
//! - P. Lallemand, L.-S. Luo, *Theory of the lattice Boltzmann method: dispersion, dissipation,
//! isotropy, Galilean invariance, and stability*, PRE 61:6546, 2000 (stability limits).
//! - H. Edelsbrunner, E. P. Mücke, *Simulation of Simplicity*, ACM TOG 9(1):66–104, 1990;
//! J. R. Shewchuk 1997 (exact predicates, through `modelcraft_mesh::predicates`).
//!
//! No code of any lattice Boltzmann package was read; everything is written from these papers and
//! the book.
//!
//! # Determinism
//! Double precision, fixed loop orders, no threads and no parallel reductions: the same model
//! gives the same result bit for bit on every platform, wasm included.
#![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)]
#![forbid(unsafe_code)]
// Population arrays are indexed by cell and direction counts fixed when they are allocated.
#![allow(clippy::needless_range_loop)]
pub mod flow;
pub mod lattice;
pub mod setup;
pub mod solve;
pub mod voxel;
pub use flow::{Flow, Group};
pub use setup::{Obstacle, Setup, fluid_step};
pub use solve::{FluidResult, FluidStats, ForceHistory, Frame, run, solve_fluid};
pub use voxel::Grid;
#[derive(Debug, Clone, PartialEq, thiserror::Error)]
pub enum Error {
#[error("invalid input: {0}")]
Invalid(String),
#[error("not found: {0}")]
NotFound(String),
#[error("limit reached: {0}")]
Limit(String),
/// The run diverged.
#[error("unstable: {0}")]
Unstable(String),
#[error(transparent)]
Model(#[from] modelcraft_model::Error),
#[error(transparent)]
Mesh(#[from] modelcraft_mesh::Error),
}
pub type Result<T> = std::result::Result<T, Error>;
/// Caps on a run's size and work; a run beyond any of them is refused before it starts.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct Limits {
pub max_cells: u64,
pub max_steps: u64,
/// Fluid cells times time steps.
pub max_cell_updates: u64,
/// The memory the output frames may take.
pub max_frame_bytes: u64,
}
impl Limits {
/// 16 million cells: about 5 GB of populations in double precision.
pub const MAX_CELLS: u64 = 16_000_000;
pub const MAX_STEPS: u64 = 10_000_000;
pub const MAX_CELL_UPDATES: u64 = 1_000_000_000_000;
pub const MAX_FRAME_BYTES: u64 = 2_000_000_000;
}
impl Default for Limits {
fn default() -> Self {
Limits {
max_cells: Limits::MAX_CELLS,
max_steps: Limits::MAX_STEPS,
max_cell_updates: Limits::MAX_CELL_UPDATES,
max_frame_bytes: Limits::MAX_FRAME_BYTES,
}
}
}
crates/lbm/src/setup.rs
+334
View File
@@ -0,0 +1,334 @@
//! From a fluid step of the model to lattice units: the grid, the time step, the relaxation
//! times, and the checks that keep the scheme in its stable and accurate range.
//!
//! # Units (Krüger et al. 2017, ch. 7)
//! Lattice units take the spacing `Δx`, the time step `Δt` and the fluid's density `ρ₀` as 1. A
//! velocity converts by `Δx/Δt`, a kinematic viscosity by `Δx²/Δt`, an acceleration by `Δx/Δt²`,
//! a pressure by `ρ₀ Δx²/Δt²` and a force by `ρ₀ Δx⁴/Δt²`. The lattice viscosity fixes the
//! symmetric relaxation time, `ν = c_s² (τ⁺ − ½)` with `c_s² = 1/3`.
//!
//! # Time step
//! Unless the step gives one, `Δt` is the largest that keeps the characteristic lattice velocity
//! at [`TARGET_LATTICE_VELOCITY`] (Mach 0.087: the weakly compressible error is O(Ma²)) and
//! `τ⁺` at most [`AUTO_TAU_MAX`] (beyond it the scheme only gets less accurate in time). The
//! characteristic velocity is the largest of the inlet and wall velocities, the initial velocity,
//! `√(2 Δp / ρ)` across pressure faces and, for a driving acceleration `a`, the smaller of the
//! viscous plane-channel peak `a H² / (8 ν)` (H the box's smallest extent) and the free-fall speed
//! `√(2 a L)` (L its largest).
//!
//! # Limits
//! A run is refused, naming the numbers and what to change, when the lattice Mach number of the
//! characteristic velocity exceeds [`MAX_MACH`] (the scheme is weakly compressible and loses
//! stability toward Ma ≈ 0.5), when `τ⁺ < `[`MIN_TAU`] (the grid cannot resolve this viscosity at
//! this velocity; linear stability analysis, Lallemand & Luo 2000, shows the stable velocity range
//! shrinks to nothing as τ → ½), or when the lattice acceleration exceeds [`MAX_LATTICE_ACCELERATION`].
//! While running, a lattice speed above [`MAX_LATTICE_SPEED`] or a non-finite or non-positive
//! density stops the run with [`crate::Error::Unstable`].
use modelcraft_mesh::surface::BoundarySurface;
use modelcraft_model::{BoxFace, FaceCondition, FluidDomain, Model, Procedure, SurfaceCondition};
use crate::voxel::Grid;
use crate::{Error, Limits, Result};
/// The characteristic lattice velocity the automatic time step aims at.
pub const TARGET_LATTICE_VELOCITY: f64 = 0.05;
/// The largest symmetric relaxation time the automatic time step allows.
pub const AUTO_TAU_MAX: f64 = 1.0;
/// The largest lattice Mach number `u Δt / (Δx c_s)` of the characteristic velocity.
pub const MAX_MACH: f64 = 0.3;
/// The smallest symmetric relaxation time τ⁺.
pub const MIN_TAU: f64 = 0.505;
/// The largest lattice acceleration `a Δt² / Δx`.
pub const MAX_LATTICE_ACCELERATION: f64 = 0.01;
/// A lattice speed above this during a run means it went unstable (Mach 0.87).
pub const MAX_LATTICE_SPEED: f64 = 0.5;
/// The default two-relaxation-time parameter Λ: bounce-back walls exactly halfway between cells
/// (Ginzburg, Verhaeghe & d'Humières 2008).
pub const DEFAULT_MAGIC: f64 = 3.0 / 16.0;
/// The surface of an obstacle body and the persistent name of each of its faces (by the
/// triangles' face tags); prepared by the caller from the geometry.
#[derive(Debug, Clone)]
pub struct Obstacle {
pub body: String,
pub surface: BoundarySurface,
/// Face names by tag; tags without a name are called `face0`, `face1`, ….
pub faces: Vec<String>,
}
/// An obstacle of the domain as the run sees it.
#[derive(Debug, Clone, PartialEq)]
pub struct ObstacleFaces {
/// The obstacle's index in the caller's obstacle list.
pub index: usize,
/// Every face tag's name.
pub names: Vec<String>,
/// Every face tag's wall velocity in lattice units (zero for walls at rest).
pub velocities: Vec<[f64; 3]>,
}
/// A fluid step in lattice units.
#[derive(Debug, Clone, PartialEq)]
pub struct Setup {
pub step: String,
pub domain: FluidDomain,
pub grid: Grid,
/// ρ₀ and the kinematic viscosity ν = μ / ρ₀, in model units.
pub density: f64,
pub viscosity: f64,
pub time_step: f64,
pub steps: u64,
pub frames: u32,
/// τ⁺ and τ⁻ (Λ = (τ⁺ − ½)(τ⁻ − ½)).
pub tau_plus: f64,
pub tau_minus: f64,
pub magic: f64,
/// The characteristic velocity in model and lattice units, and its lattice Mach number.
pub velocity: f64,
pub lattice_velocity: f64,
pub mach: f64,
/// The body acceleration and initial velocity in lattice units.
pub acceleration: [f64; 3],
pub initial_velocity: [f64; 3],
/// The domain's obstacles, in its order.
pub obstacles: Vec<ObstacleFaces>,
pub warnings: Vec<String>,
}
fn norm(v: [f64; 3]) -> f64 {
(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt()
}
/// The fluid step named `step` (or the only fluid step when `step` is `None`) and its domain.
pub fn fluid_step<'m>(model: &'m Model, step: Option<&str>) -> Result<(&'m modelcraft_model::Step, &'m FluidDomain)> {
let s = model.select_step(step, |p| p.fluid_domain().is_some(), "fluid")?;
let name = s.procedure.fluid_domain().unwrap_or_default();
let domain = model.fluid_domain(name).ok_or_else(|| Error::NotFound(format!("fluid domain `{name}`")))?;
Ok((s, domain))
}
impl Setup {
/// Convert a fluid step, checking it against the stability limits and `limits`.
pub fn new(model: &Model, step: Option<&str>, obstacles: &[Obstacle], limits: &Limits) -> Result<Setup> {
// The model checked the domain and its material when they were defined.
let (s, domain) = fluid_step(model, step)?;
let domain = domain.clone();
let Procedure::Fluid { end_time, steps, time_step, frames, magic, .. } = &s.procedure else {
return Err(Error::Invalid(format!("step `{}` is not a fluid step", s.name)));
};
let material = model.material(&domain.material);
let (Some(density), Some(fluid)) = (material.and_then(|m| m.density), material.and_then(|m| m.fluid)) else {
return Err(Error::Invalid(format!("material `{}` needs a density and fluid behaviour", domain.material)));
};
let viscosity = fluid.viscosity() / density;
if !(viscosity > 0.0 && viscosity.is_finite()) {
return Err(Error::Invalid(format!("the kinematic viscosity μ/ρ of material `{}` is {viscosity}", domain.material)));
}
// The grid.
let counts = domain.cells()?;
let mut cells = [0usize; 3];
let mut total: u64 = 1;
for a in 0..3 {
cells[a] = usize::try_from(counts[a]).map_err(|_| Error::Limit("the grid is too large".into()))?;
total = total.saturating_mul(counts[a]);
}
if total > limits.max_cells {
return Err(Error::Limit(format!(
"the domain has {total} cells ({} × {} × {}), above the cap of {}; coarsen the spacing or raise `max_cells`",
cells[0], cells[1], cells[2], limits.max_cells
)));
}
for face in BoxFace::ALL {
let open = matches!(domain.faces.get(face), FaceCondition::Velocity { .. } | FaceCondition::Pressure { .. });
if open && cells[face.axis()] < 2 {
return Err(Error::Invalid(format!("box face {} is an inlet or outlet, so the box needs at least 2 cells across it", face.name())));
}
}
let grid = Grid { origin: domain.min, spacing: domain.spacing, cells };
let h = domain.spacing;
// The characteristic velocity.
let mut u_char: f64 = norm(domain.initial_velocity);
let mut pressures = Vec::new();
for face in BoxFace::ALL {
match *domain.faces.get(face) {
FaceCondition::MovingWall { velocity } | FaceCondition::Velocity { velocity, .. } => u_char = u_char.max(norm(velocity)),
FaceCondition::Pressure { pressure } => pressures.push(pressure),
FaceCondition::Wall {} | FaceCondition::Slip {} | FaceCondition::Periodic {} => {}
}
}
for sb in &domain.surfaces {
if let SurfaceCondition::MovingWall { velocity } = sb.condition {
u_char = u_char.max(norm(velocity));
}
}
let lo = pressures.iter().copied().fold(f64::INFINITY, f64::min);
let hi = pressures.iter().copied().fold(f64::NEG_INFINITY, f64::max);
if hi > lo {
u_char = u_char.max((2.0 * (hi - lo) / density).sqrt());
}
let a = norm(domain.acceleration);
if a > 0.0 {
let extents: Vec<f64> = (0..3).map(|k| domain.max[k] - domain.min[k]).collect();
let small = extents.iter().copied().fold(f64::INFINITY, f64::min);
let large = extents.iter().copied().fold(0.0, f64::max);
u_char = u_char.max((a * small * small / (8.0 * viscosity)).min((2.0 * a * large).sqrt()));
}
if !u_char.is_finite() {
return Err(Error::Invalid("the flow's velocity scale is not a finite number".into()));
}
// The time step and the number of steps.
let diffusive = (AUTO_TAU_MAX - 0.5) * h * h / (3.0 * viscosity);
let auto = if u_char > 0.0 { (TARGET_LATTICE_VELOCITY * h / u_char).min(diffusive) } else { diffusive };
let (dt, n) = match (*time_step, *end_time, *steps) {
(Some(dt), Some(t), _) => {
let n = (t / dt).round();
if !(n >= 1.0) {
return Err(Error::Invalid(format!("the end time {t} is shorter than the time step {dt}")));
}
(dt, n)
}
(Some(dt), None, Some(n)) => (dt, n as f64),
(None, Some(t), _) => {
let n = (t / auto).ceil().max(1.0);
(t / n, n)
}
(None, None, Some(n)) => (auto, n as f64),
(_, None, None) => return Err(Error::Invalid(format!("step `{}` needs `end_time` or `steps`", s.name))),
};
if !(dt > 0.0 && dt.is_finite()) {
return Err(Error::Invalid(format!("the time step {dt} is not a positive number")));
}
if !(n <= limits.max_steps as f64) {
return Err(Error::Limit(format!(
"the step needs {n} time steps of {dt}, above the cap of {}; shorten it or raise `max_steps`",
limits.max_steps
)));
}
let n = n as u64;
let updates = (total as f64) * (n as f64);
if !(updates <= limits.max_cell_updates as f64) {
return Err(Error::Limit(format!(
"{total} cells for {n} time steps is {updates:e} cell updates, above the cap of {:e}; coarsen, shorten or raise `max_cell_updates`",
limits.max_cell_updates as f64
)));
}
// Lattice parameters and the stability limits.
let nu = viscosity * dt / (h * h);
let tau_plus = 3.0 * nu + 0.5;
let lattice_velocity = u_char * dt / h;
let mach = lattice_velocity / crate::lattice::CS2.sqrt();
if mach > MAX_MACH {
return Err(Error::Invalid(format!(
"the lattice Mach number is {mach:.3} (velocity {u_char} with time step {dt} and spacing {h}), above {MAX_MACH}: lower the time step"
)));
}
if !(tau_plus >= MIN_TAU) {
return Err(Error::Invalid(format!(
"the relaxation time τ = {tau_plus:.5} is below {MIN_TAU}: the grid cannot resolve viscosity {viscosity} at this time step (cell Reynolds number {:.1}); refine the spacing or raise the time step",
if nu > 0.0 { lattice_velocity / nu } else { f64::INFINITY }
)));
}
let magic = magic.unwrap_or(DEFAULT_MAGIC);
let tau_minus = 0.5 + magic / (tau_plus - 0.5);
let acceleration = domain.acceleration.map(|c| c * dt * dt / h);
if norm(acceleration) > MAX_LATTICE_ACCELERATION {
return Err(Error::Invalid(format!(
"the lattice acceleration is {:.3e}, above {MAX_LATTICE_ACCELERATION}: lower the time step",
norm(acceleration)
)));
}
// Obstacles, their face names and their faces' wall velocities.
let mut domain_obstacles = Vec::with_capacity(domain.obstacles.len());
let mut warnings = Vec::new();
for body in &domain.obstacles {
let (index, o) = obstacles
.iter()
.enumerate()
.find(|(_, o)| &o.body == body)
.ok_or_else(|| Error::NotFound(format!("obstacle body `{body}` has no surface (open the design or import its surface)")))?;
let tags = o.surface.tri_face().iter().copied().max().map_or(1, |m| m as usize + 1).max(o.faces.len());
let names: Vec<String> = (0..tags).map(|t| o.faces.get(t).cloned().unwrap_or_else(|| format!("face{t}"))).collect();
let mut velocities = vec![[0.0; 3]; tags];
for sb in domain.surfaces.iter().filter(|sb| &sb.face.body == body) {
let mut hit = false;
for (tag, name) in names.iter().enumerate() {
if sb.face.selects(&modelcraft_model::GeometryRef::face(body.as_str(), name.as_str())?) {
hit = true;
if let (Some(slot), SurfaceCondition::MovingWall { velocity }) = (velocities.get_mut(tag), sb.condition) {
*slot = velocity.map(|c| c * dt / h);
}
}
}
if !hit {
return Err(Error::NotFound(format!("the {} is not a face of the obstacle's surface", sb.face)));
}
}
domain_obstacles.push(ObstacleFaces { index, names, velocities });
}
if mach > 0.1 + 1e-12 {
warnings.push(format!(
"lattice Mach number {mach:.3}: the compressibility error, of order Ma² ≈ {:.1}%, may be visible",
100.0 * mach * mach
));
}
if tau_plus < 0.51 {
warnings.push(format!("τ = {tau_plus:.4} is close to ½: under-resolved flows may go unstable; a finer spacing is safer"));
}
Ok(Setup {
step: s.name.clone(),
grid,
density,
viscosity,
time_step: dt,
steps: n,
frames: *frames,
tau_plus,
tau_minus,
magic,
velocity: u_char,
lattice_velocity,
mach,
acceleration,
initial_velocity: domain.initial_velocity.map(|c| c * dt / h),
obstacles: domain_obstacles,
warnings,
domain,
})
}
/// Model units per lattice unit of velocity, pressure and force.
pub fn velocity_scale(&self) -> f64 {
self.grid.spacing / self.time_step
}
pub fn pressure_scale(&self) -> f64 {
self.density * self.grid.spacing * self.grid.spacing / (self.time_step * self.time_step)
}
pub fn force_scale(&self) -> f64 {
self.pressure_scale() * self.grid.spacing * self.grid.spacing
}
/// A velocity in lattice units.
pub fn to_lattice(&self, v: [f64; 3]) -> [f64; 3] {
let k = 1.0 / self.velocity_scale();
v.map(|c| c * k)
}
/// The lattice density of a gauge pressure: `1 + p / (c_s² ρ₀ Δx²/Δt²)`.
pub fn lattice_density(&self, pressure: f64) -> f64 {
1.0 + pressure / (crate::lattice::CS2 * self.pressure_scale())
}
/// The gauge pressure of a lattice density (the inverse of [`Setup::lattice_density`]).
pub fn gauge_pressure(&self, rho: f64) -> f64 {
(rho - 1.0) * crate::lattice::CS2 * self.pressure_scale()
}
}
crates/lbm/src/solve.rs
+190
View File
@@ -0,0 +1,190 @@
//! Running a fluid step: the time loop, output frames and force histories, in model units.
use modelcraft_model::Model;
use crate::flow::Flow;
use crate::setup::{Obstacle, Setup};
use crate::voxel::Grid;
use crate::{Error, Limits, Result};
/// Most samples a force history keeps; longer runs record every k-th step.
pub const MAX_HISTORY: u64 = 10_000;
/// A relative density variation above this at the end of a run is reported.
pub const COMPRESSIBILITY_WARNING: f64 = 0.01;
/// The velocity and pressure field at one time, per cell (index `i + nx (j + ny k)`; zero in
/// solid cells). Stored in single precision: output for display and probing.
#[derive(Debug, Clone, PartialEq)]
pub struct Frame {
pub time: f64,
pub step: u64,
pub velocity: Vec<[f32; 3]>,
/// Gauge pressure.
pub pressure: Vec<f32>,
/// The total fluid mass (density times volume).
pub mass: f64,
}
/// The force the fluid exerts on one wall over time.
#[derive(Debug, Clone, PartialEq)]
pub struct ForceHistory {
/// `x_min` … for box faces, the face's persistent name for obstacle faces.
pub name: String,
/// The obstacle body, for obstacle faces.
pub body: Option<String>,
/// One value per [`FluidResult::history_times`].
pub force: Vec<[f64; 3]>,
}
/// Run statistics.
#[derive(Debug, Clone, PartialEq)]
pub struct FluidStats {
pub cells: usize,
pub fluid_cells: usize,
pub steps: u64,
pub time_step: f64,
pub end_time: f64,
pub tau_plus: f64,
pub tau_minus: f64,
pub magic: f64,
/// The characteristic velocity (model units), its lattice value and its lattice Mach number.
pub velocity: f64,
pub lattice_velocity: f64,
pub mach: f64,
/// Fluid cells times steps.
pub cell_updates: u64,
}
/// The result of a fluid step.
#[derive(Debug, Clone, PartialEq)]
pub struct FluidResult {
pub step: String,
pub grid: Grid,
/// Per cell: not fluid (inside an obstacle).
pub solid: Vec<bool>,
/// The initial field and `frames` more, evenly spaced to the end.
pub frames: Vec<Frame>,
pub history_times: Vec<f64>,
/// Every box wall and every obstacle face.
pub forces: Vec<ForceHistory>,
pub stats: FluidStats,
pub warnings: Vec<String>,
}
impl FluidResult {
/// The cell containing a point and the frame's velocity and pressure there (`None` outside
/// the grid; zeros in a solid cell).
pub fn probe(&self, frame: usize, p: [f64; 3]) -> Option<(usize, [f32; 3], f32)> {
let cell = self.grid.cell_at(p)?;
let f = self.frames.get(frame)?;
Some((cell, *f.velocity.get(cell)?, *f.pressure.get(cell)?))
}
}
/// Solve the fluid step named `step` (or the model's only fluid step) with the obstacles' surfaces
/// prepared by the caller.
pub fn solve_fluid(model: &Model, step: Option<&str>, obstacles: &[Obstacle], limits: &Limits) -> Result<FluidResult> {
let setup = Setup::new(model, step, obstacles, limits)?;
let mut flow = Flow::new(&setup, obstacles)?;
run(&setup, &mut flow, limits)
}
fn frame(setup: &Setup, flow: &Flow, step: u64) -> Result<Frame> {
let n = setup.grid.len();
let mut velocity = Vec::new();
let mut pressure = Vec::new();
velocity.try_reserve_exact(n).map_err(|_| Error::Limit("not enough memory for an output frame".into()))?;
pressure.try_reserve_exact(n).map_err(|_| Error::Limit("not enough memory for an output frame".into()))?;
let vs = setup.velocity_scale();
for cell in 0..n {
match flow.moments(cell) {
Some((rho, u)) => {
velocity.push(u.map(|c| (c * vs) as f32));
pressure.push(setup.gauge_pressure(rho) as f32);
}
None => {
velocity.push([0.0; 3]);
pressure.push(0.0);
}
}
}
let mass = flow.mass() * setup.density * setup.grid.spacing.powi(3);
Ok(Frame { time: step as f64 * setup.time_step, step, velocity, pressure, mass })
}
/// Run a prepared flow (fresh from [`Flow::new`], possibly with its own initial state) through
/// the setup's steps, recording frames and force histories.
pub fn run(setup: &Setup, flow: &mut Flow, limits: &Limits) -> Result<FluidResult> {
let n = setup.grid.len();
let total = setup.steps;
let frames = u64::from(setup.frames).min(total).max(1);
let frame_bytes = (frames + 1) as f64 * n as f64 * 16.0;
if !(frame_bytes <= limits.max_frame_bytes as f64) {
return Err(Error::Limit(format!(
"{} frames of {n} cells need {:.0} MB, above the cap of {:.0} MB; write fewer frames or raise `max_frame_bytes`",
frames + 1,
frame_bytes / 1e6,
limits.max_frame_bytes as f64 / 1e6
)));
}
let every = total.div_ceil(MAX_HISTORY).max(1);
let fs = setup.force_scale();
let mut out_frames = vec![frame(setup, flow, 0)?];
let mut history_times = Vec::new();
let mut forces: Vec<Vec<[f64; 3]>> = vec![Vec::new(); flow.groups().len()];
let mut next_frame = 1u64;
for k in 1..=total {
flow.step()?;
if k % every == 0 || k == total {
history_times.push(k as f64 * setup.time_step);
for (h, f) in forces.iter_mut().zip(flow.forces()) {
h.push(f.map(|c| c * fs));
}
}
// Frame j (1..=frames) at step round(j · total / frames).
while next_frame <= frames && k == (next_frame * total + frames / 2) / frames {
out_frames.push(frame(setup, flow, k)?);
next_frame += 1;
}
}
let mut warnings = setup.warnings.clone();
// The weakly compressible scheme turns pressure differences into density differences of
// Δp / (ρ c_s²) in lattice units; large ones show the time step is too long for the flow.
let spread = (0..n).filter_map(|c| flow.moments(c)).map(|(rho, _)| (rho - 1.0).abs()).fold(0.0, f64::max);
if spread > COMPRESSIBILITY_WARNING {
warnings.push(format!(
"the density varies by {:.1}% at the end: the weakly compressible error is visible; a shorter time step reduces it (quadratically)",
100.0 * spread
));
}
for body in flow.empty_obstacles() {
warnings.push(format!("obstacle `{body}` covers no cell centre: it is smaller than the spacing or outside the box"));
}
let solid: Vec<bool> = (0..n).map(|c| !flow.is_fluid(c)).collect();
Ok(FluidResult {
step: setup.step.clone(),
grid: setup.grid,
solid,
frames: out_frames,
history_times,
forces: flow.groups().iter().zip(forces).map(|(g, force)| ForceHistory { name: g.name.clone(), body: g.body.clone(), force }).collect(),
stats: FluidStats {
cells: n,
fluid_cells: flow.fluid_cells(),
steps: total,
time_step: setup.time_step,
end_time: total as f64 * setup.time_step,
tau_plus: setup.tau_plus,
tau_minus: setup.tau_minus,
magic: setup.magic,
velocity: setup.velocity,
lattice_velocity: setup.lattice_velocity,
mach: setup.mach,
cell_updates: (flow.fluid_cells() as u64).saturating_mul(total),
},
warnings,
})
}
crates/lbm/src/voxel.rs
+453
View File
@@ -0,0 +1,453 @@
//! Voxelisation: which cell centres of a grid lie inside a closed, outward-oriented triangle
//! surface, decided with exact predicates, and which triangle a lattice link meets.
//!
//! # Inside or outside
//! A cell centre `p` is inside when the winding number of the surface around it is positive. The
//! winding number is counted along the ray from `p` in +x: every triangle the ray crosses adds
//! the sign of its normal's x component (+1 leaving the solid, −1 entering it). All centres on one
//! grid line (fixed y, z) share the triangles their ray can meet, so each triangle is visited once
//! per line it covers:
//! - whether the line meets the triangle is the sign of three [`orient2d`] tests of the line's
//! (y, z) against the triangle's projected edges;
//! - which centres lie before the triangle's plane is the sign of [`orient3d`] of the triangle and
//! the centre, searched from a floating-point estimate of the crossing (the sign is monotone
//! along the line, so the estimate only decides where the exact search starts).
//!
//! Degenerate positions (a line through a projected edge or vertex, a centre on a triangle's
//! plane) are resolved by one symbolic perturbation of the centre, `p + (ε, δ, δ²)` with
//! infinitesimals `1 ≫ ε ≫ δ` (Edelsbrunner & Mücke, *Simulation of Simplicity*, ACM TOG 9(1),
//! 1990): the line moves by `(δ, δ²)` in (y, z), so a line through a shared edge or vertex meets
//! exactly one of the triangles around it, and a centre on a plane the line crosses moves off it
//! by `ε` in +x. Every test sees the same perturbed point, so a watertight surface gives a
//! consistent count and the result never depends on rounding. Triangles parallel to x are met by
//! no perturbed line and are skipped.
//!
//! # Link faces
//! Momentum-exchange forces are reported per B-rep face, so every wall link (a fluid centre next
//! to a solid one) is given the face of the triangle nearest its midpoint ([`TriangleTree`]): a
//! bounding-volume tree with branch-and-bound nearest search and the closest point on a triangle
//! after C. Ericson, *Real-Time Collision Detection*, Morgan Kaufmann 2005, §5.1.5. Ties go to the
//! lower triangle index, so the assignment is deterministic. This is a classification for
//! reporting, not a meshing decision, so floating-point distances are adequate.
use modelcraft_mesh::predicates::{Point2, Point3, Sign, orient2d, orient3d, snap_to_domain};
use modelcraft_mesh::surface::BoundarySurface;
use crate::{Error, Result};
/// A uniform grid of cubic cells: cell `(i, j, k)` has its centre at
/// `origin + (i + ½, j + ½, k + ½) · spacing` and index `i + nx (j + ny k)`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Grid {
pub origin: [f64; 3],
pub spacing: f64,
pub cells: [usize; 3],
}
impl Grid {
/// The number of cells.
pub fn len(&self) -> usize {
self.cells[0] * self.cells[1] * self.cells[2]
}
pub fn is_empty(&self) -> bool {
self.len() == 0
}
/// The index of cell `(i, j, k)`.
#[inline(always)]
pub fn index(&self, [i, j, k]: [usize; 3]) -> usize {
i + self.cells[0] * (j + self.cells[1] * k)
}
/// The `(i, j, k)` of a cell index.
#[inline(always)]
pub fn coords(&self, c: usize) -> [usize; 3] {
let [nx, ny, _] = self.cells;
[c % nx, (c / nx) % ny, c / (nx * ny)]
}
/// The centre coordinate of cell `i` along `axis`, snapped into the exact predicates' domain.
#[inline]
pub fn centre_along(&self, axis: usize, i: usize) -> f64 {
snap_to_domain(self.origin[axis] + (i as f64 + 0.5) * self.spacing)
}
/// The centre of a cell.
pub fn centre(&self, c: usize) -> [f64; 3] {
let ijk = self.coords(c);
[self.centre_along(0, ijk[0]), self.centre_along(1, ijk[1]), self.centre_along(2, ijk[2])]
}
/// The cell containing a point, if it is in the grid.
pub fn cell_at(&self, p: [f64; 3]) -> Option<usize> {
let mut ijk = [0usize; 3];
for a in 0..3 {
let t = (p[a] - self.origin[a]) / self.spacing;
if !(t >= 0.0 && t < self.cells[a] as f64) {
return None;
}
ijk[a] = (t.floor() as usize).min(self.cells[a] - 1);
}
Some(self.index(ijk))
}
}
/// The sign of `orient2d(a, b, p)` for `p` perturbed by `(δ, δ²)` in its two coordinates.
/// `orient2d` is affine in `p` with gradient `(a_v − b_v, b_u − a_u)`, so when it is zero the
/// sign is that of the first nonzero gradient component; both are zero only when `a = b`.
#[inline]
fn edge_sign(a: &Point2, b: &Point2, p: &Point2) -> Sign {
let s = orient2d(a, b, p);
if !s.is_zero() {
return s;
}
let ([au, av], [bu, bv]) = (a.coords(), b.coords());
if av != bv { Sign::of(av - bv) } else { Sign::of(bu - au) }
}
fn point3(surface: &BoundarySurface, v: u32) -> Result<Point3> {
surface.point3(v).ok_or_else(|| Error::Invalid(format!("surface vertex {v} is missing")))
}
/// Mark every cell of `grid` whose centre is inside `surface` (a positive winding number):
/// `inside[c]` is set to true for those. Returns how many cells were marked that were not marked
/// before.
pub fn voxelise(surface: &BoundarySurface, grid: &Grid, inside: &mut [bool]) -> Result<usize> {
let [nx, ny, nz] = grid.cells;
if inside.len() != grid.len() {
return Err(Error::Invalid(format!("{} flags for {} cells", inside.len(), grid.len())));
}
let lines = ny * nz;
// Winding-number changes along each line: +s at the line's first cell and −s at the first
// centre beyond the triangle's plane.
let mut steps = vec![0i32; lines * (nx + 1)];
let h = grid.spacing;
let range = |axis: usize, lo: f64, hi: f64, n: usize| -> (usize, usize) {
// Centres with lo ≤ centre ≤ hi, widened by one cell on each side: the exact tests decide.
let first = ((lo - grid.origin[axis]) / h - 0.5).floor() - 1.0;
let last = ((hi - grid.origin[axis]) / h - 0.5).ceil() + 1.0;
let clamp = |t: f64| if t.is_nan() || t < 0.0 { 0usize } else { (t as usize).min(n - 1) };
(clamp(first), clamp(last))
};
for tri in surface.triangles() {
let [a, b, c] = [point3(surface, tri[0])?, point3(surface, tri[1])?, point3(surface, tri[2])?];
let proj = |p: &Point3| -> Result<Point2> {
let [_, y, z] = p.coords();
Ok(Point2::new(y, z)?)
};
let (a2, b2, c2) = (proj(&a)?, proj(&b)?, proj(&c)?);
let s = orient2d(&a2, &b2, &c2);
if s.is_zero() {
continue;
}
let ([ax, ay, az], [bx, by, bz], [cx, cy, cz]) = (a.coords(), b.coords(), c.coords());
let (j0, j1) = range(1, ay.min(by).min(cy), ay.max(by).max(cy), ny);
let (k0, k1) = range(2, az.min(bz).min(cz), az.max(bz).max(cz), nz);
// The plane's normal, for the floating-point estimate of the crossing only.
let (u, v) = ([bx - ax, by - ay, bz - az], [cx - ax, cy - ay, cz - az]);
let n = [u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]];
for k in k0..=k1 {
let z = grid.centre_along(2, k);
for j in j0..=j1 {
let y = grid.centre_along(1, j);
let p2 = Point2::new(y, z)?;
if edge_sign(&a2, &b2, &p2) != s || edge_sign(&b2, &c2, &p2) != s || edge_sign(&c2, &a2, &p2) != s {
continue;
}
// Centres before the plane (the +x ray from them crosses the triangle):
// orient3d(a, b, c, p) = −(p − a)·n has the sign s there. A centre on the plane is
// perturbed by +ε in x (n_x has the sign s) and counts as beyond it.
let behind = |i: usize| -> Result<bool> { Ok(orient3d(&a, &b, &c, &Point3::new(grid.centre_along(0, i), y, z)?) == s) };
let x = if n[0] != 0.0 { ax - (n[1] * (y - ay) + n[2] * (z - az)) / n[0] } else { ax };
let t = ((x - grid.origin[0]) / h - 0.5).ceil();
let mut split = if t.is_nan() || t < 0.0 { 0 } else { (t as usize).min(nx) };
while split > 0 && !behind(split - 1)? {
split -= 1;
}
while split < nx && behind(split)? {
split += 1;
}
let line = (j + ny * k) * (nx + 1);
let w = s as i32;
if let Some(first) = steps.get_mut(line) {
*first += w;
}
if let Some(end) = steps.get_mut(line + split) {
*end -= w;
}
}
}
}
let mut marked = 0;
for line in 0..lines {
let mut w = 0i32;
for i in 0..nx {
w += steps.get(line * (nx + 1) + i).copied().unwrap_or(0);
if w > 0
&& let Some(flag) = inside.get_mut(line * nx + i)
&& !*flag
{
*flag = true;
marked += 1;
}
}
}
Ok(marked)
}
/// A bounding-volume tree over a surface's triangles for nearest-triangle queries.
#[derive(Debug, Clone)]
pub struct TriangleTree {
triangles: Vec<[[f64; 3]; 3]>,
order: Vec<u32>,
nodes: Vec<Node>,
}
#[derive(Debug, Clone, Copy)]
struct Node {
lo: [f64; 3],
hi: [f64; 3],
/// Leaf: the range `first..first + count` of `order`; inner: `count == 0` and the children
/// are `first` and `first + 1`.
first: u32,
count: u32,
}
const LEAF: usize = 4;
fn bounds(tris: &[[[f64; 3]; 3]], ids: &[u32]) -> ([f64; 3], [f64; 3]) {
let mut lo = [f64::INFINITY; 3];
let mut hi = [f64::NEG_INFINITY; 3];
for &t in ids {
if let Some(tri) = tris.get(t as usize) {
for p in tri {
for a in 0..3 {
lo[a] = lo[a].min(p[a]);
hi[a] = hi[a].max(p[a]);
}
}
}
}
(lo, hi)
}
fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
}
fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}
/// The squared distance from `p` to triangle `abc` (Ericson 2005, §5.1.5: the closest point by
/// Voronoi regions of the vertices, edges and face).
fn distance2(p: [f64; 3], [a, b, c]: &[[f64; 3]; 3]) -> f64 {
let (ab, ac, ap) = (sub(*b, *a), sub(*c, *a), sub(p, *a));
let at = |s: f64, t: f64| [a[0] + s * ab[0] + t * ac[0], a[1] + s * ab[1] + t * ac[1], a[2] + s * ab[2] + t * ac[2]];
let d2 = |q: [f64; 3]| {
let d = sub(p, q);
dot(d, d)
};
let (d1, d2v) = (dot(ab, ap), dot(ac, ap));
if d1 <= 0.0 && d2v <= 0.0 {
return d2(*a);
}
let bp = sub(p, *b);
let (d3, d4) = (dot(ab, bp), dot(ac, bp));
if d3 >= 0.0 && d4 <= d3 {
return d2(*b);
}
let vc = d1 * d4 - d3 * d2v;
if vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0 {
return d2(at(d1 / (d1 - d3), 0.0));
}
let cp = sub(p, *c);
let (d5, d6) = (dot(ab, cp), dot(ac, cp));
if d6 >= 0.0 && d5 <= d6 {
return d2(*c);
}
let vb = d5 * d2v - d1 * d6;
if vb <= 0.0 && d2v >= 0.0 && d6 <= 0.0 {
return d2(at(0.0, d2v / (d2v - d6)));
}
let va = d3 * d6 - d5 * d4;
if va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0 {
let w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
return d2([b[0] + w * (c[0] - b[0]), b[1] + w * (c[1] - b[1]), b[2] + w * (c[2] - b[2])]);
}
let denom = 1.0 / (va + vb + vc);
d2(at(vb * denom, vc * denom))
}
fn box_distance2(p: [f64; 3], n: &Node) -> f64 {
let mut d = 0.0;
for a in 0..3 {
let e = (n.lo[a] - p[a]).max(0.0).max(p[a] - n.hi[a]);
d += e * e;
}
d
}
impl TriangleTree {
pub fn new(surface: &BoundarySurface) -> TriangleTree {
let pts = surface.points();
let at = |v: u32| pts.get(v as usize).copied().unwrap_or([0.0; 3]);
let triangles: Vec<[[f64; 3]; 3]> = surface.triangles().iter().map(|t| [at(t[0]), at(t[1]), at(t[2])]).collect();
let centroid = |t: u32, a: usize| triangles.get(t as usize).map_or(0.0, |tri| tri[0][a] + tri[1][a] + tri[2][a]);
let mut order: Vec<u32> = (0..triangles.len() as u32).collect();
let mut nodes = Vec::with_capacity(2 * triangles.len() / LEAF + 1);
let (lo, hi) = bounds(&triangles, &order);
nodes.push(Node { lo, hi, first: 0, count: order.len() as u32 });
// Split nodes breadth-first; each split halves a range, so this ends after fewer than
// 2 n / LEAF nodes.
let mut next = 0;
while next < nodes.len() {
let Some(node) = nodes.get(next).copied() else { break };
next += 1;
let (first, count) = (node.first as usize, node.count as usize);
if count <= LEAF {
continue;
}
let Some(ids) = order.get_mut(first..first + count) else { continue };
let extent = sub(node.hi, node.lo);
let axis = if extent[0] >= extent[1] && extent[0] >= extent[2] {
0
} else if extent[1] >= extent[2] {
1
} else {
2
};
ids.sort_by(|&p, &q| centroid(p, axis).total_cmp(&centroid(q, axis)).then(p.cmp(&q)));
let half = count / 2;
let (left, right) = ids.split_at(half);
let (llo, lhi) = bounds(&triangles, left);
let (rlo, rhi) = bounds(&triangles, right);
let child = nodes.len() as u32;
nodes.push(Node { lo: llo, hi: lhi, first: first as u32, count: half as u32 });
nodes.push(Node { lo: rlo, hi: rhi, first: (first + half) as u32, count: (count - half) as u32 });
if let Some(n) = nodes.get_mut(next - 1) {
n.first = child;
n.count = 0;
}
}
TriangleTree { triangles, order, nodes }
}
/// The triangle nearest `p` (the lowest index among equally near ones).
pub fn nearest(&self, p: [f64; 3]) -> Option<u32> {
let mut best: Option<(f64, u32)> = None;
let better = |d: f64, t: u32, best: &Option<(f64, u32)>| best.is_none_or(|(bd, bt)| d < bd || (d == bd && t < bt));
let mut stack = vec![0u32];
// Every node is pushed at most once.
while let Some(i) = stack.pop() {
let Some(node) = self.nodes.get(i as usize) else { continue };
if best.is_some_and(|(bd, _)| box_distance2(p, node) > bd) {
continue;
}
if node.count > 0 {
let range = node.first as usize..(node.first + node.count) as usize;
for &t in self.order.get(range).unwrap_or_default() {
if let Some(tri) = self.triangles.get(t as usize) {
let d = distance2(p, tri);
if better(d, t, &best) {
best = Some((d, t));
}
}
}
} else {
let (l, r) = (node.first, node.first + 1);
let dist = |c: u32| self.nodes.get(c as usize).map_or(f64::INFINITY, |n| box_distance2(p, n));
// Visit the nearer child first (pushed last).
if dist(l) <= dist(r) {
stack.extend([r, l]);
} else {
stack.extend([l, r]);
}
}
}
best.map(|(_, t)| t)
}
}
#[cfg(test)]
mod tests {
use super::*;
use modelcraft_mesh::surface::SurfaceOptions;
/// The surface of the box [lo, hi], outward, two triangles per side, side k tagged k.
pub(crate) fn cube(lo: [f64; 3], hi: [f64; 3]) -> BoundarySurface {
let p = |i: usize| [if i & 1 == 0 { lo[0] } else { hi[0] }, if i & 2 == 0 { lo[1] } else { hi[1] }, if i & 4 == 0 { lo[2] } else { hi[2] }];
let positions: Vec<[f64; 3]> = (0..8).map(p).collect();
// Quads counter-clockwise seen from outside.
let quads = [[0, 4, 6, 2], [1, 3, 7, 5], [0, 1, 5, 4], [2, 6, 7, 3], [0, 2, 3, 1], [4, 5, 7, 6]];
let mut triangles = Vec::new();
let mut tags = Vec::new();
for (k, q) in quads.iter().enumerate() {
triangles.push([q[0], q[1], q[2]]);
triangles.push([q[0], q[2], q[3]]);
tags.extend([k as u32; 2]);
}
BoundarySurface::new(&positions, &triangles, &tags, &SurfaceOptions::default()).unwrap()
}
#[test]
fn projected_orientation_is_the_normal_x_sign() {
let pt = |x, y, z| Point2::new(y, z).map(|p| (x, p)).unwrap().1;
// Triangle (0,0,0), (0,1,0), (0,0,1): normal +x; its (y, z) projection is counter-clockwise.
assert_eq!(orient2d(&pt(0.0, 0.0, 0.0), &pt(0.0, 1.0, 0.0), &pt(0.0, 0.0, 1.0)), Sign::Positive);
}
#[test]
fn a_box_is_voxelised_exactly_even_with_faces_on_centres() {
// Cells of size 1 from 0 to 8; the box's faces pass through cell centres (x = 2.5, 5.5)
// and between them (y = 2, 6), and its edges lie on grid lines.
let grid = Grid { origin: [0.0; 3], spacing: 1.0, cells: [8, 8, 8] };
let s = cube([2.5, 2.0, 2.5], [5.5, 6.0, 5.5]);
let mut inside = vec![false; grid.len()];
let n = voxelise(&s, &grid, &mut inside).unwrap();
// Centres on a face are decided by the perturbation p + (ε, δ, δ²): beyond it in +x, y and z.
let mut count = 0;
for c in 0..grid.len() {
let [x, y, z] = grid.centre(c);
let strictly = x > 2.5 && x < 5.5 && y > 2.0 && y < 6.0 && z > 2.5 && z < 5.5;
let maybe = (2.5..=5.5).contains(&x) && (2.0..=6.0).contains(&y) && (2.5..=5.5).contains(&z);
if strictly {
assert!(inside[c], "{x} {y} {z}");
}
if !maybe {
assert!(!inside[c], "{x} {y} {z}");
}
count += usize::from(inside[c]);
}
assert_eq!(n, count);
// x: 2.5 (moved inside), 3.5, 4.5 (5.5 moves outside); y: 2.5 … 5.5; z as x.
assert_eq!(n, 3 * 4 * 3, "{n}");
}
#[test]
fn nearest_triangle_and_its_face() {
let s = cube([0.0; 3], [1.0; 3]);
let tree = TriangleTree::new(&s);
// Just outside the +y side, near the -x edge: side 3 (y = hi).
let t = tree.nearest([0.1, 1.2, 0.5]).unwrap();
assert_eq!(s.tri_face()[t as usize], 3);
let t = tree.nearest([0.5, 0.5, -0.3]).unwrap();
assert_eq!(s.tri_face()[t as usize], 4);
// Against brute force on scattered points.
for i in 0..200 {
let f = |k: u32| ((i * 7919 + k * 104_729) % 1000) as f64 / 400.0 - 0.75;
let p = [f(1), f(2), f(3)];
let brute = (0..s.triangles().len() as u32)
.min_by(|&a, &b| {
let tri = |t: u32| {
let v = s.triangles()[t as usize];
[s.points()[v[0] as usize], s.points()[v[1] as usize], s.points()[v[2] as usize]]
};
distance2(p, &tri(a)).total_cmp(&distance2(p, &tri(b))).then(a.cmp(&b))
})
.unwrap();
assert_eq!(tree.nearest(p), Some(brute), "{p:?}");
}
}
}