mesh: split a projection block through its whole sheet; name twisted elements' block junction (#399)

Topology::split cuts a block at given fractions and every block chained to it
through parallel edges, orienting the fractions through the chain, so the
topology stays conforming. A certified-invalid element is reported with its
block, cell, centre and the block vertex and blocks meeting there; the mesh is
never emitted.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

crates/mesh/src/project/mod.rs

@@ -516,12 +516,12 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
 
     // Certified quality, the gap to the geometry and the tags.
     let (mut min_sj, mut sum_sj, mut volume) = (f64::INFINITY, 0.0, 0.0);
-    let mut invalid = 0usize;
-    for h in &lat.hexes {
+    let mut invalid: Vec<usize> = Vec::new();
+    for (e, h) in lat.hexes.iter().enumerate() {
         let q = element_quality(ElementKind::Hex8, &corners_of(&nodes, h))?;
         let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0);
         if q.validity().is_none_or(|v| v != crate::quality::Validity::Valid) {
-            invalid += 1;
+            invalid.push(e);
         }
         min_sj = min_sj.min(sj);
         sum_sj += sj;
@@ -568,11 +568,12 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
                 format!("worst in block {} near ({:.3}, {:.3}, {:.3})", c.0, x[0], x[1], x[2])
             })
     };
-    if invalid > 0 {
+    if !invalid.is_empty() {
         return Err(Error::Failed(format!(
-            "{invalid} of {} elements are inverted or degenerate after projection (smallest scaled Jacobian {min_sj:.3}, {}): the blocks do not fit the part — add an O-grid or a template where a block corner lands on a curved face",
+            "{} of {} elements are inverted, twisted or degenerate after projection (certified Jacobian; smallest corner scaled Jacobian {min_sj:.3}): {}. The blocks do not fit the part there — add an O-grid or a template where a block corner lands on a curved face, or split the blocks so the faces meeting at that junction are not squeezed",
+            invalid.len(),
             lat.hexes.len(),
-            worst()
+            where_invalid(t, &lat, &nodes, &invalid)
         )));
     }
     if min_sj < opts.min_quality {
@@ -595,6 +596,36 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o
     Ok(ProjectedMesh { nodes, hexes: lat.hexes, boundary, report })
 }
 
+/// Most invalid elements an error names.
+const NAMED_INVALID: usize = 4;
+
+/// Where the invalid elements `invalid` are: per element its block, cell and centre, and when it
+/// touches a block vertex (a junction), the vertex and the blocks that meet there.
+fn where_invalid(t: &Topology, lat: &lattice::Lattice, nodes: &[P], invalid: &[usize]) -> String {
+    let mut parts: Vec<String> = Vec::new();
+    for &e in invalid.iter().take(NAMED_INVALID) {
+        let (Some(h), Some((b, cell))) = (lat.hexes.get(e), lat.hex_cell.get(e)) else { continue };
+        let x = scale(h.iter().fold([0.0; 3], |a, &v| add(a, nodes.get(v as usize).copied().unwrap_or_default())), 0.125);
+        let mut s = format!("block {b} cell {cell:?} near ({:.3}, {:.3}, {:.3})", x[0], x[1], x[2]);
+        let vertices: BTreeSet<u32> = h
+            .iter()
+            .filter_map(|&v| match lat.entity.get(v as usize) {
+                Some(Entity::Vertex(k)) => Some(*k),
+                _ => None,
+            })
+            .collect();
+        for k in vertices {
+            let around: Vec<String> = t.blocks.iter().enumerate().filter(|(_, blk)| blk.contains(&k)).map(|(i, _)| i.to_string()).collect();
+            s.push_str(&format!(" at block vertex {k} where blocks {} meet", around.join(", ")));
+        }
+        parts.push(s);
+    }
+    if invalid.len() > NAMED_INVALID {
+        parts.push(format!("and {} more", invalid.len() - NAMED_INVALID));
+    }
+    parts.join("; ")
+}
+
 /// The polyline's points redistributed evenly by arc length (ends kept).
 fn even(line: &[P]) -> Vec<P> {
     let m = line.len().saturating_sub(1);

crates/mesh/src/project/tests.rs

@@ -193,6 +193,38 @@ fn o_grid_keeps_neighbours_and_opens_faces() {
     t.check().unwrap();
 }
 
+/// A split runs through the whole sheet of blocks: here a 2×2 grid split across its first
+/// direction, and the reversed neighbour counts its fractions from the shared end.
+#[test]
+fn split_cuts_the_whole_sheet_and_stays_conforming() {
+    let mut t = Topology::grid([&[0.0, 1.0, 2.0], &[0.0, 1.0, 2.0], &[0.0, 1.0]]).unwrap();
+    // Block 1 relabelled so its direction 0 runs the other way (still right-handed).
+    t.blocks[1] = [t.blocks[1][2], t.blocks[1][3], t.blocks[1][0], t.blocks[1][1], t.blocks[1][6], t.blocks[1][7], t.blocks[1][4], t.blocks[1][5]];
+    t.check().unwrap();
+    // Splitting block 0 along Y cuts block 1 (its X neighbour shares the Y edges) but not 2 or 3.
+    let made = t.split(0, 1, &[0.25]).unwrap();
+    assert_eq!(made.len(), 2);
+    assert_eq!(t.blocks.len(), 6);
+    t.check().unwrap();
+    let uses = lattice::face_uses(&t).unwrap();
+    assert!(uses.values().all(|&u| u <= 2));
+    // Every new vertex is a quarter of the way up in Y; the cut faces are shared.
+    assert!(t.vertices.iter().any(|v| (v[1] - 0.25).abs() < 1e-12));
+    assert!(t.vertices.iter().all(|v| [0.0, 0.25, 1.0, 2.0].iter().any(|y| (v[1] - y).abs() < 1e-12)));
+    // Splitting along X goes through the whole grid's X sheet: four more blocks.
+    let made = t.split(2, 0, &[1.0 / 3.0, 2.0 / 3.0]).unwrap();
+    assert_eq!(t.blocks.len(), 6 + 3 * 2);
+    assert!(!made.is_empty());
+    t.check().unwrap();
+    assert!(lattice::face_uses(&t).unwrap().values().all(|&u| u <= 2));
+    let mut u = t.clone();
+    for bad in [vec![], vec![0.0], vec![0.5, 0.4], vec![1.2], vec![f64::NAN]] {
+        assert!(u.split(0, 0, &bad).is_err(), "{bad:?}");
+    }
+    assert!(u.split(0, 3, &[0.5]).is_err() && u.split(99, 0, &[0.5]).is_err());
+    assert_eq!(u, t);
+}
+
 #[test]
 fn hostile_topologies_are_errors() {
     let (s, g) = sphere(1.0, 32);

crates/mesh/src/project/topology.rs

@@ -31,6 +31,17 @@ pub(crate) fn face_number(d: usize, s: usize) -> u8 {
     FACE_SIDES.iter().position(|x| *x == (d, s)).map_or(1, |i| i as u8 + 1)
 }
 
+/// The four edges of a block along direction `d`, each from its side-0 corner to its side-1 corner.
+fn edges_along(block: &[u32; 8], d: usize) -> [[u32; 2]; 4] {
+    let mut out = [[0; 2]; 4];
+    for (slot, (k, bits)) in out.iter_mut().zip(CORNERS.iter().enumerate().filter(|(_, bits)| bits[d] == 0)) {
+        let mut e = *bits;
+        e[d] = 1;
+        *slot = [block[k], block[corner_at(e)]];
+    }
+    out
+}
+
 /// Blocks over shared vertices.
 #[derive(Debug, Clone, Default, PartialEq)]
 pub struct Topology {
@@ -322,6 +333,24 @@ impl Topology {
         Ok(best.1)
     }
 
+    /// The direction of block `b` (0 to 2) most nearly along `dir`, and whether it runs against it.
+    pub fn direction_along(&self, b: usize, dir: P) -> Result<(usize, bool)> {
+        let c = self.corners(b)?;
+        let mut best = (f64::NEG_INFINITY, 0, false);
+        for d in 0..3 {
+            let mut v = [0.0; 3];
+            for [p, q] in edges_along(&[0, 1, 2, 3, 4, 5, 6, 7], d) {
+                v = add(v, sub(c[q as usize], c[p as usize]));
+            }
+            let l = crate::vec3::norm(v);
+            let cos = if l > 0.0 { crate::vec3::dot(v, dir) / l } else { 0.0 };
+            if cos.abs() > best.0 {
+                best = (cos.abs(), d, cos < 0.0);
+            }
+        }
+        Ok((best.1, best.2))
+    }
+
     /// The corners of a block face in cyclic order.
     fn face_quad(&self, b: usize, f: u8) -> Result<[u32; 4]> {
         let (d, s) = FACE_SIDES.get(usize::from(f).wrapping_sub(1)).copied().ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face")))?;
@@ -391,6 +420,119 @@ impl Topology {
         Ok(made)
     }
 
+    /// Split block `b` along its direction `d` (0 to 2) at the fractions `at` (increasing, between
+    /// 0 and 1, measured from its side-0 face), and with it the whole sheet of blocks the cut runs
+    /// through: every block holding an edge parallel to the cut edges (chained through shared
+    /// edges) is split the same way, so the topology stays conforming. Each new vertex lies on its
+    /// edge at the fraction, counted from the end `b`'s side-0 face reaches through the chain. A
+    /// split block keeps its index for its first slab (from its own side 0); the other slabs are
+    /// appended. Returns the appended blocks.
+    pub fn split(&mut self, b: usize, d: usize, at: &[f64]) -> Result<Vec<usize>> {
+        if d > 2 {
+            return Err(Error::Invalid(format!("a block has directions 0 to 2, not {d}")));
+        }
+        if at.is_empty() || at.len() > 1000 || !at.iter().all(|t| *t > 0.0 && *t < 1.0) || at.windows(2).any(|w| !(w[1] > w[0])) {
+            return Err(Error::Invalid(format!("split fractions must be increasing numbers between 0 and 1, got {at:?}")));
+        }
+        let seed = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
+        let key = |[p, q]: [u32; 2]| [p.min(q), p.max(q)];
+        // The sheet: each cut edge (ascending key) with the vertex its fractions count from, and
+        // each block in it with its cut direction, grown through an edge → blocks index.
+        let mut by_edge: BTreeMap<[u32; 2], Vec<usize>> = BTreeMap::new();
+        for (i, block) in self.blocks.iter().enumerate() {
+            for dd in 0..3 {
+                for e in edges_along(block, dd) {
+                    by_edge.entry(key(e)).or_default().push(i);
+                }
+            }
+        }
+        let mut start: BTreeMap<[u32; 2], u32> = BTreeMap::new();
+        for e in edges_along(&seed, d) {
+            start.insert(key(e), e[0]);
+        }
+        let mut cut: BTreeMap<usize, usize> = BTreeMap::from([(b, d)]);
+        let mut queue = vec![(b, d)];
+        while let Some((i, dd)) = queue.pop() {
+            let block = self.blocks.get(i).copied().unwrap_or([0; 8]);
+            let along = edges_along(&block, dd);
+            // The block runs the sheet's way, or against it: read off an edge already in the
+            // sheet (the one it was reached through), and orient its other edges to match.
+            let against = along.iter().find_map(|e| start.get(&key(*e)).map(|s| *s != e[0])).unwrap_or(false);
+            for e in along {
+                let k = key(e);
+                start.entry(k).or_insert(if against { e[1] } else { e[0] });
+                for &n in by_edge.get(&k).map_or(&[][..], Vec::as_slice) {
+                    if cut.contains_key(&n) {
+                        continue;
+                    }
+                    let nb = self.blocks.get(n).copied().unwrap_or([0; 8]);
+                    if let Some(nd) = (0..3).find(|&x| edges_along(&nb, x).iter().any(|f| key(*f) == k)) {
+                        cut.insert(n, nd);
+                        queue.push((n, nd));
+                    }
+                }
+            }
+        }
+        // Each cut block's four edges must all count from its side-0 face or all from its side-1
+        // face; a sheet that comes back on itself reversed cannot take uneven fractions.
+        let mut flipped: BTreeMap<usize, bool> = BTreeMap::new();
+        for (&i, &dd) in &cut {
+            let block = self.blocks.get(i).copied().unwrap_or([0; 8]);
+            let ways = edges_along(&block, dd).map(|e| start.get(&key(e)) != Some(&e[0]));
+            if ways.iter().any(|w| *w != ways[0]) {
+                return Err(Error::Invalid(format!(
+                    "block {i} would be cut from both ends: the sheet through block {b} comes back on itself reversed; split another block"
+                )));
+            }
+            flipped.insert(i, ways[0]);
+        }
+        if self.blocks.len() + cut.len() * at.len() > MAX_BLOCKS {
+            return Err(Error::Limit(format!("the split makes more than {MAX_BLOCKS} blocks")));
+        }
+        // The vertices along each cut edge, from its start.
+        let mut rows: BTreeMap<[u32; 2], Vec<u32>> = BTreeMap::new();
+        for (k, &s) in &start {
+            let e = if k[0] == s { k[1] } else { k[0] };
+            let (Some(ps), Some(pe)) = (self.vertices.get(s as usize).copied(), self.vertices.get(e as usize).copied()) else {
+                return Err(Error::Invalid(format!("vertex {s} or {e} does not exist")));
+            };
+            let mut row = vec![s];
+            for &t in at {
+                row.push(self.vertices.len() as u32);
+                self.vertices.push(crate::vec3::lerp(ps, pe, t));
+            }
+            row.push(e);
+            rows.insert(*k, row);
+        }
+        let n = at.len() + 1;
+        let mut made = Vec::new();
+        for (&i, &dd) in &cut {
+            let block = self.blocks.get(i).copied().unwrap_or([0; 8]);
+            let rev = flipped.get(&i).copied().unwrap_or(false);
+            // Slab `j`, counted from the block's own side 0.
+            let slab = |j: usize| -> [u32; 8] {
+                let mut out = block;
+                for (k, bits) in CORNERS.iter().enumerate() {
+                    let (mut low, mut high) = (*bits, *bits);
+                    low[dd] = 0;
+                    high[dd] = 1;
+                    let row = rows.get(&key([block[corner_at(low)], block[corner_at(high)]])).map_or(&[][..], Vec::as_slice);
+                    let local = j + bits[dd];
+                    out[k] = row.get(if rev { n - local } else { local }).copied().unwrap_or(block[k]);
+                }
+                out
+            };
+            if let Some(first) = self.blocks.get_mut(i) {
+                *first = slab(0);
+            }
+            for j in 1..n {
+                made.push(self.blocks.len());
+                self.blocks.push(slab(j));
+            }
+        }
+        Ok(made)
+    }
+
     /// Remove blocks (the others keep their order and are renumbered).
     pub fn remove(&mut self, blocks: &[usize]) -> Result<()> {
         if let Some(b) = blocks.iter().find(|b| **b >= self.blocks.len()) {