Merge pull request #2453 from hannobraun/geom

Add `SurfaceGeom::triangle_at`; rewrite point/vector conversion on top of it

crates/fj-core/src/algorithms/approx/circle.rs

@@ -59,11 +59,13 @@ pub fn approx_circle<const D: usize>(
     points
 }
 
-struct PathApproxParams {
+/// Path approximation parameters for a circle
+pub struct PathApproxParams {
     increment: Scalar,
 }
 
 impl PathApproxParams {
+    /// Compute path approximation parameters for the given circle and tolerance
     pub fn for_circle<const D: usize>(
         circle: &Circle<D>,
         tolerance: impl Into<Tolerance>,
@@ -82,10 +84,12 @@ impl PathApproxParams {
         Self { increment }
     }
 
+    /// Return the increment
     pub fn increment(&self) -> Scalar {
         self.increment
     }
 
+    /// Generate points to approximate the circle within the boundary
     pub fn points(
         &self,
         boundary: impl Into<CurveBoundary<Point<1>>>,

crates/fj-core/src/algorithms/approx/curve.rs

@@ -76,6 +76,8 @@ fn approx_circle_on_straight_surface(
     surface: &SurfaceGeom,
     tolerance: impl Into<Tolerance>,
 ) -> Vec<ApproxPoint<1>> {
+    let tolerance = tolerance.into();
+
     approx_circle(circle, boundary, tolerance)
         .into_iter()
         .map(|(point_curve, point_surface)| {
@@ -93,7 +95,8 @@ fn approx_circle_on_straight_surface(
             //    point available, so it needs to be computed later anyway, in
             //    the general case.
 
-            let point_global = surface.point_from_surface_coords(point_surface);
+            let point_global =
+                surface.point_from_surface_coords(point_surface, tolerance);
             ApproxPoint::new(point_curve, point_global)
         })
         .collect()
@@ -105,6 +108,8 @@ fn approx_line_on_any_surface(
     surface: &SurfaceGeom,
     tolerance: impl Into<Tolerance>,
 ) -> Vec<ApproxPoint<1>> {
+    let tolerance = tolerance.into();
+
     let range_u = CurveBoundary::from(
         boundary
             .inner
@@ -121,7 +126,8 @@ fn approx_line_on_any_surface(
     for (u, _) in approx_u {
         let t = (u.t - line.origin().u) / line.direction().u;
         let point_surface = line.point_from_line_coords([t]);
-        let point_global = surface.point_from_surface_coords(point_surface);
+        let point_global =
+            surface.point_from_surface_coords(point_surface, tolerance);
         points.push(ApproxPoint::new(u, point_global));
     }
 
@@ -258,7 +264,7 @@ mod tests {
                     .layers
                     .geometry
                     .of_surface(&surface)
-                    .point_from_surface_coords(point_surface);
+                    .point_from_surface_coords(point_surface, tolerance);
                 ApproxPoint::new(point_local, point_global)
             })
             .collect::<Vec<_>>();
@@ -286,7 +292,7 @@ mod tests {
                     .layers
                     .geometry
                     .of_surface(&surface)
-                    .point_from_surface_coords(point_surface);
+                    .point_from_surface_coords(point_surface, tolerance);
                 ApproxPoint::new(point_local, point_global)
             })
             .collect::<Vec<_>>();

crates/fj-core/src/algorithms/approx/cycle.rs

@@ -53,6 +53,7 @@ pub fn approx_cycle(
                 half_edge.curve(),
                 surface,
                 start_position_curve,
+                tolerance,
                 &mut cache.vertex,
                 geometry,
             );

crates/fj-core/src/algorithms/approx/mod.rs

@@ -25,7 +25,10 @@ use vertex::VertexApproxCache;
 
 use crate::geometry::Geometry;
 
-pub use self::tolerance::{InvalidTolerance, Tolerance};
+pub use self::{
+    circle::PathApproxParams,
+    tolerance::{InvalidTolerance, Tolerance},
+};
 
 /// Approximate an object
 pub trait Approx: Sized {

crates/fj-core/src/algorithms/approx/vertex.rs

@@ -8,7 +8,7 @@ use crate::{
     topology::{Curve, Surface, Vertex},
 };
 
-use super::ApproxPoint;
+use super::{ApproxPoint, Tolerance};
 
 /// # Approximate a vertex position
 pub fn approx_vertex(
@@ -16,6 +16,7 @@ pub fn approx_vertex(
     curve: &Handle<Curve>,
     surface: &Handle<Surface>,
     position_curve: Point<1>,
+    tolerance: impl Into<Tolerance>,
     cache: &mut VertexApproxCache,
     geometry: &Geometry,
 ) -> ApproxPoint<1> {
@@ -32,7 +33,7 @@ pub fn approx_vertex(
         None => {
             let position_global = geometry
                 .of_surface(surface)
-                .point_from_surface_coords(position_surface);
+                .point_from_surface_coords(position_surface, tolerance);
             cache.insert(vertex, position_global)
         }
     };

crates/fj-core/src/algorithms/bounding_volume/face.rs

@@ -1,8 +1,9 @@
 use std::ops::Deref;
 
-use fj_math::Aabb;
+use fj_math::{Aabb, Vector};
 
 use crate::{
+    algorithms::approx::Tolerance,
     geometry::{Geometry, GlobalPath, SurfaceGeom},
     topology::Face,
 };
@@ -29,10 +30,28 @@ impl super::BoundingVolume<3> for &Face {
 
                         aabb_bottom.merged(&aabb_top)
                     }
-                    GlobalPath::Line(_) => Aabb {
-                        min: surface.point_from_surface_coords(aabb2.min),
-                        max: surface.point_from_surface_coords(aabb2.max),
-                    },
+                    GlobalPath::Line(_) => {
+                        // A bounding volume must include the body it bounds,
+                        // but does not need to match it precisely. So it's
+                        // okay, if it's a bit larger.
+                        //
+                        // Let's just choose a reasonable tolerance value here,
+                        // then make sure we enlarge the AABB accordingly, to
+                        // make sure it fits.
+                        let tolerance_f64 = 0.001;
+                        let tolerance = Tolerance::from_scalar(tolerance_f64)
+                            .expect("Tolerance provided is larger than zero");
+                        let offset = Vector::from([tolerance_f64; 3]);
+
+                        Aabb {
+                            min: surface.point_from_surface_coords(
+                                aabb2.min, tolerance,
+                            ) - offset,
+                            max: surface.point_from_surface_coords(
+                                aabb2.max, tolerance,
+                            ) + offset,
+                        }
+                    }
                 }
             })
     }

crates/fj-core/src/algorithms/triangulate/mod.rs

@@ -180,32 +180,32 @@ mod tests {
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(a);
+            .point_from_surface_coords(a, core.tolerance());
         let b = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(b);
+            .point_from_surface_coords(b, core.tolerance());
         let e = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(e);
+            .point_from_surface_coords(e, core.tolerance());
         let f = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(f);
+            .point_from_surface_coords(f, core.tolerance());
         let g = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(g);
+            .point_from_surface_coords(g, core.tolerance());
         let h = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(h);
+            .point_from_surface_coords(h, core.tolerance());
 
         // Let's test that some correct triangles are present. We don't need to
         // test them all.
@@ -275,27 +275,27 @@ mod tests {
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(a);
+            .point_from_surface_coords(a, core.tolerance());
         let b = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(b);
+            .point_from_surface_coords(b, core.tolerance());
         let c = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(c);
+            .point_from_surface_coords(c, core.tolerance());
         let d = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(d);
+            .point_from_surface_coords(d, core.tolerance());
         let e = core
             .layers
             .geometry
             .of_surface(&surface)
-            .point_from_surface_coords(e);
+            .point_from_surface_coords(e, core.tolerance());
 
         assert!(triangles.contains_triangle([a, b, d]));
         assert!(triangles.contains_triangle([a, d, e]));

crates/fj-core/src/geometry/surface.rs

@@ -1,6 +1,8 @@
 //! The geometry that defines a surface
 
-use fj_math::{Line, Point, Transform, Vector};
+use fj_math::{Point, Scalar, Transform, Triangle, Vector};
+
+use crate::algorithms::approx::{PathApproxParams, Tolerance};
 
 use super::GlobalPath;
 
@@ -26,31 +28,121 @@ pub enum SurfaceGeom {
 }
 
 impl SurfaceGeom {
+    /// # Access the origin of the surface
+    pub fn origin(&self) -> Point<3> {
+        let Self::Basic { u, .. } = self;
+        match u {
+            GlobalPath::Circle(circle) => circle.center(),
+            GlobalPath::Line(line) => line.origin(),
+        }
+    }
+
+    /// # Return the triangle at the provided point on the surface
+    ///
+    /// Select a triangle of the surface's triangle mesh representation, the one
+    /// at the provided surface point. Return that triangle, as well as the
+    /// barycentric coordinates of the provided point on the triangle.
+    ///
+    /// ## Triangle Size and Validity
+    ///
+    /// If a surface is curved along both axes, the triangle's size is chosen
+    /// such, that it approximates the surface, with the maximum allowed
+    /// deviation of the actual surface defined by the provided tolerance
+    /// argument.
+    ///
+    /// Otherwise, the size of the returned triangle is at least partially
+    /// arbitrary. Take the extreme case of a plane: Since it is not curved at
+    /// all, the returned triangle can be arbitrarily large.
+    ///
+    /// However, since surfaces are infinite, and we can't represent infinite
+    /// triangles, there is no sensible upper bound for the size. Instead, to
+    /// prevent an arbitrary choice for the size of triangles, which would imply
+    /// properties of the surface that are not true, and might therefore be
+    /// confusing, the triangles returned by this function have a length of zero
+    /// along axes that do not require approximation.
+    ///
+    /// The most extreme case would be a plane, for which the returned triangle
+    /// is collapsed to a point. For a cylinder, the triangle would have the
+    /// appropriate width to approximate the curved axis given the provided
+    /// tolerance, while having zero height.
+    ///
+    /// ## Implementation Note
+    ///
+    /// At the time this was written, there was no dedicated type to represent
+    /// barycentric coordinates. Nor any other code that used them, I think.
+    ///
+    /// If this changes, and a special type for barycentric coordinates is
+    /// added, it would make sense to return that here.
+    pub fn triangle_at(
+        &self,
+        point_surface: impl Into<Point<2>>,
+        tolerance: impl Into<Tolerance>,
+    ) -> (Triangle<3>, [Scalar; 3]) {
+        let point_surface = point_surface.into();
+
+        let Self::Basic { u, v } = self;
+        match u {
+            GlobalPath::Circle(circle) => {
+                let params = PathApproxParams::for_circle(circle, tolerance);
+
+                let a = point_surface.u - params.increment();
+                let b = point_surface.u + params.increment();
+                let c = a; // triangle is degenerate, as per function docs
+
+                let triangle_points_in_circle_space = [a, b, c];
+                let triangle_points_in_global_space =
+                    triangle_points_in_circle_space
+                        .map(|point_circle| {
+                            circle.point_from_circle_coords([point_circle])
+                        })
+                        .map(|point_global| {
+                            point_global + *v * point_surface.v
+                        });
+
+                let triangle = Triangle::from(triangle_points_in_global_space);
+                let barycentric_coords = [0.5, 0.5, 0.0].map(Into::into);
+
+                (triangle, barycentric_coords)
+            }
+            GlobalPath::Line(line) => {
+                let a = line.direction();
+                let b = *v;
+
+                let point_global =
+                    line.origin() + a * point_surface.u + b * point_surface.v;
+
+                // We don't need to approximate a plane, so our triangle can be
+                // arbitrarily large or small. Here we choose the smallest
+                // possible size (it is collapsed to a point), as per the
+                // documentation of this function.
+                let triangle = Triangle::from([point_global; 3]);
+                let barycentric_coords = [1. / 3.; 3].map(Into::into);
+
+                (triangle, barycentric_coords)
+            }
+        }
+    }
+
     /// Convert a point in surface coordinates to model coordinates
     pub fn point_from_surface_coords(
         &self,
         point: impl Into<Point<2>>,
+        tolerance: impl Into<Tolerance>,
     ) -> Point<3> {
-        let point = point.into();
-        let Self::Basic { u, .. } = self;
-        u.point_from_path_coords([point.u])
-            + self.path_to_line().vector_from_line_coords([point.v])
+        let (triangle, barycentric_coords) = self.triangle_at(point, tolerance);
+        triangle.point_from_barycentric_coords(barycentric_coords)
     }
 
     /// Convert a vector in surface coordinates to model coordinates
     pub fn vector_from_surface_coords(
         &self,
         vector: impl Into<Vector<2>>,
+        tolerance: impl Into<Tolerance>,
     ) -> Vector<3> {
         let vector = vector.into();
-        let Self::Basic { u, .. } = self;
-        u.vector_from_path_coords([vector.u])
-            + self.path_to_line().vector_from_line_coords([vector.v])
-    }
-
-    fn path_to_line(&self) -> Line<3> {
-        let Self::Basic { u, v } = self;
-        Line::from_origin_and_direction(u.origin(), *v)
+        let point =
+            self.point_from_surface_coords(Point { coords: vector }, tolerance);
+        point - self.origin()
     }
 
     /// Transform the surface geometry
@@ -69,7 +161,10 @@ mod tests {
     use fj_math::{Line, Point, Vector};
     use pretty_assertions::assert_eq;
 
-    use crate::geometry::{GlobalPath, SurfaceGeom};
+    use crate::{
+        algorithms::approx::Tolerance,
+        geometry::{GlobalPath, SurfaceGeom},
+    };
 
     #[test]
     fn point_from_surface_coords() {
@@ -81,8 +176,11 @@ mod tests {
             v: Vector::from([0., 0., 2.]),
         };
 
+        // Value doesn't matter; we're dealing with a plane.
+        let tolerance = Tolerance::from_scalar(1.).unwrap();
+
         assert_eq!(
-            surface.point_from_surface_coords([2., 4.]),
+            surface.point_from_surface_coords([2., 4.], tolerance),
             Point::from([1., 5., 9.]),
         );
     }
@@ -97,8 +195,11 @@ mod tests {
             v: Vector::from([0., 0., 2.]),
         };
 
+        // Value doesn't matter; we're dealing with a plane.
+        let tolerance = Tolerance::from_scalar(1.).unwrap();
+
         assert_eq!(
-            surface.vector_from_surface_coords([2., 4.]),
+            surface.vector_from_surface_coords([2., 4.], tolerance),
             Vector::from([0., 4., 8.]),
         );
     }

crates/fj-core/src/operations/holes.rs

@@ -93,7 +93,10 @@ impl AddHole for Shell {
                 core.layers
                     .geometry
                     .of_surface(location.face.surface())
-                    .point_from_surface_coords(location.position)
+                    .point_from_surface_coords(
+                        location.position,
+                        core.tolerance(),
+                    )
             };
 
             let entry_point = point(&entry_location);

crates/fj-core/src/operations/sweep/path.rs

@@ -66,18 +66,26 @@ impl SweepSurfacePath for SurfacePath {
 
         let u = match self {
             SurfacePath::Circle(circle) => {
-                let center = surface.point_from_surface_coords(circle.center());
-                let a = surface.vector_from_surface_coords(circle.a());
-                let b = surface.vector_from_surface_coords(circle.b());
+                let center = surface.point_from_surface_coords(
+                    circle.center(),
+                    core.tolerance(),
+                );
+                let a = surface
+                    .vector_from_surface_coords(circle.a(), core.tolerance());
+                let b = surface
+                    .vector_from_surface_coords(circle.b(), core.tolerance());
 
                 let circle = Circle::new(center, a, b);
 
                 GlobalPath::Circle(circle)
             }
             SurfacePath::Line(line) => {
-                let origin = surface.point_from_surface_coords(line.origin());
-                let direction =
-                    surface.vector_from_surface_coords(line.direction());
+                let origin = surface
+                    .point_from_surface_coords(line.origin(), core.tolerance());
+                let direction = surface.vector_from_surface_coords(
+                    line.direction(),
+                    core.tolerance(),
+                );
 
                 let line = Line::from_origin_and_direction(origin, direction);
 

crates/fj-core/src/validate/solid.rs

@@ -125,6 +125,7 @@ impl SolidValidationError {
                                 .unwrap()
                                 .position,
                         ),
+                        config.tolerance,
                     ),
                     h.start_vertex().clone(),
                 ))

crates/fj-core/src/validation/checks/coincident_half_edges_are_not_siblings.rs

@@ -3,6 +3,7 @@ use std::fmt;
 use fj_math::{Point, Scalar};
 
 use crate::{
+    algorithms::approx::Tolerance,
     geometry::{CurveBoundary, Geometry},
     queries::{
         AllHalfEdgesWithSurface, BoundingVerticesOfHalfEdge, CycleOfHalfEdge,
@@ -117,24 +118,29 @@ impl ValidationCheck<Shell> for CoincidentHalfEdgesAreNotSiblings {
                 }
 
                 let Some(points_and_distances) = distances(
-                    half_edge_a.clone(),
-                    object
-                        .find_cycle_of_half_edge(half_edge_a)
-                        .unwrap()
-                        .half_edges()
-                        .after(half_edge_a)
-                        .unwrap()
-                        .start_vertex(),
-                    surface_a,
-                    half_edge_b.clone(),
-                    object
-                        .find_cycle_of_half_edge(half_edge_b)
-                        .unwrap()
-                        .half_edges()
-                        .after(half_edge_b)
-                        .unwrap()
-                        .start_vertex(),
-                    surface_b,
+                    (
+                        half_edge_a.clone(),
+                        object
+                            .find_cycle_of_half_edge(half_edge_a)
+                            .unwrap()
+                            .half_edges()
+                            .after(half_edge_a)
+                            .unwrap()
+                            .start_vertex(),
+                        surface_a,
+                    ),
+                    (
+                        half_edge_b.clone(),
+                        object
+                            .find_cycle_of_half_edge(half_edge_b)
+                            .unwrap()
+                            .half_edges()
+                            .after(half_edge_b)
+                            .unwrap()
+                            .start_vertex(),
+                        surface_b,
+                    ),
+                    config.tolerance,
                     geometry,
                 ) else {
                     // The geometry to compute the distances is not available,
@@ -179,12 +185,17 @@ impl ValidationCheck<Shell> for CoincidentHalfEdgesAreNotSiblings {
 ///
 /// Returns an [`Iterator`] of the distance at each sample.
 fn distances(
-    half_edge_a: Handle<HalfEdge>,
-    end_vertex_a: &Handle<Vertex>,
-    surface_a: &Handle<Surface>,
-    half_edge_b: Handle<HalfEdge>,
-    end_vertex_b: &Handle<Vertex>,
-    surface_b: &Handle<Surface>,
+    (half_edge_a, end_vertex_a, surface_a): (
+        Handle<HalfEdge>,
+        &Handle<Vertex>,
+        &Handle<Surface>,
+    ),
+    (half_edge_b, end_vertex_b, surface_b): (
+        Handle<HalfEdge>,
+        &Handle<Vertex>,
+        &Handle<Surface>,
+    ),
+    tolerance: Tolerance,
     geometry: &Geometry,
 ) -> Option<Vec<([Point<3>; 2], Scalar)>> {
     fn sample(
@@ -192,6 +203,7 @@ fn distances(
         half_edge: &Handle<HalfEdge>,
         end_vertex: &Handle<Vertex>,
         surface: &Handle<Surface>,
+        tolerance: Tolerance,
         geometry: &Geometry,
     ) -> Option<Point<3>> {
         let [start, end] = [
@@ -217,7 +229,7 @@ fn distances(
         Some(
             geometry
                 .of_surface(surface)
-                .point_from_surface_coords(surface_coords),
+                .point_from_surface_coords(surface_coords, tolerance),
         )
     }
 
@@ -230,13 +242,20 @@ fn distances(
     let mut distances = Vec::new();
     for i in 0..sample_count {
         let percent = i as f64 * step;
-        let sample1 =
-            sample(percent, &half_edge_a, end_vertex_a, surface_a, geometry)?;
+        let sample1 = sample(
+            percent,
+            &half_edge_a,
+            end_vertex_a,
+            surface_a,
+            tolerance,
+            geometry,
+        )?;
         let sample2 = sample(
             1.0 - percent,
             &half_edge_b,
             end_vertex_b,
             surface_b,
+            tolerance,
             geometry,
         )?;
         distances.push(([sample1, sample2], sample1.distance_to(&sample2)))

crates/fj-core/src/validation/checks/curve_geometry_mismatch.rs

@@ -162,10 +162,16 @@ impl ValidationCheck<Shell> for CurveGeometryMismatch {
                             .path
                             .point_from_path_coords(point_curve);
 
-                        let a_global =
-                            surface_geom_a.point_from_surface_coords(a_surface);
-                        let b_global =
-                            surface_geom_b.point_from_surface_coords(b_surface);
+                        let a_global = surface_geom_a
+                            .point_from_surface_coords(
+                                a_surface,
+                                config.tolerance,
+                            );
+                        let b_global = surface_geom_b
+                            .point_from_surface_coords(
+                                b_surface,
+                                config.tolerance,
+                            );
 
                         let distance = (a_global - b_global).magnitude();