hoomd_geometry/shape/
hyperbolic_convex_polytope.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 Overlap Check for Hyperbolic Surfaces
5
6use crate::{BoundingSphereRadius, IntersectsAtGlobal};
7use hoomd_manifold::{Hyperbolic, Minkowski};
8use hoomd_utility::valid::PositiveReal;
9use hoomd_vector::{Angle, Metric};
10use robust::{Coord, orient2d};
11use std::f64::consts::PI;
12
13/// A faceted hyperbolic solid defined by the convex hull of its vertices.
14///
15/// # Examples
16///
17/// Construction and basic methods:
18/// ```
19/// use hoomd_geometry::{
20///     BoundingSphereRadius, shape::HyperbolicConvexPolytope,
21/// };
22///
23/// # fn main() -> Result<(), hoomd_geometry::Error> {
24/// let hyperbolic_square = HyperbolicConvexPolytope::<3>::regular(4, 0.5);
25///
26/// let bounding_radius = hyperbolic_square.bounding_sphere_radius();
27///
28/// assert_eq!(bounding_radius.get(), 0.5);
29/// # Ok(())
30/// # }
31/// ```
32///
33/// Overlap check:
34/// ```
35/// use hoomd_geometry::{IntersectsAtGlobal, shape::HyperbolicConvexPolytope};
36/// use hoomd_manifold::Hyperbolic;
37/// use hoomd_vector::Angle;
38///
39/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
40/// let square = HyperbolicConvexPolytope::<3>::regular(4, 0.5);
41/// assert!(square.intersects_at_global(
42///     &square,
43///     &Hyperbolic::<3>::default(),
44///     &Angle::default(),
45///     &Hyperbolic::<3>::from_polar_coordinates(0.49, 2.3),
46///     &Angle::from(0.4)
47/// ));
48///
49/// assert!(!square.intersects_at_global(
50///     &square,
51///     &Hyperbolic::<3>::default(),
52///     &Angle::default(),
53///     &Hyperbolic::<3>::from_polar_coordinates(3.0, 2.3),
54///     &Angle::from(0.2)
55/// ));
56/// # Ok(())
57/// # }
58/// ```
59#[derive(Clone, Debug, PartialEq)]
60pub struct HyperbolicConvexPolytope<const N: usize> {
61    /// The vertices of the shape
62    vertices: Vec<Hyperbolic<N>>,
63    /// The radius of the bounding sphere of the shape in the Hyperbolic metric.
64    bounding_radius: f64,
65}
66
67impl<const N: usize> HyperbolicConvexPolytope<N> {
68    /// Get the vertices of the shape.
69    #[inline]
70    #[must_use]
71    pub fn vertices(&self) -> &[Hyperbolic<N>] {
72        &self.vertices
73    }
74}
75
76impl<const N: usize> BoundingSphereRadius for HyperbolicConvexPolytope<N> {
77    #[inline]
78    fn bounding_sphere_radius(&self) -> PositiveReal {
79        self.bounding_radius
80            .try_into()
81            .expect("hard coded constant should be positive")
82    }
83}
84
85/// A two-dimensional hyperbolic convex polytope.
86/// ```rust
87/// use hoomd_geometry::shape::HyperbolicConvexPolygon;
88///
89/// # fn main() -> Result<(), hoomd_geometry::Error> {
90/// let hyperbolic_pentagon = HyperbolicConvexPolygon::regular(5, 1.0_f64);
91/// # Ok(())
92/// # }
93/// ```
94pub type HyperbolicConvexPolygon = HyperbolicConvexPolytope<3>;
95
96impl HyperbolicConvexPolytope<3> {
97    /// Create a regular *n*-gon with *n* vertices and a given circumradius in
98    /// hyperbolic space.
99    #[inline]
100    #[must_use]
101    pub fn regular(n: usize, circumradius: f64) -> HyperbolicConvexPolytope<3> {
102        HyperbolicConvexPolytope {
103            vertices: (0..n)
104                .map(|x| {
105                    let theta = 2.0 * PI * (x as f64) / (n as f64);
106                    Hyperbolic::<3>::from_polar_coordinates(circumradius, theta)
107                })
108                .collect::<Vec<_>>(),
109            bounding_radius: circumradius,
110        }
111    }
112    /// Calculate the distance between the center of a `HyperbolicConvexPolytope<3>`
113    /// centered at the origin and the edge at angle phi. The calculation works by
114    /// computing the intersection of the the hyperboloid with a plane passing
115    /// through the origin and the two adjacent vertices of the polygon.
116    ///
117    /// # Example
118    /// ```
119    /// use approxim::assert_relative_eq;
120    /// use hoomd_geometry::shape::HyperbolicConvexPolytope;
121    /// use std::f64::consts::PI;
122    ///
123    /// let bounding_radius = 0.5;
124    /// let hyperbolic_square =
125    ///     HyperbolicConvexPolytope::<3>::regular(4, bounding_radius);
126    ///
127    /// // calculation using hyperbolic trigonometry
128    /// let pi_fourths_distance =
129    ///     (((PI / 4.0).cos()) * (bounding_radius.tanh())).atanh();
130    ///
131    /// let edge_distance = hyperbolic_square.edge_distance(PI / 4.0);
132    ///
133    /// assert_relative_eq!(pi_fourths_distance, edge_distance);
134    /// ```
135    #[inline]
136    #[must_use]
137    pub fn edge_distance(&self, phi: f64) -> f64 {
138        let n = self.vertices.len() as f64;
139        let phi_mod = phi.rem_euclid(2.0 * PI / n) - PI / n;
140        let eta_tanh = self.bounding_radius.tanh();
141        let arg = (eta_tanh * ((2.0 * PI / n).sin()))
142            / ((PI / n - phi_mod).sin() + (PI / n + phi_mod).sin());
143        arg.atanh()
144    }
145    #[inline]
146    /// Transform `points` to the frame where the `vertex_num`-th vertex of an
147    /// oriented hyperbolic polygon with with bounding radius `bounding_radius`
148    /// is at the origin.
149    fn to_vertex_frame_oriented(
150        body_position: &Hyperbolic<3>,
151        body_orientation: Angle,
152        vertex_num: usize,
153        bounding_radius: f64,
154        points: &[Hyperbolic<3>],
155        num_of_sides: usize,
156    ) -> Vec<Hyperbolic<3>> {
157        let theta = body_position.coordinates()[1].atan2(body_position.coordinates()[0]);
158        let eta = (body_position.coordinates()[2]).acosh();
159        let tau_over_two = -eta / 2.0;
160        let poincare = body_position.to_poincare();
161        let body_angle_body = -2.0
162            * (((tau_over_two.sinh())
163                * ((theta.cos()) * poincare[1] - (theta.sin()) * poincare[0]))
164                .atan2(
165                    (tau_over_two.cosh())
166                        + (tau_over_two.sinh())
167                            * ((theta.cos()) * poincare[0] + (theta.sin()) * poincare[1]),
168                ))
169            + body_orientation.theta;
170
171        let alpha =
172            body_angle_body - theta + 2.0 * (vertex_num as f64) * PI / (num_of_sides as f64);
173        let r = bounding_radius;
174        let vertex_translate = |point: &Hyperbolic<3>| -> Hyperbolic<3> {
175            let pt = point.point().coordinates;
176            let (eta_sinh, eta_cosh, r_sinh, r_cosh) = (eta.sinh(), eta.cosh(), r.sinh(), r.cosh());
177            let (alpha_cos, alpha_sin, theta_cos, theta_sin) =
178                (alpha.cos(), alpha.sin(), theta.cos(), theta.sin());
179            let translated = Minkowski::from([
180                (eta_cosh * r_cosh * alpha_cos * theta_cos - r_cosh * alpha_sin * theta_sin
181                    + eta_sinh * r_sinh * theta_cos)
182                    * pt[0]
183                    + (eta_cosh * r_cosh * alpha_cos * theta_sin
184                        + r_cosh * alpha_sin * theta_cos
185                        + eta_sinh * r_sinh * theta_sin)
186                        * pt[1]
187                    - (eta_sinh * r_cosh * alpha_cos + eta_cosh * r_sinh) * pt[2],
188                -(eta_cosh * alpha_sin * theta_cos + alpha_cos * theta_sin) * pt[0]
189                    + (-eta_cosh * alpha_sin * theta_sin + alpha_cos * theta_cos) * pt[1]
190                    + eta_sinh * alpha_sin * pt[2],
191                (-eta_cosh * r_sinh * alpha_cos * theta_cos + r_sinh * alpha_sin * theta_sin
192                    - eta_sinh * r_cosh * theta_cos)
193                    * pt[0]
194                    - (eta_cosh * r_sinh * alpha_cos * theta_sin
195                        + r_sinh * alpha_sin * theta_cos
196                        + eta_sinh * r_cosh * theta_sin)
197                        * pt[1]
198                    + (eta_sinh * r_sinh * alpha_cos + eta_cosh * r_cosh) * pt[2],
199            ]);
200            Hyperbolic::from_minkowski_coordinates(translated)
201        };
202        points.iter().map(vertex_translate).collect::<Vec<_>>()
203    }
204    /// Transform a vertex of a hyperbolic polygon into the system frame.
205    #[inline]
206    fn vertex_to_system_frame(
207        vertex: &Hyperbolic<3>,
208        body_orientation: Angle,
209        body_position: &Hyperbolic<3>,
210    ) -> Hyperbolic<3> {
211        let theta = body_position.coordinates()[1].atan2(body_position.coordinates()[0]);
212        let nu = (body_position.coordinates()[2]).acosh();
213        let tau_over_two = -nu / 2.0;
214        let poincare = body_position.to_poincare();
215        let body_angle_body = (-2.0
216            * (((tau_over_two.sinh())
217                * ((theta.cos()) * poincare[1] - (theta.sin()) * poincare[0]))
218                .atan2(
219                    (tau_over_two.cosh())
220                        + (tau_over_two.sinh())
221                            * ((theta.cos()) * poincare[0] + (theta.sin()) * poincare[1]),
222                ))
223            + body_orientation.theta)
224            .rem_euclid(2.0 * PI);
225        let phi = body_angle_body - theta;
226        let pt = vertex.point().coordinates;
227        let (nu_sinh, nu_cosh) = (nu.sinh(), nu.cosh());
228        let (theta_sin, theta_cos, phi_sin, phi_cos) =
229            (theta.sin(), theta.cos(), phi.sin(), phi.cos());
230        let transformed = Minkowski::from([
231            (nu_cosh * phi_cos * theta_cos - phi_sin * theta_sin) * pt[0]
232                - (nu_cosh * phi_sin * theta_cos + phi_cos * theta_sin) * pt[1]
233                + nu_sinh * theta_cos * pt[2],
234            (nu_cosh * phi_cos * theta_sin + phi_sin * theta_cos) * pt[0]
235                + (-nu_cosh * phi_sin * theta_sin + phi_cos * theta_cos) * pt[1]
236                + nu_sinh * theta_sin * pt[2],
237            nu_sinh * phi_cos * pt[0] - nu_sinh * phi_sin * pt[1] + nu_cosh * pt[2],
238        ]);
239        Hyperbolic::from_minkowski_coordinates(transformed)
240    }
241}
242
243impl IntersectsAtGlobal<HyperbolicConvexPolytope<3>, Hyperbolic<3>, Angle>
244    for HyperbolicConvexPolytope<3>
245{
246    #[inline]
247    #[allow(clippy::too_many_lines, reason = "complicated function")]
248    fn intersects_at_global(
249        &self,
250        other: &HyperbolicConvexPolytope<3>,
251        r_self: &Hyperbolic<3>,
252        o_self: &Angle,
253        r_other: &Hyperbolic<3>,
254        o_other: &Angle,
255    ) -> bool {
256        let d = r_self.distance(r_other);
257        if d > 2.0 * self.bounding_radius {
258            return false;
259        }
260        let mut result = true;
261        let mut v_count = 0_usize;
262        let n_self = self.vertices.len();
263        let n_other = other.vertices.len();
264        while result && (v_count < n_self + n_other) {
265            if v_count < n_self {
266                let v_num = (v_count) % n_self;
267                let v_next = (v_num + 1) % n_self;
268                let v_next_next = (v_num + 2) % n_self;
269                // translate all vertices
270                // need to do this for every other vertex
271                let v_1 = Self::vertex_to_system_frame(&self.vertices[v_num], *o_self, r_self);
272                let v_2 = Self::vertex_to_system_frame(&self.vertices[v_next], *o_self, r_self);
273                let v_3 =
274                    Self::vertex_to_system_frame(&self.vertices[v_next_next], *o_self, r_self);
275                let other_vertices = self
276                    .vertices
277                    .iter()
278                    .map(|vertex| Self::vertex_to_system_frame(vertex, *o_other, r_other))
279                    .collect::<Vec<Hyperbolic<3>>>();
280                let self_translated = Self::to_vertex_frame_oriented(
281                    r_self,
282                    *o_self,
283                    v_next,
284                    self.bounding_radius,
285                    &[v_1, v_2, v_3],
286                    n_self,
287                );
288                let other_translated = Self::to_vertex_frame_oriented(
289                    r_self,
290                    *o_self,
291                    v_next,
292                    self.bounding_radius,
293                    &other_vertices,
294                    n_self,
295                );
296                // convert to poincare coordinates to perform orientation checks.
297                let self_coord = self_translated
298                    .iter()
299                    .map(|pt: &Hyperbolic<3>| {
300                        let poincare = pt.to_poincare();
301                        Coord {
302                            x: poincare[0],
303                            y: poincare[1],
304                        }
305                    })
306                    .collect::<Vec<Coord<f64>>>();
307                let other_coord = other_translated
308                    .iter()
309                    .map(|pt: &Hyperbolic<3>| {
310                        let poincare = pt.to_poincare();
311                        Coord {
312                            x: poincare[0],
313                            y: poincare[1],
314                        }
315                    })
316                    .collect::<Vec<Coord<f64>>>();
317                // then do edge check
318                let mut overlap = false;
319                let mut counter = 0_usize;
320                while !overlap && (counter < other_vertices.len()) {
321                    if orient2d(self_coord[0], self_coord[1], other_coord[counter]) >= 0.0 {
322                        // break out of loop once one of the vertices is on the wrong side of the line
323                        overlap = true;
324                        break;
325                    }
326                    counter += 1;
327                }
328                counter = 0_usize;
329                // if overlap is false, then no need to check next edge FIX THIS!
330                if overlap {
331                    overlap = false;
332                    while !overlap && (counter < other_vertices.len()) {
333                        if orient2d(self_coord[1], self_coord[2], other_coord[counter]) >= 0.0 {
334                            // break out of loop once one of the vertices is on the wrong side of the line
335                            overlap = true;
336                            break;
337                        }
338                        counter += 1;
339                    }
340                }
341                result = overlap;
342                v_count += 2;
343            } else {
344                let v_num = (v_count) % n_other;
345                let v_next = (v_num + 1) % n_other;
346                let v_next_next = (v_num + 2) % n_other;
347                // translate all vertices
348                let v_1 = Self::vertex_to_system_frame(&other.vertices[v_num], *o_other, r_other);
349                let v_2 = Self::vertex_to_system_frame(&other.vertices[v_next], *o_other, r_other);
350                let v_3 =
351                    Self::vertex_to_system_frame(&other.vertices[v_next_next], *o_other, r_other);
352                let other_vertices = self
353                    .vertices
354                    .iter()
355                    .map(|vertex| Self::vertex_to_system_frame(vertex, *o_self, r_self))
356                    .collect::<Vec<Hyperbolic<3>>>();
357                let self_translated = Self::to_vertex_frame_oriented(
358                    r_other,
359                    *o_other,
360                    v_next,
361                    other.bounding_radius,
362                    &[v_1, v_2, v_3],
363                    n_other,
364                );
365                let other_translated = Self::to_vertex_frame_oriented(
366                    r_other,
367                    *o_other,
368                    v_next,
369                    other.bounding_radius,
370                    &other_vertices,
371                    n_other,
372                );
373                // convert to poincare coordinates to perform orientation checks.
374                let self_coord = self_translated
375                    .iter()
376                    .map(|pt: &Hyperbolic<3>| {
377                        let poincare = pt.to_poincare();
378                        Coord {
379                            x: poincare[0],
380                            y: poincare[1],
381                        }
382                    })
383                    .collect::<Vec<Coord<f64>>>();
384                let other_coord = other_translated
385                    .iter()
386                    .map(|pt: &Hyperbolic<3>| {
387                        let poincare = pt.to_poincare();
388                        Coord {
389                            x: poincare[0],
390                            y: poincare[1],
391                        }
392                    })
393                    .collect::<Vec<Coord<f64>>>();
394                // then do edge check
395                let mut overlap = false;
396                let mut counter = 0_usize;
397                while !overlap && (counter < other_vertices.len()) {
398                    if orient2d(self_coord[0], self_coord[1], other_coord[counter]) >= 0.0 {
399                        // break out of loop once one of the vertices is on the wrong side of the line
400                        overlap = true;
401                        break;
402                    }
403                    counter += 1;
404                }
405                counter = 0_usize;
406                if overlap {
407                    overlap = false;
408                    while !overlap && (counter < other_vertices.len()) {
409                        if orient2d(self_coord[1], self_coord[2], other_coord[counter]) >= 0.0 {
410                            // break out of loop once one of the vertices is on the wrong side of the line
411                            overlap = true;
412                        }
413                        counter += 1;
414                    }
415                }
416                result = overlap;
417                v_count += 2;
418            }
419        }
420        result
421    }
422}
423
424#[cfg(test)]
425mod tests {
426    use super::*;
427    use crate::shape::EightEight;
428    use approxim::assert_relative_eq;
429    use hoomd_manifold::{Hyperbolic, Minkowski};
430    use hoomd_vector::Angle;
431
432    #[test]
433    fn octagon_edges() {
434        let center_dist = 1.528_570_919_480_998;
435        let quarter_dist = 1.643_866_837_922_488;
436        let octagon = HyperbolicConvexPolytope::<3>::regular(8, EightEight::EIGHTEIGHT);
437        assert_relative_eq!(
438            center_dist,
439            octagon.edge_distance(-3.0 * PI / 8.0),
440            epsilon = 1e-12
441        );
442        assert_relative_eq!(
443            quarter_dist,
444            octagon.edge_distance(PI / 16.0),
445            epsilon = 1e-12
446        );
447        assert_relative_eq!(
448            EightEight::EIGHTEIGHT,
449            octagon.edge_distance(PI / 4.0),
450            epsilon = 1e-12
451        );
452    }
453    #[test]
454    fn square_edges() {
455        let center_dist = 0.602_080_559_268_716;
456        let quarter_dist = 0.666_842_324_123_307;
457        let square = HyperbolicConvexPolytope::<3>::regular(4, 1.0);
458        assert_relative_eq!(center_dist, square.edge_distance(PI / 4.0), epsilon = 1e-12);
459        assert_relative_eq!(
460            quarter_dist,
461            square.edge_distance(-PI / 8.0),
462            epsilon = 1e-12
463        );
464        assert_relative_eq!(1.0, square.edge_distance(PI / 2.0), epsilon = 1e-12);
465    }
466    #[test]
467    fn center_at_oriented_vertex() {
468        let square = HyperbolicConvexPolytope::<3>::regular(4, 0.5);
469        let (boost, rotation, orientation) = (0.5, PI / 4.0, 0.4);
470        let body_position = Hyperbolic::<3>::from_polar_coordinates(boost, rotation);
471        let square_system = square
472            .vertices()
473            .iter()
474            .map(|v| {
475                HyperbolicConvexPolygon::vertex_to_system_frame(
476                    v,
477                    Angle::from(orientation),
478                    &body_position,
479                )
480            })
481            .collect::<Vec<Hyperbolic<3>>>();
482        let translated = HyperbolicConvexPolygon::to_vertex_frame_oriented(
483            &body_position,
484            Angle::from(orientation),
485            2_usize,
486            0.5,
487            &square_system,
488            4_usize,
489        );
490        assert_relative_eq!(0.0, translated[2].coordinates()[0], epsilon = 1e-12);
491        assert_relative_eq!(0.0, translated[2].coordinates()[1], epsilon = 1e-12);
492        assert_relative_eq!(1.0, translated[2].coordinates()[2], epsilon = 1e-12);
493    }
494    #[test]
495    fn no_square_overlap() {
496        let square = HyperbolicConvexPolytope::<3>::regular(4, 0.5);
497        let boost: f64 = 3.0;
498        let rotation: f64 = 2.3;
499        let orientation: f64 = 0.4;
500        let x_j = Hyperbolic::<3>::from_polar_coordinates(boost, rotation);
501        assert!(!square.intersects_at_global(
502            &square,
503            &Hyperbolic::<3>::default(),
504            &Angle::default(),
505            &x_j,
506            &Angle::from(orientation)
507        ));
508    }
509    #[test]
510    fn square_overlap() {
511        let square = HyperbolicConvexPolytope::<3>::regular(4, 0.5);
512        let boost: f64 = 0.49;
513        let rotation: f64 = 2.3;
514        let orientation: f64 = 0.4;
515        let x_j = Hyperbolic::<3>::from_polar_coordinates(boost, rotation);
516        assert!(square.intersects_at_global(
517            &square,
518            &Hyperbolic::<3>::default(),
519            &Angle::default(),
520            &x_j,
521            &Angle::from(orientation)
522        ));
523    }
524    #[test]
525    fn overlap_translation_check() {
526        let r_0 = 0.5;
527        let square = HyperbolicConvexPolytope::<3>::regular(4, r_0);
528        let com = Hyperbolic::<3>::from_polar_coordinates(1.0, 0.0);
529        let distance = 2.0;
530        let num_spaces: usize = 10;
531        let num_trials: usize = 15;
532        let spacing = 1.321_592_891_727_355;
533        let trials = (0..num_trials)
534            .map(|n| -(n as f64) * spacing / (num_spaces as f64))
535            .collect::<Vec<f64>>();
536
537        let nudged_centers = trials
538            .iter()
539            .map(|inch| {
540                let original_center = [
541                    (1.0_f64 + distance).sinh(),
542                    0.0_f64,
543                    (1.0_f64 + distance).cosh(),
544                ];
545                let translated = Minkowski::from([
546                    original_center[0] * (inch.cosh()) + original_center[2] * (inch.sinh()),
547                    original_center[1],
548                    original_center[0] * (inch.sinh()) + original_center[2] * (inch.cosh()),
549                ]);
550                Hyperbolic::from_minkowski_coordinates(translated)
551            })
552            .collect::<Vec<Hyperbolic<3>>>();
553        // Check over overlaps
554        for translated in nudged_centers.iter().take(num_spaces) {
555            assert!(!square.intersects_at_global(
556                &square,
557                &com,
558                &Angle::from(PI / 4.0),
559                translated,
560                &Angle::from(PI / 4.0)
561            ));
562        }
563        for translated in nudged_centers.iter().take(num_trials).skip(num_spaces) {
564            assert!(square.intersects_at_global(
565                &square,
566                &com,
567                &Angle::from(PI / 4.0),
568                translated,
569                &Angle::from(PI / 4.0)
570            ));
571        }
572    }
573    #[test]
574    fn overlap_rotation_check() {
575        let r_0 = 0.5;
576        let boost: f64 = 0.339_203_554_136_322;
577        let distance: f64 = 0.45;
578        let square = HyperbolicConvexPolytope::<3>::regular(4, r_0);
579        let num_spaces: usize = 10;
580        let num_trials: usize = 15;
581        let spacing = 0.365_106_058_818_114;
582        let trials = (0..num_trials)
583            .map(|n| (n as f64) * spacing / (num_spaces as f64))
584            .collect::<Vec<f64>>();
585        let center_1 = Hyperbolic::<3>::from_polar_coordinates(-boost, 0.0);
586        let center_2 = Hyperbolic::<3>::from_polar_coordinates(distance, 0.0);
587        // Check over overlaps
588        for ep in trials.iter().take(num_spaces) {
589            assert!(!square.intersects_at_global(
590                &square,
591                &center_1,
592                &Angle::from(PI / 4.0),
593                &center_2,
594                &Angle::from(ep + PI / 4.0)
595            ));
596        }
597        for ep in trials.iter().take(num_trials).skip(num_spaces) {
598            assert!(square.intersects_at_global(
599                &square,
600                &center_1,
601                &Angle::from(PI / 4.0),
602                &center_2,
603                &Angle::from(ep + PI / 4.0)
604            ));
605        }
606    }
607}
hoomd_geometry/shape/hyperbolic_convex_polytope.rs

hoomd_geometry/shape/
hyperbolic_convex_polytope.rs
