surface_p10: a leaf whose surface cannot bound a part on its own is not halved

SurfaceCache::boxes_over gave a surface carried by no B-spline net and not a plane
(an offset surface, a procedural blend, a curve-bounded or STEP surface) the whole
surface's box and normal map for every part. Halving such a leaf dropped no part,
narrowed no interval and cleared no perpendicular, so parts_crossed_once went the
full 20 levels: 2^20 parts for one leaf of the shelled wall in
issue_1795_extended_edges_mated (1590 s, 0.4 s before the part solve), and
issue_1795_shell_hole_through_arc took 3555 s (0.6 s).

boxes_over is now None for such a surface, and its leaf is one part, solved as
upstream solves it. A leaf is also split into at most MAX_PARTS (256) parts.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

solidcore/crates/goldens/expected/surface.surface_p10.jsonl

@@ -68,7 +68,7 @@
 {"case":"gen0019","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":{"blocks":[[29.04339384379934,37.04339384379933,561.842978165522]],"n":47}}}
 {"case":"gen0020","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":{"blocks":[[8.931789963706066,12.931789963706066,38.27314715573294]],"n":17}}}
 {"case":"gen0021","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":[0.0,0.0,0.0]}}
-{"case":"gen0022","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":[0.0,0.0,0.0]}}
+{"case":"gen0022","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":{"blocks":[[11.831537595960805,15.831537595960805,67.68312101060808]],"n":18}}}
 {"case":"gen0023","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":[0.0,0.0,0.0]}}
 {"case":"gen0024","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":{"blocks":[[9.936154451873268,13.936154451873268,50.43372848535808]],"n":17}}}
 {"case":"gen0025","deferred":null,"error":null,"ints":{},"topology":null,"values":{"rec":{"blocks":[[9.925515262295008,13.925515262295008,49.551620172494786]],"n":17}}}

solidcore/crates/solidcore-ops/tests/classify_point_grazing_ray.rs

@@ -0,0 +1,79 @@
+//! A point just inside a solid cylinder, near the top of its circle, classifies as inside.
+//! The classification fires a ray along +x first; from such a point the ray leaves through the
+//! side face at a shallow angle (a few degrees), and that crossing must be found.
+//! SolidCraft #12539: an R1 fillet beside a 1 mm chamfer round an R30.475 outline was refused
+//! as cutting the part, because the middle of the round's section, 0.08 mm inside the wall,
+//! classified as outside.
+#![allow(clippy::unwrap_used, clippy::expect_used, clippy::panic)] // test
+
+mod common;
+
+use common::built;
+use solidcore::math::Point3;
+use solidcore::primitive::api_primitives::create_cylinder;
+use solidcore_ops::{Body, PointClass};
+use solidcore_record::brep::Brep as Record;
+
+const R: f64 = 30.475;
+const AXIS: [f64; 2] = [-8.0, 0.0];
+
+/// The point at `angle` degrees round the axis, `radial` from it, at height `z`.
+#[allow(clippy::disallowed_methods)] // the test's sample points, not ops computing geometry
+fn at(angle: f64, radial: f64, z: f64) -> [f64; 3] {
+    let a = angle.to_radians();
+    [AXIS[0] + radial * a.cos(), AXIS[1] + radial * a.sin(), z]
+}
+
+/// The points either side of the wall round the top of the circle that classify wrongly.
+fn misclassified(cyl: &Body) -> Vec<String> {
+    let mut wrong = Vec::new();
+    // Round the top of the circle, where a +x ray meets the wall at a shallow angle, and below
+    // it, where it meets the wall steeply. 0.075 mm either side of the wall: inside exactly when
+    // nearer the axis than the radius.
+    for step in 0..=60 {
+        let angle = 60.0 + 0.5 * f64::from(step);
+        for (radial, want) in [
+            (R - 0.075, PointClass::Inside),
+            (R + 0.075, PointClass::Outside),
+        ] {
+            for z in [1.03, 10.0, 44.0] {
+                let p = at(angle, radial, z);
+                let got = cyl.classify_point(p).unwrap();
+                if got != want {
+                    wrong.push(format!(
+                        "{angle} deg, {radial} from the axis, z {z}: {got:?}, want {want:?}"
+                    ));
+                }
+            }
+        }
+    }
+    wrong
+}
+
+/// SolidCraft's cylinder: the extruded circle, its seam at 0.6131 rad from x, its rational
+/// spline in four quarter-circle spans.
+#[test]
+fn points_beside_an_extruded_circle_wall_classify_by_their_distance_to_its_axis() {
+    let record: Record =
+        serde_json::from_str(include_str!("data/extruded_circle_r30_475.json")).unwrap();
+    let cyl = Body::from_record(&record, true).unwrap().0;
+    let wrong = misclassified(&cyl);
+    assert!(
+        wrong.is_empty(),
+        "{} of 366 points misclassified:\n{}",
+        wrong.len(),
+        wrong.join("\n")
+    );
+}
+
+#[test]
+fn points_beside_a_cylinder_wall_classify_by_their_distance_to_its_axis() {
+    let cyl = built(|k| create_cylinder(k, Point3::new(AXIS[0], AXIS[1], 0.0), R, 45.0));
+    let wrong = misclassified(&cyl);
+    assert!(
+        wrong.is_empty(),
+        "{} of 366 points misclassified:\n{}",
+        wrong.len(),
+        wrong.join("\n")
+    );
+}

solidcore/crates/solidcore/src/surface/surface/p10.rs

@@ -44,8 +44,19 @@
 //!   runs once (nothing clears the flag) and is a single pass here.
 //! - **Debug only code** (`SM_DEBUG_CODE` drawing and dumps) and the assertions are not
 //!   ported.
-//! - **Bounds.** The leaf solve takes at most six guesses (upstream's fixed arrays). No
-//!   recursion.
+//! - **A leaf the line may cross twice** (not upstream's, WR-758, #3120). Upstream solves a
+//!   leaf from one set of guesses and, where its normal map has a direction
+//!   perpendicular to the line, looks for a second crossing only next to the first; a
+//!   first crossing behind a ray's start was dropped before that, so a line just under
+//!   a cylinder's top missed the crossing ahead. Such a leaf is solved part by part
+//!   ([`LsiGlobalSolver::parts_crossed_once`]): halved until each part's normal map has
+//!   no perpendicular to the line, so each holds at most one crossing, or until its
+//!   normals lie within the angular tolerance. A leaf without one is solved as upstream
+//!   solves it.
+//! - **Bounds.** Each leaf, or part of one, takes at most six guesses (upstream's fixed
+//!   arrays). A leaf is halved into at most [`MAX_PARTS`] parts, [`MAX_PART_DEPTH`]
+//!   deep, with an explicit stack; a leaf whose surface cannot bound a part on its own
+//!   (an offset surface) is not halved. No recursion.
 
 use super::p1::{is_node_valid, is_surface_cache_owner, surface_to_evaluate};
 use super::p2::{original_surface, LsiGlobalSolver, Surface, TsectCurveType};
@@ -61,11 +72,12 @@ use crate::geom::extent2d::Extent2d;
 use crate::geom::extent3d::LineHit as BoxHit;
 use crate::geom::geom_utility::p2::ContinuityType;
 use crate::geom::global_solver::{Branch, GlobalSolve, SolverOperation};
+use crate::geom::polar_box::PolarBox;
 use crate::geom::pseudo_box::LineHit as PseudoHit;
 use crate::geom::solution_array::{
     Solution, SolutionArray, SolutionEnd, SolutionObject, SolutionType, SortKey,
 };
-use crate::geom::tol::{are_same, scaled_zero2};
+use crate::geom::tol::{are_same, scaled_zero2, ANG_TOL_RAD};
 use crate::geom::tree::{AuxDataType, NodeId, Tree};
 use crate::geom::vector2d::{Point2, Vec2};
 use crate::geom::vector3d::{Point3, Vec3};
@@ -302,6 +314,77 @@ fn pseudo_hit(hit: &Result<PseudoHit>) -> Option<(u32, f64, f64)> {
 
 // ---- SmLSIGlobalSolver::LocalSolve -------------------------------------------------------------
 
+/// The most halvings [`LsiGlobalSolver::parts_crossed_once`] makes of a leaf: a part's
+/// side is then at least 2^-20 of the leaf's (about a millionth), below which two
+/// crossings of a line are one near-tangent crossing for the part's guesses.
+const MAX_PART_DEPTH: u32 = 20;
+
+/// The most parts [`LsiGlobalSolver::parts_crossed_once`] makes of a leaf: a part left to
+/// halve when the parts made and waiting reach it is kept as it is, crossed perhaps more
+/// than once, so a leaf's solve is bounded whatever its normal map does.
+const MAX_PARTS: usize = 256;
+
+/// A part of a leaf's domain the line may cross ([`LsiGlobalSolver::parts_crossed_once`]).
+struct LeafPart {
+    /// Its uv domain.
+    domain: Extent2d,
+    /// The line interval its box allows.
+    interval: Extent1d,
+    /// The line crosses it at most once (its normal map has no direction perpendicular to
+    /// the line).
+    once: bool,
+}
+
+/// The infinite line through `point` along `vector` (`sTmpLine`).
+fn line_curve(point: Point3, vector: Vec3) -> Curve {
+    Curve::from(Line::infinite(point, vector, 3, false))
+}
+
+/// The widest angle between two corners of `polar` (pi when it has no corners: unbounded or
+/// uninitialised).
+fn normal_spread(polar: &PolarBox) -> f64 {
+    let Ok(corners) = polar.boundary_vectors() else {
+        return std::f64::consts::PI;
+    };
+    let mut spread = 0.0_f64;
+    for (i, a) in corners.iter().enumerate() {
+        for b in corners.iter().skip(i + 1) {
+            spread = spread.max(
+                a.unitize_or_z()
+                    .dot(b.unitize_or_z())
+                    .clamp(-1.0, 1.0)
+                    .acos(),
+            );
+        }
+    }
+    spread
+}
+
+/// How far the unit normals of `surface` turn across `domain` in `u` (`along_u`) or `v`:
+/// the angles between them at the ends of the domain's two sides in that direction,
+/// summed. A corner whose normal cannot be evaluated (a pole) adds nothing.
+fn normal_turn(k: &Kernel, surface: &Surface, domain: Extent2d, along_u: bool) -> Result<f64> {
+    let (lo, hi) = (domain.min(), domain.max());
+    let normal = |u: f64, v: f64| -> Result<Option<Vec3>> {
+        Ok(
+            ignored(surface.evaluate_normal(k, Point2::new(u, v), Side::Left, Side::Left))?
+                .and_then(|n| n.unitize().ok()),
+        )
+    };
+    let sides = if along_u {
+        [((lo.x, lo.y), (hi.x, lo.y)), ((lo.x, hi.y), (hi.x, hi.y))]
+    } else {
+        [((lo.x, lo.y), (lo.x, hi.y)), ((hi.x, lo.y), (hi.x, hi.y))]
+    };
+    let mut turn = 0.0;
+    for ((u0, v0), (u1, v1)) in sides {
+        if let (Some(a), Some(b)) = (normal(u0, v0)?, normal(u1, v1)?) {
+            turn += a.dot(b).clamp(-1.0, 1.0).acos();
+        }
+    }
+    Ok(turn)
+}
+
 /// Read-only state of one leaf solve.
 struct LsiLeaf<'s, O> {
     cache: &'s SurfaceCache<O>,
@@ -338,6 +421,45 @@ impl<O: SolutionObject> LsiLeaf<'_, O> {
         solution.nodes[0] = self.nodes[0];
         solution.nodes[1] = self.nodes[1];
     }
+
+    /// This leaf over `domain`, part of its own, with the line interval `interval` and
+    /// the line `line`: a side the part shares with the leaf keeps the leaf's
+    /// continuity; a side inside the leaf is a discontinuity, so a crossing on it is kept
+    /// by the parts either side (and merged as a duplicate), never left to a neighbouring
+    /// branch that is not there.
+    fn part(&self, domain: Extent2d, interval: Extent1d, line: Curve) -> Self {
+        let [u_conts, v_conts, t_conts] = self.conts;
+        let side = |shared: bool, cont: ContinuityType| {
+            if shared {
+                cont
+            } else {
+                ContinuityType::Discontinuous
+            }
+        };
+        let (lo, hi) = (self.domain.min(), self.domain.max());
+        let (plo, phi) = (domain.min(), domain.max());
+        Self {
+            cache: self.cache,
+            owner: self.owner,
+            parent: self.parent,
+            eval: self.eval,
+            domain,
+            line,
+            interval,
+            conts: [
+                [
+                    side(plo.x == lo.x, u_conts[0]),
+                    side(phi.x == hi.x, u_conts[1]),
+                ],
+                [
+                    side(plo.y == lo.y, v_conts[0]),
+                    side(phi.y == hi.y, v_conts[1]),
+                ],
+                t_conts,
+            ],
+            nodes: self.nodes,
+        }
+    }
 }
 
 impl<'t, 's, O: SolutionObject + 't> LsiGlobalSolver<'t, 's, O> {
@@ -392,6 +514,254 @@ impl<'t, 's, O: SolutionObject + 't> LsiGlobalSolver<'t, 's, O> {
         Ok(keep)
     }
 
+    /// Solves `leaf` (a leaf, or a part of one) from its guesses (upstream's guess
+    /// points and their refinement by the surface's local curve intersection), adding the
+    /// crossings kept to `base.solutions`. `once`: the line crosses `leaf` at most once
+    /// (its normal map has no direction perpendicular to the line), so its first crossing
+    /// kept ends the solve, and true is returned; otherwise every guess is followed,
+    /// with the search for a near-tangent second crossing, and false is returned.
+    fn solve_part(&mut self, k: &Kernel, leaf: &LsiLeaf<'s, O>, once: bool) -> Result<bool> {
+        let domain = leaf.domain;
+        let interval = leaf.interval;
+        // let guess = LineInterval/SurfaceDomain CenterPoints
+        let mut first_t = interval.evaluate(0.5);
+        let first_uv_guess = domain.evaluate(0.5, 0.5);
+        let mut uv_guesses = Vec::new();
+        let mut t_guesses = Vec::new();
+
+        // Let SurfaceGuessPoint = Drop LinePoint down to SurfacePoint Tangent Plane
+        let drop = |uv: Point2, t: f64| -> Result<Point2> {
+            let point = self.line_point + t * self.line_vector;
+            super::p7::drop_points_to_tangent_plane(
+                k,
+                leaf.eval,
+                uv,
+                Side::Left,
+                Side::Left,
+                &[point],
+            )?
+            .first()
+            .copied()
+            .ok_or(SolidError::sm_err(OP_LSI))
+        };
+        let mut first_uv = drop(first_uv_guess, first_t)?;
+
+        // When the updated SurfaceGuessPoint is outside the leafNodeIntersectionDomain
+        let uv_eps = EFF_ZERO * (1.0 + domain.size().max_dimension());
+        if !domain.contains_point(first_uv, uv_eps) {
+            // move guess to SurfaceMidPoint and near LineInterval beginPoint.
+            first_t = interval.evaluate(0.1);
+            first_uv = drop(domain.evaluate(0.5, 0.5), first_t)?;
+
+            // When SurfaceMid/LineStart Drop is also outside the leafNodeIntersectionDomain
+            if !domain.contains_point(first_uv, 0.0) {
+                // move guess to SurfaceMidPoint and near LineInterval endPoint.
+                first_t = interval.evaluate(0.9);
+                first_uv = drop(domain.evaluate(0.5, 0.5), first_t)?;
+
+                // when this 3rd point is also outside the leafNodeIntersectionDomain
+                if !domain.contains_point(first_uv, 0.0) {
+                    // clamp the guess to the leafNodeIntersectionDomain
+                    first_uv = domain.clamp_point(first_uv);
+
+                    // add in 3 more guesses: SurfLeafCenterPoint/(Interval Start, Mid, End)
+                    for fraction in [0.5, 0.1, 0.9] {
+                        uv_guesses.push(domain.evaluate(0.5, 0.5));
+                        t_guesses.push(interval.evaluate(fraction));
+                    }
+                }
+            }
+        }
+        uv_guesses.insert(0, first_uv);
+        t_guesses.insert(0, first_t);
+
+        // for every Guess (more than 1 when LineGuess does not drop to SurfaceDomain)
+        let mut i = 0;
+        while let (Some(&uv_guess), Some(&t_guess)) = (uv_guesses.get(i), t_guesses.get(i)) {
+            let index = i;
+            i += 1;
+
+            // pass the call along to local curve intersection
+            let Some(hit) = surface_local_curve_intersect(
+                k,
+                leaf.eval,
+                &leaf.domain,
+                &leaf.line,
+                &leaf.interval,
+                self.base.tolerance_3d,
+                uv_guess,
+                t_guess,
+            )?
+            else {
+                continue;
+            };
+
+            // process intersections
+            let solution = leaf.solution(hit.t, hit.uv, hit.deviation)?;
+
+            // skip negative parameter ray/surface intersections
+            if self.ray_fire {
+                let scaled_zero = EFF_ZERO * (1.0 + self.line_point.max_dimension());
+                if hit.t < -scaled_zero {
+                    continue;
+                }
+            }
+
+            // keep exact solutions - check near solutions
+            if !filter_solution(
+                k,
+                leaf.eval,
+                &leaf.domain,
+                &leaf.line,
+                &leaf.interval,
+                leaf.conts,
+                &solution,
+            )? {
+                continue;
+            }
+
+            // for points within the surface's trim boundary
+            if self.is_point_valid(k, leaf.owner, leaf.cache, hit.uv)? {
+                self.add_sorted_solution(&solution, SortKey::ByAllParameters, true)?;
+            }
+
+            // Try to eliminate additional intersections with the patch's normal map.
+            if once {
+                return Ok(true);
+            }
+
+            // look for another intersection very close to this one (near tangent)
+            if index < 1 {
+                let class = classify_curve_surface_xsect_point(
+                    k,
+                    leaf.eval,
+                    &leaf.domain,
+                    &leaf.line,
+                    &leaf.interval,
+                    hit.uv,
+                    hit.t,
+                    None,
+                );
+                // status ignored upstream (SE)
+                if let Some(XsectPointClass {
+                    guess: Some((new_uv, new_t)),
+                    ..
+                }) = ignored(class)?
+                {
+                    if uv_guesses.len() >= MAX_GUESSES {
+                        return Err(SolidError::invalid_input(
+                            OP_LSI,
+                            "intersection guess limit exceeded",
+                        ));
+                    }
+                    uv_guesses.push(new_uv);
+                    t_guesses.push(new_t);
+                }
+            }
+        }
+
+        Ok(false)
+    }
+
+    /// The parts of `leaf` the line may cross, each crossed at most once where it can be
+    /// told: `leaf`'s normal map has a direction perpendicular to the line, so the line
+    /// may cross it more than once. Its domain is halved across the direction its
+    /// normals turn most in, until a part's normal map has no direction perpendicular to
+    /// the line (the normal-cone condition of Sederberg and Meyers, 1988: the line then
+    /// crosses the part at most once, `once`). A part whose box (grown by the tolerance)
+    /// the line misses within `leaf`'s interval, or only behind a ray's start, is
+    /// dropped; each part kept carries the line interval its box allows. A part that
+    /// still has a perpendicular after [`MAX_PART_DEPTH`] halvings is a near-tangent
+    /// spot, kept with `once` false, as is one whose normals cannot be bounded and one
+    /// whose normals lie within the angular tolerance ([`ANG_TOL_RAD`]): halving a flat
+    /// part does not narrow its normals, and a line lying in it is a coincidence, which
+    /// the leaf's own solve takes. A leaf whose surface cannot bound a part on its own
+    /// ([`SurfaceCache::boxes_over`] is `None`) is one part, solved as upstream solves it.
+    /// An explicit stack, at most [`MAX_PART_DEPTH`] deep and [`MAX_PARTS`] parts.
+    fn parts_crossed_once(
+        &self,
+        k: &Kernel,
+        leaf: &LsiLeaf<'s, O>,
+        line_vec: Vec3,
+        pass: bool,
+    ) -> Result<Vec<LeafPart>> {
+        let tol = self.base.tolerance_3d;
+        let mut parts = Vec::new();
+        let mut stack = vec![(leaf.domain, 0_u32)];
+        while let Some((domain, depth)) = stack.pop() {
+            // a surface with no bounds of a part's own (an offset surface, say) bounds
+            // every part as the whole leaf: halving drops nothing and narrows nothing, so
+            // the leaf is solved whole, as upstream solves it
+            let Some((bbox, polar)) = leaf.cache.boxes_over(k, &domain)? else {
+                parts.push(LeafPart {
+                    domain,
+                    interval: leaf.interval,
+                    once: false,
+                });
+                continue;
+            };
+            // the line interval within the part's box (the leaf's where it has no box)
+            let mut interval = leaf.interval;
+            if let Some(mut bbox) = bbox {
+                bbox.expand_absolute(tol);
+                let Some((enter, exit)) = box_hit(&passed_on(
+                    bbox.intersect_line(self.line_point, self.line_vector),
+                    pass,
+                )?) else {
+                    continue;
+                };
+                let mut hit = Extent1d::from_min_max(enter, exit);
+                hit.expand_absolute(2.0 * tol);
+                if hit.are_disjoint(&leaf.interval, 0.0) {
+                    continue;
+                }
+                interval = hit.intersect(&leaf.interval)?;
+            }
+            if self.ray_fire && interval.max() <= 0.0 {
+                continue;
+            }
+            // a part whose normals cannot be bounded is solved as it is: halving it tells
+            // nothing more. Nor does halving a part flat to the angular tolerance: its halves'
+            // normals are as close to the line's normal plane, and a line lying in it is a
+            // coincidence, which the leaf's own solve takes.
+            let once = polar.is_some_and(|b| !b.has_perpendicular_to_vector(line_vec));
+            let flat = polar.is_some_and(|b| normal_spread(&b) <= ANG_TOL_RAD);
+            if once
+                || flat
+                || polar.is_none()
+                || depth >= MAX_PART_DEPTH
+                || parts.len() + stack.len() + 2 > MAX_PARTS
+            {
+                parts.push(LeafPart {
+                    domain,
+                    interval,
+                    once,
+                });
+                continue;
+            }
+            let (lo, hi) = (domain.min(), domain.max());
+            let halves = if normal_turn(k, leaf.eval, domain, true)?
+                >= normal_turn(k, leaf.eval, domain, false)?
+            {
+                let mid = 0.5 * (lo.x + hi.x);
+                [
+                    Extent2d::from_bounds(lo.x, lo.y, mid, hi.y),
+                    Extent2d::from_bounds(mid, lo.y, hi.x, hi.y),
+                ]
+            } else {
+                let mid = 0.5 * (lo.y + hi.y);
+                [
+                    Extent2d::from_bounds(lo.x, lo.y, hi.x, mid),
+                    Extent2d::from_bounds(lo.x, mid, hi.x, hi.y),
+                ]
+            };
+            for half in halves.into_iter().rev() {
+                stack.push((half, depth + 1));
+            }
+        }
+        Ok(parts)
+    }
+
     /// `SmLSIGlobalSolver::LocalSolve`: intersects the line with the Bezier patch of
     /// the surface leaf node `branch[0]` (guess points refined by the surface's local
     /// curve intersection, then the line ends and the patch boundaries at
@@ -551,67 +921,14 @@ impl<'t, 's, O: SolutionObject + 't> LsiGlobalSolver<'t, 's, O> {
             None => bound_ivl,
         };
 
-        // let guess = LineInterval/SurfaceDomain CenterPoints
-        let mut first_t = ivl.evaluate(0.5);
-        let first_uv_guess = uv_domain.evaluate(0.5, 0.5);
-        let mut uv_guesses = Vec::new();
-        let mut t_guesses = Vec::new();
-
-        // Let SurfaceGuessPoint = Drop LinePoint down to SurfacePoint Tangent Plane
-        let drop = |uv: Point2, t: f64| -> Result<Point2> {
-            let point = self.line_point + t * self.line_vector;
-            super::p7::drop_points_to_tangent_plane(
-                k,
-                srf_eval,
-                uv,
-                Side::Left,
-                Side::Left,
-                &[point],
-            )?
-            .first()
-            .copied()
-            .ok_or(SolidError::sm_err(OP_LSI))
-        };
-        let mut first_uv = drop(first_uv_guess, first_t)?;
-
-        // When the updated SurfaceGuessPoint is outside the leafNodeIntersectionDomain
-        let uv_eps = EFF_ZERO * (1.0 + uv_domain.size().max_dimension());
-        if !uv_domain.contains_point(first_uv, uv_eps) {
-            // move guess to SurfaceMidPoint and near LineInterval beginPoint.
-            first_t = ivl.evaluate(0.1);
-            first_uv = drop(uv_domain.evaluate(0.5, 0.5), first_t)?;
-
-            // When SurfaceMid/LineStart Drop is also outside the leafNodeIntersectionDomain
-            if !uv_domain.contains_point(first_uv, 0.0) {
-                // move guess to SurfaceMidPoint and near LineInterval endPoint.
-                first_t = ivl.evaluate(0.9);
-                first_uv = drop(uv_domain.evaluate(0.5, 0.5), first_t)?;
-
-                // when this 3rd point is also outside the leafNodeIntersectionDomain
-                if !uv_domain.contains_point(first_uv, 0.0) {
-                    // clamp the guess to the leafNodeIntersectionDomain
-                    first_uv = uv_domain.clamp_point(first_uv);
-
-                    // add in 3 more guesses: SurfLeafCenterPoint/(Interval Start, Mid, End)
-                    for fraction in [0.5, 0.1, 0.9] {
-                        uv_guesses.push(uv_domain.evaluate(0.5, 0.5));
-                        t_guesses.push(ivl.evaluate(fraction));
-                    }
-                }
-            }
-        }
-        uv_guesses.insert(0, first_uv);
-        t_guesses.insert(0, first_t);
-
         // Make the line object (sTmpIvl = the bounded interval of interest).
-        let line = Curve::from(Line::infinite(self.line_point, self.line_vector, 3, false));
         let leaf = LsiLeaf {
             cache,
             owner: tree.owner(),
             parent,
             eval: srf_eval,
             domain: uv_domain,
-            line,
+            line: line_curve(self.line_point, self.line_vector),
             interval: ivl,
             conts: [u_conts, v_conts, t_conts],
             nodes: [
@@ -620,93 +937,27 @@ impl<'t, 's, O: SolutionObject + 't> LsiGlobalSolver<'t, 's, O> {
             ],
         };
 
-        // for every Guess (more than 1 when LineGuess does not drop to SurfaceDomain)
-        let mut i = 0;
-        while let (Some(&uv_guess), Some(&t_guess)) = (uv_guesses.get(i), t_guesses.get(i)) {
-            let index = i;
-            i += 1;
-
-            // pass the call along to local curve intersection
-            let Some(hit) = surface_local_curve_intersect(
-                k,
-                leaf.eval,
-                &leaf.domain,
-                &leaf.line,
-                &leaf.interval,
-                self.base.tolerance_3d,
-                uv_guess,
-                t_guess,
-            )?
-            else {
-                continue;
-            };
-
-            // process intersections
-            let solution = leaf.solution(hit.t, hit.uv, hit.deviation)?;
-
-            // skip negative parameter ray/surface intersections
-            if self.ray_fire {
-                let scaled_zero = EFF_ZERO * (1.0 + self.line_point.max_dimension());
-                if hit.t < -scaled_zero {
-                    continue;
-                }
-            }
-
-            // keep exact solutions - check near solutions
-            if !filter_solution(
-                k,
-                leaf.eval,
-                &leaf.domain,
-                &leaf.line,
-                &leaf.interval,
-                leaf.conts,
-                &solution,
-            )? {
-                continue;
-            }
-
-            // for points within the surface's trim boundary
-            if self.is_point_valid(k, leaf.owner, cache, hit.uv)? {
-                self.add_sorted_solution(&solution, SortKey::ByAllParameters, true)?;
-            }
-
-            // Try to eliminate additional intersections with the patch's normal map.
-            let line_vec = self.line_vector.unitize()?;
-            if !cache
-                .polar_box_unstored(k, surface_id)?
-                .has_perpendicular_to_vector(line_vec)
-            {
-                return Ok(false);
-            }
-
-            // look for another intersection very close to this one (near tangent)
-            if index < 1 {
-                let class = classify_curve_surface_xsect_point(
-                    k,
-                    leaf.eval,
-                    &leaf.domain,
-                    &leaf.line,
-                    &leaf.interval,
-                    hit.uv,
-                    hit.t,
-                    None,
+        // A leaf whose normal map has no direction perpendicular to the line is crossed
+        // at most once: its crossing ends the solve, as upstream's. One whose normal map
+        // has such a direction may be crossed twice (a line just under a cylinder's top
+        // crosses its leaf on either side of the top): it is solved part by part, each
+        // part crossed at most once, so a crossing behind a ray's start cannot hide one
+        // ahead of it (WR-758, #3120).
+        let line_vec = self.line_vector.unitize()?;
+        if cache
+            .polar_box_unstored(k, surface_id)?
+            .has_perpendicular_to_vector(line_vec)
+        {
+            for part in self.parts_crossed_once(k, &leaf, line_vec, pass)? {
+                let part_leaf = leaf.part(
+                    part.domain,
+                    part.interval,
+                    line_curve(self.line_point, self.line_vector),
                 );
-                // status ignored upstream (SE)
-                if let Some(XsectPointClass {
-                    guess: Some((new_uv, new_t)),
-                    ..
-                }) = ignored(class)?
-                {
-                    if uv_guesses.len() >= MAX_GUESSES {
-                        return Err(SolidError::invalid_input(
-                            OP_LSI,
-                            "intersection guess limit exceeded",
-                        ));
-                    }
-                    uv_guesses.push(new_uv);
-                    t_guesses.push(new_t);
-                }
+                self.solve_part(k, &part_leaf, part.once)?;
             }
+        } else if self.solve_part(k, &leaf, true)? {
+            return Ok(false);
         }
 
         // the extra work required at discontinuities of the line ends and patch boundaries

solidcore/crates/solidcore/src/surface/surface_cache/p1.rs

@@ -2617,6 +2617,62 @@ impl<O: Copy + PartialEq> SurfaceCache<O> {
         Ok(last_polar)
     }
 
+    /// The axis aligned box and the normal map (polar box) of the cache's surface over
+    /// `domain`, any part of its domain, not grown by a tolerance (not upstream's: the line solve splits a leaf into parts the line crosses at
+    /// most once, [`crate::surface::surface::p2::LsiGlobalSolver`]). A B-spline's own
+    /// boxes are its whole net's whatever the domain, so the box is the net's of the
+    /// sub-patch over `domain` (`CopySubPatch`), its normals the surface's own; a plane's
+    /// is its domain's. `None` for a surface carried by no B-spline net and not a plane
+    /// (an offset, procedural blend, curve-bounded or STEP surface): its bounds are the
+    /// whole surface's whatever the domain, so they tell nothing of one part. A box or
+    /// normal map the surface cannot bound is `None` in the pair.
+    ///
+    /// # Errors
+    ///
+    /// [`SolidError::Unported`] and cancellations from the bounding.
+    pub(crate) fn boxes_over(
+        &self,
+        k: &Kernel,
+        domain: &Extent2d,
+    ) -> Result<Option<(Option<crate::geom::extent3d::Extent3d>, Option<PolarBox>)>> {
+        let req = BoundingBoxRequest {
+            normal: true,
+            pseudo: false,
+            polar: true,
+            pseudo_guess: None,
+            polar_guess: None,
+            // (the caller grows the box by its own tolerance: a sub-patch has no face to
+            // read a zone tolerance from)
+            expand_by_zone_tol: false,
+            expand_tolerance: None,
+        };
+        let surface = self.surface();
+        let net = match surface.kind_with_nurb(k)? {
+            SurfaceKind::BSpline(b) => Some(b),
+            SurfaceKind::BallPatch(bp) => Some(bp.net()),
+            SurfaceKind::Cone(s) => Some(s.rev().base()),
+            SurfaceKind::Sphere(s) => Some(s.base().base()),
+            SurfaceKind::Torus(s) => Some(s.base().base()),
+            SurfaceKind::SurfOfRevolution(s) => Some(s.base()),
+            SurfaceKind::SurfOfExtrusion(s) => Some(s.base()),
+            SurfaceKind::Plane(_) => None,
+            _ => return Ok(None),
+        };
+        let boxes = match net {
+            Some(net) => match ignored(net.copy_sub_patch(domain))? {
+                Some(sub) => {
+                    ignored(sub.calculate_bounding_box_of(k, domain, &req, Some(surface)))?
+                }
+                None => None,
+            },
+            None => ignored(calculate_bounding_box(k, surface, domain, &req))?,
+        };
+        Ok(Some(match boxes {
+            Some(b) => (b.normal, b.polar),
+            None => (None, None),
+        }))
+    }
+
     /// `mBA_pSurface->CalculateBoundingBox(m_sUVDomain, NULL, NULL, &m_sPolarBox,
     /// NULL, guess)`; `None` when the (ignored) status is a failure.
     fn compute_polar_box(

solidcore/plan/upstream-issues/WR-758--lsi-localsolve-keeps-one-crossing-of-a-leaf-the-line-crosses-twice.md

@@ -0,0 +1,9 @@
+#### WR-758. `SmLSIGlobalSolver::LocalSolve` keeps one crossing of a leaf the line crosses twice
+- **Severity:** Wrong result
+- **Area:** Intersection (line/surface), point classification
+- **Location:** `source/SMLib/src/SmSurface.cpp` (`SmLSIGlobalSolver::LocalSolve`: the guess loop, its test for ray hits behind the start, and the normal-map test after a kept crossing).
+- **What happens:** A leaf is solved from one set of guesses: the middle of the leaf with the middle of the line interval, and further guesses only when that one does not drop onto the leaf. A crossing behind a ray's start is skipped (`continue`) before the polar box test and before the search for a second crossing next to the first. Where a leaf's normal map has a direction perpendicular to the line, the line can cross the leaf twice. When Newton converges to the crossing behind the start, the crossing ahead is never searched for. That search would also look only next to the first crossing, not at a crossing elsewhere in the leaf. A ray that leaves a solid through such a leaf then hits nothing, and `Point3DClassify` puts a point inside the solid in the infinite region.
+- **Trigger:** A ray, or any line, that crosses one leaf of a surface's tree twice. Example: a +x classification ray from just inside a cylinder whose spline spans start away from the top of the circle (SolidCraft's extruded circle, seam at 0.6131 rad, spans 35.1 to 125.1 deg): the line just under the top crosses the leaf from 80.1 to 102.6 deg on either side of 90 deg.
+- **How found:** SolidCraft #12539: an R1 fillet beside a 1 mm chamfer was refused as cutting the part, because the middle of the round's section, 0.08 mm inside an R30.475 wall, classified as outside (solidcore #3120).
+- **Confirmed:** yes, on the port (`surface.surface_p10` matches the C++ bit for bit in difftest). `crates/solidcore-ops/tests/classify_point_grazing_ray.rs`: 57 of 366 points beside the wall classified wrongly before the fix, none after.
+- **Rust port:** Fixed. A leaf whose normal map has a perpendicular to the line is halved until each part's has none, so each part holds at most one crossing (Sederberg and Meyers' normal-cone condition). Each part is solved from its own guesses. Parts the line misses, or reaches only behind a ray's start, are dropped. A leaf without such a direction is solved as upstream solves it.