hoomd_geometry/shape/
cylinder.rs

1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement [`Cylinder`]
5
6use serde::{Deserialize, Serialize};
7
8use super::Circle;
9use crate::Volume;
10use hoomd_utility::valid::PositiveReal;
11use hoomd_vector::{Cartesian, Versor};
12
13/// A circle with normal `[0 0 1]` swept by `h/2` in the `+z` and `-z` directions.
14///
15/// # Example
16///
17/// [`Cylinder`] implements the [`Volume`] trait, which is equivalent to
18/// $` \pi r^2 h `$.
19/// ```
20/// use hoomd_geometry::{Volume, shape::Cylinder};
21/// use std::f64::consts::PI;
22///
23/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
24/// let cyl = Cylinder {
25///     radius: 2.0.try_into()?,
26///     height: 3.0.try_into()?,
27/// };
28/// assert_eq!(cyl.volume(), PI * (2.0 * 2.0) * 3.0);
29/// # Ok(())
30/// # }
31/// ```
32#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
33pub struct Cylinder {
34    /// Radius of the [`Cylinder`]
35    pub radius: PositiveReal,
36    /// Height of the [`Cylinder`]
37    pub height: PositiveReal,
38}
39
40impl Volume for Cylinder {
41    #[inline]
42    fn volume(&self) -> f64 {
43        Circle {
44            radius: self.radius,
45        }
46        .volume()
47            * self.height.get()
48    }
49}
50
51impl Cylinder {
52    /// Determine whether two cylinders would intersect if they both had infinite length.
53    ///
54    /// This implementation is based on [Collision detection of cylindrical rigid bodies for motion planning](https://doi.org/10.1109/ROBOT.2006.1641925), and provides a
55    /// useful alternative to bounding sphere-based checks when the underlying bodies
56    /// are highly elongated.
57    ///
58    /// # Example
59    ///
60    /// This example is based on the scenario of two highly anisotropic shapes,
61    /// whose bounding spheres have a large amount of volume relative to a tighter
62    /// fitting bounding shape. Note that the effective volume checked for the infinite
63    /// bounding cylinder is actually the intersection of that cylinder and the ball
64    /// queried in the neighbor list, so it is guaranteed to be finite.
65    ///
66    /// ```
67    /// # use hoomd_geometry::shape::Cylinder;
68    /// # use hoomd_vector::{Cartesian, Versor};
69    /// // Two parallel cylinders that wrap a dipyramid of height:width ratio 5:1.
70    /// let c1 = Cylinder {
71    ///     radius: 1.0.try_into().unwrap(),
72    ///     height: 10.0.try_into().unwrap(),
73    /// };
74    /// let c2 = Cylinder {
75    ///     radius: 1.0.try_into().unwrap(),
76    ///     height: 10.0.try_into().unwrap(),
77    /// };
78    /// let o_ij = Versor::default();
79    ///
80    /// // Case 1: Bounding spheres would intersect, but the bounding cylinders do not.
81    /// let v_ij_no_intersect = Cartesian::from([2.1, 0.0, 0.0]);
82    /// assert!(!c1.intersects_at_infinite(&c2, &v_ij_no_intersect, &o_ij));
83    ///
84    /// // Case 2: Infinite cylinders intersect, meaning we need to further check the
85    /// // underlying shapes for collision.
86    /// let v_ij_intersect = Cartesian::from([1.9, 0.0, 0.0]);
87    /// assert!(c1.intersects_at_infinite(&c2, &v_ij_intersect, &o_ij));
88    /// ```
89    #[must_use]
90    #[inline]
91    pub fn intersects_at_infinite(&self, other: &Self, v_ij: &Cartesian<3>, o_ij: &Versor) -> bool {
92        let [vx, vy, _vz] = v_ij.coordinates;
93
94        // Calculate sx and sy directly from the Versor components
95        // for rotating (0, 0, 1) by quaternion q = (w, x, y, z):
96        let q = o_ij.get();
97        let sx = 2.0 * (q.vector[0] * q.vector[2] + q.scalar * q.vector[1]);
98        let sy = 2.0 * (q.vector[1] * q.vector[2] - q.scalar * q.vector[0]);
99
100        let n_magnitude_sq = sx.powi(2) + sy.powi(2);
101
102        let dot_product = vy * sx - vx * sy;
103        let mut distance_sq = dot_product.powi(2) / n_magnitude_sq;
104        if !distance_sq.is_finite() {
105            distance_sq = vx.powi(2) + vy.powi(2);
106        }
107
108        // A collision occurs if the shortest distance between the axes is <= r1+r2
109        distance_sq <= (self.radius.get() + other.radius.get()).powi(2)
110    }
111}
112
113#[cfg(test)]
114mod tests {
115    use super::*;
116    use hoomd_vector::{Cartesian, InnerProduct, Versor};
117    use rstest::rstest;
118    use std::f64::consts::PI;
119
120    #[rstest]
121    #[case::parallel_intersecting(1.0, 1.0, [1.5, 0.0, 0.0], Versor::default(), true)]
122    #[case::parallel_touching(1.0, 1.000, [2.0, 0.0, 0.0], Versor::default(), true)]
123    #[case::parallel_barely_not_touching(1.0, 0.999_999, [2.0, 0.0, 0.0], Versor::default(), false)]
124    #[case::parallel_not_intersecting(1.0, 1.0, [2.000_001, 0.0, 0.0], Versor::default(), false)]
125    #[case::parallel_different_radii_touching(1.0, 0.5, [1.5, 0.0, 0.0], Versor::default(), true)]
126    #[case::parallel_different_radii_not_intersecting(1.0, 0.5, [1.500_001, 0.0, 0.0], Versor::default(), false)]
127    #[case::perpendicular_intersecting_at_origin(1.0, 1.0, [0.0, 0.0, 0.0], Versor::from_axis_angle(Cartesian::from([1., 0., 0.]).to_unit_unchecked().0, PI / 2.0), true)]
128    #[case::perpendicular_skew_intersecting(1.0, 1.0, [1.5, 0.0, 0.0], Versor::from_axis_angle(Cartesian::from([1., 0., 0.]).to_unit_unchecked().0, PI / 2.0), true)]
129    #[case::perpendicular_skew_touching(1.0, 1.0, [0.0, 2.0, 0.0], Versor::from_axis_angle(Cartesian::from([0., 1., 0.]).to_unit_unchecked().0, PI / 2.0), true)]
130    #[case::perpendicular_skew_barely_not_touching(1.0, 0.999_999, [0.0, 2.0, 0.0], Versor::from_axis_angle(Cartesian::from([0., 1., 0.]).to_unit_unchecked().0, PI / 2.0), false)]
131    #[case::perpendicular_skew_not_intersecting(1.0, 1.0, [2.000_001, 0.0, 5.0], Versor::from_axis_angle(Cartesian::from([1., 0., 0.]).to_unit_unchecked().0, PI / 2.0), false)]
132    #[case::skew_intersecting(1.0, 1.0, [1.0, 1.0, 0.0], Versor::from_axis_angle(Cartesian::from([0., 1., 0.]).to_unit_unchecked().0, PI / 5.0), true)]
133    #[case::skew_touching(1.0, 1.0, [0.0, 2.0, 0.0], Versor::from_axis_angle(Cartesian::from([0., 1., 0.]).to_unit_unchecked().0, PI / 5.0), true)]
134    #[case::skew_not_intersecting(1.0, 1.0, [0.0, 2.000_001, 0.0], Versor::from_axis_angle(Cartesian::from([0., 1., 0.]).to_unit_unchecked().0, PI / 5.0), false)]
135    fn test_intersects_at_infinite(
136        #[case] r1: f64,
137        #[case] r2: f64,
138        #[case] v_ij: impl Into<Cartesian<3>>,
139        #[case] o_ij: Versor,
140        #[case] expected: bool,
141    ) -> Result<(), Box<dyn std::error::Error>> {
142        let c1 = Cylinder {
143            radius: r1.try_into()?,
144            height: 1.0.try_into()?,
145        };
146        let c2 = Cylinder {
147            radius: r2.try_into()?,
148            height: 1.0.try_into()?,
149        };
150        let v_ij_cartesian = v_ij.into();
151
152        assert_eq!(
153            c1.intersects_at_infinite(&c2, &v_ij_cartesian, &o_ij),
154            expected,
155        );
156
157        Ok(())
158    }
159}
hoomd_geometry/shape/cylinder.rs

hoomd_geometry/shape/
cylinder.rs
