mirror of
https://github.com/bherbruck/solidcraft.git
synced 2026-11-03 02:25:55 +00:00
solidcraft-smlib: revolved surfaces into SMLib as exact rational B-splines
A PlacedRevolution record goes to SMLib as its exact 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, so the directions swap. The B-spline matches the record's surface to 1e-12 with the same normal, forward or not, mirrored or not; truck's sphere and torus now reach SMLib with their analytic volume (1e-6) as manifold solids. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Brennen Herbruck
committed by
bherbruck
bherbruck
parent
bfb238e48b
commit
da7ccba9e4
@@ -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}");
|
||||
}
|
||||
|
||||
+106
-1
@@ -276,11 +276,116 @@ fn surface(a: &mut Arrays, s: &b::Surface, domain: [[f64; 2]; 2]) -> R<bool> {
|
||||
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<f64>, 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<([f64; 3], f64)>> = 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<Item = ([f64; 3], f64)>) {
|
||||
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<()> {
|
||||
|
||||
@@ -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}");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user