STEP import: an edge through a face's pole is split there

A corner ball's edge may run from its equator over the sphere's pole to a vertex just past
it (Poppy Humanoid's creative_commons_license: 1.8° past, 0.016 mm off the pole). On the
sphere's parameters the pole is a line, and the edge's pcurve would have to jump half a turn
along it: "pcurve starts 2 from the end of the use before it, on another period". Each
edge that passes through a pole of a face it bounds (a sphere's, a revolved face's on its
axis) inside its interval is now split at the nearest point to the pole (golden-section
search from 16 samples), the new vertex placed on the pole; the loop's half turn is then at
a vertex on the pole, where a loop may run along it. Each solid's notes say how many edges
were split. The geometry check takes a vertex at a pole to within the vertex's tolerance, as
the reader already does (SurfaceEval::pole_within).

Without the split the edge's jump read as a 0.024 mm gap and was covered by a tolerance:
the measured gap now only covers what the file's geometry really has.

Test: a spherical triangle (radius 0.5) with an edge over the pole to 1.8° past it reads
with its spherical-excess area and its edge split.

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

crates/geom/src/brep/geometry.rs

@@ -166,11 +166,17 @@ impl Brep {
                             ends.push(if u.forward { [uv[0], uv[2]] } else { [uv[2], uv[0]] });
                         }
                         // Each use continues the one before it (the last closes the loop), except
-                        // along a direction that collapses there (a pole).
+                        // along a direction that collapses there (a pole): at the vertex where
+                        // they meet, to within the vertex's tolerance.
                         for (position, pair) in ends.iter().enumerate() {
                             let before = ends.get(position.checked_sub(1).unwrap_or(ends.len().saturating_sub(1)));
                             let (Some([_, end]), [start, _]) = (before, pair) else { continue };
-                            let pole = surface.pole(*start);
+                            let reach = edges
+                                .get(position)
+                                .and_then(|u| self.edges.get(u.edge as usize).map(|e| e.vertices[if u.forward { 0 } else { 1 }]))
+                                .and_then(|v| self.vertices.get(v as usize))
+                                .map_or(0.0, |v| v.tolerance);
+                            let pole = surface.pole_within(*start, reach);
                             for d in 0..2 {
                                 let (Some(Some(p)), Some(false)) = (period.get(d), pole.get(d)) else { continue };
                                 let gap = (start[d] - end[d]).abs();

crates/kernel/src/step_in/tests.rs

@@ -2039,3 +2039,60 @@ fn spline_edges_off_their_plane_by_less_than_the_file_states_get_their_closed_fo
     let m = measure(&imp.bodies[0].body).unwrap();
     assert!(rel(m.area, 900.0) < 1e-6, "{}", m.area);
 }
+
+/// A spherical triangle (radius 0.5) whose edge runs from the equator over the sphere's pole to
+/// a vertex 1.8° past it, as a corner ball meets its rounds (Poppy Humanoid writes one): on the
+/// sphere's parameters the pole is a line where that edge's pcurve would jump half a turn. The
+/// edge is split at the pole, and the face reads with its area by the spherical excess.
+#[test]
+fn an_edge_over_a_spheres_pole_is_split_there() {
+    let r = 0.5f64;
+    let over = (91.8f64).to_radians();
+    let a = [r, 0.0, 0.0];
+    let b = [0.0, r, 0.0];
+    let p = [r * over.cos(), 0.0, r * over.sin()];
+    let f = |x: f64| format!("{x:.17E}");
+    let pt = |id: u64, q: [f64; 3]| format!("#{id}=CARTESIAN_POINT('',({},{},{}));\n", f(q[0]), f(q[1]), f(q[2]));
+    let dir = |id: u64, q: [f64; 3]| format!("#{id}=DIRECTION('',({},{},{}));\n", f(q[0]), f(q[1]), f(q[2]));
+    let mut t = String::from(HEAD);
+    t += &pt(1, [0.0; 3]);
+    t += &dir(2, [0.0, 0.0, 1.0]);
+    t += &dir(3, [1.0, 0.0, 0.0]);
+    t += "#4=AXIS2_PLACEMENT_3D('',#1,#2,#3);\n#5=SPHERICAL_SURFACE('',#4,0.5);\n";
+    t += &pt(10, a);
+    t += &pt(11, b);
+    t += &pt(12, p);
+    t += "#13=VERTEX_POINT('',#10);\n#14=VERTEX_POINT('',#11);\n#15=VERTEX_POINT('',#12);\n";
+    // A to P in the x–z plane over the pole (normal −y, from x towards z); P to B (normal P × B,
+    // from P); B to A on the equator (normal −z, from y towards x).
+    let n = |u: [f64; 3], v: [f64; 3]| {
+        let c = [u[1] * v[2] - u[2] * v[1], u[2] * v[0] - u[0] * v[2], u[0] * v[1] - u[1] * v[0]];
+        let l = (c[0] * c[0] + c[1] * c[1] + c[2] * c[2]).sqrt();
+        [c[0] / l, c[1] / l, c[2] / l]
+    };
+    let unit = |u: [f64; 3]| [u[0] / r, u[1] / r, u[2] / r];
+    t += &dir(20, [0.0, -1.0, 0.0]);
+    t += "#21=AXIS2_PLACEMENT_3D('',#1,#20,#3);\n#22=CIRCLE('',#21,0.5);\n#23=EDGE_CURVE('',#13,#15,#22,.T.);\n";
+    t += &dir(30, n(p, b));
+    t += &dir(31, unit(p));
+    t += "#32=AXIS2_PLACEMENT_3D('',#1,#30,#31);\n#33=CIRCLE('',#32,0.5);\n#34=EDGE_CURVE('',#15,#14,#33,.T.);\n";
+    t += &dir(40, [0.0, 0.0, -1.0]);
+    t += &dir(41, [0.0, 1.0, 0.0]);
+    t += "#42=AXIS2_PLACEMENT_3D('',#1,#40,#41);\n#43=CIRCLE('',#42,0.5);\n#44=EDGE_CURVE('',#14,#13,#43,.T.);\n";
+    // Counter-clockwise seen from outside: A, B, P, so every edge the other way round.
+    t += "#50=ORIENTED_EDGE('',*,*,#44,.F.);\n#51=ORIENTED_EDGE('',*,*,#34,.F.);\n#52=ORIENTED_EDGE('',*,*,#23,.F.);\n";
+    t += "#53=EDGE_LOOP('',(#50,#51,#52));\n#54=FACE_OUTER_BOUND('',#53,.T.);\n#55=ADVANCED_FACE('',(#54),#5,.T.);\n";
+    t += "#60=OPEN_SHELL('',(#55));\n#61=SHELL_BASED_SURFACE_MODEL('Ball corner',(#60));\n#62=MANIFOLD_SURFACE_SHAPE_REPRESENTATION('',(#61),#70);\n";
+    t += &context(70, "(LENGTH_UNIT()NAMED_UNIT(*)SI_UNIT(.MILLI.,.METRE.))");
+    t += TAIL;
+    let imp = step_import(&t).unwrap();
+    assert_eq!(imp.bodies.len(), 1, "{:?}", imp.warnings);
+    assert!(imp.warnings.iter().any(|w| w.contains("1 edges split where they pass through a face's pole")), "{:?}", imp.warnings);
+    // The spherical excess of the unit vectors' triangle, tan(E/2) = |a·(b×c)| / (1 + a·b + b·c + c·a).
+    let (ua, ub, up) = (unit(a), unit(b), unit(p));
+    let dot = |u: [f64; 3], v: [f64; 3]| u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
+    let triple = dot(ua, n(ub, up)) * (1.0 - dot(ub, up).powi(2)).sqrt();
+    let excess = 2.0 * triple.abs().atan2(1.0 + dot(ua, ub) + dot(ub, up) + dot(up, ua));
+    let m = measure(&imp.bodies[0].body).unwrap();
+    assert!(rel(m.area, excess * r * r) < 1e-3, "{} vs {}", m.area, excess * r * r);
+}

crates/step/src/topology.rs

@@ -794,6 +794,10 @@ impl<'b, 'a> Builder<'b, 'a> {
     fn finish(mut self, item: u64, name: String, shells: Vec<Vec<(u32, bool)>>, kind: Kind) -> R<Solid> {
         let _span = solidcraft_profile::span("import", "step.solid_record");
         let mut notes = Vec::new();
+        let at_poles = self.split_at_poles()?;
+        if at_poles > 0 {
+            notes.push(format!("{at_poles} edges split where they pass through a face's pole"));
+        }
         let seams = self.add_seams()?;
         if seams > 0 {
             notes.push(format!("{seams} seam edges added to faces that wind around their surface"));

crates/step/src/topology/seams.rs

@@ -207,6 +207,69 @@ impl Builder<'_, '_> {
         Ok(added)
     }
 
+    /// Split every edge that passes through a pole of a face it bounds (a sphere's, a revolved
+    /// face's on its axis) inside its interval: a great circle from the equator over the pole to a
+    /// vertex just past it. On the face's parameters the pole is a line, along which the edge's
+    /// pcurve would have to jump half a turn; split there, the jump is at a vertex on the pole,
+    /// where a loop may run along the pole's line. Returns how many edges were split.
+    pub(super) fn split_at_poles(&mut self) -> R<usize> {
+        let mut splits = 0;
+        for fi in 0..self.faces.len() {
+            let Some(face) = self.faces.get(fi) else { continue };
+            let Some(ax) = axis(&face.surface) else { continue };
+            let poles: Vec<Vec3> = ax.poles.iter().flatten().copied().collect();
+            if poles.is_empty() {
+                continue;
+            }
+            let mut edges: Vec<u32> =
+                face.loops.iter().filter_map(|l| if let LoopB::Edges(es) = l { Some(es) } else { None }).flatten().map(|u| u.0).collect();
+            edges.sort_unstable();
+            edges.dedup();
+            for e in edges {
+                for pole in &poles {
+                    let Some(t) = self.through_at(e, *pole)? else { continue };
+                    splits += 1;
+                    if splits > MAX_SPLITS {
+                        return Err(format!("more than {MAX_SPLITS} edges split at poles"));
+                    }
+                    let v = self.split_edge(e, t)?;
+                    // The new vertex is the pole, as the face's parameters have it.
+                    if let Some(x) = self.vertices.get_mut(v as usize) {
+                        let gap = Vec3::from(x.point).dist(*pole);
+                        x.point = arr(*pole);
+                        x.tolerance = x.tolerance.max(1.5 * gap);
+                    }
+                }
+            }
+        }
+        Ok(splits)
+    }
+
+    /// Where edge `e` passes through `pole` inside its interval (within the tolerance, away from
+    /// its ends), if it does: its nearest point to the pole, refined by golden-section search.
+    fn through_at(&self, e: u32, pole: Vec3) -> R<Option<f64>> {
+        let edge = self.edges.get(e as usize).ok_or("edge reference")?;
+        let c = CurveEval::new(&edge.curve).map_err(|m| m.to_string())?;
+        let [lo, hi] = edge.interval;
+        let reach = self.tol.max(edge.tolerance);
+        let d = |t: f64| c.point(t).map_or(f64::INFINITY, |p| p.dist(pole));
+        let step = (hi - lo) / SAMPLES as f64;
+        let k = (0..=SAMPLES).min_by(|&i, &j| d(lo + step * i as f64).total_cmp(&d(lo + step * j as f64))).unwrap_or(0);
+        let (mut a, mut b) = (lo + step * k.saturating_sub(1) as f64, (lo + step * (k + 1) as f64).min(hi));
+        let g = 0.5 * (5f64.sqrt() - 1.0);
+        for _ in 0..100 {
+            let (x, y) = (b - g * (b - a), a + g * (b - a));
+            if d(x) < d(y) {
+                b = y;
+            } else {
+                a = x;
+            }
+        }
+        let t = 0.5 * (a + b);
+        let ends = [lo, hi].map(|s| c.point(s).map_or(0.0, |p| p.dist(pole)));
+        Ok((d(t) <= reach && ends.iter().all(|x| *x > reach) && t > lo && t < hi).then_some(t))
+    }
+
     /// Turns (in u, in v) a loop makes around a periodic surface.
     fn winding(&self, eval: &SurfaceEval, es: &[(u32, bool)], periods: [Option<f64>; 2]) -> R<[i64; 2]> {
         let mut last: Option<[f64; 2]> = None;