kernel: helical sweeps and lofts to a point on SMLib only

sweep_helix is SMLib's screw sweep (smlib_screw_sweep) and loft_to_point SMLib's
loft to a point (smlib_loft_to_point), with their checks and history as before.
Deleted: the B-spline walls fitted through sampled stations (smooth_helix), the
ruled walls between stations and the face-by-face pyramid, and the Pappus
estimate only they used.

(not built: base in progress)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt
This commit is contained in:
jyviko
2026-10-22 00:33:33 -04:00
parent 2d3d0a081b
commit 2b77f09f86
+9 -152
View File
@@ -1,14 +1,9 @@
//! Helical sweeps (coils, springs) and lofts to a point (pyramids, cones).
use mt::builder;
use solidcraft_geom::{Grid, Loop2, Mesh, Plane, Region2, Seg2, Vec2, Vec3};
use truck_modeling as mt;
use solidcraft_geom::{Loop2, Mesh, Plane, Region2, Seg2, Vec2, Vec3};
use crate::body::{Body, Solid, p3};
use crate::build::{NetFit, cap_and_close, interpolate_net, uniform_wire};
use crate::numbered::NumberedLoop;
use crate::record::{self, Side};
use crate::{KernelError, Result, guard};
use crate::body::Body;
use crate::{KernelError, Result};
/// Rotate `v` about the unit axis `k` by `ang` (Rodrigues).
fn rot(v: Vec3, k: Vec3, ang: f64) -> Vec3 {
@@ -17,9 +12,8 @@ fn rot(v: Vec3, k: Vec3, ang: f64) -> Vec3 {
/// Sweep a planar region's outer loop along a helix about `axis` through `axis_origin`:
/// `turns` revolutions (right-handed about `axis`), rising `pitch` per turn. The profile keeps
/// its attitude relative to the axis, as a screw motion does. The walls are B-spline surfaces
/// (one per profile segment); if those cannot be built, ruled walls between stations 3° apart
/// stand in.
/// its attitude relative to the axis, as a screw motion does: SMLib's screw sweep, one wall per
/// profile segment.
pub fn sweep_helix(plane: &Plane, region: &Region2, axis_origin: Vec3, axis: Vec3, pitch: f64, turns: f64) -> Result<Body> {
let k = axis.normalized().ok_or_else(|| KernelError::Invalid("helix axis".into()))?;
if !(pitch.is_finite() && turns.is_finite() && turns > 0.0 && turns <= 200.0 && pitch.abs() < 1e7 && axis_origin.is_finite()) {
@@ -38,34 +32,7 @@ pub fn sweep_helix(plane: &Plane, region: &Region2, axis_origin: Vec3, axis: Vec
if turns > 1.0 && pitch.abs() <= extent {
return Err(KernelError::Invalid(format!("the coil runs into itself: the pitch must exceed the section ({extent:.3} mm)")));
}
if crate::backend()? == crate::Backend::Smlib {
return crate::smlib_sweeps::screw(plane, &crate::smlib_sweeps::pieces([&profile]), axis_origin, k, pitch, turns);
}
// One smooth B-spline wall per profile segment; the ruled walls below are the fallback.
if let Ok(b) = smooth_helix(plane, &profile, axis_origin, k, pitch, turns) {
return Ok(b);
}
let n = ((turns * 120.0).ceil() as usize).clamp(4, 24_000);
let total = std::f64::consts::TAU * turns;
let mut wires = Vec::with_capacity(n + 1);
for i in 0..=n {
let a = total * i as f64 / n as f64;
let lift = k * (pitch * a / std::f64::consts::TAU);
let o = axis_origin + rot(plane.origin - axis_origin, k, a) + lift;
let pl = Plane { origin: o, x: rot(plane.x, k, a), y: rot(plane.y, k, a) };
wires.push(uniform_wire(&pl, lp)?);
}
guard("helix sweep", || {
// Every station's wire has the profile's segments, edge for edge.
let mut faces: Vec<Side> = Vec::new();
for w in wires.windows(2) {
let (Some(a), Some(b)) = (w.first(), w.get(1)) else { continue };
let sides = builder::try_wire_homotopy(a, b).map_err(|e| KernelError::Failed(format!("helix: {e}")))?;
faces.extend(record::ruled(&sides, a, b, &profile.from, &profile.from));
}
let (Some(first), Some(last)) = (wires.first(), wires.last()) else { return Err(KernelError::Failed("helix".into())) };
cap_and_close(faces, std::slice::from_ref(first), std::slice::from_ref(last), false)
})
crate::smlib_sweeps::screw(plane, &crate::smlib_sweeps::pieces([&profile]), axis_origin, k, pitch, turns)
}
/// A screw sweep of a convex polygon (in `plane`'s coordinates) as a closed triangle mesh,
@@ -107,84 +74,8 @@ pub fn sweep_helix_mesh(plane: &Plane, polygon: &[Vec2], axis_origin: Vec3, axis
Ok(mesh)
}
/// Samples per profile segment and stations per turn of the smooth helical walls.
const SEG_SAMPLES: usize = 12;
const STATIONS_PER_TURN: f64 = 12.0;
/// Helical walls as B-spline surfaces: each profile segment, sampled at every station of the
/// screw motion, interpolated by cubics across the profile and then along the helix.
fn smooth_helix(plane: &Plane, profile: &NumberedLoop, o: Vec3, k: Vec3, pitch: f64, turns: f64) -> Result<Body> {
let lp = &profile.lp;
let nseg = lp.segs.len();
if nseg == 0 || nseg > 64 {
return Err(KernelError::Invalid("helix profile".into()));
}
let n = ((turns * STATIONS_PER_TURN).ceil() as usize).clamp(4, 4800);
let total = std::f64::consts::TAU * turns;
let at = |i: usize, q: solidcraft_geom::Vec2| {
let a = total * i as f64 / n as f64;
let w = plane.to_world(q);
o + rot(w - o, k, a) + k * (pitch * a / std::f64::consts::TAU)
};
let fail = |w: &str| KernelError::Failed(format!("smooth helix: {w}"));
// Per segment: the wall surface through its samples at every station.
let mut faces_data = Vec::with_capacity(nseg);
for s in &lp.segs {
// u along the helix (the stations), v along the profile segment.
let samples =
Grid::from_fn([n + 1, SEG_SAMPLES], |i, j| at(i, s.point_at(j as f64 / (SEG_SAMPLES - 1) as f64))).ok_or_else(|| fail("samples"))?;
faces_data.push(interpolate_net(&samples).ok_or_else(|| fail("wall surface"))?);
}
guard("smooth helix", || {
let start: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(at(0, s.start())))).collect();
let end: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(at(n, s.start())))).collect();
let mut e0 = Vec::with_capacity(nseg);
let mut e1 = Vec::with_capacity(nseg);
let mut seams = Vec::with_capacity(nseg);
for (j, NetFit { u_ends: [c0, c1], v_ends: [seam, _], .. }) in faces_data.iter().enumerate() {
let (Some(a0), Some(b0), Some(a1), Some(b1)) = (start.get(j), start.get((j + 1) % nseg), end.get(j), end.get((j + 1) % nseg)) else {
return Err(fail("vertices"));
};
e0.push(mt::Edge::try_new(a0, b0, mt::Curve::BSplineCurve(c0.clone())).map_err(|e| fail(&e.to_string()))?);
e1.push(mt::Edge::try_new(a1, b1, mt::Curve::BSplineCurve(c1.clone())).map_err(|e| fail(&e.to_string()))?);
seams.push(mt::Edge::try_new(a0, a1, mt::Curve::BSplineCurve(seam.clone())).map_err(|e| fail(&e.to_string()))?);
}
let mut walls = Vec::with_capacity(nseg);
for (j, NetFit { surface: surf, .. }) in faces_data.iter().enumerate() {
let (Some(a), Some(b), Some(sa), Some(sb)) = (e0.get(j), e1.get(j), seams.get(j), seams.get((j + 1) % nseg)) else {
return Err(fail("edges"));
};
// Counter-clockwise in the surface's (u along the helix, v along the profile) space.
let w: mt::Wire = vec![sa.clone(), b.clone(), sb.inverse(), a.inverse()].into();
let wall = mt::Face::try_new(vec![w], mt::Surface::BSplineSurface(surf.clone())).map_err(|e| fail(&e.to_string()))?;
walls.push((wall, Some(record::profile(profile.numbers(j)))));
}
let w0: mt::Wire = e0.iter().cloned().collect();
let w1: mt::Wire = e1.iter().cloned().collect();
// Either orientation of the walls; keep the one that measures as a proper solid.
let expected = lp_area_path(plane, lp, o, k, turns);
// A coarse mesh (within ~10%) is enough to tell a proper solid from an empty or inside-out one.
let section_tol = (lp.signed_area().abs().sqrt() * 0.05).max(1e-3);
for ws in [walls.clone(), walls.iter().map(|(f, prov)| (f.inverse(), prov.clone())).collect()] {
if let Ok(b) = cap_and_close(ws, std::slice::from_ref(&w0), std::slice::from_ref(&w1), false)
&& b.mesh(&crate::MeshTolerance::chord(section_tol)).is_ok_and(|m| (m.measure().volume - expected).abs() <= 0.15 * expected.abs())
{
return Ok(b);
}
}
Err(fail("the walls do not close into the expected solid"))
})
}
/// Pappus estimate of a screw sweep's volume: section area × the centroid's path around the axis.
fn lp_area_path(plane: &Plane, lp: &Loop2, o: Vec3, k: Vec3, turns: f64) -> f64 {
let c = plane.to_world(lp.centroid());
let d = c - o;
let r = (d - k * d.dot(k)).len();
lp.signed_area().abs() * std::f64::consts::TAU * r * turns
}
/// A solid from a planar loop of straight segments to an apex point off its plane (a pyramid).
/// A solid from a planar loop of straight segments to an apex point off its plane (a pyramid):
/// SMLib's loft to a point.
pub fn loft_to_point(plane: &Plane, lp: &Loop2, apex: Vec3) -> Result<Body> {
if !apex.is_finite() || lp.segs.is_empty() || lp.segs.len() > 10_000 {
return Err(KernelError::Invalid("loft to a point: bad section or apex".into()));
@@ -198,46 +89,12 @@ pub fn loft_to_point(plane: &Plane, lp: &Loop2, apex: Vec3) -> Result<Body> {
// Wind the base so its normal points away from the apex.
let above = (apex - plane.origin).dot(plane.normal()) > 0.0;
let base = crate::numbered::number_loops([lp]).pop().ok_or_else(|| KernelError::Invalid("loft to a point: bad section".into()))?.oriented(above);
if crate::backend()? == crate::Backend::Smlib {
return crate::smlib_sweeps::loft_to_point(plane, &crate::smlib_sweeps::pieces([&base]), apex);
}
let lp = &base.lp;
guard("loft to point", || {
let verts: Vec<mt::Vertex> = lp.segs.iter().map(|s| builder::vertex(p3(plane.to_world(s.start())))).collect();
let top = builder::vertex(p3(apex));
let m = verts.len();
let edges: Vec<mt::Edge> = (0..m).filter_map(|i| Some(builder::line(verts.get(i)?, verts.get((i + 1) % m)?))).collect();
let sides: Vec<mt::Edge> = verts.iter().map(|v| builder::line(v, &top)).collect();
let mut faces = Vec::with_capacity(m + 1);
// History as the faces are made: the base is the section (the start cap), each
// triangle comes from its base segment.
let mut history = Vec::with_capacity(m + 1);
let bw: mt::Wire = edges.iter().cloned().collect();
faces.push(builder::try_attach_plane(&[bw]).map_err(|e| KernelError::Failed(format!("loft base: {e}")))?);
history.push(record::cap(false));
for i in 0..m {
let (Some(e), Some(up), Some(down)) = (edges.get(i), sides.get((i + 1) % m), sides.get(i)) else { continue };
// Triangle i: along the base edge, up to the apex, back down — opposite to the base's winding.
let w: mt::Wire = vec![e.inverse(), down.clone(), up.inverse()].into();
faces.push(builder::try_attach_plane(&[w]).map_err(|e| KernelError::Failed(format!("loft side: {e}")))?);
history.push(record::profile(base.numbers(i)));
}
let shell: mt::Shell = faces.into();
let body = match Solid::try_new(vec![shell.clone()]) {
Ok(s) => Body::new(s)?,
Err(_) => {
let inv: mt::Shell = shell.face_iter().map(|f| f.inverse()).collect();
Body::new(Solid::try_new(vec![inv]).map_err(|e| KernelError::Failed(format!("loft to point: {e}")))?)?
}
};
record::attach(body, &[], Some(history), Default::default())
})
crate::smlib_sweeps::loft_to_point(plane, &crate::smlib_sweeps::pieces([&base]), apex)
}
#[cfg(test)]
mod tests {
use super::*;
use solidcraft_geom::Vec2;
#[test]
fn coil_volume_is_section_times_helix_length() {