From 7225d3f2b696d2eb2c61cd149eb2cd9507bc884e Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 22 Oct 2026 00:53:26 -0400 Subject: [PATCH] kernel: Split Face and Split Body on SMLib, the tool as exact record surface patches A tool is now a set of schema 4 record surfaces, each over a parameter box, unbounded where the tool is: a plane; a face's surface over its own box; a sketch curve as one patch per segment swept along the direction (a line a plane strip, an arc swept square to its sketch a cylinder strip, any other segment its exact rational Bezier extruded), an open curve continued along its end tangents past the body. SMLib cuts the faces where the patches meet them (Brep::split_faces: its intersection curves trimmed to each face by MergeCurvesOnSurface) and slices a body with one patch (Brep::split_body); the history is SMLib's own tags. Nothing is traced on a mesh, sampled or fitted. Which side of the tool a point is on is its offset along the normal of the nearest patch, each patch's surface prepared once with its projection grid (ToolSide), never per point. ToolSide::of and side_of_point take the body the tool must reach past; face_side becomes face_sides (every face of a body, one mesh); analytic_field becomes analytic_face (a face of a body, on its record); split_any returns the face map instead of the truck imprint. split_edges_at and fit (a Catmull-Rom fit through samples) are deleted with the truck code. ssi::Surface::from_record turns a record surface over a box into SMLib's. (not built: base in progress) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_017itGh5fpMEousnZ91cYcVt --- crates/kernel/src/splitface.rs | 1943 +++++++++----------------------- crates/smlib/src/ssi.rs | 53 + 2 files changed, 588 insertions(+), 1408 deletions(-) diff --git a/crates/kernel/src/splitface.rs b/crates/kernel/src/splitface.rs index e511edfc1f..89bbb6c421 100644 --- a/crates/kernel/src/splitface.rs +++ b/crates/kernel/src/splitface.rs @@ -1,21 +1,32 @@ -//! Split Face: divide chosen faces of a body along where a tool meets them (a plane, another -//! face's surface, or sketch curves swept along the sketch normal), leaving the solid itself -//! unchanged. Edges the split line ends on are cut in every face that uses them. +//! Split Face and Split Body: divide chosen faces of a body along where a tool meets them (a +//! plane, another face's surface, or sketch curves swept along a direction), leaving the solid +//! unchanged, or cut a body in two along the tool. Edges the split line ends on are cut in every +//! face that uses them. //! -//! The tool is a signed field (which side of the tool a point is on). Its zero set on a face is -//! traced on the face's triangulation, every crossing is placed exactly on the surface by -//! bisection in the face's parameters, and the line ends exactly on the boundary edges. +//! A tool is a set of exact surface patches: schema 4 record surfaces, each over a parameter box, +//! unbounded where the tool is (a plane, a face's surface extended). A plane is one patch, so is a +//! face's surface; a sketch curve is one patch per segment swept along the direction (a line a +//! plane strip, an arc swept square to its sketch a cylinder strip, any other segment its exact +//! rational Bézier extruded), an open curve's ends continued by lines past the body. SMLib cuts +//! the faces where the patches meet them ([`solidcraft_smlib::Brep::split_faces`]), its own +//! intersection curves trimmed to each face, and slices a body with one patch +//! ([`solidcraft_smlib::Brep::split_body`]); the history is SMLib's own. +//! +//! Which side of a tool a point is on is its offset along the normal of the nearest patch +//! ([`ToolSide`]): each patch's surface is prepared once (its projection grid kept), never per +//! point. Every patch's normal points to the tool's positive side: a plane's normal, a face's +//! outward side, the left of a sketch curve (inside a closed one, counter-clockwise). -use std::collections::HashMap; +use std::f64::consts::TAU; -use mt::{BoundedCurve, InnerSpace, MetricSpace, ParametricCurve, ParametricSurface, ParametricSurface3D, SearchNearestParameter}; -use solidcraft_geom::{Plane, Seg2, UnitVec3, Vec2, Vec3}; -use truck_meshalgo::tessellation::RobustMeshableShape; -use truck_modeling as mt; +use solidcraft_geom::brep::eval::{CurveEval, SurfaceEval}; +use solidcraft_geom::brep::{self, Frame, Surface}; +use solidcraft_geom::tess::uv::Starts; +use solidcraft_geom::{Aabb3, Plane, Seg2, UnitVec3, Vec2, Vec3}; -use crate::body::{Body, from_p3, p3}; -use crate::record::Rebuild; -use crate::{KernelError, Result, guard}; +use crate::body::Body; +use crate::history::{EdgeOrigin, FaceOrigin, History}; +use crate::{KernelError, Result, guard, record}; /// What splits the faces. #[derive(Clone, Debug)] @@ -25,1111 +36,491 @@ pub enum SplitTool<'a> { /// The surface of another body's face (`face` in that body's face order), extended. Face { body: &'a Body, face: usize }, /// Sketch curves (chains of segments in `plane`) swept along `along` (none: the plane's - /// normal); each curve splits on its own, an open one extended past its ends. A single line - /// splits exactly (a plane), and so does a single arc swept along the normal (a cylinder); - /// longer chains, splines and arcs swept aslant follow the curve. + /// normal); each curve splits on its own, an open one extended past its ends. Curves { plane: Plane, curves: Vec>, along: Option }, } -/// Which side of the tool a point is on (zero on the tool). -enum Field { - Plane { - o: Vec3, - n: Vec3, - }, - Surface(mt::Surface), - /// A circular cylinder (radius `r` about the line through `o` along unit `axis`); `out` is - /// +1 when the tool's normal points away from the axis. - Cylinder { - o: Vec3, - axis: UnitVec3, - r: f64, - out: f64, - }, - /// A chain of points in `plane`, swept along `along` (the plane's normal when none). A - /// closed chain is positive inside it; an open one on its left, its end segments extended. - Curve { - plane: Plane, - pts: Vec, - closed: bool, - along: Option, - }, +/// Most segments a sketch tool may hold. +const MAX_SEGMENTS: usize = 4096; + +fn invalid(what: impl Into) -> KernelError { + KernelError::Invalid(what.into()) } -impl Field { - fn eval(&self, p: Vec3) -> f64 { - match self { - Field::Plane { o, n } => (p - *o).dot(*n), - Field::Surface(s) => match s.search_nearest_parameter(p3(p), None, 50) { - Some((u, v)) => { - let q = from_p3(s.subs(u, v)); - let n = s.normal(u, v); - (p - q).dot(Vec3::new(n.x, n.y, n.z)) - } - None => f64::NAN, - }, - Field::Cylinder { o, axis, r, out } => { - let v = p - *o; - let radial = v - *axis * v.dot(*axis); - *out * (radial.len() - *r) - } - Field::Curve { plane, pts, closed, along } => { - // Where p's line along the sweep direction meets the plane. - let p = match along { - Some(d) => p - *d * ((p - plane.origin).dot(plane.normal()) / d.dot(plane.normal())), - None => p, - }; - let q = plane.to_local(p); - let q = Vec2::new(q.x, q.y); - if *closed { - // Distance to the loop, positive inside it. - let d = (0..pts.len().saturating_sub(1)) - .filter_map(|i| Some((*pts.get(i)?, *pts.get(i + 1)?))) - .map(|(a, b)| { - let d = b - a; - let l2 = d.dot(d); - let t = if l2 > 0.0 { ((q - a).dot(d) / l2).clamp(0.0, 1.0) } else { 0.0 }; - (q - (a + d * t)).len() - }) - .fold(f64::INFINITY, f64::min); - let uv: Vec<(f64, f64)> = pts.iter().map(|p| (p.x, p.y)).collect(); - return if inside(&uv, (q.x, q.y)) { d } else { -d }; - } - // Signed distance to the polyline, its end segments extended as lines. - let n = pts.len(); - let mut best = (f64::INFINITY, 0.0); - for i in 0..n.saturating_sub(1) { - let (Some(&a), Some(&b)) = (pts.get(i), pts.get(i + 1)) else { continue }; - let d = b - a; - let l2 = d.dot(d); - if !(l2 > 0.0) { - continue; - } - let mut t = (q - a).dot(d) / l2; - if i > 0 { - t = t.max(0.0); - } - if i + 2 < n { - t = t.min(1.0); - } - let c = a + d * t; - let dist = (q - c).len(); - if dist < best.0 { - let side = d.x * (q.y - a.y) - d.y * (q.x - a.x); - best = (dist, if side >= 0.0 { dist } else { -dist }); - } - } - best.1 - } - } - } +fn failed(what: &str, e: impl std::fmt::Display) -> KernelError { + KernelError::Failed(format!("SMLib {what}: {e}")) } -/// A crossing of the split line with a face boundary edge. -#[derive(Clone)] -struct Hit { - edge: mt::Edge, - t: f64, - p: mt::Point3, +fn a3(v: Vec3) -> [f64; 3] { + [v.x, v.y, v.z] } -/// One traced split line on a face: points (exact) with their surface parameters, and whether -/// it closes on itself. -struct Chain { - pts: Vec<(mt::Point3, (f64, f64))>, - closed: bool, +fn v3(p: [f64; 3]) -> Vec3 { + Vec3::new(p[0], p[1], p[2]) } -/// Mesh tolerance for tracing a face. -fn trace_tol(face: &mt::Face, size: f64) -> f64 { - let _ = face; - trace_tol_of(size) +/// One piece of a tool: a record surface over `domain` (`±∞` where the tool runs on without +/// end), `flip` when the tool's positive side is against the surface's normal. +#[derive(Clone, Debug)] +struct Patch { + surface: Surface, + domain: [[f64; 2]; 2], + flip: bool, } -/// Mesh tolerance for tracing on a body of this size. -fn trace_tol_of(size: f64) -> f64 { - (size * 2e-3).max(1e-4) -} +const ALL: [f64; 2] = [f64::NEG_INFINITY, f64::INFINITY]; -/// Trace the zero set of `g` over a face: chains of exact points. -fn trace(face: &mt::Face, g: &Field, size: f64) -> Result> { - if let Some(chains) = trace_flat(face, g, size) { - return Ok(chains); +impl Patch { + fn frame(origin: Vec3, x: Vec3, y: Vec3) -> Frame { + Frame { origin: a3(origin), x_axis: a3(x), y_axis: a3(y) } } - let tol = trace_tol(face, size); - let shell: mt::Shell = vec![face.clone()].into(); - let meshed = shell.robust_triangulation(tol); - let Some(mf) = meshed.face_iter().next() else { return Ok(Vec::new()) }; - let Some(pm) = mf.surface() else { return Err(KernelError::Failed("the face could not be meshed".into())) }; - let surface = face.surface(); - let (pos, uvs) = (pm.positions(), pm.uv_coords()); - // The face's triangles subdivided in parameter space, finely enough that a small split line - // (a hole's circle on a big face) crosses some triangle sides. Corners shared between - // triangles are found by their (rounded) parameters. - let step = size / 100.0; - let mut ids: HashMap<(i64, i64), usize> = HashMap::new(); - let mut corner_uv: Vec<(f64, f64)> = Vec::new(); - let (mut umin, mut umax, mut vmin, mut vmax) = (f64::INFINITY, f64::NEG_INFINITY, f64::INFINITY, f64::NEG_INFINITY); - for uv in uvs { - umin = umin.min(uv.x); - umax = umax.max(uv.x); - vmin = vmin.min(uv.y); - vmax = vmax.max(uv.y); - } - let qu = ((umax - umin).abs() * 1e-9).max(1e-12); - let qv = ((vmax - vmin).abs() * 1e-9).max(1e-12); - let mut corner = |uv: (f64, f64)| -> usize { - let k = ((uv.0 / qu).round() as i64, (uv.1 / qv).round() as i64); - *ids.entry(k).or_insert_with(|| { - corner_uv.push(uv); - corner_uv.len() - 1 - }) - }; - let mut base: Vec<[((f64, f64), mt::Point3); 3]> = Vec::new(); - for t in pm.faces().triangle_iter() { - let mut c = [((0.0, 0.0), mt::Point3::new(0.0, 0.0, 0.0)); 3]; - for (k, v) in t.iter().enumerate() { - let Some(uv) = v.uv.and_then(|i| uvs.get(i)) else { return Err(KernelError::Failed("mesh without parameters".into())) }; - let p = pos.get(v.pos).copied().unwrap_or(mt::Point3::new(0.0, 0.0, 0.0)); - c[k] = ((uv.x, uv.y), p); - } - base.push(c); - } - // One division for every triangle, so shared sides divide alike. - let longest = base.iter().flat_map(|c| (0..3).map(move |k| c[k].1.distance(c[(k + 1) % 3].1))).fold(0.0, f64::max); - let n = ((longest / step).ceil() as usize).clamp(1, 48); - let mut tris: Vec<[(usize, (f64, f64)); 3]> = Vec::new(); - for c in base { - let at = |i: usize, j: usize| { - let (a, b, cc) = (c[0].0, c[1].0, c[2].0); - let (s, t) = (i as f64 / n as f64, j as f64 / n as f64); - (a.0 + (b.0 - a.0) * s + (cc.0 - a.0) * t, a.1 + (b.1 - a.1) * s + (cc.1 - a.1) * t) - }; - for i in 0..n { - for j in 0..n - i { - let (p0, p1, p2) = (at(i, j), at(i + 1, j), at(i, j + 1)); - tris.push([(corner(p0), p0), (corner(p1), p1), (corner(p2), p2)]); - if i + j + 1 < n { - let p3 = at(i + 1, j + 1); - tris.push([(corner(p1), p1), (corner(p3), p3), (corner(p2), p2)]); - } - } - } - if tris.len() > 4_000_000 { - return Err(KernelError::Failed("the face is too finely divided to trace".into())); - } - } - let gv: Vec = corner_uv.iter().map(|uv| g.eval(from_p3(surface.subs(uv.0, uv.1)))).collect(); - if gv.iter().any(|x| !x.is_finite()) { - return Err(KernelError::Failed("the splitting tool cannot be evaluated on the face".into())); - } - let side = |i: usize| gv.get(i).copied().unwrap_or(0.0) >= 0.0; - // Mesh edges crossed by the zero set, and which triangles link them. - let key = |a: usize, b: usize| (a.min(b), a.max(b)); - let mut node_uv: HashMap<(usize, usize), ((f64, f64), (f64, f64))> = HashMap::new(); - let mut links: HashMap<(usize, usize), Vec<(usize, usize)>> = HashMap::new(); - let mut edge_tris: HashMap<(usize, usize), usize> = HashMap::new(); - for t in &tris { - for k in 0..3 { - let (a, b) = (t[k], t[(k + 1) % 3]); - *edge_tris.entry(key(a.0, b.0)).or_default() += 1; - } - let crossed: Vec<(usize, usize)> = (0..3) - .filter_map(|k| { - let (a, b) = (t[k], t[(k + 1) % 3]); - (side(a.0) != side(b.0)).then(|| { - let kk = key(a.0, b.0); - let (ua, ub) = if a.0 < b.0 { (a.1, b.1) } else { (b.1, a.1) }; - node_uv.entry(kk).or_insert((ua, ub)); - kk - }) + + /// The patch with every unbounded range closed where it has passed all of `reach`: the + /// span of the box's corners in that parameter, which the surface is straight along (a + /// plane's two, a cylinder's or a cone's height). + fn bounded(&self, reach: &Aabb3) -> Result { + let corners: Vec = (0..8) + .map(|k| { + let pick = |bit: usize, lo: f64, hi: f64| if k & bit == 0 { lo } else { hi }; + Vec3::new(pick(1, reach.min.x, reach.max.x), pick(2, reach.min.y, reach.max.y), pick(4, reach.min.z, reach.max.z)) }) .collect(); - if let [x, y] = crossed[..] { - links.entry(x).or_default().push(y); - links.entry(y).or_default().push(x); - } - } - // Exact point on a crossed mesh edge: bisection along its parameter segment. - let exact = |k: (usize, usize)| -> Option<(mt::Point3, (f64, f64))> { - let (ua, ub) = *node_uv.get(&k)?; - let f = |s: f64| { - let uv = (ua.0 + (ub.0 - ua.0) * s, ua.1 + (ub.1 - ua.1) * s); - (g.eval(from_p3(surface.subs(uv.0, uv.1))), uv) + // The parameters of a corner along the surface's straight directions. + let params = |c: Vec3| -> Option<[f64; 2]> { + match &self.surface { + Surface::Plane { origin, u_axis, v_axis } => { + let (du, dv, r) = (v3(*u_axis), v3(*v_axis), c - v3(*origin)); + let (a, b, d) = (du.dot(du), du.dot(dv), dv.dot(dv)); + let det = a * d - b * b; + (det.abs() > 0.0).then(|| [(d * du.dot(r) - b * dv.dot(r)) / det, (a * dv.dot(r) - b * du.dot(r)) / det]) + } + Surface::Cylinder { frame, .. } | Surface::Cone { frame, .. } => { + let z = v3(frame.x_axis).cross(v3(frame.y_axis)); + Some([0.0, (c - v3(frame.origin)).dot(z)]) + } + _ => None, + } }; - let (mut lo, mut hi) = (0.0, 1.0); - let s0 = f(0.0).0 >= 0.0; - for _ in 0..48 { - let mid = 0.5 * (lo + hi); - if (f(mid).0 >= 0.0) == s0 { - lo = mid; - } else { - hi = mid; + let mut domain = self.domain; + for (k, range) in domain.iter_mut().enumerate() { + if range.iter().all(|x| x.is_finite()) { + continue; + } + let along: Option> = corners.iter().map(|c| params(*c).map(|p| p[k])).collect(); + let along = along.ok_or_else(|| invalid("not supported yet: a split tool unbounded on this surface"))?; + let (lo, hi) = along.iter().fold((f64::INFINITY, f64::NEG_INFINITY), |(lo, hi), x| (lo.min(*x), hi.max(*x))); + let pad = (hi - lo) * 0.05 + 1.0; + *range = [range[0].max(lo - pad), range[1].min(hi + pad)]; + if !(range[0] < range[1]) { + // The tool ends before the body: an empty patch, which cuts nothing. + *range = [range[0], range[0]]; } } - let (val, uv) = f(0.5 * (lo + hi)); - // A side whose ends were on the same side after all (a corner exactly on the tool). - if !(val.abs() <= size * 1e-7) { - return None; - } - Some((surface.subs(uv.0, uv.1), uv)) - }; - // Walk the links into chains. - let mut seen: std::collections::HashSet<(usize, usize)> = std::collections::HashSet::new(); - let mut chains = Vec::new(); - let mut nodes: Vec<(usize, usize)> = links.keys().copied().collect(); - nodes.sort(); - // Open chains start at mesh boundary edges (one triangle). - nodes.sort_by_key(|k| edge_tris.get(k).copied().unwrap_or(0) != 1); - for start in nodes { - if seen.contains(&start) { - continue; - } - let mut path = vec![start]; - seen.insert(start); - let mut cur = start; - let closed; - loop { - let next = links.get(&cur).into_iter().flatten().find(|n| !seen.contains(*n)).copied(); - match next { - Some(n) => { - seen.insert(n); - path.push(n); - cur = n; - } - None => { - closed = path.len() > 2 && links.get(&cur).is_some_and(|l| l.contains(&start)); - break; - } - } - if path.len() > 1_000_000 { - return Err(KernelError::Failed("split line too long".into())); - } - } - // Points about a step apart (the ends kept). - let all: Vec<_> = path.iter().filter_map(|k| exact(*k)).collect(); - let mut pts: Vec<(mt::Point3, (f64, f64))> = Vec::new(); - for (i, p) in all.iter().enumerate() { - let keep = i == 0 || i + 1 == all.len() || pts.last().is_none_or(|l| l.0.distance(p.0) >= step * 0.5); - if keep { - pts.push(*p); - } - } - if pts.len() >= 2 { - chains.push(Chain { pts, closed }); - } + Ok(Patch { domain, ..self.clone() }) } - Ok(chains) } -/// A plane across a flat face (a plane, or a B-spline surface that is one) with straight edges, -/// exactly: the line the two planes meet in, inside the face between pairs of its crossings with -/// the boundary, taken in order along the line. Which side a corner is on is decided once, at -/// the corner, so a plane through a corner crosses the boundary there once or not at all, as it -/// should. `None` (traced on a mesh instead) for any other face or tool, and where an edge lies -/// in the plane. -fn trace_flat(face: &mt::Face, g: &Field, size: f64) -> Option> { - let Field::Plane { o, n } = g else { return None }; - let surface = face.surface(); - let Field::Plane { n: m, .. } = surface_field(&surface) else { return None }; - let along = m.cross(*n).normalized()?; - let tol = size * 1e-7 + 1e-9; - let mut hits: Vec = Vec::new(); - for e in face.boundaries().iter().flat_map(|w| w.edge_iter()) { - let (a, b) = (from_p3(e.front().point()), from_p3(e.back().point())); - let (ga, gb) = ((a - *o).dot(*n), (b - *o).dot(*n)); - // Straight: a line, or a curve (a B-spline of degree one) whose points all lie on its - // chord; such an edge crosses the plane once at most. - let straight = matches!(e.curve(), mt::Curve::Line(_)) || { - let c = e.curve(); - let (t0, t1) = c.range_tuple(); - (1..8).all(|k| from_p3(c.subs(t0 + (t1 - t0) * k as f64 / 8.0)).dist_to_segment(a, b) <= tol) - }; - if !straight || (ga.abs() <= tol && gb.abs() <= tol) || !ga.is_finite() || !gb.is_finite() { - return None; - } - if (ga >= 0.0) != (gb >= 0.0) { - hits.push(p3(a + (b - a) * (ga / (ga - gb)).clamp(0.0, 1.0))); - } - } - if !hits.len().is_multiple_of(2) { - return None; - } - hits.sort_by(|a, b| from_p3(*a).dot(along).total_cmp(&from_p3(*b).dot(along))); - let uv = |p: mt::Point3| surface.search_nearest_parameter(p, None, 50); - let mut chains = Vec::new(); - for pair in hits.chunks(2) { - let [a, b] = pair else { return None }; - // A corner the plane only touches: no line. - if a.distance(*b) <= tol { - continue; - } - chains.push(Chain { pts: vec![(*a, uv(*a)?), (*b, uv(*b)?)], closed: false }); - } - Some(chains) +/// A box past the body (by its own size), which every unbounded patch is cut to. +fn reach_of(body: &Body) -> Aabb3 { + let b = body.outline_bounds(); + let pad = Vec3::new(1.0, 1.0, 1.0) * (b.diagonal() + 1.0); + Aabb3 { min: b.min - pad, max: b.max + pad } } -/// Where the zero set crosses a boundary edge of the face (exact, by bisection on the curve). -fn boundary_hits(face: &mt::Face, g: &Field) -> Vec { - let mut out = Vec::new(); - for w in face.absolute_boundaries() { - for e in w.edge_iter() { - let c = e.curve(); - let (t0, t1) = c.range_tuple(); - let n = 64; - let val = |t: f64| g.eval(from_p3(c.subs(t))); - let mut prev = (t0, val(t0)); - for k in 1..=n { - let t = t0 + (t1 - t0) * k as f64 / n as f64; - let cur = (t, val(t)); - if (prev.1 >= 0.0) != (cur.1 >= 0.0) { - let (mut lo, mut hi) = (prev.0, cur.0); - let s0 = prev.1 >= 0.0; - for _ in 0..60 { - let mid = 0.5 * (lo + hi); - if (val(mid) >= 0.0) == s0 { - lo = mid; - } else { - hi = mid; - } - } - let t = 0.5 * (lo + hi); - out.push(Hit { edge: e.clone(), t, p: c.subs(t) }); - } - prev = cur; +/// The tool as patches, one group for each field (a sketch tool: one for each curve). +fn cutters(tool: &SplitTool, reach: &Aabb3) -> Result>> { + match tool { + SplitTool::Plane(p) => { + let (x, y) = (p.x.normalized(), p.normal().cross(p.x).normalized()); + let (Some(x), Some(y)) = (x, y) else { return Err(invalid("a split plane with no normal")) }; + let surface = Surface::Plane { origin: a3(p.origin), u_axis: a3(x), v_axis: a3(y) }; + Ok(vec![vec![Patch { surface, domain: [ALL, ALL], flip: false }]]) + } + SplitTool::Face { body, face } => Ok(vec![vec![face_patch(body, *face)?]]), + SplitTool::Curves { plane, curves, along } => { + let sweep = sweep_of(plane, *along)?; + if curves.iter().map(Vec::len).sum::() > MAX_SEGMENTS { + return Err(invalid("too many split curves")); } + curves.iter().map(|c| curve_patches(plane, c, sweep, reach)).collect() } } - out } -/// The exact boundary crossings each open chain runs between (`None` for closed chains and -/// for chains whose ends have no crossing of their own). -/// -/// A chain ends where the zero set crosses the boundary of the face's mesh, whose sides are -/// chords up to `mesh_tol` inside a curved boundary: the true crossing lies further along the -/// zero set, by the chord's sag over the sine of the angle the two meet at. A fixed snap -/// distance below the sag drops splits depending on where the mesh's vertices fall (a wall -/// crossing a disc's rim 0.06 mm past the last chord on a 50 mm part). Ends are matched within -/// a few sags (crossings at 15° or steeper), one crossing to one end, nearest pairs first, so -/// a wider reach never lets two ends share a crossing. -fn chain_ends(chains: &[Chain], hits: &[Hit], mesh_tol: f64) -> Vec> { - let reach = 4.0 * mesh_tol; - // (distance, end slot: chain * 2 + 0 for the first point, 1 for the last; hit index) - let mut pairs: Vec<(f64, usize, usize)> = Vec::new(); - for (c, ch) in chains.iter().enumerate() { - if ch.closed { - continue; - } - for (k, end) in [ch.pts.first(), ch.pts.last()].into_iter().enumerate() { - let Some(&(p, _)) = end else { continue }; - for (h, hit) in hits.iter().enumerate() { - let d = hit.p.distance(p); - if d <= reach { - pairs.push((d, c * 2 + k, h)); - } - } - } - } - pairs.sort_by(|x, y| x.0.total_cmp(&y.0)); - let mut slot: Vec> = vec![None; chains.len() * 2]; - let mut taken = vec![false; hits.len()]; - for (_, end, h) in pairs { - if slot.get(end).is_some_and(Option::is_none) && taken.get(h).is_some_and(|t| !*t) { - if let Some(s) = slot.get_mut(end) { - *s = Some(h); - } - if let Some(t) = taken.get_mut(h) { - *t = true; - } - } - } - (0..chains.len()) - .map(|c| match (slot.get(c * 2).copied().flatten(), slot.get(c * 2 + 1).copied().flatten()) { - (Some(a), Some(b)) => Some((hits.get(a)?.clone(), hits.get(b)?.clone())), +/// The surface of face `face` of `body`, over its surface's whole parameter box (unbounded +/// where the surface is), facing out of the body. +fn face_patch(body: &Body, face: usize) -> Result { + body.require_brep("a split tool")?; + let record = crate::native::brep_record(body)?; + let f = record.faces.get(face).ok_or_else(|| invalid(format!("the tool body has no face {face}")))?; + let forward = record + .shells + .iter() + .find_map(|s| match s { + brep::Shell::Faces { faces, .. } => faces.iter().find(|o| o.face as usize == face).map(|o| o.forward), _ => None, }) - .collect() + .ok_or_else(|| invalid(format!("face {face} of the tool body lies in no shell")))?; + // The surface's own parameter box, unbounded where it is. + let dom = SurfaceEval::new(&f.surface).map_err(|e| invalid(format!("the split tool's surface: {e}")))?.domain(); + let range = |k: usize| dom.range[k].map_or(ALL, |(a, b)| [a, b]); + Ok(Patch { surface: f.surface.clone(), domain: [range(0), range(1)], flip: !forward }) } -/// A curve through points: a line when they are collinear, else a cubic through every point -/// (Catmull-Rom tangents; `closed` when the last point repeats the first). -pub(crate) fn fit(pts: &[mt::Point3], scale: f64) -> Result { - let (Some(a), Some(b)) = (pts.first(), pts.last()) else { return Err(KernelError::Failed("empty split line".into())) }; - let d = b - a; - let len = d.magnitude(); - let straight = len > 0.0 && pts.iter().all(|p| (p - a).cross(d).magnitude() / len <= scale * 1e-9); - if straight { - return Ok(mt::Curve::Line(mt::Line(*a, *b))); +/// The unit sweep direction of a sketch tool, on the sketch normal's side (`None`: along the +/// normal). +fn sweep_of(plane: &Plane, along: Option) -> Result { + let n = plane.normal().normalized().ok_or_else(|| invalid("a sketch plane with no normal"))?; + let Some(d) = along else { return Ok(n) }; + let d = d.normalized().filter(|d| d.is_finite()).ok_or_else(|| invalid("the projection direction must be a nonzero vector"))?; + if d.dot(n).abs() < 1e-6 { + return Err(invalid("the projection direction runs along the sketch plane")); } - let m = pts.len().saturating_sub(1); - if m < 1 { - return Err(KernelError::Failed("split line too short".into())); + Ok(if d.dot(n) < 0.0 { -d } else { d }) +} + +/// One sketch curve swept along `sweep`: a patch for each segment, unbounded along the sweep, +/// and an open curve's ends continued along their tangents past `reach`. Each patch's normal +/// points to the curve's left. +fn curve_patches(plane: &Plane, segs: &[Seg2], sweep: Vec3, reach: &Aabb3) -> Result> { + let (Some(first), Some(last)) = (segs.first(), segs.last()) else { return Err(invalid("an empty split curve")) }; + let ok = |p: Vec2| p.x.is_finite() && p.y.is_finite(); + if !segs.iter().all(|s| ok(s.start()) && ok(s.end())) { + return Err(invalid("split curves need finite segments")); } - let closed = a.distance(*b) <= scale * 1e-9 && m >= 3; - let at = |i: isize| -> mt::Point3 { - let i = if closed { i.rem_euclid(m as isize) as usize } else { i.clamp(0, m as isize) as usize }; - pts.get(i).copied().unwrap_or(*a) + let closed = (first.start() - last.end()).len() <= 1e-9 * (1.0 + first.start().len()); + let far = reach.diagonal() + (plane.origin - (reach.min + reach.max) * 0.5).len(); + let mut out = Vec::new(); + let line = |a: Vec2, b: Vec2, out: &mut Vec| -> Result<()> { + let (wa, wb) = (plane.to_world(a), plane.to_world(b)); + let d = wb - wa; + let len = d.len(); + let d = d.normalized().ok_or_else(|| invalid("a zero-length split line"))?; + // u along the sweep, v along the line: the normal sweep × line is the line's left. + let surface = Surface::Plane { origin: a3(wa), u_axis: a3(sweep), v_axis: a3(d) }; + out.push(Patch { surface, domain: [ALL, [0.0, len]], flip: false }); + Ok(()) }; - let mut ctrl = vec![*a]; - let mut knots = vec![0.0; 4]; - for i in 0..m as isize { - let (p0, p1) = (at(i), at(i + 1)); - let t0 = (at(i + 1) - at(i - 1)) / 6.0; - let t1 = (at(i + 2) - at(i)) / 6.0; - ctrl.extend([p0 + t0, p1 - t1, p1]); - let k = (i + 1) as f64; - knots.extend(if i + 1 == m as isize { [k; 4].to_vec() } else { [k; 3].to_vec() }); - } - let bsp = mt::BSplineCurve::try_new(mt::KnotVec::from(knots), ctrl).map_err(|e| KernelError::Failed(format!("split line: {e}")))?; - Ok(mt::Curve::BSplineCurve(bsp)) -} - -/// The curve of a split line: an exact circular arc when a cylinder meets a plane square to its -/// axis (a sketch arc splitting a planar face), else [`fit`]. -fn fit_on(g: &Field, face: &mt::Face, pts: &[mt::Point3], scale: f64) -> Result { - if let (Field::Cylinder { axis, .. }, mt::Surface::Plane(pl)) = (g, face.surface()) - && let (Some(a), Some(b), Some(m)) = (pts.first(), pts.last(), pts.get(pts.len() / 2)) - && a.distance(*b) > scale * 1e-9 - { - let n = pl.normal(); - if n.cross(crate::body::v3(axis.vec())).magnitude() < 1e-9 { - let (va, vb) = (mt::Vertex::new(*a), mt::Vertex::new(*b)); - let e = guard("split arc", || Ok(mt::builder::circle_arc(&va, &vb, *m)))?; - return Ok(e.curve()); + // An open curve runs on along its end tangents: an end line lengthened, else a line added. + let mut segs = segs.to_vec(); + let (mut lead, mut tail) = (None, None); + if !closed { + let t0 = first.start_tangent().normalized().ok_or_else(|| invalid("a split curve with no tangent"))?; + let t1 = last.end_tangent().normalized().ok_or_else(|| invalid("a split curve with no tangent"))?; + let (s0, s1) = (first.start(), last.end()); + match segs.first_mut() { + Some(Seg2::Line { a, .. }) => *a = s0 - t0 * far, + _ => lead = Some((s0 - t0 * far, s0)), + } + match segs.last_mut() { + Some(Seg2::Line { b, .. }) => *b = s1 + t1 * far, + _ => tail = Some((s1, s1 + t1 * far)), } } - fit(pts, scale) -} - -/// Cut an edge at a parameter (none when the point is already a vertex): the vertex and the -/// pieces in the edge's absolute direction. -fn cut_edge(e: &mt::Edge, t: f64, p: mt::Point3, tol: f64) -> Result<(mt::Vertex, Option<(mt::Edge, mt::Edge)>)> { - let abs = if e.orientation() { e.clone() } else { e.inverse() }; - if abs.front().point().distance2(p) <= tol * tol { - return Ok((abs.front().clone(), None)); + if let Some((a, b)) = lead { + line(a, b, &mut out)?; } - if abs.back().point().distance2(p) <= tol * tol { - return Ok((abs.back().clone(), None)); - } - let v = mt::Vertex::new(p); - // `t` is on the edge as it was: when an earlier crossing already cut it (a split line - // crossing it twice), the parameter is found again on the piece. - let cut = abs.cut_with_parameter(&v, t).or_else(|| { - use mt::SearchNearestParameter; - let t = abs.curve().search_nearest_parameter(p, None, 100)?; - abs.cut_with_parameter(&v, t) - }); - let (a, b) = cut.ok_or_else(|| KernelError::Failed("cannot cut an edge at the split line".into()))?; - Ok((v, Some((a, b)))) -} - -/// Replace cut edges (absolute edge → its two pieces) in a face's loops. -fn recut(face: &mt::Face, cuts: &[(mt::Edge, (mt::Edge, mt::Edge))]) -> mt::Face { - let touches = face.absolute_boundaries().iter().flat_map(|w| w.edge_iter()).any(|e| cuts.iter().any(|(c, _)| c.id() == e.id())); - if !touches { - return face.clone(); - } - let wires: Vec = face - .absolute_boundaries() - .iter() - .map(|w| { - let mut out: Vec = Vec::new(); - for e in w.edge_iter() { - expand(e, cuts, &mut out, 0); + let square = sweep.cross(plane.normal()).len() <= 1e-12; + for s in &segs { + match *s { + Seg2::Line { a, b } => line(a, b, &mut out)?, + Seg2::Arc { center, radius, start, sweep: turn } if square && radius.is_finite() && radius > 0.0 && turn.abs() <= TAU => { + // Angles about the sketch normal from the sketch's x; a counter-clockwise arc has + // its left inside, against the cylinder's outward normal. + let n = plane.normal().normalized().ok_or_else(|| invalid("a sketch plane with no normal"))?; + let surface = Surface::Cylinder { frame: Patch::frame(plane.to_world(center), plane.x, n.cross(plane.x)), radius }; + let (lo, hi) = if turn >= 0.0 { (start, start + turn) } else { (start + turn, start) }; + out.push(Patch { surface, domain: [[lo, hi], ALL], flip: turn > 0.0 }); } - out.into() - }) - .collect(); - let mut f = mt::Face::new_unchecked(wires, face.surface()); - if !face.orientation() { - f.invert(); - } - f -} - -/// `e` as the pieces it was cut into, in its direction; a piece cut again (a line crossing an -/// edge twice) as its own pieces. -fn expand(e: &mt::Edge, cuts: &[(mt::Edge, (mt::Edge, mt::Edge))], out: &mut Vec, depth: usize) { - match cuts.iter().find(|(c, _)| c.id() == e.id()).filter(|_| depth < 64) { - Some((c, (a, b))) if c.orientation() == e.orientation() => { - expand(a, cuts, out, depth + 1); - expand(b, cuts, out, depth + 1); - } - Some((_, (a, b))) => { - expand(&b.inverse(), cuts, out, depth + 1); - expand(&a.inverse(), cuts, out, depth + 1); - } - None => out.push(e.clone()), - } -} - -/// The parameter of the point of `c` nearest `q`, searched from the nearest of points along it (a -/// long arc's search from its middle can settle on the far side). -fn nearest(c: &mt::Curve, q: mt::Point3) -> Option { - let (t0, t1) = c.range_tuple(); - let hint = (0..=16).map(|k| t0 + (t1 - t0) * k as f64 / 16.0).min_by(|x, y| c.subs(*x).distance2(q).total_cmp(&c.subs(*y).distance2(q))); - c.search_nearest_parameter(q, hint, 100) -} - -/// `body` with its edges cut at the given points where they pass through them (not at their -/// ends): the faces either side keep the same region, bounded by the pieces. Used to give two -/// faces that cover the same region the same vertices before they are glued (#674). -pub(crate) fn split_edges_at(body: &Body, points: &[Vec3]) -> Result { - let size = body.size(); - let tol = size * 1e-7 + 1e-9; - guard("split edges", || { - let mut cuts: Vec<(mt::Edge, (mt::Edge, mt::Edge))> = Vec::new(); - let mut rec = Rebuild::default(); - // Each current piece and the body edge it is a piece of. - let mut current: Vec<(mt::Edge, mt::Edge)> = body - .solid - .edge_iter() - .map(|e| { - let abs = if e.orientation() { e.clone() } else { e.inverse() }; - (abs.clone(), abs) - }) - .collect(); - for p in points { - let q = p3(*p); - let hit = current.iter().position(|(e, _)| { - if e.front().point().distance(q) <= tol * 10.0 || e.back().point().distance(q) <= tol * 10.0 { - return false; + ref s => { + for (pts, w) in beziers(s)? { + out.push(extruded(plane, &pts, &w, sweep, reach)?); } - let c = e.curve(); - let (t0, t1) = c.range_tuple(); - nearest(&c, q).is_some_and(|t| t > t0.min(t1) && t < t0.max(t1) && c.subs(t).distance(q) <= tol * 10.0) - }); - let Some(k) = hit else { continue }; - let (e, root) = current.swap_remove(k); - let Some(t) = nearest(&e.curve(), q) else { continue }; - let (_, cut) = cut_edge(&e, t, q, tol)?; - if let Some((a, b)) = cut { - rec.piece(&a, &root); - rec.piece(&b, &root); - current.push((a.clone(), root.clone())); - current.push((b.clone(), root.clone())); - cuts.push((e, (a, b))); - } else { - current.push((e, root)); } } - if cuts.is_empty() { - return Ok(body.clone()); + } + if let Some((a, b)) = tail { + line(a, b, &mut out)?; + } + Ok(out) +} + +/// A segment as rational Bézier pieces (control points and weights), exactly: an arc in pieces +/// of a quarter turn at most. +fn beziers(s: &Seg2) -> Result, Vec)>> { + Ok(match *s { + Seg2::Line { a, b } => vec![(vec![a, b], vec![1.0; 2])], + Seg2::Arc { center, radius, start, sweep } => { + if !(radius.is_finite() && radius > 0.0 && sweep.is_finite() && sweep.abs() <= TAU && sweep != 0.0) { + return Err(invalid("a split arc needs a positive radius and a turn of at most one revolution")); + } + let n = (sweep.abs() / std::f64::consts::FRAC_PI_2).ceil().max(1.0) as usize; + let h = sweep / (2 * n) as f64; + let at = |t: f64| center + Vec2::new(t.cos(), t.sin()) * radius; + (0..n) + .map(|i| { + let t0 = start + 2.0 * h * i as f64; + let apex = center + Vec2::new((t0 + h).cos(), (t0 + h).sin()) * (radius / h.cos()); + (vec![at(t0), apex, at(t0 + 2.0 * h)], vec![1.0, h.cos(), 1.0]) + }) + .collect() + } + Seg2::Cubic { p0, p1, p2, p3 } => vec![(vec![p0, p1, p2, p3], vec![1.0; 4])], + Seg2::Conic { a, apex, b, w } => vec![(vec![a, apex, b], vec![1.0, w, 1.0])], + Seg2::Bezier { n, p } => { + let pts = p.get(..=usize::from(n)).ok_or_else(|| invalid("Bézier degree"))?.to_vec(); + let w = vec![1.0; pts.len()]; + vec![(pts, w)] } - let shells: Vec = body - .solid - .boundaries() - .iter() - .map(|sh| { - sh.face_iter() - .map(|f| { - let nf = recut(f, &cuts); - rec.whole(&nf, f); - nf - }) - .collect::>() - .into() - }) - .collect(); - let solid = mt::Solid::new_unchecked(shells); - let split = Body { - solid: std::sync::Arc::new(solid), - size: 0.0, - boundary_tol: body.boundary_tol, - volume: Some(std::sync::Arc::new(std::sync::Mutex::new(None))), - coarse_missing: None, - mesh: None, - color: body.color, - paint: None, - splits: body.splits.clone(), - face_sources: None, - history: None, - display_of: None, - meshes: None, - phased: false, - tolerances: None, - regularization_protected: true, - }; - Ok(crate::record::rebuilt_body(split, &[&body.solid], &rec)) }) } -/// A face rebuilt from loops (absolute) on the same surface and orientation as `like`. -fn face_like(like: &mt::Face, wires: Vec) -> mt::Face { - let mut f = mt::Face::new_unchecked(wires, like.surface()); - if !like.orientation() { - f.invert(); +/// A rational Bézier piece in `plane` swept along `sweep` across `reach`, exactly: a spline of +/// degree one along the sweep (u) and the piece's along the curve (v), so its normal is the +/// curve's left. +fn extruded(plane: &Plane, pts: &[Vec2], w: &[f64], sweep: Vec3, reach: &Aabb3) -> Result { + let n = plane.normal().normalized().ok_or_else(|| invalid("a sketch plane with no normal"))?; + // The sweep heights that take the sketch plane through the whole box. + let (mut lo, mut hi) = (f64::INFINITY, f64::NEG_INFINITY); + for k in 0..8 { + let pick = |bit: usize, a: f64, b: f64| if k & bit == 0 { a } else { b }; + let c = Vec3::new(pick(1, reach.min.x, reach.max.x), pick(2, reach.min.y, reach.max.y), pick(4, reach.min.z, reach.max.z)); + let h = (c - plane.origin).dot(n) / sweep.dot(n); + (lo, hi) = (lo.min(h), hi.max(h)); } - f -} - -/// Polygon of a loop in the surface's parameters (for inside tests). -fn loop_uv(s: &mt::Surface, w: &mt::Wire) -> Vec<(f64, f64)> { - let mut out = Vec::new(); - let mut hint = None; - for e in w.edge_iter() { - let c = e.oriented_curve(); - let (t0, t1) = c.range_tuple(); - for k in 0..8 { - let p = c.subs(t0 + (t1 - t0) * k as f64 / 8.0); - if let Some(uv) = s.search_nearest_parameter(p, hint, 50) { - hint = Some(uv); - out.push(uv); - } - } - } - out -} - -fn inside(poly: &[(f64, f64)], q: (f64, f64)) -> bool { - let mut c = false; - let n = poly.len(); - for i in 0..n { - let (Some(a), Some(b)) = (poly.get(i), poly.get((i + 1) % n)) else { continue }; - if (a.1 > q.1) != (b.1 > q.1) && q.0 < a.0 + (b.0 - a.0) * (q.1 - a.1) / (b.1 - a.1) { - c = !c; - } - } - c -} - -/// Split one face (index `fi` in `faces`) by one field; `faces` is updated (the face replaced -/// by its pieces, boundary edges cut in every face). -/// Returns how many lines split it and how many faces now stand in its place (at `fi`). The -/// pieces of cut edges are recorded in `rec`. -fn split_one(faces: &mut Vec, fi: usize, g: &Field, size: f64, new_edges: &mut Vec, rec: &mut Rebuild) -> Result<(usize, usize)> { - let Some(face) = faces.get(fi).cloned() else { return Ok((0, 1)) }; - let chains = trace(&face, g, size)?; - if chains.is_empty() { - return Ok((0, 1)); - } - let tol = size * 1e-7 + 1e-9; - let hits = boundary_hits(&face, g); - let ends = chain_ends(&chains, &hits, trace_tol(&face, size)); - // Pieces of the original face as they are split. - let mut pieces: Vec = vec![face.clone()]; - let mut all_cuts: Vec<(mt::Edge, (mt::Edge, mt::Edge))> = Vec::new(); - let mut made = 0; - for (ch, end) in chains.into_iter().zip(ends) { - if ch.closed { - // A closed line inside the face: a new face inside it and a hole in the piece that holds it. - let mut pts: Vec = ch.pts.iter().map(|x| x.0).collect(); - if let Some(first) = pts.first().copied() { - pts.push(first); - } - let curve = fit(&pts, size)?; - let Some(&(p0, uv0)) = ch.pts.first() else { continue }; - let v = mt::Vertex::new(p0); - let e = mt::Edge::new_unchecked(&v, &v, curve); - new_edges.push(e.clone()); - let uv: Vec<(f64, f64)> = ch.pts.iter().map(|x| x.1).collect(); - let area: f64 = (0..uv.len()) - .map(|i| { - let (a, b) = (uv[i], uv[(i + 1) % uv.len()]); - a.0 * b.1 - b.0 * a.1 + let degree = pts.len().saturating_sub(1); + let knots = |d: usize| std::iter::repeat_n(0.0, d + 1).chain(std::iter::repeat_n(1.0, d + 1)).collect::>(); + let homogeneous: Vec> = [lo, hi] + .iter() + .map(|h| { + pts.iter() + .zip(w) + .map(|(p, wt)| { + let q = plane.to_world(*p) + sweep * *h; + [q.x * wt, q.y * wt, q.z * wt, *wt] }) - .sum(); - let s = face.surface(); - let Some(pi) = pieces.iter().position(|f| { - let ws = f.absolute_boundaries(); - ws.first().is_some_and(|o| inside(&loop_uv(&s, o), uv0)) && !ws.iter().skip(1).any(|h| inside(&loop_uv(&s, h), uv0)) - }) else { - continue; - }; - let Some(host) = pieces.get(pi).cloned() else { continue }; - // The new loop counter-clockwise in parameters (an outer loop); its reverse a hole. - let ccw: mt::Wire = if area > 0.0 { vec![e.clone()].into() } else { vec![e.inverse()].into() }; - let hole: mt::Wire = ccw.inverse(); - let ws = host.absolute_boundaries().clone(); - let mut outer_ws = vec![]; - let mut inner_ws = vec![ccw.clone()]; - for (k, w) in ws.iter().enumerate() { - let q = loop_uv(&s, w).first().copied(); - if k > 0 && q.is_some_and(|q| inside(&uv, q)) { - inner_ws.push(w.clone()); - } else { - outer_ws.push(w.clone()); - } - } - outer_ws.push(hole); - let a = face_like(&host, outer_ws); - let b = face_like(&host, inner_ws); - pieces.splice(pi..=pi, [a, b]); - made += 1; - continue; - } - // An open line: it must end on the boundary at two crossings. - let Some((ha, hb)) = end else { continue }; - // Cut the boundary edges (an edge may already have been cut by an earlier line). - let vert = |h: &Hit, cuts: &mut Vec<(mt::Edge, (mt::Edge, mt::Edge))>| -> Result { - // The current edge holding the hit: the original, or a piece of it. - let mut e = h.edge.clone(); - for (c, (x, y)) in cuts.iter() { - if c.id() == e.id() || c.id() == e.inverse().id() { - // The piece the hit lies on (a cut line's pieces are parametrised afresh, - // so by where the point is, not by the original parameter). - let off = |pc: &mt::Edge| { - use mt::SearchNearestParameter; - let cc = pc.curve(); - let (t0, t1) = cc.range_tuple(); - let t = cc.search_nearest_parameter(h.p, None, 100).unwrap_or(h.t).clamp(t0.min(t1), t0.max(t1)); - cc.subs(t).distance(h.p) - }; - e = if off(x) <= off(y) { x.clone() } else { y.clone() }; - } - } - let (v, cut) = cut_edge(&e, h.t, h.p, tol)?; - if let Some(pair) = cut { - let abs = if e.orientation() { e.clone() } else { e.inverse() }; - cuts.push((abs, pair)); - } - Ok(v) - }; - let va = vert(&ha, &mut all_cuts)?; - let vb = vert(&hb, &mut all_cuts)?; - if va.id() == vb.id() { - continue; - } - pieces = pieces.iter().map(|f| recut(f, &all_cuts)).collect(); - let mut pts: Vec = ch.pts.iter().map(|x| x.0).collect(); - if let Some(f) = pts.first_mut() { - *f = va.point(); - } - if let Some(l) = pts.last_mut() { - *l = vb.point(); - } - // A line from one loop of a piece to another (into a hole): the two loops become one, - // joined by the line run both ways; a later line divides it. - let in_loop = |w: &mt::Wire, v: &mt::Vertex| w.edge_iter().any(|e| e.front().id() == v.id()); - if let Some(pi) = pieces.iter().position(|f| { - let ws = f.absolute_boundaries(); - ws.iter().any(|w| in_loop(w, &va)) && ws.iter().any(|w| in_loop(w, &vb)) && !ws.iter().any(|w| in_loop(w, &va) && in_loop(w, &vb)) - }) && let Some(host) = pieces.get(pi).cloned() - { - let ws = host.absolute_boundaries().clone(); - let (Some(wa_i), Some(wb_i)) = (ws.iter().position(|w| in_loop(w, &va)), ws.iter().position(|w| in_loop(w, &vb))) else { continue }; - let from = |w: &mt::Wire, v: &mt::Vertex| -> Vec { - let edges: Vec = w.edge_iter().cloned().collect(); - let k = edges.iter().position(|e| e.front().id() == v.id()).unwrap_or(0); - edges.iter().skip(k).chain(edges.iter().take(k)).cloned().collect() - }; - let (Some(w1), Some(w2)) = (ws.get(wa_i), ws.get(wb_i)) else { continue }; - let slit = mt::Edge::new_unchecked(&va, &vb, fit_on(g, &face, &pts, size)?); - new_edges.push(slit.clone()); - let mut joined = from(w1, &va); - joined.push(slit.clone()); - joined.extend(from(w2, &vb)); - joined.push(slit.inverse()); - // The outer loop first. - let others = ws.iter().enumerate().filter(|(k, _)| *k != wa_i && *k != wb_i).map(|(_, w)| w.clone()); - let joined: mt::Wire = joined.into(); - let loops: Vec = if wa_i == 0 || wb_i == 0 { - std::iter::once(joined).chain(others).collect() - } else { - others.chain(std::iter::once(joined)).collect() - }; - pieces.splice(pi..=pi, [face_like(&host, loops)]); - continue; - } - // The piece whose loop passes both vertices. - let Some(pi) = pieces.iter().position(|f| { - f.absolute_boundaries().iter().any(|w| { - let ids: Vec<_> = w.edge_iter().map(|e| e.front().id()).collect(); - ids.contains(&va.id()) && ids.contains(&vb.id()) - }) - }) else { - continue; - }; - let Some(host) = pieces.get(pi).cloned() else { continue }; - let ws = host.absolute_boundaries().clone(); - let Some(wi) = ws.iter().position(|w| { - let ids: Vec<_> = w.edge_iter().map(|e| e.front().id()).collect(); - ids.contains(&va.id()) && ids.contains(&vb.id()) - }) else { - continue; - }; - let Some(w) = ws.get(wi) else { continue }; - let edges: Vec = w.edge_iter().cloned().collect(); - let n = edges.len(); - let (Some(ia), Some(ib)) = (edges.iter().position(|e| e.front().id() == va.id()), edges.iter().position(|e| e.front().id() == vb.id())) - else { - continue; - }; - let split = mt::Edge::new_unchecked(&va, &vb, fit_on(g, &face, &pts, size)?); - new_edges.push(split.clone()); - let run = |from: usize, to: usize| -> Vec { - let mut out = Vec::new(); - let mut i = from; - while i != to { - if let Some(e) = edges.get(i) { - out.push(e.clone()); - } - i = (i + 1) % n; - } - out - }; - // Piece A: a → b along the loop, back along the split. Piece B: b → a, then the split. - let mut wa = run(ia, ib); - wa.push(split.inverse()); - let mut wb = run(ib, ia); - wb.push(split.clone()); - let (wa, wb): (mt::Wire, mt::Wire) = (wa.into(), wb.into()); - if !wa.is_closed() || !wb.is_closed() { - return Err(KernelError::Failed("split loops do not close".into())); - } - // Other loops go to the piece on their side. - let s = face.surface(); - let (pa_poly, pb_poly) = (loop_uv(&s, &wa), loop_uv(&s, &wb)); - let (mut la, mut lb) = (vec![wa], vec![wb]); - for (k, other) in ws.iter().enumerate() { - if k == wi { - continue; - } - let q = loop_uv(&s, other).first().copied(); - match q { - Some(q) if inside(&pa_poly, q) => la.push(other.clone()), - Some(q) if inside(&pb_poly, q) => lb.push(other.clone()), - _ => return Err(KernelError::Failed("a hole meets the split line".into())), - } - } - let (a, b) = (face_like(&host, la), face_like(&host, lb)); - pieces.splice(pi..=pi, [a, b]); - made += 1; - } - if made == 0 { - return Ok((0, 1)); - } - // In the order they were cut: a piece cut again is a piece of the edge it came from. - for (c, (a, b)) in &all_cuts { - rec.piece(a, c); - rec.piece(b, c); - } - let count = pieces.len(); - // Neighbours take the cut edges. - let rest: Vec = faces.iter().enumerate().filter(|(k, _)| *k != fi).map(|(_, f)| recut(f, &all_cuts)).collect(); - let mut out = Vec::with_capacity(rest.len() + pieces.len()); - out.extend(rest.iter().take(fi).cloned()); - out.extend(pieces); - out.extend(rest.iter().skip(fi).cloned()); - *faces = out; - Ok((made, count)) + .collect() + }) + .collect(); + let surface = Surface::RationalSpline { u_knots: knots(1), v_knots: knots(degree), homogeneous }; + Ok(Patch { surface, domain: [[0.0, 1.0], [0.0, 1.0]], flip: false }) } -/// The plane or circular cylinder a (B-spline) surface is, as a field with the surface's -/// normal; other surfaces as they are. -fn surface_field(s: &mt::Surface) -> Field { - let (Some((u0, u1)), Some((v0, v1))) = s.try_range_tuple() else { return Field::Surface(s.clone()) }; - if !matches!(s, mt::Surface::Plane(_) | mt::Surface::BSplineSurface(_) | mt::Surface::NurbsSurface(_)) { - return Field::Surface(s.clone()); - } - let at = |a: f64, b: f64| from_p3(s.subs(u0 + (u1 - u0) * a, v0 + (v1 - v0) * b)); - let nrm = |a: f64, b: f64| { - let n = s.normal(u0 + (u1 - u0) * a, v0 + (v1 - v0) * b); - Vec3::new(n.x, n.y, n.z) - }; - let grid: Vec<(f64, f64)> = (0..=6).flat_map(|i| (0..=6).map(move |j| (i as f64 / 6.0, j as f64 / 6.0))).collect(); - let size = grid.iter().map(|(a, b)| at(*a, *b).dist(at(0.0, 0.0))).fold(1e-9, f64::max); - let tol = size * 1e-7; - // A plane: one normal, every point on it. - let n0 = nrm(0.5, 0.5); - let o = at(0.5, 0.5); - if grid.iter().all(|(a, b)| (at(*a, *b) - o).dot(n0).abs() <= tol) { - return Field::Plane { o, n: n0 }; - } - // A cylinder: straight along one parameter, a circle along the other. - for along_v in [true, false] { - let pt = |t: f64, w: f64| if along_v { at(t, w) } else { at(w, t) }; - let d = pt(0.0, 1.0) - pt(0.0, 0.0); - let Some(axis) = d.unit() else { continue }; - let straight = (0..=6).all(|k| { - let t = k as f64 / 6.0; - (0..=4).all(|m| (pt(t, m as f64 / 4.0) - pt(t, 0.0) - d * (m as f64 / 4.0)).len() <= tol) - }); - if !straight { - continue; - } - let prof: Vec = (0..=8).map(|k| pt(k as f64 / 8.0, 0.0)).collect(); - let (a, b, c) = (prof[0], prof[3], prof[6]); - let (ab, ac) = (b - a, c - a); - let nn = ab.cross(ac); - let n2 = nn.dot(nn); - if !(n2 > 1e-24) { - continue; - } - let centre = a + (nn.cross(ab) * ac.dot(ac) + ac.cross(nn) * ab.dot(ab)) * (0.5 / n2); - let r = (a - centre).len(); - if prof.iter().all(|p| ((*p - centre).len() - r).abs() <= tol) && nn.normalized().is_some_and(|nu| nu.cross(axis).len() < 1e-6) { - let mid = pt(0.5, 0.5); - let v = mid - centre; - let radial = v - axis * v.dot(axis); - let out = if nrm(0.5, 0.5).dot(radial) >= 0.0 { 1.0 } else { -1.0 }; - return Field::Cylinder { o: centre, axis, r, out }; - } - } - Field::Surface(s.clone()) +/// One patch ready to measure points against: its surface prepared once. +struct Measured { + eval: SurfaceEval, + starts: Starts, + domain: [[f64; 2]; 2], + sign: f64, } -/// A plane or circular cylinder behind a surface (for offsetting). +impl Measured { + fn new(p: &Patch) -> Result { + let eval = SurfaceEval::new(&p.surface).map_err(|e| invalid(format!("the split tool's surface: {e}")))?; + Ok(Measured { eval, starts: Starts::default(), domain: p.domain, sign: if p.flip { -1.0 } else { 1.0 } }) + } + + /// The nearest point of the patch to `p` (its parameters kept in the patch's box, a + /// periodic one taken round to it), and `p`'s offset along the patch's side there. + fn near(&self, p: Vec3) -> Option<(f64, f64)> { + let uv = self.eval.project_with(&self.starts, p, None, f64::INFINITY)?; + let period = self.eval.domain().period; + let mut q = [0.0; 2]; + for k in 0..2 { + let [lo, hi] = self.domain[k]; + let mut x = uv[k]; + if let Some(t) = period[k] { + x = lo + (x - lo).rem_euclid(t); + // Past the end of a strip: whichever end is nearer round the turn. + if x > hi && x - hi > lo + t - x { + x -= t; + } + } + q[k] = x.clamp(lo.min(hi), hi.max(lo)); + } + let (at, su, sv) = self.eval.derivs(q[0], q[1])?; + let n = su.cross(sv).normalized()?; + Some(((p - at).len(), (p - at).dot(n) * self.sign)) + } +} + +/// Which side of a tool a point is on, the tool's surfaces prepared once for many points. +pub(crate) struct ToolSide(Vec); + +impl ToolSide { + /// The side field of `tool` (one field: a sketch tool with one curve), its patches reaching + /// past `near`. + pub(crate) fn of(tool: &SplitTool, near: &Body) -> Result { + let mut groups = cutters(tool, &reach_of(near))?; + match (groups.pop(), groups.is_empty()) { + (Some(g), true) => Ok(ToolSide(g.iter().map(Measured::new).collect::>()?)), + _ => Err(invalid("one tool at a time")), + } + } + + /// The side of `p`: its sign says which, zero on the tool (NaN where it can't be told): the + /// offset from the nearest patch along its side. + pub(crate) fn at(&self, p: Vec3) -> f64 { + self.0.iter().filter_map(|m| m.near(p)).min_by(|a, b| a.0.total_cmp(&b.0)).map_or(f64::NAN, |(_, side)| side) + } +} + +/// Which side of a tool a point is on. +pub(crate) fn side_of_point(tool: &SplitTool, near: &Body, p: Vec3) -> Result { + Ok(ToolSide::of(tool, near)?.at(p)) +} + +/// How many fields a tool is (a sketch tool: one per curve). +pub(crate) fn field_count(tool: &SplitTool) -> Result { + Ok(match tool { + SplitTool::Curves { curves, .. } => curves.len(), + _ => 1, + }) +} + +/// A plane or circular cylinder behind a face's surface, facing as the face does. pub(crate) enum Analytic { Plane { n: Vec3 }, Cylinder { o: Vec3, axis: UnitVec3, r: f64, out: f64 }, } -pub(crate) fn analytic_field(s: &mt::Surface) -> Option { - match surface_field(s) { - Field::Plane { n, .. } => Some(Analytic::Plane { n }), - Field::Cylinder { o, axis, r, out } => Some(Analytic::Cylinder { o, axis, r, out }), +/// The plane or circular cylinder face `face` of `body` lies on (a B-spline whose control net is +/// flat is a plane), with the face's orientation; none for other surfaces. +pub(crate) fn analytic_face(body: &Body, face: usize) -> Result> { + let p = face_patch(body, face)?; + let sign = if p.flip { -1.0 } else { 1.0 }; + let flat = |pts: &[Vec3]| -> Option { + let o = *pts.first()?; + let size = pts.iter().map(|q| q.dist(o)).fold(0.0, f64::max); + let n = pts.iter().flat_map(|a| pts.iter().map(move |b| (*a - o).cross(*b - o))).max_by(|a, b| a.len().total_cmp(&b.len()))?.normalized()?; + pts.iter().all(|q| (*q - o).dot(n).abs() <= size * 1e-9).then_some(n) + }; + Ok(match &p.surface { + Surface::Plane { u_axis, v_axis, .. } => v3(*u_axis).cross(v3(*v_axis)).normalized().map(|n| Analytic::Plane { n: n * sign }), + Surface::PlaneThrough { origin, u_point, v_point } => { + (v3(*u_point) - v3(*origin)).cross(v3(*v_point) - v3(*origin)).normalized().map(|n| Analytic::Plane { n: n * sign }) + } + Surface::Cylinder { frame, radius } => { + v3(frame.x_axis).cross(v3(frame.y_axis)).unit().map(|axis| Analytic::Cylinder { o: v3(frame.origin), axis, r: *radius, out: sign }) + } + Surface::Spline { points, .. } => { + let pts: Vec = points.iter().flatten().map(|q| v3(*q)).collect(); + // The net's normal, turned to the surface's own. + let m = Measured::new(&p)?; + flat(&pts).and_then(|n| { + let (_, su, sv) = m.eval.derivs(p.domain[0][0], p.domain[1][0])?; + let own = su.cross(sv); + Some(Analytic::Plane { n: if own.dot(n) < 0.0 { -n } else { n } * sign }) + }) + } _ => None, - } -} - -/// Which side of a tool a point is on. -pub(crate) fn side_of_point(tool: &SplitTool, p: Vec3) -> Result { - Ok(ToolSide::of(tool)?.at(p)) -} - -/// Which side of a tool points are on, the tool's field built once for many points. -pub(crate) struct ToolSide(Field); - -impl ToolSide { - pub(crate) fn of(tool: &SplitTool) -> Result { - let mut f = fields(tool)?; - match (f.pop(), f.is_empty()) { - (Some(g), true) => Ok(ToolSide(g)), - _ => Err(KernelError::Invalid("one tool at a time".into())), - } - } - - /// The side of `p`: its sign says which, zero on the tool (NaN where it can't be told). - pub(crate) fn at(&self, p: Vec3) -> f64 { - self.0.eval(p) - } -} - -/// Which side of a tool a face lies on (its field well inside it). -pub(crate) fn face_side(face: &mt::Face, tool: &SplitTool, size: f64) -> Result { - let f = fields(tool)?; - let [g] = f.as_slice() else { return Err(KernelError::Invalid("one tool at a time".into())) }; - side_of(face, g, size) -} - -fn fields(tool: &SplitTool) -> Result> { - Ok(match tool { - SplitTool::Plane(p) => vec![Field::Plane { o: p.origin, n: p.normal() }], - SplitTool::Face { body: other, face } => { - other.require_brep("a split tool")?; - let f = other.solid.face_iter().nth(*face).ok_or_else(|| KernelError::Invalid(format!("the tool body has no face {face}")))?; - vec![surface_field(&f.oriented_surface())] - } - SplitTool::Curves { plane, curves, along } => { - let n = plane.normal(); - // The sweep direction, turned to the plane's side (a direction along the plane sweeps - // nothing). - let along = match along { - None => None, - Some(d) => { - let d = d - .normalized() - .filter(|d| d.is_finite()) - .ok_or_else(|| KernelError::Invalid("the projection direction must be a nonzero vector".into()))?; - if d.dot(n).abs() < 1e-6 { - return Err(KernelError::Invalid("the projection direction runs along the sketch plane".into())); - } - let d = if d.dot(n) < 0.0 { -d } else { d }; - (d.cross(n).len() > 1e-12).then_some(d) - } - }; - let sweep = along.unwrap_or(n); - let mut out = Vec::new(); - for c in curves { - out.push(match c.as_slice() { - [Seg2::Line { a, b }] => { - let (wa, wb) = (plane.to_world(*a), plane.to_world(*b)); - let n = sweep.cross(wb - wa).normalized().ok_or_else(|| KernelError::Invalid("a zero-length split line".into()))?; - Field::Plane { o: wa, n } - } - [Seg2::Arc { center, radius, .. }] if along.is_none() && radius.is_finite() && *radius > 0.0 => Field::Cylinder { - o: plane.to_world(*center), - axis: n.unit().ok_or_else(|| KernelError::Invalid("a sketch plane with no normal".into()))?, - r: *radius, - out: 1.0, - }, - segs => { - let mut pts: Vec = Vec::new(); - for sg in segs { - for p in sg.polyline(1e-3) { - if pts.last().is_none_or(|q| (*q - p).len() > 1e-9) { - pts.push(p); - } - } - } - if pts.len() < 2 || pts.iter().any(|p| !(p.x.is_finite() && p.y.is_finite())) || pts.len() > 1_000_000 { - return Err(KernelError::Invalid("split curves need finite segments".into())); - } - let span = pts.iter().fold(0.0f64, |m, p| m.max(p.len())); - let closed = pts.len() > 3 && pts.first().zip(pts.last()).is_some_and(|(a, b)| (*a - *b).len() <= span * 1e-9 + 1e-9); - Field::Curve { plane: *plane, pts, closed, along } - } - }); - } - out - } }) } -/// The body's shells with the chosen faces split by each field in turn, the new edges, and -/// how many lines split something. -/// What an imprint made: the shells' faces, for each face the body face it came from, the new -/// edges, and how many lines split each chosen face. -pub(crate) struct Imprint { - shells: Vec>, - origin: Vec, - new_edges: Vec, - made: HashMap, - /// The faces and edges made from the body's (each face over the same region as the body - /// face it came from, or a piece of it; cut edges' pieces). - rec: Rebuild, +/// Which side of a tool each face of `body` lies on (in face order): the field at the centre of +/// the face's triangle where it is largest, so well inside the face. +pub(crate) fn face_sides(body: &Body, tool: &SplitTool) -> Result> { + let side = ToolSide::of(tool, body)?; + let mesh = body.mesh(&crate::MeshTolerance::chord((body.size() * 2e-3).max(1e-4)))?; + let mut best = vec![0.0f64; body.face_count()]; + for (t, f) in mesh.triangles.iter().zip(&mesh.tri_face) { + let Some([a, b, c]) = mesh.tri(t) else { continue }; + let x = side.at((a + b + c) * (1.0 / 3.0)); + if let Some(slot) = best.get_mut(*f as usize) + && x.is_finite() + && x.abs() > slot.abs() + { + *slot = x; + } + } + Ok(best) } -fn imprint(body: &Body, faces: &[usize], fields: &[Field], size: f64) -> Result { - let mut out = Imprint { shells: Vec::new(), origin: Vec::new(), new_edges: Vec::new(), made: HashMap::new(), rec: Rebuild::default() }; - let body_faces: Vec<&mt::Face> = body.solid.face_iter().collect(); - let mut base = 0; - for sh in body.solid.boundaries() { - let mut fs: Vec = sh.face_iter().cloned().collect(); - // Which faces (by position) are to be split: the chosen ones and, later, their pieces; - // and the body face each came from. - let mut target: Vec = (0..fs.len()).map(|k| faces.contains(&(base + k))).collect(); - let mut origin: Vec = (0..fs.len()).map(|k| base + k).collect(); - base += fs.len(); - for g in fields { - let mut k = 0; - while k < fs.len() { - if !target.get(k).copied().unwrap_or(false) { - k += 1; - continue; - } - let src = origin.get(k).copied().unwrap_or(0); - let (m, pieces) = split_one(&mut fs, k, g, size, &mut out.new_edges, &mut out.rec)?; - *out.made.entry(src).or_default() += m; - // The pieces stand at k..k + pieces; positions after shift. - target = target.iter().take(k).copied().chain(std::iter::repeat_n(true, pieces)).chain(target.iter().skip(k + 1).copied()).collect(); - origin = origin.iter().take(k).copied().chain(std::iter::repeat_n(src, pieces)).chain(origin.iter().skip(k + 1).copied()).collect(); - k += pieces; +/// The regions of a body cut by a tool: for each face (in face order), the region it lies in, +/// and how many there are. Faces are in one region when they meet at an edge the tool does not +/// run along; an edge along the tool (a split line) parts them. +pub(crate) fn regions(split: &Body, tool: &SplitTool) -> Result<(Vec, usize)> { + let record = crate::native::brep_record(split)?; + let fields: Vec = cutters(tool, &reach_of(split))? + .iter() + .map(|g| g.iter().map(Measured::new).collect::>>().map(ToolSide)) + .collect::>()?; + let tol = (split.size() * 2e-3).max(1e-4); + let along = |e: &brep::Edge| -> bool { + let Ok(c) = CurveEval::new(&e.curve) else { return false }; + let [t0, t1] = e.interval; + let at: Vec = [0.25, 0.5, 0.75].iter().filter_map(|s| c.point(t0 + (t1 - t0) * s)).collect(); + at.len() == 3 && fields.iter().any(|g| at.iter().all(|p| g.at(*p).abs() <= tol)) + }; + let n = record.faces.len(); + let mut parent: Vec = (0..n).collect(); + fn root(parent: &mut [usize], mut i: usize) -> usize { + while let Some(&p) = parent.get(i).filter(|p| **p != i) { + i = p; + } + i + } + for e in &record.edges { + if along(e) { + continue; + } + let faces: Vec = e.uses.iter().map(|u| u.face as usize).collect(); + for pair in faces.windows(2) { + let (a, b) = (root(&mut parent, pair[0]), root(&mut parent, pair[1])); + if let Some(slot) = parent.get_mut(a.max(b)) { + *slot = a.min(b); } } - for (f, src) in fs.iter().zip(&origin) { - if let Some(from) = body_faces.get(*src) { - out.rec.whole(f, from); - } + } + let mut ids = std::collections::BTreeMap::new(); + let region: Vec = (0..n) + .map(|i| { + let r = root(&mut parent, i); + let next = ids.len(); + *ids.entry(r).or_insert(next) + }) + .collect(); + Ok((region, ids.len())) +} + +/// SMLib's surfaces of the tool's patches, in the tool's order, each reaching past `reach`. +fn smlib_tools(groups: &[Vec], reach: &Aabb3) -> Result> { + let mut out = Vec::new(); + for p in groups.iter().flatten() { + let b = p.bounded(reach)?; + if b.domain.iter().any(|r| !(r[0] < r[1])) { + continue; } - out.shells.push(fs); - out.origin.extend(origin); + out.push(solidcraft_smlib::ssi::Surface::from_record(&b.surface, b.domain).map_err(|e| failed("split tool", e))?); } Ok(out) } +/// The record's position of SMLib's face `i` (identity when `order` is empty). +fn canonical(order: &[usize], i: usize) -> Option { + if order.is_empty() { Some(i) } else { order.get(i).copied() } +} + +/// A result checked by SMLib and brought back with its history over `[body]`. +fn brought_back(name: &str, body: &Body, ob: &solidcraft_smlib::history::Order, r: &solidcraft_smlib::Brep) -> Result<(Body, History)> { + let check = r.check(solidcraft_smlib::CheckLevel::All).map_err(|e| failed(name, e))?; + if !check.is_valid() { + return Err(failed(name, format!("SMLib's check refuses the result: {check}"))); + } + let (out, or) = crate::native::body_from_smlib_ordered(r)?; + let or = or.ok_or_else(|| failed(name, "the result's record is numbered unlike SMLib's body"))?; + let h = r.history().map_err(|e| failed(name, e))?.history(std::slice::from_ref(ob), &or).map_err(|e| failed(name, e))?; + let out = out.with_color(body.color()).with_split_lines(body.split_lines().cloned()); + Ok((out, h)) +} + +/// `body` in SMLib, tagged as operand 0, with its record's order. +fn tagged(name: &str, body: &Body) -> Result<(solidcraft_smlib::Brep, solidcraft_smlib::history::Order)> { + let (sb, ob) = crate::native::smlib_body_ordered(body)?; + let ob = ob.ok_or_else(|| failed(name, "SMLib numbers the body's faces unlike its record"))?; + Ok((sb.tagged(0).map_err(|e| failed(name, e))?, ob)) +} + /// Split the given faces (indices in `body`'s face order) where the tool meets them. The solid /// is unchanged. A listed face the tool misses is an error naming it. The result's history is /// over `[body]`: each piece of a split face is Split from it, every other face Kept. @@ -1139,387 +530,122 @@ pub fn split_faces(body: &Body, faces: &[usize], tool: &SplitTool) -> Result Result<(Body, Vec)> { - let (b, imp) = split_any(body, faces, tool)?; - // The new edges, remembered so that export keeps the pieces apart. - let mut lines: Vec> = body.splits.as_ref().map(|l| l.to_vec()).unwrap_or_default(); - lines.extend(imp.new_edges.iter().filter_map(|e| { - let c = e.curve(); - let (t0, t1) = c.range_tuple(); - let n = if matches!(c, mt::Curve::Line(_)) { 1 } else { 64 }; - let pts: Vec = (0..=n).map(|k| from_p3(c.subs(t0 + (t1 - t0) * k as f64 / n as f64))).collect(); - pts.iter().all(|p| p.is_finite()).then_some(pts) - })); - let b = b.with_split_lines(Some(std::sync::Arc::new(lines))); - if let Some(missed) = faces.iter().find(|f| imp.made.get(f).copied().unwrap_or(0) == 0) { - return Err(KernelError::Invalid(format!("the tool does not cross face {missed}"))); + let (b, map, made) = split(body, faces, tool)?; + if let Some((missed, _)) = faces.iter().zip(&made).find(|(_, m)| **m == 0) { + return Err(invalid(format!("the tool does not cross face {missed}"))); } - Ok((b, imp.origin)) + Ok((b, map)) } -/// Split where the tool meets the listed faces (missing some is fine; missing all is an error). -pub(crate) fn split_any(body: &Body, faces: &[usize], tool: &SplitTool) -> Result<(Body, Imprint)> { - body.require_brep("split face")?; - let size = body.size(); - let fields = fields(tool)?; - let total = body.solid.face_iter().count(); - if faces.iter().any(|f| *f >= total) { - return Err(KernelError::Invalid("no such face".into())); +/// Split where the tool meets the listed faces (missing some is fine; missing all is an error), +/// with the body face each face of the result came from. +pub(crate) fn split_any(body: &Body, faces: &[usize], tool: &SplitTool) -> Result<(Body, Vec)> { + split(body, faces, tool).map(|(b, map, _)| (b, map)) +} + +/// The split, the body face each result face came from, and for each listed face how many faces +/// it became besides itself. +fn split(body: &Body, faces: &[usize], tool: &SplitTool) -> Result<(Body, Vec, Vec)> { + const NAME: &str = "split face"; + body.require_brep(NAME)?; + let total = body.face_count(); + if faces.is_empty() || faces.iter().any(|f| *f >= total) { + return Err(invalid("no such face")); } - guard("split face", || { - let mut imp = imprint(body, faces, &fields, size)?; - if imp.made.values().all(|m| *m == 0) { - return Err(KernelError::Invalid("the tool does not cross the selected faces".into())); - } - let shells = std::mem::take(&mut imp.shells); - let solid = mt::Solid::new_unchecked(shells.into_iter().map(mt::Shell::from).collect()); - let split = Body { - solid: std::sync::Arc::new(solid), - size: 0.0, - boundary_tol: body.boundary_tol, - volume: Some(std::sync::Arc::new(std::sync::Mutex::new(None))), - coarse_missing: None, - mesh: None, - color: body.color, - paint: None, - splits: body.splits.clone(), - face_sources: None, - history: None, - display_of: None, - meshes: None, - phased: false, - tolerances: None, - regularization_protected: true, - }; - Ok((crate::record::rebuilt_body(split, &[&body.solid], &imp.rec), imp)) + if faces.iter().enumerate().any(|(i, f)| faces.get(..i).is_some_and(|before| before.contains(f))) { + return Err(invalid("a face is listed twice")); + } + let reach = reach_of(body); + let tools = smlib_tools(&cutters(tool, &reach)?, &reach)?; + let missed = || invalid("the tool does not cross the selected faces"); + if tools.is_empty() { + return Err(missed()); + } + guard(NAME, || { + let (sb, ob) = tagged(NAME, body)?; + let at: Vec = faces + .iter() + .map(|f| if ob.faces.is_empty() { Some(*f) } else { ob.faces.iter().position(|c| c == f) }) + .collect::>() + .ok_or_else(|| failed(NAME, "a listed face is not in SMLib's body"))?; + let refs: Vec<&solidcraft_smlib::ssi::Surface> = tools.iter().collect(); + let (r, made) = sb.split_faces(&at, &refs).map_err(|e| failed(NAME, e))?.ok_or_else(missed)?; + let (out, h) = brought_back(NAME, body, &ob, &r)?; + let map = h + .faces + .iter() + .map(|f| f.origins.iter().find_map(|o| if let FaceOrigin::Face { face, .. } = o { Some(*face) } else { None })) + .collect::>>() + .ok_or_else(|| failed(NAME, "a face of the result comes from no face of the body"))?; + // The new edges, remembered so that export keeps the pieces apart. + let mut lines: Vec> = body.split_lines().map(|l| l.to_vec()).unwrap_or_default(); + lines.extend(new_lines(&out, &h)?); + let mut out = out.with_split_lines(Some(std::sync::Arc::new(lines))).with_checked_history(Some(h), &record::counts(&[body])); + // An authored boundary: later rebuilds keep it. + out.regularization_protected = true; + Ok((out, map, made)) }) } -/// The field at a point well inside a face (a mesh triangle's centre where the field is -/// largest; for a plane, the face's corner farthest from it). -fn side_of(face: &mt::Face, g: &Field, size: f64) -> Result { - // A face the plane does not cross is on the side of whichever corner lies off it (the - // farthest, against rounding); a face with every corner on the plane is meshed. - if let Field::Plane { .. } = g { - let far = face - .boundaries() - .iter() - .flat_map(|w| w.vertex_iter()) - .map(|v| g.eval(from_p3(v.point()))) - .fold(0.0, |b: f64, x| if x.abs() > b.abs() { x } else { b }); - if far.abs() > size * 1e-6 { - return Ok(far); - } - } - let shell: mt::Shell = vec![face.clone()].into(); - let meshed = shell.robust_triangulation((size * 2e-3).max(1e-4)); - let pm = meshed.face_iter().next().and_then(|f| f.surface()).ok_or_else(|| KernelError::Failed("a face could not be meshed".into()))?; - let pos = pm.positions(); - let mut best: f64 = 0.0; - for t in pm.faces().triangle_iter() { - let c = t.iter().filter_map(|v| pos.get(v.pos)).fold(Vec3::ZERO, |a, p| a + from_p3(*p) / 3.0); - let x = g.eval(c); - if x.is_finite() && x.abs() > best.abs() { - best = x; - } - } - Ok(best) -} - -/// Loops from edges (each used once), chained end to start; edges are turned round as needed. -fn chain_loops(edges: &[mt::Edge]) -> Result>> { - let mut pool: Vec = edges.to_vec(); +/// The edges of `split` that no edge of the body was before (its split lines), as polylines +/// through their exact points. +fn new_lines(split: &Body, h: &History) -> Result>> { + let record = crate::native::brep_record(split)?; let mut out = Vec::new(); - while let Some(first) = pool.pop() { - let start = first.front().id(); - let mut at = first.back().clone(); - let mut lp = vec![first]; - while at.id() != start { - let Some(k) = pool.iter().position(|e| e.front().id() == at.id() || e.back().id() == at.id()) else { - return Err(KernelError::Failed("the section through the body does not close".into())); - }; - let e = pool.remove(k); - let e = if e.front().id() == at.id() { e } else { e.inverse() }; - at = e.back().clone(); - lp.push(e); - if lp.len() > 1_000_000 { - return Err(KernelError::Failed("section loop too long".into())); - } + for (e, prov) in record.edges.iter().zip(&h.edges) { + if prov.origins.iter().any(|o| matches!(o, EdgeOrigin::Edge { .. })) { + continue; + } + let c = CurveEval::new(&e.curve).map_err(|err| failed("split face", err))?; + let n = if matches!(e.curve, brep::Curve::Line { .. }) { 1 } else { 64 }; + let [t0, t1] = e.interval; + let pts: Option> = (0..=n).map(|k| c.point(t0 + (t1 - t0) * k as f64 / n as f64)).collect(); + if let Some(pts) = pts.filter(|p| p.iter().all(|q| q.is_finite())) { + out.push(pts); } - out.push(lp); } Ok(out) } -/// The regions of a body cut by a tool: for each face (in face order), the region it lies in, -/// and how many there are. Faces are in one region when they meet at an edge the tool does not -/// run along; an edge along the tool (a split line) parts them. -pub(crate) fn regions(split: &Body, tool: &SplitTool) -> Result<(Vec, usize)> { - let fields = fields(tool)?; - let tol = trace_tol_of(split.size()); - let faces: Vec<&mt::Face> = split.solid.face_iter().collect(); - let mut along_tool: HashMap = HashMap::new(); - let mut by_edge: HashMap> = HashMap::new(); - for (i, f) in faces.iter().enumerate() { - for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { - by_edge.entry(e.id()).or_default().push(i); - along_tool.entry(e.id()).or_insert_with(|| { - let c = e.curve(); - let (t0, t1) = c.range_tuple(); - let at: Vec = [0.25, 0.5, 0.75].iter().map(|s| from_p3(c.subs(t0 + (t1 - t0) * s))).collect(); - fields.iter().any(|g| at.iter().all(|p| g.eval(*p).abs() <= tol)) - }); - } - } - let mut region = vec![usize::MAX; faces.len()]; - let mut count = 0; - for s in 0..faces.len() { - if region.get(s) != Some(&usize::MAX) { - continue; - } - if let Some(r) = region.get_mut(s) { - *r = count; - } - let mut stack = vec![s]; - while let Some(i) = stack.pop() { - let Some(f) = faces.get(i) else { continue }; - for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { - if along_tool.get(&e.id()).copied().unwrap_or(false) { - continue; - } - for &j in by_edge.get(&e.id()).into_iter().flatten() { - if let Some(r) = region.get_mut(j) - && *r == usize::MAX - { - *r = count; - stack.push(j); - } - } - } - } - count += 1; - } - Ok((region, count)) -} - -/// How many fields a tool is (a sketch tool: one per curve). -pub(crate) fn field_count(tool: &SplitTool) -> Result { - Ok(fields(tool)?.len()) -} - -/// Faces grouped by shared edges (in face order). -fn components(faces: Vec) -> Vec> { - let n = faces.len(); - let mut by_edge: HashMap> = HashMap::new(); - for (i, f) in faces.iter().enumerate() { - for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { - by_edge.entry(e.id()).or_default().push(i); - } - } - let mut comp = vec![usize::MAX; n]; - let mut count = 0; - for s in 0..n { - if comp[s] != usize::MAX { - continue; - } - comp[s] = count; - let mut stack = vec![s]; - while let Some(i) = stack.pop() { - let Some(f) = faces.get(i) else { continue }; - for e in f.boundaries().iter().flat_map(|w| w.edge_iter()) { - for &j in by_edge.get(&e.id()).into_iter().flatten() { - if comp[j] == usize::MAX { - comp[j] = count; - stack.push(j); - } - } - } - } - count += 1; - } - let mut out: Vec> = vec![Vec::new(); count]; - for (f, c) in faces.into_iter().zip(comp) { - if let Some(v) = out.get_mut(c) { - v.push(f); - } - } - out -} - -/// Split a body into the parts on either side of a tool: a plane, or another face's surface -/// (extended). Faces are divided where the tool meets them and every part is closed with the -/// piece of the tool's surface inside the body; parts come in face order, those on the tool's -/// back first. A body the tool misses comes back whole. Each part's history is over `[body]`: +/// Split a body into the parts on either side of a tool: a plane, another face's surface +/// (extended), or one sketch line swept along its direction. Every part is closed with the piece +/// of the tool's surface inside the body; parts on the tool's back come first, each side's in +/// face order. A body the tool misses comes back whole. Each part's history is over `[body]`: /// its faces the body's (or pieces of them), the section new. pub fn split_body(body: &Body, tool: &SplitTool) -> Result> { - body.require_brep("split body")?; - if body.solid.boundaries().len() != 1 { - return Err(KernelError::Invalid("not supported yet: splitting a body with voids".into())); - } - let size = body.size(); - let fields = fields(tool)?; - let [g] = fields.as_slice() else { return Err(KernelError::Invalid("a body is split by a plane or a face".into())) }; - // The tool surface, with the field's normal. - // The tool's whole surface for the section (a sheet that is part of a plane or cylinder - // is extended to all of it), with the field's normal. - let surface = match g { - Field::Plane { o, n } => { - let x = (if n.x.abs() < 0.9 { Vec3::X } else { Vec3::Y }).cross(*n).normalized().unwrap_or(Vec3::X); - let y = n.cross(x); - let o = p3(*o); - mt::Surface::Plane(mt::Plane::new(o, o + crate::body::v3(x), o + crate::body::v3(y))) - } - Field::Cylinder { o, axis, r, out } => { - // A line along the axis at radius r, long enough for the body, revolved. - let mut lo = f64::INFINITY; - let mut hi = f64::NEG_INFINITY; - for v in body.solid.vertex_iter() { - let t = (from_p3(v.point()) - *o).dot(*axis); - lo = lo.min(t); - hi = hi.max(t); - } - let (lo, hi) = (lo - size, hi + size); - let x = (if axis.x.abs() < 0.9 { Vec3::X } else { Vec3::Y }).cross(*axis).normalized().unwrap_or(Vec3::X); - let (a, b) = (*o + *axis * lo + x * *r, *o + *axis * hi + x * *r); - let line = mt::Curve::Line(mt::Line(p3(a), p3(b))); - crate::step_in::geom::revolved_oriented(line, p3(*o), crate::body::v3(axis.vec()), 0.5, crate::body::v3(x * *out)) - } - Field::Surface(s) => s.clone(), - Field::Curve { .. } => return Err(KernelError::Invalid("a body is split by a plane or a face".into())), - }; - let total = body.solid.face_iter().count(); - let all: Vec = (0..total).collect(); - guard("split body", || { - let imp = imprint(body, &all, &fields, size)?; - let new_edges = imp.new_edges; - if imp.made.values().all(|m| *m == 0) || new_edges.is_empty() { - return Ok(vec![crate::boolean_record::operand(body, 0, &crate::record::counts(&[body]))]); - } - let mut rec = imp.rec; - let faces: Vec = imp.shells.into_iter().flatten().collect(); - // The section: loops of the new edges on the tool surface, outer loops counter-clockwise - // about its normal and holes inside them. - let loops = chain_loops(&new_edges)?; - let polys: Vec> = loops.iter().map(|l| loop_uv(&surface, &l.clone().into())).collect(); - let signed = |p: &[(f64, f64)]| { - (0..p.len()) - .map(|i| { - let (a, b) = (p[i], p[(i + 1) % p.len()]); - a.0 * b.1 - b.0 * a.1 - }) - .sum::() - }; - let depth: Vec = (0..loops.len()) - .map(|i| { - let q = polys.get(i).and_then(|p| p.first()).copied().unwrap_or((0.0, 0.0)); - (0..loops.len()).filter(|&j| j != i && polys.get(j).is_some_and(|p| inside(p, q))).count() - }) - .collect(); - let wire = |i: usize, ccw: bool| -> mt::Wire { - let w: mt::Wire = loops.get(i).cloned().unwrap_or_default().into(); - let is_ccw = polys.get(i).is_some_and(|p| signed(p) > 0.0); - if is_ccw == ccw { w } else { w.inverse() } - }; - let mut section = Vec::new(); - // On a surface of revolution, loops that run once round it (a cylinder through a box: - // its circles at the top and bottom) bound bands: taken in order along the profile, - // each pair is one face, the lower loop running back round and the upper forward. - let mut wrapping = vec![false; loops.len()]; - if let mt::Surface::RevolutedCurve(rc) = &surface { - let (o, axis) = (rc.entity().origin(), rc.entity().axis().normalize()); - let x0 = if axis.x.abs() < 0.9 { mt::Vector3::unit_x() } else { mt::Vector3::unit_y() }; - let ex = (x0 - axis * x0.dot(axis)).normalize(); - let ey = axis.cross(ex); - let tau = std::f64::consts::TAU; - let samples = |w: &mt::Wire| -> Vec { - w.edge_iter() - .flat_map(|e| { - let c = e.oriented_curve(); - let (t0, t1) = c.range_tuple(); - (0..32).map(move |k| c.subs(t0 + (t1 - t0) * k as f64 / 32.0)).collect::>() - }) - .collect() - }; - let mut bands: Vec<(f64, usize)> = Vec::new(); - for (i, l) in loops.iter().enumerate() { - let pts = samples(&l.clone().into()); - let ang = |p: &mt::Point3| (p - o).dot(ey).atan2((p - o).dot(ex)); - let turn: f64 = (0..pts.len()) - .map(|k| { - let d = ang(&pts[(k + 1) % pts.len()]) - ang(&pts[k]); - d - tau * (d / tau).round() - }) - .sum(); - if (turn.abs() - tau).abs() < 0.1 { - let along = pts.iter().map(|p| (p - o).dot(axis)).sum::() / pts.len().max(1) as f64; - bands.push((along, i)); - if let Some(w) = wrapping.get_mut(i) { - *w = true; - } - } - } - bands.sort_by(|a, b| a.0.total_cmp(&b.0)); - // A loop runs so that the band (towards `toward` along the axis) is on its left about - // the surface normal. - let oriented = |i: usize, toward: f64| -> mt::Wire { - let w: mt::Wire = loops.get(i).cloned().unwrap_or_default().into(); - let pts = samples(&w); - let (Some(p), Some(q)) = (pts.first(), pts.get(1)) else { return w }; - let n = surface.search_nearest_parameter(*p, None, 50).map(|(u, v)| surface.normal(u, v)).unwrap_or(axis); - let left = n.cross(q - p); - if left.dot(axis) * toward >= 0.0 { w } else { w.inverse() } - }; - for pair in bands.chunks(2) { - let [(_, lo), (_, hi)] = pair else { return Err(KernelError::Failed("an unmatched band edge in the section".into())) }; - section.push(mt::Face::new_unchecked(vec![oriented(*lo, 1.0), oriented(*hi, -1.0)], surface.clone())); - } - } - for (i, d) in depth.iter().enumerate() { - if d % 2 != 0 || wrapping.get(i).copied().unwrap_or(false) { + const NAME: &str = "split body"; + body.require_brep(NAME)?; + let reach = reach_of(body); + let groups = cutters(tool, &reach)?; + let [group] = groups.as_slice() else { return Err(invalid("a body is split by a plane or a face")) }; + let [patch] = group.as_slice() else { return Err(invalid("a body is split by a plane, a face or one sketch line")) }; + let counts = record::counts(&[body]); + let whole = || crate::boolean_record::operand(body, 0, &counts); + let Some(surface) = smlib_tools(std::slice::from_ref(group), &reach)?.pop() else { return Ok(vec![whole()]) }; + guard(NAME, || { + let (sb, ob) = tagged(NAME, body)?; + let Some((back, front)) = sb.split_body(&surface).map_err(|e| failed(NAME, e))? else { return Ok(vec![whole()]) }; + // SMLib's sides are the surface's; the tool's are the face's. + let sides = if patch.flip { [front, back] } else { [back, front] }; + let (mut parts, mut cut) = (Vec::new(), false); + for side in &sides { + let (part, h) = brought_back(NAME, body, &ob, side)?; + // The section comes from no face of the body. + cut |= h.faces.iter().any(|f| f.origins.is_empty()); + let part = part.with_checked_history(Some(h), &counts); + let lumps = part.lumps()?; + if lumps.len() < 2 { + parts.push(part); continue; } - let mut ws = vec![wire(i, true)]; - for (j, dj) in depth.iter().enumerate() { - let q = polys.get(j).and_then(|p| p.first()).copied().unwrap_or((0.0, 0.0)); - if *dj == d + 1 && polys.get(i).is_some_and(|p| inside(p, q)) { - ws.push(wire(j, false)); - } - } - section.push(mt::Face::new_unchecked(ws, surface.clone())); - } - // Each face goes to the side its inside is on. - let (mut back, mut front) = (Vec::new(), Vec::new()); - for f in faces { - let x = side_of(&f, g, size)?; - if x.abs() <= size * 1e-9 { - return Err(KernelError::Invalid("not supported yet: the tool lies along a face of the body".into())); - } - if x < 0.0 { - back.push(f); - } else { - front.push(f); + for lump in lumps { + let h = crate::boolean_record::through(&lump, 0, part.history()); + parts.push(lump.with_checked_history(h, &counts)); } } - // The section is new: it comes from nothing of the body. - for f in §ion { - rec.made(f, Vec::new()); - } - // The section closes the back side facing the tool's normal, the front side facing away. - back.extend(section.iter().cloned()); - front.extend(section.iter().map(|f| f.inverse())); - let mut parts = Vec::new(); - for side in [back, front] { - for comp in components(side) { - let shell = mt::Shell::from(comp); - if shell.shell_condition() != mt::ShellCondition::Closed { - return Err(KernelError::Failed("a part of the split does not close".into())); - } - parts.push(crate::record::rebuilt(mt::Solid::new_unchecked(vec![shell]), &[&body.solid], &rec)?.with_color(body.color)); - } - } - Ok(parts) + Ok(if cut { parts } else { vec![whole()] }) }) } - #[cfg(test)] mod tests { use super::*; @@ -1688,7 +814,8 @@ mod tests { let quarter = std::f64::consts::PI * 36.0 / 4.0; // (Face areas come from meshes.) assert!(rel(tops[0], quarter) < 1e-3 && rel(tops[1], 200.0 - quarter) < 1e-3, "{tops:?}"); - assert!(s.solid.edge_iter().any(|e| matches!(e.curve(), mt::Curve::NurbsCurve(_)))); + let record = crate::native::brep_record(&s).unwrap(); + assert!(record.edges.iter().any(|e| matches!(e.curve, brep::Curve::Circle { .. } | brep::Curve::RationalSpline { .. }))); // A listed face the tool misses is named: the plane z = 10 crosses the sides, not the top. let mid = Plane { origin: Vec3::new(0.0, 0.0, 10.0), ..Plane::XY }; let side = b.faces(0.01).unwrap().iter().position(|f| (f.centroid.x - 10.0).abs() < 1e-6).unwrap(); diff --git a/crates/smlib/src/ssi.rs b/crates/smlib/src/ssi.rs index 628b6dece7..cd31c95e14 100644 --- a/crates/smlib/src/ssi.rs +++ b/crates/smlib/src/ssi.rs @@ -86,6 +86,59 @@ impl Surface { ffi::ssi::nurbs([int(degree[0])?, int(degree[1])?], int(nu)?, int(nv)?, knots_u, knots_v, &flat, weights).map(Surface).map_err(error) } + /// A schema 4 record surface over the parameter box `domain` (`[[u0, u1], [v0, v1]]`, the + /// record's own parameters: radians for angles), with the record's parameterisation and + /// normal. A spline is taken over its whole knot domain. A cone only as a whole turn, widening + /// along its axis; a revolution not at all ([`Error::Unsupported`]). + pub fn from_record(surface: &solidcraft_geom::brep::Surface, domain: [[f64; 2]; 2]) -> Result { + use solidcraft_geom::brep::Surface as S; + let [u, v] = domain; + let deg = |r: [f64; 2]| r.map(f64::to_degrees); + let sub = |a: [f64; 3], b: [f64; 3]| [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; + let len = |a: [f64; 3]| (a[0] * a[0] + a[1] * a[1] + a[2] * a[2]).sqrt(); + let cross = |a: [f64; 3], b: [f64; 3]| [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]; + // SMLib's plane runs along its axes made unit: the record's parameters scale by their lengths. + let plane = |origin: [f64; 3], du: [f64; 3], dv: [f64; 3]| { + let (a, b) = (len(du), len(dv)); + Surface::plane(origin, du, dv, [u[0] * a, u[1] * a], [v[0] * b, v[1] * b]) + }; + match surface { + S::Plane { origin, u_axis, v_axis } => plane(*origin, *u_axis, *v_axis), + S::PlaneThrough { origin, u_point, v_point } => plane(*origin, sub(*u_point, *origin), sub(*v_point, *origin)), + S::Cylinder { frame, radius } => Surface::cylinder(frame.origin, frame.x_axis, frame.y_axis, *radius, deg(u), v), + S::Cone { frame, radius_at_origin, half_angle } => { + let t = half_angle.tan(); + if !(t > 0.0) || !(u[1] - u[0] >= std::f64::consts::TAU * (1.0 - 1e-12)) { + return Err(Error::Unsupported); + } + // The apex lies where the radius is 0; heights run from it along the axis. + let at = -radius_at_origin / t; + let z = cross(frame.x_axis, frame.y_axis); + let apex = [0, 1, 2].map(|k| frame.origin[k] + z[k] * at); + Surface::cone(apex, frame.x_axis, frame.y_axis, half_angle.to_degrees(), [(v[0] - at).max(0.0), v[1] - at]) + } + S::Sphere { frame, radius } => Surface::sphere(frame.origin, frame.x_axis, frame.y_axis, *radius, deg(u), deg(v)), + S::Torus { frame, major_radius, minor_radius } => { + Surface::torus(frame.origin, frame.x_axis, frame.y_axis, *major_radius, *minor_radius, deg(u), deg(v)) + } + S::Spline { u_knots, v_knots, points } => { + let (nu, nv) = (points.len(), points.first().map_or(0, Vec::len)); + let flat: Vec<[f64; 3]> = points.iter().flatten().copied().collect(); + let degree = [u_knots.len().checked_sub(nu + 1).ok_or(Error::Argument)?, v_knots.len().checked_sub(nv + 1).ok_or(Error::Argument)?]; + Surface::nurbs(degree, nu, nv, u_knots, v_knots, &flat, None) + } + S::RationalSpline { u_knots, v_knots, homogeneous } => { + let (nu, nv) = (homogeneous.len(), homogeneous.first().map_or(0, Vec::len)); + let all: Vec<[f64; 4]> = homogeneous.iter().flatten().copied().collect(); + let flat: Vec<[f64; 3]> = all.iter().map(|h| [h[0] / h[3], h[1] / h[3], h[2] / h[3]]).collect(); + let weights: Vec = all.iter().map(|h| h[3]).collect(); + let degree = [u_knots.len().checked_sub(nu + 1).ok_or(Error::Argument)?, v_knots.len().checked_sub(nv + 1).ok_or(Error::Argument)?]; + Surface::nurbs(degree, nu, nv, u_knots, v_knots, &flat, Some(&weights)) + } + S::PlacedRevolution { .. } => Err(Error::Unsupported), + } + } + /// The UV domain the intersection runs within, as `[u0, u1, v0, v1]`. pub fn domain(&self) -> Result<[f64; 4]> { ffi::ssi::domain(&self.0).map_err(error)