Fix triangulation of sharp, concave faces

crates/fj-kernel/src/algorithms/approx/face.rs

@@ -34,6 +34,7 @@ impl Approx for &FaceSet {
         let min_distance = ValidationConfig::default().distinct_min_distance;
         let mut all_points: BTreeSet<ApproxPoint<2>> = BTreeSet::new();
 
+        // Run some validation code on the approximation.
         for approx in &approx {
             let approx: &FaceApprox = approx;
 

crates/fj-kernel/src/algorithms/triangulate/delaunay.rs

@@ -1,17 +1,48 @@
+use std::collections::BTreeMap;
+
 use fj_math::{Point, Scalar, Triangle, Winding};
 use spade::HasPosition;
 
-use crate::objects::Handedness;
+use crate::{algorithms::approx::cycle::CycleApprox, objects::Handedness};
 
 /// Create a Delaunay triangulation of all points
 pub fn triangulate(
-    points: Vec<TriangulationPoint>,
+    cycles: impl IntoIterator<Item = CycleApprox>,
     coord_handedness: Handedness,
 ) -> Vec<[TriangulationPoint; 3]> {
     use spade::Triangulation as _;
 
-    let triangulation = spade::DelaunayTriangulation::<_>::bulk_load(points)
-        .expect("Inserted invalid values into triangulation");
+    let mut triangulation = spade::ConstrainedDelaunayTriangulation::<_>::new();
+
+    let mut points = BTreeMap::new();
+
+    for cycle_approx in cycles {
+        let mut handle_prev = None;
+
+        for point in cycle_approx.points() {
+            let handle = match points.get(&point) {
+                Some(handle) => *handle,
+                None => {
+                    let handle = triangulation
+                        .insert(TriangulationPoint {
+                            point_surface: point.local_form,
+                            point_global: point.global_form,
+                        })
+                        .expect("Inserted invalid point into triangulation");
+
+                    points.insert(point, handle);
+
+                    handle
+                }
+            };
+
+            if let Some(handle_prev) = handle_prev {
+                triangulation.add_constraint(handle_prev, handle);
+            }
+
+            handle_prev = Some(handle);
+        }
+    }
 
     let mut triangles = Vec::new();
     for triangle in triangulation.inner_faces() {

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

@@ -6,7 +6,7 @@ mod polygon;
 use fj_interop::mesh::Mesh;
 use fj_math::Point;
 
-use self::{delaunay::TriangulationPoint, polygon::Polygon};
+use self::polygon::Polygon;
 
 use super::approx::{face::FaceApprox, Approx, Tolerance};
 
@@ -44,27 +44,20 @@ where
 
 impl Triangulate for FaceApprox {
     fn triangulate_into_mesh(self, mesh: &mut Mesh<Point<3>>) {
-        let points: Vec<_> = self
-            .points()
-            .into_iter()
-            .map(|point| TriangulationPoint {
-                point_surface: point.local_form,
-                point_global: point.global_form,
-            })
-            .collect();
         let face_as_polygon = Polygon::new()
             .with_exterior(
                 self.exterior
                     .points()
                     .into_iter()
                     .map(|point| point.local_form),
             )
-            .with_interiors(self.interiors.into_iter().map(|interior| {
+            .with_interiors(self.interiors.iter().map(|interior| {
                 interior.points().into_iter().map(|point| point.local_form)
             }));
 
+        let cycles = [self.exterior].into_iter().chain(self.interiors);
         let mut triangles =
-            delaunay::triangulate(points, self.coord_handedness);
+            delaunay::triangulate(cycles, self.coord_handedness);
         triangles.retain(|triangle| {
             face_as_polygon
                 .contains_triangle(triangle.map(|point| point.point_surface))
@@ -176,29 +169,26 @@ mod tests {
         Ok(())
     }
 
-    #[ignore]
     #[test]
     fn sharp_concave_shape() -> anyhow::Result<()> {
         let objects = Objects::new();
 
-        //
-        //                c
-        //               /|
-        //   e          / |
-        //   |\       /   |
-        //   | |     /    |
-        //   | \   /      |
-        //   |  \ /       |
-        //   |   d        |
-        //   a ---------- b
-        //
+        //   e       c
+        //   |\     /|
+        //   \ \   / b
+        //    \ \ / /
+        //     \ d /
+        //      \a/
 
-        let a = [0., 0.];
-        let b = [0.4, 0.];
-        //let b = [0.5, 0.]; // test passes with this change
-        let c = [0.4, 1.0];
+        // Naive Delaunay triangulation will create a triangle (c, d, e), which
+        // is not part of the polygon. The higher-level triangulation will
+        // filter that out, but it will result in missing triangles.
+
+        let a = [0.1, 0.0];
+        let b = [0.2, 0.9];
+        let c = [0.2, 1.0];
         let d = [0.1, 0.1];
-        let e = [0., 0.8];
+        let e = [0.0, 1.0];
 
         let surface = objects.surfaces.xy_plane();
         let face = Face::partial()
@@ -209,17 +199,15 @@ mod tests {
 
         let triangles = triangulate(face)?;
 
-        let a3 = surface.geometry().point_from_surface_coords(a);
-        let b3 = surface.geometry().point_from_surface_coords(b);
-        let c3 = surface.geometry().point_from_surface_coords(c);
-        let d3 = surface.geometry().point_from_surface_coords(d);
-        let e3 = surface.geometry().point_from_surface_coords(e);
-
-        assert!(triangles.contains_triangle([a3, b3, d3]));
-        assert!(triangles.contains_triangle([b3, c3, d3]));
-        assert!(triangles.contains_triangle([a3, d3, e3]));
+        let a = surface.geometry().point_from_surface_coords(a);
+        let b = surface.geometry().point_from_surface_coords(b);
+        let c = surface.geometry().point_from_surface_coords(c);
+        let d = surface.geometry().point_from_surface_coords(d);
+        let e = surface.geometry().point_from_surface_coords(e);
 
-        assert!(!triangles.contains_triangle([b3, e3, d3]));
+        assert!(triangles.contains_triangle([a, b, d]));
+        assert!(triangles.contains_triangle([a, d, e]));
+        assert!(triangles.contains_triangle([b, c, d]));
 
         Ok(())
     }