Skip to main content

axiolid_construct/
revolve.rs

1//! Revolution of a profile about an axis.
2//!
3//! The profile is flattened to rings in its own plane, then each ring point
4//! is rotated about the axis in steps chosen by the chord budget. A full
5//! turn closes on itself; a partial turn is capped by the profile at each
6//! end.
7//!
8//! Pappus gives the oracle: a full revolution has volume `2*pi*R*A`, where
9//! `R` is the centroid's distance from the axis and `A` the profile area.
10
11use axiolid_contracts::{GeomError, GeomResult};
12use axiolid_core::{Point3, Scalar, Tolerance, Vec3};
13use axiolid_mesh::TriMesh;
14
15use crate::profile::Rings;
16
17/// Rotation of `p` about the axis through `origin` along unit `dir`.
18///
19/// Rodrigues' formula. Written out rather than pulled from a matrix type
20/// so the axis stays arbitrary: a revolution is not a Z-up operation.
21pub(crate) fn rotate(p: Point3, origin: Point3, dir: Vec3, angle: Scalar) -> Point3 {
22    let v = p - origin;
23    let (s, c) = angle.sin_cos();
24    origin + v * c + dir.cross(v) * s + dir * (dir.dot(v) * (1.0 - c))
25}
26
27/// Angular steps so the swept arc meets the chord budget.
28///
29/// The widest point of the profile governs: a ring point at distance `r`
30/// from the axis traces a circle of that radius, and its sagitta is
31/// `r(1 - cos(dtheta/2))`. Using the maximum radius means every other
32/// point is sampled at least as finely.
33pub(crate) fn steps(max_radius: Scalar, angle: Scalar, tol: Scalar) -> usize {
34    if !(max_radius.is_finite() && max_radius > 0.0 && tol.is_finite() && tol > 0.0) {
35        return 8;
36    }
37    let ratio = (1.0 - (tol / max_radius).min(1.0)).clamp(-1.0, 1.0);
38    let per = 2.0 * ratio.acos().max(1e-9);
39    ((angle.abs() / per).ceil() as usize).clamp(2, 4096)
40}
41
42/// Revolve a profile about an axis into a closed solid.
43///
44/// A full turn wraps its rings by index, so the seam shares vertices by
45/// construction rather than by two samplings agreeing numerically. That is
46/// the same structure `tessellate_primitive` uses for a cylinder, and it
47/// avoids the trim-based path recorded as broken in issue #2.
48pub fn revolve(
49    rings: &Rings,
50    axis_origin: Point3,
51    axis_direction: Vec3,
52    angle: Scalar,
53    tolerance: Tolerance,
54) -> GeomResult<TriMesh> {
55    if !angle.is_finite() || angle == 0.0 {
56        return Err(GeomError::InvalidInput(format!(
57            "revolution angle must be finite and non-zero, got {angle}"
58        )));
59    }
60    let len = axis_direction.length();
61    if !len.is_finite() || len <= 0.0 {
62        return Err(GeomError::InvalidInput(
63            "revolution axis must be a finite non-zero direction".to_owned(),
64        ));
65    }
66    let dir = axis_direction / len;
67    // Widest distance from the axis governs the angular step: a point at
68    // radius r traces a circle of that radius.
69    let mut max_r: Scalar = 0.0;
70    for p in rings.outer.iter().chain(rings.holes.iter().flatten()) {
71        let v = Point3::new(p.x, p.y, 0.0) - axis_origin;
72        max_r = max_r.max((v - dir * dir.dot(v)).length());
73    }
74    let full = (angle.abs() - core::f64::consts::TAU).abs() <= 1e-9;
75    let n = steps(max_r, angle, tolerance.linear());
76    // A full turn emits n stations and wraps onto the first; a partial turn
77    // emits n + 1 so both ends exist to be capped.
78    let count = if full { n } else { n + 1 };
79    let stations: Vec<crate::loft::Station> = (0..count)
80        .map(|s| {
81            let t = angle * (s as Scalar) / (n as Scalar);
82            crate::loft::place(rings, |p| {
83                rotate(Point3::new(p.x, p.y, 0.0), axis_origin, dir, t)
84            })
85        })
86        .collect();
87    let stations: Vec<_> = stations.into_iter().rev().collect();
88    crate::loft::loft(rings, &stations, full)
89}