diff --git a/ROADMAP.md b/ROADMAP.md index a866bac0..9fa5b3c9 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -83,8 +83,8 @@ headline counts in-scope tabs only and lists the deferred ones separately. - No persistent naming (fillet edges are re-found by position; projected sketch geometry too). - Free-form curves cut at a tangent touch (no clean crossing) enter that profile as polylines. Sweeps along 3D sketch curves follow their polyline (smooth walls when the path is smooth). -- Shell and offset faces work on bodies of planes and cylinders (convex or not: filleted boxes, - bosses); draft on convex planar bodies and on walls (planes, cylinders → cones) between caps; loft is ruled. +- Shell and offset faces work on bodies of planes, cylinders, cones and spheres (convex or not: + filleted boxes, bosses, domes, countersinks); draft on convex planar bodies and on walls (planes, cylinders → cones) between caps; loft is ruled. ## Kernel diff --git a/crates/kernel/src/offset.rs b/crates/kernel/src/offset.rs index 1fd9f939..414d3ca9 100644 --- a/crates/kernel/src/offset.rs +++ b/crates/kernel/src/offset.rs @@ -20,6 +20,11 @@ enum Surf { /// Axis point and direction, radius; `convex` when the outward normal points away from /// the axis. Cylinder { o: Vec3, a: Vec3, r: f64, convex: bool }, + /// Centre and radius; `convex` when the outward normal points away from the centre. + Sphere { c: Vec3, r: f64, convex: bool }, + /// Apex, unit axis (into the cone, away from the apex) and half-angle; `convex` when the + /// outward normal points away from the axis. + Cone { v: Vec3, a: Vec3, half: f64, convex: bool }, } impl Surf { @@ -32,6 +37,16 @@ impl Surf { let radial = (v - a * v.dot(a)).normalized()?; Some(if convex { radial } else { -radial }) } + Surf::Sphere { c, convex, .. } => { + let radial = (p - c).normalized()?; + Some(if convex { radial } else { -radial }) + } + Surf::Cone { v, a, half, convex } => { + let w = p - v; + let q = (w - a * w.dot(a)).normalized()?; + let n = q * half.cos() - a * half.sin(); + Some(if convex { n } else { -n }) + } } } fn moved(&self, s: f64) -> Option { @@ -41,6 +56,16 @@ impl Surf { let r2 = if convex { r + s } else { r - s }; (r2 > 1e-9).then_some(Surf::Cylinder { o, a, r: r2, convex }) } + Surf::Sphere { c, r, convex } => { + let r2 = if convex { r + s } else { r - s }; + (r2 > 1e-9).then_some(Surf::Sphere { c, r: r2, convex }) + } + // Moving a cone along its normal moves its apex along the axis. + Surf::Cone { v, a, half, convex } => { + let s = if convex { s } else { -s }; + let sh = half.sin(); + (sh > 1e-9).then_some(Surf::Cone { v: v - a * (s / sh), a, half, convex }) + } } } fn dist(&self, p: Vec3) -> f64 { @@ -50,11 +75,26 @@ impl Surf { let v = p - o; (v - a * v.dot(a)).len() - r } + Surf::Sphere { c, r, .. } => p.dist(c) - r, + Surf::Cone { v, a, half, .. } => { + let w = p - v; + let h = w.dot(a); + (w - a * h).len() * half.cos() - h * half.sin() + } } } /// The affine map taking this surface onto `to` (same kind, same axis). fn map_to(&self, to: &Surf) -> Option { match (*self, *to) { + (Surf::Sphere { c, r, .. }, Surf::Sphere { r: r2, .. }) => { + let k = r2 / r; + let t = c * (1.0 - k); + Some(mt::Matrix4::from_translation(mt::Vector3::new(t.x, t.y, t.z)) * mt::Matrix4::from_scale(k)) + } + (Surf::Cone { v, .. }, Surf::Cone { v: v2, .. }) => { + let t = v2 - v; + Some(mt::Matrix4::from_translation(mt::Vector3::new(t.x, t.y, t.z))) + } (Surf::Plane { n, d }, Surf::Plane { d: d2, .. }) => { let t = n * (d2 - d); Some(mt::Matrix4::from_translation(mt::Vector3::new(t.x, t.y, t.z))) @@ -109,7 +149,15 @@ fn samples(f: &mt::Face) -> Vec<(Vec3, Vec3)> { let c = e.oriented_curve(); let (t0, t1) = c.range_tuple(); for k in 0..5 { - if let Some(x) = surf.search_nearest_parameter(c.subs(t0 + (t1 - t0) * k as f64 / 5.0), uv.last().copied(), 50) { + // The previous point's parameters as a hint, unless that lands elsewhere (a + // singular point such as a pole holds the search). + let q = c.subs(t0 + (t1 - t0) * k as f64 / 5.0); + let near = |x: (f64, f64)| from_p3(surf.subs(x.0, x.1)).dist(from_p3(q)) < 1e-6 * (1.0 + from_p3(q).len()); + let found = surf + .search_nearest_parameter(q, uv.last().copied(), 50) + .filter(|x| near(*x)) + .or_else(|| surf.search_nearest_parameter(q, None, 100).filter(|x| near(*x))); + if let Some(x) = found { uv.push(x); } } @@ -143,13 +191,14 @@ fn surf_of(f: &mt::Face, tol: f64) -> Option { } let s = samples(f); let (p0, n0) = *s.first()?; - let far = s.iter().map(|x| x.1).min_by(|a, b| a.dot(n0).total_cmp(&b.dot(n0)))?; - let a = n0.cross(far).normalized()?; + // The normal most square to the first (a half cylinder's two ends are opposite). + let far = s.iter().map(|x| x.1).min_by(|a, b| a.dot(n0).abs().total_cmp(&b.dot(n0).abs()))?; + let Some(a) = n0.cross(far).normalized() else { return sphere_or_cone(&s, tol) }; if s.iter().any(|(_, n)| n.dot(a).abs() > 1e-6) { - return None; + return sphere_or_cone(&s, tol); } // The axis: where the normal lines through two samples meet (seen along a). - let (p1, n1) = *s.iter().max_by(|x, y| (x.1 - n0).len().total_cmp(&(y.1 - n0).len()))?; + let (p1, n1) = *s.iter().min_by(|x, y| x.1.dot(n0).abs().total_cmp(&y.1.dot(n0).abs()))?; // p0 + t0 n0 = p1 + t1 n1 (projected): solve in the plane ⟂ a. let w = p1 - p0; let c = n0.dot(n1); @@ -168,6 +217,72 @@ fn surf_of(f: &mt::Face, tol: f64) -> Option { ok.then_some(cyl) } +/// A sphere (normal lines all through one point) or a cone (tangent planes all through one +/// point, normals at one angle to an axis) through the samples. +fn sphere_or_cone(s: &[(Vec3, Vec3)], tol: f64) -> Option { + // Least squares point nearest all the normal lines: Σ (I − n nᵀ) x = Σ (I − n nᵀ) p. + let mut m = [[0.0f64; 3]; 3]; + let mut r = [0.0f64; 3]; + for (p, n) in s { + let nv = [n.x, n.y, n.z]; + let pv = [p.x, p.y, p.z]; + for i in 0..3 { + for j in 0..3 { + let a = (if i == j { 1.0 } else { 0.0 }) - nv[i] * nv[j]; + m[i][j] += a; + r[i] += a * pv[j]; + } + } + } + if let Some(x) = crate::polyhedron::solve3_pub(m, r) { + let c = Vec3::new(x[0], x[1], x[2]); + let (p0, n0) = *s.first()?; + let rad = p0.dist(c); + let convex = (p0 - c).dot(n0) > 0.0; + let sp = Surf::Sphere { c, r: rad, convex }; + if rad > tol && s.iter().all(|(p, n)| sp.dist(*p).abs() < tol && sp.normal(*p).is_some_and(|q| q.dot(*n) > 1.0 - 1e-6)) { + return Some(sp); + } + } + // Cone: the apex lies on every tangent plane, n·(v − p) = 0. + let mut m = [[0.0f64; 3]; 3]; + let mut r = [0.0f64; 3]; + for (p, n) in s { + let nv = [n.x, n.y, n.z]; + let k = n.dot(*p); + for i in 0..3 { + for j in 0..3 { + m[i][j] += nv[i] * nv[j]; + } + r[i] += nv[i] * k; + } + } + let x = crate::polyhedron::solve3_pub(m, r)?; + let v = Vec3::new(x[0], x[1], x[2]); + // Axis: from the apex toward the samples' centroid (the cone opens that way). + let cen = s.iter().fold(Vec3::ZERO, |acc, (p, _)| acc + *p) * (1.0 / s.len() as f64); + // Normals make a constant angle with the axis: fit it from three spread normals. + let n0 = s.first()?.1; + let n1 = s.iter().map(|x| x.1).max_by(|a, b| (*a - n0).len().total_cmp(&(*b - n0).len()))?; + let n2 = s.iter().map(|x| x.1).max_by(|a, b| (*a - n0).cross(n1 - n0).len().total_cmp(&(*b - n0).cross(n1 - n0).len()))?; + let mut a = (n1 - n0).cross(n2 - n0).normalized()?; + if a.dot(cen - v) < 0.0 { + a = -a; + } + let half = n0.dot(a).abs().clamp(0.0, 1.0).asin(); + if !(half > 1e-3 && half < std::f64::consts::FRAC_PI_2 - 1e-3) { + return None; + } + let (p0, n0) = *s.first()?; + let q = { + let w = p0 - v; + (w - a * w.dot(a)).normalized()? + }; + let convex = n0.dot(q) > 0.0; + let cone = Surf::Cone { v, a, half, convex }; + s.iter().all(|(p, n)| cone.dist(*p).abs() < tol && cone.normal(*p).is_some_and(|m| m.dot(*n) > 1.0 - 1e-6)).then_some(cone) +} + /// Move `p` (on all of `old`) onto the moved surfaces `new`: the smallest displacement δ with /// n_i · δ = s_i for each face's normal n_i at p, then checked against the moved surfaces. fn move_point(p: Vec3, old: &[Surf], shifts: &[f64], new: &[Surf], tol: f64) -> Option { @@ -242,10 +357,11 @@ pub(crate) fn offset_body(b: &Body, shift: impl Fn(usize, Vec3) -> f64) -> Resul let mut new = Vec::with_capacity(faces.len()); let mut shifts = Vec::with_capacity(faces.len()); for (i, f) in faces.iter().enumerate() { - let sf = surf_of(f, tol * 100.0).ok_or_else(|| fail("only bodies with planar and cylindrical faces"))?; + let sf = surf_of(f, tol * 100.0).ok_or_else(|| fail("only bodies with planes, cylinders, cones and spheres"))?; let n = match sf { Surf::Plane { n, .. } => n, - Surf::Cylinder { a, .. } => a.any_perp(), + Surf::Cylinder { a, .. } | Surf::Cone { a, .. } => a.any_perp(), + Surf::Sphere { .. } => Vec3::Z, }; let s = shift(i, n); old.push(sf); @@ -356,6 +472,7 @@ pub(crate) fn draft_walls(b: &Body, chosen: &[usize], neutral: &solvecraft_geom: let wall = match sf { Surf::Plane { n, .. } => n.dot(pull).abs() < 1e-9, Surf::Cylinder { a, .. } => a.dot(pull).abs() > 1.0 - 1e-9, + Surf::Sphere { .. } | Surf::Cone { .. } => false, }; let cap = matches!(sf, Surf::Plane { n, .. } if n.dot(pull).abs() > 1.0 - 1e-9); if chosen.contains(&i) && !wall { @@ -446,6 +563,7 @@ pub(crate) fn draft_walls(b: &Body, chosen: &[usize], neutral: &solvecraft_geom: let rc = mt::RevolutedCurve::by_revolution(line, p3(o), mt::Vector3::new(a.x, a.y, a.z)); mt::Surface::RevolutedCurve(mt::Processor::new(rc)) } + Surf::Sphere { .. } | Surf::Cone { .. } => return Err(fail("draft: only planes and cylinders along the pull can be drafted")), }; if let Some(slot) = surfaces.get_mut(i) { *slot = Some(surface); diff --git a/crates/kernel/src/polyhedron.rs b/crates/kernel/src/polyhedron.rs index 1a0f0bd1..d333c994 100644 --- a/crates/kernel/src/polyhedron.rs +++ b/crates/kernel/src/polyhedron.rs @@ -362,8 +362,18 @@ fn shell_planar(b: &Body, open: &[Vec3], thickness: f64) -> Result { } let margin = thickness + size * 0.05; let cavity = offset_planar(&healed, |fi, _| if opened.contains(&fi) { margin } else { -thickness })?; - let r = crate::ops::boolean(&healed, &cavity, crate::BoolOp::Cut)?.ok_or_else(|| KernelError::Failed("the shell removed everything".into()))?; - Ok(r) + match crate::ops::boolean(&healed, &cavity, crate::BoolOp::Cut) { + Ok(Some(r)) => Ok(r), + Ok(None) => Err(KernelError::Failed("the shell removed everything".into())), + // Healing can leave faces the boolean dislikes (revolved patches): the body as it was. + Err(e) if healed.face_count() == b.face_count() => { + let cavity = offset_planar(b, |fi, _| if opened.contains(&fi) { margin } else { -thickness }).map_err(|_| e.clone())?; + crate::ops::boolean(b, &cavity, crate::BoolOp::Cut) + .map_err(|_| e)? + .ok_or_else(|| KernelError::Failed("the shell removed everything".into())) + } + Err(e) => Err(e), + } } fn point_tri(p: Vec3, a: Vec3, b: Vec3, c: Vec3) -> f64 { diff --git a/crates/kernel/src/tests.rs b/crates/kernel/src/tests.rs index 775b714e..98924243 100644 --- a/crates/kernel/src/tests.rs +++ b/crates/kernel/src/tests.rs @@ -1211,3 +1211,30 @@ fn fillet_three_edges_at_a_corner() { let m = measure(&f).unwrap().merged; assert_eq!(m.face_types.get("sphere"), Some(&1), "{m:?}"); } + +#[test] +fn shell_with_spherical_and_conical_faces() { + // A cylinder capped by a hemisphere (a dome), revolved, opened at the bottom. + use solvecraft_geom::Seg2; + let xz = Plane::new(Vec3::ZERO, Vec3::X, Vec3::new(0.0, 0.0, -1.0)).unwrap(); + let segs = vec![ + Seg2::Line { a: Vec2::new(0.0, 0.0), b: Vec2::new(10.0, 0.0) }, + Seg2::Line { a: Vec2::new(10.0, 0.0), b: Vec2::new(10.0, -20.0) }, + Seg2::Arc { center: Vec2::new(0.0, -20.0), radius: 10.0, start: 0.0, sweep: -PI / 2.0 }, + Seg2::Line { a: Vec2::new(0.0, -30.0), b: Vec2::new(0.0, 0.0) }, + ]; + let b = revolve(&xz, &[Region2 { outer: Loop2 { segs }, holes: vec![] }], Vec2::ZERO, Vec2::new(0.0, -1.0), 2.0 * PI).unwrap().pop().unwrap(); + let s = shell(&b, &[Vec3::new(0.0, 0.0, 0.0)], 1.0).unwrap(); + let outer = PI * 100.0 * 20.0 + 2.0 / 3.0 * PI * 1000.0; + let inner = PI * 81.0 * 20.0 + 2.0 / 3.0 * PI * 729.0; + let v = measure(&s).unwrap().volume; + assert!(rel(v, outer - inner) < 2e-3, "{v} vs {}", outer - inner); + // A cone frustum (a revolved trapezoid), opened at its wide base. + let prof = + Region2 { outer: Loop2::polygon(&[Vec2::new(0.0, 0.0), Vec2::new(10.0, 0.0), Vec2::new(5.0, -10.0), Vec2::new(0.0, -10.0)]), holes: vec![] }; + let fr = revolve(&xz, &[prof], Vec2::ZERO, Vec2::new(0.0, -1.0), 2.0 * PI).unwrap().pop().unwrap(); + let v0 = measure(&fr).unwrap().volume; + let s = shell(&fr, &[Vec3::new(0.0, 0.0, 0.0)], 1.0).unwrap(); + let v = measure(&s).unwrap().volume; + assert!(v > 0.1 * v0 && v < 0.6 * v0, "{v} of {v0}"); +}