hoomd_geometry/shape/
simplex3.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//! A tetrahedron in three dimensions. This struct should be viewed as a prototype for
5//! more complex geometries in addition to its standalone functionality.
6use hoomd_utility::valid::PositiveReal;
7use itertools::Itertools;
8use serde::{Deserialize, Serialize};
9use std::{array, fmt};
10
11use hoomd_vector::{Cartesian, Cross, InnerProduct, Rotate, Rotation, RotationMatrix};
12
13use crate::{BoundingSphereRadius, IntersectsAt, IntersectsAtGlobal, SupportMapping, Volume};
14
15/// The hull of any 4 noncoplanar points in three dimensions.
16///
17/// # Examples
18///
19/// A [`Simplex3`], or equivalently a (nonuniform) tetrahedron, is the simplest faceted
20/// three-dimensional geometry. Because a simplex has an intrinsic idea of its position in
21/// space (defined by its barycenter), this struct provides implementations for spatial
22/// translation and a center of mass.
23///
24/// A uniform tetrahedron with edge length 2*sqrt(2), centered at the origin
25/// Vertices are defined by 3-length positive-signed permutations of [1,-1]
26///
27/// ```rust
28/// use hoomd_geometry::{IntersectsAt, Volume, shape::Simplex3};
29/// use hoomd_vector::{Cartesian, Rotation, Versor};
30///
31/// let tet = Simplex3::default();
32///
33/// assert_eq!(tet.vertices()[0], [1.0; 3].into());
34/// assert_eq!(tet.centroid(), [0.0; 3].into());
35///
36/// assert_eq!(tet.volume(), 8.0 / 3.0);
37///
38/// assert_eq!(
39///     tet.get_edge_vectors()[0],
40///     tet.vertices()[1] - tet.vertices()[0]
41/// );
42/// ```
43///
44/// [`Simplex3`] instances support a fast intersection detection algorithm:
45///
46/// ```rust
47/// # use hoomd_geometry::{shape::Simplex3, IntersectsAt};
48/// # use hoomd_vector::{Cartesian, Rotation, Versor};
49/// # let tet = Simplex3::default();
50/// let other_tet = Simplex3::default();
51/// let displacement = [2.0, 2.0, 0.0].into();
52/// let q_ij = Versor::identity();
53///
54/// assert!(tet.intersects_at(&other_tet, &displacement, &q_ij));
55///
56/// assert!(!tet.intersects_at(&other_tet, &[2.001, 2.0, 0.0].into(), &q_ij));
57/// ```
58///
59/// Although generally advised, Simplex3 tetrahedra are not required to be convex:
60/// ```rust
61/// # use hoomd_geometry::{shape::Simplex3, Volume};
62/// # use hoomd_vector::Cartesian;
63/// let planar_tetrahedron = Simplex3::from(
64///     [[0.0; 3], [0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
65///         .map(Cartesian::from),
66/// );
67///
68/// assert_eq!(planar_tetrahedron.volume(), 0.0);
69/// assert_eq!(planar_tetrahedron.centroid(), [1.0, 1.0, 0.0].into());
70/// ```
71#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
72pub struct Simplex3 {
73    /// Vertices of the simplex
74    vertices: [Cartesian<3>; 4], // NOT public, to force orientation on construction
75}
76
77impl SupportMapping<Cartesian<3>> for Simplex3 {
78    #[inline]
79    fn support_mapping(&self, n: &Cartesian<3>) -> Cartesian<3> {
80        let dots = self.vertices.map(|v| v.dot(n));
81        self.vertices[dots
82            .iter()
83            // position_max_by consumes ~9% of Xenocollide runtime!
84            .position_max_by(|x, y| x.partial_cmp(y).unwrap_or(std::cmp::Ordering::Equal))
85            .expect("dot product results should always be comparable")]
86    }
87}
88
89impl Volume for Simplex3 {
90    #[inline]
91    fn volume(&self) -> f64 {
92        let (ba, ca, da) = (
93            self.b() - self.a(),
94            self.c() - self.a(),
95            self.d() - self.a(),
96        );
97        1.0 / 6.0 * da.dot(&ba.cross(&ca)).abs()
98    }
99}
100
101impl From<[Cartesian<3>; 4]> for Simplex3 {
102    #[inline]
103    fn from(vertices: [Cartesian<3>; 4]) -> Self {
104        Simplex3 { vertices }.orient()
105    }
106}
107impl From<[[f64; 3]; 4]> for Simplex3 {
108    #[inline]
109    fn from(arrs: [[f64; 3]; 4]) -> Self {
110        Simplex3 {
111            vertices: arrs.map(Cartesian::from),
112        }
113        .orient()
114    }
115}
116
117impl Default for Simplex3 {
118    /// A uniform tetrahedron, centered on the origin with edges bisecting half of
119    /// the faces of a cube on their diagonals.
120    #[inline]
121    fn default() -> Self {
122        // Oriented by construction.
123        Simplex3 {
124            vertices: [
125                [1.0, 1.0, 1.0].into(),
126                [1.0, -1.0, -1.0].into(),
127                [-1.0, 1.0, -1.0].into(),
128                [-1.0, -1.0, 1.0].into(),
129            ],
130        }
131    }
132}
133
134impl fmt::Display for Simplex3 {
135    #[inline]
136    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
137        write!(
138            f,
139            "Simplex {{ [{}] }}",
140            self.vertices
141                .iter()
142                .map(Cartesian::to_string)
143                .collect::<Vec<String>>()
144                .join(", ")
145        )
146    }
147}
148impl Simplex3 {
149    /// The vertices of the simplex.
150    #[inline]
151    #[must_use]
152    pub fn vertices(&self) -> [Cartesian<3>; 4] {
153        self.vertices
154    }
155
156    /// Translate a simplex via rowwise addition of a Cartesian3
157    #[inline]
158    #[must_use]
159    pub fn translate_by(&mut self, rhs: &Cartesian<3>) -> Self {
160        Self {
161            vertices: self.vertices.map(|v| v + *rhs),
162        }
163    }
164
165    /// The center of mass of the tetrahedron.
166    #[inline]
167    #[must_use]
168    pub fn centroid(&self) -> Cartesian<3> {
169        array::from_fn::<_, 3, _>(|i| self.vertices.iter().fold(0.0, |acc, v| acc + v[i])).into()
170    }
171
172    /// Get the edges of the tetrahedron as edge endpoint coordinates. In vertex index
173    /// form, this returns values in the order [(1, 0), (2, 0), (3, 0), (2, 1), (3, 2)]
174    #[inline]
175    #[must_use]
176    pub fn get_edges(&self) -> [[Cartesian<3>; 2]; 5] {
177        [
178            [self.b(), self.a()],
179            [self.c(), self.a()],
180            [self.d(), self.a()],
181            [self.c(), self.b()],
182            [self.d(), self.b()],
183        ]
184    }
185
186    /// Edge vectors, in the same order as `get_edges` and pointing left to right.
187    #[inline]
188    #[must_use]
189    pub fn get_edge_vectors(&self) -> [Cartesian<3>; 5] {
190        self.get_edges().map(|[l, r]| l - r)
191    }
192
193    #[inline]
194    #[must_use]
195    /// 0th vertex of the tetrahedron
196    pub(crate) fn a(&self) -> Cartesian<3> {
197        self.vertices[0]
198    }
199    #[inline]
200    #[must_use]
201    /// 1st vertex of the tetrahedron
202    pub(crate) fn b(&self) -> Cartesian<3> {
203        self.vertices[1]
204    }
205    #[inline]
206    #[must_use]
207    /// 2nd vertex of the tetrahedron
208    pub(crate) fn c(&self) -> Cartesian<3> {
209        self.vertices[2]
210    }
211    #[inline]
212    #[must_use]
213    /// 3rd vertex of the tetrahedron
214    pub(crate) fn d(&self) -> Cartesian<3> {
215        self.vertices[3]
216    }
217    /// Return the vertices of an oriented tetrahedron. Users should call ``orient_in_place``
218    #[inline]
219    pub(crate) fn orient(&self) -> Simplex3 {
220        let dp = (self.d() - self.a()).dot(&((self.b() - self.a()).cross(&(self.c() - self.a()))));
221        if dp < 0.0 {
222            Simplex3 {
223                vertices: self.vertices,
224            }
225        } else {
226            Simplex3 {
227                vertices: [self.a(), self.c(), self.b(), self.d()],
228            }
229        }
230    }
231}
232/// Check if plane ``P_i`` defined by 4 coordinates and containing face ``i`` is a
233/// separating plane (or conversely, if the face normal is a separating axis)
234#[inline]
235#[must_use]
236fn check_face_on_p_is_separating(aff: &[f64; 4]) -> u8 {
237    aff.iter().enumerate().fold(
238        0u8,
239        |acc, (i, &x)| {
240            if x > 0.0 { acc | (1 << i) } else { acc }
241        },
242    )
243}
244
245/// Check whether plane ``P_j`` parallel to a face on q is a separating plane.
246#[inline]
247#[must_use]
248fn check_face_on_q_is_separating(aff: &[f64; 4]) -> bool {
249    aff.iter().all(|&x| x > 0.0)
250}
251
252/// Check that no point in our projected vertices lies in the (-, -) quadrant
253/// If a point satisfying that condition is found, there exists a separating line.
254#[expect(
255    clippy::too_many_arguments,
256    reason = "Internal function not exposed to users."
257)]
258#[inline]
259#[must_use]
260fn edge_test(
261    ma: u8,
262    mb: u8,
263    a: u8,
264    b: u8,
265    i: usize,
266    j: usize,
267    ea: &[f64; 4],
268    eb: &[f64; 4],
269) -> bool {
270    let cp = (ea[j] * eb[i]) - (ea[i] * eb[j]);
271    // NOTE: original paper used cp > | < 0.0, but we want to detect exact edge-edge overlaps.
272    ((ma & a) != 0 && (mb & b) != 0 && (cp >= 0.0)) || (ma & b) != 0 && (mb & a) != 0 && (cp <= 0.0)
273}
274
275/// Separating edge cases used in ``check_edge_is_separating``
276const _SEPARATING_EDGE_CASES: [(u8, u8, usize, usize); 6] = [
277    (1, 2, 0, 1),
278    (1, 4, 0, 2),
279    (1, 8, 0, 3),
280    (2, 4, 1, 2),
281    (2, 8, 1, 3),
282    (4, 8, 2, 3),
283];
284
285/// Check if there exists a separating plane containing the edge e shared by faces
286/// f0 and f1 of the tetrahedron.
287#[inline]
288#[must_use]
289fn check_edge_is_separating(aff_a: &[f64; 4], aff_b: &[f64; 4], ma: u8, mb: u8) -> bool {
290    if (ma | mb) != 15 {
291        return false; // If there is a vertex of b contained in the (-, -) quadrant
292    }
293    // exclude the vertices in (+,+) quadrant
294    let xa = ma & (ma ^ mb);
295    let xb = mb & (ma ^ mb);
296
297    if _SEPARATING_EDGE_CASES
298        .iter()
299        .any(|&(a, b, i, j)| edge_test(xa, xb, a, b, i, j, aff_a, aff_b))
300    {
301        return false;
302    }
303
304    // There exists a separating plane supported by the edge shared by f0 and f1
305    true
306}
307
308impl<R> IntersectsAtGlobal<Simplex3, Cartesian<3>, R> for Simplex3
309where
310    R: Rotation + Rotate<Cartesian<3>>,
311    RotationMatrix<3>: From<R>,
312{
313    #[inline]
314    fn intersects_at_global(
315        &self,
316        other: &Simplex3,
317        r_self: &Cartesian<3>,
318        o_self: &R,
319        r_other: &Cartesian<3>,
320        o_other: &R,
321    ) -> bool {
322        let (v_ij, o_ij) = hoomd_vector::pair_system_to_local(r_self, o_self, r_other, o_other);
323
324        self.intersects_at(other, &v_ij, &o_ij)
325    }
326}
327
328impl<R> IntersectsAt<Simplex3, Cartesian<3>, R> for Simplex3
329where
330    R: Copy,
331    RotationMatrix<3>: From<R>,
332{
333    /// Original C code of algorithm:
334    /// <https://web.archive.org/web/20030907000716/http://www.acm.org/jgt/papers/GanovelliPonchioRocchini02/tet_a_tet.html>
335    ///
336    /// Recent Clojure reimplementation that orients shapes:
337    /// <https://gist.github.com/postspectacular/9021724>
338    #[inline]
339    fn intersects_at(&self, other: &Simplex3, v_ij: &Cartesian<3>, o_ij: &R) -> bool {
340        let r = RotationMatrix::from(*o_ij);
341        let q = Simplex3::from(other.vertices.map(|v| r.rotate(&v))).translate_by(v_ij);
342
343        let q_deltas = q.vertices.map(|v| v - self.vertices[0]);
344
345        // Edge difference vectors for tetrahedron p
346        let mut edge_vectors_p = [Cartesian::<3>::default(); 5];
347        let mut masks = [0u8; 4];
348        let mut affs = [[0_f64; 4]; 4];
349
350        edge_vectors_p[0] = self.vertices[1] - self.vertices[0];
351        edge_vectors_p[1] = self.vertices[2] - self.vertices[0];
352
353        // Handle all possible SAT cases, in a performance-optimized order
354
355        // f 0 1
356        let n = edge_vectors_p[0].cross(&edge_vectors_p[1]);
357        affs[0] = q_deltas.map(|v| v.dot(&n));
358        let mask = check_face_on_p_is_separating(&affs[0]);
359        if mask == 15 {
360            return false;
361        }
362        masks[0] = mask;
363
364        // f 2 0
365        edge_vectors_p[2] = self.vertices[3] - self.vertices[0];
366        let n = edge_vectors_p[2].cross(&edge_vectors_p[0]);
367
368        affs[1] = q_deltas.map(|v| v.dot(&n));
369        let mask = check_face_on_p_is_separating(&affs[1]);
370        if mask == 15 {
371            return false;
372        }
373        masks[1] = mask;
374
375        // e 0 1
376        if check_edge_is_separating(&affs[0], &affs[1], masks[0], masks[1]) {
377            return false;
378        }
379
380        // f 1 2
381        let n = edge_vectors_p[1].cross(&edge_vectors_p[2]);
382        affs[2] = q_deltas.map(|v| v.dot(&n));
383        let mask = check_face_on_p_is_separating(&affs[2]);
384        if mask == 15 {
385            return false;
386        }
387        masks[2] = mask;
388
389        // e 0 2
390        if check_edge_is_separating(&affs[0], &affs[2], masks[0], masks[2]) {
391            return false;
392        }
393
394        // e 1 2
395        if check_edge_is_separating(&affs[1], &affs[2], masks[1], masks[2]) {
396            return false;
397        }
398
399        // f* 4 3
400        edge_vectors_p[3] = self.vertices[2] - self.vertices[1];
401        edge_vectors_p[4] = self.vertices[3] - self.vertices[1];
402
403        let n = edge_vectors_p[4].cross(&edge_vectors_p[3]);
404        affs[3] = q.vertices.map(|v| (v - self.vertices[1]).dot(&n));
405        let mask = check_face_on_p_is_separating(&affs[3]);
406        if mask == 15 {
407            return false;
408        }
409        masks[3] = mask;
410
411        // e 0|1|2 3
412        if [(0, 3), (1, 3), (2, 3)]
413            .iter()
414            .any(|&(i, j)| check_edge_is_separating(&affs[i], &affs[j], masks[i], masks[j]))
415        {
416            return false;
417        }
418
419        // If (at least) a vertex of q is inside tetrahedron p
420        // (vertex bounded by all 4 halfspaces)
421        if masks.iter().fold(0, |acc, &m| acc | m) != 15 {
422            return true;
423        }
424
425        // From now on, if there is a separating plane it is parallel to a face of q
426
427        // Difference between vertices on p and a vertex on q
428        let p_deltas = self.vertices.map(|v| v - q.vertices[0]);
429        let mut edge_vectors_q = [Cartesian::<3>::default(); 5];
430        edge_vectors_q[0] = q.vertices[1] - q.vertices[0];
431        edge_vectors_q[1] = q.vertices[2] - q.vertices[0];
432
433        // f 0 1
434        let n = edge_vectors_q[0].cross(&edge_vectors_q[1]);
435        if check_face_on_q_is_separating(&p_deltas.map(|v| v.dot(&n))) {
436            return false;
437        }
438
439        edge_vectors_q[2] = q.vertices[3] - q.vertices[0];
440
441        // f 2 0
442        let n = edge_vectors_q[2].cross(&edge_vectors_q[0]);
443        if check_face_on_q_is_separating(&p_deltas.map(|v| v.dot(&n))) {
444            return false;
445        }
446
447        // f 1 2
448        let n = edge_vectors_q[1].cross(&edge_vectors_q[2]);
449        if check_face_on_q_is_separating(&p_deltas.map(|v| v.dot(&n))) {
450            return false;
451        }
452        edge_vectors_q[3] = q.vertices[2] - q.vertices[1];
453        edge_vectors_q[4] = q.vertices[3] - q.vertices[1];
454
455        // f* 4 3
456        let n = edge_vectors_q[4].cross(&edge_vectors_q[3]);
457        let aff = self.vertices.map(|v| (v - other.vertices[1]).dot(&n));
458        if check_face_on_q_is_separating(&aff) {
459            return false;
460        }
461        true // No separating planes -> intersection!
462    }
463}
464
465impl BoundingSphereRadius for Simplex3 {
466    /// Radius of a sphere that bounds the shape.
467    ///
468    /// The sphere has the same local origin as the shape `self`.
469    ///
470    /// # Example
471    ///
472    /// ```
473    /// use approxim::assert_relative_eq;
474    /// use hoomd_geometry::{BoundingSphereRadius, shape::Simplex3};
475    ///
476    /// # fn main() -> Result<(), hoomd_geometry::Error> {
477    /// let tet = Simplex3::default();
478    ///
479    /// let bounding_radius = tet.bounding_sphere_radius();
480    ///
481    /// assert_relative_eq!(bounding_radius.get(), 3.0_f64.sqrt());
482    /// # Ok(())
483    /// # }
484    /// ```
485    #[inline]
486    fn bounding_sphere_radius(&self) -> PositiveReal {
487        self.vertices
488            .iter()
489            .map(Cartesian::norm_squared)
490            .fold(0.0, &f64::max)
491            .sqrt()
492            .try_into()
493            .expect("All norms are zero or NaN -- check your inputs!")
494    }
495}
496
497#[cfg(test)]
498mod tests {
499    use crate::xenocollide::collide3d;
500
501    use super::*;
502    use hoomd_vector::{Quaternion, Rotation, Versor};
503
504    use rstest::rstest;
505
506    #[test]
507    fn test_compute_mask() {
508        let arrays = (0u8..=15).map(|i| {
509            (
510                i,
511                [
512                    f64::from(i & 1),
513                    f64::from((i >> 1) & 1),
514                    f64::from((i >> 2) & 1),
515                    f64::from((i >> 3) & 1),
516                ],
517            )
518        });
519        arrays.for_each(|(i, arr)| assert_eq!(check_face_on_p_is_separating(&arr), i));
520    }
521
522    #[rstest(
523        ma, mb, ea, eb,
524        case(
525            0, 2,
526            [0.0, -100_000.0, 0.0, -500_000.0],
527            [0.0, 50_000.0, -500_000.0, 0.0]
528        ),
529        case(
530            6, 2,
531            [0.0, 1_025_000.0, 500_000.0, -500_000.0],
532            [0.0, 50_000.0, -500_000.0, 0.0]
533        ),
534        case(
535            0, 8,
536            [0.0, -100_000.0, 0.0, -500_000.0],
537            [-500_000.0, -1_475_000.0, -500_000.0, 500_000.0]
538        ),
539        case(
540            6, 8,
541            [0.0, 1_025_000.0, 500_000.0, -500_000.0],
542            [-500_000.0, -1_475_000.0, -500_000.0, 500_000.0]
543        ),
544        case(
545            2, 8,
546            [0.0, 50_000.0, -500_000.0, 0.0],
547            [-500_000.0, -1_475_000.0, -500_000.0, 500_000.0]
548        ),
549        case(
550            0, 0,
551            [0.0, 0.0, 0.0, -50_000.0],
552            [-500_005.0, -500_000.0, -1_000_000.0, -612_500.0]
553        ),
554        case(
555            0, 0,
556            [0.0, 0.0, 0.0, -50_000.0],
557            [0.0, -500_000.0, 0.0, -225_000.0]
558        ),
559        case(
560            0, 0,
561            [-500_005.0, -500_000.0, -1_000_000.0, -612_500.0],
562            [0.0, -500_000.0, 0.0, -225_000.0]
563        )
564    )]
565    fn test_edge_a(ma: u8, mb: u8, ea: [f64; 4], eb: [f64; 4]) {
566        // let t = Simplex3::default();
567        assert!(!check_edge_is_separating(&ea, &eb, ma, mb));
568    }
569    #[rstest(
570        n=> [[0.0, 0.0, -5000.0], [-5000.0, 5000.0, 1250.0], [0.0, -10000.0, 2500.0]],
571    )]
572    fn test_face_a(n: [f64; 3]) {
573        let deltas: [Cartesian<3>; 4] = [
574            [0.0, 0.0, 0.0],
575            [-200.0, 0.0, 20.0],
576            [-50.0, 50.0, 0.0],
577            [150.0, 25.0, 100.0],
578        ]
579        .map(Cartesian::from);
580        let aff = deltas.map(|v| v.dot(&n.into()));
581        let result = check_face_on_p_is_separating(&aff);
582        let expected_result = match n {
583            [0.0, 0.0, -5000.0] => (0, [0.0, -100_000.0, 0.0, -500_000.0]),
584            [-5000.0, 5000.0, 1250.0] => (6, [0.0, 1_025_000.0, 500_000.0, -500_000.0]),
585            [0.0, -10_000.0, 2500.0] => (2, [0.0, 50_000.0, -500_000.0, 0.0]),
586            _ => unreachable!(),
587        };
588        assert_eq!((result, aff), expected_result);
589    }
590
591    #[test]
592    fn test_tetisect() {
593        let mut p = Simplex3::from([
594            [0.0, 0.0, 0.0],
595            [50.0, 50.0, 0.0],
596            [100.0, 0.0, 0.0],
597            [50.0, 25.0, 100.0],
598        ]);
599        let mut q = Simplex3::from([
600            [0.0, 0.0, 0.0],
601            [-50.0, 50.0, 0.0],
602            [-200.0, 0.0, 20.0],
603            [150.0, 25.0, 100.0],
604        ]);
605
606        // First test case is constructed to intersect
607        assert!(p.intersects_at(&q, &Cartesian::default(), &Versor::identity()));
608
609        let p_centered = p.translate_by(&-p.centroid());
610        let q_centered = q.translate_by(&-q.centroid());
611
612        assert_eq!(
613            p.intersects_at(&q, &Cartesian::default(), &Versor::identity()),
614            collide3d(
615                &p_centered,
616                &q_centered,
617                &(q.centroid() - p.centroid()),
618                &Versor::identity()
619            )
620        );
621
622        let mut q_nooverlap = Simplex3::from([
623            [100.001, 0.0, 0.0],
624            [150.0, 50.0, 0.0],
625            [200.0, 0.0, 0.0],
626            [150.0, 25.0, 10.0],
627        ]);
628
629        // Second test case is constructed to NOT intersect
630        assert!(!p.intersects_at(&q_nooverlap, &Cartesian::default(), &Versor::identity()));
631        let q_nooverlap_centered = q_nooverlap.translate_by(&-p.centroid());
632        assert_eq!(
633            p.intersects_at(&q_nooverlap, &Cartesian::default(), &Versor::identity()),
634            collide3d(
635                &p_centered,
636                &q_nooverlap_centered,
637                &(q_nooverlap.centroid() - p.centroid()),
638                &Versor::identity()
639            )
640        );
641    }
642
643    #[rstest(
644        v_ij, o_ij, overlaps,
645        case::perfect_overlap(
646            [0.0, 0.0, 0.0].into(),
647            Versor::identity(),
648            true,
649        ),
650        case::particle_at_infinity(
651            [f64::INFINITY, 0.0, 0.0].into(),
652            Versor::identity(),
653            false,
654        ),
655        case::particle_at_negative_infinity(
656            [f64::NEG_INFINITY, 0.0, 0.0].into(),
657            Versor::identity(),
658            false,
659        ),
660        case::tip_tip_intersection_exact(
661            [2.0, 2.0, 2.0].into(),
662            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
663            true,
664        ),
665        case::tip_tip_intersection_imprecise(
666            [1.999_999_999, 1.999_999_999, 1.999_999_999].into(),
667            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
668            true,
669        ),
670        case::tip_tip_intersection_nooverlap(
671            [2.000_000_001, 2.000_000_001, 2.000_000_001].into(),
672            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
673            false,
674        ),
675        case::unrotated_tip_tip_intersection_exact(
676            [2.0, 2.0, 0.0].into(),
677            Versor::identity(),
678            true,
679        ),
680        case::unrotated_tip_tip_intersection_imprecise(
681            [1.999_999_999, 1.999_999_999, 0.0].into(),
682            Versor::identity(),
683            true,
684        ),
685        case::unrotated_tip_tip_intersection_nooverlap(
686            [2.000_000_001, 2.000_000_001, 0.0].into(),
687            Versor::identity(),
688            false,
689        ),
690        // NOTE: as above, this type of overlap was not counted in the original paper
691        case::tip_edge_intersection_exact(
692            [1.0, 1.0, 2.0].into(),
693            Versor::default(),
694            true,
695        ),
696        case::tip_edge_intersection_imprecise(
697            [1.0, 1.0, 1.999_999_999].into(),
698            Versor::default(),
699            true,
700        ),
701        case::tip_edge_intersection_nooverlap(
702            [1.0, 1.0, 2.000_000_001].into(),
703            Versor::default(),
704            false,
705        ),
706        case::parallel_edge_edge_intersection_exact(
707            [1.0, 1.0, 2.0].into(),
708            // Hard-coded quaternion required to prevent error accumulation
709            Quaternion::from([2.0_f64.sqrt() / 2.0, 2.0_f64.sqrt()/2.0, 0.0, 0.0]).to_versor_unchecked(),
710            true,
711        ),
712        case::parallel_edge_edge_intersection_imprecise(
713            [1.0, 1.0, 1.999_999_999].into(),
714            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
715            true,
716        ),
717        case::parallel_edge_edge_intersection_nooverlap(
718            [1.0, 1.0, 2.000_000_001].into(),
719            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
720            false,
721        ),
722        case::orthogonal_edge_edge_intersection_exact(
723            [1.0, 0.0, 2.0].into(),
724            Versor::identity(),
725            true,
726        ),
727        case::orthogonal_edge_edge_intersection_imprecise(
728            [1.0, 0.0, 1.999_999_999].into(),
729            Versor::identity(),
730            true,
731        ),
732        case::orthogonal_edge_edge_intersection_nooverlap(
733            [1.0, 0.0, 2.000_000_001].into(),
734            Versor::identity(),
735            false,
736        ),
737        case::nonorthogonal_edge_edge_intersection_exact(
738            [1.0, 0.0, 2.0].into(),
739            Versor::from_axis_angle([0.0, 0.0, 1.0].try_into().unwrap(), std::f64::consts::PI / 3.1),
740            true,
741        ),
742        case::nonorthogonal_edge_edge_intersection_imprecise(
743            [1.0, 0.0, 1.999_999_999].into(),
744            Versor::from_axis_angle([0.0, 0.0, 1.0].try_into().unwrap(), std::f64::consts::PI / 3.1),
745            true,
746        ),
747        case::nonorthogonal_edge_edge_intersection_nooverlap(
748            [1.0, 0.0, 2.000_000_001].into(),
749            Versor::from_axis_angle([0.0, 0.0, 1.0].try_into().unwrap(), std::f64::consts::PI / 3.1),
750            false,
751        ),
752        // Overlapping portions of the shape exactly align faces and edges
753        case::partial_aligned_overlap_exact(
754            [0.0, 1.0, -1.0].into(),
755            Versor::identity(),
756            true,
757        ),
758        case::partial_aligned_overlap_imprecise(
759            [0.0, 1.0, -0.999_999_999].into(),
760            Versor::identity(),
761            true,
762        ),
763        // Top edges are parallel, shapes partially within one another
764        case::partial_parallel_overlap(
765            [0.0, 0.0, -1.0].into(),
766            Versor::identity(),
767            true,
768        ),
769        // Vertex of p passes through q and lies on an edge
770        case::vertex_into_edge_shallow_exact(
771            [0.0, 1.0, 2.0].into(),
772            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_4),
773            true,
774        ),
775        case::vertex_into_edge_shallow_imprecise(
776            [0.0, 0.999_999_999, 2.0].into(),
777            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_4),
778            true,
779        ),
780        case::vertex_into_edge_deep_exact(
781            [0.0, 1.0, 1.0].into(),
782            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_4),
783            true,
784        ),
785        case::vertex_into_edge_deep_imprecise(
786            [0.0, 0.999_999_999, 1.0].into(),
787            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_4),
788            true,
789        ),
790        // Degenerate case that fails for both HOOMD-Blue and tet-a-tet
791        /*
792        case::vertex_face_imprecise(
793            [1.0, 1.0, 2.0].into(),
794            Versor::from_axis_angle([1.0, 0.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_4),
795            true,
796        ),
797        */
798        case::vertex_face_nooverlap(
799            [1.2765, -1.2765, 1.2765].into(),
800            Versor::from_axis_angle([1.0, 1.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
801            false,
802        ),
803        case::vertex_face_near_exact(
804            [1.275, -1.275, 1.275].into(),
805            Versor::from_axis_angle([1.0, 1.0, 0.0].try_into().unwrap(), std::f64::consts::FRAC_PI_2),
806            true,
807        ),
808    )]
809    fn test_tetrahedron_overlap_param(
810        #[values("intersects_at", "xenocollide")] method: &str,
811        v_ij: Cartesian<3>,
812        o_ij: Versor,
813        overlaps: bool,
814    ) {
815        let p = Simplex3::default();
816        let q = Simplex3::default();
817
818        let result = if method == "intersects_at" {
819            p.intersects_at(&q, &v_ij, &o_ij)
820        } else {
821            collide3d(&p, &q, &v_ij, &o_ij)
822        };
823
824        assert_eq!(
825            result, overlaps,
826            "Method `{method}` gave wrong overlap result.",
827        );
828    }
829
830    #[rstest]
831    fn test_tetrahedron_volume() {
832        let p = Simplex3::default();
833        assert_eq!(p.volume(), 8.0 / 3.0);
834    }
835}
hoomd_geometry/shape/simplex3.rs

hoomd_geometry/shape/
simplex3.rs
