From 700a7e3d10a4cd94b1d03fc8c918ed749c09b5ed Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Fri, 16 Oct 2026 17:12:34 -0400 Subject: [PATCH] mesh: platform-independent maths (libm) so meshes are bit-identical on macOS and Linux MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The system maths library's sin, cos, atan2, log … are not correctly rounded and differ between platforms: mesh.hex.generate gave the plate of the strong-relax test a different mesh on Linux than on macOS (nodes from 2809 on), so the test passed on a Mac and failed in CI. Every transcendental call in fieldcraft-mesh now goes through the libm crate (a pure-Rust port of musl's), which gives the same bits everywhere. Part of #438. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- Cargo.lock | 1 + crates/mesh/Cargo.toml | 4 ++++ crates/mesh/src/block.rs | 8 ++++---- crates/mesh/src/hex/brep.rs | 2 +- crates/mesh/src/hex/decompose.rs | 13 ++++++++----- crates/mesh/src/hex/extrude_to.rs | 2 +- crates/mesh/src/hex/geom.rs | 6 +++--- crates/mesh/src/hex/layered.rs | 6 +++--- crates/mesh/src/hex/mod.rs | 2 +- crates/mesh/src/hex/region.rs | 9 ++++++--- crates/mesh/src/octree/mod.rs | 6 +++--- crates/mesh/src/planar.rs | 2 +- crates/mesh/src/project/mod.rs | 6 +++--- crates/mesh/src/relax/strong.rs | 2 +- crates/mesh/src/surface/topology.rs | 2 +- 15 files changed, 41 insertions(+), 30 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 68587a2..606c60b 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -1418,6 +1418,7 @@ name = "fieldcraft-mesh" version = "0.6.0" dependencies = [ "fieldcraft-model", + "libm", "num-bigint", "num-rational", "num-traits", diff --git a/crates/mesh/Cargo.toml b/crates/mesh/Cargo.toml index 948375e..806c14f 100644 --- a/crates/mesh/Cargo.toml +++ b/crates/mesh/Cargo.toml @@ -10,6 +10,10 @@ workspace = true [dependencies] fieldcraft-model = { workspace = true } +# Pure-Rust maths (a port of musl's libm) for sin, cos, log, exp, pow, hypot, acos: the same +# bits on every platform, where the system library's are not correctly rounded and differ (the +# plate's hex mesh differed between macOS and Linux, #438). atan2 and cbrt are crate::vec3's. +libm = "0.2" thiserror = { workspace = true } [dev-dependencies] diff --git a/crates/mesh/src/block.rs b/crates/mesh/src/block.rs index daf8afd..ee69caa 100644 --- a/crates/mesh/src/block.rs +++ b/crates/mesh/src/block.rs @@ -146,7 +146,7 @@ impl Path { if norm(sub(back, p)) > 1e-9 * la.max(lb) || ((c * c + s * s).sqrt() - 1.0).abs() > 1e-9 { return Err(Error::Invalid(format!("corner {p:?} is not on the edge's ellipse or arc about {center:?}"))); } - Ok(s.atan2(c)) + Ok(crate::vec3::atan2(s, c)) }; let start = angle(from)?; let mut sweep = angle(to)? - start; @@ -170,7 +170,7 @@ impl Path { Path::Line { from, to } => add_scaled(from, s, sub(to, from)), Path::Ellipse { center, a, b, start, sweep, .. } => { let t = start + s * sweep; - add_scaled(add_scaled(center, t.cos(), a), t.sin(), b) + add_scaled(add_scaled(center, libm::cos(t), a), libm::sin(t), b) } } } @@ -179,7 +179,7 @@ impl Path { /// The parameters of `n` elements along one direction with last/first size ratio `bias`, at /// `order` points per element (`order · n + 1` values from 0 to 1; mid-side values halfway). fn parameters(n: usize, bias: f64, order: usize) -> Vec { - let q = if n > 1 { bias.powf(1.0 / (n - 1) as f64) } else { 1.0 }; + let q = if n > 1 { libm::pow(bias, 1.0 / (n - 1) as f64) } else { 1.0 }; let mut ends = Vec::with_capacity(n + 1); let (mut sum, mut size) = (0.0, 1.0); ends.push(0.0); @@ -645,7 +645,7 @@ mod tests { ]; let m = block(&Block { kind: ElementKind::Hex20, corners, divisions: [3, 6, 2], bias: [1.0; 3], curves }).unwrap(); for p in &m.points { - let r = p[0].hypot(p[1]); + let r = libm::hypot(p[0], p[1]); assert!(r > a - 1e-12 && r < b + 1e-12, "{p:?}"); } let exact = std::f64::consts::PI / 4.0 * (b * b - a * a) * 0.5; diff --git a/crates/mesh/src/hex/brep.rs b/crates/mesh/src/hex/brep.rs index 3318d1b..a3f3435 100644 --- a/crates/mesh/src/hex/brep.rs +++ b/crates/mesh/src/hex/brep.rs @@ -372,7 +372,7 @@ fn composite_tags(surface: &BoundarySurface, surfaces: &[Option]) -> Ve unit(cross(sub(b, a), sub(c, a))) }) }; - let cos_min = SMOOTH_ANGLE.cos(); + let cos_min = libm::cos(SMOOTH_ANGLE); // Per pair of faces (low tag, high tag): whether every segment between them is smooth. let mut pairs: BTreeMap<(u32, u32), bool> = BTreeMap::new(); for fe in surface.feature_edges().iter().filter(|e| e.kind == FeatureKind::FaceBoundary) { diff --git a/crates/mesh/src/hex/decompose.rs b/crates/mesh/src/hex/decompose.rs index b777866..ee29dea 100644 --- a/crates/mesh/src/hex/decompose.rs +++ b/crates/mesh/src/hex/decompose.rs @@ -166,7 +166,7 @@ fn sub(a: Q, b: Q) -> Q { } fn len(a: Q) -> f64 { - a[0].hypot(a[1]) + libm::hypot(a[0], a[1]) } fn cross(a: Q, b: Q) -> f64 { @@ -186,14 +186,17 @@ fn unit(a: Q) -> Option { /// region on the left: π on a straight line, less at a convex corner, more at a concave one. fn interior(a: Q, o: Q, b: Q) -> f64 { let (u, v) = (sub(o, a), sub(b, o)); - std::f64::consts::PI - cross(u, v).atan2(dot(u, v)) + std::f64::consts::PI - crate::vec3::atan2(cross(u, v), dot(u, v)) } /// The inward bisector at a corner of interior angle `theta` whose outgoing edge runs along unit /// `w`: `w` turned left by half the angle. fn bisector(w: Q, theta: f64) -> Q { let half = 0.5 * theta; - [w[0] * half.cos() - w[1] * half.sin(), w[0] * half.sin() + w[1] * half.cos()] + { + let (s, c) = (libm::sin(half), libm::cos(half)); + [w[0] * c - w[1] * s, w[0] * s + w[1] * c] + } } fn centroid(pts: &[Q]) -> Q { @@ -1320,7 +1323,7 @@ mod tests { /// The half crank's cross-section (y ≤ 0) as the extrusion projects it: the outline, the /// bosses' circles where they rise above the arm, and the pocket's outline. fn crank_half() -> (Profile, Vec>) { - let deg = |p: Q, c: Q| (p[1] - c[1]).atan2(p[0] - c[0]).to_degrees(); + let deg = |p: Q, c: Q| crate::vec3::atan2(p[1] - c[1], p[0] - c[0]).to_degrees(); let spindle = [0.0, 0.0]; let pedal = [170.0, 0.0]; let j1 = [13.455063072993756, -18.65372021077189]; @@ -1344,7 +1347,7 @@ mod tests { }; let outer = vec![ a(spindle, 23.0, [-23.0, 0.0], j1, true), - a(c1, (j1[0] - c1[0]).hypot(j1[1] - c1[1]), j1, r1, false), + a(c1, libm::hypot(j1[0] - c1[0], j1[1] - c1[1]), j1, r1, false), line(r1, r2), a(c2, 6.0, r2, j2, false), a(pedal, 15.0, j2, [185.0, 0.0], true), diff --git a/crates/mesh/src/hex/extrude_to.rs b/crates/mesh/src/hex/extrude_to.rs index 9b062e7..7132d83 100644 --- a/crates/mesh/src/hex/extrude_to.rs +++ b/crates/mesh/src/hex/extrude_to.rs @@ -410,7 +410,7 @@ fn section_boundary(quads: &[[u32; 4]]) -> Result { /// The layer parameters `0 = t_0 < … < t_n = 1`, thicknesses in geometric progression with the /// last over the first equal to `bias`. pub(crate) fn parameters(n: usize, bias: f64) -> Vec { - let r = if n > 1 { bias.powf(1.0 / (n - 1) as f64) } else { 1.0 }; + let r = if n > 1 { libm::pow(bias, 1.0 / (n - 1) as f64) } else { 1.0 }; let mut out = Vec::with_capacity(n + 1); let (mut sum, mut h) = (0.0, 1.0); out.push(0.0); diff --git a/crates/mesh/src/hex/geom.rs b/crates/mesh/src/hex/geom.rs index 0c1e98c..26f2306 100644 --- a/crates/mesh/src/hex/geom.rs +++ b/crates/mesh/src/hex/geom.rs @@ -101,7 +101,7 @@ impl Surface { let h = dot(w, axis); let radial = unit(sub(w, scale(axis, h)))?; let rho = norm(sub(w, scale(axis, h))); - let (c, s) = (half_angle.cos(), half_angle.sin()); + let (c, s) = (libm::cos(half_angle), libm::sin(half_angle)); let t = (h * c + rho * s).max(0.0); add(apex, add(scale(axis, t * c), scale(radial, t * s))) } @@ -128,7 +128,7 @@ impl Surface { Surface::Cone { apex, axis, half_angle } => { let w = sub(p, apex); let radial = unit(sub(w, scale(axis, dot(w, axis))))?; - unit(sub(scale(radial, half_angle.cos()), scale(axis, half_angle.sin()))) + unit(sub(scale(radial, libm::cos(half_angle)), scale(axis, libm::sin(half_angle)))) } Surface::Sphere { center, .. } => unit(sub(p, center)), Surface::Torus { center, axis, major, .. } => { @@ -383,7 +383,7 @@ impl<'a> Inside<'a> { let (a, b, c) = (sub(*a, x), sub(*b, x), sub(*c, x)); let (la, lb, lc) = (norm(a), norm(b), norm(c)); let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la; - sum += 2.0 * dot(a, cross(b, c)).atan2(den); + sum += 2.0 * crate::vec3::atan2(dot(a, cross(b, c)), den); } Ok(sum / (4.0 * std::f64::consts::PI)) } diff --git a/crates/mesh/src/hex/layered.rs b/crates/mesh/src/hex/layered.rs index 167bb41..0af417a 100644 --- a/crates/mesh/src/hex/layered.rs +++ b/crates/mesh/src/hex/layered.rs @@ -98,8 +98,8 @@ impl Frame { let w = sub(p, o); let z = dot(w, a); let (x, y) = (dot(w, e1), dot(w, e2)); - let t = y.atan2(x); - ([z, x.hypot(y)], if t < 0.0 { t + TAU } else { t }) + let t = crate::vec3::atan2(y, x); + ([z, libm::hypot(x, y)], if t < 0.0 { t + TAU } else { t }) } } } @@ -107,7 +107,7 @@ impl Frame { pub(crate) fn point(&self, q: Q, h: f64) -> P { match *self { Frame::Translate { o, u, v, d } => add(o, add(add(scale(u, q[0]), scale(v, q[1])), scale(d, h))), - Frame::Rotate { o, a, e1, e2 } => add(o, add(scale(a, q[0]), scale(add(scale(e1, h.cos()), scale(e2, h.sin())), q[1]))), + Frame::Rotate { o, a, e1, e2 } => add(o, add(scale(a, q[0]), scale(add(scale(e1, libm::cos(h)), scale(e2, libm::sin(h))), q[1]))), } } diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index 87e7564..4db9c5b 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -479,7 +479,7 @@ impl CellGrid { } (lo, hi) }); - let per_axis = ((ctxs.len() as f64).cbrt().ceil() as usize * 2).clamp(1, MAX_GRID_DIM); + let per_axis = (crate::vec3::cbrt(ctxs.len() as f64).ceil() as usize * 2).clamp(1, MAX_GRID_DIM); let extent = (0..3).map(|k| hi[k] - lo[k]).fold(0.0, f64::max); let step = if extent > 0.0 && extent.is_finite() { extent / per_axis as f64 } else { 1.0 }; let dim = |k: usize| (((hi[k] - lo[k]) / step).ceil() as usize).clamp(1, MAX_GRID_DIM); diff --git a/crates/mesh/src/hex/region.rs b/crates/mesh/src/hex/region.rs index df026b3..c84d611 100644 --- a/crates/mesh/src/hex/region.rs +++ b/crates/mesh/src/hex/region.rs @@ -68,7 +68,7 @@ fn angles(lp: &[Vec], ccw: bool) -> Option> { let next = lp.get(k)?; let (a, o, b) = (*prev.get(prev.len().checked_sub(2)?)?, *next.first()?, *next.get(1)?); let (u, v) = ([o[0] - a[0], o[1] - a[1]], [b[0] - o[0], b[1] - o[1]]); - let turn = (u[0] * v[1] - u[1] * v[0]).atan2(u[0] * v[0] + u[1] * v[1]); + let turn = crate::vec3::atan2(u[0] * v[1] - u[1] * v[0], u[0] * v[0] + u[1] * v[1]); out.push(PI - if ccw { turn } else { -turn }); } Some(out) @@ -443,8 +443,11 @@ fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result P { /// The grid axes (rows): the coordinate axes, or the normals of the two largest orthogonal planar /// patches when more boundary area is aligned with them. fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] { - let tol = 1f64.to_radians().sin(); + let tol = libm::sin(1f64.to_radians()); let normal = |t: u32| g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); let mut planar: Vec<(f64, usize, P)> = Vec::new(); for (i, p) in g.patches.iter().enumerate() { @@ -501,7 +501,7 @@ fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] { /// that only measures how near the edge it started; the faces of a wall face each other. fn triangle_thickness(g: &Geometry) -> Vec { let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0; - let opposed = -WALL_OPENING.to_radians().cos(); + let opposed = -libm::cos(WALL_OPENING.to_radians()); (0..g.surface.triangles().len() as u32) .map(|t| { let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); @@ -625,7 +625,7 @@ impl Context<'_> { return Ok(true); } let reach = 0.5 * 3f64.sqrt() * e * 1.05; - let cos = self.options.curvature_angle.to_radians().cos(); + let cos = libm::cos(self.options.curvature_angle.to_radians()); let mut first: BTreeMap = BTreeMap::new(); for t in self.g.overlapping(&bx) { let [a, b, c] = self.g.triangle(t); diff --git a/crates/mesh/src/planar.rs b/crates/mesh/src/planar.rs index 8a69cea..d2566b6 100644 --- a/crates/mesh/src/planar.rs +++ b/crates/mesh/src/planar.rs @@ -266,7 +266,7 @@ pub fn face_boundary(positions: &[[f64; 3]], triangles: &[[u32; 3]], tri_face: & .map(|i| { let (p, q, r) = (at(i + n - 1), at(i), at(i + 1)); let (u, w) = ([q[0] - p[0], q[1] - p[1]], [r[0] - q[0], r[1] - q[1]]); - let turn = (u[0] * w[1] - u[1] * w[0]).atan2(u[0] * w[0] + u[1] * w[1]).abs(); + let turn = crate::vec3::atan2(u[0] * w[1] - u[1] * w[0], u[0] * w[0] + u[1] * w[1]).abs(); ids[i] != ids[(i + n - 1) % n] || turn > CORNER_ANGLE }) .collect(); diff --git a/crates/mesh/src/project/mod.rs b/crates/mesh/src/project/mod.rs index 3e32313..cc79340 100644 --- a/crates/mesh/src/project/mod.rs +++ b/crates/mesh/src/project/mod.rs @@ -551,7 +551,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o *m = m.max(measured); } if let Some(l) = class_log.get_mut(cls[d]).filter(|_| len > tol) { - l.0 += measured.ln(); + l.0 += libm::log(measured); l.1 += 1; } } @@ -598,7 +598,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o let Some(nx) = g.normal(fa, qx) else { continue }; let h = dist(pa, qx); if h > 0.5 * probe { - kappa = kappa.max(crate::vec3::dot(na, nx).clamp(-1.0, 1.0).acos() / h); + kappa = kappa.max(libm::acos(crate::vec3::dot(na, nx).clamp(-1.0, 1.0)) / h); } } } @@ -616,7 +616,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o .zip(&class_max) .zip(&class_arc) .map(|((&(sum, k), &longest), &arc)| { - let mean = if k == 0 { 1.0 } else { (sum / k as f64).exp().round() }; + let mean = if k == 0 { 1.0 } else { libm::exp(sum / k as f64).round() }; (mean.max((longest / MAX_STRETCH).ceil()) as usize).max(arc).clamp(1, crate::block::MAX_DIVISIONS) }) .collect(); diff --git a/crates/mesh/src/relax/strong.rs b/crates/mesh/src/relax/strong.rs index 316fe11..8cf401f 100644 --- a/crates/mesh/src/relax/strong.rs +++ b/crates/mesh/src/relax/strong.rs @@ -65,7 +65,7 @@ enum Phase { /// The objective's term for one element of quality `q` whose size went from `v0` to `v`. fn phi(phase: Phase, q: f64, valid: bool, v: f64, v0: f64) -> f64 { if valid && q > 0.0 { - let g = if v > 0.0 && v0 > 0.0 { (v / v0).ln() } else { 0.0 }; + let g = if v > 0.0 && v0 > 0.0 { libm::log(v / v0) } else { 0.0 }; let shape = match phase { Phase::Smooth => (1.0 / q - 1.0).powi(2), Phase::MinMax => q.powi(-MIN_MAX_POWER) - 1.0, diff --git a/crates/mesh/src/surface/topology.rs b/crates/mesh/src/surface/topology.rs index a8447cf..576322b 100644 --- a/crates/mesh/src/surface/topology.rs +++ b/crates/mesh/src/surface/topology.rs @@ -236,7 +236,7 @@ fn solid_angle(s: &Surface, t: usize, x: [f64; 3]) -> f64 { let [a, b, c] = s.tri_xyz(t).map(|p| sub(p, x)); let (la, lb, lc) = (norm(a), norm(b), norm(c)); let den = la * lb * lc + dot(a, b) * lc + dot(a, c) * lb + dot(b, c) * la; - 2.0 * det(a, b, c).atan2(den) + 2.0 * crate::vec3::atan2(det(a, b, c), den) } /// Edges between different B-rep faces, and edges inside a face whose normals turn by more than