mesh: smoothing lifts elements up to 0.7 and visits mirror twins alike (#402)

The compass search after the sweeps only lifted elements below 0.5 (the quality floor), so a
mesh whose worst elements sat just above it stayed there: every corner of an enclosure base
stopped at a scaled Jacobian of 0.50-0.52 where its walls meet the corner bosses. It now lifts
every element below 0.7, never lowering another below the lower of its own value and 0.5.

The sweeps and the search are Gauss-Seidel, so the node order shaped the result and mirror
image corners smoothed differently. Nodes are now visited by a key that a mirror or a quarter
turn about the nodes' centroid keeps (distance from it, then the sorted offsets along the
axes). What spread is left between mirrored corners (0.686 vs 0.702 on the test base) is most likely the projection onto the body's
tessellation, which the kernel does not make mirror symmetric (not yet proven).

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

crates/mesh/src/project/smooth.rs

@@ -43,6 +43,10 @@ const SWEEP_FLOOR: f64 = 0.3;
 /// Scaled Jacobians closer than this count as the same when two grids are compared.
 const SAME_QUALITY: f64 = 1e-6;
 
+/// The compass search lifts every element below this scaled Jacobian, as far as its nodes can
+/// move without lowering another below the lower of its own value and [`QUALITY_FLOOR`].
+const LIFT_BELOW: f64 = 0.7;
+
 /// Rounds of the compass search on the worst elements.
 const OPTIMIZE_ROUNDS: usize = 20;
 
@@ -339,6 +343,27 @@ pub(crate) fn smooth(m: Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Sm
     report
 }
 
+/// The order nodes are visited in, a sweep being Gauss–Seidel (each node sees the moves made
+/// before it): by distance from the nodes' centroid, then by height along each axis's absolute
+/// offset from it, so a mesh and its mirror image (or a quarter turn of it) about the centroid
+/// visit corresponding nodes alike and smooth alike. Keys are rounded so mirror twins tie; twins
+/// lie apart, so their order among themselves does not matter.
+fn sweep_order(x: &[P]) -> Vec<u32> {
+    let n = x.len().max(1) as f64;
+    let c = x.iter().fold([0.0; 3], |a, p| add(a, scale(*p, 1.0 / n)));
+    let span = x.iter().map(|p| dist(*p, c)).fold(0.0, f64::max).max(f64::MIN_POSITIVE);
+    let q = |v: f64| (v / span * 1e7).round() as i64;
+    let mut order: Vec<u32> = (0..x.len() as u32).collect();
+    order.sort_by_key(|&v| {
+        let p = at(x, v);
+        let d = sub(p, c);
+        let mut offs = [q(d[0].abs()), q(d[1].abs()), q(d[2].abs())];
+        offs.sort_unstable();
+        (q(norm(d)), offs, v)
+    });
+    order
+}
+
 /// One smoothing run; sweep moves keep every element at or above the lower of its scaled
 /// Jacobian before and `sweep_floor`.
 fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing, iterations: usize, tol: f64, sweep_floor: f64) -> SmoothReport {
@@ -395,11 +420,21 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing,
     if method == Smoothing::None {
         return report;
     }
+    let order = sweep_order(m.nodes);
+    let rank: Vec<usize> = {
+        let mut r = vec![0; n];
+        for (i, &v) in order.iter().enumerate() {
+            if let Some(x) = r.get_mut(v as usize) {
+                *x = i;
+            }
+        }
+        r
+    };
     let ptol = 1e-3 * tol;
     for it in 0..iterations {
         report.iterations = it + 1;
         let mut moved = 0.0f64;
-        for v in 0..n as u32 {
+        for &v in &order {
             let c = cons(v);
             let p = at(m.nodes, v);
             let mean = |nbrs: &mut dyn Iterator<Item = u32>| -> Option<P> {
@@ -461,8 +496,10 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing,
     // IJNME 40:3979–4002, 1997, on Knupp's hexahedral metric), boundary nodes kept on their
     // targets.
     for _ in 0..OPTIMIZE_ROUNDS {
-        let poor: BTreeSet<u32> =
-            m.hexes.iter().filter(|h| hex_scaled_jacobian(&h.map(|q| at(m.nodes, q))) < QUALITY_FLOOR).flat_map(|h| h.iter().copied()).collect();
+        let mut poor: Vec<u32> =
+            m.hexes.iter().filter(|h| hex_scaled_jacobian(&h.map(|q| at(m.nodes, q))) < LIFT_BELOW).flat_map(|h| h.iter().copied()).collect();
+        poor.sort_unstable_by_key(|&v| rank.get(v as usize).copied().unwrap_or(usize::MAX));
+        poor.dedup();
         if poor.is_empty() {
             break;
         }