diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index 345faf6d66..001e8d0a40 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -1876,14 +1876,14 @@ pub fn loft(sections: &[(Plane, Loop2)]) -> Result { /// sections between. Without rails, the plain loft. pub fn loft_with_rails(sections: &[(Plane, Loop2)], rails: &[Vec], fit: RailFit) -> Result { if crate::backend()? == crate::Backend::Smlib { - if !rails.is_empty() { - return Err(crate::smlib_sweeps::refuse("loft along guide rails")); - } let reversed = loft_reversed(sections)?; let numbered = sections.iter().map(|(p, _)| *p).zip(numbered::number_loops(sections.iter().map(|(_, l)| l))).collect::>(); let loops = loft_numbered_loops(&numbered, reversed)?; let first = if reversed { loops.len().saturating_sub(1) } else { 0 }; let pieces: Vec<(Plane, Vec)> = loops.iter().map(|(p, l)| (*p, crate::smlib_sweeps::pieces([l]))).collect(); + if !rails.is_empty() { + return crate::smlib_sweeps::loft_rails(&pieces, first, rails, fit); + } return crate::smlib_sweeps::loft(&pieces, first); } let (wires, edge_sources, reversed) = loft_wires(sections)?; @@ -2349,31 +2349,27 @@ pub(crate) fn rotate_about(p: Vec3, center: Vec3, k: Vec3, angle: f64) -> Vec3 { center + v * c + k.cross(v) * s + k * (k.dot(v) * (1.0 - c)) } -/// [`sweep_with`] through SMLib's curve sweep: the profile turning with the path or only moving -/// along it (orientation parallel), along a path without corners (SMLib's sweep follows smooth -/// paths). Twists, end scales, guide rails and orientation directions have no SMLib sweep. +/// [`sweep_with`] through SMLib, along a path without corners (SMLib's sweeps follow smooth +/// paths): its curve sweep for the profile turning with the path or only moving along it +/// (orientation parallel); a twist, an end scale, guide rails or an orientation direction place +/// the profile at stations along the path ([`crate::smlib_sweeps::sweep_controlled`]). fn smlib_sweep(plane: &Plane, region: &Region2, path: &[PathSeg], ctl: &crate::SweepControl) -> Result { ctl.check()?; region_ok(region)?; if path.is_empty() || path.len() > 1000 { return Err(KernelError::Invalid("the path needs 1…1000 segments".into())); } - let translate = match ctl.orientation { - crate::SweepOrientation::Perpendicular => false, - crate::SweepOrientation::Parallel => true, - crate::SweepOrientation::Direction(_) => return Err(crate::smlib_sweeps::refuse("sweep with an orientation direction")), - }; - if ctl.twist != 0.0 || ctl.scale != 1.0 { - return Err(crate::smlib_sweeps::refuse("sweep with a twist or an end scale")); - } - if !ctl.rails.is_empty() { - return Err(crate::smlib_sweeps::refuse("sweep along guide rails")); - } let ends: Vec<(Vec3, Vec3)> = path.iter().map(PathSeg::tangents).collect::>().ok_or_else(|| KernelError::Invalid("a path segment has no direction".into()))?; if ends.windows(2).any(|w| w[0].1.dot(w[1].0) < (1e-3f64).cos()) { return Err(crate::smlib_sweeps::refuse("sweep along a path with corners (SMLib's sweep follows smooth paths)")); } + let translate = match ctl.orientation { + _ if ctl.varies() => return crate::smlib_sweeps::sweep_controlled(plane, &oriented(region)?, path, ctl), + crate::SweepOrientation::Perpendicular => false, + crate::SweepOrientation::Parallel => true, + crate::SweepOrientation::Direction(_) => return crate::smlib_sweeps::sweep_controlled(plane, &oriented(region)?, path, ctl), + }; crate::smlib_sweeps::sweep(plane, &crate::smlib_sweeps::pieces(oriented(region)?.loops()), path, translate) } diff --git a/crates/kernel/src/loft_rails.rs b/crates/kernel/src/loft_rails.rs index 72acd5e6e6..84200ce63d 100644 --- a/crates/kernel/src/loft_rails.rs +++ b/crates/kernel/src/loft_rails.rs @@ -18,9 +18,9 @@ use crate::sweep_path::RailFit; use crate::{KernelError, Result}; /// Most wires a railed loft is skinned through, sections and stations together. -const MAX_WIRES: usize = 33; +pub(crate) const MAX_WIRES: usize = 33; /// Most stations between two sections. -const MAX_PER_SPAN: usize = 15; +pub(crate) const MAX_PER_SPAN: usize = 15; type Curve = mt::BSplineCurve; @@ -42,14 +42,14 @@ struct Rail { } /// An affine move: `x` goes to `to + m·(x − from)` (`m` by columns). -struct Affine { +pub(crate) struct Affine { from: Vec3, to: Vec3, m: [Vec3; 3], } impl Affine { - fn apply(&self, x: Vec3) -> Vec3 { + pub(crate) fn apply(&self, x: Vec3) -> Vec3 { let d = x - self.from; self.to + self.m[0] * d.x + self.m[1] * d.y + self.m[2] * d.z } @@ -259,7 +259,7 @@ fn nearest_on(section: &[Curve], p: Vec3) -> Option<(f64, Vec3, f64)> { /// A wire parameter `f` of the way from `g0` to `g1`; round a closed wire of `n` edges the /// shorter way. -fn interpolate_g(g0: f64, g1: f64, f: f64, n: usize, closed: bool) -> f64 { +pub(crate) fn interpolate_g(g0: f64, g1: f64, f: f64, n: usize, closed: bool) -> f64 { let n = n as f64; if !closed { return (g0 + f * (g1 - g0)).clamp(0.0, n); @@ -294,7 +294,7 @@ fn centroid(ctrl: &[Vec]) -> Vec3 { /// rail's point), the station's middle `middle`, the loft's direction there `dir` (unit); `size` /// the loft's size. One rail: turned about `dir` and sized about the middle, shifted along `dir`; /// two: carried with the first, turned and sized with the line between them, leaning along `dir`. -fn station_map(pairs: &[(Vec3, Vec3)], middle: Vec3, dir: Vec3, fit: RailFit, size: f64) -> Result { +pub(crate) fn station_map(pairs: &[(Vec3, Vec3)], middle: Vec3, dir: Vec3, fit: RailFit, size: f64) -> Result { let (from, to, a, r, two) = match pairs { [(p, q)] => (middle, middle, *p - middle, *q - middle, false), [(p1, q1), (p2, q2)] => (*p1, *q1, *p2 - *p1, *q2 - *q1, true), diff --git a/crates/kernel/src/smlib_sweeps.rs b/crates/kernel/src/smlib_sweeps.rs index 83782aa789..9e19b3202c 100644 --- a/crates/kernel/src/smlib_sweeps.rs +++ b/crates/kernel/src/smlib_sweeps.rs @@ -4,6 +4,11 @@ //! record as a kernel body. A result SMLib's own check refuses is an error, never a body. An //! option SMLib has no sweep for is "not supported yet", never truck's sweep in its place. //! +//! A sweep that twists, scales, keeps an orientation direction or follows guide rails, and a +//! loft that follows guide rails, place the profile at stations by the laws the kernel's own +//! sweeps and lofts define ([`crate::sweep_path`], [`crate::loft_rails`]); each section is the +//! profile's exact image, and SMLib skins the walls through them. +//! //! History is in the builders' vocabulary ([`FaceOrigin::Profile`], [`FaceOrigin::Cap`]), found //! from the result's geometry ([`crate::provenance`] does the same for operations without exact //! history): a wall comes from the profile segments whose curves its edge at the start runs @@ -63,7 +68,7 @@ pub(crate) struct Draft { } #[cfg(feature = "smlib")] -pub(crate) use imp::{extrude, loft, revolve, screw, sweep}; +pub(crate) use imp::{extrude, loft, loft_rails, revolve, screw, sweep, sweep_controlled}; #[cfg(not(feature = "smlib"))] mod stubs { @@ -87,9 +92,15 @@ mod stubs { pub(crate) fn screw(_: &Plane, _: &[Piece], _: Vec3, _: Vec3, _: f64, _: f64) -> Result { Err(not_built()) } + pub(crate) fn sweep_controlled(_: &Plane, _: &NumberedRegion, _: &[crate::PathSeg], _: &crate::SweepControl) -> Result { + Err(not_built()) + } + pub(crate) fn loft_rails(_: &[(Plane, Vec)], _: usize, _: &[Vec], _: crate::RailFit) -> Result { + Err(not_built()) + } } #[cfg(not(feature = "smlib"))] -pub(crate) use stubs::{extrude, loft, revolve, screw, sweep}; +pub(crate) use stubs::{extrude, loft, loft_rails, revolve, screw, sweep, sweep_controlled}; /// Refuses what no SMLib sweep does, before any work. pub(crate) fn refuse(what: &str) -> KernelError { @@ -100,7 +111,7 @@ pub(crate) fn refuse(what: &str) -> KernelError { mod imp { use std::collections::BTreeSet; - use solidcraft_smlib::sweeps::{Curve, SweepFrame}; + use solidcraft_smlib::sweeps::{Affine, Curve, Ends, SweepFrame}; use solidcraft_smlib::{Brep, CheckLevel}; use solidcraft_geom::Vec2; @@ -423,6 +434,110 @@ mod imp { with_history(body, &roles, pieces) } + /// The length of a path segment. + fn seg_len(s: &crate::PathSeg) -> f64 { + match *s { + crate::PathSeg::Line { a, b } => a.dist(b), + crate::PathSeg::Arc { a, center, angle, .. } => angle.abs() * a.dist(center), + } + } + + /// The point `f` (0…1) of the way along a path segment, and the unit direction of travel there. + fn seg_at(s: &crate::PathSeg, f: f64) -> Option<(Vec3, Vec3)> { + match *s { + crate::PathSeg::Line { a, b } => Some((a.lerp(b, f), (b - a).normalized()?)), + crate::PathSeg::Arc { a, center, axis, angle } => { + let k = axis.normalized()?; + let q = crate::build::rotate_about(a, center, k, angle * f); + Some((q, (k.cross(q - center) * angle.signum()).normalized()?)) + } + } + } + + /// `n + 1` points evenly spaced by length along `path` (exact on its lines and arcs), with the + /// unit directions of travel there. + fn path_samples(path: &[crate::PathSeg], n: usize) -> Result<(Vec, Vec)> { + let lens: Vec = path.iter().map(seg_len).collect(); + let total: f64 = lens.iter().sum(); + if !(total > 1e-9 && total.is_finite()) || n == 0 { + return Err(KernelError::Invalid("the sweep path has no length".into())); + } + let bad = || KernelError::Invalid("a path segment has no direction".into()); + let (mut pts, mut tangents) = (Vec::with_capacity(n + 1), Vec::with_capacity(n + 1)); + let (mut k, mut before) = (0usize, 0.0f64); + for i in 0..=n { + let at = total * i as f64 / n as f64; + while k + 1 < path.len() && lens.get(k).is_some_and(|l| before + l < at) { + before += lens.get(k).copied().unwrap_or(0.0); + k += 1; + } + let (seg, len) = path.get(k).zip(lens.get(k)).ok_or_else(bad)?; + let f = if *len > 0.0 { ((at - before) / len).clamp(0.0, 1.0) } else { 0.0 }; + let (q, t) = seg_at(seg, f).ok_or_else(bad)?; + pts.push(q); + tangents.push(t); + } + Ok((pts, tangents)) + } + + /// Fewest and most stations of a sweep that twists, scales or follows rails. + const MIN_STATIONS: usize = 16; + const MAX_STATIONS: usize = 2048; + + /// The profile `region` on `plane` swept along the smooth `path` the way `ctl` asks: twisted, + /// scaled, kept at an orientation or guided by rails. The profile is placed at stations along + /// the path by [`crate::sweep_path`]'s laws (the same the polyline sweep uses), each section + /// exactly the profile's image, and SMLib skins the walls through them. Round a closed path the + /// last section has to land on the first, and the ends are joined. + pub(crate) fn sweep_controlled(plane: &Plane, region: &NumberedRegion, path: &[crate::PathSeg], ctl: &crate::SweepControl) -> Result { + let pieces = pieces(region.loops()); + let c = curves(plane, &pieces)?; + let start = match path.first() { + Some(crate::PathSeg::Line { a, .. } | crate::PathSeg::Arc { a, .. }) => *a, + None => return Err(KernelError::Invalid("the path needs 1…1000 segments".into())), + }; + // Each point of the profile turns with the path and the twist, and SMLib's skin is a cubic + // interpolation through the stations with free ends, whose error near the ends falls only + // as the square of the step: a square of reach 7.1 twisted a quarter turn over 40 mm is + // 2e-7 off in volume at 1.4 degree steps and 2e-6 at 2.8 (measured against V = A L). A + // station per degree of turn. + let turn: f64 = path.iter().map(|s| if let crate::PathSeg::Arc { angle, .. } = s { angle.abs() } else { 0.0 }).sum::() + ctl.twist.abs(); + let step = std::f64::consts::PI / 180.0; + let rail_points = ctl.rails.iter().map(Vec::len).max().unwrap_or(0); + let n = ((turn / step).ceil() as usize).max(rail_points.saturating_sub(1)).clamp(MIN_STATIONS, MAX_STATIONS); + let (pts, tangents) = path_samples(path, n)?; + let st = crate::sweep_path::controlled_stations(plane, &pts, &tangents, &[], ctl)?; + let (Some(first), Some(last)) = (st.first(), st.last()) else { return Err(KernelError::Invalid("sweep path".into())) }; + let s0 = first.at; + let size = pts.windows(2).map(|w| w[0].dist(w[1])).sum::() + region.outer.lp.signed_area().abs().sqrt(); + let closed = pts.last().is_some_and(|e| e.dist(start) <= 1e-7 * (1.0 + size)); + if closed { + for l in region.loops() { + crate::sweep_path::closing_offset(plane, l, &st)?; + } + } + let maps: Vec = st + .iter() + .map(|s| { + let [a, b, c] = s.m; + Affine { m: [p(a), p(b), p(c)], t: p(s.at - s.apply(s0)) } + }) + .collect(); + let brep = Brep::sweep_stations(&c, &maps, if closed { Ends::Joined } else { Ends::Capped }).map_err(|e| failed("sweep", e))?; + let body = body_of(&brep, "sweep")?; + let mid = st.get(st.len() / 2).copied().unwrap_or(*last); + let inside = seeds(plane, &pieces, |q| mid.at + mid.apply(q - s0)); + let roles = if closed { + // Round a closed path there are no caps, and the start holds seams only. + roles(&body, &plane.offset(size * 10.0), None, &[], &inside, false)? + } else { + let end = Plane::new(last.at + last.apply(plane.origin - s0), last.apply(plane.x), last.apply(plane.y)) + .ok_or_else(|| KernelError::Invalid("the sweep's end section is degenerate".into()))?; + roles(&body, plane, Some(&end), &seeds(plane, &pieces, |q| q), &inside, false)? + }; + with_history(body, &roles, &pieces) + } + /// A loft through `sections` (each a closed loop of pieces, all re-split to the same count /// and lined up), `first` the index of the section given first (its cap is the start cap). pub(crate) fn loft(sections: &[(Plane, Vec)], first: usize) -> Result { @@ -433,7 +548,13 @@ mod imp { 3 => 2, _ => 3, }; - let brep = Brep::loft(&curves, degree, true).map_err(|e| failed("loft", e))?; + loft_through(sections, first, &curves, degree) + } + + /// The loft skinned through `curves` (the sections', or with stations between them), `degree` + /// across them; `sections` and `first` as [`loft`] takes them name the faces. + fn loft_through(sections: &[(Plane, Vec)], first: usize, curves: &[Vec], degree: u32) -> Result { + let brep = Brep::loft(curves, degree, true).map_err(|e| failed("loft", e))?; let body = body_of(&brep, "loft")?; let (Some((start, lead)), Some(last)) = (sections.get(first), sections.get(if first == 0 { sections.len().saturating_sub(1) } else { 0 })) else { @@ -457,6 +578,281 @@ mod imp { with_history(body, &roles, &columns) } + /// A curve as a NURBS net: degree, control points (x, y, z) and weights, clamped knots. A line + /// is of degree 1; an arc is rational quadratic, one span per quarter turn or less. + struct Net { + degree: usize, + points: Vec, + weights: Vec, + knots: Vec, + } + + impl Net { + fn of(c: &Curve) -> Option { + let v = |q: [f64; 3]| Vec3::new(q[0], q[1], q[2]); + Some(match c { + Curve::Line { start, end } => { + Net { degree: 1, points: vec![v(*start), v(*end)], weights: vec![1.0; 2], knots: vec![0.0, 0.0, 1.0, 1.0] } + } + Curve::Arc { center, axis, reference, radius, start, end } => { + let (c, k, x) = (v(*center), v(*axis), v(*reference)); + let y = k.cross(x); + let sweep = end - start; + let m = (sweep / (std::f64::consts::PI / 2.0) - 1e-9).ceil().max(1.0) as usize; + let d = sweep / m as f64; + let at = |a: f64, r: f64| c + (x * a.cos() + y * a.sin()) * r; + let (mut points, mut weights, mut knots) = (vec![at(*start, *radius)], vec![1.0], vec![0.0; 3]); + for j in 0..m { + let a0 = start + d * j as f64; + points.push(at(a0 + d / 2.0, radius / (d / 2.0).cos())); + weights.push((d / 2.0).cos()); + points.push(at(a0 + d, *radius)); + weights.push(1.0); + let u = (j + 1) as f64 / m as f64; + knots.extend(if j + 1 < m { vec![u, u] } else { vec![1.0; 3] }); + } + Net { degree: 2, points, weights, knots } + } + Curve::Nurbs { degree, points, knots } => Net { + degree: *degree, + points: points.iter().map(|q| Vec3::new(q[0], q[1], q[2])).collect(), + weights: points.iter().map(|q| q[3]).collect(), + knots: knots.clone(), + }, + }) + } + + /// The same curve with other control points. + fn with(&self, points: Vec) -> Net { + Net { degree: self.degree, points, weights: self.weights.clone(), knots: self.knots.clone() } + } + + /// Whether `o` has the same degree, knots and weights (its control points pair up). + fn pairs_with(&self, o: &Net) -> bool { + let near = |a: &[f64], b: &[f64]| a.len() == b.len() && a.iter().zip(b).all(|(x, y)| (x - y).abs() <= 1e-9 * (1.0 + x.abs())); + self.degree == o.degree && self.points.len() == o.points.len() && near(&self.knots, &o.knots) && near(&self.weights, &o.weights) + } + + fn curve(&self) -> Curve { + let points = self.points.iter().zip(&self.weights).map(|(q, w)| [q.x, q.y, q.z, *w]).collect(); + Curve::Nurbs { degree: self.degree, points, knots: self.knots.clone() } + } + + /// The point `f` (0…1) of the way along the parameter range. + fn at(&self, f: f64) -> Option { + let ctrl = self.points.iter().zip(&self.weights).map(|(q, w)| [q.x * w, q.y * w, q.z * w, *w]).collect(); + let c = solidcraft_geom::bspline::CurveN::<4> { degree: self.degree, knots: self.knots.clone(), ctrl }; + let (t0, t1) = (*self.knots.first()?, *self.knots.last()?); + let ([x, y, z, w], _) = c.eval(t0 + (t1 - t0) * f.clamp(0.0, 1.0))?; + (w.abs() > 1e-300).then(|| Vec3::new(x / w, y / w, z / w)) + } + } + + /// The point of a section (its curves' nets) nearest `p`: as a wire parameter (the curve's + /// index plus how far along its parameter range), the point, and how far it is. + fn nearest_on(section: &[Net], p: Vec3) -> Option<(f64, Vec3, f64)> { + const SAMPLES: usize = 256; + let mut best: Option<(f64, Vec3, f64)> = None; + for (k, c) in section.iter().enumerate() { + let d = |f: f64| c.at(f).map_or(f64::INFINITY, |q| q.dist(p)); + let Some(i) = (0..=SAMPLES).min_by(|a, b| d(*a as f64 / SAMPLES as f64).total_cmp(&d(*b as f64 / SAMPLES as f64))) else { continue }; + // Golden-section search about the nearest sample. + let (mut lo, mut hi) = ((i as f64 - 1.0).max(0.0) / SAMPLES as f64, (i as f64 + 1.0).min(SAMPLES as f64) / SAMPLES as f64); + let g = 0.5 * (5f64.sqrt() - 1.0); + for _ in 0..80 { + let (a, b) = (hi - g * (hi - lo), lo + g * (hi - lo)); + if d(a) < d(b) { + hi = b; + } else { + lo = a; + } + } + let f = 0.5 * (lo + hi); + let q = c.at(f)?; + if best.as_ref().is_none_or(|b| q.dist(p) < b.2) { + best = Some((k as f64 + f, q, q.dist(p))); + } + } + best + } + + /// The point at wire parameter `g` of a section. + fn wire_point(section: &[Net], g: f64) -> Result { + let k = (g.floor().max(0.0) as usize).min(section.len().saturating_sub(1)); + section.get(k).and_then(|c| c.at(g - k as f64)).ok_or_else(|| KernelError::Failed("loft rails: section".into())) + } + + fn centroid(section: &[Net]) -> Vec3 { + let (sum, n) = section.iter().flat_map(|c| &c.points).fold((Vec3::ZERO, 0usize), |(s, n), q| (s + *q, n + 1)); + sum * (1.0 / n.max(1) as f64) + } + + /// A loft through `sections` (as [`loft`] takes them) guided by one or two `rails` (polylines + /// from the first section to the last, through every section's curve), by the law + /// [`crate::loft_rails`] gives a loft's rails: between each two sections stations are the + /// plain loft there (the sections' control points straight between two, the cubic through + /// them across more) moved so that the points the rails guide land on the rails. Each section + /// curve has to pair with the others' (the same kind, so their control points pair up). + pub(crate) fn loft_rails(sections: &[(Plane, Vec)], first: usize, rails: &[Vec], fit: crate::RailFit) -> Result { + if rails.is_empty() || rails.len() > 2 { + return Err(KernelError::Invalid("a loft takes one or two guide rails".into())); + } + let fail = |m: &str| KernelError::Failed(format!("loft rails: {m}")); + let count = sections.len(); + if count < 2 { + return Err(fail("sections")); + } + let curves: Vec> = sections.iter().map(|(pl, ps)| curves(pl, ps)).collect::>()?; + let nets: Vec> = curves + .iter() + .map(|cs| cs.iter().map(Net::of).collect::>>()) + .collect::>() + .ok_or_else(|| fail("section curves"))?; + let lead = nets.first().ok_or_else(|| fail("sections"))?; + let n = lead.len(); + if n == 0 || nets.iter().any(|s| s.len() != n || s.iter().zip(lead).any(|(a, b)| !a.pairs_with(b))) { + return Err(refuse("loft along guide rails between sections whose curves differ in kind (a line against an arc)")); + } + let mut bounds = solidcraft_geom::Aabb3::EMPTY; + for q in nets.iter().flatten().flat_map(|c| &c.points) { + bounds.add(*q); + } + let size = bounds.diagonal(); + let tol = (size * 1e-4).max(1e-6); + // The sections' parameters along the loft: in proportion to the mean distance between + // their curves' starts. + let mut along = vec![0.0]; + for w in nets.windows(2) { + let step = w[0].iter().zip(&w[1]).filter_map(|(a, b)| Some(a.points.first()?.dist(*b.points.first()?))).sum::() / n as f64; + along.push(along.last().copied().unwrap_or(0.0) + step); + } + let total = along.last().copied().unwrap_or(0.0); + if !(total > 1e-9) || along.windows(2).any(|w| !(w[1] - w[0] > 1e-9 * total)) { + return Err(KernelError::Invalid("two loft sections sit on top of each other".into())); + } + along.iter_mut().for_each(|u| *u /= total); + // The plain loft at v: each control point straight between two sections, the cubic (or + // quadratic) through them across more. + let degree_across = (count - 1).min(3); + let across: Vec>> = if count == 2 { + Vec::new() + } else { + (0..n) + .map(|k| { + (0..lead.get(k).map_or(0, |c| c.points.len())) + .map(|j| { + let pts: Vec<[f64; 3]> = nets.iter().filter_map(|s| s.get(k)?.points.get(j).map(|q| [q.x, q.y, q.z])).collect(); + solidcraft_geom::bspline::CurveN::interpolate(&along, &pts, degree_across).map_err(|e| fail(&e.to_string())) + }) + .collect::>>() + }) + .collect::>()? + }; + let plain = |v: f64| -> Result> { + let (Some(a), Some(b)) = (nets.first(), nets.last()) else { return Err(fail("sections")) }; + if count == 2 { + return Ok(a.iter().zip(b).map(|(ca, cb)| ca.with(ca.points.iter().zip(&cb.points).map(|(p, q)| p.lerp(*q, v)).collect())).collect()); + } + a.iter() + .zip(&across) + .map(|(c, cs)| { + let pts = cs.iter().map(|x| x.eval(v).map(|([px, py, pz], _)| Vec3::new(px, py, pz))).collect::>>(); + pts.map(|p| c.with(p)).ok_or_else(|| fail("the plain loft")) + }) + .collect() + }; + // The plain loft's middle and its direction there. + let middle = |v: f64| -> Result<(Vec3, Vec3)> { + let h = 1e-4; + let (a, b) = ((v - h).max(0.0), (v + h).min(1.0)); + let d = centroid(&plain(b)?) - centroid(&plain(a)?); + let d = d.normalized().ok_or_else(|| KernelError::Invalid("two loft sections sit on top of each other".into()))?; + Ok((centroid(&plain(v)?), d)) + }; + // Each rail from its end nearer the first section's middle, and where it meets each + // section: the guided point (wire parameter), that point, the nudge onto it from the + // crossing, and how far along the rail the crossing is. + let start = centroid(lead); + type Anchor = (f64, Vec3, Vec3, (usize, f64)); + let mut guided: Vec<(Vec, Vec)> = Vec::with_capacity(rails.len()); + for r in rails { + let mut pts = + crate::sweep_path::clean(r).map_err(|_| KernelError::Invalid("a guide rail needs 2…100000 finite points along its length".into()))?; + if let (Some(a), Some(b)) = (pts.first(), pts.last()) + && b.dist(start) < a.dist(start) + { + pts.reverse(); + } + let mut cursor = (0usize, 0.0f64); + let mut anchors = Vec::with_capacity(count); + for (i, ((plane, _), section)) in sections.iter().zip(&nets).enumerate() { + let Some((hit, at)) = crate::sweep_path::rail_crossing(&pts, plane.origin, plane.normal(), cursor, tol) else { + return Err(KernelError::Invalid(format!( + "a guide rail does not reach section {}: a rail has to run from the first section to the last, through every one", + i + 1 + ))); + }; + cursor = at; + let (g, on, d) = nearest_on(section, hit).ok_or_else(|| fail("section"))?; + if d > tol { + return Err(KernelError::Invalid(format!( + "a guide rail crosses section {}'s plane {d:.3} mm from its curves: a rail has to pass through every section", + i + 1 + ))); + } + anchors.push((g, on, on - hit, at)); + } + guided.push((pts, anchors)); + } + if let [a, b] = &guided[..] + && a.1.iter().zip(&b.1).any(|(p, q)| p.1.dist(q.1) <= tol) + { + return Err(KernelError::Invalid("the two guide rails meet a section at the same point".into())); + } + // The sections' normals, turned to agree with the first's. + let mut normals: Vec = Vec::with_capacity(count); + for (p, _) in sections { + let nn = p.normal(); + let prev = normals.last().copied().unwrap_or(nn); + normals.push(if nn.dot(prev) < 0.0 { nn * -1.0 } else { nn }); + } + let spans = count - 1; + let per = (crate::loft_rails::MAX_WIRES.saturating_sub(count) / spans).clamp(1, crate::loft_rails::MAX_PER_SPAN); + let mut all: Vec> = Vec::with_capacity(count + spans * per); + for i in 0..spans { + all.push(nets.get(i).ok_or_else(|| fail("sections"))?.iter().map(Net::curve).collect()); + let (Some(&v0), Some(&v1), Some(&n0), Some(&n1)) = (along.get(i), along.get(i + 1), normals.get(i), normals.get(i + 1)) else { + return Err(fail("sections")); + }; + let mut cursors: Vec<(usize, f64)> = guided.iter().map(|(_, a)| a.get(i).map_or((0, 0.0), |x| x.3)).collect(); + for q in 1..=per { + let f = q as f64 / (per + 1) as f64; + let v = v0 + f * (v1 - v0); + let station = plain(v)?; + let (mid, dir) = middle(v)?; + let normal = n0.lerp(n1, f).normalized().ok_or_else(|| fail("section planes"))?; + let mut pairs = Vec::with_capacity(guided.len()); + for ((pts, anchors), cursor) in guided.iter().zip(cursors.iter_mut()) { + let (Some(a0), Some(a1)) = (anchors.get(i), anchors.get(i + 1)) else { return Err(fail("rail")) }; + let anchor = wire_point(&station, crate::loft_rails::interpolate_g(a0.0, a1.0, f, n, true))?; + let Some((hit, at)) = crate::sweep_path::rail_crossing(pts, anchor, normal, *cursor, tol) else { + return Err(KernelError::Invalid(format!( + "a guide rail turns back between sections {} and {}: it has to run on from one section to the next", + i + 1, + i + 2 + ))); + }; + *cursor = at; + pairs.push((anchor, hit + a0.2 * (1.0 - f) + a1.2 * f)); + } + let map = crate::loft_rails::station_map(&pairs, mid, dir, fit, size)?; + all.push(station.iter().map(|c| c.with(c.points.iter().map(|q| map.apply(*q)).collect()).curve()).collect()); + } + } + all.push(nets.last().ok_or_else(|| fail("sections"))?.iter().map(Net::curve).collect()); + loft_through(sections, first, &all, 3) + } + /// The profile `pieces` on `plane` carried by a screw motion about the axis through `origin` /// along `axis`: `turns` right-handed turns, rising `pitch` per turn. pub(crate) fn screw(plane: &Plane, pieces: &[Piece], origin: Vec3, axis: Vec3, pitch: f64, turns: f64) -> Result { diff --git a/crates/kernel/src/smlib_sweeps_tests.rs b/crates/kernel/src/smlib_sweeps_tests.rs index 52bb97467d..f0489bb023 100644 --- a/crates/kernel/src/smlib_sweeps_tests.rs +++ b/crates/kernel/src/smlib_sweeps_tests.rs @@ -117,3 +117,50 @@ fn a_sweep_keeps_a_profile_off_its_path_where_it_lies() { let b = smlib_sweeps::sweep(&plane, &pieces_of(&disc), &path, false).unwrap(); close(volume(&b), PI * 4.0 * (PI / 2.0 * 25.0), 1e-3, "bent tube"); } + +fn numbered(r: &Region2) -> crate::numbered::NumberedRegion { + number_regions(std::slice::from_ref(r)).pop().unwrap().oriented() +} + +#[test] +fn a_sweep_twists_scales_and_follows_rails_through_stations() { + let square = region(Loop2::polygon(&[Vec2::new(-5.0, -5.0), Vec2::new(5.0, -5.0), Vec2::new(5.0, 5.0), Vec2::new(-5.0, 5.0)])); + let up = [crate::PathSeg::Line { a: Vec3::ZERO, b: Vec3::new(0.0, 0.0, 40.0) }]; + // Turning about the section's centre keeps its area: V = 100 * 40, the caps and the four + // walls named. + let ctl = crate::SweepControl { twist: PI / 2.0, ..Default::default() }; + let b = smlib_sweeps::sweep_controlled(&Plane::XY, &numbered(&square), &up, &ctl).unwrap(); + close(volume(&b), 4000.0, 1e-3, "twisted square"); + let o = origins(&b); + assert!(o.contains(&cap(false)) && o.contains(&cap(true)) && (0..4).all(|k| o.contains(&profile(k))), "{o:?}"); + // Scaled from 5 to 10 over 30, or between rails that widen as much: a cone frustum. + let disc = region(Loop2::circle(Vec2::ZERO, 5.0)); + let frustum = |h: f64| h / 3.0 * PI * (25.0 + 100.0 + 50.0); + let short = [crate::PathSeg::Line { a: Vec3::ZERO, b: Vec3::new(0.0, 0.0, 30.0) }]; + let ctl = crate::SweepControl { scale: 2.0, ..Default::default() }; + close(volume(&smlib_sweeps::sweep_controlled(&Plane::XY, &numbered(&disc), &short, &ctl).unwrap()), frustum(30.0), 1e-3, "end scale"); + let rails = [vec![Vec3::new(5.0, 0.0, 0.0), Vec3::new(10.0, 0.0, 30.0)], vec![Vec3::new(-5.0, 0.0, 0.0), Vec3::new(-10.0, 0.0, 30.0)]]; + for n in 1..=2 { + let ctl = crate::SweepControl { rails: rails[..n].to_vec(), ..Default::default() }; + close(volume(&smlib_sweeps::sweep_controlled(&Plane::XY, &numbered(&disc), &short, &ctl).unwrap()), frustum(30.0), 1e-3, "rails"); + } +} + +#[test] +fn a_twisted_sweep_round_a_closed_path_joins_its_ends() { + // An 8 x 8 square round a circle of radius 30, a quarter turn of twist (onto itself): Pappus. + let xz = Plane { origin: Vec3::ZERO, x: Vec3::X, y: Vec3::Z }; + let square = region(Loop2::polygon(&[Vec2::new(26.0, -4.0), Vec2::new(34.0, -4.0), Vec2::new(34.0, 4.0), Vec2::new(26.0, 4.0)])); + let a = Vec3::new(30.0, 0.0, 0.0); + let round = [ + crate::PathSeg::Arc { a, center: Vec3::ZERO, axis: Vec3::Z, angle: PI }, + crate::PathSeg::Arc { a: -a, center: Vec3::ZERO, axis: Vec3::Z, angle: PI }, + ]; + let ctl = crate::SweepControl { twist: PI / 2.0, ..Default::default() }; + let b = smlib_sweeps::sweep_controlled(&xz, &numbered(&square), &round, &ctl).unwrap(); + close(volume(&b), 64.0 * TAU * 30.0, 1e-3, "twisted ring"); + assert!(origins(&b).iter().all(|o| o.contains("Profile")), "{:?}", origins(&b)); + // A twist that does not bring the square back onto itself cannot close. + let ctl = crate::SweepControl { twist: PI / 3.0, ..Default::default() }; + assert!(smlib_sweeps::sweep_controlled(&xz, &numbered(&square), &round, &ctl).is_err()); +} diff --git a/crates/kernel/src/sweep_path.rs b/crates/kernel/src/sweep_path.rs index b0720affc9..53490c4af1 100644 --- a/crates/kernel/src/sweep_path.rs +++ b/crates/kernel/src/sweep_path.rs @@ -104,7 +104,7 @@ impl SweepControl { } /// Whether the profile changes along the path (it twists, scales or follows rails). - fn varies(&self) -> bool { + pub(crate) fn varies(&self) -> bool { self.twist != 0.0 || self.scale != 1.0 || !self.rails.is_empty() } @@ -141,13 +141,13 @@ impl SweepControl { /// A linear map of world vectors (columns) plus the station it is centred at. #[derive(Clone, Copy)] -struct Station { - at: Vec3, - m: [Vec3; 3], +pub(crate) struct Station { + pub at: Vec3, + pub m: [Vec3; 3], } impl Station { - fn apply(&self, v: Vec3) -> Vec3 { + pub(crate) fn apply(&self, v: Vec3) -> Vec3 { self.m[0] * v.x + self.m[1] * v.y + self.m[2] * v.z } } @@ -177,7 +177,7 @@ pub(crate) fn closed_path(pts: &[Vec3]) -> bool { /// The cyclic correspondence of a loop's last section to its first. Matching whole curves /// (rather than vertices alone) prevents sewing an incompatible twist or scale shut. -fn closing_offset(plane: &Plane, profile: &NumberedLoop, st: &[Station]) -> Result { +pub(crate) fn closing_offset(plane: &Plane, profile: &NumberedLoop, st: &[Station]) -> Result { let fail = || KernelError::Invalid("a sweep along a closed path must return to its starting profile; adjust the twist, end scale or rails".into()); let (Some(first), Some(last)) = (st.first(), st.last()) else { return Err(fail()) }; @@ -301,7 +301,13 @@ fn frames_toward(tangents: &[Vec3], d: Vec3) -> Result> { /// at its own (twisted and scaled about the profile's normal through the first station, then /// turned by the orientation's frames), optionally stretched across a corner (mitre). Twist and /// scale grow with the length run along `pts`. -fn controlled_stations(plane: &Plane, pts: &[Vec3], tangents: &[Vec3], mitre: &[Option<(Vec3, f64)>], ctl: &SweepControl) -> Result> { +pub(crate) fn controlled_stations( + plane: &Plane, + pts: &[Vec3], + tangents: &[Vec3], + mitre: &[Option<(Vec3, f64)>], + ctl: &SweepControl, +) -> Result> { if !ctl.rails.is_empty() { return railed_stations(pts, tangents, mitre, ctl); } @@ -524,7 +530,7 @@ pub fn sweep_path(plane: &Plane, region: &Region2, path: &[Vec3]) -> Result Result { ctl.check()?; if crate::backend()? == crate::Backend::Smlib { - return Err(crate::smlib_sweeps::refuse("sweep along a path of points (a 3D sketch curve)")); + return smlib_sweep_points(plane, region, path, ctl); } if region.area().abs() < 1e-9 { return Err(KernelError::Invalid("profile has no area".into())); @@ -560,6 +566,24 @@ pub fn sweep_path_with(plane: &Plane, region: &Region2, path: &[Vec3], ctl: &Swe ruled(plane, &profile, &pts, ctl) } +/// [`sweep_path_with`] through SMLib: a path of points is the polyline through them, its +/// straight runs the lines between them. A line, or a polyline with corners, is swept as those +/// lines ([`crate::sweep_with`], which takes a path with corners only where SMLib has a sweep for +/// it); a smooth polyline of more points stands for a curve (a drawn 3D spline) that SMLib has no +/// path for yet. +fn smlib_sweep_points(plane: &Plane, region: &Region2, path: &[Vec3], ctl: &SweepControl) -> Result { + let pts = clean(path)?; + let smooth = pts.len() > 2 && !pts.windows(3).any(|w| turn(w[1] - w[0], w[2] - w[1]) > CORNER); + if smooth { + return Err(crate::smlib_sweeps::refuse("sweep along a smooth path of points (a 3D sketch spline)")); + } + let segs: Vec = pts.windows(2).map(|w| crate::PathSeg::Line { a: w[0], b: w[1] }).collect(); + if segs.len() > 1000 { + return Err(KernelError::Invalid("the path needs 1…1000 segments".into())); + } + crate::sweep_with(plane, region, &segs, ctl) +} + fn central_tangents(p: &[Vec3]) -> Result> { let n = p.len(); let closed = closed_path(p);