hoomd_geometry/shape/
twelvetwelve.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 [`TwelveTwelve`]
5
6use serde::{Deserialize, Serialize};
7
8use crate::IsPointInside;
9use hoomd_manifold::{Hyperbolic, Minkowski};
10use hoomd_vector::Metric;
11use std::f64::consts::PI;
12
13/// A regular dodecagon in two-dimensional hyperbolic space.
14///
15/// [`TwelveTwelve`] implements a single regular dodecagon in the {12,12} tiling
16/// of two-dimensional hyperbolic space. The scaling of the octagon is set such
17/// that each of the angles is $` \frac{2\pi}{12} `$ so that twelve equivalent
18/// dodecagons meet at each vertex.
19#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
20pub struct TwelveTwelve {}
21
22impl IsPointInside<Hyperbolic<3>> for TwelveTwelve {
23    /// Checks if a given Hyperbolic point is inside [`TwelveTwelve`].
24    ///
25    /// # Example
26    /// ```
27    /// use hoomd_geometry::{IsPointInside, shape::TwelveTwelve};
28    /// use hoomd_manifold::Hyperbolic;
29    /// use std::f64::consts::PI;
30    ///
31    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
32    /// let twelve_twelve = TwelveTwelve {};
33    ///
34    /// let point = Hyperbolic::<3>::from_polar_coordinates(1.0, PI / 8.0);
35    /// assert!(twelve_twelve.is_point_inside(&point));
36    /// # Ok(())
37    /// # }
38    /// ```
39    #[inline]
40    fn is_point_inside(&self, point: &Hyperbolic<3>) -> bool {
41        TwelveTwelve::distance_to_boundary(point) >= 0.0
42    }
43}
44
45impl TwelveTwelve {
46    /// Computes the shortest distance between a given point and the boundary
47    /// of `TwelveTwelve`.
48    ///
49    /// The shortest distance is computed by finding the arclength of the geodesic
50    /// which passes through the query point and intersects the boundary at a
51    /// right angle.
52    ///
53    /// # Example
54    /// ```
55    /// use approxim::assert_relative_eq;
56    /// use hoomd_geometry::shape::TwelveTwelve;
57    /// use hoomd_manifold::{Hyperbolic, Minkowski};
58    /// use std::f64::consts::PI;
59    ///
60    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
61    /// let v: f64 = TwelveTwelve::CUSP_TO_EDGE - 0.4;
62    /// let theta: f64 = PI / 12.0;
63    /// let x = Hyperbolic::from_minkowski_coordinates(
64    ///     [
65    ///         (v.sinh()) * (theta.cos()),
66    ///         (v.sinh()) * (theta.sin()),
67    ///         (v.cosh()),
68    ///     ]
69    ///     .into(),
70    /// );
71    /// assert_relative_eq!(
72    ///     TwelveTwelve::distance_to_boundary(&x),
73    ///     0.4,
74    ///     epsilon = 1e-12
75    /// );
76    /// # Ok(())
77    /// # }
78    /// ```
79    #[inline]
80    #[must_use]
81    pub fn distance_to_boundary(point: &Hyperbolic<3>) -> f64 {
82        let theta =
83            (point.coordinates()[1].atan2(point.coordinates()[0])).rem_euclid(PI / 6.0) - PI / 12.0;
84        let boost = (point.coordinates()[2]).acosh();
85        let (b_sinh, b_cosh) = (boost.sinh(), point.coordinates()[2]);
86        let xi = Self::CUSP_TO_EDGE;
87        let (xi_sinh, xi_cosh) = (xi.sinh(), xi.cosh());
88        // boost into frame where edge is the vertical diameter
89        let edge_as_diameter: Hyperbolic<3> =
90            Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([
91                xi_cosh * b_sinh * (theta.cos()) - xi_sinh * b_cosh,
92                b_sinh * (theta.sin()),
93                -xi_sinh * b_sinh * (theta.cos()) + xi_cosh * b_cosh,
94            ]));
95        let flipped = Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([
96            -edge_as_diameter.coordinates()[0],
97            edge_as_diameter.coordinates()[1],
98            edge_as_diameter.coordinates()[2],
99        ]));
100        let sign = -(edge_as_diameter.coordinates()[0]).signum();
101        sign * (edge_as_diameter.distance(&flipped)) / 2.0
102    }
103    /// Apply a lattice transformation to a point.
104    #[inline]
105    #[must_use]
106    pub fn gamma(eta: f64, theta: f64, point: &[f64; 3]) -> [f64; 3] {
107        let (eta_sinh_squared, two_eta_sinh, theta_sin, theta_cos) = (
108            (eta.sinh()).powi(2),
109            (2.0 * eta).sinh(),
110            theta.sin(),
111            theta.cos(),
112        );
113        [
114            (2.0 * (eta_sinh_squared) * ((theta_cos).powi(2)) + 1.0) * point[0]
115                + ((2.0 * theta).sin()) * (eta_sinh_squared) * point[1]
116                + (two_eta_sinh) * (theta_cos) * point[2],
117            ((2.0 * theta).sin()) * (eta_sinh_squared) * point[0]
118                + (2.0 * (eta_sinh_squared) * ((theta_sin).powi(2)) + 1.0) * point[1]
119                + (two_eta_sinh) * (theta_sin) * point[2],
120            (two_eta_sinh) * (theta_cos) * point[0]
121                + (two_eta_sinh) * (theta_sin) * point[1]
122                + ((2.0 * eta).cosh()) * point[2],
123        ]
124    }
125    /// Calculate the change in angle in the tangent bundle associated with a lattice
126    /// transformation.
127    #[inline]
128    #[must_use]
129    pub fn reorient(theta: f64, point: &[f64; 3]) -> f64 {
130        let (q_u, q_v) = (point[0] / (1.0 + point[2]), point[1] / (1.0 + point[2]));
131        let alpha = (1.0 + (PI / 6.0).cos()).sqrt();
132        let beta = (theta.cos()) * ((2.0 * ((PI / 6.0).cos())).sqrt());
133        let gamma = (theta.sin()) * ((2.0 * ((PI / 6.0).cos())).sqrt());
134        let p_x = alpha + beta * q_u + gamma * q_v;
135        let p_y = beta * q_v - gamma * q_u;
136        -2.0 * (p_y.atan2(p_x))
137    }
138    /// Cusp-to-vertex distance for the {12,12} tiling for Gauss curvature K = -1.
139    pub const TWELVETWELVE: f64 = 3.325_771_782_117_242;
140    /// Cusp-to-middle-of-edge distance for the {12,12} tiling for Gauss curvature K = -1.
141    pub const CUSP_TO_EDGE: f64 = 1.991_652_391_049_437;
142    /// Length of one of the sides of the {12,12} tiling for Gauss curvature K = -1.
143    pub const EDGE_LENGTH: f64 = 3.983_304_782_098_874;
144}
145
146#[cfg(test)]
147mod tests {
148    use super::*;
149    use approxim::assert_relative_eq;
150    use hoomd_manifold::{Hyperbolic, HyperbolicDisk, Minkowski};
151    use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
152    use std::ops::Not;
153
154    #[test]
155    fn boundary_distance() {
156        // Distance to the edge of the {12,12} fundamental domain
157        let e = Hyperbolic::<3>::from_polar_coordinates(1.0, 0.1);
158        let e_edge_distance = TwelveTwelve::distance_to_boundary(&e);
159        let e_edge_distance_numeric = 1.028_489_242_583_61;
160        assert_relative_eq!(e_edge_distance, e_edge_distance_numeric, epsilon = 1e-12);
161
162        let f = Hyperbolic::<3>::from_polar_coordinates(0.6, 0.2 + PI / 6.0);
163        let f_edge_distance = TwelveTwelve::distance_to_boundary(&f);
164        let f_edge_distance_numeric = 1.393_774_804_611_314;
165        assert_relative_eq!(f_edge_distance, f_edge_distance_numeric, epsilon = 1e-12);
166    }
167
168    #[test]
169    fn inside_is_inside() {
170        let twelve_twelve = TwelveTwelve {};
171        let r = TwelveTwelve::CUSP_TO_EDGE;
172        let mut rng = StdRng::seed_from_u64(239);
173        let disk = HyperbolicDisk {
174            disk_radius: r.try_into().expect("hard-coded positive number"),
175            point: Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([0.0, 0.0, 1.0])),
176        };
177        let random_point: Hyperbolic<3> = disk.sample(&mut rng);
178        assert!(twelve_twelve.is_point_inside(&random_point));
179
180        let point_1 = Hyperbolic::<3>::from_polar_coordinates(1.99, PI / 12.0);
181        assert!(twelve_twelve.is_point_inside(&point_1));
182
183        let point_2 = Hyperbolic::<3>::from_polar_coordinates(3.32, PI / 6.0);
184        assert!(twelve_twelve.is_point_inside(&point_2));
185    }
186
187    #[test]
188    fn outside_is_outside() {
189        let twelve_twelve = TwelveTwelve {};
190        let point_1 = Hyperbolic::<3>::from_polar_coordinates(2.0, PI / 12.0);
191        assert!((twelve_twelve.is_point_inside(&point_1)).not());
192
193        let point_2 = Hyperbolic::<3>::from_polar_coordinates(3.33, PI / 6.0);
194        assert!((twelve_twelve.is_point_inside(&point_2)).not());
195    }
196}