Tetrahedron sampling (#13430)

# Objective

Add interior and boundary sampling for the `Tetrahedron` primitive. This
is part of ongoing work to bring the primitives to parity with each
other in terms of their capabilities.

## Solution

`Tetrahedron` implements the `ShapeSample` trait. To support this, there
is a new public method `Tetrahedron::faces` which gets the faces of a
tetrahedron as `Triangle3d`s. There are more sophisticated ideas for
getting the faces we might want to consider in the future (e.g.
adjusting according to the orientation), but this method gives the most
mathematically straightforward answer, giving the faces the orientation
induced by the tetrahedron itself.

crates/bevy_math/src/primitives/dim2.rs

@@ -528,11 +528,19 @@ impl Triangle2d {
     }
 
     /// Reverse the [`WindingOrder`] of the triangle
-    /// by swapping the first and last vertices
+    /// by swapping the first and last vertices.
     #[inline(always)]
     pub fn reverse(&mut self) {
         self.vertices.swap(0, 2);
     }
+
+    /// This triangle but reversed.
+    #[inline(always)]
+    #[must_use]
+    pub fn reversed(mut self) -> Self {
+        self.reverse();
+        self
+    }
 }
 
 impl Measured2d for Triangle2d {

crates/bevy_math/src/primitives/dim3.rs

@@ -786,6 +786,14 @@ impl Triangle3d {
         self.vertices.swap(0, 2);
     }
 
+    /// This triangle but reversed.
+    #[inline(always)]
+    #[must_use]
+    pub fn reversed(mut self) -> Triangle3d {
+        self.reverse();
+        self
+    }
+
     /// Get the centroid of the triangle.
     ///
     /// This function finds the geometric center of the triangle by averaging the vertices:
@@ -915,6 +923,22 @@ impl Tetrahedron {
     pub fn centroid(&self) -> Vec3 {
         (self.vertices[0] + self.vertices[1] + self.vertices[2] + self.vertices[3]) / 4.0
     }
+
+    /// Get the triangles that form the faces of this tetrahedron.
+    ///
+    /// Note that the orientations of the faces are determined by that of the tetrahedron; if the
+    /// signed volume of this tetrahedron is positive, then the triangles' normals will point
+    /// outward, and if the signed volume is negative they will point inward.
+    #[inline(always)]
+    pub fn faces(&self) -> [Triangle3d; 4] {
+        let [a, b, c, d] = self.vertices;
+        [
+            Triangle3d::new(b, c, d),
+            Triangle3d::new(a, c, d).reversed(),
+            Triangle3d::new(a, b, d),
+            Triangle3d::new(a, b, c).reversed(),
+        ]
+    }
 }
 
 impl Measured3d for Tetrahedron {

crates/bevy_math/src/sampling/shape_sampling.rs

@@ -215,6 +215,63 @@ impl ShapeSample for Triangle3d {
     }
 }
 
+impl ShapeSample for Tetrahedron {
+    type Output = Vec3;
+
+    fn sample_interior<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
+        let [v0, v1, v2, v3] = self.vertices;
+
+        // Generate a random point in a cube:
+        let mut coords: [f32; 3] = [
+            rng.gen_range(0.0..1.0),
+            rng.gen_range(0.0..1.0),
+            rng.gen_range(0.0..1.0),
+        ];
+
+        // The cube is broken into six tetrahedra of the form 0 <= c_0 <= c_1 <= c_2 <= 1,
+        // where c_i are the three euclidean coordinates in some permutation. (Since 3! = 6,
+        // there are six of them). Sorting the coordinates folds these six tetrahedra into the
+        // tetrahedron 0 <= x <= y <= z <= 1 (i.e. a fundamental domain of the permutation action).
+        coords.sort_by(|x, y| x.partial_cmp(y).unwrap());
+
+        // Now, convert a point from the fundamental tetrahedron into barycentric coordinates by
+        // taking the four successive differences of coordinates; note that these telescope to sum
+        // to 1, and this transformation is linear, hence preserves the probability density, since
+        // the latter comes from the Lebesgue measure.
+        //
+        // (See https://en.wikipedia.org/wiki/Lebesgue_measure#Properties — specifically, that
+        // Lebesgue measure of a linearly transformed set is its original measure times the
+        // determinant.)
+        let (a, b, c, d) = (
+            coords[0],
+            coords[1] - coords[0],
+            coords[2] - coords[1],
+            1. - coords[2],
+        );
+
+        // This is also a linear mapping, so probability density is still preserved.
+        v0 * a + v1 * b + v2 * c + v3 * d
+    }
+
+    fn sample_boundary<R: Rng + ?Sized>(&self, rng: &mut R) -> Self::Output {
+        let triangles = self.faces();
+        let areas = triangles.iter().map(|t| t.area());
+
+        if areas.clone().sum::<f32>() > 0.0 {
+            // There is at least one triangle with nonzero area, so this unwrap succeeds.
+            let dist = WeightedIndex::new(areas).unwrap();
+
+            // Get a random index, then sample the interior of the associated triangle.
+            let idx = dist.sample(rng);
+            triangles[idx].sample_interior(rng)
+        } else {
+            // In this branch the tetrahedron has zero surface area; just return a point that's on
+            // the tetrahedron.
+            self.vertices[0]
+        }
+    }
+}
+
 impl ShapeSample for Cylinder {
     type Output = Vec3;