hoomd_geometry/shape/cuboid.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 [`Hypercuboid`]
5
6use itertools::multizip;
7use rand::{
8 Rng,
9 distr::{Distribution, Uniform},
10};
11use serde::{Deserialize, Serialize};
12use serde_with::serde_as;
13use std::{array, ops::Mul};
14
15use crate::{BoundingSphereRadius, Error, IsPointInside, MapPoint, Scale, SupportMapping, Volume};
16use hoomd_utility::valid::PositiveReal;
17use hoomd_vector::Cartesian;
18
19/// A shape with with all perpendicular angles made from axis-aligned edges.
20///
21/// A [`Hypercuboid`] is the N-dimensional analog of a rectangle, and is defined by
22/// its edge lengths. Each perpendicular edge of the cuboid is aligned along the
23/// corresponding Cartesian axis. The hypercuboid is placed with its centroid at the
24/// origin.
25///
26/// # Example
27///
28/// Construction and basic methods:
29/// ```
30/// use hoomd_geometry::{Volume, shape::Hypercuboid};
31///
32/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
33/// let unit_cube = Hypercuboid {
34/// edge_lengths: [1.0.try_into()?; 3],
35/// };
36/// assert_eq!(unit_cube.volume(), 1.0);
37///
38/// let min_extents = unit_cube.minimal_extents();
39/// let max_extents = unit_cube.maximal_extents();
40/// assert_eq!(min_extents, [-0.5; 3]);
41/// assert_eq!(max_extents, [0.5; 3]);
42///
43/// let rectangular_prism = Hypercuboid {
44/// edge_lengths: [1.0.try_into()?, 1.0.try_into()?, 9.0.try_into()?],
45/// };
46///
47/// assert_eq!(rectangular_prism.volume(), 9.0);
48/// # Ok(())
49/// # }
50/// ```
51///
52/// Perform a fast AABB intersection tests:
53/// ```
54/// use hoomd_geometry::shape::Hypercuboid;
55///
56/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
57/// let unit_cube = Hypercuboid {
58/// edge_lengths: [1.0.try_into()?; 3],
59/// };
60/// let rectangular_prism = Hypercuboid {
61/// edge_lengths: [1.0.try_into()?, 1.0.try_into()?, 9.0.try_into()?],
62/// };
63///
64/// assert!(unit_cube.intersects_aligned(&rectangular_prism, &[1.0; 3].into()));
65/// assert!(
66/// !unit_cube.intersects_aligned(&rectangular_prism, &[1.1; 3].into())
67/// );
68/// # Ok(())
69/// # }
70/// ```
71///
72/// Wrap with [`Convex`](crate::Convex) to check intersections of oriented cuboids:
73///
74/// ```
75/// use hoomd_geometry::{Convex, IntersectsAt, shape::Rectangle};
76/// use hoomd_vector::Angle;
77/// use std::f64::consts::PI;
78///
79/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
80/// let square = Convex(Rectangle {
81/// edge_lengths: [1.0.try_into()?; 2],
82/// });
83///
84/// assert!(!square.intersects_at(
85/// &square,
86/// &[1.1, 0.0].into(),
87/// &Angle::default()
88/// ));
89/// assert!(square.intersects_at(
90/// &square,
91/// &[1.1, 0.0].into(),
92/// &Angle::from(PI / 4.0)
93/// ));
94/// # Ok(())
95/// # }
96/// ```
97#[serde_as]
98#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
99pub struct Hypercuboid<const N: usize> {
100 /// The lengths of each edge of the cuboid.
101 #[serde_as(as = "[_; N]")]
102 pub edge_lengths: [PositiveReal; N],
103}
104
105/// An axis-aligned rectangle.
106///
107/// # Examples
108///
109/// Basic construction and methods:
110/// ```
111/// use hoomd_geometry::{Volume, shape::Rectangle};
112///
113/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
114/// let rectangle = Rectangle {
115/// edge_lengths: [2.0.try_into()?, 4.0.try_into()?],
116/// };
117/// assert_eq!(rectangle.volume(), 8.0);
118/// # Ok(())
119/// # }
120/// ```
121///
122/// Intersection tests:
123/// ```
124/// use hoomd_geometry::{Convex, IntersectsAt, shape::Rectangle};
125/// use hoomd_vector::{Angle, Cartesian};
126/// use std::f64::consts::PI;
127///
128/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
129/// let rectangle = Rectangle {
130/// edge_lengths: [4.0.try_into()?, 2.0.try_into()?],
131/// };
132/// let rectangle = Convex(rectangle);
133///
134/// assert!(!rectangle.intersects_at(
135/// &rectangle,
136/// &[0.0, 2.1].into(),
137/// &Angle::default()
138/// ));
139/// assert!(rectangle.intersects_at(
140/// &rectangle,
141/// &[0.0, 2.1].into(),
142/// &Angle::from(PI / 2.0)
143/// ));
144/// # Ok(())
145/// # }
146/// ```
147pub type Rectangle = Hypercuboid<2>;
148
149/// An axis-aligned cuboid.
150///
151/// # Examples
152///
153/// Basic construction and methods:
154/// ```
155/// use hoomd_geometry::{Volume, shape::Cuboid};
156///
157/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
158/// let cuboid = Cuboid {
159/// edge_lengths: [2.0.try_into()?, 4.0.try_into()?, 7.0.try_into()?],
160/// };
161/// assert_eq!(cuboid.volume(), 56.0);
162/// # Ok(())
163/// # }
164/// ```
165///
166/// Intersection tests:
167/// ```
168/// use hoomd_geometry::{Convex, IntersectsAt, shape::Cuboid};
169/// use hoomd_vector::{Cartesian, Versor};
170///
171/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
172/// let cuboid = Cuboid {
173/// edge_lengths: [4.0.try_into()?, 2.0.try_into()?, 7.0.try_into()?],
174/// };
175/// let cuboid = Convex(cuboid);
176///
177/// assert!(!cuboid.intersects_at(
178/// &cuboid,
179/// &[0.0, 2.1, 0.0].into(),
180/// &Versor::default()
181/// ));
182/// assert!(cuboid.intersects_at(
183/// &cuboid,
184/// &[2.0, 1.0, 6.5].into(),
185/// &Versor::default()
186/// ));
187/// # Ok(())
188/// # }
189/// ```
190pub type Cuboid = Hypercuboid<3>;
191
192impl Hypercuboid<3> {
193 /// Length of the `Hypercuboid` edge along the x axis
194 #[inline]
195 #[must_use]
196 pub fn a(&self) -> PositiveReal {
197 self.edge_lengths[0]
198 }
199 /// Length of the `Hypercuboid` edge along the y axis
200 #[inline]
201 #[must_use]
202 pub fn b(&self) -> PositiveReal {
203 self.edge_lengths[1]
204 }
205 /// Length of the `Hypercuboid` edge along the z axis
206 #[inline]
207 #[must_use]
208 pub fn c(&self) -> PositiveReal {
209 self.edge_lengths[2]
210 }
211}
212
213impl<const N: usize> Hypercuboid<N> {
214 /// Construct a cuboid with all edge lengths equal.
215 ///
216 /// # Example
217 /// ```
218 /// use hoomd_geometry::shape::Rectangle;
219 ///
220 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
221 /// let square = Rectangle::with_equal_edges(10.0.try_into()?);
222 /// # Ok(())
223 /// # }
224 #[inline]
225 #[must_use]
226 pub fn with_equal_edges(l: PositiveReal) -> Self {
227 Self {
228 edge_lengths: [l; N],
229 }
230 }
231
232 /// Test for intersections between two *axis-aligned* cuboids.
233 ///
234 /// This test is much faster than a general oriented cuboid (OBB) intersection, which
235 /// can be achieved by wrapping with the [`Convex`](crate::Convex) newtype.
236 ///
237 /// # Example
238 ///
239 /// ```
240 /// use hoomd_geometry::shape::Hypercuboid;
241 ///
242 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
243 /// let unit_cube = Hypercuboid {
244 /// edge_lengths: [1.0.try_into()?; 3],
245 /// };
246 /// let rectangular_prism = Hypercuboid {
247 /// edge_lengths: [1.0.try_into()?, 1.0.try_into()?, 9.0.try_into()?],
248 /// };
249 ///
250 /// assert!(unit_cube.intersects_aligned(&rectangular_prism, &[1.0; 3].into()));
251 /// assert!(
252 /// !unit_cube.intersects_aligned(&rectangular_prism, &[1.1; 3].into())
253 /// );
254 /// # Ok(())
255 /// # }
256 /// ```
257 #[must_use]
258 #[inline]
259 pub fn intersects_aligned(&self, other: &Hypercuboid<N>, v_ij: &Cartesian<N>) -> bool {
260 let b_mins = Cartesian::from(other.minimal_extents()) + *v_ij;
261 let b_maxs = Cartesian::from(other.maximal_extents()) + *v_ij;
262 multizip((
263 self.minimal_extents(),
264 b_maxs,
265 self.maximal_extents(),
266 b_mins,
267 ))
268 .all(|(a_min, b_max, a_max, b_min)| (a_min <= b_max) && (a_max >= b_min))
269 }
270}
271
272impl<const N: usize> Volume for Hypercuboid<N> {
273 #[inline]
274 fn volume(&self) -> f64 {
275 self.edge_lengths
276 .iter()
277 .map(PositiveReal::get)
278 .reduce(f64::mul)
279 .expect("N should be >= 1")
280 }
281}
282
283impl<const N: usize> BoundingSphereRadius for Hypercuboid<N> {
284 #[inline]
285 fn bounding_sphere_radius(&self) -> PositiveReal {
286 f64::sqrt(
287 self.edge_lengths
288 .iter()
289 .map(PositiveReal::get)
290 .map(|x| (x / 2.0).powi(2))
291 .sum(),
292 )
293 .try_into()
294 .expect("expression evaluates to a positive real")
295 }
296}
297
298impl<const N: usize> SupportMapping<Cartesian<N>> for Hypercuboid<N> {
299 #[inline]
300 fn support_mapping(&self, n: &Cartesian<N>) -> Cartesian<N> {
301 let mut iter = n
302 .into_iter()
303 .zip(self.edge_lengths)
304 .map(|(n_i, l_i)| l_i.get() / 2.0 * n_i.signum());
305 array::from_fn(|_| iter.next().unwrap_or_default()).into()
306 }
307}
308
309impl<const N: usize> Hypercuboid<N> {
310 #[inline]
311 #[must_use]
312 /// Determine the maximal extents of the cuboid along each Cartesian axis.
313 ///
314 /// # Example
315 ///
316 /// ```
317 /// use hoomd_geometry::shape::Hypercuboid;
318 ///
319 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
320 /// let unit_cube = Hypercuboid {
321 /// edge_lengths: [1.0.try_into()?; 3],
322 /// };
323 ///
324 /// let max_extents = unit_cube.maximal_extents();
325 /// assert_eq!(max_extents, [0.5; 3]);
326 /// # Ok(())
327 /// # }
328 /// ```
329 pub fn maximal_extents(&self) -> [f64; N] {
330 array::from_fn(|i| self.edge_lengths[i].get() / 2.0)
331 }
332
333 #[inline]
334 #[must_use]
335 /// Determine the minimal extents of the cuboid along each Cartesian axis.
336 ///
337 /// # Example
338 ///
339 /// ```
340 /// use hoomd_geometry::shape::Hypercuboid;
341 ///
342 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
343 /// let unit_cube = Hypercuboid {
344 /// edge_lengths: [1.0.try_into()?; 3],
345 /// };
346 ///
347 /// let min_extents = unit_cube.minimal_extents();
348 /// assert_eq!(min_extents, [-0.5; 3]);
349 /// # Ok(())
350 /// # }
351 /// ```
352 pub fn minimal_extents(&self) -> [f64; N] {
353 array::from_fn(|i| -self.edge_lengths[i].get() / 2.0)
354 }
355}
356
357impl Hypercuboid<2> {
358 /// Represent the shape in the GSD box format.
359 ///
360 /// # Example
361 ///
362 /// ```
363 /// use hoomd_geometry::shape::Rectangle;
364 ///
365 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
366 /// let rectangle = Rectangle {
367 /// edge_lengths: [10.0.try_into()?, 15.0.try_into()?],
368 /// };
369 ///
370 /// let gsd_box = rectangle.to_gsd_box();
371 /// assert_eq!(gsd_box, [10.0, 15.0, 0.0, 0.0, 0.0, 0.0]);
372 /// # Ok(())
373 /// # }
374 /// ```
375 #[inline]
376 #[must_use]
377 pub fn to_gsd_box(&self) -> [f64; 6] {
378 [
379 self.edge_lengths[0].get(),
380 self.edge_lengths[1].get(),
381 0.0,
382 0.0,
383 0.0,
384 0.0,
385 ]
386 }
387}
388
389impl Hypercuboid<3> {
390 /// Represent the shape in the GSD box format.
391 ///
392 /// # Example
393 ///
394 /// ```
395 /// use hoomd_geometry::shape::Cuboid;
396 ///
397 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
398 /// let rectangle = Cuboid {
399 /// edge_lengths: [10.0.try_into()?, 15.0.try_into()?, 20.0.try_into()?],
400 /// };
401 ///
402 /// let gsd_box = rectangle.to_gsd_box();
403 /// assert_eq!(gsd_box, [10.0, 15.0, 20.0, 0.0, 0.0, 0.0]);
404 /// # Ok(())
405 /// # }
406 /// ```
407 #[inline]
408 #[must_use]
409 pub fn to_gsd_box(&self) -> [f64; 6] {
410 [
411 self.edge_lengths[0].get(),
412 self.edge_lengths[1].get(),
413 self.edge_lengths[2].get(),
414 0.0,
415 0.0,
416 0.0,
417 ]
418 }
419}
420
421impl<const N: usize> IsPointInside<Cartesian<N>> for Hypercuboid<N> {
422 /// Check if a cartesian vector is inside a cuboid.
423 ///
424 /// By conventions typically used in periodic boundary conditions, points
425 /// exactly at the minimal extent are inside the shape but points exactly
426 /// on the maximal extent are not:
427 /// ```math
428 /// -\frac{L_x}{2} \le x \lt \frac{L_x}{2}
429 /// ```
430 /// ```math
431 /// -\frac{L_y}{2} \le y \lt \frac{L_y}{2}
432 /// ```
433 /// ... and so on
434 ///
435 /// ```
436 /// use hoomd_geometry::{IsPointInside, shape::Hypercuboid};
437 ///
438 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
439 /// let cuboid = Hypercuboid {
440 /// edge_lengths: [6.0.try_into()?, 8.0.try_into()?],
441 /// };
442 ///
443 /// assert!(cuboid.is_point_inside(&[2.5, -3.5].into()));
444 /// assert!(!cuboid.is_point_inside(&[4.0, -3.5].into()));
445 /// # Ok(())
446 /// # }
447 /// ```
448 #[inline]
449 fn is_point_inside(&self, point: &Cartesian<N>) -> bool {
450 point
451 .into_iter()
452 .zip(&self.edge_lengths)
453 .all(|(x, l)| -l.get() / 2.0 <= x && x < l.get() / 2.0)
454 }
455}
456
457impl<const N: usize> Scale for Hypercuboid<N> {
458 /// Construct a scaled hypercuboid.
459 ///
460 /// The resulting hypercuboid's edge lengths $` L_\mathrm{new} `$ are
461 /// the original's $` L `$ scaled by $` v `$:
462 /// ```math
463 /// L_\mathrm{new} = v L
464 /// ```
465 ///
466 /// The centroid remains at the origin.
467 ///
468 /// # Example
469 ///
470 /// ```
471 /// use hoomd_geometry::{Scale, shape::Rectangle};
472 ///
473 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
474 /// let rectangle = Rectangle::with_equal_edges(10.0.try_into()?);
475 ///
476 /// let scaled_rectangle = rectangle.scale_length(0.5.try_into()?);
477 ///
478 /// assert_eq!(scaled_rectangle.edge_lengths[0].get(), 5.0);
479 /// # Ok(())
480 /// # }
481 /// ```
482 #[inline]
483 fn scale_length(&self, v: PositiveReal) -> Self {
484 Self {
485 edge_lengths: array::from_fn(|i| self.edge_lengths[i] * v),
486 }
487 }
488
489 /// Construct a scaled hypercuboid.
490 ///
491 /// The resulting hypercuboid's edge lengths $` L_\mathrm{new} `$ are
492 /// the original's $` L `$ scaled by $` v^\frac{1}{N} `$:
493 /// ```math
494 /// L_\mathrm{new} = v^\frac{1}{N} L
495 /// ```
496 ///
497 /// The centroid remains at the origin.
498 ///
499 /// # Examples
500 ///
501 /// ```
502 /// use hoomd_geometry::{Scale, shape::Rectangle};
503 ///
504 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
505 /// let rectangle = Rectangle::with_equal_edges(10.0.try_into()?);
506 ///
507 /// let scaled_rectangle = rectangle.scale_volume(4.0.try_into()?);
508 ///
509 /// assert_eq!(scaled_rectangle.edge_lengths[0].get(), 20.0);
510 /// # Ok(())
511 /// # }
512 /// ```
513 #[inline]
514 fn scale_volume(&self, v: PositiveReal) -> Self {
515 let v = v.get().powf(1.0 / N as f64);
516 self.scale_length(v.try_into().expect("v^{1/N} should be a positive real"))
517 }
518}
519
520impl<const N: usize> MapPoint<Cartesian<N>> for Hypercuboid<N> {
521 /// Map a point from one hypercuboid to another.
522 ///
523 /// Given a point P *inside `self`*, map it to the other shape
524 /// by scaling.
525 ///
526 /// # Errors
527 ///
528 /// Returns [`Error::PointOutsideShape`] when `point` is outside the shape
529 /// `self`.
530 ///
531 /// # Example
532 ///
533 /// ```
534 /// use hoomd_geometry::{MapPoint, shape::Rectangle};
535 /// use hoomd_vector::Cartesian;
536 ///
537 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
538 /// let rectangle_a = Rectangle::with_equal_edges(10.0.try_into()?);
539 /// let rectangle_b = Rectangle::with_equal_edges(20.0.try_into()?);
540 ///
541 /// let mapped_point =
542 /// rectangle_a.map_point(Cartesian::from([-1.0, 1.0]), &rectangle_b);
543 ///
544 /// assert_eq!(mapped_point, Ok(Cartesian::from([-2.0, 2.0])));
545 /// assert_eq!(
546 /// rectangle_a.map_point(Cartesian::from([-100.0, 1.0]), &rectangle_b),
547 /// Err(hoomd_geometry::Error::PointOutsideShape)
548 /// );
549 /// # Ok(())
550 /// # }
551 /// ```
552 #[inline]
553 fn map_point(&self, point: Cartesian<N>, other: &Self) -> Result<Cartesian<N>, Error> {
554 if !self.is_point_inside(&point) {
555 return Err(Error::PointOutsideShape);
556 }
557
558 let scale: [_; N] = array::from_fn(|i| other.edge_lengths[i] / self.edge_lengths[i]);
559 Ok(Cartesian::from(array::from_fn(|i| {
560 (scale[i].get() * point[i]).clamp(
561 -other.edge_lengths[i].get() / 2.0,
562 (other.edge_lengths[i].get() / 2.0).next_down(),
563 )
564 })))
565 }
566}
567
568impl<const N: usize> Distribution<Cartesian<N>> for Hypercuboid<N> {
569 /// Generate points uniformly distributed in the cuboid.
570 ///
571 /// # Example
572 ///
573 /// ```
574 /// use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
575 ///
576 /// use hoomd_geometry::{IsPointInside, shape::Hypercuboid};
577 ///
578 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
579 /// let cuboid = Hypercuboid {
580 /// edge_lengths: [6.0.try_into()?, 8.0.try_into()?],
581 /// };
582 /// let mut rng = StdRng::seed_from_u64(1);
583 ///
584 /// let point = cuboid.sample(&mut rng);
585 /// assert!(cuboid.is_point_inside(&point));
586 /// # Ok(())
587 /// # }
588 /// ```
589 #[inline]
590 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Cartesian<N> {
591 let minimal_extents = self.minimal_extents();
592 let maximal_extents = self.maximal_extents();
593
594 array::from_fn(|i| {
595
596 let uniform = Uniform::new(minimal_extents[i], maximal_extents[i])
597 .expect("cuboid should always have real valued extents where the minimum is less than the maximum");
598 uniform.sample(rng)}).into()
599 }
600}
601
602#[cfg(test)]
603#[expect(clippy::used_underscore_binding, reason = "Required for const tests.")]
604mod tests {
605 use super::*;
606 use approxim::assert_relative_eq;
607 use assert2::check;
608 use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
609 use rstest::*;
610 use std::marker::PhantomData;
611
612 /// Number of random samples to test.
613 const N: usize = 1024;
614
615 #[rstest(
616 edges0 => [[2.0.try_into().expect("test value is a positive real"), 2.0.try_into().expect("test value is a positive real"), 2.0.try_into().expect("test value is a positive real")]],
617 edges1 => [[1.0.try_into().expect("test value is a positive real"), 1.0.try_into().expect("test value is a positive real"), 1.0.try_into().expect("test value is a positive real")]],
618 )]
619 fn test_box_intersections_aligned(edges0: [PositiveReal; 3], edges1: [PositiveReal; 3]) {
620 let (s0, s1) = (
621 Hypercuboid {
622 edge_lengths: edges0,
623 },
624 Hypercuboid {
625 edge_lengths: edges1,
626 },
627 );
628 // Should all be false (no intersection), which we invert to true
629 check!(!s0.intersects_aligned(&s1, &[10.0, 10.0, 10.0].into()));
630 // Boundaries are aligned
631 check!(s0.intersects_aligned(&s1, &[1.5, 1.5, 1.5].into()));
632 // Both at origin - will always intersect for any cuboids
633 check!(s0.intersects_aligned(&s1, &[0.0, 0.0, 0.0].into()));
634 }
635 #[rstest(
636 edges0 => [[2.0.try_into().expect("test value is a positive real"), 2.0.try_into().expect("test value is a positive real")]],
637 edges1 => [[1.0.try_into().expect("test value is a positive real"), 1.0.try_into().expect("test value is a positive real")]],
638 )]
639 fn test_box_intersections_2d_aligned(edges0: [PositiveReal; 2], edges1: [PositiveReal; 2]) {
640 let (c0, c1) = (
641 Hypercuboid {
642 edge_lengths: edges0,
643 },
644 Hypercuboid {
645 edge_lengths: edges1,
646 },
647 );
648 // Should all be false (no intersection), which we invert to true
649 check!(!c0.intersects_aligned(&c1, &[10.0, 10.0].into()));
650 // Boundaries are aligned
651 check!(c0.intersects_aligned(&c1, &[1.5, 1.5].into()));
652 // Both at origin - will always intersect for any cuboids
653 check!(c0.intersects_aligned(&c1, &[0.0, 0.0].into()));
654 }
655
656 #[rstest(
657 _n => [
658 PhantomData::<Hypercuboid<1>>,
659 PhantomData::<Hypercuboid<2>>,
660 PhantomData::<Hypercuboid<3>>,
661 PhantomData::<Hypercuboid<4>>
662 ],
663 l => [1e-6, 1.0, 3.456, 99_999_999.9],
664 )]
665 fn test_box_extents<const N: usize>(_n: PhantomData<Hypercuboid<N>>, l: f64) {
666 let c = Hypercuboid {
667 edge_lengths: [l.try_into().expect("test value is a positive real"); N],
668 };
669 check!(c.maximal_extents() == [l / 2.0; N]);
670 check!(c.minimal_extents() == [-l / 2.0; N]);
671 }
672
673 #[rstest(
674 _n => [
675 PhantomData::<Hypercuboid<1>>,
676 PhantomData::<Hypercuboid<2>>,
677 PhantomData::<Hypercuboid<3>>,
678 PhantomData::<Hypercuboid<4>>
679 ],
680 l => [1e-6, 1.0, 3.456, 99_999_999.9],
681 )]
682 fn test_box_volume<const N: usize>(
683 _n: PhantomData<Hypercuboid<N>>,
684 l: f64,
685 ) -> anyhow::Result<()> {
686 let c = Hypercuboid {
687 edge_lengths: [l.try_into()?; N],
688 };
689 assert_relative_eq!(
690 c.volume(),
691 if N != 0 {
692 l.powi(i32::try_from(N).unwrap())
693 } else {
694 0.0
695 }
696 );
697
698 Ok(())
699 }
700
701 #[rstest(
702 l => [1e-6, 1.0, 3.456, 99_999_999.9],
703 )]
704 fn test_box_abc(l: f64) -> anyhow::Result<()> {
705 let c = Hypercuboid {
706 edge_lengths: [l.try_into()?; 3],
707 };
708 check!([c.a(), c.b(), c.c()] == [l.try_into()?; 3]);
709
710 Ok(())
711 }
712
713 #[test]
714 fn bounding_sphere_radius_2d() -> anyhow::Result<()> {
715 let cuboid = Hypercuboid {
716 edge_lengths: [1.0.try_into()?, 1.0.try_into()?],
717 };
718 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 2.0_f64.sqrt() / 2.0);
719
720 let cuboid = Hypercuboid {
721 edge_lengths: [2.0.try_into()?, 2.0.try_into()?],
722 };
723 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 2.0_f64.sqrt());
724
725 let cuboid = Hypercuboid {
726 edge_lengths: [6.0.try_into()?, 8.0.try_into()?],
727 };
728 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 5.0);
729
730 Ok(())
731 }
732
733 #[test]
734 fn bounding_sphere_radius_3d() -> anyhow::Result<()> {
735 let cuboid = Hypercuboid {
736 edge_lengths: [1.0.try_into()?, 1.0.try_into()?, 1.0.try_into()?],
737 };
738 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 3.0_f64.sqrt() / 2.0);
739
740 let cuboid = Hypercuboid {
741 edge_lengths: [2.0.try_into()?, 2.0.try_into()?, 2.0.try_into()?],
742 };
743 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 3.0_f64.sqrt());
744
745 let cuboid = Hypercuboid {
746 edge_lengths: [2.0.try_into()?, 4.0.try_into()?, 6.0.try_into()?],
747 };
748 assert_relative_eq!(cuboid.bounding_sphere_radius().get(), 14.0_f64.sqrt());
749
750 Ok(())
751 }
752
753 #[test]
754 fn support_mapping_2d() -> anyhow::Result<()> {
755 let cuboid = Hypercuboid {
756 edge_lengths: [2.0.try_into()?, 4.0.try_into()?],
757 };
758
759 assert_relative_eq!(
760 cuboid.support_mapping(&Cartesian::from([1.0, 0.1])),
761 [1.0, 2.0].into()
762 );
763 assert_relative_eq!(
764 cuboid.support_mapping(&Cartesian::from([1.0, -0.1])),
765 [1.0, -2.0].into()
766 );
767 assert_relative_eq!(
768 cuboid.support_mapping(&Cartesian::from([-0.1, 1.0])),
769 [-1.0, 2.0].into()
770 );
771 assert_relative_eq!(
772 cuboid.support_mapping(&Cartesian::from([-0.1, -1.0])),
773 [-1.0, -2.0].into()
774 );
775
776 Ok(())
777 }
778
779 #[test]
780 fn support_mapping_3d() -> anyhow::Result<()> {
781 let cuboid = Hypercuboid {
782 edge_lengths: [2.0.try_into()?, 4.0.try_into()?, 6.0.try_into()?],
783 };
784
785 assert_relative_eq!(
786 cuboid.support_mapping(&Cartesian::from([1.0, 0.1, 0.1])),
787 [1.0, 2.0, 3.0].into()
788 );
789 assert_relative_eq!(
790 cuboid.support_mapping(&Cartesian::from([1.0, 0.1, -0.1])),
791 [1.0, 2.0, -3.0].into()
792 );
793 assert_relative_eq!(
794 cuboid.support_mapping(&Cartesian::from([1.0, -0.1, 0.1])),
795 [1.0, -2.0, 3.0].into()
796 );
797 assert_relative_eq!(
798 cuboid.support_mapping(&Cartesian::from([1.0, -0.1, -0.1])),
799 [1.0, -2.0, -3.0].into()
800 );
801 assert_relative_eq!(
802 cuboid.support_mapping(&Cartesian::from([-1.0, 0.1, 0.1])),
803 [-1.0, 2.0, 3.0].into()
804 );
805 assert_relative_eq!(
806 cuboid.support_mapping(&Cartesian::from([-1.0, 0.1, -0.1])),
807 [-1.0, 2.0, -3.0].into()
808 );
809 assert_relative_eq!(
810 cuboid.support_mapping(&Cartesian::from([-1.0, -0.1, 0.1])),
811 [-1.0, -2.0, 3.0].into()
812 );
813 assert_relative_eq!(
814 cuboid.support_mapping(&Cartesian::from([-1.0, -0.1, -0.1])),
815 [-1.0, -2.0, -3.0].into()
816 );
817
818 Ok(())
819 }
820
821 #[test]
822 fn is_point_inside() -> anyhow::Result<()> {
823 let cuboid = Hypercuboid {
824 edge_lengths: [2.0.try_into()?, 4.0.try_into()?],
825 };
826
827 check!(cuboid.is_point_inside(&Cartesian::from([0.0, 0.0])));
828 check!(cuboid.is_point_inside(&Cartesian::from([-1.0, 0.0])));
829 check!(cuboid.is_point_inside(&Cartesian::from([0.0, -2.0])));
830 check!(cuboid.is_point_inside(&Cartesian::from([-1.0, -2.0])));
831 check!(cuboid.is_point_inside(&Cartesian::from([0.5, -1.0])));
832
833 check!(!cuboid.is_point_inside(&Cartesian::from([1.0, 0.0])));
834 check!(!cuboid.is_point_inside(&Cartesian::from([0.0, 2.0])));
835 check!(!cuboid.is_point_inside(&Cartesian::from([1.0, 2.0])));
836 check!(!cuboid.is_point_inside(&Cartesian::from([10.0, -20.0])));
837
838 Ok(())
839 }
840
841 #[test]
842 fn distribution() -> anyhow::Result<()> {
843 let cuboid = Hypercuboid {
844 edge_lengths: [6.0.try_into()?, 10.0.try_into()?],
845 };
846 let mut rng = StdRng::seed_from_u64(3);
847
848 let points: Vec<_> = (&cuboid).sample_iter(&mut rng).take(N).collect();
849 check!(&points.iter().all(|p| cuboid.is_point_inside(p)));
850 check!(&points.iter().any(|p| p[0] < -2.8));
851 check!(&points.iter().any(|p| p[0] > 2.8));
852 check!(&points.iter().any(|p| p[1] < -4.8));
853 check!(&points.iter().any(|p| p[1] > 4.8));
854
855 Ok(())
856 }
857
858 #[test]
859 fn test_scale_length() -> anyhow::Result<()> {
860 let cuboid = Hypercuboid {
861 edge_lengths: [6.0.try_into()?, 10.0.try_into()?],
862 };
863
864 let scaled_cuboid = cuboid.scale_length(2.0.try_into()?);
865 check!(scaled_cuboid.edge_lengths[0].get() == 12.0);
866 check!(scaled_cuboid.edge_lengths[1].get() == 20.0);
867
868 let scaled_cuboid = cuboid.scale_length(0.5.try_into()?);
869 check!(scaled_cuboid.edge_lengths[0].get() == 3.0);
870 check!(scaled_cuboid.edge_lengths[1].get() == 5.0);
871
872 Ok(())
873 }
874
875 #[test]
876 fn test_scale_volume() -> anyhow::Result<()> {
877 let cuboid = Hypercuboid {
878 edge_lengths: [6.0.try_into()?, 10.0.try_into()?],
879 };
880
881 let scaled_cuboid = cuboid.scale_volume(4.0.try_into()?);
882 check!(scaled_cuboid.edge_lengths[0].get() == 12.0);
883 check!(scaled_cuboid.edge_lengths[1].get() == 20.0);
884
885 let scaled_cuboid = cuboid.scale_volume(0.25.try_into()?);
886 check!(scaled_cuboid.edge_lengths[0].get() == 3.0);
887 check!(scaled_cuboid.edge_lengths[1].get() == 5.0);
888
889 Ok(())
890 }
891
892 #[test]
893 fn test_map_basic() -> anyhow::Result<()> {
894 let cuboid_a = Hypercuboid {
895 edge_lengths: [6.0.try_into()?, 10.0.try_into()?],
896 };
897
898 let cuboid_b = Hypercuboid {
899 edge_lengths: [12.0.try_into()?, 5.0.try_into()?],
900 };
901
902 check!(
903 cuboid_a.map_point(Cartesian::from([0.0, 0.0]), &cuboid_b)
904 == Ok(Cartesian::from([0.0, 0.0]))
905 );
906 check!(
907 cuboid_b.map_point(Cartesian::from([0.0, 0.0]), &cuboid_a)
908 == Ok(Cartesian::from([0.0, 0.0]))
909 );
910
911 check!(
912 cuboid_a.map_point(Cartesian::from([100.0, 0.0]), &cuboid_b)
913 == Err(Error::PointOutsideShape)
914 );
915 check!(
916 cuboid_b.map_point(Cartesian::from([0.0, -200.0]), &cuboid_a)
917 == Err(Error::PointOutsideShape)
918 );
919
920 check!(
921 cuboid_a.map_point(Cartesian::from([2.0, 1.0]), &cuboid_b)
922 == Ok(Cartesian::from([4.0, 0.5]))
923 );
924 check!(
925 cuboid_b.map_point(Cartesian::from([-4.0, 0.5]), &cuboid_a)
926 == Ok(Cartesian::from([-2.0, 1.0]))
927 );
928
929 check!(
930 cuboid_a.map_point(Cartesian::from([-3.0, -5.0]), &cuboid_b)
931 == Ok(Cartesian::from([-6.0, -2.5]))
932 );
933 check!(
934 cuboid_b.map_point(Cartesian::from([-6.0, -2.5]), &cuboid_a)
935 == Ok(Cartesian::from([-3.0, -5.0]))
936 );
937
938 Ok(())
939 }
940
941 #[test]
942 fn test_map_corner() -> anyhow::Result<()> {
943 let mut rng = StdRng::seed_from_u64(3);
944 let uniform = Uniform::new(1.0, 1000.0)?;
945
946 for _ in 0..65536 {
947 let a = uniform.sample(&mut rng);
948 let b = uniform.sample(&mut rng);
949 let cuboid_a = Hypercuboid::<2>::with_equal_edges(a.try_into()?);
950 let cuboid_b = Hypercuboid::<2>::with_equal_edges(b.try_into()?);
951
952 // Test that points right on the boundary of one shape remain inside
953 // the other shape. If not implemented correctly, map_point might
954 // round and place a point just outside the shape. This test fails
955 // when the `.clamp` call in `map_point` is commented out.
956
957 let lower_left =
958 cuboid_a.map_point(Cartesian::from([-a / 2.0, -a / 2.0]), &cuboid_b)?;
959 check!(
960 cuboid_b.is_point_inside(&lower_left),
961 "{lower_left:?} should be inside {cuboid_b:?}"
962 );
963
964 let upper_right = cuboid_a.map_point(
965 Cartesian::from([(a / 2.0).next_down(), (a / 2.0).next_down()]),
966 &cuboid_b,
967 )?;
968 check!(
969 cuboid_b.is_point_inside(&upper_right),
970 "{upper_right:?} should be inside {cuboid_b:?}"
971 );
972 }
973
974 Ok(())
975 }
976}