diff --git a/crates/mesh/src/hex/decompose.rs b/crates/mesh/src/hex/decompose.rs index ee29dea..44f023b 100644 --- a/crates/mesh/src/hex/decompose.rs +++ b/crates/mesh/src/hex/decompose.rs @@ -916,8 +916,8 @@ fn assignable(loops: &[Poly], cuts: &[Vec], s: &Settings) -> bool { fn o_grid(loops: &[Poly], s: &Settings, budget: &mut usize) -> Vec> { let (Some(outer), Some(hole)) = (loops.first(), loops.get(1)) else { return Vec::new() }; let c = centroid(&hole.pts); - let corner = |p: &Poly, i: usize| p.is_end(i) && p.angle(i) <= 0.75 * std::f64::consts::PI; - let concave = |p: &Poly, i: usize| p.is_end(i) && p.angle(i) >= 1.25 * std::f64::consts::PI; + let corner = |p: &Poly, i: usize| p.is_end(i) && region::is_corner(p.angle(i)); + let concave = |p: &Poly, i: usize| p.is_end(i) && region::is_reflex(p.angle(i)); // A cut from loop point (l, i) along d that must end on loop `target`. let ray = |l: usize, i: usize, d: Q, target: usize| -> Option> { let o = loops.get(l)?.at(i); @@ -1134,8 +1134,9 @@ mod tests { let n = ((a1 - a0).abs() / 5.0).ceil().max(1.0) as usize; (0..=n) .map(|k| { + // The mesher's portable trig, so every platform gets the same points (#556). let t = (a0 + (a1 - a0) * k as f64 / n as f64).to_radians(); - [c[0] + r * t.cos(), c[1] + r * t.sin()] + [c[0] + r * libm::cos(t), c[1] + r * libm::sin(t)] }) .collect() } @@ -1308,6 +1309,10 @@ mod tests { ]); let m = mesh_profile(&profile, &[], 2.5); check(&m, &profile, 0.5); + // The intended O-grid on every platform (#556): Windows' last-bit-different arcs used to + // give 10 regions and 2063 quads with a 0.375 corner. + assert_eq!(m.section.schemes.len(), 9, "{:?}", m.section.schemes); + assert_eq!(m.section.quads.iter().map(Vec::len).sum::(), 848); } #[test] diff --git a/crates/mesh/src/hex/region.rs b/crates/mesh/src/hex/region.rs index c84d611..6e6af54 100644 --- a/crates/mesh/src/hex/region.rs +++ b/crates/mesh/src/hex/region.rs @@ -47,6 +47,24 @@ const MAX_ANNULUS_NODES: usize = 4096; /// Smoothing sweeps over the interior of a triangle primitive. const TRIANGLE_SWEEPS: usize = 20; +/// How close (radians) an interior angle must come to a classification threshold to count as on +/// it. Exact geometry lands on the thresholds (an O-grid's 45° cuts meet straight sides at +/// exactly 135°), and the computed angle is then rounded either way by the last bits of its +/// points: those differ between platforms' `sin`/`cos` and between tessellations, and the +/// scheme must not (#556). An angle within this of a threshold is taken as exactly on it, and +/// on a threshold is smooth, neither a corner nor reflex: that keeps the slot's intended O-grid. +const ANGLE_TOLERANCE: f64 = 1e-9; + +/// A corner: an interior angle clearly below 135° (exactly 135° is not one). +pub(crate) fn is_corner(angle: f64) -> bool { + angle < 0.75 * PI - ANGLE_TOLERANCE +} + +/// A reflex turn no mapped or triangle region can take: clearly above 225°. +pub(crate) fn is_reflex(angle: f64) -> bool { + angle > 1.25 * PI + ANGLE_TOLERANCE +} + #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub(crate) enum Scheme { /// Corner positions (curve uses) in the single loop. @@ -89,10 +107,10 @@ pub(crate) fn classify(loops: &[Vec>]) -> Scheme { match loops { [lp] => { let Some(a) = angles(lp, ccw(lp)) else { return Scheme::Free }; - if a.iter().any(|&x| x >= 1.25 * PI) { + if a.iter().any(|&x| is_reflex(x)) { return Scheme::Free; } - let corners: Vec = (0..a.len()).filter(|&k| a.get(k).is_some_and(|&x| x <= 0.75 * PI)).collect(); + let corners: Vec = (0..a.len()).filter(|&k| a.get(k).is_some_and(|&x| is_corner(x))).collect(); if corners.is_empty() { return Scheme::Disk; } @@ -103,7 +121,7 @@ pub(crate) fn classify(loops: &[Vec>]) -> Scheme { } } [outer, inner] => { - let smooth = |lp: &Vec>| angles(lp, ccw(lp)).is_some_and(|a| a.iter().all(|&x| x > 0.75 * PI && x < 1.25 * PI)); + let smooth = |lp: &Vec>| angles(lp, ccw(lp)).is_some_and(|a| a.iter().all(|&x| !is_corner(x) && !is_reflex(x))); if smooth(outer) && smooth(inner) { Scheme::Annulus } else { Scheme::Free } } _ => Scheme::Free, @@ -568,6 +586,30 @@ mod tests { use super::*; use crate::size::UniformSize; + /// A corner of exactly 135°, the threshold, is not a corner however the last bits of its + /// points round (Windows' `sin` against glibc's, or two tessellations): a trapezoid with one + /// stays a triangle rather than flipping to mapped when that vertex moves by an ulp (#556). + #[test] + fn a_corner_on_the_threshold_is_one_whichever_way_it_rounds() { + // (0, 0) 90°, (3, 0) 45°, (2, 1) 135°, (0, 1) 90°: three corners and a 135° bend. + let trapezoid = |p: Q| { + let pts = [[0.0, 0.0], [3.0, 0.0], p, [0.0, 1.0]]; + let curves: Vec> = (0..4).map(|k| vec![pts[k], pts[(k + 1) % 4]]).collect(); + classify(&[curves]) + }; + let steps: [fn(f64) -> f64; 3] = [|v| v, f64::next_up, f64::next_down]; + for nx in steps { + for ny in steps { + let mut p = [nx(2.0), ny(1.0)]; + // A few ulps further each way as well. + for _ in 0..3 { + assert_eq!(trapezoid(p), Scheme::Triangle([0, 1, 3]), "{p:?}"); + p = [nx(p[0]), ny(p[1])]; + } + } + } + } + fn circle(r: f64, n: usize, ccw: bool) -> Vec { (0..n) .map(|k| {