mesh: projection smoother keeps the best sweep grid and lifts from the better start (#489)

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

crates/mesh/src/project/smooth.rs

@@ -319,6 +319,33 @@ fn mesh_quality(x: &[P], hexes: &[[u32; 8]]) -> Local {
     l
 }
 
+/// A grid's smallest and mean scaled Jacobian, to pick the best grid the sweeps passed through.
+#[derive(Debug, Clone, Copy)]
+struct SweepScore {
+    min: f64,
+    mean: f64,
+}
+
+impl SweepScore {
+    /// A larger smallest scaled Jacobian, or the same one (to rounding) with a larger mean.
+    fn better_than(self, other: SweepScore) -> bool {
+        if (self.min - other.min).abs() > SCORE_TIE { self.min > other.min } else { self.mean > other.mean + SCORE_TIE }
+    }
+}
+
+/// Two smallest scaled Jacobians this close count as equal.
+const SCORE_TIE: f64 = 1e-9;
+
+fn sweep_score(x: &[P], hexes: &[[u32; 8]]) -> SweepScore {
+    let (mut min, mut sum) = (f64::INFINITY, 0.0);
+    for hex in hexes {
+        let sj = hex_scaled_jacobian(&hex.map(|q| at(x, q)));
+        min = min.min(sj);
+        sum += sj;
+    }
+    SweepScore { min, mean: if hexes.is_empty() { 0.0 } else { sum / hexes.len() as f64 } }
+}
+
 /// Smooth in place: at most `iterations` sweeps, stopping when no node moves more than `tol`.
 ///
 /// The sweeps let a move take an element down to [`SWEEP_FLOOR`] (a sheared block straightens
@@ -403,6 +430,10 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing,
         return report;
     }
     let ptol = 1e-3 * tol;
+    // The sweeps do not raise the smallest scaled Jacobian monotonically (a move may take an
+    // element down to the sweep floor), so where they stop decides the result; keep the best grid
+    // seen (smallest scaled Jacobian, then the mean) and lift the worst elements from that.
+    let mut best = (sweep_score(m.nodes, hexes), m.nodes.clone());
     for it in 0..iterations {
         report.iterations = it + 1;
         let mut moved = 0.0f64;
@@ -458,10 +489,25 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing,
                 }
             }
         }
+        let score = sweep_score(m.nodes, hexes);
+        if score.better_than(best.0) {
+            best = (score, m.nodes.clone());
+        }
         if moved <= tol {
             break;
         }
     }
+    // Where the sweeps stopped and the best grid they passed through each go through the lift
+    // and the straightening; the better result is kept (the lift can raise an element the best
+    // sweep grid left low, so neither start always wins).
+    let end = m.nodes.clone();
+    let mut starts = vec![end];
+    if best.1 != starts[0] {
+        starts.push(best.1);
+    }
+    let mut finished: Option<(SweepScore, Local, Vec<P>)> = None;
+    for start in starts {
+        *m.nodes = start;
         // The worst elements: each of their nodes searched for the place that raises the smallest
         // scaled Jacobian around it (a compass search, the optimization-based smoothing of L. A.
         // Freitag and C. Ollivier-Gooch, *Tetrahedral mesh improvement using swapping and smoothing*,
@@ -560,5 +606,14 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing,
                 break;
             }
         }
+        let (score, local) = (sweep_score(m.nodes, hexes), mesh_quality(m.nodes, hexes));
+        let better = finished.as_ref().is_none_or(|(s, l, _)| score.better_than(*s) && local.ratio <= l.ratio.max(EDGE_RATIO_CEILING));
+        if better {
+            finished = Some((score, local, m.nodes.clone()));
+        }
+    }
+    if let Some((_, _, nodes)) = finished {
+        *m.nodes = nodes;
+    }
     report
 }