M2: loft (arc-length ruled sections) and sweep along sketch paths; three-level measurement voting

Oracle: 23/29 parts match Fusion.

Co-Authored-By: Claude Opus 5.5 <[email protected]>

crates/kernel/src/build.rs

@@ -395,3 +395,246 @@ fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2)
         Body::new(solid).map(Some)
     })
 }
+
+/// Close a ruled shell (side faces between matching wires) with planar caps into a body.
+fn cap_and_close(mut faces: Vec<mt::Face>, w0: &[mt::Wire], w1: &[mt::Wire]) -> Result<Body> {
+    let bottom = builder::try_attach_plane(&w0.iter().map(|w| w.inverse()).collect::<Vec<_>>()).map_err(|e| KernelError::Failed(format!("cap: {e}")))?;
+    let top = builder::try_attach_plane(w1).map_err(|e| KernelError::Failed(format!("cap: {e}")))?;
+    let mut last = String::new();
+    for (b, t) in [(bottom.clone(), top.clone()), (bottom.inverse(), top.inverse()), (bottom.clone(), top.inverse()), (bottom.inverse(), top)] {
+        let mut all = std::mem::take(&mut faces);
+        let keep = all.clone();
+        all.push(b);
+        all.push(t);
+        let shell: mt::Shell = all.into();
+        match Solid::try_new(vec![shell]) {
+            Ok(s) => return Body::new(s),
+            Err(e) => last = e.to_string(),
+        }
+        faces = keep;
+    }
+    Err(KernelError::Failed(format!("closing the shell: {last}")))
+}
+
+/// Re-split a loop into exactly `n` segments: a full circle gets breaks at `angles` (or evenly),
+/// other loops split their longest segments in half.
+fn resample(lp: &Loop2, n: usize, angles: &[f64]) -> Loop2 {
+    let circle = lp.segs.iter().all(|s| matches!(s, Seg2::Arc { .. }))
+        && lp.segs.windows(2).all(|w| match (w[0], w[1]) {
+            (Seg2::Arc { center: c1, radius: r1, .. }, Seg2::Arc { center: c2, radius: r2, .. }) => c1.dist(c2) < 1e-9 && (r1 - r2).abs() < 1e-9,
+            _ => false,
+        });
+    if circle && let Some(Seg2::Arc { center, radius, .. }) = lp.segs.first().copied() {
+        let mut a: Vec<f64> = if angles.len() == n { angles.to_vec() } else { (0..n).map(|i| std::f64::consts::TAU * i as f64 / n as f64).collect() };
+        // Keep the given order but make them increase.
+        for i in 1..a.len() {
+            while a[i] <= a[i - 1] {
+                a[i] += std::f64::consts::TAU;
+            }
+        }
+        let segs = (0..n)
+            .map(|i| {
+                let s = a[i];
+                let e = if i + 1 < n { a[i + 1] } else { a[0] + std::f64::consts::TAU };
+                Seg2::Arc { center, radius, start: s, sweep: e - s }
+            })
+            .collect();
+        return Loop2 { segs };
+    }
+    let mut segs = lp.segs.clone();
+    while segs.len() < n {
+        let Some((k, _)) = segs.iter().enumerate().max_by(|a, b| a.1.length().total_cmp(&b.1.length())) else { break };
+        let s = segs.remove(k);
+        let (p, q) = match s {
+            Seg2::Line { a, b } => {
+                let m = a.lerp(b, 0.5);
+                (Seg2::Line { a, b: m }, Seg2::Line { a: m, b })
+            }
+            Seg2::Arc { center, radius, start, sweep } => {
+                (Seg2::Arc { center, radius, start, sweep: sweep / 2.0 }, Seg2::Arc { center, radius, start: start + sweep / 2.0, sweep: sweep / 2.0 })
+            }
+        };
+        segs.insert(k, q);
+        segs.insert(k, p);
+    }
+    Loop2 { segs }
+}
+
+/// Ruled loft through two or more planar sections (outer loops only).
+pub fn loft(sections: &[(Plane, Loop2)]) -> Result<Body> {
+    if sections.len() < 2 || sections.len() > 50 {
+        return Err(KernelError::Invalid("a loft needs 2…50 sections".into()));
+    }
+    let n = sections.iter().map(|(_, l)| l.segs.len()).max().unwrap_or(0);
+    if n == 0 || n > MAX_SEGS {
+        return Err(KernelError::Invalid("empty section".into()));
+    }
+    // Orient every section counter-clockwise about the first section's normal, then match
+    // segment counts and starting points to the section with the most segments.
+    let n0 = sections[0].0.normal();
+    let oriented: Vec<(Plane, Loop2)> = sections
+        .iter()
+        .map(|(p, l)| {
+            let l = if p.normal().dot(n0) >= 0.0 { l.ccw() } else { l.ccw().reversed() };
+            (*p, l)
+        })
+        .collect();
+    let lead = oriented.iter().max_by_key(|(_, l)| l.segs.len()).cloned().ok_or_else(|| KernelError::Invalid("no sections".into()))?;
+    let lead_c = lead.0.to_world(lead.1.centroid());
+    let lead_dirs: Vec<Vec3> = lead.1.segs.iter().map(|s| lead.0.to_world(s.start()) - lead_c).collect();
+    let mut wires: Vec<mt::Wire> = Vec::new();
+    for (p, l) in &oriented {
+        let c = l.centroid();
+        // Break angles: the lead section's vertex directions seen in this plane.
+        let angles: Vec<f64> = lead_dirs.iter().map(|d| Vec2::new(d.dot(p.x), d.dot(p.y)).angle()).collect();
+        let mut r = resample(l, n, &angles);
+        // Start at the vertex closest in direction to the lead's first vertex.
+        if let Some(d0) = lead_dirs.first() {
+            let target = Vec2::new(d0.dot(p.x), d0.dot(p.y)).angle();
+            let best = (0..r.segs.len()).min_by(|a, b| {
+                let ang = |k: usize| {
+                    let v = r.segs[k].start() - c;
+                    let d = (v.angle() - target).rem_euclid(std::f64::consts::TAU);
+                    d.min(std::f64::consts::TAU - d)
+                };
+                ang(*a).total_cmp(&ang(*b))
+            });
+            if let Some(k) = best {
+                r.segs.rotate_left(k);
+            }
+        }
+        wires.push(uniform_wire(p, &r)?);
+    }
+    guard("loft", || {
+        let mut faces: Vec<mt::Face> = Vec::new();
+        for w in wires.windows(2) {
+            let sides = builder::try_wire_homotopy(&w[0], &w[1]).map_err(|e| KernelError::Failed(format!("loft: {e}")))?;
+            faces.extend(sides.face_iter().cloned());
+        }
+        let (Some(first), Some(last)) = (wires.first(), wires.last()) else { return Err(KernelError::Failed("loft".into())) };
+        cap_and_close(faces, std::slice::from_ref(first), std::slice::from_ref(last))
+    })
+}
+
+/// A 3D path segment for sweeps.
+#[derive(Clone, Copy, Debug)]
+pub enum PathSeg {
+    Line { a: Vec3, b: Vec3 },
+    /// Arc around `center`, rotating by `angle` about `axis` (right hand) from `a`.
+    Arc { a: Vec3, center: Vec3, axis: Vec3, angle: f64 },
+}
+
+/// Sweep a planar region (outer loop) along a chain of path segments, keeping the profile
+/// perpendicular to the path (lines translate it, arcs rotate it).
+pub fn sweep(plane: &Plane, region: &Region2, path: &[PathSeg]) -> Result<Body> {
+    region_ok(region)?;
+    if path.is_empty() || path.len() > 1000 {
+        return Err(KernelError::Invalid("the path needs 1…1000 segments".into()));
+    }
+    let w0 = wire(plane, &region.outer.ccw())?;
+    guard("sweep", || {
+        let mut shell = mt::Shell::new();
+        let mut cur = w0.clone();
+        for seg in path {
+            let mut part: mt::Shell = match *seg {
+                PathSeg::Line { a, b } => builder::tsweep(&cur, v3(b - a)),
+                PathSeg::Arc { center, axis, angle, .. } => {
+                    let ax = axis.normalized().ok_or_else(|| KernelError::Invalid("arc axis".into()))?;
+                    builder::rsweep(&cur, p3(center), v3(ax), mt::Rad(angle))
+                }
+            };
+            let next = part.extract_boundaries().into_iter().find(|w| w.edge_iter().all(|e| !cur.edge_iter().any(|c| c.id() == e.id()))).ok_or_else(|| KernelError::Failed("sweep: lost the profile".into()))?;
+            shell.append(&mut part);
+            cur = next.inverse();
+        }
+        let faces: Vec<mt::Face> = shell.face_iter().cloned().collect();
+        cap_and_close(faces, std::slice::from_ref(&w0), std::slice::from_ref(&cur))
+    })
+}
+
+/// Cox–de Boor basis values N_{j,3}(u) for all j.
+fn cubic_basis(knots: &[f64], n_ctrl: usize, u: f64) -> Vec<f64> {
+    let p = 3;
+    let m = knots.len();
+    let mut n = vec![0.0; m.saturating_sub(1)];
+    let last = n_ctrl.saturating_sub(1);
+    for (j, v) in n.iter_mut().enumerate() {
+        let (Some(&a), Some(&b)) = (knots.get(j), knots.get(j + 1)) else { continue };
+        // Half-open spans, with the last span closed at the end.
+        if (a <= u && u < b) || (u >= b && (b - a) > 0.0 && j == last && (u - b).abs() < 1e-12) {
+            *v = 1.0;
+        }
+    }
+    for k in 1..=p {
+        for j in 0..m.saturating_sub(1 + k) {
+            let (Some(&tj), Some(&tjk), Some(&tj1), Some(&tjk1)) = (knots.get(j), knots.get(j + k), knots.get(j + 1), knots.get(j + k + 1)) else { continue };
+            let left = if tjk - tj > 0.0 { (u - tj) / (tjk - tj) * n.get(j).copied().unwrap_or(0.0) } else { 0.0 };
+            let right = if tjk1 - tj1 > 0.0 { (tjk1 - u) / (tjk1 - tj1) * n.get(j + 1).copied().unwrap_or(0.0) } else { 0.0 };
+            if let Some(v) = n.get_mut(j) {
+                *v = left + right;
+            }
+        }
+    }
+    n.truncate(n_ctrl);
+    n
+}
+
+/// Cubic B-spline through points at uniform parameters (clamped, knots by averaging).
+fn interpolate_cubic(pts: &[Vec3]) -> Option<mt::BSplineCurve<mt::Point3>> {
+    let m = pts.len();
+    if m < 4 {
+        return None;
+    }
+    let u: Vec<f64> = (0..m).map(|i| i as f64 / (m - 1) as f64).collect();
+    let mut knots = vec![0.0; 4];
+    for j in 1..(m - 3) {
+        knots.push((j..j + 3).filter_map(|i| u.get(i)).sum::<f64>() / 3.0);
+    }
+    knots.extend([1.0; 4]);
+    // Dense solve of N · P = Q for each coordinate.
+    let mut a: Vec<Vec<f64>> = u.iter().map(|ui| cubic_basis(&knots, m, *ui)).collect();
+    let mut q: Vec<[f64; 3]> = pts.iter().map(|p| [p.x, p.y, p.z]).collect();
+    for col in 0..m {
+        let piv = (col..m).max_by(|x, y| a[*x][col].abs().total_cmp(&a[*y][col].abs()))?;
+        if a[piv][col].abs() < 1e-14 {
+            return None;
+        }
+        a.swap(col, piv);
+        q.swap(col, piv);
+        for r in 0..m {
+            if r == col {
+                continue;
+            }
+            let f = a[r][col] / a[col][col];
+            if f == 0.0 {
+                continue;
+            }
+            for c in col..m {
+                let v = a[col][c];
+                a[r][c] -= f * v;
+            }
+            let qc = q[col];
+            for k in 0..3 {
+                q[r][k] -= f * qc[k];
+            }
+        }
+    }
+    let ctrl: Vec<mt::Point3> = (0..m).map(|i| mt::Point3::new(q[i][0] / a[i][i], q[i][1] / a[i][i], q[i][2] / a[i][i])).collect();
+    mt::BSplineCurve::try_new(mt::KnotVec::from(knots), ctrl).ok()
+}
+
+/// A closed wire whose edges are cubic B-splines parameterised (nearly) by arc length, so that
+/// ruled surfaces between sections pair points proportionally along each edge.
+fn uniform_wire(plane: &Plane, lp: &Loop2) -> Result<mt::Wire> {
+    const SAMPLES: usize = 24;
+    let n = lp.segs.len();
+    let verts: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(plane.to_world(s.start())))).collect();
+    let mut edges = Vec::with_capacity(n);
+    for (i, s) in lp.segs.iter().enumerate() {
+        let (Some(a), Some(b)) = (verts.get(i), verts.get((i + 1) % n)) else { continue };
+        let pts: Vec<Vec3> = (0..SAMPLES).map(|k| plane.to_world(s.point_at(k as f64 / (SAMPLES - 1) as f64))).collect();
+        let curve = interpolate_cubic(&pts).ok_or_else(|| KernelError::Failed("loft curve".into()))?;
+        edges.push(mt::Edge::new(a, b, mt::Curve::BSplineCurve(curve)));
+    }
+    Ok(edges.into())
+}

crates/kernel/src/lib.rs

@@ -23,7 +23,7 @@ mod topo;
 
 pub use blend::{chamfer, fillet};
 pub use body::{Body, EdgeInfo, FaceInfo};
-pub use build::{box_solid, cylinder, extrude, extrude_tapered, revolve, sphere, torus};
+pub use build::{PathSeg, box_solid, cylinder, extrude, extrude_tapered, loft, revolve, sphere, sweep, torus};
 pub use measure::{BodyMeasure, measure};
 pub use ops::{BoolOp, boolean, split_by_plane, transform, transform_matrix};
 pub use polyhedron::{HalfSpace, convex_polyhedron, draft, shell};

crates/kernel/src/measure.rs

@@ -40,13 +40,21 @@ fn face_sums(m: &Mesh, nf: usize) -> Vec<(f64, f64, Vec3)> {
 pub fn measure(b: &Body) -> Result<BodyMeasure> {
     let size = b.size();
     let fine = b.tessellate((size * 5e-5).max(1e-4))?;
+    let medium = b.tessellate((size * 2.5e-4).max(2e-4))?;
     let coarse = b.tessellate((size * 1e-3).max(1e-3))?;
     let nf = b.face_count();
-    let (sf, sc) = (face_sums(&fine, nf), face_sums(&coarse, nf));
+    let (sf, sm, sc) = (face_sums(&fine, nf), face_sums(&medium, nf), face_sums(&coarse, nf));
+    let close = |x: f64, y: f64| (x - y).abs() <= 0.01 * x.abs().max(y.abs()) + 1e-9;
     let (mut area, mut v6, mut mom) = (0.0, 0.0, Vec3::ZERO);
-    for (f, c) in sf.iter().zip(&sc) {
-        let bad = (f.0 - c.0).abs() > 0.01 * c.0.max(f.0) + 1e-9;
-        let pick = if bad { c } else { f };
+    for ((f, m), c) in sf.iter().zip(&sm).zip(&sc) {
+        // Trust the finest tessellation that agrees with the next coarser one.
+        let pick = if close(f.0, m.0) {
+            f
+        } else if close(m.0, c.0) {
+            m
+        } else {
+            f
+        };
         area += pick.0;
         v6 += pick.1;
         mom += pick.2;

crates/kernel/src/tests.rs

@@ -192,6 +192,30 @@ fn polyhedra_shell_and_draft() {
     assert!(shell(&cyl, &[Vec3::new(0.0, 0.0, 5.0)], 1.0).is_err(), "curved faces");
 }
 
+#[test]
+fn loft_and_sweep() {
+    let sq = Loop2::polygon(&[Vec2::new(-20.0, -20.0), Vec2::new(20.0, -20.0), Vec2::new(20.0, 20.0), Vec2::new(-20.0, 20.0)]);
+    let top = Plane::XY.offset(40.0);
+    let l = loft(&[(Plane::XY, sq), (top, Loop2::circle(Vec2::ZERO, 10.0))]).unwrap();
+    let m = measure(&l).unwrap();
+    println!("loft volume {} faces {:?}", m.volume, m.merged);
+    assert!(m.volume > 314.0 * 40.0 && m.volume < 1600.0 * 40.0);
+    // A tube: up 30, quarter bend of radius 20, across 30.
+    let disc = Region2 { outer: Loop2::circle(Vec2::ZERO, 3.0), holes: vec![] };
+    let path = [
+        PathSeg::Line { a: Vec3::ZERO, b: Vec3::new(0.0, 0.0, 30.0) },
+        PathSeg::Arc { a: Vec3::new(0.0, 0.0, 30.0), center: Vec3::new(20.0, 0.0, 30.0), axis: Vec3::new(0.0, 1.0, 0.0), angle: std::f64::consts::FRAC_PI_2 },
+        PathSeg::Line { a: Vec3::new(20.0, 0.0, 50.0), b: Vec3::new(50.0, 0.0, 50.0) },
+    ];
+    let t = sweep(&Plane::XY, &disc, &path).unwrap();
+    for tol in [0.1, 0.01, 0.003, 0.001] {
+        println!("tube tess {tol}: {}", t.tessellate(tol).unwrap().measure().volume);
+    }
+    let area = PI * 9.0;
+    let expect = area * (30.0 + 30.0 + 20.0 * std::f64::consts::FRAC_PI_2);
+    assert!(rel(measure(&t).unwrap().volume, expect) < 1e-3, "{} vs {expect}", measure(&t).unwrap().volume);
+}
+
 #[test]
 #[ignore]
 fn debug_convergence() {