diff --git a/crates/doc/src/eval.rs b/crates/doc/src/eval.rs index 6e14cac8..e903aaf4 100644 --- a/crates/doc/src/eval.rs +++ b/crates/doc/src/eval.rs @@ -373,6 +373,73 @@ fn region_samples(r: &Region2) -> Vec { out } +/// Grow profile edges that lie on a face of a target body away from the profile, by a small +/// amount, when the space they grow into is already material (join) or empty (cut). `probe` +/// maps a sketch point to the world points where the tool would be (along the extrude, around +/// the revolve). Only loops of straight edges are adjusted. +fn grow_profile_for_coplanar(st: &ModelState, r: &Region2, op: Operation, targets: &[String], probe: &dyn Fn(Vec2) -> Vec) -> Region2 { + if !matches!(op, Operation::Cut | Operation::Join) { + return r.clone(); + } + let bodies: Vec<&ModelBody> = st.bodies.iter().filter(|b| targets.is_empty() || targets.contains(&b.name)).collect(); + if bodies.is_empty() { + return r.clone(); + } + let meshes: Vec> = bodies.iter().map(|b| b.mesh()).collect(); + let size = bodies.iter().map(|b| b.body.size()).fold(0.0, f64::max); + let delta = (size * 0.005).max(1e-3); + let lp = r.outer.ccw(); + let lines: Option> = + lp.segs.iter().map(|s| if let solvecraft_geom::Seg2::Line { a, b } = *s { Some((a, b)) } else { None }).collect(); + let Some(lines) = lines else { return r.clone() }; + let inside_any = |p: Vec3| meshes.iter().any(|m| m.contains(p)); + let mut grow = vec![0.0; lines.len()]; + let mut any = false; + for (i, (a, b)) in lines.iter().enumerate() { + let Some(dir) = (*b - *a).normalized() else { continue }; + let out = Vec2::new(dir.y, -dir.x); + // Is the edge on a body face? Material on exactly one side of it. + let mids = [a.lerp(*b, 0.25), a.lerp(*b, 0.5), a.lerp(*b, 0.75)]; + let side = |k: f64| mids.iter().flat_map(|m| probe(*m + out * (delta * k))).map(inside_any).collect::>(); + let (outer, inner) = (side(0.5), side(-0.5)); + let want_outer = op == Operation::Join; + let on_face = outer.iter().all(|x| *x == want_outer) && inner.iter().any(|x| *x != want_outer); + let room = side(1.0).iter().all(|x| *x == want_outer); + if on_face && room { + if let Some(g) = grow.get_mut(i) { + *g = delta; + } + any = true; + } + } + if !any { + return r.clone(); + } + // Offset the chosen lines and re-intersect neighbours. + let n = lines.len(); + let shifted: Vec<(Vec2, Vec2)> = lines + .iter() + .zip(&grow) + .map(|((a, b), g)| { + let dir = (*b - *a).normalized().unwrap_or(Vec2::X); + let out = Vec2::new(dir.y, -dir.x); + (*a + out * *g, dir) + }) + .collect(); + let mut pts = Vec::with_capacity(n); + for i in 0..n { + let (Some(&(p0, d0)), Some(&(p1, d1))) = (shifted.get((i + n - 1) % n), shifted.get(i)) else { return r.clone() }; + let den = d0.cross(d1); + if den.abs() < 1e-12 { + pts.push(p1); + continue; + } + let t = (p1 - p0).cross(d1) / den; + pts.push(p0 + d0 * t); + } + Region2 { outer: solvecraft_geom::Loop2::polygon(&pts), holes: r.holes.clone() } +} + /// Booleans struggle with coincident faces (a cut starting on the face it cuts, a boss /// starting on the face it joins). Extending the tool past such an end changes nothing when the /// space beyond is empty (cut) or already material (join), so do exactly that. @@ -498,7 +565,10 @@ fn feature_tools(vals: &BTreeMap, f: &Feature, st: &ModelState) - let mut tools = Vec::new(); for r in ®ions { let (l2, h2) = extend_for_coplanar(st, &ss.plane, r, lo + off, hi + off, *operation, targets); - tools.extend(kernel::extrude(&ss.plane, std::slice::from_ref(r), l2, h2)?); + let n = ss.plane.normal(); + let probe = |p: Vec2| (1..4).map(|k| ss.plane.to_world(p) + n * (lo + off + (hi - lo) * k as f64 / 4.0)).collect::>(); + let r = grow_profile_for_coplanar(st, r, *operation, targets, &probe); + tools.extend(kernel::extrude(&ss.plane, std::slice::from_ref(&r), l2, h2)?); } Ok(tools) } @@ -507,6 +577,20 @@ fn feature_tools(vals: &BTreeMap, f: &Feature, st: &ModelState) - let regions = solvecraft_sketch::merge_regions(&select_profiles(&ss, profiles)?); let (o, d) = revolve_axis(&ss, axis)?; let ang = val(vals, angle, Kind::Angle)?; + let (operation, targets) = feature_op(f); + let axis_w = ss.plane.dir_to_world(d); + let origin_w = ss.plane.to_world(o); + let probe = |p: Vec2| { + let q = ss.plane.to_world(p); + (1..4) + .map(|k| { + let a = ang * k as f64 / 4.0; + let v = q - origin_w; + origin_w + v * a.cos() + axis_w.cross(v) * a.sin() + axis_w * (axis_w.dot(v) * (1.0 - a.cos())) + }) + .collect::>() + }; + let regions: Vec = regions.iter().map(|r| grow_profile_for_coplanar(st, r, operation, &targets, &probe)).collect(); Ok(kernel::revolve(&ss.plane, ®ions, o, d, ang)?) } FeatureKind::Box { corner, length, width, height, .. } => { diff --git a/crates/kernel/src/build.rs b/crates/kernel/src/build.rs index 3eed4222..fa50c0f0 100644 --- a/crates/kernel/src/build.rs +++ b/crates/kernel/src/build.rs @@ -140,7 +140,22 @@ pub fn revolve(plane: &Plane, regions: &[Region2], axis_origin: Vec2, axis_dir: .collect(), }; let r = &Region2 { outer: snap_loop(&r.outer), holes: r.holes.iter().map(snap_loop).collect() }; + // A full revolution of a profile with one edge on the axis: sweep the rest of the loop + // as a cone-like shell (no zero-area faces, which booleans can't handle). + if angle.abs() >= std::f64::consts::TAU - 1e-9 + && r.holes.is_empty() + && let Some(b) = revolve_touching_axis(plane, &r.outer, axis_origin, d)? + { + out.push(b); + continue; + } let f = face(plane, r)?; + // Partial revolves sweep toward the side the sketch normal points to (as features + // extrude toward it). + let c = plane.to_world(r.centroid()); + let radial = (c - o) - axis * (c - o).dot(axis); + let toward = axis.cross(radial).dot(plane.normal()); + let angle = if toward < 0.0 { -angle } else { angle }; let solid: Solid = guard("revolve", || Ok(builder::rsweep(&f, p3(o), v3(axis), mt::Rad(angle))))?; out.push(Body::new(solid)?); } @@ -296,3 +311,52 @@ pub fn extrude_tapered(plane: &Plane, region: &Region2, length: f64, dir_sign: f Err(KernelError::Failed(format!("tapered extrude: {last}"))) }) } + +/// Revolve 360° a loop that has exactly one line segment on the axis, without degenerate faces. +fn revolve_touching_axis(plane: &Plane, lp: &Loop2, axis_origin: Vec2, d: Vec2) -> Result> { + let on_axis = |p: Vec2| d.cross(p - axis_origin).abs() < 1e-9; + let n = lp.segs.len(); + let axis_segs: Vec = + (0..n).filter(|i| lp.segs.get(*i).is_some_and(|s| matches!(s, Seg2::Line { a, b } if on_axis(*a) && on_axis(*b)))).collect(); + let [k] = axis_segs[..] else { return Ok(None) }; + // The open chain after the axis segment, back around to it. + let mut chain: Vec = (1..n).filter_map(|j| lp.segs.get((k + j) % n).copied()).collect(); + // A single arc (a sphere) is split off-centre: a seam on the equator would lie in the very + // planes other bodies are often cut with. + if let [Seg2::Arc { center, radius, start, sweep }] = chain[..] { + let f = 0.37; + chain = vec![ + Seg2::Arc { center, radius, start, sweep: sweep * f }, + Seg2::Arc { center, radius, start: start + sweep * f, sweep: sweep * (1.0 - f) }, + ]; + } + if chain.is_empty() || chain.iter().skip(1).take(chain.len().saturating_sub(2)).any(|s| on_axis(s.start()) && on_axis(s.end())) { + return Ok(None); + } + // Built with the axis through the world origin (truck tessellates such cones reliably), then + // moved into place. + let shift = plane.to_world(axis_origin); + let local = Plane { origin: plane.origin - shift, ..*plane }; + let plane = &local; + guard("revolve", || { + let verts: Vec = chain + .iter() + .map(|s| builder::vertex(p3(plane.to_world(s.start())))) + .chain(chain.last().map(|s| builder::vertex(p3(plane.to_world(s.end()))))) + .collect(); + let mut edges: Vec = Vec::new(); + for (i, s) in chain.iter().enumerate() { + let (Some(a), Some(b)) = (verts.get(i), verts.get(i + 1)) else { continue }; + edges.push(match *s { + Seg2::Line { .. } => builder::line(a, b), + Seg2::Arc { .. } => builder::circle_arc(a, b, p3(plane.to_world(s.mid()))), + }); + } + let w: mt::Wire = edges.into(); + let axis = plane.dir_to_world(d); + let shell = builder::cone(&w, v3(axis), mt::Rad(std::f64::consts::TAU)); + let solid = Solid::try_new(vec![shell]).map_err(|e| KernelError::Failed(format!("revolve: {e}")))?; + let solid = if shift.len() > 0.0 { builder::translated(&solid, v3(shift)) } else { solid }; + Body::new(solid).map(Some) + }) +} diff --git a/crates/kernel/src/measure.rs b/crates/kernel/src/measure.rs index d7359791..58b7cd02 100644 --- a/crates/kernel/src/measure.rs +++ b/crates/kernel/src/measure.rs @@ -1,5 +1,5 @@ use serde::Serialize; -use solvecraft_geom::{Aabb3, Vec3}; +use solvecraft_geom::{Aabb3, Mesh, Vec3}; use crate::Result; use crate::body::Body; @@ -22,21 +22,46 @@ pub struct BodyMeasure { pub merged: TopoCounts, } -/// Measure a body from a fine tessellation (chord tolerance 5e-5 of the body size). Curved -/// faces are slightly under-estimated (relative error typically below 2e-4). +/// Per-face sums: area, 6 × signed volume and 24 × first moments (divergence theorem, origin 0). +fn face_sums(m: &Mesh, nf: usize) -> Vec<(f64, f64, Vec3)> { + let mut out = vec![(0.0, 0.0, Vec3::ZERO); nf]; + for (t, f) in m.triangles.iter().zip(&m.tri_face) { + let (Some([a, b, c]), Some(slot)) = (m.tri(t), out.get_mut(*f as usize)) else { continue }; + let d = a.dot(b.cross(c)); + slot.0 += (b - a).cross(c - a).len() * 0.5; + slot.1 += d; + slot.2 += (a + b + c) * d; + } + out +} + +/// Measure a body from a fine tessellation (chord tolerance 5e-5 of the body size). Faces whose +/// fine tessellation disagrees with a coarse one (a meshing failure) use the coarse result. pub fn measure(b: &Body) -> Result { let size = b.size(); - let m2 = b.tessellate((size * 5e-5).max(1e-4))?; - let a2 = m2.measure(); + let fine = b.tessellate((size * 5e-5).max(1e-4))?; + let coarse = b.tessellate((size * 1e-3).max(1e-3))?; + let nf = b.face_count(); + let (sf, sc) = (face_sums(&fine, nf), face_sums(&coarse, nf)); + let (mut area, mut v6, mut mom) = (0.0, 0.0, Vec3::ZERO); + for (f, c) in sf.iter().zip(&sc) { + let bad = (f.0 - c.0).abs() > 0.01 * c.0.max(f.0) + 1e-9; + let pick = if bad { c } else { f }; + area += pick.0; + v6 += pick.1; + mom += pick.2; + } + let volume = v6 / 6.0; + let centroid = if v6.abs() > 1e-12 { mom / (4.0 * v6) } else { fine.bounds().center() }; Ok(BodyMeasure { - volume: a2.volume, - area: a2.area, - centroid: a2.centroid, - bbox: m2.bounds(), + volume, + area, + centroid, + bbox: coarse.bounds().union(&fine.bounds()), faces: b.face_count(), edges: b.edge_count(), vertices: b.vertex_count(), shells: b.shell_count(), - merged: merged_topology(b, &m2)?, + merged: merged_topology(b, &coarse)?, }) } diff --git a/crates/kernel/src/ops.rs b/crates/kernel/src/ops.rs index 32e486bf..fa2a698c 100644 --- a/crates/kernel/src/ops.rs +++ b/crates/kernel/src/ops.rs @@ -24,6 +24,51 @@ fn volume(b: &Body) -> f64 { b.tessellate(b.size() * 2e-3).map(|m| m.measure().volume).unwrap_or(f64::NAN) } +/// Halton low-discrepancy value (base `b`) for index `i`. +fn halton(mut i: usize, b: usize) -> f64 { + let (mut f, mut r) = (1.0, 0.0); + while i > 0 { + f /= b as f64; + r += f * (i % b) as f64; + i /= b; + } + r +} + +/// Fraction of sample points where `result` disagrees with `op` applied to `a` and `b` +/// (membership by ray parity). Samples cover the overlap box densely and the union box lightly. +fn mismatch(a: &solvecraft_geom::Mesh, b: &solvecraft_geom::Mesh, result: &solvecraft_geom::Mesh, op: BoolOp) -> Option { + if a.triangles.len() + b.triangles.len() + result.triangles.len() > 60_000 { + return None; + } + let (ba, bb) = (a.bounds(), b.bounds()); + let lo = ba.min.max(bb.min); + let hi = ba.max.min(bb.max); + let all = ba.union(&bb); + let mut boxes = vec![(all, 300usize)]; + if lo.x < hi.x && lo.y < hi.y && lo.z < hi.z { + boxes.push((solvecraft_geom::Aabb3 { min: lo, max: hi }, 900)); + } + let (mut n, mut bad) = (0usize, 0usize); + for (bx, count) in boxes { + let s = bx.size(); + for i in 1..=count { + let p = bx.min + Vec3::new(s.x * halton(i, 2), s.y * halton(i, 3), s.z * halton(i, 5)); + let (ia, ib) = (a.contains(p), b.contains(p)); + let want = match op { + BoolOp::Union => ia || ib, + BoolOp::Cut => ia && !ib, + BoolOp::Intersect => ia && ib, + }; + n += 1; + if result.contains(p) != want { + bad += 1; + } + } + } + Some(bad as f64 / n.max(1) as f64) +} + /// Boolean of two bodies. The result may be empty (`Ok(None)`) for a cut that removes /// everything or an intersection of disjoint bodies. Results are checked against volume bounds /// and retried with shifted copies and other tolerances when they fail or look wrong. @@ -31,10 +76,20 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { let size = a.size().max(b.size()); let (va, vb) = (volume(a), volume(b)); let slack = 2e-3 * (va + vb) + 1e-9; - let plausible = |v: f64| match op { - BoolOp::Union => v >= va.max(vb) - slack && v <= va + vb + slack, - BoolOp::Cut => v >= va - vb - slack && v <= va + slack, - BoolOp::Intersect => v <= va.min(vb) + slack, + let tol_m = size * 5e-4; + let meshes = (a.tessellate(tol_m).ok(), b.tessellate(tol_m).ok()); + let plausible = |v: f64, r: &Body| { + let bounds = match op { + BoolOp::Union => v >= va.max(vb) - slack && v <= va + vb + slack, + BoolOp::Cut => v >= va - vb - slack && v <= va + slack, + BoolOp::Intersect => v <= va.min(vb) + slack, + }; + // Membership check against the operands (catches misclassified pieces). + let consistent = match (&meshes, r.tessellate(tol_m)) { + ((Some(ma), Some(mb)), Ok(mr)) => mismatch(ma, mb, &mr, op).is_none_or(|f| f < 0.004), + _ => true, + }; + bounds && consistent }; let mut last = String::new(); let mut empty_votes = 0; @@ -67,7 +122,7 @@ pub fn boolean(a: &Body, b: &Body, op: BoolOp) -> Result> { Ok(s) => match Body::new(s) { Ok(body) => { let v = volume(&body); - if v > 0.0 && plausible(v) { + if v > 0.0 && plausible(v, &body) { return Ok(Some(body)); } last = format!("implausible result volume {v:.4}"); diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index 2f1d3dea..4cc936c7 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -295,3 +295,38 @@ fn debug_chamfer_cut() { println!("tess {tol}: {}", r.tessellate(tol).unwrap().measure().volume); } } + +#[test] +#[ignore] +fn debug_revolved_hole_tool() { + let b = box_solid(Vec3::ZERO, Vec3::new(60.0, 40.0, 15.0)).unwrap(); + // Drill point profile in the XZ plane through x = 50: revolve about the vertical line x = 50. + let pl = Plane { origin: Vec3::new(50.0, 20.0, 0.0), x: Vec3::X, y: Vec3::Z }; + let r = 3.0; + let tip = r / 59f64.to_radians().tan(); + let prof = + Region2 { outer: Loop2::polygon(&[Vec2::new(0.0, 16.0), Vec2::new(0.0, 5.0 - tip), Vec2::new(r, 5.0), Vec2::new(r, 16.0)]), holes: vec![] }; + let tool = revolve(&pl, &[prof], Vec2::ZERO, Vec2::Y, std::f64::consts::TAU).unwrap().pop().unwrap(); + println!("tool vol {}", measure(&tool).unwrap().volume); + let res = boolean(&b, &tool, BoolOp::Cut); + println!("cut: {:?}", res.map(|o| o.map(|x| measure(&x).unwrap().volume))); + let expect = 36000.0 - std::f64::consts::PI * 9.0 * 10.0 - std::f64::consts::PI * 9.0 * tip / 3.0; + println!("expect {expect}"); +} + +#[test] +#[ignore] +fn debug_sphere() { + let s = sphere(Vec3::new(1.0, 2.0, 3.0), 10.0).unwrap(); + let m = measure(&s).unwrap(); + println!("vol {} area {} faces {} bbox {:?}", m.volume, m.area, m.faces, m.bbox); + for f in s.faces(0.05).unwrap() { + println!("face {} area {:.3} c {:?}", f.index, f.area, f.centroid); + } + let b = box_solid(Vec3::new(-25.0, -25.0, -20.0), Vec3::new(25.0, 25.0, 0.0)).unwrap(); + for (name, c) in [("origin", Vec3::ZERO), ("up", Vec3::new(0.0, 0.0, 3.0)), ("off", Vec3::new(1.3, 2.1, 0.7))] { + let sp = sphere(c, 15.0).unwrap(); + println!("{name} cut {:?}", boolean(&b, &sp, BoolOp::Cut).map(|o| o.map(|x| measure(&x).unwrap().volume))); + println!("{name} and {:?}", boolean(&b, &sp, BoolOp::Intersect).map(|o| o.map(|x| measure(&x).unwrap().volume))); + } +}