diff --git a/crates/mesh/src/relax.rs b/crates/mesh/src/relax.rs index bef6be7..1cb51dd 100644 --- a/crates/mesh/src/relax.rs +++ b/crates/mesh/src/relax.rs @@ -23,7 +23,8 @@ //! An element's quality is its scaled Jacobian (a triangle's: its smallest angle over 60°), from //! [`crate::quality::element_quality`]; a solid is valid when its certified Jacobian check is, //! a shell when its normal keeps its side. So no move inverts an element, and none makes the -//! worst element around a node worse. +//! worst element around a node worse, and none makes an element's edge ratio exceed the +//! region's largest before the relaxation. //! //! **Constraints** ([`Freedom`]): an interior node moves freely; a node on a geometry face moves //! in the face's parameters and is evaluated back on it ([`Constraints::slide_face`]); a node on a @@ -378,6 +379,12 @@ impl State<'_> { (worst, valid) } + /// The largest edge ratio of the elements around `n`. + fn local_ratio(&self, n: u32) -> f64 { + let around = self.around.get(&n).map(Vec::as_slice).unwrap_or(&[]); + around.iter().filter_map(|&i| self.elems.get(i)).map(|e| edge_ratio(e.kind, &self.coords(e))).fold(0.0, f64::max) + } + /// The worst element around `n`: its index and quality. fn worst_around(&self, n: u32) -> Option<(usize, f64)> { let mut best: Option<(usize, f64)> = None; @@ -428,7 +435,7 @@ impl State<'_> { } let undo = self.place(n, q); let (worst, valid) = self.local(n); - let ok = worst >= before.0 + gain && before.1.iter().zip(&valid).all(|(was, is)| !was || *is); + let ok = worst >= before.0 + gain && before.1.iter().zip(&valid).all(|(was, is)| !was || *is) && self.local_ratio(n) <= self.ratio_cap; if !ok { self.restore(undo); } @@ -606,6 +613,9 @@ pub fn relax(model: &Model, weights: &BTreeMap, iterations: usize, geo let mut state = State::new(model, &movers, geometry, false)?; let original = state.pos.clone(); let before = region_quality(&state); + // A move may not stretch an element past the region's longest one (the scaled Jacobian alone + // lets an element flatten: the crank's worst edge ratio went from 10.5 to 2759). + state.ratio_cap = before.max_edge_ratio.max(1.0); let mut done = 0; for _ in 0..iterations { fieldcraft_model::progress::stage("relax sweep", Some(iterations as f64), "sweeps"); diff --git a/crates/mesh/src/relax/tests.rs b/crates/mesh/src/relax/tests.rs index c450292..4760345 100644 --- a/crates/mesh/src/relax/tests.rs +++ b/crates/mesh/src/relax/tests.rs @@ -433,3 +433,19 @@ fn strong_relaxation_keeps_a_good_mesh_and_refuses_bad_options() { } assert!(morph(&m, &BTreeSet::new(), MorphOptions::default(), &geo).unwrap().positions.is_empty()); } + +/// The guarded mode never stretches an element past the region's longest one before (the +/// scaled Jacobian alone does not see an element flattening: on the crank the worst edge ratio +/// went from 10.5 to 2759, #522). +#[test] +fn guarded_relaxation_keeps_the_largest_edge_ratio() { + for (amount, radius) in [(0.04, 0.5), (0.08, 0.7), (0.12, 1.0)] { + let original = hex_annulus(8); + let mut m = original.clone(); + perturb(&mut m, CENTRE, radius, amount); + let weights = sphere_weights(&m, CENTRE, radius + 0.2, 1.0, Falloff::None, None).unwrap(); + let r = relax(&m, &weights, 50, &Geo { model: original }).unwrap(); + assert!(r.after.max_edge_ratio <= r.before.max_edge_ratio, "{amount}: {} → {}", r.before.max_edge_ratio, r.after.max_edge_ratio); + assert!(r.after.invalid <= r.before.invalid); + } +}