mesh: strong relax objective never ranks an inverted element below a valid one

In the max-min phase a valid element's term q^-16 reached 1e96 as q -> 0+, far above
an invalid element's 2e8, so the descent kept elements inverted rather than pass
through flat: the distorted plate stayed tangled (4 invalid) once its mesh was the
same on every platform. The shape term now continues as its tangent line below
q = 0.3 (finite, steep) and an invalid element costs 1e11, above any valid one.

Adds a test that the plate's mesh, its perturbation and the strong relaxation are
the same bits with 1 and 4 threads and match a fingerprint taken on macOS (CI checks
it on Linux).

Part of #438.

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

crates/engine/src/cmd/relax/tests.rs

@@ -357,3 +357,37 @@ fn strong_mode_parameters_are_checked() {
     let r = s.execute("mesh.node.relax", &json!({"mode": "strong", "below": 0.01})).unwrap();
     assert_eq!(r["moved"], 0, "{r}");
 }
+
+/// The mesh and the strong relaxation are the same bits on every platform and thread count:
+/// the plate's hex mesh differed between macOS and Linux (the system maths library's `atan2`,
+/// `sin` … are not correctly rounded), so the distorted-plate test passed on one and failed on
+/// the other (#438). The fingerprint is of the nodes after mesh.hex.generate, after the
+/// perturbation and after the strong relaxation; CI checks it on Linux, it was taken on macOS.
+#[test]
+fn mesh_and_strong_relaxation_are_the_same_bits_everywhere() {
+    let fingerprint = |s: &Session, h: &mut u64| {
+        // FNV-1a over every node's label and coordinate bits.
+        for (n, p) in s.model.nodes() {
+            for word in [u64::from(*n), p.xyz[0].to_bits(), p.xyz[1].to_bits(), p.xyz[2].to_bits()] {
+                for b in word.to_le_bytes() {
+                    *h = (*h ^ u64::from(b)).wrapping_mul(0x0000_0100_0000_01b3);
+                }
+            }
+        }
+    };
+    let run = || {
+        let mut h = 0xcbf2_9ce4_8422_2325;
+        let mut s = plate();
+        s.execute("mesh.hex.generate", &json!({"size": 2.0})).unwrap();
+        fingerprint(&s, &mut h);
+        perturb(&mut s, ON_HOLE, 6.0, 0.35);
+        fingerprint(&s, &mut h);
+        s.execute("mesh.node.relax", &json!({"mode": "strong", "below": 0.5, "reach": 4.0})).unwrap();
+        fingerprint(&s, &mut h);
+        h
+    };
+    let one = fieldcraft_parallel::with_threads(1, run);
+    let many = fieldcraft_parallel::with_threads(4, run);
+    assert_eq!(one, many, "1 and 4 threads differ");
+    assert_eq!(one, 0x7445_e04c_0d65_c1d0, "fingerprint {one:#018x}");
+}

crates/mesh/src/relax/strong.rs

@@ -12,9 +12,10 @@
 //! 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
+//! φ continuing as its tangent line below q = 0.3, so it rises steeply but stays finite as an
+//! element degenerates, while the invalid branch falls as an inverted element untangles and lies
+//! above every valid element's term (a barrier: no step that lowers `F` inverts an element unless
+//! it untangles worse ones; and passing through q → 0⁺ never costs more than staying inverted). 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
@@ -43,8 +44,14 @@ use fieldcraft_model::Model;
 use super::{Constraints, Elem, Freedom, P, RegionQuality, State, add, dist, dot, measure, region_quality, scale, sub};
 use crate::{Error, Result};
 
-/// Weight of an invalid element in the objective.
-const INVALID: f64 = 1e8;
+/// Weight of an invalid element in the objective: above every valid element's term (the shape
+/// term is at most [`shape`] at 0, a few 1e9, plus the size term), so no step that untangles an
+/// element ever looks worse than staying inverted.
+const INVALID: f64 = 1e11;
+
+/// Below this quality the shape term continues as its tangent line, so it stays finite (and
+/// below [`INVALID`]) as an element degenerates instead of growing without bound.
+const KNEE: f64 = 0.3;
 
 /// Weight of the size term.
 const SIZE_WEIGHT: f64 = 0.5;
@@ -66,16 +73,30 @@ enum Phase {
 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 { libm::log(v / v0) } 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
+        shape(phase, q) + SIZE_WEIGHT * g * g
     } else {
         INVALID * (2.0 - q.min(1.0))
     }
 }
 
+/// The shape term of a valid element of quality `q > 0`: `φ(q)` from [`KNEE`] up, its tangent
+/// line below. Unbounded, `q⁻¹⁶` reached 1e96 as an element neared flat, so far above an
+/// invalid element's term that the descent kept elements inverted rather than pass through
+/// q → 0⁺ (#438: the distorted plate stayed tangled).
+fn shape(phase: Phase, q: f64) -> f64 {
+    // φ and its derivative at `x`.
+    let at = |x: f64| match phase {
+        Phase::Smooth => ((1.0 / x - 1.0).powi(2), -2.0 * (1.0 / x - 1.0) / (x * x)),
+        Phase::MinMax => (x.powi(-MIN_MAX_POWER) - 1.0, -f64::from(MIN_MAX_POWER) * x.powi(-MIN_MAX_POWER - 1)),
+    };
+    if q >= KNEE {
+        at(q).0
+    } else {
+        let (v, d) = at(KNEE);
+        v + d * (q - KNEE)
+    }
+}
+
 /// 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.