Surface: close T-junctions left by independently tessellated faces

A face boundary vertex lying (within the weld distance) on an undivided edge of
the neighbouring face splits that edge's triangle. Only edges no other triangle
uses are looked at, so closed surfaces are unchanged. The half of a crank arm
cut on its symmetry plane tessellates with such a junction where the cut face
meets the top face; it now meshes (tet and swept hex).

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

crates/mesh/src/surface/tests.rs

@@ -392,3 +392,23 @@ fn orient2d_sign(a: &Point3, b: &Point3, c: &Point3) -> i32 {
     let p = |q: &Point3| crate::predicates::Point2::new(q.coords()[0], q.coords()[1]).unwrap();
     crate::predicates::orient2d(&p(a), &p(b), &p(c)) as i32
 }
+
+/// A vertex of one face on an undivided edge of its neighbour (a T-junction, as independent
+/// face tessellations leave where a face is split) is closed by splitting that edge.
+#[test]
+fn t_junctions_are_closed() {
+    let (mut p, mut t, mut f) = cube([0.0; 3], 1.0);
+    // Edge 0–2 (x = 0, z = 0) gets a vertex at its middle on the x = 0 side only.
+    let m = p.len() as u32;
+    p.push([0.0, 0.5, 0.0]);
+    let i = t.iter().position(|x| *x == [0, 6, 2]).unwrap();
+    t[i] = [0, 6, m];
+    t.push([m, 6, 2]);
+    f.push(f[i]);
+    let (sp, st, sf) = soup(&(p, t, f));
+    let s = BoundarySurface::new(&sp, &st, &sf, &SurfaceOptions { weld_tolerance: 0.0, weld_distance: 1e-9, feature_angle: 30.0 }).unwrap();
+    assert!((s.volume() - 1.0).abs() < 1e-12, "{}", s.volume());
+    assert_eq!(s.triangles().len(), 14);
+    // Without a weld distance nothing is split, and the open edge is reported.
+    assert!(BoundarySurface::new(&sp, &st, &sf, &SurfaceOptions { weld_tolerance: 0.0, weld_distance: 0.0, feature_angle: 30.0 }).is_err());
+}

crates/mesh/src/surface/weld.rs

@@ -78,9 +78,95 @@ pub(crate) fn weld(positions: &[[f64; 3]], triangles: &[[u32; 3]], relative: f64
             source_triangle.push(index(t));
         }
     }
+    if tol > 0.0 {
+        split_t_junctions(&points, &mut out, &mut source_triangle, tol);
+    }
     Ok(Welded { points, triangles: out, source_vertex, source_triangle })
 }
 
+/// Rounds of T-junction splitting, and the most open edges looked at.
+const T_JUNCTION_ROUNDS: usize = 8;
+const MAX_OPEN_EDGES: usize = 100_000;
+
+/// Close T-junctions: where a vertex of one face's boundary lies on an edge of the neighbouring
+/// face (within `tol`) that the tessellation did not divide there, split that edge's triangle at
+/// it. Only edges no other triangle uses are looked at, so a closed surface is left as it is.
+/// The standard repair of tessellations whose faces are triangulated independently (Barequet and
+/// Kumar 1997, *Repairing CAD models*).
+fn split_t_junctions(points: &[[f64; 3]], triangles: &mut Vec<[u32; 3]>, source: &mut Vec<u32>, tol: f64) {
+    use std::collections::{BTreeMap, BTreeSet};
+    for _ in 0..T_JUNCTION_ROUNDS {
+        let mut used: BTreeMap<(u32, u32), u32> = BTreeMap::new();
+        for tri in triangles.iter() {
+            for k in 0..3 {
+                *used.entry((tri[k], tri[(k + 1) % 3])).or_default() += 1;
+            }
+        }
+        // Edges used once, in one direction only.
+        let open: BTreeSet<(u32, u32)> = used.iter().filter(|(e, n)| **n == 1 && !used.contains_key(&(e.1, e.0))).map(|(e, _)| *e).collect();
+        if open.is_empty() || open.len() > MAX_OPEN_EDGES {
+            return;
+        }
+        let ends: BTreeSet<u32> = open.iter().flat_map(|&(a, b)| [a, b]).collect();
+        let at = |v: u32| points.get(v as usize).copied();
+        let mut changed = false;
+        let mut next = Vec::with_capacity(triangles.len() + 8);
+        let mut next_source = Vec::with_capacity(triangles.len() + 8);
+        for (tri, &src) in triangles.iter().zip(source.iter()) {
+            let mut split = None;
+            for k in 0..3 {
+                let (a, b, c) = (tri[k], tri[(k + 1) % 3], tri[(k + 2) % 3]);
+                if !open.contains(&(a, b)) {
+                    continue;
+                }
+                let (Some(pa), Some(pb)) = (at(a), at(b)) else { continue };
+                let ab = sub(pb, pa);
+                let len2 = dot(ab, ab);
+                if !(len2 > 0.0) {
+                    continue;
+                }
+                let mut on: Vec<(f64, u32)> = ends
+                    .iter()
+                    .filter(|&&v| v != a && v != b && v != c)
+                    .filter_map(|&v| {
+                        let p = at(v)?;
+                        let t = dot(sub(p, pa), ab) / len2;
+                        let q = [pa[0] + ab[0] * t, pa[1] + ab[1] * t, pa[2] + ab[2] * t];
+                        let d = sub(p, q);
+                        let inside = t * len2.sqrt() > tol && (1.0 - t) * len2.sqrt() > tol;
+                        (inside && dot(d, d) <= tol * tol).then_some((t, v))
+                    })
+                    .collect();
+                if !on.is_empty() {
+                    on.sort_by(|x, y| x.0.total_cmp(&y.0).then(x.1.cmp(&y.1)));
+                    split = Some((a, b, c, on));
+                    break;
+                }
+            }
+            match split {
+                Some((a, b, c, on)) => {
+                    let mut prev = a;
+                    for (_, v) in on.iter().copied().chain(std::iter::once((1.0, b))) {
+                        next.push([prev, v, c]);
+                        next_source.push(src);
+                        prev = v;
+                    }
+                    changed = true;
+                }
+                None => {
+                    next.push(*tri);
+                    next_source.push(src);
+                }
+            }
+        }
+        *triangles = next;
+        *source = next_source;
+        if !changed {
+            return;
+        }
+    }
+}
+
 /// Representatives found so far.
 enum Lookup {
     /// Identical coordinates only (keyed by bits; `+ 0.0` turns -0 into +0).