diff --git a/crates/geom/src/brep/geometry.rs b/crates/geom/src/brep/geometry.rs index 6c0b92fc35..3d6f4e0076 100644 --- a/crates/geom/src/brep/geometry.rs +++ b/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(); diff --git a/crates/kernel/src/step_in/tests.rs b/crates/kernel/src/step_in/tests.rs index 7fdd21d72b..aa623ce210 100644 --- a/crates/kernel/src/step_in/tests.rs +++ b/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); +} diff --git a/crates/step/src/topology.rs b/crates/step/src/topology.rs index 4d1b9f1a36..50a9ad41d7 100644 --- a/crates/step/src/topology.rs +++ b/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>, kind: Kind) -> R { 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")); diff --git a/crates/step/src/topology/seams.rs b/crates/step/src/topology/seams.rs index ef2afae54c..46dd4f9b35 100644 --- a/crates/step/src/topology/seams.rs +++ b/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 { + 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 = ax.poles.iter().flatten().copied().collect(); + if poles.is_empty() { + continue; + } + let mut edges: Vec = + 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> { + 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; 2]) -> R<[i64; 2]> { let mut last: Option<[f64; 2]> = None;