diff --git a/crates/render/src/contour.rs b/crates/render/src/contour.rs index af1a5ac..e46fcfa 100644 --- a/crates/render/src/contour.rs +++ b/crates/render/src/contour.rs @@ -5,7 +5,7 @@ //! sharp iso-lines on any renderer, with no shader support: the GPU viewport and the CPU //! rasterizer draw the same flat-coloured pieces. -use crate::{Rgb, Vec3}; +use crate::{Rgb, SceneTri, Vec3}; /// Equal bands from `min` to `max`. #[derive(Clone, Copy, Debug, PartialEq)] @@ -94,6 +94,19 @@ pub fn split_triangle(p: [Vec3; 3], v: [f64; 3], bands: &Bands, out: &mut impl F } } +/// [`split_triangle`] with each piece fanned into scene triangles of `normal`, coloured by `color` +/// of its band (`None`: a value that is not a number). +pub fn banded(p: [Vec3; 3], v: [f64; 3], normal: Vec3, bands: &Bands, color: &impl Fn(Option) -> Rgb, out: &mut Vec) { + split_triangle(p, v, bands, &mut |piece, band| { + let color = color(band); + for k in 1..piece.len().saturating_sub(1) { + if let (Some(&a), Some(&b), Some(&c)) = (piece.first(), piece.get(k), piece.get(k + 1)) { + out.push(SceneTri { p: [a, b, c], normal, color }); + } + } + }); +} + /// The part of a convex polygon where the value is at least (`keep_above`) or at most `level` /// (Sutherland–Hodgman against one level). fn clip(poly: &[(Vec3, f64)], level: f64, keep_above: bool) -> Vec<(Vec3, f64)> { @@ -189,6 +202,12 @@ mod tests { assert_eq!(k, None); }); assert_eq!(n, 2); + // Fanned into scene triangles, keeping the area and the bands. + let mut tris = Vec::new(); + banded(p, [0.0, 10.0, 0.0], Vec3::Z, &b, &|k| Rgb(k.map_or(0, |k| k as u8), 0, 0), &mut tris); + assert!((tris.iter().map(|t| area(&t.p)).sum::() - 50.0).abs() < 1e-9); + assert_eq!(tris.iter().map(|t| t.color.0).max(), Some(4)); + assert!(tris.iter().all(|t| t.normal == Vec3::Z)); } /// A clamped range puts the values above it in the "over" band, split exactly at the diff --git a/crates/render/src/flow.rs b/crates/render/src/flow.rs new file mode 100644 index 0000000..2da5802 --- /dev/null +++ b/crates/render/src/flow.rs @@ -0,0 +1,615 @@ +//! Display data of a flow field on a regular grid of cells (a lattice Boltzmann result): a slice +//! plane through the cells carrying their values, a field sampled at a surface beside the cells, +//! and streamlines traced through the velocity field. +//! +//! Values are cell-centred; a cell with no value (inside an obstacle) holds NaN. A slice joins the +//! cell centres of the layers either side of its plane (interpolated linearly between them) into +//! triangles carrying their corners' values, so [`crate::contour::split_triangle`] bands them as it +//! bands a nodal field; the outermost centres are moved out to the domain's faces, so the slice +//! covers the whole domain. A triangle touching cells with no value takes the mean of its other +//! corners there, and one with none is left out (the obstacle is drawn there). +//! +//! Streamlines follow the trilinear interpolation of the cell velocities between the cell centres +//! (solid cells at rest: the no-slip value), integrated along their arc length by the midpoint +//! (second-order Runge–Kutta) method with a step of half a cell. Each stops where it leaves the +//! domain, enters a solid cell or reaches still fluid, or at its step cap; all of them together +//! stop at a cap on points, so hostile counts cost a bounded time. + +use crate::{Aabb3, Vec3}; + +/// Most quadrilaterals a slice is made of: a larger layer is sampled every few cells. +pub const MAX_SLICE_QUADS: usize = 1 << 20; +/// Most streamlines traced at once. +pub const MAX_STREAMLINES: usize = 2_000; +/// Most steps of one streamline. +pub const MAX_STREAMLINE_STEPS: usize = 100_000; +/// Most points of all the streamlines traced at once. +pub const MAX_STREAMLINE_POINTS: usize = 2_000_000; +/// A streamline's step, in cells. +const STEP: f64 = 0.5; +/// A speed below this fraction of the field's largest is still fluid, where a streamline ends. +const STILL: f64 = 1e-6; + +/// A regular grid of cubic cells: cell `(i, j, k)` spans `origin + spacing · (i, j, k)` to one +/// spacing further, and its values sit at index `i + nx (j + ny k)`. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct CellGrid { + pub origin: [f64; 3], + pub spacing: f64, + pub cells: [usize; 3], +} + +impl CellGrid { + /// The number of cells, or `None` for a grid that cannot hold values: no cells, a spacing + /// that is not a positive number, an origin that is not finite, a count that overflows. + pub fn cell_count(&self) -> Option { + let ok = self.spacing.is_finite() && self.spacing > 0.0 && self.origin.iter().all(|c| c.is_finite()) && self.cells.iter().all(|&n| n > 0); + if !ok { + return None; + } + self.cells[0].checked_mul(self.cells[1])?.checked_mul(self.cells[2]) + } + + fn index(&self, [i, j, k]: [usize; 3]) -> usize { + i + self.cells[0] * (j + self.cells[1] * k) + } + + /// The `(i, j, k)` of cell index `c`. + pub fn coords(&self, c: usize) -> [usize; 3] { + let [nx, ny, _] = self.cells.map(|n| n.max(1)); + [c % nx, (c / nx) % ny, c / (nx * ny)] + } + + /// The centre of cell `(i, j, k)`. + pub fn centre(&self, ijk: [usize; 3]) -> Vec3 { + Vec3::new(self.centre_along(0, ijk[0]), self.centre_along(1, ijk[1]), self.centre_along(2, ijk[2])) + } + + fn centre_along(&self, axis: usize, i: usize) -> f64 { + self.start(axis) + (i as f64 + 0.5) * self.spacing + } + + /// The near end of the grid along `axis` (its origin there). + pub fn start(&self, axis: usize) -> f64 { + self.origin.get(axis).copied().unwrap_or(0.0) + } + + /// The far end of the grid along `axis`. + pub fn end(&self, axis: usize) -> f64 { + self.start(axis) + self.cells.get(axis).copied().unwrap_or(0) as f64 * self.spacing + } + + /// The box the cells fill. + pub fn bounds(&self) -> Aabb3 { + let mut b = Aabb3::EMPTY; + b.add(Vec3::from(self.origin)); + b.add(Vec3::new(self.end(0), self.end(1), self.end(2))); + b + } + + /// The cell holding `p`, if it is in the grid. + pub fn cell_at(&self, p: Vec3) -> Option { + let mut ijk = [0usize; 3]; + for (axis, slot) in ijk.iter_mut().enumerate() { + let n = self.cells.get(axis).copied().unwrap_or(0); + let t = (coord(p, axis) - self.start(axis)) / self.spacing; + if !(t >= 0.0 && t < n as f64) { + return None; + } + *slot = (t.floor() as usize).min(n.saturating_sub(1)); + } + Some(self.index(ijk)) + } + + /// The cells around `p` for trilinear interpolation between their centres: per axis the lower + /// index, the upper one and the upper one's weight. `None` when `p` lies more than `margin` + /// cells outside the grid. + fn stencil(&self, p: Vec3, margin: f64) -> Option<[(usize, usize, f64); 3]> { + let mut out = [(0, 0, 0.0); 3]; + for (axis, slot) in out.iter_mut().enumerate() { + let n = self.cells.get(axis).copied().unwrap_or(0); + let t = (coord(p, axis) - self.start(axis)) / self.spacing - 0.5; + if !(t >= -0.5 - margin && t <= n as f64 - 0.5 + margin) { + return None; + } + let t = t.clamp(0.0, n.saturating_sub(1) as f64); + let lo = t.floor() as usize; + let hi = (lo + 1).min(n.saturating_sub(1)); + *slot = (lo, hi, t - lo as f64); + } + Some(out) + } +} + +/// The eight cells of a trilinear stencil ([`CellGrid::stencil`]) with their weights. +fn corners(st: &[(usize, usize, f64); 3]) -> impl Iterator + '_ { + (0..8).map(move |corner| { + let mut ijk = [0; 3]; + let mut w = 1.0; + for (axis, (lo, hi, f)) in st.iter().enumerate() { + let up = (corner >> axis) & 1 == 1; + ijk[axis] = if up { *hi } else { *lo }; + w *= if up { *f } else { 1.0 - f }; + } + (ijk, w) + }) +} + +/// The coordinate of `p` along `axis` (0, 1, 2: x, y, z). +fn coord(p: Vec3, axis: usize) -> f64 { + match axis { + 0 => p.x, + 1 => p.y, + _ => p.z, + } +} + +/// A triangle carrying a value at each corner. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct ValueTri { + pub p: [Vec3; 3], + pub v: [f64; 3], +} + +/// A slice plane through a grid: its triangles with their corners' values, and its normal. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Slice { + pub normal: Vec3, + pub tris: Vec, +} + +/// The cell indices a slice samples along an axis of `n` cells, every `stride`-th and the last, +/// each with the coordinate its vertex is drawn at: the cell centre, the first and last moved out +/// to the grid's faces (`lo`, `hi`). One cell gives two vertices, one at each face. +fn samples(n: usize, stride: usize, lo: f64, hi: f64, centre: impl Fn(usize) -> f64) -> Vec<(usize, f64)> { + let mut out: Vec<(usize, f64)> = (0..n).step_by(stride.max(1)).map(|i| (i, centre(i))).collect(); + if let Some(last) = n.checked_sub(1) + && out.last().is_some_and(|(i, _)| *i != last) + { + out.push((last, centre(last))); + } + if out.len() == 1 + && let Some(&(i, _)) = out.first() + { + out.push((i, hi)); + } + if let Some(first) = out.first_mut() { + first.1 = lo; + } + if let Some(last) = out.last_mut() { + last.1 = hi; + } + out +} + +/// The slice of `values` (one per cell, NaN where there is none) on the plane normal to `axis` +/// (0, 1, 2: x, y, z) at coordinate `at`. `None` when the plane misses the grid, the axis is not +/// one, or the values do not match the grid. +pub fn slice(grid: &CellGrid, values: &[f64], axis: usize, at: f64) -> Option { + let n = grid.cell_count()?; + if values.len() != n || axis > 2 || !at.is_finite() { + return None; + } + let lo = grid.start(axis); + if !(at >= lo && at <= grid.end(axis)) { + return None; + } + let (u, v) = ((axis + 1) % 3, (axis + 2) % 3); + let (nw, nu, nv) = (grid.cells.get(axis).copied()?, grid.cells.get(u).copied()?, grid.cells.get(v).copied()?); + // The layers either side of the plane, and the upper one's weight. + let t = ((at - lo) / grid.spacing - 0.5).clamp(0.0, nw.saturating_sub(1) as f64); + let k0 = t.floor() as usize; + let k1 = (k0 + 1).min(nw.saturating_sub(1)); + let w = t - k0 as f64; + let value = |a: usize, b: usize| -> f64 { + let layer = |k: usize| { + let mut ijk = [0; 3]; + ijk[axis] = k; + ijk[u] = a; + ijk[v] = b; + values.get(grid.index(ijk)).copied().filter(|x| x.is_finite()).unwrap_or(f64::NAN) + }; + let (x0, x1) = (layer(k0), layer(k1)); + match (x0.is_nan(), x1.is_nan()) { + (false, false) => x0 + (x1 - x0) * w, + // Beside an obstacle: the nearer layer's value, if it has one. + (false, true) if w < 0.5 => x0, + (true, false) if w >= 0.5 => x1, + _ => f64::NAN, + } + }; + let stride = (nu.saturating_mul(nv) as f64 / MAX_SLICE_QUADS as f64).sqrt().ceil().max(1.0) as usize; + let su = samples(nu, stride, grid.start(u), grid.end(u), |i| grid.centre_along(u, i)); + let sv = samples(nv, stride, grid.start(v), grid.end(v), |i| grid.centre_along(v, i)); + let grid_values: Vec> = su.iter().map(|&(a, _)| sv.iter().map(|&(b, _)| value(a, b)).collect()).collect(); + let at_vertex = |a: usize, b: usize| -> Option<(Vec3, f64)> { + let (cu, cv) = (su.get(a)?.1, sv.get(b)?.1); + let mut c = [0.0; 3]; + c[axis] = at; + c[u] = cu; + c[v] = cv; + Some((Vec3::from(c), *grid_values.get(a)?.get(b)?)) + }; + let mut tris = Vec::new(); + for a in 0..su.len().saturating_sub(1) { + for b in 0..sv.len().saturating_sub(1) { + let (Some(c00), Some(c10), Some(c11), Some(c01)) = (at_vertex(a, b), at_vertex(a + 1, b), at_vertex(a + 1, b + 1), at_vertex(a, b + 1)) + else { + continue; + }; + let corners = [c00, c10, c11, c01]; + let (sum, count) = corners.iter().filter(|c| !c.1.is_nan()).fold((0.0, 0usize), |(s, n), c| (s + c.1, n + 1)); + if count == 0 { + continue; + } + let mean = sum / count as f64; + let [p0, p1, p2, p3] = corners.map(|c| c.0); + let [v0, v1, v2, v3] = corners.map(|c| if c.1.is_nan() { mean } else { c.1 }); + tris.push(ValueTri { p: [p0, p1, p2], v: [v0, v1, v2] }); + tris.push(ValueTri { p: [p0, p2, p3], v: [v0, v2, v3] }); + } + } + let mut normal = [0.0; 3]; + normal[axis] = 1.0; + Some(Slice { normal: Vec3::from(normal), tris }) +} + +/// The value of `values` (one per cell, NaN where there is none) at `p`, interpolated +/// trilinearly between the cell centres around it from those with a value; with none of them, the +/// nearest cell with a value within two cells. NaN outside the grid (more than a cell beyond it), +/// with nothing near, or when the values do not match the grid. For colouring an obstacle's +/// surface by the fluid beside it. +pub fn sample(grid: &CellGrid, values: &[f64], p: Vec3) -> f64 { + if grid.cell_count() != Some(values.len()) { + return f64::NAN; + } + let Some(st) = grid.stencil(p, 1.0) else { return f64::NAN }; + let at = |ijk: [usize; 3]| values.get(grid.index(ijk)).copied().filter(|x| x.is_finite()); + let (mut sum, mut weight) = (0.0, 0.0); + for (ijk, w) in corners(&st) { + if let Some(x) = at(ijk).filter(|_| w > 0.0) { + sum += w * x; + weight += w; + } + } + if weight > 1e-9 { + return sum / weight; + } + // The nearest cell with a value, two cells around the nearest centre at most. + let near = st.map(|(lo, hi, f)| if f < 0.5 { lo } else { hi }); + let span = |axis: usize| { + let n = near.get(axis).copied().unwrap_or(0); + n.saturating_sub(2)..=(n + 2).min(grid.cells.get(axis).copied().unwrap_or(1).saturating_sub(1)) + }; + let mut best: Option<(f64, f64)> = None; + for k in span(2) { + for j in span(1) { + for i in span(0) { + let Some(x) = at([i, j, k]) else { continue }; + let d = grid.centre([i, j, k]).dist(p); + if best.is_none_or(|(bd, _)| d < bd) { + best = Some((d, x)); + } + } + } + } + best.map_or(f64::NAN, |(_, x)| x) +} + +/// A traced streamline: its points from the seed downstream and the speed at each. +#[derive(Clone, Debug, Default, PartialEq)] +pub struct Streamline { + pub points: Vec, + pub speed: Vec, +} + +/// The velocity at `p` interpolated trilinearly between the cell centres (solid cells at rest); +/// `None` outside the grid. +fn velocity_at(grid: &CellGrid, velocity: &[[f32; 3]], solid: &[bool], p: Vec3) -> Option { + let st = grid.stencil(p, 0.0)?; + let mut u = Vec3::ZERO; + for (ijk, w) in corners(&st) { + let c = grid.index(ijk); + if w == 0.0 || solid.get(c).copied().unwrap_or(false) { + continue; + } + let Some(v) = velocity.get(c) else { continue }; + let v = Vec3::new(f64::from(v[0]), f64::from(v[1]), f64::from(v[2])); + if v.is_finite() { + u = u + v * w; + } + } + Some(u) +} + +/// Streamlines of `velocity` (one per cell; `solid` flags the cells inside obstacles, or is empty) +/// from each seed (at most [`MAX_STREAMLINES`]) downstream, each at most `max_steps` steps of half +/// a cell (capped at [`MAX_STREAMLINE_STEPS`]), all together at most [`MAX_STREAMLINE_POINTS`] +/// points. A seed outside the fluid gives no streamline; nothing is traced when the fields do not +/// match the grid. +pub fn streamlines(grid: &CellGrid, velocity: &[[f32; 3]], solid: &[bool], seeds: &[Vec3], max_steps: usize) -> Vec { + let Some(n) = grid.cell_count() else { return Vec::new() }; + if velocity.len() != n || !(solid.is_empty() || solid.len() == n) { + return Vec::new(); + } + let fastest = + velocity.iter().map(|v| Vec3::new(f64::from(v[0]), f64::from(v[1]), f64::from(v[2])).len()).filter(|s| s.is_finite()).fold(0.0, f64::max); + if !(fastest > 0.0) { + return Vec::new(); + } + let still = STILL * fastest; + let h = STEP * grid.spacing; + // The flow's direction and speed at a point in the fluid. + let dir = |p: Vec3| -> Option<(Vec3, f64)> { + let cell = grid.cell_at(p)?; + if solid.get(cell).copied().unwrap_or(false) { + return None; + } + let u = velocity_at(grid, velocity, solid, p)?; + let s = u.len(); + (s > still).then(|| (u * (1.0 / s), s)) + }; + let max_steps = max_steps.min(MAX_STREAMLINE_STEPS); + let mut budget = MAX_STREAMLINE_POINTS; + let mut out = Vec::new(); + for &seed in seeds.iter().take(MAX_STREAMLINES) { + if budget == 0 { + break; + } + let Some((mut d, s)) = dir(seed) else { continue }; + let mut line = Streamline { points: vec![seed], speed: vec![s] }; + budget -= 1; + let mut p = seed; + for _ in 0..max_steps { + if budget == 0 { + break; + } + let Some((mid, _)) = dir(p + d * (0.5 * h)) else { break }; + let q = p + mid * h; + let Some((next, s)) = dir(q) else { break }; + line.points.push(q); + line.speed.push(s); + budget -= 1; + (p, d) = (q, next); + } + if line.points.len() >= 2 { + out.push(line); + } + } + out +} + +/// `count` seeds (at most [`MAX_STREAMLINES`]) on a lattice across the grid's face at the low +/// end of `axis` (`high`: the high end), half a cell inside it: over the solid cells seen along the +/// axis, the box around them grown by half its size on every side and kept on the face, so the +/// streamlines pass over and around the obstacles; over the whole face without solid cells. +/// Empty for an axis that is not one or a grid that holds no cells. +pub fn face_seeds(grid: &CellGrid, solid: &[bool], axis: usize, high: bool, count: usize) -> Vec { + let Some(n) = grid.cell_count() else { return Vec::new() }; + let count = count.min(MAX_STREAMLINES); + if axis > 2 || count == 0 { + return Vec::new(); + } + let (u, v) = ((axis + 1) % 3, (axis + 2) % 3); + let face = |a: usize| (grid.start(a), grid.end(a)); + let (mut ru, mut rv) = (face(u), face(v)); + // The solid cells' extent across the face, by cell index. + let mut box_ = [usize::MAX, 0, usize::MAX, 0]; + for (c, _) in solid.iter().enumerate().take(n).filter(|(_, s)| **s) { + let ijk = grid.coords(c); + let (a, b) = (ijk[u], ijk[v]); + box_ = [box_[0].min(a), box_[1].max(a), box_[2].min(b), box_[3].max(b)]; + } + if box_[0] <= box_[1] { + let grow = |(lo, hi): (f64, f64), first: usize, last: usize, axis: usize| { + let a = grid.start(axis) + first as f64 * grid.spacing; + let b = grid.start(axis) + (last + 1) as f64 * grid.spacing; + let m = 0.5 * (b - a); + ((a - m).max(lo), (b + m).min(hi)) + }; + ru = grow(ru, box_[0], box_[1], u); + rv = grow(rv, box_[2], box_[3], v); + } + let (wu, wv) = (ru.1 - ru.0, rv.1 - rv.0); + if !(wu > 0.0 && wv > 0.0) { + return Vec::new(); + } + let cu = ((count as f64 * wu / wv).sqrt().round() as usize).clamp(1, count); + let cv = (count / cu).max(1); + let at = if high { grid.end(axis) - 0.5 * grid.spacing } else { grid.start(axis) + 0.5 * grid.spacing }; + let mut out = Vec::with_capacity(cu * cv); + for a in 0..cu { + for b in 0..cv { + let mut c = [0.0; 3]; + c[axis] = at; + c[u] = ru.0 + (a as f64 + 0.5) / cu as f64 * wu; + c[v] = rv.0 + (b as f64 + 0.5) / cv as f64 * wv; + out.push(Vec3::from(c)); + } + } + out +} + +#[cfg(test)] +mod tests { + use super::*; + + fn grid(n: [usize; 3]) -> CellGrid { + CellGrid { origin: [0.0; 3], spacing: 1.0, cells: n } + } + + /// A field linear in x (the value is the cell centre's x) slices to a plane with exact values + /// at its corners, covering the whole domain. + #[test] + fn a_slice_covers_the_domain_and_interpolates() { + let g = grid([4, 3, 2]); + let values: Vec = (0..g.cell_count().unwrap()).map(|c| (c % 4) as f64 + 0.5).collect(); + let s = slice(&g, &values, 2, 1.0).unwrap(); + assert_eq!(s.normal, Vec3::Z); + assert_eq!(s.tris.len(), 2 * 3 * 2, "4 × 3 samples, 3 × 2 quads"); + let mut b = Aabb3::EMPTY; + for t in &s.tris { + for (p, v) in t.p.iter().zip(t.v) { + b.add(*p); + assert_eq!(p.z, 1.0); + // Inside, the value is x; at the faces, the outer cells' values. + let expect = p.x.clamp(0.5, 3.5); + assert!((v - expect).abs() < 1e-12, "{p:?} {v}"); + } + } + assert_eq!((b.min, b.max), (Vec3::new(0.0, 0.0, 1.0), Vec3::new(4.0, 3.0, 1.0))); + // Between two layers the values are interpolated: z-dependent field, plane at z = 1. + let vz: Vec = (0..g.cell_count().unwrap()).map(|c| (c / 12) as f64 * 10.0).collect(); + let s = slice(&g, &vz, 2, 1.0).unwrap(); + assert!(s.tris.iter().all(|t| t.v.iter().all(|v| (v - 5.0).abs() < 1e-12))); + // Along x and y too; one cell along an axis still makes a quad across it. + assert!(slice(&g, &values, 0, 2.0).is_some_and(|s| s.normal == Vec3::X && s.tris.len() == 2 * 2)); + assert!(slice(&grid([1, 1, 1]), &[3.0], 1, 0.5).is_some_and(|s| s.tris.len() == 2)); + } + + /// Cells with no value: a quad touching some takes its other corners' mean, one with none is + /// left out; hostile planes and fields are refused. + #[test] + fn a_slice_skips_solid_cells_and_refuses_bad_input() { + let g = grid([3, 3, 1]); + let mut values = vec![1.0; 9]; + values[4] = f64::NAN; + let s = slice(&g, &values, 2, 0.5).unwrap(); + assert_eq!(s.tris.len(), 8); + assert!(s.tris.iter().all(|t| t.v == [1.0; 3])); + let all_solid = slice(&g, &[f64::NAN; 9], 2, 0.5).unwrap(); + assert!(all_solid.tris.is_empty()); + for (axis, at) in [(3, 0.5), (2, -0.1), (2, 1.1), (2, f64::NAN), (2, f64::INFINITY)] { + assert!(slice(&g, &values, axis, at).is_none(), "{axis} {at}"); + } + assert!(slice(&g, &values[..8], 2, 0.5).is_none()); + let bad = [ + CellGrid { spacing: 0.0, ..g }, + CellGrid { spacing: f64::NAN, ..g }, + CellGrid { origin: [f64::INFINITY, 0.0, 0.0], ..g }, + CellGrid { cells: [0, 3, 1], ..g }, + CellGrid { cells: [usize::MAX, usize::MAX, 2], ..g }, + ]; + for b in bad { + assert!(slice(&b, &values, 2, 0.5).is_none(), "{b:?}"); + assert!(sample(&b, &values, Vec3::ZERO).is_nan()); + assert!(streamlines(&b, &[[1.0, 0.0, 0.0]; 9], &[], &[Vec3::new(0.5, 0.5, 0.5)], 10).is_empty()); + assert!(face_seeds(&b, &[], 0, false, 10).is_empty()); + } + } + + /// A huge layer is sampled every few cells, keeping its extent. + #[test] + fn a_large_slice_is_sampled() { + let g = grid([2048, 1024, 1]); + let values = vec![1.0; g.cell_count().unwrap()]; + let s = slice(&g, &values, 2, 0.5).unwrap(); + assert!(s.tris.len() <= 2 * MAX_SLICE_QUADS, "{}", s.tris.len()); + let xs = s.tris.iter().flat_map(|t| t.p).map(|p| p.x).fold(0.0, f64::max); + assert_eq!(xs, 2048.0); + } + + /// Sampling beside a solid block takes the fluid cells' values only. + #[test] + fn sampling_beside_solid_cells_reads_the_fluid() { + let g = grid([4, 1, 1]); + let values = [f64::NAN, f64::NAN, 2.0, 4.0]; + // On the block's face at x = 2: cells 1 (solid) and 2. + assert_eq!(sample(&g, &values, Vec3::new(2.0, 0.5, 0.5)), 2.0); + assert_eq!(sample(&g, &values, Vec3::new(3.0, 0.5, 0.5)), 3.0); + // Deep inside the block, the nearest fluid cell; beyond the grid, nothing. + assert_eq!(sample(&g, &values, Vec3::new(0.2, 0.5, 0.5)), 2.0); + assert!(sample(&g, &values, Vec3::new(9.0, 0.5, 0.5)).is_nan()); + assert!(sample(&g, &[f64::NAN; 4], Vec3::new(1.0, 0.5, 0.5)).is_nan()); + assert!(sample(&g, &values, Vec3::new(f64::NAN, 0.5, 0.5)).is_nan()); + } + + /// A uniform flow along x carries a streamline straight to the outlet in half-cell steps; a + /// solid cell stops it; still fluid and seeds outside give none. + #[test] + fn streamlines_follow_the_flow() { + let g = grid([10, 3, 3]); + let n = g.cell_count().unwrap(); + let u = vec![[2.0f32, 0.0, 0.0]; n]; + let seed = Vec3::new(0.5, 1.5, 1.5); + let lines = streamlines(&g, &u, &[], &[seed], 1000); + assert_eq!(lines.len(), 1); + let l = &lines[0]; + assert!(l.points.iter().all(|p| (p.y - 1.5).abs() < 1e-12 && (p.z - 1.5).abs() < 1e-12)); + let last = *l.points.last().unwrap(); + assert!(last.x > 9.4 && last.x < 10.0, "{last:?}"); + assert_eq!(l.points.len(), 19); + assert!(l.speed.iter().all(|s| (s - 2.0).abs() < 1e-6)); + // A step cap. + assert_eq!(streamlines(&g, &u, &[], &[seed], 3)[0].points.len(), 4); + // A solid cell across the way at x = 5. + let mut solid = vec![false; n]; + solid[g.index([5, 1, 1])] = true; + let blocked = streamlines(&g, &u, &solid, &[seed], 1000); + assert!(blocked[0].points.iter().all(|p| p.x < 5.0), "{:?}", blocked[0].points.last()); + // Outside, in a solid cell, in still fluid: no streamline. + let none = [Vec3::new(-1.0, 1.5, 1.5), Vec3::new(5.5, 1.5, 1.5), Vec3::new(f64::NAN, 0.0, 0.0)]; + assert!(streamlines(&g, &u, &solid, &none, 1000).is_empty()); + assert!(streamlines(&g, &vec![[0.0; 3]; n], &[], &[seed], 1000).is_empty()); + assert!(streamlines(&g, &u[..5], &[], &[seed], 1000).is_empty()); + assert!(streamlines(&g, &u, &solid[..5], &[seed], 1000).is_empty()); + } + + /// A rigid rotation about z: streamlines are circles, the midpoint rule keeps their radius. + #[test] + fn streamlines_of_a_vortex_close() { + let g = CellGrid { origin: [-10.0, -10.0, 0.0], spacing: 0.5, cells: [40, 40, 2] }; + let u: Vec<[f32; 3]> = (0..g.cell_count().unwrap()) + .map(|c| { + let p = g.centre([c % 40, (c / 40) % 40, c / 1600]); + [-p.y as f32, p.x as f32, 0.0] + }) + .collect(); + let lines = streamlines(&g, &u, &[], &[Vec3::new(5.0, 0.0, 0.5)], 2000); + let l = &lines[0]; + assert_eq!(l.points.len(), 2001, "a closed streamline runs to its step cap"); + for p in &l.points { + assert!(((p.x * p.x + p.y * p.y).sqrt() - 5.0).abs() < 0.05, "{p:?}"); + } + } + + /// Hostile counts are capped: seeds, steps and points. + #[test] + fn streamline_work_is_bounded() { + let g = grid([64, 64, 64]); + let n = g.cell_count().unwrap(); + let u = vec![[0.0f32, 0.0, 1.0]; n]; + let seeds = face_seeds(&g, &[], 2, false, usize::MAX); + assert_eq!(seeds.len(), 44 * 45, "capped at {MAX_STREAMLINES}"); + assert!(seeds.len() <= MAX_STREAMLINES); + let lines = streamlines(&g, &u, &[], &seeds, usize::MAX); + let points: usize = lines.iter().map(|l| l.points.len()).sum(); + assert!(points <= MAX_STREAMLINE_POINTS); + assert!(lines.iter().all(|l| l.points.len() == 127), "each runs up the 64 cells in 126 steps"); + } + + /// Seeds cover the face, or the obstacle's shadow on it grown by half its size each way. + #[test] + fn seeds_lie_upstream_of_the_obstacle() { + let g = grid([10, 10, 10]); + let all = face_seeds(&g, &[], 0, false, 25); + assert_eq!(all.len(), 25); + assert!(all.iter().all(|p| p.x == 0.5 && (0.0..10.0).contains(&p.y) && (0.0..10.0).contains(&p.z))); + assert_eq!(all[0], Vec3::new(0.5, 1.0, 1.0)); + // A 2 × 2 block of cells at y, z in [4, 6]: seeds over y, z in [3, 7]. + let mut solid = vec![false; 1000]; + for (j, k) in [(4, 4), (4, 5), (5, 4), (5, 5)] { + solid[g.index([7, j, k])] = true; + } + let near = face_seeds(&g, &solid, 0, false, 16); + assert_eq!(near.len(), 16); + assert!(near.iter().all(|p| (3.0..=7.0).contains(&p.y) && (3.0..=7.0).contains(&p.z)), "{near:?}"); + let high = face_seeds(&g, &solid, 0, true, 4); + assert!(high.iter().all(|p| p.x == 9.5)); + assert!(face_seeds(&g, &solid, 3, false, 4).is_empty()); + assert!(face_seeds(&g, &solid, 0, false, 0).is_empty()); + // A wide face gets more seeds across its width. + let wide = face_seeds(&grid([1, 40, 10]), &[], 0, false, 16); + assert_eq!(wide.len(), 8 * 2); + } +} diff --git a/crates/render/src/lib.rs b/crates/render/src/lib.rs index 4301f79..6cbedf3 100644 --- a/crates/render/src/lib.rs +++ b/crates/render/src/lib.rs @@ -1,7 +1,8 @@ //! Display data for FieldCraft, with no GPU dependency: what a finite element model looks like //! on screen ([`DisplayMesh`]: the visible faces, element edges, line elements and nodes), an //! orbit camera with perspective or orthographic projection, picking ([`pick`]: what a ray through -//! a pixel hits), and a CPU rasterizer for headless snapshots. The desktop viewport draws the +//! a pixel hits), flow fields on a grid of cells ([`flow`]: slices, surface values and +//! streamlines), and a CPU rasterizer for headless snapshots. The desktop viewport draws the //! same [`Scene`] on the GPU. #![deny(clippy::unwrap_used, clippy::expect_used, clippy::panic, clippy::unimplemented, clippy::todo, clippy::unreachable)] #![forbid(unsafe_code)] @@ -9,6 +10,7 @@ mod camera; pub mod contour; mod display; +pub mod flow; mod font; pub mod glyphs; pub mod pick; diff --git a/docs/architecture.md b/docs/architecture.md index 589dec7..ef2cf26 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -24,7 +24,7 @@ crates in strictly lower layers. | L1 | `fieldcraft-mesh` | Surface preparation (weld, watertightness, feature edges), robust predicates, Delaunay tetrahedralization and refinement, Voronoi-based optimization, hex meshing, quality metrics. | started | | L1 | `fieldcraft-fem` | Element kernels shared by every structural solver: the element table (tet4, tet10, wedge6, hex8 full and incompatible modes, hex20 full and reduced), shape functions, quadrature, element stiffness, consistent and lumped mass, body and pressure loads, stresses and their extrapolation, isotropic and orthotropic elasticity, the uniform-strain (one-point) hex8 and tet4 kinematics, hourglass shape vectors and element frequency bounds of the explicit solver, objective stress rates and J2/Johnson–Cook return mapping. Pure functions, no solver loop. | started | | L1 | `fieldcraft-decks` | Abaqus input, Nastran bulk data and LS-DYNA keyword readers and writers, to and from `model`. Text out, no file system; `deck.write` saves it. | started | -| L1 | `fieldcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots, geometry surfaces with their face boundaries, deformed shapes and banded contour plots of nodal results, picking. | started | +| L1 | `fieldcraft-render` | Display data without a GPU dependency: the visible faces, element edges and nodes of a model, the orbit camera, a CPU rasterizer for headless snapshots, geometry surfaces with their face boundaries, deformed shapes and banded contour plots of nodal results, slices, surface values and streamlines of flow fields on grids of cells, picking. | started | | L2 | `fieldcraft-solver` | Implicit solver: assembly from the resolved model, boundary conditions, sparse direct solution, linear static and modal analysis, stresses and reactions. | started | | L2 | `fieldcraft-explicit` | Explicit dynamics: central-difference time integration with lumped mass, one-point hex8 with hourglass control, nodal-pressure tet4, bulk viscosity, rigid walls, large deformation with objective rates, plasticity, damage and erosion, energy account, histories and field frames. | started | | L2 | `fieldcraft-lbm` | Lattice Boltzmann flow solver on the CPU: D3Q19 two-relaxation-time scheme, bounce-back and slip walls, Zou–He inlets and outlets, obstacles voxelised from the face-tagged surface with exact predicates, momentum-exchange forces per face. Uses `mesh` for the predicates and the boundary surface. | started |