diff --git a/crates/mesh/src/relax.rs b/crates/mesh/src/relax.rs index 6f725fd..1b9ae8e 100644 --- a/crates/mesh/src/relax.rs +++ b/crates/mesh/src/relax.rs @@ -274,6 +274,8 @@ struct Elem { nodes: Vec, /// A shell's normal as it was: its side, which a move must keep. normal: Option

, + /// Its size as it was ([`strong::size`]): the strong mode keeps sizes near it. + size0: f64, } /// The shell normal from its corners (`None` for solids). @@ -327,6 +329,10 @@ struct State<'a> { mids_of: BTreeMap>, skin: BTreeSet, geometry: &'a dyn Constraints, + /// Hold nodes on geometry edges (the strong mode keeps feature lines as they are). + pin_curves: bool, + /// The largest edge ratio the strong mode may make. + ratio_cap: f64, } impl State<'_> { @@ -341,6 +347,7 @@ impl State<'_> { /// How `n` may move: a skin node the caller left free is held. fn freedom(&self, n: u32) -> Freedom { match self.geometry.freedom(n) { + Some(Freedom::Curve(_)) if self.pin_curves => Freedom::Fixed, Some(f) => f, None if self.skin.contains(&n) => Freedom::Fixed, None => Freedom::Free, @@ -500,6 +507,68 @@ impl State<'_> { } } +impl<'a> State<'a> { + /// The elements with a corner among `movers`, their nodes' positions and adjacency; with + /// `pin_curves`, nodes on geometry edges are held. + fn new(model: &Model, movers: &BTreeSet, geometry: &'a dyn Constraints, pin_curves: bool) -> Result> { + let mut state = State { + pos: BTreeMap::new(), + elems: Vec::new(), + around: BTreeMap::new(), + neighbours: BTreeMap::new(), + mids: BTreeMap::new(), + mids_of: BTreeMap::new(), + skin: skin_nodes(model, movers), + geometry, + pin_curves, + ratio_cap: f64::INFINITY, + }; + for (&label, e) in model.elements() { + if !quality::is_measured(e.kind) || e.nodes.len() != e.kind.node_count() { + continue; + } + let corners = e.kind.corner_count(); + if !e.nodes.iter().take(corners).any(|n| movers.contains(n)) { + continue; + } + let mut x = Vec::with_capacity(e.nodes.len()); + for &n in &e.nodes { + let p = model.node(n).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("element {label} refers to missing node {n}")))?; + state.pos.insert(n, p); + x.push(p); + } + let index = state.elems.len(); + for &n in e.nodes.iter().take(corners) { + if movers.contains(&n) { + state.around.entry(n).or_default().push(index); + } + } + for edge in e.kind.edges() { + let (Some(&a), Some(&b)) = (edge.first().and_then(|&i| e.nodes.get(i)), edge.get(1).and_then(|&i| e.nodes.get(i))) else { continue }; + for (u, v) in [(a, b), (b, a)] { + if movers.contains(&u) { + state.neighbours.entry(u).or_default().insert(v); + } + } + if let Some(&m) = edge.get(2).and_then(|&i| e.nodes.get(i)) { + let (pa, pb, pm) = (x.get(edge[0]).copied(), x.get(edge[1]).copied(), x.get(edge[2]).copied()); + if let (Some(pa), Some(pb), Some(pm)) = (pa, pb, pm) { + state.mids.entry(m).or_insert((a, b, sub(pm, scale(add(pa, pb), 0.5)))); + for c in [a, b] { + let list = state.mids_of.entry(c).or_default(); + if !list.contains(&m) { + list.push(m); + } + } + } + } + } + state.elems.push(Elem { label, kind: e.kind, normal: shell_normal(e.kind, &x), size0: strong::size(e.kind, &x), nodes: e.nodes.clone() }); + } + Ok(state) + } +} + /// Quality over the region. fn region_quality(s: &State<'_>) -> RegionQuality { let mut r = RegionQuality { elements: s.elems.len(), min_quality: f64::INFINITY, max_edge_ratio: 0.0, invalid: 0, worst: Vec::new() }; @@ -534,58 +603,7 @@ pub fn relax(model: &Model, weights: &BTreeMap, iterations: usize, geo return Err(Error::Invalid(format!("node {n} has weight {w}; weights are in (0, 1]"))); } let movers: BTreeSet = weights.keys().copied().collect(); - let mut state = State { - pos: BTreeMap::new(), - elems: Vec::new(), - around: BTreeMap::new(), - neighbours: BTreeMap::new(), - mids: BTreeMap::new(), - mids_of: BTreeMap::new(), - skin: skin_nodes(model, &movers), - geometry, - }; - for (&label, e) in model.elements() { - if !quality::is_measured(e.kind) || e.nodes.len() != e.kind.node_count() { - continue; - } - let corners = e.kind.corner_count(); - if !e.nodes.iter().take(corners).any(|n| movers.contains(n)) { - continue; - } - let mut x = Vec::with_capacity(e.nodes.len()); - for &n in &e.nodes { - let p = model.node(n).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("element {label} refers to missing node {n}")))?; - state.pos.insert(n, p); - x.push(p); - } - let index = state.elems.len(); - for &n in e.nodes.iter().take(corners) { - if movers.contains(&n) { - state.around.entry(n).or_default().push(index); - } - } - for edge in e.kind.edges() { - let (Some(&a), Some(&b)) = (edge.first().and_then(|&i| e.nodes.get(i)), edge.get(1).and_then(|&i| e.nodes.get(i))) else { continue }; - for (u, v) in [(a, b), (b, a)] { - if movers.contains(&u) { - state.neighbours.entry(u).or_default().insert(v); - } - } - if let Some(&m) = edge.get(2).and_then(|&i| e.nodes.get(i)) { - let (pa, pb, pm) = (x.get(edge[0]).copied(), x.get(edge[1]).copied(), x.get(edge[2]).copied()); - if let (Some(pa), Some(pb), Some(pm)) = (pa, pb, pm) { - state.mids.entry(m).or_insert((a, b, sub(pm, scale(add(pa, pb), 0.5)))); - for c in [a, b] { - let list = state.mids_of.entry(c).or_default(); - if !list.contains(&m) { - list.push(m); - } - } - } - } - } - state.elems.push(Elem { label, kind: e.kind, normal: shell_normal(e.kind, &x), nodes: e.nodes.clone() }); - } + let mut state = State::new(model, &movers, geometry, false)?; let original = state.pos.clone(); let before = region_quality(&state); let mut done = 0; @@ -865,7 +883,7 @@ pub fn apply_field( let mut reverted = false; for (label, kind, nodes) in &touched { let x0 = at(model, nodes); - let e = Elem { label: *label, kind: *kind, nodes: nodes.clone(), normal: shell_normal(*kind, &x0) }; + let e = Elem { label: *label, kind: *kind, nodes: nodes.clone(), normal: shell_normal(*kind, &x0), size0: 0.0 }; if measure(&e, &x0).1 && !measure(&e, &at(&work, nodes)).1 { for n in nodes { if moved.remove(n).is_some() @@ -907,5 +925,8 @@ pub fn apply_field( Ok(out) } +mod strong; +pub use strong::{MAX_SWEEPS, Optimized, Phase, Sweep, optimize}; + #[cfg(test)] mod tests; diff --git a/crates/mesh/src/relax/strong.rs b/crates/mesh/src/relax/strong.rs new file mode 100644 index 0000000..0aa2c48 --- /dev/null +++ b/crates/mesh/src/relax/strong.rs @@ -0,0 +1,373 @@ +//! Strong relaxation: an optimisation of the node positions of a whole region toward elements of +//! scaled Jacobian 1 (`mesh.node.relax` with `"mode": "strong"`). +//! +//! The unknowns are the positions of the region's corner nodes: interior nodes in space, nodes on +//! a geometry face in the face's parameters ([`Constraints::slide_face`]); nodes on geometry edges +//! and vertices are held, so feature lines stay as they are. The objective is +//! +//! ```text +//! F = Σₑ φ(qₑ) + ½ ln²(Vₑ/V₀ₑ) for a valid element (q its scaled Jacobian, V its size), +//! φ(q) = (1/q − 1)² phase 1 (smooth): every element toward 1, the worst most, +//! φ(q) = q⁻¹⁶ − 1 phase 2 (max–min): a smooth stand-in for the smallest q, +//! B (2 − q) for an invalid element (B large), +//! ``` +//! +//! growing without bound as an element degenerates (a barrier: no step that lowers `F` inverts an +//! element), while the invalid branch falls as an inverted element untangles and lies far above +//! any valid configuration. The size term is there because the scaled Jacobian does not see an +//! element flattening along one direction: without it the optimiser squeezed elements against +//! the boundary to square up their neighbours (the Pi 5 base's largest edge ratio went from 5 to +//! 34). For the same reason an element longer than the region's largest edge ratio before costs +//! as much as an invalid one. +//! +//! `F` is minimised by block coordinate descent (nonlinear Gauss–Seidel): each sweep visits the +//! movable nodes in label order and lowers the part of `F` that depends on the node, first by +//! trying the centroid of its edge neighbours (Laplacian), then by a gradient step (central +//! differences, backtracking line search, the step length adapted per node). A move is kept when +//! it lowers the local objective, so individual elements may get worse on the way. L. Freitag, +//! P. Plassmann, *Local optimization-based simplicial mesh untangling and improvement*, IJNME 49, +//! 2000; P. Knupp, *Hexahedral and tetrahedral mesh untangling*, Eng. Comput. 17, 2001 (the +//! barrier objective on the Jacobian); J. Nocedal, S. Wright, *Numerical Optimization*, ch. 3 +//! (backtracking). +//! +//! Phase 1 has up to half the sweeps, phase 2 (from phase 1's best state) the rest; a phase stops +//! when a sweep lowers its objective by less than a relative 1e-5, and after the first sweep only +//! nodes whose neighbourhood moved are visited. The sweeps keep the best state they pass through (highest smallest quality, then highest +//! mean). The result is that state if it has no invalid element and its smallest quality is not +//! below the starting one; otherwise nothing moves. Deterministic: no randomness, fixed order. + +use std::collections::{BTreeMap, BTreeSet}; + +use fieldcraft_model::Model; + +use super::{Constraints, Elem, Freedom, P, RegionQuality, State, add, dist, dot, measure, region_quality, scale, sub}; +use crate::{Error, Result}; + +/// Most sweeps of a strong relaxation. +pub const MAX_SWEEPS: usize = 2000; +/// Weight of an invalid element in the objective. +const INVALID: f64 = 1e8; +/// A phase stops when a sweep lowers its objective by less than this fraction. +const CONVERGED: f64 = 1e-5; + +/// One sweep's numbers. +#[derive(Debug, Clone, Copy, PartialEq, Default)] +pub struct Sweep { + pub phase: Phase, + /// The phase's objective. + pub objective: f64, + pub min_quality: f64, + pub mean_quality: f64, + pub invalid: usize, +} + +/// What [`optimize`] did. +#[derive(Debug, Clone, PartialEq, Default)] +pub struct Optimized { + /// The new position of every node that moved. + pub positions: BTreeMap, + /// Corner nodes that moved. + pub moved: usize, + pub max_displacement: f64, + pub before: RegionQuality, + pub after: RegionQuality, + pub mean_before: f64, + pub mean_after: f64, + /// Before the first sweep, then after each. + pub history: Vec, + /// The sweep whose state was kept (0: the start, nothing moved). + pub kept: usize, +} + +/// Weight of the size term. +const SIZE_WEIGHT: f64 = 0.5; + +/// Exponent of the max–min phase. +const MIN_MAX_POWER: i32 = 16; + +/// What a phase minimises. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum Phase { + /// `(1/q − 1)²`: every element toward 1, the worst weighted most. + #[default] + Smooth, + /// `q⁻¹⁶`: a smooth stand-in for the smallest quality, which it lifts. + MinMax, +} + +/// The objective's term for one element of quality `q` whose size went from `v0` to `v`. +fn phi(phase: Phase, q: f64, valid: bool, v: f64, v0: f64) -> f64 { + if valid && q > 0.0 { + let g = if v > 0.0 && v0 > 0.0 { (v / v0).ln() } else { 0.0 }; + let shape = match phase { + Phase::Smooth => (1.0 / q - 1.0).powi(2), + Phase::MinMax => q.powi(-MIN_MAX_POWER) - 1.0, + }; + shape + SIZE_WEIGHT * g * g + } else { + INVALID * (2.0 - q.min(1.0)) + } +} + +/// A linear hexahedron's scaled Jacobian at its corners and centre (the samples +/// [`crate::quality::element_quality`] takes), straight from the trilinear map: the sweeps' +/// inner loop, where the certified check would dominate. A state is certified before it is kept. +fn hex8_scaled_jacobian(x: &[P]) -> (f64, f64) { + let corners = &crate::quality::jacobian::HEX_CORNERS; + let at = |k: usize| x.get(k).copied().unwrap_or([0.0; 3]); + let (mut worst, mut centre) = (f64::INFINITY, 0.0); + for u in corners.iter().map(|c| c.map(|v| v as f64)).chain(std::iter::once([0.5; 3])) { + let mut col = [[0.0; 3]; 3]; + for (k, bits) in corners.iter().enumerate() { + for (d, c) in col.iter_mut().enumerate() { + let w: f64 = (0..3) + .map(|e| { + let b = bits[e] as f64; + if e == d { 2.0 * b - 1.0 } else { b * u[e] + (1.0 - b) * (1.0 - u[e]) } + }) + .product(); + *c = add(*c, scale(at(k), w)); + } + } + let l = dot(col[0], col[0]).sqrt() * dot(col[1], col[1]).sqrt() * dot(col[2], col[2]).sqrt(); + let det = dot(col[0], super::cross(col[1], col[2])); + let q = if l > 0.0 && l.is_finite() { det / l } else { -1.0 }; + worst = worst.min(q); + centre = det; + } + (worst, centre) +} + +/// An element's size for the size term: a linear hexahedron's Jacobian at its centre (its +/// volume, near enough), else the size [`crate::quality::element_quality`] reports. +pub(super) fn size(kind: fieldcraft_model::ElementKind, x: &[P]) -> f64 { + if kind == fieldcraft_model::ElementKind::Hex8 { + hex8_scaled_jacobian(x).1 + } else { + crate::quality::element_quality(kind, x).map_or(0.0, |q| q.size) + } +} + +/// An element's objective term, quality and validity; `certified` takes the certified check +/// for every kind. +fn term_of(phase: Phase, e: &Elem, s: &State<'_>, certified: bool) -> (f64, f64, bool) { + let x = s.coords(e); + let hex = e.kind == fieldcraft_model::ElementKind::Hex8; + let (q, valid, v) = if !certified && hex { + let (q, v) = hex8_scaled_jacobian(&x); + (q, q > 0.0, v) + } else { + let (q, valid) = measure(e, &x); + (q, valid, size(e.kind, &x)) + }; + // An element stretched past the cap is as bad as an invalid one. + let ratio = super::edge_ratio(e.kind, &x); + let stretch = if ratio > s.ratio_cap { INVALID * ratio / s.ratio_cap } else { 0.0 }; + (phi(phase, q, valid, v, e.size0) + stretch, q, valid) +} + +impl State<'_> { + /// The objective over the elements at `n`. + fn local_objective(&self, phase: Phase, n: u32) -> f64 { + let mut f = 0.0; + for &i in self.around.get(&n).map(Vec::as_slice).unwrap_or(&[]) { + if let Some(e) = self.elems.get(i) { + f += term_of(phase, e, self, false).0; + } + } + f + } + + /// `n` moved by `d` under `freedom`, kept when the local objective drops below `f0`; the new + /// local objective when kept. + fn try_lower(&mut self, phase: Phase, n: u32, freedom: Freedom, d: P, f0: f64) -> Option { + let p = self.at(n); + let q = self.moved(freedom, p, d)?; + if dist(p, q) == 0.0 { + return None; + } + let undo = self.place(n, q); + let f = self.local_objective(phase, n); + if f < f0 { + Some(f) + } else { + self.restore(undo); + None + } + } + + /// One optimisation of node `n`; `step` is its adaptive step (relative to the local edge + /// length). Whether it moved. + fn optimize_node(&mut self, phase: Phase, n: u32, step: &mut f64) -> bool { + let freedom = self.freedom(n); + if freedom == Freedom::Fixed { + return false; + } + let mut f0 = self.local_objective(phase, n); + if f0 == 0.0 { + return false; + } + let mut moved = false; + // Laplacian: the centroid of the edge neighbours, kept when it lowers the objective. + let nb: Vec = self.neighbours.get(&n).map(|s| s.iter().copied().collect()).unwrap_or_default(); + if !nb.is_empty() { + let c = scale(nb.iter().fold([0.0; 3], |s, &m| add(s, self.at(m))), 1.0 / nb.len() as f64); + let d = sub(c, self.at(n)); + for k in [1.0, 0.5, 0.25] { + if let Some(f) = self.try_lower(phase, n, freedom, scale(d, k), f0) { + f0 = f; + moved = true; + break; + } + } + } + // Gradient step (central differences in space; on a face the slide maps the difference + // onto the surface, so the gradient is the surface gradient). + let h = self.local_length(n); + if !(h > 0.0) { + return moved; + } + let here = self.at(n); + let eps = 1e-5 * h; + let mut g = [0.0; 3]; + for k in 0..3 { + let mut f = [0.0; 2]; + for (j, sign) in [1.0, -1.0].into_iter().enumerate() { + let mut d = [0.0; 3]; + d[k] = sign * eps; + f[j] = match self.moved(freedom, here, d) { + Some(q) => { + let undo = self.place(n, q); + let v = self.local_objective(phase, n); + self.restore(undo); + v + } + None => f0, + }; + } + g[k] = (f[0] - f[1]) / (2.0 * eps); + } + let gn = dot(g, g).sqrt(); + if !(gn > 0.0 && gn.is_finite()) { + return moved; + } + let dir = scale(g, -1.0 / gn); + let mut t = *step; + for attempt in 0..8 { + if let Some(_f) = self.try_lower(phase, n, freedom, scale(dir, t * h), f0) { + // Grow after a first-try success, keep the shortened step otherwise. + *step = if attempt == 0 { (t * 2.0).min(0.5) } else { t }; + return true; + } + t *= 0.5; + } + *step = (t * 2.0).max(1e-6); + moved + } + + fn sweep_numbers(&self, phase: Phase) -> Sweep { + let mut s = Sweep { phase, min_quality: f64::INFINITY, ..Sweep::default() }; + for e in &self.elems { + let (f, q, valid) = term_of(phase, e, self, true); + s.objective += f; + s.min_quality = s.min_quality.min(q); + s.mean_quality += q; + s.invalid += usize::from(!valid); + } + if self.elems.is_empty() { + s.min_quality = 0.0; + } else { + s.mean_quality /= self.elems.len() as f64; + } + s + } +} + +/// Whether `a` is a better state than `b`: fewer invalid elements, then a higher smallest +/// quality, then a higher mean. +fn better(a: &Sweep, b: &Sweep) -> bool { + a.invalid < b.invalid + || (a.invalid == b.invalid && (a.min_quality > b.min_quality || (a.min_quality == b.min_quality && a.mean_quality > b.mean_quality))) +} + +/// Optimise the corner nodes `nodes` (and the mid-side nodes following them) for up to `sweeps` +/// sweeps, faces' nodes sliding on them and edges' and vertices' nodes held, as `geometry` says. +/// The model is not changed: the result has the new positions. +pub fn optimize(model: &Model, nodes: &BTreeSet, sweeps: usize, geometry: &dyn Constraints) -> Result { + if sweeps == 0 || sweeps > MAX_SWEEPS { + return Err(Error::Invalid(format!("sweeps must be 1 to {MAX_SWEEPS}, got {sweeps}"))); + } + if nodes.len() > super::MAX_NODES { + return Err(Error::Limit(format!("more than {} nodes to relax", super::MAX_NODES))); + } + let mut state = State::new(model, nodes, geometry, true)?; + let original = state.pos.clone(); + let before = region_quality(&state); + // No element may end up more stretched than the most stretched one was. + state.ratio_cap = before.max_edge_ratio.max(1.0); + let start = state.sweep_numbers(Phase::Smooth); + let mut history = vec![start]; + let (mut best, mut best_pos, mut kept) = (start, None::>, 0); + let movers: Vec = nodes.iter().copied().filter(|n| state.around.contains_key(n) && state.freedom(*n) != Freedom::Fixed).collect(); + let mut steps: BTreeMap = movers.iter().map(|&n| (n, 0.1)).collect(); + // Phase 1 has up to half the sweeps (the rest when it stops early); phase 2 the rest. Each + // sweep visits the nodes whose neighbourhood changed in the previous one. + let mut k = 0; + for (phase, budget) in [(Phase::Smooth, sweeps.div_ceil(2)), (Phase::MinMax, sweeps)] { + let mut active: BTreeSet = movers.iter().copied().collect(); + let mut last = state.sweep_numbers(phase).objective; + while k < budget && !active.is_empty() { + k += 1; + let mut next = BTreeSet::new(); + for &n in &movers { + if !active.contains(&n) { + continue; + } + let mut step = steps.get(&n).copied().unwrap_or(0.1); + if state.optimize_node(phase, n, &mut step) { + next.insert(n); + next.extend(state.neighbours.get(&n).into_iter().flatten().copied()); + } + steps.insert(n, step); + } + active = next; + let now = state.sweep_numbers(phase); + history.push(now); + if better(&now, &best) { + best = now; + best_pos = Some(state.pos.clone()); + kept = k; + } + // Converged: the objective no longer drops. + if last - now.objective <= CONVERGED * last.max(f64::MIN_POSITIVE) { + break; + } + last = now.objective; + } + // Phase 2 starts from the best state. + if let Some(p) = &best_pos { + state.pos = p.clone(); + } + } + let accept = best.invalid == 0 && best.min_quality >= start.min_quality; + match best_pos { + Some(p) if accept => state.pos = p, + _ => { + state.pos = original.clone(); + kept = 0; + } + } + let after = region_quality(&state); + let end = state.sweep_numbers(Phase::Smooth); + let mut out = Optimized { before, after, mean_before: start.mean_quality, mean_after: end.mean_quality, history, kept, ..Default::default() }; + for (n, p) in &state.pos { + let Some(o) = original.get(n) else { continue }; + let d = dist(*o, *p); + if d > 0.0 { + out.max_displacement = out.max_displacement.max(d); + out.moved += usize::from(nodes.contains(n)); + out.positions.insert(*n, *p); + } + } + Ok(out) +} diff --git a/crates/mesh/src/relax/tests.rs b/crates/mesh/src/relax/tests.rs index 26accc5..0e767e5 100644 --- a/crates/mesh/src/relax/tests.rs +++ b/crates/mesh/src/relax/tests.rs @@ -145,7 +145,7 @@ fn region_min(m: &Model, labels: &[u32]) -> (f64, usize) { for l in labels { let e = m.element(*l).unwrap(); let x: Vec

= e.nodes.iter().map(|n| m.node(*n).unwrap().xyz).collect(); - let el = Elem { label: *l, kind: e.kind, nodes: e.nodes.clone(), normal: None }; + let el = Elem { label: *l, kind: e.kind, nodes: e.nodes.clone(), normal: None, size0: 0.0 }; let (q, v) = measure(&el, &x); worst = worst.min(q); invalid += usize::from(!v); @@ -342,3 +342,72 @@ fn hostile_parameters_are_refused() { let r = relax(&m, &w, 3, &NoGeometry).unwrap(); assert!(r.positions.is_empty() || r.positions.keys().all(|n| Annulus::faces(m.node(*n).unwrap().xyz).is_empty())); } + +/// Every node of the annulus pushed about (surface nodes on their faces). +fn distorted_annulus(n: usize, amount: f64) -> (Model, Model) { + let original = hex_annulus(n); + let mut m = original.clone(); + perturb(&mut m, [1.0, 1.0, 0.5], 10.0, amount); + (original, m) +} + +/// Strong mode on a distorted block: the region optimised as a whole lifts the smallest and the +/// mean scaled Jacobian well above where they were; nothing inverts; nodes on faces stay on them +/// (the cylinders exactly); nodes on edges and vertices do not move; the same input gives the +/// same output. +#[test] +fn strong_relaxation_of_a_distorted_block_keeps_faces_and_edges() { + let (original, m) = distorted_annulus(6, 0.06); + let geo = Geo { model: original.clone() }; + let nodes: BTreeSet = m.nodes().keys().copied().collect(); + let r = optimize(&m, &nodes, 200, &geo).unwrap(); + assert!(r.after.min_quality > r.before.min_quality + 0.2, "{:?} → {:?}", r.before.min_quality, r.after.min_quality); + assert!(r.mean_after > r.mean_before, "mean {} → {}", r.mean_before, r.mean_after); + assert_eq!(r.after.invalid, 0); + assert!(r.kept > 0 && r.history.len() > 2); + // Within a phase the objective falls sweep by sweep (the descent keeps only lowering moves). + assert!(r.history.windows(2).filter(|w| w[0].phase == w[1].phase).all(|w| w[1].objective <= w[0].objective), "{:?}", r.history); + assert!(r.history.iter().any(|h| h.phase == Phase::MinMax), "{:?}", r.history); + let out = applied(&m, &Relaxed { positions: r.positions.clone(), ..Relaxed::default() }); + let elements: Vec = m.elements().keys().copied().collect(); + let (after, invalid) = region_min(&out, &elements); + assert!((after - r.after.min_quality).abs() < 1e-12 && invalid == 0); + for (n, p) in out.nodes() { + let faces = Annulus::faces(original.node(*n).unwrap().xyz); + assert!(Annulus::on(p.xyz, &faces), "node {n} left its faces {faces:?}"); + if faces.len() >= 2 { + assert_eq!(p.xyz, m.node(*n).unwrap().xyz, "node {n} on an edge moved"); + } + } + assert_eq!(optimize(&m, &nodes, 200, &geo).unwrap(), r, "not deterministic"); +} + +/// Strong mode untangles: an interior node pushed through its neighbours inverts elements; the +/// optimisation brings every element back to valid. +#[test] +fn strong_relaxation_untangles_inverted_elements() { + let (original, mut m) = distorted_annulus(4, 0.0); + // An interior node (i = j = k = 2 of 4) pushed past its neighbour. + let inside = (1 + 2 + 5 * (2 + 5 * 2)) as u32; + let p = m.node(inside).unwrap().xyz; + m.set_node(inside, [p[0] + 0.4, p[1] + 0.4, p[2] + 0.3]).unwrap(); + let geo = Geo { model: original }; + let nodes: BTreeSet = m.nodes().keys().copied().collect(); + let r = optimize(&m, &nodes, 100, &geo).unwrap(); + assert!(r.before.invalid > 0, "{:?}", r.before); + assert_eq!(r.after.invalid, 0, "{:?}", r.after); + assert!(r.after.min_quality > 0.5, "{:?}", r.after); +} + +/// Nothing better found leaves the mesh as it is; bad parameters are errors. +#[test] +fn strong_relaxation_keeps_a_perfect_mesh_and_refuses_bad_sweeps() { + let m = hex_annulus(3); + let geo = Geo { model: m.clone() }; + let nodes: BTreeSet = m.nodes().keys().copied().collect(); + let r = optimize(&m, &nodes, 20, &geo).unwrap(); + assert!(r.after.min_quality >= r.before.min_quality && r.after.invalid == 0); + assert!(optimize(&m, &nodes, 0, &geo).is_err()); + assert!(optimize(&m, &nodes, MAX_SWEEPS + 1, &geo).is_err()); + assert!(optimize(&m, &BTreeSet::new(), 5, &geo).unwrap().positions.is_empty()); +}