diff --git a/crates/bevy_math/src/primitives/dim2.rs b/crates/bevy_math/src/primitives/dim2.rs
index c69491fac2..d7c47e65dc 100644
--- a/crates/bevy_math/src/primitives/dim2.rs
+++ b/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 {
diff --git a/crates/bevy_math/src/primitives/dim3.rs b/crates/bevy_math/src/primitives/dim3.rs
index e1b5ae5128..1c2e217ede 100644
--- a/crates/bevy_math/src/primitives/dim3.rs
+++ b/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 {
diff --git a/crates/bevy_math/src/sampling/shape_sampling.rs b/crates/bevy_math/src/sampling/shape_sampling.rs
index 3728034413..ed3359b832 100644
--- a/crates/bevy_math/src/sampling/shape_sampling.rs
+++ b/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;