hoomd_geometry/shape/
sphere.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 [`Hypersphere`]
5
6use rand::{Rng, distr::Distribution};
7use serde::{Deserialize, Serialize};
8use std::{array, f64::consts::PI, ops::Mul};
9
10use crate::{
11    BoundingSphereRadius, Error, IntersectsAt, IntersectsAtGlobal, IsPointInside, MapPoint, Scale,
12    SupportMapping, Volume,
13};
14use hoomd_utility::valid::PositiveReal;
15use hoomd_vector::{Cartesian, InnerProduct, Rotate, Rotation, distribution::Ball};
16
17/// The (single, double, ...)-factorial function
18pub fn factorial(n: usize, ntuple: usize) -> usize {
19    assert!(ntuple > 0);
20    if n == 0 {
21        1
22    } else {
23        (1..=n)
24            .rev()
25            .step_by(ntuple)
26            .reduce(usize::mul)
27            .expect("1..=(n!=0) is never empty")
28    }
29}
30
31/// Compute the volume prefactor for the volume of a rounded shape
32pub(crate) fn sphere_volume_prefactor(n: usize) -> f64 {
33    // FUTURE: replace with std::f64::gamma when its in main
34    let dim_factor = (if n.rem_euclid(2) == 0 { n } else { n - 1 } / 2) as f64;
35    if n.rem_euclid(2) == 0 {
36        PI.powf(dim_factor) / (factorial(n / 2, 1) as f64)
37    } else {
38        2.0 * (2.0 * PI).powf(dim_factor) / (factorial(n, 2) as f64)
39    }
40}
41
42/// All points within a given `radius` from the origin.
43///
44/// # Examples
45///
46/// Basic construction and methods:
47/// ```
48/// use hoomd_geometry::{SupportMapping, Volume, shape::Hypersphere};
49/// use hoomd_vector::Cartesian;
50/// use std::f64::consts::PI;
51///
52/// let unit_sphere = Hypersphere::<3>::default();
53/// let volume = unit_sphere.volume();
54///
55/// assert_eq!(unit_sphere.radius.get(), 1.0);
56/// assert_eq!(volume, 4.0 * PI / 3.0);
57///
58/// assert_eq!(
59///     unit_sphere.support_mapping(&Cartesian::from([1.0; 3])),
60///     [1.0 / f64::sqrt(3.0); 3].into()
61/// )
62/// ```
63///
64/// Test for intersections:
65/// ```
66/// use hoomd_geometry::{IntersectsAt, shape::Hypersphere};
67/// use hoomd_vector::{Cartesian, Versor};
68///
69/// let unit_sphere = Hypersphere::<3>::default();
70///
71/// assert!(!unit_sphere.intersects_at(
72///     &unit_sphere,
73///     &Cartesian::from([2.1, 0.0, 0.0]),
74///     &Versor::default()
75/// ));
76/// assert!(unit_sphere.intersects_at(
77///     &unit_sphere,
78///     &Cartesian::from([0.0, 1.9, 0.0]),
79///     &Versor::default()
80/// ));
81/// ```
82#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
83pub struct Hypersphere<const N: usize> {
84    /// Radius of the sphere
85    pub radius: PositiveReal,
86}
87
88/// A circle in two dimensions.
89///
90/// # Examples
91///
92/// Basic construction and methods:
93/// ```
94/// use hoomd_geometry::{Volume, shape::Circle};
95/// use hoomd_vector::Cartesian;
96/// use std::f64::consts::PI;
97///
98/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
99/// let circle = Circle {
100///     radius: 2.0.try_into()?,
101/// };
102/// let volume = circle.volume();
103///
104/// assert_eq!(volume, PI * 4.0);
105/// # Ok(())
106/// # }
107/// ```
108///
109/// Test for intersections:
110/// ```
111/// use hoomd_geometry::{IntersectsAt, shape::Circle};
112/// use hoomd_vector::{Angle, Cartesian};
113///
114/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
115/// let circle = Circle {
116///     radius: 1.0.try_into()?,
117/// };
118///
119/// assert_eq!(
120///     circle.intersects_at(
121///         &circle,
122///         &Cartesian::from([2.1, 0.0]),
123///         &Angle::default()
124///     ),
125///     false
126/// );
127/// assert_eq!(
128///     circle.intersects_at(
129///         &circle,
130///         &Cartesian::from([0.0, 1.9]),
131///         &Angle::default()
132///     ),
133///     true
134/// );
135/// # Ok(())
136/// # }
137/// ```
138pub type Circle = Hypersphere<2>;
139
140/// A sphere in three dimensions.
141///
142/// # Examples
143///
144/// Basic construction and methods:
145/// ```
146/// use hoomd_geometry::{SupportMapping, Volume, shape::Sphere};
147/// use hoomd_vector::Cartesian;
148/// use std::f64::consts::PI;
149///
150/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
151/// let unit_sphere = Sphere {
152///     radius: 1.0.try_into()?,
153/// };
154/// let volume = unit_sphere.volume();
155///
156/// assert_eq!(unit_sphere.radius.get(), 1.0);
157/// assert_eq!(volume, 4.0 * PI / 3.0);
158///
159/// assert_eq!(
160///     unit_sphere.support_mapping(&Cartesian::from([1.0; 3])),
161///     [1.0 / f64::sqrt(3.0); 3].into()
162/// );
163/// # Ok(())
164/// # }
165/// ```
166///
167/// Test for intersections:
168/// ```
169/// use hoomd_geometry::{IntersectsAt, shape::Sphere};
170/// use hoomd_vector::{Cartesian, Versor};
171///
172/// let unit_sphere = Sphere::default();
173///
174/// assert!(!unit_sphere.intersects_at(
175///     &unit_sphere,
176///     &Cartesian::from([2.1, 0.0, 0.0]),
177///     &Versor::default()
178/// ));
179/// assert!(unit_sphere.intersects_at(
180///     &unit_sphere,
181///     &Cartesian::from([0.0, 1.9, 0.0]),
182///     &Versor::default()
183/// ));
184/// ```
185pub type Sphere = Hypersphere<3>;
186
187impl<const N: usize> Default for Hypersphere<N> {
188    #[inline]
189    fn default() -> Self {
190        Hypersphere {
191            radius: 1.0.try_into().expect("1.0 is a positive real"),
192        }
193    }
194}
195
196impl<const N: usize> Hypersphere<N> {
197    /// Create a sphere with a given positive real radius.
198    #[must_use]
199    #[inline]
200    pub fn with_radius(radius: PositiveReal) -> Self {
201        Hypersphere { radius }
202    }
203
204    /// Test whether one sphere intersects with another.
205    ///
206    /// The vector `v_ij` points from the local origin of `self` to the local
207    /// origin of `other`.
208    #[inline]
209    pub fn intersects<V>(&self, other: &Hypersphere<N>, v_ij: &V) -> bool
210    where
211        V: InnerProduct,
212    {
213        (v_ij).norm_squared() <= (other.radius.get() + self.radius.get()).powi(2)
214    }
215}
216
217// TRAITS
218
219impl<const N: usize, V: InnerProduct> SupportMapping<V> for Hypersphere<N> {
220    #[inline]
221    fn support_mapping(&self, n: &V) -> V {
222        *n / n.norm() * self.radius.get()
223    }
224}
225
226impl<const N: usize> Volume for Hypersphere<N> {
227    #[inline]
228    fn volume(&self) -> f64 {
229        sphere_volume_prefactor(N)
230            * self
231                .radius
232                .get()
233                .powi(N.try_into().expect("Dimension should not overflow i32!"))
234    }
235}
236
237impl<const N: usize, R> IntersectsAtGlobal<Hypersphere<N>, Cartesian<N>, R> for Hypersphere<N>
238where
239    R: Rotation + Rotate<Cartesian<N>>,
240{
241    #[inline]
242    fn intersects_at_global(
243        &self,
244        other: &Hypersphere<N>,
245        r_self: &Cartesian<N>,
246        o_self: &R,
247        r_other: &Cartesian<N>,
248        o_other: &R,
249    ) -> bool {
250        let (v_ij, o_ij) = hoomd_vector::pair_system_to_local(r_self, o_self, r_other, o_other);
251
252        self.intersects_at(other, &v_ij, &o_ij)
253    }
254}
255
256impl<const N: usize, V, R> IntersectsAt<Hypersphere<N>, V, R> for Hypersphere<N>
257where
258    V: InnerProduct,
259    R: Rotation + Rotate<V>,
260{
261    #[inline]
262    fn intersects_at(&self, other: &Hypersphere<N>, v_ij: &V, _o_ij: &R) -> bool {
263        (v_ij).norm_squared() <= (other.radius.get() + self.radius.get()).powi(2)
264    }
265}
266
267impl<const N: usize> BoundingSphereRadius for Hypersphere<N> {
268    #[inline]
269    fn bounding_sphere_radius(&self) -> PositiveReal {
270        self.radius
271    }
272}
273
274impl<const N: usize, V> IsPointInside<V> for Hypersphere<N>
275where
276    V: InnerProduct,
277{
278    /// Check if a vector is inside a hypersphere.
279    ///
280    /// ```
281    /// use hoomd_geometry::{IsPointInside, shape::Sphere};
282    /// use hoomd_vector::Cartesian;
283    ///
284    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
285    /// let sphere = Sphere {
286    ///     radius: 3.0.try_into()?,
287    /// };
288    ///
289    /// assert!(sphere.is_point_inside(&Cartesian::from([2.5, 0.0, 0.0])));
290    /// assert!(!sphere.is_point_inside(&Cartesian::from([3.0, -3.0, 2.0])));
291    /// # Ok(())
292    /// # }
293    /// ```
294    #[inline]
295    fn is_point_inside(&self, point: &V) -> bool {
296        point.dot(point) < self.radius.get().powi(2)
297    }
298}
299
300impl<const N: usize> Scale for Hypersphere<N> {
301    /// Construct a scaled hypersphere.
302    ///
303    /// The resulting hypersphere's radious $` r_\mathrm{new} `$ is
304    /// the original's $` r `$ scaled by $` v `$:
305    /// ```math
306    /// r_\mathrm{new} = v r
307    /// ```
308    ///
309    /// The centroid remains at the origin.
310    ///
311    /// # Example
312    ///
313    /// ```
314    /// use hoomd_geometry::{Scale, shape::Sphere};
315    ///
316    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
317    /// let sphere = Sphere {
318    ///     radius: 5.0.try_into()?,
319    /// };
320    ///
321    /// let scaled_sphere = sphere.scale_length(0.5.try_into()?);
322    ///
323    /// assert_eq!(scaled_sphere.radius.get(), 2.5);
324    /// # Ok(())
325    /// # }
326    /// ```
327    #[inline]
328    fn scale_length(&self, v: PositiveReal) -> Self {
329        Self {
330            radius: self.radius * v,
331        }
332    }
333
334    /// Construct a scaled hypersphere.
335    ///
336    /// The resulting hypersphere's radius $` r_\mathrm{new} `$ is
337    /// the original's $` r `$ scaled by $` v^\frac{1}{N} `$:
338    /// ```math
339    /// r_\mathrm{new} = v^\frac{1}{N} r
340    /// ```
341    ///
342    /// The centroid remains at the origin.
343    ///
344    /// # Example
345    ///
346    /// # Example
347    ///
348    /// ```
349    /// use hoomd_geometry::{Scale, shape::Circle};
350    ///
351    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
352    /// let sphere = Circle {
353    ///     radius: 5.0.try_into()?,
354    /// };
355    ///
356    /// let scaled_sphere = sphere.scale_volume(0.25.try_into()?);
357    ///
358    /// assert_eq!(scaled_sphere.radius.get(), 2.5);
359    /// # Ok(())
360    /// # }
361    /// ```
362    #[inline]
363    fn scale_volume(&self, v: PositiveReal) -> Self {
364        let v = v.get().powf(1.0 / N as f64);
365        self.scale_length(v.try_into().expect("v^{1/N} should be a positive real"))
366    }
367}
368
369impl<const N: usize> MapPoint<Cartesian<N>> for Hypersphere<N> {
370    /// Map a point from one hypersphere to another.
371    ///
372    /// Given a point P *inside `self`*, map it to the other shape
373    /// by scaling.
374    ///
375    /// # Errors
376    ///
377    /// Returns [`Error::PointOutsideShape`] when `point` is outside the shape
378    /// `self`.
379    ///
380    /// # Example
381    ///
382    /// ```
383    /// use hoomd_geometry::{MapPoint, shape::Circle};
384    /// use hoomd_vector::Cartesian;
385    ///
386    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
387    /// let closed_a = Circle {
388    ///     radius: 10.0.try_into()?,
389    /// };
390    /// let closed_b = Circle {
391    ///     radius: 20.0.try_into()?,
392    /// };
393    ///
394    /// let mapped_point =
395    ///     closed_a.map_point(Cartesian::from([-1.0, 1.0]), &closed_b);
396    ///
397    /// assert_eq!(mapped_point, Ok(Cartesian::from([-2.0, 2.0])));
398    /// assert_eq!(
399    ///     closed_a.map_point(Cartesian::from([-100.0, 1.0]), &closed_b),
400    ///     Err(hoomd_geometry::Error::PointOutsideShape)
401    /// );
402    /// # Ok(())
403    /// # }
404    /// ```
405    #[inline]
406    fn map_point(&self, point: Cartesian<N>, other: &Self) -> Result<Cartesian<N>, Error> {
407        if !self.is_point_inside(&point) {
408            return Err(Error::PointOutsideShape);
409        }
410
411        // When multiplying by the scale factor, the floating point multiply
412        // might round up in a way that places the mapped point outside the
413        // other shape. Progressively make the scale smaller until the
414        // check passes.
415        let mut scale = other.radius / self.radius;
416        loop {
417            let point = Cartesian::from(array::from_fn(|i| scale.get() * point[i]));
418            if other.is_point_inside(&point) {
419                return Ok(point);
420            }
421
422            scale = scale
423                .get()
424                .next_down()
425                .try_into()
426                .expect("scale should remain a positive real");
427        }
428    }
429}
430
431impl<const N: usize> Distribution<Cartesian<N>> for Hypersphere<N> {
432    /// Generate points uniformly distributed in the hypersphere.
433    ///
434    /// # Example
435    ///
436    /// ```
437    /// use hoomd_geometry::{IsPointInside, shape::Sphere};
438    /// use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
439    ///
440    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
441    /// let sphere = Sphere {
442    ///     radius: 5.0.try_into()?,
443    /// };
444    /// let mut rng = StdRng::seed_from_u64(1);
445    ///
446    /// let point = sphere.sample(&mut rng);
447    /// assert!(sphere.is_point_inside(&point));
448    /// # Ok(())
449    /// # }
450    /// ```
451    #[inline]
452    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Cartesian<N> {
453        let ball = Ball {
454            radius: self.radius,
455        };
456        ball.sample(rng)
457    }
458}
459
460#[cfg(test)]
461#[expect(
462    clippy::used_underscore_binding,
463    reason = "Used in test parameterization."
464)]
465#[expect(
466    clippy::unreadable_literal,
467    reason = "exact test results need not be readable"
468)]
469mod tests {
470    use super::*;
471    use crate::Convex;
472    use approxim::assert_relative_eq;
473    use assert2::check;
474    use hoomd_vector::{Cartesian, Versor};
475    use rand::{
476        SeedableRng,
477        distr::{Distribution, Uniform},
478        rngs::StdRng,
479    };
480    use rstest::*;
481    use std::marker::PhantomData;
482
483    /// Number of random samples to test.
484    const N: usize = 1024;
485
486    fn volume_map(n: usize, r: f64) -> f64 {
487        match n {
488            0 => 1.0 * r.powi(i32::try_from(n).unwrap()),
489            1 => 2.0 * r.powi(i32::try_from(n).unwrap()),
490            2 => PI * r.powi(i32::try_from(n).unwrap()),
491            3 => 4.0 / 3.0 * PI * r.powi(i32::try_from(n).unwrap()),
492            4 => PI.powi(2) / 2.0 * r.powi(i32::try_from(n).unwrap()),
493            5 => 8.0 * PI.powi(2) / 15.0 * r.powi(i32::try_from(n).unwrap()),
494            _ => unreachable!(),
495        }
496    }
497
498    #[rstest]
499    #[case(PhantomData::<Hypersphere<0>>)]
500    #[case(PhantomData::<Hypersphere<1>>)]
501    #[case(PhantomData::<Hypersphere<2>>)]
502    #[case(PhantomData::<Hypersphere<3>>)]
503    #[case(PhantomData::<Hypersphere<4>>)]
504    #[case(PhantomData::<Hypersphere<5>>)]
505    fn test_volume_and_radius<const N: usize>(
506        #[case] _n: PhantomData<Hypersphere<N>>,
507        #[values(0.01, 1.0, 33.3, 1e6)] radius: f64,
508    ) -> anyhow::Result<()> {
509        let s = Hypersphere::<N> {
510            radius: radius.try_into()?,
511        };
512
513        if radius == 1.0 {
514            check!(s.radius.get() == 1.0);
515            check!(s == Hypersphere::<N>::default());
516        } else {
517            check!(s.radius.get() == radius);
518        }
519
520        assert_relative_eq!(s.volume(), volume_map(N, radius));
521
522        Ok(())
523    }
524
525    #[rstest]
526    fn test_n_factorial(#[values(1, 2, 3, 4)] m: usize) {
527        check!(factorial(m, m) == m);
528    }
529    #[rstest]
530    fn test_single_double_factorial(#[values(1, 5, 10, 18, 20)] n: usize) {
531        check!(factorial(n, 1) == factorial(n, 2) * factorial(n - 1, 2));
532    }
533
534    #[rstest]
535    #[case(PhantomData::<Hypersphere<0>>)]
536    #[case(PhantomData::<Hypersphere<1>>)]
537    #[case(PhantomData::<Hypersphere<2>>)]
538    #[case(PhantomData::<Hypersphere<3>>)]
539    #[case(PhantomData::<Hypersphere<4>>)]
540    #[case(PhantomData::<Hypersphere<5>>)]
541    fn test_support_fn<const N: usize>(
542        #[case] _n: PhantomData<Hypersphere<N>>,
543        #[values(0.1, 1.0, 33.3)] radius: f64,
544    ) -> anyhow::Result<()> {
545        let s = Hypersphere::<N> {
546            radius: radius.try_into()?,
547        };
548        let v = Cartesian::<N>::from([radius.powi(2) / 1.8; N]);
549        check!(v / v.norm() * radius == s.support_mapping(&v));
550
551        Ok(())
552    }
553
554    #[test]
555    fn support_mapping() -> anyhow::Result<()> {
556        let sphere = Sphere::with_radius(2.0.try_into()?);
557
558        assert_relative_eq!(
559            sphere.support_mapping(&Cartesian::from([0.0, 0.0, 1.0])),
560            [0.0, 0.0, 2.0].into()
561        );
562        assert_relative_eq!(
563            sphere.support_mapping(&Cartesian::from([0.0, 0.0, 01.0])),
564            [0.0, 0.0, 2.0].into()
565        );
566        assert_relative_eq!(
567            sphere.support_mapping(&Cartesian::from([0.0, 1.0, 0.0])),
568            [0.0, 2.0, 0.0].into()
569        );
570        assert_relative_eq!(
571            sphere.support_mapping(&Cartesian::from([0.0, -1.0, 0.0])),
572            [0.0, -2.0, 0.0].into()
573        );
574        assert_relative_eq!(
575            sphere.support_mapping(&Cartesian::from([1.0, 0.0, 0.0])),
576            [2.0, 0.0, 0.0].into()
577        );
578        assert_relative_eq!(
579            sphere.support_mapping(&Cartesian::from([-1.0, 0.0, 0.0])),
580            [-2.0, 0.0, 0.0].into()
581        );
582
583        assert_relative_eq!(
584            sphere.support_mapping(&Cartesian::from([1.0, 1.0, 1.0])),
585            [1.1547005383792517, 1.1547005383792517, 1.1547005383792517].into()
586        );
587
588        Ok(())
589    }
590
591    #[test]
592    fn intersects_at() -> anyhow::Result<()> {
593        let sphere0 = Sphere::with_radius(2.0.try_into()?);
594        let sphere1 = Sphere::with_radius(4.0.try_into()?);
595        let identity = Versor::default();
596
597        check!(sphere0.intersects_at(&sphere1, &[0.0, 0.0, 5.9].into(), &identity));
598        check!(sphere0.intersects_at(&sphere1, &[0.0, 5.9, 0.0].into(), &identity));
599        check!(sphere0.intersects_at(&sphere1, &[5.9, 0.0, 0.0].into(), &identity));
600        check!(sphere0.intersects_at(&sphere1, &[3.4, 3.4, 3.4].into(), &identity));
601
602        check!(!sphere0.intersects_at(&sphere1, &[0.0, 0.0, 6.1].into(), &identity));
603        check!(!sphere0.intersects_at(&sphere1, &[0.0, 6.1, 0.0].into(), &identity));
604        check!(!sphere0.intersects_at(&sphere1, &[6.1, 0.0, 0.0].into(), &identity));
605        check!(!sphere0.intersects_at(&sphere1, &[3.52, 3.52, 3.52].into(), &identity));
606
607        let sphere0 = Convex(sphere0);
608        let sphere1 = Convex(sphere1);
609
610        check!(sphere0.intersects_at(&sphere1, &[0.0, 0.0, 5.9].into(), &identity));
611        check!(sphere0.intersects_at(&sphere1, &[0.0, 5.9, 0.0].into(), &identity));
612        check!(sphere0.intersects_at(&sphere1, &[5.9, 0.0, 0.0].into(), &identity));
613        check!(sphere0.intersects_at(&sphere1, &[3.4, 3.4, 3.4].into(), &identity));
614
615        check!(!sphere0.intersects_at(&sphere1, &[0.0, 0.0, 6.1].into(), &identity));
616        check!(!sphere0.intersects_at(&sphere1, &[0.0, 6.1, 0.0].into(), &identity));
617        check!(!sphere0.intersects_at(&sphere1, &[6.1, 0.0, 0.0].into(), &identity));
618        check!(!sphere0.intersects_at(&sphere1, &[3.52, 3.52, 3.52].into(), &identity));
619
620        Ok(())
621    }
622
623    #[test]
624    fn is_point_inside() -> anyhow::Result<()> {
625        let circle = Circle::with_radius(2.0.try_into()?);
626
627        check!(circle.is_point_inside(&Cartesian::from([0.0, 0.0])));
628        check!(circle.is_point_inside(&Cartesian::from([0.0, 1.0])));
629        check!(circle.is_point_inside(&Cartesian::from([0.0, -1.0])));
630        check!(circle.is_point_inside(&Cartesian::from([1.0, 0.0])));
631        check!(circle.is_point_inside(&Cartesian::from([-1.0, 0.0])));
632        check!(circle.is_point_inside(&Cartesian::from([2.0_f64.next_down(), 0.0])));
633        check!(circle.is_point_inside(&Cartesian::from([0.0, 2.0_f64.next_down()])));
634
635        check!(!circle.is_point_inside(&Cartesian::from([2.0, 0.0])));
636        check!(!circle.is_point_inside(&Cartesian::from([0.0, 2.0])));
637        check!(!circle.is_point_inside(&Cartesian::from([1.5, 1.5])));
638
639        Ok(())
640    }
641
642    #[test]
643    fn distribution() -> anyhow::Result<()> {
644        let circle = Circle::with_radius(4.0.try_into()?);
645        let mut rng = StdRng::seed_from_u64(4);
646
647        let points: Vec<_> = (&circle).sample_iter(&mut rng).take(N).collect();
648        check!(&points.iter().all(|p| circle.is_point_inside(p)));
649        check!(&points.iter().any(|p| p.dot(p) > 3.9));
650
651        Ok(())
652    }
653
654    #[test]
655    fn test_scale_length() -> anyhow::Result<()> {
656        let circle = Circle::with_radius(4.0.try_into()?);
657
658        let scaled_circle = circle.scale_length(2.0.try_into()?);
659        check!(scaled_circle.radius.get() == 8.0);
660
661        let scaled_circle = circle.scale_length(0.5.try_into()?);
662        check!(scaled_circle.radius.get() == 2.0);
663
664        Ok(())
665    }
666
667    #[test]
668    fn test_scale_volume() -> anyhow::Result<()> {
669        let circle = Circle::with_radius(4.0.try_into()?);
670
671        let scaled_circle = circle.scale_volume(4.0.try_into()?);
672        check!(scaled_circle.radius.get() == 8.0);
673
674        let scaled_circle = circle.scale_volume(0.25.try_into()?);
675        check!(scaled_circle.radius.get() == 2.0);
676
677        Ok(())
678    }
679
680    #[test]
681    fn test_map_basic() -> anyhow::Result<()> {
682        let circle_a = Circle::with_radius(4.0.try_into()?);
683        let circle_b = Circle::with_radius(8.0.try_into()?);
684
685        check!(
686            circle_a.map_point(Cartesian::from([0.0, 0.0]), &circle_b)
687                == Ok(Cartesian::from([0.0, 0.0]))
688        );
689        check!(
690            circle_b.map_point(Cartesian::from([0.0, 0.0]), &circle_a)
691                == Ok(Cartesian::from([0.0, 0.0]))
692        );
693
694        check!(
695            circle_a.map_point(Cartesian::from([100.0, 0.0]), &circle_b)
696                == Err(Error::PointOutsideShape)
697        );
698        check!(
699            circle_b.map_point(Cartesian::from([0.0, -200.0]), &circle_a)
700                == Err(Error::PointOutsideShape)
701        );
702
703        check!(
704            circle_a.map_point(Cartesian::from([3.0, 0.0]), &circle_b)
705                == Ok(Cartesian::from([6.0, 0.0]))
706        );
707        check!(
708            circle_b.map_point(Cartesian::from([-6.0, 0.0]), &circle_a)
709                == Ok(Cartesian::from([-3.0, 0.0]))
710        );
711
712        check!(
713            circle_a.map_point(Cartesian::from([-1.0, 2.0]), &circle_b)
714                == Ok(Cartesian::from([-2.0, 4.0]))
715        );
716        check!(
717            circle_b.map_point(Cartesian::from([2.0, -4.0]), &circle_a)
718                == Ok(Cartesian::from([1.0, -2.0]))
719        );
720
721        Ok(())
722    }
723
724    #[test]
725    fn test_map_surface() -> anyhow::Result<()> {
726        let mut rng = StdRng::seed_from_u64(3);
727        let uniform_radius = Uniform::new(1.0, 1000.0)?;
728        let uniform_angle = Uniform::new(0.0, 2.0 * PI)?;
729
730        for _ in 0..16384 {
731            let a = uniform_radius.sample(&mut rng);
732            let b = uniform_radius.sample(&mut rng);
733            let circle_a = Circle::with_radius(a.try_into()?);
734            let circle_b = Circle::with_radius(b.try_into()?);
735
736            // Test that points right on the boundary of one shape remain inside
737            // the other shape. If not implemented correctly, map_point might
738            // round  and place a point just outside the shape. This test fails
739            // without corner case handling in `map_point`.
740
741            let left = circle_a.map_point(Cartesian::from([(-a).next_up(), 0.0]), &circle_b)?;
742            check!(
743                circle_b.is_point_inside(&left),
744                "{left:?} should be inside {circle_b:?}"
745            );
746
747            let right = circle_a.map_point(Cartesian::from([a.next_down(), 0.0]), &circle_b)?;
748            check!(
749                circle_b.is_point_inside(&right),
750                "{right:?} should be inside {circle_b:?}"
751            );
752
753            let bottom = circle_a.map_point(Cartesian::from([0.0, (-a).next_up()]), &circle_b)?;
754            check!(
755                circle_b.is_point_inside(&bottom),
756                "{bottom:?} should be inside {circle_b:?}"
757            );
758
759            let top = circle_a.map_point(Cartesian::from([0.0, a.next_down()]), &circle_b)?;
760            check!(
761                circle_b.is_point_inside(&top),
762                "{top:?} should be inside {circle_b:?}"
763            );
764
765            for _ in 0..32 {
766                let theta = uniform_angle.sample(&mut rng);
767                let point = Cartesian::from([a * theta.cos(), b * theta.sin()]);
768
769                // `point` may be rounded in such a way that it falls outside the shape.
770                // Skip these cases. 32 iterations is sufficient to produce many interior
771                // points.
772                if !circle_a.is_point_inside(&point) {
773                    continue;
774                }
775
776                let mapped_point = circle_a.map_point(point, &circle_b)?;
777                check!(
778                    circle_b.is_point_inside(&mapped_point),
779                    "{mapped_point:?} should be inside {circle_b:?}"
780                );
781            }
782        }
783
784        Ok(())
785    }
786}
hoomd_geometry/shape/sphere.rs

hoomd_geometry/shape/
sphere.rs
