brep: a face's points search one grid of its surface, made once

uv::project, given no hint that converges within `accept`, evaluated a grid
of up to 16384 points of the whole surface for every point it placed, then
sorted all of them by distance. The grid is the same for every point of a
surface, and a face's pcurves place hundreds of points on one surface: on a
modelled thread's helical faces that grid was most of a boolean's time on
SMLib, in the operands' records, SMLib's export of each body and the
result's (threads/2526, the ladder's threads).

uv::Starts keeps the grid, filled on first use; uv::project_with takes it,
and a face's pcurve fit (pcurve::On) keeps one per face. The 32 nearest
starts are picked with a partial select instead of a full sort, in the same
order (nearest first, grid order where equally near), so every point gets
the parameters it got before. (The tessellation's boundary projection can
take the same grid; crates/kernel/src/tess is the mesh agent's.)

Not built: base in progress (rebased onto kernel-smlib/no-truck 41adb372).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt

crates/geom/src/brep/eval.rs

@@ -179,8 +179,14 @@ impl SurfaceEval {
     /// The parameters of `p`, a point within `accept` of the surface, from `near` when given
     /// (not unwrapped to it).
     pub fn project(&self, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
+        self.project_with(&uv::Starts::default(), p, near, accept)
+    }
+
+    /// [`SurfaceEval::project`] with the search grid of a spline or revolution kept in `starts`
+    /// ([`uv::project_with`]), which must serve this surface only.
+    pub fn project_with(&self, starts: &uv::Starts, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
         let q = match self {
-            SurfaceEval::Eval(e, dom) => return uv::project(e, dom, p, near, accept),
+            SurfaceEval::Eval(e, dom) => return uv::project_with(e, dom, starts, p, near, accept),
             SurfaceEval::Cylinder { axes, .. } => {
                 let [a, b, h] = axes.local(p);
                 [angle(a, b, near.map(|n| n[0])), h]

crates/geom/src/brep/pcurve.rs

@@ -148,9 +148,11 @@ fn own_pcurve(c: &Curve2, [t0, t1]: [f64; 2], budget: &mut Budget) -> Result<Cur
     Ok(c.clone())
 }
 
-/// A face's surface: its evaluator and parameter box.
+/// A face's surface: its evaluator and parameter box, and the grid a point's parameters on it
+/// are searched from, made once for all the face's points.
 struct On {
     surface: SurfaceEval,
+    starts: uv::Starts,
     dom: Domain,
     /// The face's own parameter box, where its kernel gives one.
     face: Option<[[f64; 2]; 2]>,
@@ -172,7 +174,7 @@ impl On {
     fn new(s: &Surface) -> Result<Self> {
         let surface = SurfaceEval::new(s)?;
         let dom = surface.domain();
-        Ok(On { surface, dom, face: None })
+        Ok(On { surface, starts: uv::Starts::default(), dom, face: None })
     }
 
     fn within(self, face: Option<[[f64; 2]; 2]>) -> Self {
@@ -293,7 +295,7 @@ impl On {
 
     /// The parameters of `p` within `accept`, from `near` (and unwrapped to it) when given.
     fn project(&self, p: Vec3, near: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
-        let mut q = self.surface.project(p, near, accept)?;
+        let mut q = self.surface.project_with(&self.starts, p, near, accept)?;
         if let Some(n) = near {
             for (d, x) in q.iter_mut().enumerate() {
                 *x = uv::unwrap(*x, n[d], self.dom.period.get(d).copied().flatten());

crates/geom/src/tess/uv.rs

@@ -169,10 +169,57 @@ fn newton(s: &EvalSurface, dom: &Domain, p: Vec3, start: [f64; 2]) -> Option<([f
     Some((uv, dist))
 }
 
+/// The grid of surface points [`project`] searches from when it has no hint that converges: up to
+/// [`MAX_STARTS`] evaluations of the surface, the same for every point, so a caller projecting
+/// many points onto one surface (a face's pcurves, its boundary) keeps one `Starts` for it and
+/// the grid is evaluated once, on first use ([`project_with`]).
+#[derive(Clone, Debug, Default)]
+pub struct Starts(std::sync::OnceLock<Vec<(Vec3, [f64; 2])>>);
+
+/// A cache of the surface's own points: any two are alike.
+impl PartialEq for Starts {
+    fn eq(&self, _: &Starts) -> bool {
+        true
+    }
+}
+
+impl Starts {
+    /// The grid over `dom` ([`starts_along`] in each direction, thinned to [`MAX_STARTS`]):
+    /// each start's point and parameters. Empty where a direction is unbounded.
+    fn of(&self, s: &EvalSurface, dom: &Domain) -> &[(Vec3, [f64; 2])] {
+        self.0.get_or_init(|| {
+            let (Some(ru), Some(rv)) = (dom.range[0], dom.range[1]) else { return Vec::new() };
+            let (mut us, mut vs) = (starts_along(s, 0, ru), starts_along(s, 1, rv));
+            // Thin the longer list (every other start, its ends kept) until the grid fits.
+            while us.len().saturating_mul(vs.len()) > MAX_STARTS {
+                let longer = if us.len() >= vs.len() { &mut us } else { &mut vs };
+                let last = longer.last().copied();
+                *longer = longer.iter().step_by(2).copied().chain(last).collect();
+                longer.dedup();
+            }
+            let mut grid = Vec::with_capacity(us.len() * vs.len());
+            for &u in &us {
+                for &v in &vs {
+                    if let Some((q, _, _)) = s.derivs(u, v) {
+                        grid.push((q, [u, v]));
+                    }
+                }
+            }
+            grid
+        })
+    }
+}
+
 /// Parameters of `p` (a point on the surface within `accept`), from `hint` when that
 /// converges, else from the best starts of a grid over the surface's domain, as fine as its
 /// knot spans ([`starts_along`]).
 pub fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
+    project_with(s, dom, &Starts::default(), p, hint, accept)
+}
+
+/// [`project`] with the grid kept in `starts`, which must have been used with this surface and
+/// domain only.
+pub fn project_with(s: &EvalSurface, dom: &Domain, starts: &Starts, p: Vec3, hint: Option<[f64; 2]>, accept: f64) -> Option<[f64; 2]> {
     if let EvalSurface::Plane { origin, du, dv } = s {
         let r = p - *origin;
         let (a, b, c) = (du.len2(), du.dot(*dv), dv.len2());
@@ -193,27 +240,21 @@ pub fn project(s: &EvalSurface, dom: &Domain, p: Vec3, hint: Option<[f64; 2]>, a
         }
         best = Some((uv, d));
     }
-    let (mut us, mut vs) = (starts_along(s, 0, dom.range[0]?), starts_along(s, 1, dom.range[1]?));
-    // Thin the longer list (every other start, its ends kept) until the grid fits.
-    while us.len().saturating_mul(vs.len()) > MAX_STARTS {
-        let longer = if us.len() >= vs.len() { &mut us } else { &mut vs };
-        let last = longer.last().copied();
-        *longer = longer.iter().step_by(2).copied().chain(last).collect();
-        longer.dedup();
+    if dom.range.iter().any(Option::is_none) {
+        return None;
     }
-    let mut starts: Vec<(f64, [f64; 2])> = Vec::with_capacity(us.len() * vs.len());
-    for &u in &us {
-        for &v in &vs {
-            if let Some((q, _, _)) = s.derivs(u, v) {
-                starts.push(((q - p).len2(), [u, v]));
-            }
-        }
+    let grid = starts.of(s, dom);
+    // The nearest starts first (in grid order where equally near), until one converges within
+    // `accept`: where turns of a surface pass close to each other (a thread's root band a tenth
+    // of a millimetre from the next turn's), the nearest few can all lie on the wrong turn.
+    let mut near: Vec<(f64, usize)> = grid.iter().enumerate().map(|(i, (q, _))| ((*q - p).len2(), i)).collect();
+    let by_distance = |a: &(f64, usize), b: &(f64, usize)| a.0.total_cmp(&b.0).then(a.1.cmp(&b.1));
+    if near.len() > NEWTON_STARTS {
+        near.select_nth_unstable_by(NEWTON_STARTS, by_distance);
+        near.truncate(NEWTON_STARTS);
     }
-    starts.sort_by(|a, b| a.0.total_cmp(&b.0));
-    // The nearest starts first, until one converges within `accept`: where turns of a surface
-    // pass close to each other (a thread's root band a tenth of a millimetre from the next
-    // turn's), the nearest few can all lie on the wrong turn.
-    for (_, start) in starts.iter().take(NEWTON_STARTS) {
+    near.sort_by(by_distance);
+    for start in near.iter().filter_map(|(_, i)| grid.get(*i).map(|(_, uv)| uv)) {
         if let Some((uv, d)) = newton(s, dom, p, *start)
             && best.is_none_or(|(_, bd)| d < bd)
         {