diff --git a/crates/engine/src/cmd/plan.rs b/crates/engine/src/cmd/plan.rs index 282918b..efbef43 100644 --- a/crates/engine/src/cmd/plan.rs +++ b/crates/engine/src/cmd/plan.rs @@ -1559,6 +1559,47 @@ mod acceptance { accept(&mut s, &r, 0.6, 6.0); } + /// The enclosure base of the test below without its standoffs: symmetric across x = 28 and + /// y = 21, so its four corners mesh alike: the smallest scaled Jacobian round each corner + /// boss well above 0.5 (the compass search lifts elements up to 0.7) and nearly the same. + #[test] + fn an_enclosure_base_meshes_its_four_corners_alike() { + let mut commands = vec![ + json!({"command": "solid.box", "params": {"length": 56, "width": 42, "height": 16, "body_name": "Base"}}), + json!({"command": "solid.fillet", "params": {"edges": [[0, 0, 8], [56, 0, 8], [56, 42, 8], [0, 42, 8]], "radius": 5}}), + json!({"command": "solid.shell", "params": {"faces": [[28, 21, 16]], "thickness": 3, "body": "Base"}}), + ]; + let corners = [(-1.0, -1.0), (57.0, -1.0), (-1.0, 43.0), (57.0, 43.0)]; + for (x, y) in corners { + commands.push(json!({"command": "solid.cylinder", "params": {"diameter": 13, "height": 16, "base": [x, y, 0], "operation": "join", "targets": ["Base"]}})); + } + for (x, y) in corners { + commands.push(json!({"command": "solid.hole", "params": {"position": [x, y, 16], "diameter": 3.3, "depth": 8}})); + } + let (_d, mut s) = open(json!(commands)); + let plan = s.execute("mesh.plan", &json!({"body": "Base", "size": 1.5, "layers": 2})).unwrap(); + let r = s.execute("mesh.plan.apply", &json!({"plan": plan["name"]})).unwrap(); + accept(&mut s, &r, 0.5, 10.0); + // Every element's scaled Jacobian, by the corner (quadrant) its centroid lies in. + let m = s.execute("mesh.quality.mark", &json!({"limits": {"scaled_jacobian": 1.01, "min_dihedral": null, "max_dihedral": null, "aspect_ratio": null, "min_angle": null, "max_angle": null}})).unwrap(); + let mut worst = [f64::INFINITY; 4]; + for e in m["elements"].as_array().unwrap() { + let at = &e["at"]; + let q = usize::from(at[0].as_f64().unwrap() > 28.0) + 2 * usize::from(at[1].as_f64().unwrap() > 21.0); + for f in e["failed"].as_array().unwrap() { + if f["metric"] == "scaled_jacobian" { + worst[q] = worst[q].min(f["value"].as_f64().unwrap()); + } + } + } + // The layout is mirror symmetric and the smoother visits mirror twins alike; what is left + // is the projection onto the body's tessellation, which the kernel does not make mirror + // symmetric (a chordal tolerance of 0.0165 mm, on 1.5 mm elements). + let (lo, hi) = worst.iter().fold((f64::INFINITY, 0.0f64), |(l, h), w| (l.min(*w), h.max(*w))); + assert!(hi - lo < 0.03, "per corner min scaled Jacobian {worst:?}: {r}"); + assert!(lo > 0.65, "per corner min scaled Jacobian {worst:?}"); + } + /// An enclosure base: a box with rounded vertical edges, shelled open at the top, a screw /// boss joined at each corner (centred outside it, nearly touching the inner fillet) with a /// blind hole, and standoffs on the floor. Each boss's ring frames the whole circle, diff --git a/crates/engine/src/cmd/project/tests.rs b/crates/engine/src/cmd/project/tests.rs index 2a5d009..ad1aeda 100644 --- a/crates/engine/src/cmd/project/tests.rs +++ b/crates/engine/src/cmd/project/tests.rs @@ -197,9 +197,13 @@ fn squeezed_junction_elements_are_reported_not_meshed() { {"op": "connect", "from": [45, "+z"], "to": [47, "-z"], "layers": 2}, {"op": "remove", "blocks": [1, 2, 3, 4, 5]} ]); + // Smoothing lifts its worst elements just over the default floor (0.21): the squeeze is + // reported at any floor an analyst would accept. + p["min_quality"] = json!(0.3); let model = s.model.clone(); let e = s.execute("mesh.hex.project", &p).unwrap_err().to_string(); - assert!(e.contains("smallest scaled Jacobian") && e.contains("worst in block 34"), "{e}"); + // Smoothing leaves the worst element in the squeezed ring block (34) or in block 36. + assert!(e.contains("smallest scaled Jacobian") && (e.contains("worst in block 34") || e.contains("worst in block 36")), "{e}"); assert_eq!(s.model, model); } diff --git a/crates/mesh/src/project/smooth.rs b/crates/mesh/src/project/smooth.rs index ac907fe..76a18db 100644 --- a/crates/mesh/src/project/smooth.rs +++ b/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], 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 { + 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 = (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], method: Smoothing, iterations: usize, tol: f64, sweep_floor: f64) -> SmoothReport { @@ -395,11 +420,21 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec], method: Smoothing, if method == Smoothing::None { return report; } + let order = sweep_order(m.nodes); + let rank: Vec = { + 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| -> Option

{ @@ -461,8 +496,10 @@ fn pass(m: &mut Mesh, g: &dyn Geometry, targets: &[Vec], 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 = - 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 = + 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; }