From 01830e5c8d96f7e23b89eb31664ef37a7fcf7447 Mon Sep 17 00:00:00 2001 From: Omar Sobh Date: Tue, 15 Sep 2026 23:08:09 -0500 Subject: [PATCH] =?UTF-8?q?PERF-2=20P1:=20exact=20chunked=20polyline=20ind?= =?UTF-8?q?ex=20for=20the=20O-grid's=20sliding=20outer=20ring=20=E2=80=94?= =?UTF-8?q?=20PolylineIndex::nearest=20returns=20nearest=5Fon=5Fpolyline's?= =?UTF-8?q?=20point=20bit=20for=20bit=20(same=20segments,=20same=20order,?= =?UTF-8?q?=20same=20strict=20minimum;=20chunks=20skipped=20only=20when=20?= =?UTF-8?q?their=20box=20is=20farther=20than=20the=20best=20by=20more=20th?= =?UTF-8?q?an=201e-12=20relative),=20used=20by=20the=20Winslow=20sweep=20a?= =?UTF-8?q?nd=20the=20two=20ring=20projections;=20pin:=200=20mismatches=20?= =?UTF-8?q?over=2039,864=20queries=20on=20a=203,166-point=20rounded=20outl?= =?UTF-8?q?ine=20including=20vertices,=20mid-points=20and=20ties?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01YJPeT6WA2e7YvAnS875AHL --- .../specialized/rtx-cfd/src/mesh/patch_gen.rs | 97 +++++++++++++++++-- .../rtx-cfd/tests/polyline_index_exact.rs | 89 +++++++++++++++++ 2 files changed, 176 insertions(+), 10 deletions(-) create mode 100644 crates/specialized/rtx-cfd/tests/polyline_index_exact.rs diff --git a/crates/specialized/rtx-cfd/src/mesh/patch_gen.rs b/crates/specialized/rtx-cfd/src/mesh/patch_gen.rs index 293b4c9..162a18f 100644 --- a/crates/specialized/rtx-cfd/src/mesh/patch_gen.rs +++ b/crates/specialized/rtx-cfd/src/mesh/patch_gen.rs @@ -280,7 +280,88 @@ fn point_at_fraction(pts: &[[f64; 2]], cum: &[f64], f: f64) -> [f64; 2] { } /// Nearest point on the closed polyline `pts` to `q` (segment projection). -fn nearest_on_polyline(pts: &[[f64; 2]], q: [f64; 2]) -> [f64; 2] { +/// Chunked bounding boxes over a closed polyline for EXACT nearest-point +/// queries (PERF-2 P1, `docs/perf2_campaign.md`): [`Self::nearest`] returns +/// the same point as [`nearest_on_polyline`] — the same segments are tested +/// in the same order with the same strict-minimum rule — but skips every +/// chunk of segments whose box is farther from the query than the best +/// distance so far (with a 1e-12 relative margin, so a box's rounded +/// lower bound can never hide a segment that could still win). +pub struct PolylineIndex<'a> { + pts: &'a [[f64; 2]], + chunk: usize, + /// Per chunk: `[xmin, xmax, ymin, ymax]` of its segments' end points. + boxes: Vec<[f64; 4]>, +} + +impl<'a> PolylineIndex<'a> { + /// Index `pts` (a closed polyline: segment `i` joins `pts[i]` and + /// `pts[(i + 1) % n]`) in chunks of 32 segments. + pub fn new(pts: &'a [[f64; 2]]) -> Self { + let n = pts.len(); + let chunk = 32; + let boxes = (0..n.div_ceil(chunk)) + .map(|c| { + let mut b = [ + f64::INFINITY, + f64::NEG_INFINITY, + f64::INFINITY, + f64::NEG_INFINITY, + ]; + let start = c * chunk; + let end = (start + chunk).min(n); + // The chunk's segments' end points: `start..=end` (wrapping). + for i in start..=end { + let p = pts[i % n]; + b[0] = b[0].min(p[0]); + b[1] = b[1].max(p[0]); + b[2] = b[2].min(p[1]); + b[3] = b[3].max(p[1]); + } + b + }) + .collect(); + Self { pts, chunk, boxes } + } + + /// The nearest point of the polyline to `q`, bit for bit the answer of + /// [`nearest_on_polyline`]. + pub fn nearest(&self, q: [f64; 2]) -> [f64; 2] { + let pts = self.pts; + let n = pts.len(); + let (mut best, mut best_d) = (pts[0], f64::INFINITY); + for (c, b) in self.boxes.iter().enumerate() { + let dx = (b[0] - q[0]).max(0.0).max(q[0] - b[1]); + let dy = (b[2] - q[1]).max(0.0).max(q[1] - b[3]); + if dx * dx + dy * dy > best_d * (1.0 + 1e-12) { + continue; + } + let start = c * self.chunk; + let end = (start + self.chunk).min(n); + for i in start..end { + let (a, b) = (pts[i], pts[(i + 1) % n]); + let (dx, dy) = (b[0] - a[0], b[1] - a[1]); + let l2 = dx * dx + dy * dy; + let t = if l2 > 0.0 { + (((q[0] - a[0]) * dx + (q[1] - a[1]) * dy) / l2).clamp(0.0, 1.0) + } else { + 0.0 + }; + let p = [a[0] + t * dx, a[1] + t * dy]; + let d = (p[0] - q[0]).powi(2) + (p[1] - q[1]).powi(2); + if d < best_d { + best_d = d; + best = p; + } + } + } + best + } +} + +/// The nearest point of the closed polyline `pts` to `q` by a scan of every +/// segment (the reference for [`PolylineIndex`]). +pub fn nearest_on_polyline(pts: &[[f64; 2]], q: [f64; 2]) -> [f64; 2] { let n = pts.len(); let (mut best, mut best_d) = (pts[0], f64::INFINITY); for i in 0..n { @@ -489,10 +570,11 @@ pub fn winslow_smooth( sweeps += 1; moved = 0.0; if let Some(curve) = outer_curve { + let index = PolylineIndex::new(curve); for i in 0..ns { let (a, b) = (x[nn - 2][i], x[nn - 1][i]); let cand = [2.0 * b[0] - a[0], 2.0 * b[1] - a[1]]; - let p = nearest_on_polyline(curve, cand); + let p = index.nearest(cand); let d = ((p[0] - x[nn][i][0]).powi(2) + (p[1] - x[nn][i][1]).powi(2)).sqrt(); moved = moved.max(d); x[nn][i] = p; @@ -1109,6 +1191,7 @@ fn o_grid_from_outline( let t0 = std::time::Instant::now(); let hull = convex_hull(hull_src); let outer_poly = offset_convex_polygon(&hull, offset, 24); + let outer_index = PolylineIndex::new(&outer_poly); regen_charge(1, t0); let t0 = std::time::Instant::now(); // Initial outer ring: the inner point pushed along its outward normal @@ -1125,10 +1208,7 @@ fn o_grid_from_outline( let (dx, dy) = (next[0] - prev[0], next[1] - prev[1]); let l = (dx * dx + dy * dy).sqrt().max(1e-300); let normal = [-dy / l, dx / l]; - nearest_on_polyline( - &outer_poly, - [cur[0] + offset * normal[0], cur[1] + offset * normal[1]], - ) + outer_index.nearest([cur[0] + offset * normal[0], cur[1] + offset * normal[1]]) }) .collect(); regen_charge(2, t0); @@ -1146,10 +1226,7 @@ fn o_grid_from_outline( // rows' converged distribution (the crude normal push would // undo it every build — measured as the interior moving // 0.05 h with the wall at rest). - rows[nn] = rows[nn] - .iter() - .map(|&p| nearest_on_polyline(&outer_poly, p)) - .collect(); + rows[nn] = rows[nn].iter().map(|&p| outer_index.nearest(p)).collect(); rows } None => eta diff --git a/crates/specialized/rtx-cfd/tests/polyline_index_exact.rs b/crates/specialized/rtx-cfd/tests/polyline_index_exact.rs new file mode 100644 index 0000000..86675f8 --- /dev/null +++ b/crates/specialized/rtx-cfd/tests/polyline_index_exact.rs @@ -0,0 +1,89 @@ +//! PERF-2 P1 (`docs/perf2_campaign.md`): the chunked polyline index used +//! by the O-grid's sliding outer ring must return the brute-force scan's +//! point BIT FOR BIT — the recorded marches are gated on digit identity. +//! The polygon is the shape the ring sees: a convex hull offset sampled +//! 24 points per corner (a rounded cylinder + flag outline, ~3,000 +//! points); the queries sit inside, on and outside it, and on a lattice +//! whose points fall exactly on vertices and segment mid-points (ties). + +use rtx_cfd::mesh::patch_gen::{PolylineIndex, nearest_on_polyline}; +use std::f64::consts::PI; + +/// A rounded "cylinder + flag" convex outline: a circle of radius 0.05 at +/// the origin joined to a rectangle to x = 0.4, corners rounded with +/// radius 0.06 sampled 24 per corner, plus fine sampling on the arcs. +fn outline() -> Vec<[f64; 2]> { + let mut pts = Vec::new(); + // Left semicircle (the cylinder side), 256 points. + for m in 0..256 { + let th = PI / 2.0 + PI * m as f64 / 255.0; + pts.push([0.11 * th.cos(), 0.11 * th.sin()]); + } + // Bottom edge to the tip, 1,200 points. + for m in 1..1200 { + pts.push([-0.0 + 0.46 * m as f64 / 1200.0, -0.11]); + } + // Right rounded end, 512 points. + for m in 0..512 { + let th = -PI / 2.0 + PI * m as f64 / 511.0; + pts.push([0.46 + 0.11 * th.cos(), 0.11 * th.sin()]); + } + // Top edge back, 1,200 points. + for m in 1..1200 { + pts.push([0.46 - 0.46 * m as f64 / 1200.0, 0.11]); + } + pts +} + +#[test] +fn indexed_nearest_point_is_the_brute_force_point_bit_for_bit() { + let pts = outline(); + let index = PolylineIndex::new(&pts); + let mut queries: Vec<[f64; 2]> = Vec::new(); + // A lattice through the whole region (inside, on, outside). + for i in 0..120 { + for j in 0..60 { + queries.push([-0.2 + 0.8 * i as f64 / 119.0, -0.25 + 0.5 * j as f64 / 59.0]); + } + } + // Every vertex, every segment mid-point, and points pushed off them. + let n = pts.len(); + for i in 0..n { + let (a, b) = (pts[i], pts[(i + 1) % n]); + queries.push(a); + queries.push([0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])]); + queries.push([a[0] * 1.3, a[1] * 1.3]); + queries.push([a[0] * 0.7 + 0.01, a[1] * 0.7]); + } + // A deterministic pseudo-random cloud. + let mut state = 0x9e37_79b9_7f4a_7c15_u64; + for _ in 0..20_000 { + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let u = (state >> 11) as f64 / (1u64 << 53) as f64; + state ^= state << 13; + state ^= state >> 7; + state ^= state << 17; + let v = (state >> 11) as f64 / (1u64 << 53) as f64; + queries.push([-0.3 + 1.0 * u, -0.3 + 0.6 * v]); + } + let mut mismatches = 0; + for &q in &queries { + let a = nearest_on_polyline(&pts, q); + let b = index.nearest(q); + if a[0].to_bits() != b[0].to_bits() || a[1].to_bits() != b[1].to_bits() { + mismatches += 1; + if mismatches < 5 { + eprintln!("mismatch at q = {q:?}: brute {a:?} vs index {b:?}"); + } + } + } + println!( + " polyline index: {} points, {} chunks tested over {} queries, {mismatches} mismatches", + pts.len(), + pts.len().div_ceil(32), + queries.len() + ); + assert_eq!(mismatches, 0); +}