diff --git a/crates/kernel/src/native/smlib_import_tests.rs b/crates/kernel/src/native/smlib_import_tests.rs index ba67214cd3..0443c3531e 100644 --- a/crates/kernel/src/native/smlib_import_tests.rs +++ b/crates/kernel/src/native/smlib_import_tests.rs @@ -27,12 +27,14 @@ fn primitives_keep_their_analytic_volume() { assert!(solid && (got - 28.0 * PI).abs() <= 1e-6 * 28.0 * PI, "tilted cylinder: SMLib volume {got}, manifold {solid}"); } -/// Truck holds spheres and tori as revolutions (`PlacedRevolution`), which the SMLib import does -/// not carry yet: a clear error, never a wrong body. +/// Truck holds spheres and tori as revolutions (`PlacedRevolution`); they go to SMLib as their +/// exact rational B-spline surfaces. #[test] -fn revolutions_are_refused_for_now() { - for (name, body) in [("sphere", sphere(v(0., 0., 0.), 5.0).unwrap()), ("torus", torus(v(0., 0., 0.), 10.0, 2.0).unwrap())] { - let e = super::smlib_body(&body).err().map(|e| e.to_string()).unwrap_or_default(); - assert!(e.contains("not supported yet: revolved surfaces"), "{name}: {e}"); - } +fn revolutions_keep_their_analytic_volume() { + let (got, solid) = in_smlib(&sphere(v(1., 2., 3.), 5.0).unwrap()).unwrap(); + let ball = 4.0 / 3.0 * PI * 125.0; + assert!(solid && (got - ball).abs() <= 1e-6 * ball, "sphere: SMLib volume {got}, manifold {solid}"); + let (got, solid) = in_smlib(&torus(v(0., 0., 0.), 10.0, 2.0).unwrap()).unwrap(); + let ring = 2.0 * PI * PI * 10.0 * 4.0; + assert!(solid && (got - ring).abs() <= 1e-6 * ring, "torus: SMLib volume {got}, manifold {solid}"); } diff --git a/crates/smlib/src/import.rs b/crates/smlib/src/import.rs index 6eb65d8426..c8e1039492 100644 --- a/crates/smlib/src/import.rs +++ b/crates/smlib/src/import.rs @@ -276,11 +276,116 @@ fn surface(a: &mut Arrays, s: &b::Surface, domain: [[f64; 2]; 2]) -> R { a.r(flat::FACE_RANGE).extend([u0, v0, u1, v1]); return Ok(true); } - b::Surface::PlacedRevolution { .. } => return Err(Unsupported::NotYet("revolved surfaces")), + b::Surface::PlacedRevolution { profile, origin, axis, matrix, forward } => revolution(a, profile, *origin, *axis, matrix, *forward, domain)?, } Ok(false) } +/// The quarter-circle arc's rational parameter `s` in [0, 1] at angle `phi` in [0, π/2]: the +/// circle NURBS below is not linear in the angle, so end values are converted pointwise. +fn quarter_param(phi: f64) -> f64 { + let c = std::f64::consts::FRAC_1_SQRT_2; + let angle = |s: f64| { + let (b0, b1, b2) = ((1.0 - s) * (1.0 - s), 2.0 * s * (1.0 - s) * c, s * s); + (b1 + b2).atan2(b0 + b1) + }; + let (mut lo, mut hi) = (0.0f64, 1.0f64); + for _ in 0..64 { + let mid = 0.5 * (lo + hi); + if angle(mid) < phi { + lo = mid; + } else { + hi = mid; + } + } + 0.5 * (lo + hi) +} + +/// The circle NURBS's parameter (knots 0..4, one quarter turn per span) at `angle` after its start. +fn circle_param(angle: f64) -> f64 { + let q = std::f64::consts::FRAC_PI_2; + let a = angle.clamp(0.0, 4.0 * q); + let k = (a / q).floor().min(3.0); + k + quarter_param(a - k * q) +} + +/// A revolution, exactly, as a rational B-spline surface: each profile control point swept round +/// the axis by a full circle of four rational quadratic quarters starting at the face's first +/// angle, then mapped by the placement. The circle's parameter is not the angle, so the face's +/// angle range is converted pointwise and its pcurves are not carried (SMLib rebuilds the trims). +/// Where the placement mirrors, the record's natural normal is the opposite of the map's +/// (`brep::eval`), so the B-spline takes the two directions the other way round. +fn revolution( + a: &mut Arrays, + profile: &b::Curve, + origin: [f64; 3], + axis: [f64; 3], + m: &[[f64; 4]; 4], + forward: bool, + domain: [[f64; 2]; 2], +) -> R<()> { + let (knots, pts): (Vec, Vec<([f64; 3], f64)>) = match profile { + b::Curve::Line { start, end } => (vec![0.0, 0.0, 1.0, 1.0], vec![(*start, 1.0), (*end, 1.0)]), + b::Curve::Spline { knots, points } => (knots.clone(), points.iter().map(|p| (*p, 1.0)).collect()), + b::Curve::RationalSpline { knots, homogeneous } => { + (knots.clone(), homogeneous.iter().map(|h| ([h[0] / h[3], h[1] / h[3], h[2] / h[3]], h[3])).collect()) + } + _ => return Err(Unsupported::NotYet("revolutions of conic, polyline or intersection profiles")), + }; + let len = dot(axis, axis).sqrt(); + if !(len > 0.0) { + return Err(Unsupported::Invalid("revolution axis")); + } + let z = scale(axis, 1.0 / len); + // The record's u and v: the profile's parameter and the angle, or the other way round. + let ([p0, p1], [a0, a1]) = if forward { (domain[0], domain[1]) } else { (domain[1], domain[0]) }; + let t1 = circle_param(a1 - a0); + // Column-major placement: point' = M · (p, 1). + let place = |p: [f64; 3]| [0, 1, 2].map(|r| m[0][r] * p[0] + m[1][r] * p[1] + m[2][r] * p[2] + m[3][r]); + let det = { + let col = |c: usize| [m[c][0], m[c][1], m[c][2]]; + dot(cross(col(0), col(1)), col(2)) + }; + let c = std::f64::consts::FRAC_1_SQRT_2; + // net[i][k]: profile point i at circle control point k. + let mut net: Vec> = Vec::with_capacity(pts.len()); + for (p, w) in &pts { + let r = sub(*p, origin); + let centre = add(origin, scale(z, dot(r, z))); + let x0 = sub(*p, centre); + let y0 = cross(z, x0); + // Start at the face's first angle. + let (s0, c0) = a0.sin_cos(); + let (x, y) = (add(scale(x0, c0), scale(y0, s0)), sub(scale(y0, c0), scale(x0, s0))); + let row = (0..9) + .map(|k| { + let (sk, ck) = (k as f64 * std::f64::consts::FRAC_PI_4).sin_cos(); + let (f, wk) = if k % 2 == 1 { (1.0 / c, c) } else { (1.0, 1.0) }; + (place(add(centre, scale(add(scale(x, ck), scale(y, sk)), f))), w * wk) + }) + .collect(); + net.push(row); + } + let circle_knots = [0.0, 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 4.0, 4.0, 4.0]; + // Profile along the B-spline's first direction when the record's u is the profile and the + // placement keeps orientation, or the angle is u and it mirrors; else the angle goes first. + let profile_first = forward == (det > 0.0); + if profile_first { + nurbs_face(a, &knots, &circle_knots, pts.len(), 9)?; + net_rows(a, net.iter().flatten().copied()); + a.r(flat::FACE_RANGE).extend([p0, 0.0, p1, t1]); + } else { + nurbs_face(a, &circle_knots, &knots, 9, pts.len())?; + net_rows(a, (0..9).flat_map(|k| net.iter().filter_map(move |row| row.get(k).copied()))); + a.r(flat::FACE_RANGE).extend([0.0, p0, t1, p1]); + } + Ok(()) +} + +fn net_rows(a: &mut Arrays, cps: impl Iterator) { + net(a, cps); +} + /// A plane `origin + u·U + v·V`. With unit, orthogonal axes it is SMLib's plane, whose (u, v) are /// distances along them; otherwise it is the exact bilinear patch over the domain. fn plane(a: &mut Arrays, origin: [f64; 3], u: [f64; 3], v: [f64; 3], domain: [[f64; 2]; 2]) -> R<()> { diff --git a/crates/smlib/src/import/tests.rs b/crates/smlib/src/import/tests.rs index 32a43f3393..8f286c3bd9 100644 --- a/crates/smlib/src/import/tests.rs +++ b/crates/smlib/src/import/tests.rs @@ -150,3 +150,73 @@ fn a_skewed_plane_goes_as_an_exact_bilinear_patch() { let b = crate::Brep::from_record(&r, true).unwrap(); assert!((b.volume(1e-7).unwrap() - 36.0).abs() < 1e-6, "{:?}", b.volume(1e-7)); } + +/// The surface `arrays` wrote for the first face, back as a schema 4 rational spline, with its range. +fn first_nurbs(a: &Arrays) -> (Surface, [f64; 4]) { + let (nu, nv) = (a.ints[&flat::FACE_NURBS_U_VERTEX_COUNT][0] as usize, a.ints[&flat::FACE_NURBS_V_VERTEX_COUNT][0] as usize); + let cvs = &a.reals[&flat::FACE_NURBS_CVS]; + let w = &a.reals[&flat::FACE_NURBS_WEIGHTS]; + let homogeneous = (0..nu) + .map(|i| { + (0..nv) + .map(|j| { + let k = i * nv + j; + [cvs[3 * k] * w[k], cvs[3 * k + 1] * w[k], cvs[3 * k + 2] * w[k], w[k]] + }) + .collect() + }) + .collect(); + let r = &a.reals[&flat::FACE_RANGE]; + ( + Surface::RationalSpline { + u_knots: a.reals[&flat::FACE_NURBS_U_KNOTS].clone(), + v_knots: a.reals[&flat::FACE_NURBS_V_KNOTS].clone(), + homogeneous, + }, + [r[0], r[1], r[2], r[3]], + ) +} + +#[test] +fn revolutions_become_the_same_surface_as_a_rational_spline() { + use solidcraft_geom::brep::eval::SurfaceEval; + // A quarter-circle profile (rational) away from the axis, so the surface is a torus piece. + let s = std::f64::consts::FRAC_1_SQRT_2; + let profile = Curve::RationalSpline { + knots: vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0], + homogeneous: vec![[3.0, 0.0, 0.0, 1.0], [3.0 * s, 0.0, 1.0 * s, s], [2.0, 0.0, 1.0, 1.0]], + }; + let identity = [[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 0.0, 1.0]]; + let mirror = [[0.0, 1.0, 0.0, 0.0], [1.0, 0.0, 0.0, 0.0], [0.0, 0.0, 2.0, 0.0], [1.0, -2.0, 3.0, 1.0]]; + for matrix in [identity, mirror] { + for forward in [true, false] { + let rev = + Surface::PlacedRevolution { profile: Box::new(profile.clone()), origin: [0.0, 0.0, 0.0], axis: [0.0, 0.0, 1.0], matrix, forward }; + let (angles, prof) = ([0.4, 0.4 + 1.9 * std::f64::consts::PI], [0.0, 1.0]); + let domain = if forward { [prof, angles] } else { [angles, prof] }; + let mut a = Arrays::default(); + surface(&mut a, &rev, domain).unwrap(); + let (nurbs, range) = first_nurbs(&a); + let (re, ne) = (SurfaceEval::new(&rev).unwrap(), SurfaceEval::new(&nurbs).unwrap()); + let det = { + let col = |c: usize| [matrix[c][0], matrix[c][1], matrix[c][2]]; + dot(cross(col(0), col(1)), col(2)) + }; + for i in 0..=6 { + for j in 0..=6 { + let (sp, ang) = (i as f64 / 6.0, angles[0] + (angles[1] - angles[0]) * j as f64 / 6.0); + let uv = if forward { [sp, ang] } else { [ang, sp] }; + let (p, du, dv) = re.derivs(uv[0], uv[1]).unwrap(); + let n = du.cross(dv) * det.signum(); + let t = circle_param(ang - angles[0]); + let profile_first = forward == (det > 0.0); + let q = if profile_first { [sp, t] } else { [t, sp] }; + assert!(q[0] >= range[0] - 1e-12 && q[0] <= range[2] + 1e-12 && q[1] >= range[1] - 1e-12 && q[1] <= range[3] + 1e-12); + let (pn, nu_, nv_) = ne.derivs(q[0], q[1]).unwrap(); + assert!((pn - p).len() < 1e-12, "forward {forward} det {det}: {p:?} vs {pn:?}"); + assert!(nu_.cross(nv_).dot(n) > 0.0, "normal flipped: forward {forward} det {det}"); + } + } + } + } +}