diff --git a/crates/mesh/src/surface/tests.rs b/crates/mesh/src/surface/tests.rs index 8203a97..7841ee8 100644 --- a/crates/mesh/src/surface/tests.rs +++ b/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()); +} diff --git a/crates/mesh/src/surface/weld.rs b/crates/mesh/src/surface/weld.rs index b899d64..77673fd 100644 --- a/crates/mesh/src/surface/weld.rs +++ b/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, 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 = 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).