Faces by name in the face commands; model.faces gives a point on each face (#1243)

crates/kernel/src/surface_kind.rs

@@ -50,3 +50,64 @@ pub fn face_surface(b: &Body, face: usize) -> Option<FaceSurface> {
         Surf::Torus { o, a, big, h, r, .. } => FaceSurface::Torus { center: o + a * h, axis: a, major: big, minor: r },
     })
 }
+
+impl FaceSurface {
+    /// The point of the whole surface (not only of its face) nearest `p`. `None` where every
+    /// point of a ring is as near (on a cylinder's or cone's axis, at a sphere's or torus's
+    /// centre), and for a point that is not finite.
+    pub fn nearest(&self, p: Vec3) -> Option<Vec3> {
+        if !p.is_finite() {
+            return None;
+        }
+        match *self {
+            FaceSurface::Plane { normal, point } => Some(p - normal * (p - point).dot(normal)),
+            FaceSurface::Cylinder { origin, axis, radius } => {
+                let w = p - origin;
+                let along = axis * w.dot(axis);
+                Some(origin + along + (w - along).normalized()? * radius)
+            }
+            FaceSurface::Cone { apex, axis, half_angle } => {
+                let w = p - apex;
+                let out = (w - axis * w.dot(axis)).normalized()?;
+                // The cone's line through the apex that p is nearest, in the plane of the axis and p.
+                let line = axis * half_angle.cos() + out * half_angle.sin();
+                Some(apex + line * w.dot(line).max(0.0))
+            }
+            FaceSurface::Sphere { center, radius } => Some(center + (p - center).normalized()? * radius),
+            FaceSurface::Torus { center, axis, major, minor } => {
+                let w = p - center;
+                let tube = center + (w - axis * w.dot(axis)).normalized()? * major;
+                Some(tube + (p - tube).normalized()? * minor)
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn near(a: Option<Vec3>, b: Vec3) -> bool {
+        a.is_some_and(|a| a.dist(b) < 1e-12)
+    }
+
+    /// Each surface's nearest point, worked out by hand.
+    #[test]
+    fn the_nearest_point_of_each_surface() {
+        let plane = FaceSurface::Plane { normal: Vec3::Z, point: Vec3::new(0.0, 0.0, 4.0) };
+        assert!(near(plane.nearest(Vec3::new(3.0, 5.0, 9.0)), Vec3::new(3.0, 5.0, 4.0)));
+        // A cylinder of radius 5 about the line x = 1, y = 2: (7, 10, 3) is 10 from it along (3, 4)/5.
+        let cylinder = FaceSurface::Cylinder { origin: Vec3::new(1.0, 2.0, 0.0), axis: Vec3::Z, radius: 5.0 };
+        assert!(near(cylinder.nearest(Vec3::new(7.0, 10.0, 3.0)), Vec3::new(4.0, 6.0, 3.0)));
+        assert!(cylinder.nearest(Vec3::new(1.0, 2.0, 3.0)).is_none(), "on the axis every direction is as near");
+        let sphere = FaceSurface::Sphere { center: Vec3::new(1.0, 2.0, 3.0), radius: 5.0 };
+        assert!(near(sphere.nearest(Vec3::new(1.0, 2.0, 13.0)), Vec3::new(1.0, 2.0, 8.0)));
+        // A torus of ring radius 10 and tube radius 2 about Z: (14, 0, 0) is 2 beyond the tube's outside.
+        let torus = FaceSurface::Torus { center: Vec3::ZERO, axis: Vec3::Z, major: 10.0, minor: 2.0 };
+        assert!(near(torus.nearest(Vec3::new(14.0, 0.0, 0.0)), Vec3::new(12.0, 0.0, 0.0)));
+        // A cone of half-angle 45 degrees: (2, 0, 0) is √2 along its line (1, 0, 1)/√2.
+        let cone = FaceSurface::Cone { apex: Vec3::ZERO, axis: Vec3::Z, half_angle: std::f64::consts::FRAC_PI_4 };
+        assert!(cone.nearest(Vec3::new(2.0, 0.0, 0.0)).is_some_and(|q| q.dist(Vec3::new(1.0, 0.0, 1.0)) < 1e-12));
+        assert!(plane.nearest(Vec3::new(f64::NAN, 0.0, 0.0)).is_none());
+    }
+}