From 34be218639556700756afbf3764308f9bf4c095d Mon Sep 17 00:00:00 2001 From: bherbruck Date: Thu, 8 Oct 2026 16:26:48 +0000 Subject: [PATCH] Kernel: split a body by a plane or another face's surface split_body(body, tool) divides every face where the tool meets it and closes each side with the piece of the tool surface inside the body. - Patches have outer loops and holes. - Bands go round a cylinder or another surface of revolution. - Each side falls apart into its connected parts. - A body the tool misses comes back whole. Faces that are B-spline sheets but really planes or circular cylinders (extruded faces) split as the whole plane or cylinder: the tool is extended, and the field is exact and cheap. Tests (volumes): a box split by a plane (2400/3600); a holed plate split through its hole into two equal halves; a box split by a cylinder's side into the core (pi*16*30) and the rest; a plane that misses. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL --- crates/kernel/src/lib.rs | 2 +- crates/kernel/src/splitface.rs | 366 +++++++++++++++++++++++++++++- crates/kernel/src/step_in/geom.rs | 2 +- 3 files changed, 364 insertions(+), 6 deletions(-) diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index e6c88946..4b7480a8 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -52,7 +52,7 @@ pub use ops::{BoolOp, boolean, split_by_plane, transform, transform_matrix}; pub use polybool::planar_boolean; pub use polyhedron::{HalfSpace, convex_polyhedron, draft, offset_faces, shell}; pub use sew::{EdgeSpec, FaceSpec, SurfSpec, sew}; -pub use splitface::{SplitTool, split_faces}; +pub use splitface::{SplitTool, split_body, split_faces}; pub use step::step_export; pub use step_in::{ImportedBody, StepImport, StepNode, step_import, step_import_shared, step_orientation_errors, step_validate}; pub use step_out::{ExportBody, ExportProduct, StepHeader, step_export_bodies, step_export_products}; diff --git a/crates/kernel/src/splitface.rs b/crates/kernel/src/splitface.rs index 694765bb..9e0ae396 100644 --- a/crates/kernel/src/splitface.rs +++ b/crates/kernel/src/splitface.rs @@ -30,9 +30,23 @@ pub enum SplitTool<'a> { /// Which side of the tool a point is on (zero on the tool). enum Field { - Plane { o: Vec3, n: Vec3 }, + Plane { + o: Vec3, + n: Vec3, + }, Surface(mt::Surface), - Curve { plane: Plane, pts: Vec }, + /// 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: Vec3, + r: f64, + out: f64, + }, + Curve { + plane: Plane, + pts: Vec, + }, } impl Field { @@ -47,6 +61,11 @@ impl Field { } 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 } => { let q = plane.to_local(p); let q = Vec2::new(q.x, q.y); @@ -214,7 +233,11 @@ fn trace(face: &mt::Face, g: &Field, size: f64) -> Result> { hi = mid; } } - let uv = f(0.5 * (lo + hi)).1; + 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. @@ -629,13 +652,67 @@ fn split_one(faces: &mut Vec, fi: usize, g: &Field, size: f64, new_edg Ok((made, count)) } +/// 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.normalized() 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()) +} + 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![Field::Surface(f.oriented_surface())] + vec![surface_field(&f.oriented_surface())] } SplitTool::Curves { plane, curves } => { if curves.iter().any(|c| c.len() < 2 || c.iter().any(|p| !(p.x.is_finite() && p.y.is_finite()))) { @@ -700,6 +777,256 @@ pub fn split_faces(body: &Body, faces: &[usize], tool: &SplitTool) -> Result Result { + 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(); + 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())); + } + } + out.push(lp); + } + Ok(out) +} + +/// 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. +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), 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 (shells, new_edges, made) = imprint(body, &all, &fields, size)?; + if made == 0 || new_edges.is_empty() { + return Ok(vec![body.clone()]); + } + let faces: Vec = 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) { + 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); + } + } + // 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(Body::new(mt::Solid::new_unchecked(vec![shell]))?.with_color(body.color)); + } + } + Ok(parts) + }) +} + #[cfg(test)] mod tests { use super::*; @@ -783,6 +1110,37 @@ mod tests { assert!(tops.iter().all(|a| rel(*a, half) < 1e-3), "{tops:?}"); } + fn volumes(parts: &[Body]) -> Vec { + let mut v: Vec = parts.iter().map(|b| measure(b).unwrap().volume).collect(); + v.sort_by(f64::total_cmp); + v + } + + #[test] + fn split_body_by_plane_and_by_a_face() { + let b = box_solid(Vec3::ZERO, Vec3::new(10.0, 20.0, 30.0)).unwrap(); + let pl = Plane::new(Vec3::new(4.0, 0.0, 0.0), Vec3::Y, Vec3::Z).unwrap(); + let v = volumes(&split_body(&b, &SplitTool::Plane(pl)).unwrap()); + assert!(v.len() == 2 && rel(v[0], 2400.0) < 1e-9 && rel(v[1], 3600.0) < 1e-9, "{v:?}"); + // A plate with a hole split through the hole: two equal halves. + let plate = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 10.0)).unwrap(); + let pin = cylinder(Vec3::new(20.0, 15.0, -1.0), Vec3::Z, 5.0, 12.0).unwrap(); + let holed = crate::boolean(&plate, &pin, crate::BoolOp::Cut).unwrap().unwrap(); + let pl = Plane::new(Vec3::new(20.0, 0.0, 0.0), Vec3::Y, Vec3::Z).unwrap(); + let v = volumes(&split_body(&holed, &SplitTool::Plane(pl)).unwrap()); + let half = (12000.0 - std::f64::consts::PI * 250.0) / 2.0; + assert!(v.len() == 2 && rel(v[0], half) < 1e-4 && rel(v[1], half) < 1e-4, "{v:?}"); + // A box split by a cylinder's side: the core and the rest. + let tool = cylinder(Vec3::new(5.0, 10.0, -5.0), Vec3::Z, 4.0, 40.0).unwrap(); + let side = tool.faces(0.01).unwrap().iter().position(|f| f.plane_normal.is_none()).unwrap(); + let v = volumes(&split_body(&b, &SplitTool::Face { body: &tool, face: side }).unwrap()); + let core = std::f64::consts::PI * 16.0 * 30.0; + assert!(v.len() == 2 && rel(v[0], core) < 1e-3 && rel(v[1], 6000.0 - core) < 1e-3, "{v:?}"); + // A plane that misses gives the body back. + let far = Plane::new(Vec3::new(50.0, 0.0, 0.0), Vec3::Y, Vec3::Z).unwrap(); + assert_eq!(split_body(&b, &SplitTool::Plane(far)).unwrap().len(), 1); + } + #[test] fn a_tool_that_misses_is_an_error() { let b = box_solid(Vec3::ZERO, Vec3::new(10.0, 20.0, 30.0)).unwrap(); diff --git a/crates/kernel/src/step_in/geom.rs b/crates/kernel/src/step_in/geom.rs index cb86728a..af767dfb 100644 --- a/crates/kernel/src/step_in/geom.rs +++ b/crates/kernel/src/step_in/geom.rs @@ -684,7 +684,7 @@ fn revolved(curve: mt::Curve, o: mt::Point3, axis: mt::Vector3) -> mt::Surface { /// A surface of revolution whose normal at the profile point `(u, 0)` points along `expected`. /// The rotation is reversed rather than the surface inverted: truck's inverted `Processor` /// swaps the parameters but not the search hints, which breaks meshing. -fn revolved_oriented(curve: mt::Curve, o: mt::Point3, axis: mt::Vector3, u: f64, expected: mt::Vector3) -> mt::Surface { +pub(crate) fn revolved_oriented(curve: mt::Curve, o: mt::Point3, axis: mt::Vector3, u: f64, expected: mt::Vector3) -> mt::Surface { let s = revolved(curve.clone(), o, axis); if s.normal(u, 0.0).dot(expected) >= 0.0 { return s;