Blends: a curved seam where the edge chain runs on holds the contact on itself (#1291)

A round or chamfer whose edge chain runs on across a curved seam (the end circle of a lofted tube whose axis bends) is now built: the section solve keeps its contact on the seam through a plane taken afresh at the seam's nearest point at each step. Straight seams are unchanged. Robustness seed 41 is re-recorded: its second fillet, refused before, now builds. Fixes #1213.

Co-authored-by: gvocale <[email protected]>

crates/engine/src/tests.rs

@@ -306,6 +306,7 @@ fn join_of_a_lobed_disc_pattern_is_not_refused() {
 mod blend_fit;
 mod boolean_diagnostics;
 mod component_frames;
+mod curved_seam_round;
 mod large_step_import;
 mod loft_cut;
 mod loft_ring;

crates/engine/src/tests/curved_seam_round.rs

@@ -0,0 +1,105 @@
+//! A round or a chamfer on the end circle of a bent lofted tube.
+//!
+//! The wall of such a tube is two faces of the loft that meet in a seam, a curve (the loft bends),
+//! and the seam ends on the circle at each of the circle's two vertices. The blend's section at such
+//! a vertex has its contact on the seam. It was refused for a seam that is not straight ("not
+//! supported yet: a curved seam where the edge chain turns"); now the contact is held on the curved
+//! seam itself.
+
+use std::f64::consts::PI;
+
+use solidcraft_geom::Vec3;
+
+use super::*;
+
+/// Radius of the tube, radius of the bend of its axis, and the turn (radians).
+const TUBE: f64 = 20.0;
+const BEND: f64 = 100.0;
+const TURN: f64 = PI / 3.0;
+
+/// Five circles of the tube's radius, square to an arc of radius `BEND` that turns `TURN` from the
+/// origin (the axis leaves along +x, the arc's centre is (0, BEND)), lofted into one body. The
+/// start circle lies in the plane x = 0.
+fn bent_tube(s: &mut Session) {
+    let mut sections = Vec::new();
+    for i in 0..5 {
+        let t = TURN * i as f64 / 4.0;
+        let name = format!("S{i}");
+        let plane = json!({"origin": [BEND * t.sin(), BEND * (1.0 - t.cos()), 0], "x_dir": [0, 0, 1], "y_dir": [t.sin(), -t.cos(), 0]});
+        run(s, "sketch.create", json!({"plane": plane, "name": name}));
+        run(s, "sketch.circle.center", json!({"center": [0, 0], "radius": TUBE}));
+        run(s, "sketch.finish", json!({}));
+        sections.push(json!({"sketch": name}));
+    }
+    run(s, "solid.loft", json!({"sections": sections, "operation": "new", "name": "Bent", "body_name": "Tube"}));
+}
+
+/// A point of the start circle, 45 degrees round from the arc's plane (nowhere near its vertices,
+/// which are in that plane's perpendicular, the seams' ends).
+fn start_circle() -> Value {
+    let a = PI / 4.0;
+    json!([[0.0, -TUBE * a.cos(), TUBE * a.sin()]])
+}
+
+/// The tube's distance from the axis to a point: how far it is from the core circle of the bend.
+fn core(p: Vec3) -> f64 {
+    ((p.x.hypot(p.y - BEND) - BEND).hypot(p.z)).abs()
+}
+
+/// What a ball of radius `round` rolling along the start circle makes: it touches the end plane
+/// (x = 0) and the wall, so its centre is on x = `round` and `TUBE - round` from the core circle.
+/// The mesh vertices of the round are all `round` from the curve of those centres. Returns how many
+/// vertices were looked at and the worst error of that distance.
+fn off_the_ball(s: &Session, round: f64) -> (usize, f64) {
+    let body = &s.model.state().bodies[0].body;
+    let mesh = body.tessellate(0.01).unwrap();
+    let centres: Vec<Vec3> = (0..720)
+        .map(|i| {
+            let a = 2.0 * PI * i as f64 / 720.0;
+            let rho = BEND + (TUBE - round) * a.cos();
+            Vec3::new(round, BEND - (rho * rho - round * round).sqrt(), (TUBE - round) * a.sin())
+        })
+        .collect();
+    let (mut count, mut worst) = (0, 0.0_f64);
+    for p in &mesh.positions {
+        // On the round: off the end plane, inside the wall, and within reach of the ball.
+        if p.x > 1e-3 && p.x < round + 0.5 && core(*p) < TUBE - 0.02 && core(*p) > TUBE - round - 0.5 {
+            let d = centres.iter().map(|c| c.dist(*p)).fold(f64::INFINITY, f64::min);
+            count += 1;
+            worst = worst.max((d - round).abs());
+        }
+    }
+    (count, worst)
+}
+
+#[test]
+fn a_round_on_the_end_circle_of_a_bent_lofted_tube_is_a_rolling_ball_blend() {
+    let mut s = Session::default();
+    bent_tube(&mut s);
+    let before = volume(&mut s);
+    // The loft is the bent tube itself: pi r^2 times the length of the axis.
+    assert!(rel(before, PI * TUBE * TUBE * BEND * TURN) < 1e-3, "{before}");
+    run(&mut s, "solid.fillet", json!({"edges": start_circle(), "radius": 3.0}));
+    assert_eq!(s.model.state().bodies.len(), 1);
+    // About what the same round takes off a straight tube (a body of revolution, Pappus: the corner
+    // between the ball and the wall, 9 (1 - pi / 4) mm2 at 0.2234 * 3 mm inside the wall).
+    let straight = 2.0 * PI * 9.0 * (1.0 - PI / 4.0) * (TUBE - 0.2234 * 3.0);
+    let removed = before - volume(&mut s);
+    assert!((removed - straight).abs() < 0.08 * straight, "removed {removed}, a straight tube loses {straight}");
+    let (count, worst) = off_the_ball(&s, 3.0);
+    assert!(count > 200, "only {count} vertices on the round");
+    assert!(worst < 0.02, "the round's surface is {worst} mm off the ball's");
+}
+
+#[test]
+fn a_chamfer_on_the_end_circle_of_a_bent_lofted_tube_builds() {
+    let mut s = Session::default();
+    bent_tube(&mut s);
+    let before = volume(&mut s);
+    run(&mut s, "solid.chamfer", json!({"edges": start_circle(), "distance": 3.0}));
+    assert_eq!(s.model.state().bodies.len(), 1);
+    // A straight tube loses a triangle of 3 x 3 / 2 mm2 with its middle a third of 3 mm inside the wall.
+    let straight = 2.0 * PI * 4.5 * (TUBE - 1.0);
+    let removed = before - volume(&mut s);
+    assert!((removed - straight).abs() < 0.08 * straight, "removed {removed}, a straight tube loses {straight}");
+}

crates/kernel/src/curveblend.rs

@@ -421,6 +421,29 @@ struct Section {
     q2: Vec3,
 }
 
+/// The third condition of a section's solve, in place of staying in the plane square to the edge: a contact held at a seam.
+#[derive(Clone, Copy)]
+enum Extra<'a> {
+    /// Side `k`'s contact lies in this plane (a point on it, its normal): beside a straight seam.
+    Plane(usize, Vec3, Vec3),
+    /// Side `k`'s contact lies on this curved seam (a lofted wall's seam between its two halves).
+    Seam(usize, &'a mt::Edge),
+}
+
+/// The seam's point nearest `x` and its unit tangent there.
+fn nearest_on_seam(seam: &mt::Edge, x: Vec3) -> Option<(Vec3, Vec3)> {
+    use mt::{BoundedCurve, ParametricCurve, SearchNearestParameter};
+    let curve = seam.curve();
+    let (lo, hi) = curve.range_tuple();
+    let at = |t: f64| from_p3(curve.subs(t));
+    // The nearest of a few samples starts the search: a long seam holds the point's neighbourhood
+    // far from its ends.
+    let hint = (0..=64).map(|i| lo + (hi - lo) * i as f64 / 64.0).min_by(|a, b| at(*a).dist(x).total_cmp(&at(*b).dist(x)))?;
+    let t = curve.search_nearest_parameter(p3(x), hint, 100)?.clamp(lo, hi);
+    let d = curve.der(t);
+    Some((at(t), Vec3::new(d.x, d.y, d.z).normalized()?))
+}
+
 /// Nearest point of a face's surface and the outward normal there.
 pub(crate) fn project(f: &mt::Face, p: Vec3, hint: Option<(f64, f64)>) -> Option<(Vec3, Vec3, (f64, f64))> {
     use mt::{ParametricSurface, ParametricSurface3D, SearchNearestParameter};
@@ -449,22 +472,12 @@ fn in_range(f: &mt::Face, uv: (f64, f64)) -> bool {
 }
 
 /// Solve for the ball touching both sides near `p`. `extra` adds a third condition (a plane the
-/// side-k contact must lie in); without it the centre stays in the plane through `p` square to
-/// `t`. `sigma` is +1 when the ball sits outside the material (concave edge), −1 inside.
+/// side-k contact must lie in, or a curved seam it must lie on); without it the centre stays in the
+/// plane through `p` square to `t`. `sigma` is +1 when the ball sits outside the material (concave edge), −1 inside.
 /// `cusp` is ±1 at a knife edge (the faces' normals opposite: both leave the edge the same way,
 /// along `cusp · n1 × t`, and the ball lies where the cusp has opened to its width), else 0.
 #[allow(clippy::too_many_arguments)]
-fn solve(
-    f1: &[&mt::Face],
-    f2: &[&mt::Face],
-    p: Vec3,
-    t: Vec3,
-    r: f64,
-    sigma: f64,
-    cusp: f64,
-    extra: Option<(usize, Vec3, Vec3)>,
-    tol: f64,
-) -> Option<Section> {
+fn solve(f1: &[&mt::Face], f2: &[&mt::Face], p: Vec3, t: Vec3, r: f64, sigma: f64, cusp: f64, extra: Option<Extra<'_>>, tol: f64) -> Option<Section> {
     solve_from(f1, f2, p, t, r, sigma, cusp, extra, tol, None)
 }
 
@@ -479,7 +492,7 @@ fn solve_from(
     r: f64,
     sigma: f64,
     cusp: f64,
-    extra: Option<(usize, Vec3, Vec3)>,
+    extra: Option<Extra<'_>>,
     tol: f64,
     start: Option<Vec3>,
 ) -> Option<Section> {
@@ -558,10 +571,19 @@ fn solve_from(
         // Third row: stay in the section plane, or put a contact on the seam plane.
         let (row3, g3) = match extra {
             None => (t, (c - p).dot(t)),
-            Some((k, o, m)) => {
+            Some(Extra::Plane(k, o, m)) => {
                 let q = if k == 1 { q1 } else { q2 };
                 (m, (q - o).dot(m))
             }
+            Some(Extra::Seam(k, seam)) => {
+                // The plane through the seam's nearest point to the contact, square to that side's
+                // surface, taken afresh each step: it holds the contact on the curved seam as the
+                // steps close in.
+                let (q, n) = if k == 1 { (q1, n1) } else { (q2, n2) };
+                let (on_seam, d) = nearest_on_seam(seam, q)?;
+                let m = d.cross(n).normalized()?;
+                (m, (q - on_seam).dot(m))
+            }
         };
         let a = [[n1.x, n1.y, n1.z], [n2.x, n2.y, n2.z], [row3.x, row3.y, row3.z]];
         let b = [-g1, -g2, -g3];
@@ -1126,15 +1148,16 @@ fn curve_blend_with_limit(solid: &Solid, ids: &[mt::EdgeID], size: f64, r: f64,
                 slanted.push((k, fi, p, t, None));
             }
         }
-        let extra = match here.first().cloned().filter(|_| !at_end) {
+        let extra = match here.first().filter(|_| !at_end) {
             None => None,
+            // A curved seam holds the contact on itself (the chain runs on smoothly across it: a
+            // lofted wall's two halves).
+            Some((side, seam)) if !straight(seam, tol * 10.0) => Some(Extra::Seam(*side, seam)),
             Some((side, seam)) => {
+                let side = *side;
                 // The plane through the seam square to that side's surface.
                 let o = from_p3(seam.front().point());
                 let d = (from_p3(seam.back().point()) - o).normalized().ok_or_else(|| KernelError::Failed("seam".into()))?;
-                if !straight(&seam, tol * 10.0) {
-                    return Err(unsupported("a curved seam where the edge chain turns"));
-                }
                 let fs = if side == 1 { &s1 } else { &s2 };
                 // The side's normal at the joint: from its face nearest the joint (any face of a
                 // longer side is no guide here).
@@ -1144,7 +1167,7 @@ fn curve_blend_with_limit(solid: &Solid, ids: &[mt::EdgeID], size: f64, r: f64,
                     .min_by(|a, b| a.0.dist(p).total_cmp(&b.0.dist(p)))
                     .map(|x| x.1)
                     .ok_or_else(|| KernelError::Failed("seam normal".into()))?;
-                Some((side, o, d.cross(nn).normalized().ok_or_else(|| KernelError::Failed("seam plane".into()))?))
+                Some(Extra::Plane(side, o, d.cross(nn).normalized().ok_or_else(|| KernelError::Failed("seam plane".into()))?))
             }
         };
         // On the faces of the edge from this vertex alone (the last edge's at an open chain's