mesh: tag each projected quad with the face it lies on, slivers included

project_trimmed accepted the exact surface's point up to a fifth of the
nearest triangle's median edge off the face: on a tall fillet piece
(triangles 30 mm long) that let a point on a sibling piece of the same
cylinder count as on it, so every quad round the corner took the first
piece's name. The slack is now the triangle's own chordal gap to the
surface, doubled.

A face narrower than the elements (the sliver a corner boss leaves of the
inner fillet) gets one row of the quads it runs through, from the
neighbour with the most quads left. Tessellation::nearest_face finds the
face nearest a point.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/mesh/src/project/mod.rs

@@ -136,6 +136,63 @@ pub struct ProjectedMesh {
     pub report: ProjectReport,
 }
 
+/// Give each target face no boundary quad's centre landed on (a face narrower than the elements:
+/// a fillet's sliver between a cut and the walls) one row of the quads it runs through (of a
+/// target holding the face, centre within half the quad's shortest side of it): one quad of
+/// those beside each other across the face (projecting onto it at nearly the same point), all
+/// along it. A quad keeps its face when it is that face's
+/// last, so no face loses its only quad to a neighbour.
+fn claim_narrow_faces(boundary: &mut [BoundaryFace], centres: &[(P, f64, usize)], targets: &[Vec<u32>], g: &dyn Geometry) {
+    let mut count: BTreeMap<u32, usize> = BTreeMap::new();
+    for b in boundary.iter() {
+        *count.entry(b.brep_face).or_insert(0) += 1;
+    }
+    let targeted: BTreeSet<u32> = targets.iter().flatten().copied().collect();
+    for f in targeted {
+        if count.contains_key(&f) {
+            continue;
+        }
+        // (distance to the face, quad, its centre, the centre's foot on the face, shortest side)
+        let mut near: Vec<(f64, usize, P, P, f64)> = centres
+            .iter()
+            .enumerate()
+            .filter(|(_, (_, _, tg))| targets.get(*tg).is_some_and(|set| set.contains(&f)))
+            .filter_map(|(i, (centre, shortest, _))| {
+                let foot = g.project(f, *centre)?;
+                let d = dist(foot, *centre);
+                (d <= 0.5 * shortest).then_some((d, i, *centre, foot, *shortest))
+            })
+            .collect();
+        near.sort_by(|a, b| a.0.total_cmp(&b.0));
+        // The rows across the face: quads beside each other that project onto it at nearly the
+        // same point.
+        let mut rows: Vec<Vec<(usize, P, P, f64)>> = Vec::new();
+        for (_, i, centre, foot, shortest) in near {
+            let across = |(_, c, q, s): &(usize, P, P, f64)| {
+                let apart = dist(*c, centre);
+                apart < 1.5 * s.max(shortest) && dist(*q, foot) < 0.5 * apart
+            };
+            match rows.iter_mut().find(|row| row.iter().any(across)) {
+                Some(row) => row.push((i, centre, foot, shortest)),
+                None => rows.push(vec![(i, centre, foot, shortest)]),
+            }
+        }
+        // One quad per row: from the neighbour with the most quads left (the nearest of equals),
+        // so neighbours on both sides give alike.
+        for row in rows {
+            let left = |i: usize| boundary.get(i).and_then(|b| count.get(&b.brep_face)).copied().unwrap_or(0);
+            let Some(&(i, ..)) = row.iter().max_by_key(|(i, ..)| left(*i)) else { continue };
+            let Some(b) = boundary.get_mut(i) else { continue };
+            let Some(left) = count.get_mut(&b.brep_face) else { continue };
+            if *left > 1 {
+                *left -= 1;
+                b.brep_face = f;
+                *count.entry(f).or_insert(0) += 1;
+            }
+        }
+    }
+}
+
 /// The face groups joined across smooth B-rep edges.
 fn smooth_groups(surface: &BoundarySurface, faces: usize) -> Vec<usize> {
     let mut parent: Vec<usize> = (0..faces).collect();
@@ -626,6 +683,7 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
         }
     }
     let mut boundary = Vec::with_capacity(quads.len());
+    let mut centres = Vec::with_capacity(quads.len());
     let mut enclosed = 0.0;
     for (q, e, f, tg) in &quads {
         let pts = q.map(|v| nodes.get(v as usize).copied().unwrap_or_default());
@@ -634,7 +692,10 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
         let set = targets.get(*tg).map_or(&[][..], Vec::as_slice);
         let brep_face = on_target(g, set, centre).map_or_else(|| set.first().copied().unwrap_or(0), |h| h.2);
         boundary.push(BoundaryFace { element: *e, face: *f, brep_face });
+        let shortest = (0..4).map(|i| dist(pts[i], pts[(i + 1) % 4])).fold(f64::INFINITY, f64::min);
+        centres.push((centre, shortest, *tg));
     }
+    claim_narrow_faces(&mut boundary, &centres, &targets, g);
     let surface_volume = surface.volume();
     let report = ProjectReport {
         blocks: t.blocks.len(),

crates/mesh/src/tet/geometry.rs

@@ -507,6 +507,19 @@ impl Tessellation {
         self.face_tree(face).nearest4(p, &accept, &dist).map(|r| ([r[0], r[1], r[2]], r[3]))
     }
 
+    /// The B-rep face whose triangles lie nearest to `p`.
+    pub fn nearest_face(&self, p: P) -> Option<u32> {
+        let dist = |t: u32| {
+            let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]);
+            let d = sub(closest_on_triangle(p, self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2])), p);
+            (dot(d, d), [0.0, 0.0, 0.0, f64::from(t)])
+        };
+        let r = self.tri_tree.nearest4(p, &|_| true, &dist)?;
+        // The triangle index travels through the f64 exactly (triangle counts are capped far
+        // below 2^53).
+        self.tri_face.get(r[3] as usize).copied()
+    }
+
     /// The tree over `face`'s triangles (the whole boundary's for a face out of range, which
     /// the queries' filter then finds nothing in).
     fn face_tree(&self, face: u32) -> &Tree {
@@ -528,38 +541,51 @@ impl Tessellation {
     }
 }
 
-/// How far off a plane face's triangles, in triangle sizes, a point of the plane may lie and
-/// still be on the face (rounding only: the triangles lie on the plane).
+/// How far off a face's triangles, in triangle sizes, a point of its surface may lie and still
+/// be on the face, beyond the triangles' own gap to the surface (rounding only).
 const PLANE_SLACK: f64 = 1e-6;
 
 impl Tessellation {
     /// The closest point to `p` on the trimmed face `face`: the exact surface's closest point
-    /// when it lies on the face (within a fifth of the local triangle size of the face's
-    /// triangles, or on them for a plane, whose triangles can span the whole face), else the closest point of the face's triangles. Unlike
-    /// [`BoundaryGeometry::project_face`] it never returns a point of the surface's extension
-    /// beyond the face's boundary, whatever the distance from `p`.
+    /// when it lies on the face (no further off the face's nearest triangle than twice that
+    /// triangle's own chordal gap to the surface, which is rounding only for a plane), else the
+    /// closest point of the face's triangles. Unlike [`BoundaryGeometry::project_face`] it never
+    /// returns a point of the surface's extension beyond the face's boundary, whatever the
+    /// distance from `p` — nor of the same surface's other faces: the pieces of one fillet split
+    /// by a cut share their cylinder, and a point on one piece is off the others.
     pub fn project_trimmed(&self, face: u32, p: P) -> Option<P> {
         let accept = |t: u32| self.tri_face.get(t as usize) == Some(&face);
-        let nearest = |x: P| {
+        let corners = |t: u32| {
+            let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]);
+            (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2]))
+        };
+        let near = {
             let dist = |t: u32| {
-                let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]);
-                let (a, b, c) = (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2]));
-                let q = closest_on_triangle(x, a, b, c);
-                let d = sub(q, x);
-                let mut e = [measure::dist(a, b), measure::dist(b, c), measure::dist(c, a)];
-                e.sort_by(f64::total_cmp);
-                (dot(d, d), [q[0], q[1], q[2], e[1]])
+                let (a, b, c) = corners(t);
+                let q = closest_on_triangle(p, a, b, c);
+                let d = sub(q, p);
+                (dot(d, d), [q[0], q[1], q[2], 0.0])
             };
-            self.face_tree(face).nearest4(x, &accept, &dist).map(|r| ([r[0], r[1], r[2]], r[3]))
+            let r = self.face_tree(face).nearest4(p, &accept, &dist)?;
+            [r[0], r[1], r[2]]
         };
-        let (near, _) = nearest(p)?;
         let Some(surface) = self.analytic(face) else { return Some(near) };
         let Some(q) = surface.project(p) else { return Some(near) };
-        let (back, size) = nearest(q)?;
-        // A plane's triangles lie on it, so a point of the plane off them is off the face; a
-        // curved face's triangles are chords, a fraction of their size off the surface.
-        let slack = if matches!(surface, AnalyticSurface::Plane { .. }) { PLANE_SLACK } else { 0.2 };
-        Some(if measure::dist(back, q) <= slack * size { q } else { near })
+        // The nearest triangle to the surface's point, with its slack: the largest gap between
+        // the triangle and the surface (at its edges' midpoints and centre: a chord's deepest
+        // points), doubled, and a millionth of its size for rounding.
+        let gap = |x: P| surface.project(x).map_or(0.0, |y| measure::dist(x, y));
+        let dist = |t: u32| {
+            let (a, b, c) = corners(t);
+            let o = closest_on_triangle(q, a, b, c);
+            let d = sub(o, q);
+            let mid = |u: P, v: P| scale(add(u, v), 0.5);
+            let deepest = [mid(a, b), mid(b, c), mid(c, a), scale(add(add(a, b), c), 1.0 / 3.0)].into_iter().map(gap).fold(0.0, f64::max);
+            let size = measure::dist(a, b).max(measure::dist(b, c)).max(measure::dist(c, a));
+            (dot(d, d), [o[0], o[1], o[2], 2.0 * deepest + PLANE_SLACK * size])
+        };
+        let r = self.face_tree(face).nearest4(q, &accept, &dist)?;
+        Some(if measure::dist([r[0], r[1], r[2]], q) <= r[3] { q } else { near })
     }
 
     /// The outward unit normal of B-rep face `face` near `p`: the exact surface's (oriented like