From a4a17331dc8b7e27c94d1c494f1235d5990dcac6 Mon Sep 17 00:00:00 2001 From: bherbruck Date: Wed, 7 Oct 2026 22:08:19 +0000 Subject: [PATCH] Kernel: concave edge fillets; whole smooth-loop blends of planar faces (cylinders and tori: pocket floors, plate outlines); rounding every edge of a convex planar body (sphere corners); healing after booleans (exact lines/arcs for intersection curves, coplanar faces merged, collinear edges joined); faces rebuilt on their stored surfaces Oracle: 28/29. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01AS5Tn3xieewUSg3uxrMVdL --- ROADMAP.md | 7 +- crates/kernel/src/blend.rs | 128 ++++++++--- crates/kernel/src/build.rs | 2 +- crates/kernel/src/heal.rs | 361 ++++++++++++++++++++++++++++++++ crates/kernel/src/lib.rs | 2 + crates/kernel/src/loopblend.rs | 306 +++++++++++++++++++++++++++ crates/kernel/src/ops.rs | 2 +- crates/kernel/src/polyhedron.rs | 140 ++++++++++++- crates/kernel/src/tests.rs | 92 ++++++++ 9 files changed, 997 insertions(+), 43 deletions(-) create mode 100644 crates/kernel/src/heal.rs create mode 100644 crates/kernel/src/loopblend.rs diff --git a/ROADMAP.md b/ROADMAP.md index c470a2e1..ae66a88d 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -9,7 +9,7 @@ features → parametric timeline) as an open, pure-Rust, agent-drivable applicat |---|---| | Where are we? | **M0 done; M1/M2 in progress.** Sketch with solver, extrude (taper, two-sided, through all), revolve, fillet/chamfer, holes, patterns, mirror, shell, draft, loft, sweep, split, construction planes, parameters that rebuild the timeline, STEP/STL export, desktop and browser (wasm) app with Fusion-style selection (hover/selected highlights, origin planes, click-to-pick sketch planes, multi-select, window/crossing box selection, animated view cube), CLI, MCP server, oracle harness. | | Command parity (SOLID + SKETCH toolbar) | **57 / 291 (20%)** — [docs/parity.md](docs/parity.md) | -| Fusion oracle | **25 / 29 parts match** Fusion's measurements — [docs/oracle.md](docs/oracle.md) | +| Fusion oracle | **28 / 29 parts match** Fusion's measurements — [docs/oracle.md](docs/oracle.md) | | Tests | 85 (solver, profiles, kernel booleans/blends/measures, expressions, timeline, file formats, camera, engine end-to-end, hostile-input fuzz over every command) | | Gates | `cargo xtask ci`: fmt, clippy -D warnings, tests, asset attribution, layering, wasm32 build — green | | Weighted parity estimate | **≈ 5%** of Fusion's Design workspace by importance (sketch + basic solids are the core, but surfaces, assemblies, sheet metal, CAM, drawings are untouched) | @@ -35,8 +35,9 @@ features → parametric timeline) as an open, pure-Rust, agent-drivable applicat ## Known limitations -- Fillets and chamfers: straight convex edges between planar faces with perpendicular planar end - faces; two blended edges may not share a corner yet. +- Fillets and chamfers: straight edges (convex or concave) between planar faces with perpendicular + planar end faces; whole smooth loops of a planar face (pocket floors, plate outlines); every edge + of a convex planar body (sphere corners). Other corner configurations are not supported yet. - Booleans of curved bodies with coincident faces can fail (planar ones have an exact fallback); fully internal voids are not supported. Spheres are built from six pole-free patches so they combine reliably. diff --git a/crates/kernel/src/blend.rs b/crates/kernel/src/blend.rs index a81b7f94..423c900c 100644 --- a/crates/kernel/src/blend.rs +++ b/crates/kernel/src/blend.rs @@ -44,7 +44,7 @@ fn plane_normal(f: &mt::Face) -> Option { /// edges where consecutive edges no longer connect. fn rebuild_face(f: &mt::Face, subst: &HashMap, connectors: &[mt::Edge]) -> Result { let mut wires = Vec::new(); - for w in f.boundaries() { + for w in f.absolute_boundaries() { let mut edges: Vec = Vec::new(); for e in w.edge_iter() { let ne = match subst.get(&e.id()) { @@ -81,8 +81,7 @@ fn rebuild_face(f: &mt::Face, subst: &HashMap, connectors: } wires.push(mt::Wire::from(out)); } - let surface = f.oriented_surface(); - mt::Face::try_new(wires, surface).map_err(|e| KernelError::Failed(format!("blend: {e}"))) + crate::heal::absolute_face(f, wires).ok_or_else(|| KernelError::Failed("blend: a face could not be rebuilt".into())) } fn blend_one(solid: &Solid, edge: &mt::Edge, size: f64, r: f64, shape: Shape) -> Result { @@ -103,22 +102,25 @@ fn blend_one(solid: &Solid, edge: &mt::Edge, size: f64, r: f64, shape: Shape) -> let (v0, v1) = (edge.absolute_front().clone(), edge.absolute_back().clone()); let (p0, p1) = (vtx(&v0), vtx(&v1)); let d = (p1 - p0).normalized().ok_or_else(|| KernelError::Failed("zero-length edge".into()))?; - // In-face directions perpendicular to the edge, pointing away from it. - let mut t1 = n1.cross(d); - if t1.dot(n2) > 0.0 { - t1 = -t1; - } - let mut t2 = n2.cross(d); - if t2.dot(n1) > 0.0 { - t2 = -t2; - } - // Convex edge: the faces' outward normals point away from each other's interior. + // In-face directions perpendicular to the edge, into each face: a face's boundary runs + // counter-clockwise about its outward normal, so its interior is to the left. + let into = |f: &mt::Face, n: Vec3| -> Option { + let e = f.boundary_iters().into_iter().flatten().find(|e| e.id() == edge.id())?; + let dir = (vtx(e.back()) - vtx(e.front())).normalized()?; + n.cross(dir).normalized() + }; + let (Some(t1), Some(t2)) = (into(f1, n1), into(f2, n2)) else { return Err(KernelError::Failed("blend: edge direction".into())) }; + // Convex edge: each face's outward normal points away from the other face; concave: toward. let cos_phi = t1.dot(t2).clamp(-1.0, 1.0); let phi = cos_phi.acos(); - if !(phi > 1e-3 && phi < std::f64::consts::PI - 1e-3) || n1.dot(t2) > 1e-9 || n2.dot(t1) > 1e-9 { - return Err(unsupported("only convex edges can be blended")); + let convex = n1.dot(t2) < -1e-9 && n2.dot(t1) < -1e-9; + let concave = n1.dot(t2) > 1e-9 && n2.dot(t1) > 1e-9; + if !(phi > 1e-3 && phi < std::f64::consts::PI - 1e-3) || !(convex || concave) { + return Err(unsupported("the faces at the edge are tangent or the edge is ambiguous")); } - blend_geometry(solid, si, &faces, (i1, i2), edge, (v0, v1), (p0, p1, d), (n1, t1, t2, phi), size, r, shape) + // The blend's centre line is inside the material for a convex edge, outside for a concave one. + let side = if convex { -1.0 } else { 1.0 }; + blend_geometry(solid, si, &faces, (i1, i2), edge, (v0, v1), (p0, p1, d), (n1, t1, t2, phi, side), size, r, shape) } /// Distance along each face from the edge to the blend boundary. @@ -138,7 +140,7 @@ fn blend_geometry( edge: &mt::Edge, (v0, v1): (mt::Vertex, mt::Vertex), (p0, p1, d): (Vec3, Vec3, Vec3), - (n1, t1, t2, phi): (Vec3, Vec3, Vec3, f64), + (n1, t1, t2, phi, side): (Vec3, Vec3, Vec3, f64, f64), size: f64, r: f64, shape: Shape, @@ -190,16 +192,17 @@ fn blend_geometry( let b0 = builder::vertex(p3(p0 + t2 * s)); let a1 = builder::vertex(p3(p1 + t1 * s)); let b1 = builder::vertex(p3(p1 + t2 * s)); - // Corner curves on the end faces (A → B) and the blend surface. - let bisector = (t1 + t2).normalized().ok_or_else(|| KernelError::Failed("bisector".into()))?; + // Corner curves on the end faces (A → B) and the blend surface. The arc centre is r from + // both faces (inside the material for a convex edge, outside for a concave one); the arc's + // middle is the point nearest the old edge. + let centre = |p: Vec3| p + t1 * s + n1 * (r * side); + let arc_mid = |p: Vec3| { + let c = centre(p); + c + (p - c).normalized().unwrap_or(-n1) * r + }; let corner = |a: &mt::Vertex, b: &mt::Vertex, p: Vec3| -> mt::Edge { match shape { - Shape::Round => { - // Arc centre is r inside both faces; its mid point is on the bisector, toward the edge. - let c = p + t1 * s - n1 * r; - let mid = c - bisector * r; - builder::circle_arc(a, b, p3(mid)) - } + Shape::Round => builder::circle_arc(a, b, p3(arc_mid(p))), Shape::Flat => builder::line(a, b), } }; @@ -248,13 +251,14 @@ fn blend_geometry( Shape::Round => { let swept: mt::Face = builder::tsweep(&c0, v3(p1 - p0)); let mut surf = swept.oriented_surface(); - // Outward normal at the arc middle points away from the arc centre. - let c = p0 + t1 * s - n1 * r; - let mid = c - bisector * r; + // Outward normal at the arc middle: away from the centre on a convex edge (the + // centre is in the material), toward it on a concave one. + let (c, mid) = (centre(p0), arc_mid(p0)); + let outward = (mid - c) * (-side); use mt::{ParametricSurface3D, SearchParameter}; if let Some((u, v)) = surf.search_parameter(p3(mid), None, 100) { let nn = surf.normal(u, v); - if Vec3::new(nn.x, nn.y, nn.z).dot(bisector) > 0.0 { + if Vec3::new(nn.x, nn.y, nn.z).dot(outward) < 0.0 { surf = mt::Invertible::inverse(&surf); } } @@ -280,7 +284,48 @@ fn blend(body: &Body, edges: &[Vec3], r: f64, shape: Shape, what: &str) -> Resul if edges.len() > 1000 { return Err(KernelError::Invalid("too many edges".into())); } - let mut cur = body.clone(); + // Coplanar neighbours as one face (as after a join), so corners look as they should. + let mut cur = Body::new(crate::heal::heal(body.deep_copy(), body.size()))?; + // Every edge of a convex planar body: rounded as a whole (sphere corners). + if shape == Shape::Round && edges.len() >= 6 { + let solid = cur.deep_copy(); + let all = Body::unique_edges(&solid); + if all.len() == edges.len() { + let size = cur.size(); + let tol = (size * 2e-3).max(1e-3); + let mut hit = std::collections::HashSet::new(); + for p in edges { + if let Some((idx, _)) = cur.nearest_edge(*p, tol)? { + hit.insert(idx); + } + } + if hit.len() == all.len() + && let Ok(b) = crate::polyhedron::round_all_edges(&cur, r) + { + return Ok(b); + } + } + } + // A whole smooth loop of a planar face is blended in one go. + if edges.len() >= 2 { + let size = cur.size(); + let tol = (size * 2e-3).max(1e-3); + let solid = cur.deep_copy(); + let all = Body::unique_edges(&solid); + let mut ids = Vec::new(); + for p in edges { + if let Some((idx, dist)) = cur.nearest_edge(*p, tol)? + && dist <= size * 0.05 + 1e-3 + && let Some(e) = all.get(idx) + && !ids.contains(&e.id()) + { + ids.push(e.id()); + } + } + if let Some(r2) = crate::loopblend::loop_blend(&solid, &ids, size, r, shape == Shape::Round) { + return Body::new(guard(what, || r2)?); + } + } for p in edges { let size = cur.size(); let tol = (size * 2e-3).max(1e-3); @@ -336,7 +381,28 @@ pub fn chamfer(body: &Body, edges: &[Vec3], distance: f64) -> Result { return Ok(r); } } - let mut cur = body.clone(); + // Coplanar neighbours as one face (as after a join), so corners look as they should. + let mut cur = Body::new(crate::heal::heal(body.deep_copy(), body.size()))?; + // A whole smooth loop of a planar face is blended in one go. + if edges.len() >= 2 { + let size = cur.size(); + let tol = (size * 2e-3).max(1e-3); + let solid = cur.deep_copy(); + let all = Body::unique_edges(&solid); + let mut ids = Vec::new(); + for p in edges { + if let Some((idx, dist)) = cur.nearest_edge(*p, tol)? + && dist <= size * 0.05 + 1e-3 + && let Some(e) = all.get(idx) + && !ids.contains(&e.id()) + { + ids.push(e.id()); + } + } + if let Some(r2) = crate::loopblend::loop_blend(&solid, &ids, size, distance, false) { + return Body::new(guard("chamfer", || r2)?); + } + } for p in edges { cur = match chamfer_tool(&cur, *p, distance, 0).and_then(|t| crate::ops::boolean(&cur, &t, crate::BoolOp::Cut)) { Ok(Some(b)) => b, diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index bfa4018d..b0d1a0a3 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -609,7 +609,7 @@ fn cubic_basis(knots: &[f64], n_ctrl: usize, u: f64) -> Vec { } /// Cubic B-spline through points at uniform parameters (clamped, knots by averaging). -fn interpolate_cubic(pts: &[Vec3]) -> Option> { +pub(crate) fn interpolate_cubic(pts: &[Vec3]) -> Option> { let m = pts.len(); if m < 4 { return None; diff --git a/crates/kernel/src/heal.rs b/crates/kernel/src/heal.rs new file mode 100644 index 00000000..b40e7fcf --- /dev/null +++ b/crates/kernel/src/heal.rs @@ -0,0 +1,361 @@ +//! Clean-up after operations: coplanar neighbouring faces become one face, and two straight +//! edges meeting in line at a vertex nothing else uses become one edge (as Fusion keeps them). + +use std::collections::HashMap; + +use solvecraft_geom::Vec3; +use truck_modeling::{self as mt, builder}; + +use crate::body::{Solid, from_p3}; + +fn plane_of(f: &mt::Face) -> Option<(Vec3, f64)> { + match f.oriented_surface() { + mt::Surface::Plane(p) => { + let n = p.normal(); + let n = Vec3::new(n.x, n.y, n.z).normalized()?; + let o = p.origin(); + Some((n, n.dot(Vec3::new(o.x, o.y, o.z)))) + } + _ => None, + } +} + +/// A face like `f` with new boundary wires given in `f`'s absolute (stored) orientation: +/// built on the stored surface and turned like `f` (rebuilding on the turned surface can make +/// the mesher pick the wrong side of a trim on periodic surfaces). +pub(crate) fn absolute_face(f: &mt::Face, wires: Vec) -> Option { + let mut nf = mt::Face::try_new(wires, f.surface()).ok()?; + if !f.orientation() { + nf.invert(); + } + Some(nf) +} + +fn find(p: &mut [usize], mut i: usize) -> usize { + while let Some(&q) = p.get(i) { + if q == i { + break; + } + i = q; + } + i +} + +/// Points along an oriented edge (for loop areas and containment). +fn samples(e: &mt::Edge) -> Vec { + use mt::{BoundedCurve, ParametricCurve}; + let c = e.oriented_curve(); + let (t0, t1) = c.range_tuple(); + (0..8).map(|k| t0 + (t1 - t0) * k as f64 / 8.0).map(|t| from_p3(c.subs(t))).collect() +} + +fn loop_area(edges: &[mt::Edge], n: Vec3) -> f64 { + let pts: Vec = edges.iter().flat_map(samples).collect(); + let k = pts.len(); + let mut s = Vec3::ZERO; + for i in 0..k { + if let (Some(a), Some(b)) = (pts.get(i), pts.get((i + 1) % k)) { + s += a.cross(*b); + } + } + s.dot(n) * 0.5 +} + +fn inside(p: Vec3, edges: &[mt::Edge], n: Vec3) -> bool { + let u = n.any_perp(); + let w = n.cross(u); + let pts: Vec<(f64, f64)> = edges.iter().flat_map(samples).map(|q| (q.dot(u), q.dot(w))).collect(); + let (px, py) = (p.dot(u), p.dot(w)); + let k = pts.len(); + let mut c = false; + for i in 0..k { + let (Some(&(ax, ay)), Some(&(bx, by))) = (pts.get(i), pts.get((i + 1) % k)) else { continue }; + if (ay > py) != (by > py) && px < ax + (py - ay) / (by - ay) * (bx - ax) { + c = !c; + } + } + c +} + +/// Merge each group of edge-connected coplanar faces into one face. Returns None when nothing +/// changed or the result would not be a valid solid. +fn merge_coplanar(shell: &mt::Shell, tol: f64) -> Option { + let faces: Vec = shell.face_iter().cloned().collect(); + let planes: Vec> = faces.iter().map(plane_of).collect(); + let mut parent: Vec = (0..faces.len()).collect(); + let mut by_edge: HashMap> = HashMap::new(); + for (i, f) in faces.iter().enumerate() { + for e in f.edge_iter() { + by_edge.entry(e.id()).or_default().push(i); + } + } + let mut merged_any = false; + for fs in by_edge.values() { + if let [a, b] = fs[..] + && a != b + && let (Some(Some((na, da))), Some(Some((nb, db)))) = (planes.get(a), planes.get(b)) + && na.dot(*nb) > 1.0 - 1e-9 + && (da - db).abs() < tol + { + let (ra, rb) = (find(&mut parent, a), find(&mut parent, b)); + if ra != rb { + if let Some(x) = parent.get_mut(ra) { + *x = rb; + } + merged_any = true; + } + } + } + if !merged_any { + return None; + } + let mut groups: HashMap> = HashMap::new(); + for i in 0..faces.len() { + let r = find(&mut parent, i); + groups.entry(r).or_default().push(i); + } + let mut out: Vec = Vec::new(); + let mut order: Vec = groups.keys().copied().collect(); + order.sort_unstable(); + for g in order { + let members = groups.get(&g)?; + if members.len() == 1 { + out.push(faces.get(*members.first()?)?.clone()); + continue; + } + let (n, _) = (*planes.get(*members.first()?)?)?; + // Oriented boundary edges of the group; edges used twice inside it are dropped. + let mut count: HashMap = HashMap::new(); + for m in members { + for e in faces.get(*m)?.boundary_iters().into_iter().flatten() { + *count.entry(e.id()).or_insert(0) += 1; + } + } + let mut rest: Vec = Vec::new(); + for m in members { + for e in faces.get(*m)?.boundary_iters().into_iter().flatten() { + if count.get(&e.id()).copied().unwrap_or(0) == 1 { + rest.push(e); + } + } + } + // Chain into loops. + let mut loops: Vec> = Vec::new(); + while let Some(first) = rest.pop() { + let start = first.front().clone(); + let mut cur = first.back().clone(); + let mut lp = vec![first]; + let mut guard = 0; + while cur != start && guard < 100_000 { + guard += 1; + let k = rest.iter().position(|e| e.front() == &cur)?; + let e = rest.remove(k); + cur = e.back().clone(); + lp.push(e); + } + loops.push(lp); + } + let (outer, holes): (Vec>, Vec>) = loops.into_iter().partition(|l| loop_area(l, n) > 0.0); + let surface = faces.get(*members.first()?)?.oriented_surface(); + for o in &outer { + let mut wires: Vec = vec![o.clone().into()]; + for h in &holes { + let probe = h.first().map(|e| from_p3(e.front().point()))?; + // The innermost outer loop around the hole owns it. + let owner = outer.iter().filter(|x| inside(probe, x, n)).min_by(|a, b| loop_area(a, n).total_cmp(&loop_area(b, n))); + if owner.is_some_and(|x| std::ptr::eq(x, o)) { + wires.push(h.clone().into()); + } + } + out.push(mt::Face::try_new(wires, surface.clone()).ok()?); + } + } + Some(out.into()) +} + +/// Join pairs of collinear straight edges at vertices only they use. +fn merge_collinear(shell: &mt::Shell) -> Option { + let mut faces: Vec = shell.face_iter().cloned().collect(); + let mut changed = false; + for _round in 0..10_000 { + let mut at: HashMap> = HashMap::new(); + let mut seen = std::collections::HashSet::new(); + for f in &faces { + for e in f.edge_iter() { + if seen.insert(e.id()) { + at.entry(e.front().id()).or_default().push(e.absolute_clone()); + at.entry(e.back().id()).or_default().push(e.absolute_clone()); + } + } + } + let pair = at.iter().find_map(|(v, es)| { + let [a, b] = &es[..] else { return None }; + if !matches!(a.curve(), mt::Curve::Line(_)) || !matches!(b.curve(), mt::Curve::Line(_)) { + return None; + } + let other = |e: &mt::Edge| if e.front().id() == *v { e.back().clone() } else { e.front().clone() }; + let mid = if a.front().id() == *v { a.front().clone() } else { a.back().clone() }; + let (pa, pm, pb) = (from_p3(other(a).point()), from_p3(mid.point()), from_p3(other(b).point())); + let (u, w) = ((pm - pa).normalized()?, (pb - pm).normalized()?); + (u.dot(w) > 1.0 - 1e-10).then(|| (a.clone(), b.clone(), other(a), other(b))) + }); + let Some((a, b, va, vb)) = pair else { break }; + let joined = builder::line(&va, &vb); + for f in &mut faces { + if !f.edge_iter().any(|e| e.id() == a.id()) { + continue; + } + let mut wires = Vec::new(); + for w in f.absolute_boundaries() { + let mut es: Vec = Vec::new(); + for e in w.edge_iter() { + if e.id() == b.id() { + continue; + } + if e.id() == a.id() { + // a+b (in either order) become the joined edge, in this wire's direction. + es.push(if e.front() == &va { joined.clone() } else { joined.inverse() }); + continue; + } + es.push(e.clone()); + } + wires.push(mt::Wire::from(es)); + } + *f = absolute_face(f, wires)?; + } + changed = true; + } + changed.then(|| faces.into()) +} + +/// The exact curve an approximate edge (a boolean's intersection curve) really is: a line or a +/// circular arc, when every sample lies on it. +fn exact_curve(e: &mt::Edge, tol: f64) -> Option { + // Only the approximate curves booleans make; exact ones stay as they are. + if !matches!(e.curve(), mt::Curve::IntersectionCurve(_)) { + return None; + } + fit_exact(e, tol) +} + +/// A line or circular arc through the edge's samples, if they all lie on one. +fn fit_exact(e: &mt::Edge, tol: f64) -> Option { + use mt::{BoundedCurve, ParametricCurve}; + let c = e.curve(); + let (t0, t1) = c.range_tuple(); + let pts: Vec = (0..=16).map(|k| from_p3(c.subs(t0 + (t1 - t0) * k as f64 / 16.0))).collect(); + let (a, b) = (from_p3(e.absolute_front().point()), from_p3(e.absolute_back().point())); + if e.absolute_front() == e.absolute_back() || a.dist(b) < tol { + return None; + } + if pts.iter().all(|p| p.dist_to_segment(a, b) < tol) { + return Some(builder::line(e.absolute_front(), e.absolute_back())); + } + // A circle through the ends and the middle sample. + let m = *pts.get(8)?; + let n = (m - a).cross(b - a).normalized()?; + let (ab, am) = (b - a, m - a); + // Centre in the plane: solve |c-a| = |c-b| = |c-m|. + let d = 2.0 * ab.cross(am).len2(); + if d < 1e-300 { + return None; + } + // Circumcentre: a + (|u|² v×w + |v|² w×u) / (2|w|²), u = b−a, v = m−a, w = u×v. + let w = ab.cross(am); + let centre = a + (am.cross(w) * ab.len2() + w.cross(ab) * am.len2()) * (1.0 / d); + let r = centre.dist(a); + if pts.iter().all(|p| (p.dist(centre) - r).abs() < tol && (*p - a).dot(n).abs() < tol) { + return Some(builder::circle_arc(e.absolute_front(), e.absolute_back(), mt::Point3::new(m.x, m.y, m.z))); + } + None +} + +/// Replace approximate edges that are really lines or arcs with exact ones. +fn exact_edges(shell: &mt::Shell, tol: f64) -> Option { + let mut subst: HashMap = HashMap::new(); + for e in shell.edge_iter() { + if !subst.contains_key(&e.id()) + && let Some(x) = exact_curve(&e, tol) + { + subst.insert(e.id(), x); + } + } + if subst.is_empty() { + return None; + } + let mut faces = Vec::new(); + for f in shell.face_iter() { + let mut wires = Vec::new(); + for w in f.absolute_boundaries() { + let es: Vec = w + .edge_iter() + .map(|e| match subst.get(&e.id()) { + // The new edge runs absolute front → back; match this use's direction. + Some(n) => { + if e.front() == e.absolute_front() { + n.clone() + } else { + n.inverse() + } + } + None => e.clone(), + }) + .collect(); + wires.push(mt::Wire::from(es)); + } + faces.push(absolute_face(f, wires)?); + } + Some(faces.into()) +} + +/// Heal a solid; returns the input unchanged when nothing applies or healing fails. +pub fn heal(solid: Solid, size: f64) -> Solid { + let tol = (size * 1e-7).max(1e-9); + let mut shells = Vec::new(); + let mut changed = false; + for sh in solid.boundaries() { + let mut cur = sh.clone(); + if let Some(m) = exact_edges(&cur, (size * 1e-5).max(1e-7)) { + cur = m; + changed = true; + } + if let Some(m) = merge_coplanar(&cur, tol) { + cur = m; + changed = true; + } + if let Some(m) = merge_collinear(&cur) { + cur = m; + changed = true; + } + shells.push(cur); + } + if !changed { + return solid; + } + Solid::try_new(shells).unwrap_or(solid) +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn approximate_arcs_become_exact() { + // A cubic (non-rational) spline through a quarter circle. + let (a, b) = (builder::vertex(mt::Point3::new(10.0, 0.0, 0.0)), builder::vertex(mt::Point3::new(0.0, 10.0, 0.0))); + let pts: Vec = (0..=40) + .map(|k| { + let t = std::f64::consts::FRAC_PI_2 * k as f64 / 40.0; + Vec3::new(10.0 * t.cos(), 10.0 * t.sin(), 0.0) + }) + .collect(); + let curve = mt::Curve::BSplineCurve(crate::build::interpolate_cubic(&pts).expect("spline")); + let e = mt::Edge::new(&a, &b, curve); + let x = fit_exact(&e, 1e-2).expect("an arc"); + use mt::{BoundedCurve, ParametricCurve}; + let c = x.curve(); + let (t0, t1) = c.range_tuple(); + let mid = c.subs((t0 + t1) / 2.0); + assert!(((mid.x * mid.x + mid.y * mid.y).sqrt() - 10.0).abs() < 1e-9, "{mid:?}"); + } +} diff --git a/crates/kernel/src/lib.rs b/crates/kernel/src/lib.rs index 32c0ef78..5e988d2a 100644 --- a/crates/kernel/src/lib.rs +++ b/crates/kernel/src/lib.rs @@ -15,6 +15,8 @@ mod blend; mod body; mod build; +mod heal; +mod loopblend; mod measure; mod meshbody; mod ops; diff --git a/crates/kernel/src/loopblend.rs b/crates/kernel/src/loopblend.rs new file mode 100644 index 00000000..29478a04 --- /dev/null +++ b/crates/kernel/src/loopblend.rs @@ -0,0 +1,306 @@ +//! Blending a whole boundary loop of a planar face whose neighbours are walls standing +//! perpendicular to it (planes, and cylinders with axes along the face normal): a pocket floor, +//! the top outline of an extruded plate. The loop must be smooth (lines and arcs meeting +//! tangentially). The face shrinks by the radius, the walls are trimmed by it, and each loop edge +//! gets a blend face: a cylinder (or plane, for a chamfer) along a line, a torus (or cone) along +//! an arc, all meeting along the profile arcs at the joints. The result is exact. + +use std::collections::HashMap; + +use solvecraft_geom::Vec3; +use truck_modeling::{self as mt, builder}; + +use crate::body::{Solid, from_p3, p3, v3}; +use crate::{KernelError, Result}; + +fn unsupported(msg: &str) -> KernelError { + KernelError::Failed(format!("not supported yet: {msg}")) +} + +#[derive(Clone, Copy, Debug)] +enum Seg { + Line { + a: Vec3, + b: Vec3, + }, + /// Arc from `a` to `b` about `c` turning by `angle` about `axis` (right hand, angle > 0). + Arc { + a: Vec3, + c: Vec3, + axis: Vec3, + angle: f64, + mid: Vec3, + }, +} + +impl Seg { + fn start(&self) -> Vec3 { + match *self { + Seg::Line { a, .. } | Seg::Arc { a, .. } => a, + } + } + /// Unit tangent (direction of travel) at a point of the segment. + fn tangent(&self, p: Vec3) -> Option { + match *self { + Seg::Line { a, b } => (b - a).normalized(), + Seg::Arc { c, axis, .. } => axis.cross(p - c).normalized(), + } + } +} + +/// The exact shape of an oriented loop edge (a line, or an arc found from samples). +fn seg_of(e: &mt::Edge, tol: f64) -> Option { + use mt::{BoundedCurve, ParametricCurve}; + let (a, b) = (from_p3(e.front().point()), from_p3(e.back().point())); + let c = e.oriented_curve(); + let (t0, t1) = c.range_tuple(); + let pts: Vec = (0..=16).map(|k| from_p3(c.subs(t0 + (t1 - t0) * k as f64 / 16.0))).collect(); + if pts.iter().all(|p| p.dist_to_segment(a, b) < tol) { + return Some(Seg::Line { a, b }); + } + let m = *pts.get(8)?; + let (u, v) = (b - a, m - a); + let w = u.cross(v); + let d = 2.0 * w.len2(); + if d < 1e-300 { + return None; + } + let centre = a + (v.cross(w) * u.len2() + w.cross(u) * v.len2()) * (1.0 / d); + let r = centre.dist(a); + if !pts.iter().all(|p| (p.dist(centre) - r).abs() < tol) { + return None; + } + // Turning direction: the travel at a is axis × (a − c). + let t_a = (*pts.get(1)? - a).normalized()?; + let mut axis = w.normalized()?; + if axis.cross(a - centre).dot(t_a) < 0.0 { + axis = -axis; + } + let (ra, rb) = ((a - centre).normalized()?, (b - centre).normalized()?); + let mut angle = ra.dot(rb).clamp(-1.0, 1.0).acos(); + if ra.cross(rb).dot(axis) < 0.0 { + angle = std::f64::consts::TAU - angle; + } + Some(Seg::Arc { a, c: centre, axis, angle, mid: m }) +} + +fn plane_normal(f: &mt::Face) -> Option { + match f.oriented_surface() { + mt::Surface::Plane(p) => { + let n = p.normal(); + Vec3::new(n.x, n.y, n.z).normalized() + } + _ => None, + } +} + +/// Outward normal of a face at a point on it. +fn face_normal(f: &mt::Face, p: Vec3) -> Option { + use mt::{ParametricSurface3D, SearchNearestParameter, SearchParameter}; + let s = f.oriented_surface(); + let (u, v) = s.search_parameter(p3(p), None, 100).or_else(|| s.search_nearest_parameter(p3(p), None, 100))?; + let n = s.normal(u, v); + Vec3::new(n.x, n.y, n.z).normalized() +} + +/// Blend the loop of `edges` (by id) if it is a whole smooth boundary loop of a planar face with +/// perpendicular walls. `None` when the edges are not such a loop. +pub(crate) fn loop_blend(solid: &Solid, edge_ids: &[mt::EdgeID], size: f64, r: f64, round: bool) -> Option> { + let tol = (size * 1e-6).max(1e-7); + for (si, shell) in solid.boundaries().iter().enumerate() { + let faces: Vec<&mt::Face> = shell.face_iter().collect(); + for (fi, f) in faces.iter().enumerate() { + let Some(n) = plane_normal(f) else { continue }; + for (wi, w) in f.boundaries().iter().enumerate() { + let ids: Vec = w.edge_iter().map(|e| e.id()).collect(); + if ids.len() != edge_ids.len() || !edge_ids.iter().all(|e| ids.contains(e)) { + continue; + } + return Some(build(solid, si, &faces, fi, wi, n, size, tol, r, round)); + } + } + } + None +} + +#[allow(clippy::too_many_arguments)] +fn build(solid: &Solid, si: usize, faces: &[&mt::Face], fi: usize, wi: usize, n: Vec3, size: f64, tol: f64, r: f64, round: bool) -> Result { + let f = faces.get(fi).ok_or_else(|| KernelError::Failed("face".into()))?; + let wire = f.boundaries().get(wi).cloned().ok_or_else(|| KernelError::Failed("loop".into()))?; + let edges: Vec = wire.edge_iter().cloned().collect(); + let k = edges.len(); + let segs: Vec = edges + .iter() + .map(|e| seg_of(e, (size * 1e-5).max(1e-6))) + .collect::>() + .ok_or_else(|| unsupported("loop edges must be lines or arcs"))?; + // Smooth joints. + for i in 0..k { + let (Some(prev), Some(cur)) = (segs.get((i + k - 1) % k), segs.get(i)) else { continue }; + let p = cur.start(); + let (Some(ta), Some(tb)) = (prev.tangent(p), cur.tangent(p)) else { return Err(unsupported("degenerate loop edge")) }; + if ta.dot(tb) < 1.0 - 1e-6 { + return Err(unsupported("blending a loop with sharp corners")); + } + } + // Walls: the other face of each edge, perpendicular to the face; all on one side. + let mut walls: Vec = Vec::with_capacity(k); + let mut side = 0.0; + for (e, s) in edges.iter().zip(&segs) { + let wf = faces + .iter() + .position(|g| !std::ptr::eq(*g, *f) && g.edge_iter().any(|x| x.id() == e.id())) + .ok_or_else(|| unsupported("an edge with one face"))?; + let mid = match *s { + Seg::Line { a, b } => (a + b) * 0.5, + Seg::Arc { mid, .. } => mid, + }; + let g = faces.get(wf).ok_or_else(|| KernelError::Failed("wall".into()))?; + let gn = face_normal(g, mid).ok_or_else(|| unsupported("wall normal"))?; + if gn.dot(n).abs() > 1e-6 { + return Err(unsupported("the walls must stand perpendicular to the face")); + } + // Into the face from the loop: left of travel. + let m = n.cross(s.tangent(mid).ok_or_else(|| unsupported("tangent"))?); + let sd = if gn.dot(m) > 0.0 { 1.0 } else { -1.0 }; + if side != 0.0 && sd != side { + return Err(unsupported("a loop that is partly convex and partly concave")); + } + side = sd; + walls.push(wf); + } + // Offsets: into the face by r; along the wall by r (up a concave wall, down a convex one). + let into = |s: &Seg, p: Vec3| -> Vec3 { s.tangent(p).map(|t| n.cross(t)).unwrap_or(Vec3::ZERO) }; + let lift = n * (r * side); + let pa: Vec = (0..k).map(|i| segs.get(i).map(|s| s.start() + into(s, s.start()) * r).unwrap_or_default()).collect(); + let pb: Vec = (0..k).map(|i| segs.get(i).map(|s| s.start() + lift).unwrap_or_default()).collect(); + // The offset arcs must not collapse. + for s in &segs { + if let Seg::Arc { a, c, .. } = *s { + let toward = into(s, a).dot(c - a) > 0.0; + if toward && a.dist(c) <= r + tol { + return Err(unsupported("the blend is larger than a corner radius")); + } + } + } + let va: Vec = pa.iter().map(|p| builder::vertex(p3(*p))).collect(); + let vb: Vec = pb.iter().map(|p| builder::vertex(p3(*p))).collect(); + let centre = |i: usize| pa.get(i).copied().unwrap_or_default() + lift; + let profile_mid = |i: usize| { + let (Some(s), c) = (segs.get(i), centre(i)) else { return Vec3::ZERO }; + c + (s.start() - c).normalized().unwrap_or(n) * r + }; + // Profiles at the joints (face side → wall side). + let profiles: Vec = (0..k) + .map(|i| { + let (Some(a), Some(b)) = (va.get(i), vb.get(i)) else { + return builder::line(&builder::vertex(p3(Vec3::ZERO)), &builder::vertex(p3(Vec3::X))); + }; + if round { builder::circle_arc(a, b, p3(profile_mid(i))) } else { builder::line(a, b) } + }) + .collect(); + let rail = |vs: &[mt::Vertex], i: usize, offset: &dyn Fn(Vec3, &Seg) -> Vec3| -> Option { + let s = segs.get(i)?; + let (a, b) = (vs.get(i)?, vs.get((i + 1) % k)?); + Some(match *s { + Seg::Line { .. } => builder::line(a, b), + Seg::Arc { mid, .. } => builder::circle_arc(a, b, p3(offset(mid, s))), + }) + }; + let face_off = |p: Vec3, s: &Seg| p + into(s, p) * r; + let wall_off = |p: Vec3, _: &Seg| p + lift; + let face_rails: Vec = + (0..k).map(|i| rail(&va, i, &face_off)).collect::>().ok_or_else(|| KernelError::Failed("rail".into()))?; + let wall_rails: Vec = + (0..k).map(|i| rail(&vb, i, &wall_off)).collect::>().ok_or_else(|| KernelError::Failed("rail".into()))?; + // Seams between neighbouring walls at the joints get shorter. + let loop_vertex: Vec = edges.iter().map(|e| e.front().clone()).collect(); + let mut seam_subst: HashMap = HashMap::new(); + let loop_ids: Vec = edges.iter().map(|e| e.id()).collect(); + for (i, v) in loop_vertex.iter().enumerate() { + for g in faces { + for e in g.edge_iter() { + if loop_ids.contains(&e.id()) || seam_subst.contains_key(&e.id()) || !(e.front() == v || e.back() == v) { + continue; + } + let (fa, fb) = (e.absolute_front().clone(), e.absolute_back().clone()); + let other = if &fa == v { &fb } else { &fa }; + if !matches!(e.curve(), mt::Curve::Line(_)) { + return Err(unsupported("the wall seams must be straight")); + } + let along = (from_p3(other.point()) - from_p3(v.point())).dot(n) * side; + if along <= r + tol { + return Err(unsupported("the blend is taller than the walls")); + } + let nb = vb.get(i).ok_or_else(|| KernelError::Failed("joint".into()))?; + let ne = if &fa == v { builder::line(nb, &fb) } else { builder::line(&fa, nb) }; + seam_subst.insert(e.id(), ne); + } + } + } + // Rebuild faces. + let idx_of = |id: mt::EdgeID| loop_ids.iter().position(|x| *x == id); + let rebuild = |g: &mt::Face, rails: &[mt::Edge]| -> Result { + let mut wires = Vec::new(); + for w in g.absolute_boundaries() { + let mut es = Vec::new(); + for e in w.edge_iter() { + let ne = if let Some(i) = idx_of(e.id()) { + let rl = rails.get(i).ok_or_else(|| KernelError::Failed("rail".into()))?; + // Rails run in loop order (joint i → i+1). + if e.front() == loop_vertex.get(i).ok_or_else(|| KernelError::Failed("joint".into()))? { rl.clone() } else { rl.inverse() } + } else if let Some(s) = seam_subst.get(&e.id()) { + if e.front() == e.absolute_front() { s.clone() } else { s.inverse() } + } else { + e.clone() + }; + es.push(ne); + } + wires.push(mt::Wire::from(es)); + } + crate::heal::absolute_face(g, wires).ok_or_else(|| KernelError::Failed("loop blend: a face could not be rebuilt".into())) + }; + let mut out: Vec = Vec::new(); + for (gi, g) in faces.iter().enumerate() { + if gi == fi { + out.push(rebuild(g, &face_rails)?); + } else if walls.contains(&gi) { + out.push(rebuild(g, &wall_rails)?); + } else if g.edge_iter().any(|e| seam_subst.contains_key(&e.id())) { + out.push(rebuild(g, &[])?); + } else { + out.push((*g).clone()); + } + } + // Blend faces. + use mt::{ParametricSurface3D, SearchNearestParameter}; + for i in 0..k { + let (Some(s), Some(pr), Some(pr_next), Some(fr), Some(wr)) = + (segs.get(i), profiles.get(i), profiles.get((i + 1) % k), face_rails.get(i), wall_rails.get(i)) + else { + continue; + }; + let swept: mt::Shell = match *s { + Seg::Line { a, b } => vec![builder::tsweep(pr, v3(b - a))].into(), + Seg::Arc { c, axis, angle, .. } => builder::rsweep(pr, p3(c), v3(axis), mt::Rad(angle)), + }; + let mut surf = swept.face_iter().next().map(|x| x.oriented_surface()).ok_or_else(|| KernelError::Failed("blend surface".into()))?; + // Outward: toward the profile centre on a concave loop, away from it on a convex one. + let (c, mid) = (centre(i), profile_mid(i)); + let outward = (c - mid) * side; + if let Some((u, v)) = surf.search_nearest_parameter(p3(mid), None, 100) { + let nn = surf.normal(u, v); + if Vec3::new(nn.x, nn.y, nn.z).dot(outward) < 0.0 { + surf = mt::Invertible::inverse(&surf); + } + } + let w: mt::Wire = vec![fr.inverse(), pr.clone(), wr.clone(), pr_next.inverse()].into(); + let face = mt::Face::try_new(vec![w.clone()], surf.clone()).map_err(|e| KernelError::Failed(format!("blend face: {e}")))?; + out.push(face); + } + let mut shells = solid.boundaries().clone(); + if let Some(sh) = shells.get_mut(si) { + *sh = out.into(); + } + Solid::try_new(shells).map_err(|e| KernelError::Failed(format!("loop blend produced an invalid solid: {e}"))) +} diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index 8ea98f69..5cd6501f 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -121,7 +121,7 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { empty_votes += 1; last = "empty result".into(); } - Ok(s) => match Body::new(s) { + Ok(s) => match Body::new(crate::heal::heal(s, size)) { Ok(body) => { let v = volume(&body); if v > 0.0 && plausible(v, &body) { diff --git a/crates/kernel/src/polyhedron.rs b/crates/kernel/src/polyhedron.rs index 9cdeeeb4..90f3b160 100644 --- a/crates/kernel/src/polyhedron.rs +++ b/crates/kernel/src/polyhedron.rs @@ -8,7 +8,7 @@ use std::collections::HashMap; use solvecraft_geom::{Plane, Vec3}; use truck_modeling::{self as mt, builder}; -use crate::body::{Body, Solid, from_p3, p3}; +use crate::body::{Body, Solid, from_p3, p3, v3}; use crate::{KernelError, Result, guard}; /// A half-space `n · x <= d` (n is the outward unit normal). @@ -22,8 +22,9 @@ fn det3(a: Vec3, b: Vec3, c: Vec3) -> f64 { a.dot(b.cross(c)) } -/// The convex polyhedron bounded by the half-spaces. -pub fn convex_polyhedron(hs: &[HalfSpace]) -> Result { +/// Vertices, and per half-space with a face its vertex loop (counter-clockwise about the +/// outward normal), of the convex polyhedron bounded by the half-spaces. +fn combinatorics(hs: &[HalfSpace]) -> Result<(Vec, Vec<(usize, Vec)>)> { if hs.len() < 4 || hs.len() > 200 { return Err(KernelError::Invalid("a polyhedron needs 4…200 planes".into())); } @@ -53,8 +54,8 @@ pub fn convex_polyhedron(hs: &[HalfSpace]) -> Result { return Err(KernelError::Invalid("the planes don't enclose a solid".into())); } // Faces: the vertices on each plane, ordered counter-clockwise about its outward normal. - let mut faces_idx: Vec> = Vec::new(); - for h in hs { + let mut faces_idx: Vec<(usize, Vec)> = Vec::new(); + for (hi, h) in hs.iter().enumerate() { let on: Vec = (0..verts.len()).filter(|i| verts.get(*i).is_some_and(|v| (h.n.dot(*v) - h.d).abs() < scale * 1e-7)).collect(); if on.len() < 3 { continue; @@ -67,13 +68,19 @@ pub fn convex_polyhedron(hs: &[HalfSpace]) -> Result { let ang = |i: usize| verts.get(i).map(|v| (*v - c).dot(w).atan2((*v - c).dot(u))).unwrap_or(0.0); ang(*a).total_cmp(&ang(*b)) }); - faces_idx.push(ordered); + faces_idx.push((hi, ordered)); } + Ok((verts, faces_idx)) +} + +/// The convex polyhedron bounded by the half-spaces. +pub fn convex_polyhedron(hs: &[HalfSpace]) -> Result { + let (verts, faces_idx) = combinatorics(hs)?; guard("polyhedron", || { let tv: Vec = verts.iter().map(|v| builder::vertex(p3(*v))).collect(); let mut edges: HashMap<(usize, usize), mt::Edge> = HashMap::new(); let mut faces: Vec = Vec::new(); - for f in &faces_idx { + for (_, f) in &faces_idx { let n = f.len(); let mut wire: Vec = Vec::with_capacity(n); for k in 0..n { @@ -176,3 +183,122 @@ pub fn draft(b: &Body, faces: &[Vec3], neutral: &Plane, pull: Vec3, angle: f64) } convex_polyhedron(&out) } + +/// Round every edge of a convex body with planar faces by `r`: faces shrink, each edge becomes a +/// cylinder and each corner a sphere patch (the body is the inner polyhedron grown by a ball). +pub fn round_all_edges(b: &Body, r: f64) -> Result { + let hs: Vec = planar_convex(b)?.into_iter().map(|(h, _)| h).collect(); + let inner: Vec = hs.iter().map(|h| HalfSpace { n: h.n, d: h.d - r }).collect(); + let (q, faces) = combinatorics(&inner)?; + if faces.len() != hs.len() { + return Err(KernelError::Failed("not supported yet: the rounding is larger than a face".into())); + } + let normal = |hi: usize| hs.get(hi).map(|h| h.n).unwrap_or(Vec3::Z); + // Edges of the inner polyhedron: (vertex a, vertex b) → (face left, face right). + let mut edge_faces: HashMap<(usize, usize), Vec> = HashMap::new(); + for (hi, lp) in &faces { + let n = lp.len(); + for k in 0..n { + let (Some(&a), Some(&bb)) = (lp.get(k), lp.get((k + 1) % n)) else { continue }; + edge_faces.entry((a.min(bb), a.max(bb))).or_default().push(*hi); + } + } + if edge_faces.values().any(|f| f.len() != 2) { + return Err(KernelError::Failed("not supported yet: rounding this shape".into())); + } + guard("round edges", || { + // Vertex of face `hi` at inner vertex `qi`. + let mut verts: HashMap<(usize, usize), mt::Vertex> = HashMap::new(); + let mut vert = |hi: usize, qi: usize| -> mt::Vertex { + verts.entry((hi, qi)).or_insert_with(|| builder::vertex(p3(q.get(qi).copied().unwrap_or_default() + normal(hi) * r))).clone() + }; + let mut lines: HashMap<(usize, usize, usize), mt::Edge> = HashMap::new(); + let mut arcs: HashMap<(usize, usize, usize), mt::Edge> = HashMap::new(); + // Straight edge of face hi between inner vertices a and b (stored low → high). + let mut line = |hi: usize, a: usize, bb: usize, vert: &mut dyn FnMut(usize, usize) -> mt::Vertex| -> mt::Edge { + let (lo, hi_v) = (a.min(bb), a.max(bb)); + let e = lines.entry((hi, lo, hi_v)).or_insert_with(|| builder::line(&vert(hi, lo), &vert(hi, hi_v))).clone(); + if a < bb { e } else { e.inverse() } + }; + // Arc at inner vertex qi from face f to face g (stored with f < g). + let mut arc = |qi: usize, f: usize, g: usize, vert: &mut dyn FnMut(usize, usize) -> mt::Vertex| -> mt::Edge { + let (lo, hi_f) = (f.min(g), f.max(g)); + let c = q.get(qi).copied().unwrap_or_default(); + let e = arcs + .entry((qi, lo, hi_f)) + .or_insert_with(|| { + let mid = c + (normal(lo) + normal(hi_f)).normalized().unwrap_or(Vec3::Z) * r; + builder::circle_arc(&vert(lo, qi), &vert(hi_f, qi), p3(mid)) + }) + .clone(); + if f < g { e } else { e.inverse() } + }; + let mut out: Vec = Vec::new(); + // Planar faces. + for (hi, lp) in &faces { + let n = lp.len(); + let mut w: Vec = Vec::new(); + for k in 0..n { + let (Some(&a), Some(&bb)) = (lp.get(k), lp.get((k + 1) % n)) else { continue }; + w.push(line(*hi, a, bb, &mut vert)); + } + out.push(builder::try_attach_plane(&[w.into()]).map_err(|e| KernelError::Failed(format!("rounded face: {e}")))?); + } + use mt::{ParametricSurface3D, SearchNearestParameter}; + let orient = |mut surf: mt::Surface, at: Vec3, outward: Vec3| -> mt::Surface { + if let Some((u, v)) = surf.search_nearest_parameter(p3(at), None, 100) { + let nn = surf.normal(u, v); + if Vec3::new(nn.x, nn.y, nn.z).dot(outward) < 0.0 { + surf = mt::Invertible::inverse(&surf); + } + } + surf + }; + // Edge cylinders: face f's loop runs a → b, so the cylinder (to f's right) runs b → a + // along f, then across to g. + for ((a, bb), fs) in &edge_faces { + let (Some(&f0), Some(&f1)) = (fs.first(), fs.get(1)) else { continue }; + // Which face runs a → b? + let runs = |hi: usize| { + faces.iter().find(|(h, _)| *h == hi).is_some_and(|(_, lp)| { + let n = lp.len(); + (0..n).any(|k| lp.get(k) == Some(a) && lp.get((k + 1) % n) == Some(bb)) + }) + }; + let (f, g) = if runs(f0) { (f0, f1) } else { (f1, f0) }; + let w: Vec = + vec![line(f, *bb, *a, &mut vert), arc(*a, f, g, &mut vert), line(g, *a, *bb, &mut vert), arc(*bb, g, f, &mut vert)]; + let (qa, qb) = (q.get(*a).copied().unwrap_or_default(), q.get(*bb).copied().unwrap_or_default()); + let profile = arc(*a, f, g, &mut vert); + let swept = builder::tsweep(&profile, v3(qb - qa)); + let out_dir = (normal(f) + normal(g)).normalized().unwrap_or(Vec3::Z); + let mid = (qa + qb) * 0.5 + out_dir * r; + let surf = orient(swept.oriented_surface(), mid, out_dir); + out.push(mt::Face::try_new(vec![w.into()], surf).map_err(|e| KernelError::Failed(format!("rounded edge: {e}")))?); + } + // Corner sphere patches: the faces around each inner vertex, counter-clockwise. + for (qi, c) in q.iter().enumerate() { + let mut around: Vec = faces.iter().filter(|(_, lp)| lp.contains(&qi)).map(|(h, _)| *h).collect(); + let axis = around.iter().fold(Vec3::ZERO, |acc, h| acc + normal(*h)).normalized().unwrap_or(Vec3::Z); + let u = axis.any_perp(); + let w2 = axis.cross(u); + around.sort_by(|x, y| { + let ang = |h: usize| normal(h).dot(w2).atan2(normal(h).dot(u)); + ang(*x).total_cmp(&ang(*y)) + }); + let k = around.len(); + let w: Vec = (0..k).filter_map(|i| Some(arc(qi, *around.get(i)?, *around.get((i + 1) % k)?, &mut vert))).collect(); + // A sphere surface whose poles and seam are away from this patch. + let side = axis.any_perp(); + let (np, sp, back) = (*c + side * r, *c - side * r, *c - axis * r); + let meridian = builder::circle_arc(&builder::vertex(p3(np)), &builder::vertex(p3(sp)), p3(back)); + let surf = + mt::Surface::RevolutedCurve(mt::Processor::new(mt::RevolutedCurve::by_revolution(meridian.oriented_curve(), p3(*c), v3(side)))); + let surf = orient(surf, *c + axis * r, axis); + out.push(mt::Face::try_new(vec![w.into()], surf).map_err(|e| KernelError::Failed(format!("rounded corner: {e}")))?); + } + let shell: mt::Shell = out.into(); + let solid = Solid::try_new(vec![shell]).map_err(|e| KernelError::Failed(format!("rounded body: {e}")))?; + Body::new(solid) + }) +} diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index 303727bd..503f820f 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -630,4 +630,96 @@ fn mesh_bodies_move_measure_and_refuse_solid_ops() { assert!(crate::mesh_body(&pos, &[[0, 1, 99]]).is_err()); assert!(crate::mesh_body(&[Vec3::new(f64::NAN, 0.0, 0.0)], &[[0, 0, 0]]).is_err()); assert!(crate::mesh_body(&pos, &[[0, 0, 1]]).is_err()); + +fn concave_fillet() { + // An L made of a base and a wall (one solid), fillet the inside corner. + let base = box_solid(Vec3::ZERO, Vec3::new(80.0, 50.0, 8.0)).unwrap(); + let wall = box_solid(Vec3::new(0.0, 42.0, 0.0), Vec3::new(80.0, 50.0, 50.0)).unwrap(); + let l = boolean(&base, &wall, BoolOp::Union).unwrap().unwrap(); + let v0 = measure(&l).unwrap().volume; + let f = fillet(&l, &[Vec3::new(40.0, 42.0, 8.0)], 5.0).unwrap(); + let m = measure(&f).unwrap(); + assert!(rel(m.volume, v0 + (25.0 - 25.0 * PI / 4.0) * 80.0) < 1e-4, "{} vs {}", m.volume, v0 + (25.0 - 25.0 * PI / 4.0) * 80.0); + assert_eq!(m.merged.face_types.get("cylinder"), Some(&1)); +} + +#[test] +#[ignore] +fn debug_pocket_corners() { + let b = box_solid(Vec3::ZERO, Vec3::new(80.0, 60.0, 20.0)).unwrap(); + let p = box_solid(Vec3::new(15.0, 15.0, 12.0), Vec3::new(65.0, 45.0, 25.0)).unwrap(); + let mut cur = boolean(&b, &p, BoolOp::Cut).unwrap().unwrap(); + for e in [Vec3::new(65.0, 15.0, 16.0), Vec3::new(15.0, 15.0, 16.0), Vec3::new(65.0, 45.0, 16.0), Vec3::new(15.0, 45.0, 16.0)] { + match fillet(&cur, &[e], 5.0) { + Ok(n) => { + println!("ok {e:?} vol {}", measure(&n).unwrap().volume); + cur = n; + } + Err(err) => { + println!("fail {e:?}: {err}"); + break; + } + } + } +} + +#[test] +fn pocket_floor_loop_fillet() { + let b = box_solid(Vec3::ZERO, Vec3::new(80.0, 60.0, 20.0)).unwrap(); + let p = box_solid(Vec3::new(15.0, 15.0, 12.0), Vec3::new(65.0, 45.0, 25.0)).unwrap(); + let mut cur = boolean(&b, &p, BoolOp::Cut).unwrap().unwrap(); + cur = fillet(&cur, &[Vec3::new(65.0, 15.0, 16.0), Vec3::new(15.0, 15.0, 16.0), Vec3::new(65.0, 45.0, 16.0), Vec3::new(15.0, 45.0, 16.0)], 5.0) + .unwrap(); + let v1 = measure(&cur).unwrap().volume; + let floor: Vec = [(40.0, 15.0), (65.0, 30.0), (40.0, 45.0), (15.0, 30.0)].iter().map(|(x, y)| Vec3::new(*x, *y, 12.0)).collect(); + let mut pts = floor.clone(); + let d = 5.0 - 5.0 / 2f64.sqrt(); + pts.extend([(15.0 + d, 15.0 + d), (65.0 - d, 15.0 + d), (65.0 - d, 45.0 - d), (15.0 + d, 45.0 - d)].iter().map(|(x, y)| Vec3::new(*x, *y, 12.0))); + let f = fillet(&cur, &pts, 2.0).unwrap(); + let m = measure(&f).unwrap(); + // Added: the fillet profile (r²(1 − π/4)) swept along the floor outline (Pappus for the + // corners: the profile's centroid sits 2 − 4·2/(3(4 − π))·… inward; check loosely). + let straight = 2.0 * (40.0 + 20.0); + let added_lines = (4.0 - PI) * straight; + assert!(m.volume > v1 + added_lines && m.volume < v1 + added_lines + 4.0 * 2.0 * PI * 5.0, "{} {}", m.volume, v1); + println!("pocket: {:?} vol {}", m.merged, m.volume); +} + +#[test] +fn round_every_edge_of_a_box() { + let b = box_solid(Vec3::ZERO, Vec3::new(40.0, 30.0, 20.0)).unwrap(); + let mut pts = Vec::new(); + for e in b.edges(0.1).unwrap() { + pts.push(e.mid); + } + let f = fillet(&b, &pts, 3.0).unwrap(); + let m = measure(&f).unwrap(); + let r = 3.0f64; + let (a, bb, c) = (40.0 - 2.0 * r, 30.0 - 2.0 * r, 20.0 - 2.0 * r); + // Inner box grown by a ball (Steiner): V = abc + 2r(ab+bc+ca) + πr²(a+b+c) + 4/3πr³. + let want = a * bb * c + 2.0 * r * (a * bb + bb * c + c * a) + PI * r * r * (a + bb + c) + 4.0 / 3.0 * PI * r.powi(3); + assert!(rel(m.volume, want) < 1e-3, "{} vs {want}", m.volume); + assert_eq!((m.merged.faces, m.merged.edges, m.merged.vertices), (26, 48, 24), "{:?}", m.merged); +} + +#[test] +#[ignore] +fn debug_heal_sphere_cut() { + let b = box_solid(Vec3::new(-25.0, -25.0, -20.0), Vec3::new(25.0, 25.0, 0.0)).unwrap(); + let s = sphere(Vec3::ZERO, 15.0).unwrap(); + let r = guard("t", || { + let mut sb = s.deep_copy(); + sb.not(); + Ok(truck_shapeops::and(&b.deep_copy(), &sb, 0.05)) + }) + .unwrap() + .unwrap(); + let raw = Body::new(r.clone()).unwrap(); + let healed = Body::new(crate::heal::heal(r, 60.0)).unwrap(); + for (name, x) in [("raw", raw), ("healed", healed)] { + println!("{name}: vol {} faces {}", measure(&x).unwrap().volume, x.face_count()); + for f in x.faces(0.05).unwrap() { + println!(" face {} area {:.2} c {:?}", f.index, f.area, f.centroid); + } + } }