diff --git a/crates/engine/src/cmd/project/tests.rs b/crates/engine/src/cmd/project/tests.rs index 924bd2b..2a5d009 100644 --- a/crates/engine/src/cmd/project/tests.rs +++ b/crates/engine/src/cmd/project/tests.rs @@ -151,7 +151,7 @@ fn cup_params(size: f64, inset: Value) -> Value { fn cup_with_a_handle_at(size: f64) { let (_d, mut s) = cup_session(); let r = s.execute("mesh.hex.project", &cup_params(size, json!(0.5))).unwrap(); - assert!(sj(&r) >= 0.3, "{r}"); + assert_thin_layers_hold(&mut s, &r); let (v, body) = (r["volume"].as_f64().unwrap(), r["body_volume"].as_f64().unwrap()); assert!((v - body).abs() < 0.01 * body, "{r}"); assert!(r["max_gap"].as_f64().unwrap() < 1e-6, "{r}"); @@ -168,14 +168,22 @@ fn cup_with_a_handle() { #[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] fn cup_with_a_handle_refined() { cup_with_a_handle_at(1.0); + cup_with_a_handle_at(0.75); +} + +/// Every element valid, the scaled Jacobian at least 0.5 and no edge ratio above 10 (#418). +fn assert_thin_layers_hold(s: &mut Session, r: &Value) { + assert!(sj(r) >= 0.5, "{r}"); + let q = s.execute("mesh.quality", &json!({})).unwrap(); + assert_eq!(q["valid"].as_u64(), Some(s.model.elements().len() as u64), "{q}"); + assert!(q["metrics"]["scaled_jacobian"]["min"].as_f64().unwrap() >= 0.5, "{q}"); + assert!(q["metrics"]["edge_ratio"]["max"].as_f64().unwrap() <= 10.0, "{q}"); } /// Without the split, the ring block under the upper arm is squeezed between the petal's edge -/// and the window-top plane and folds: the mesher refuses the mesh and names the twisted -/// elements' block and the block junction they touch. (Smoothing now unfolds this grid to a -/// scaled Jacobian of 0.22, so the check runs without it.) +/// and the window-top plane: the mesher refuses the mesh and names the worst element's block. #[test] -fn twisted_junction_elements_are_reported_not_meshed() { +fn squeezed_junction_elements_are_reported_not_meshed() { let (_d, mut s) = cup_session(); let mut p = cup_params(2.0, json!([0.5, 0.3, 0.3])); p["operations"] = json!([ @@ -189,47 +197,39 @@ fn twisted_junction_elements_are_reported_not_meshed() { {"op": "connect", "from": [45, "+z"], "to": [47, "-z"], "layers": 2}, {"op": "remove", "blocks": [1, 2, 3, 4, 5]} ]); - p["smoothing"] = json!("none"); let model = s.model.clone(); let e = s.execute("mesh.hex.project", &p).unwrap_err().to_string(); - assert!(e.contains("twisted") && e.contains("where blocks 5, 8, 12, 29, 30, 33 meet"), "{e}"); + assert!(e.contains("smallest scaled Jacobian") && e.contains("worst in block 34"), "{e}"); assert_eq!(s.model, model); } /// The cup without its handle: the floor's O-grid inside the wall's butterfly, cavity carved. -#[test] -fn cup_floor_double_o_grid() { +/// The wall (4 mm) shares its radial divisions with the floor's longer radial edges; counted +/// from the longest, the wall was cut into slivers that sheared at the cavity's rim, and +/// smoothing collapsed their edges (edge ratio 172 at size 1, #418). +fn cup_floor_double_o_grid_at(size: f64) { let (_d, mut s) = session(json!([ {"command": "solid.cylinder", "params": {"radius": 20, "height": 30, "body_name": "Cup"}}, {"command": "solid.cylinder", "params": {"base": [0, 0, 4], "radius": 16, "height": 27, "operation": "cut", "targets": ["Cup"]}} ])); - let mut p = cup_params(2.0, json!(0.5)); + let mut p = cup_params(size, json!(0.5)); p["operations"] = json!([{"op": "o_grid", "block": 0, "open": [1, 2], "inset": 0.6}, {"op": "remove", "blocks": [1, 2, 3, 4, 5]}]); let r = s.execute("mesh.hex.project", &p).unwrap(); - assert!(sj(&r) >= 0.5, "{r}"); + assert_thin_layers_hold(&mut s, &r); let (v, body) = (r["volume"].as_f64().unwrap(), r["body_volume"].as_f64().unwrap()); assert!((v - body).abs() < 0.02 * body, "{r}"); } -/// Every element valid and no edge ratio above 10: smoothing must not collapse edges (#418). -fn assert_no_collapsed_edges(s: &mut Session) { - let q = s.execute("mesh.quality", &json!({})).unwrap(); - assert_eq!(q["valid"].as_u64(), Some(s.model.elements().len() as u64), "{q}"); - assert!(q["metrics"]["edge_ratio"]["max"].as_f64().unwrap() <= 10.0, "{q}"); +#[test] +fn cup_floor_double_o_grid() { + cup_floor_double_o_grid_at(2.0); + cup_floor_double_o_grid_at(1.0); } -/// The handle-less cup at size 1: smoothing collapsed edges at the cavity's rim to an edge ratio -/// of 172 (#418). #[test] -fn cup_floor_smoothing_keeps_edges() { - let (_d, mut s) = session(json!([ - {"command": "solid.cylinder", "params": {"radius": 20, "height": 30, "body_name": "Cup"}}, - {"command": "solid.cylinder", "params": {"base": [0, 0, 4], "radius": 16, "height": 27, "operation": "cut", "targets": ["Cup"]}} - ])); - let mut p = cup_params(1.0, json!(0.5)); - p["operations"] = json!([{"op": "o_grid", "block": 0, "open": [1, 2], "inset": 0.6}, {"op": "remove", "blocks": [1, 2, 3, 4, 5]}]); - s.execute("mesh.hex.project", &p).unwrap(); - assert_no_collapsed_edges(&mut s); +#[ignore = "heavy: cargo xtask test --heavy (CI heavy job)"] +fn cup_floor_double_o_grid_refined() { + cup_floor_double_o_grid_at(0.75); } /// Thin layers like a sealed enclosure's: a 1.6 mm floor slab and 1 mm walls round a pocket. @@ -248,7 +248,7 @@ fn thin_enclosure_layers() { &json!({"body": "Case", "size": size, "grid": {"planes": {"x": [0, 1, 39, 40], "y": [0, 1, 29, 30], "z": [0, 1.6, 12]}, "carve": [{"i": [1, 1], "j": [1, 1], "k": [1, 1]}]}}), ) .unwrap(); - assert_no_collapsed_edges(&mut s); + assert_thin_layers_hold(&mut s, &r); assert!(sj(&r) > 0.99, "{r}"); } } diff --git a/crates/mesh/src/project/mod.rs b/crates/mesh/src/project/mod.rs index e380276..29a888d 100644 --- a/crates/mesh/src/project/mod.rs +++ b/crates/mesh/src/project/mod.rs @@ -62,11 +62,16 @@ pub const DEFAULT_MAX_ARC_ANGLE: f64 = std::f64::consts::PI / 6.0; /// Samples per side of an outer block face when measuring how much its target turns. const CURVATURE_SAMPLES: usize = 32; +/// The longest an element edge may be, in element sizes, when a class of parallel block edges +/// takes its divisions from the mean of its lengths. +const MAX_STRETCH: f64 = 2.0; + /// How to project. #[derive(Debug, Clone, PartialEq)] pub struct ProjectOptions { /// The element size: each class of parallel block edges gets `round(length / size)` - /// divisions (at least one), from its longest projected edge. + /// divisions (at least one), `length` the geometric mean of its projected edges' lengths, + /// and enough that no element is longer than twice the size. pub size: f64, /// The largest angle (radians) the surface may turn through across one element on a curved /// target face: a class takes enough divisions that its elements are no longer than the @@ -360,9 +365,16 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o } } } - // 2. Divisions per class of parallel edges, from the longest projected edge. + // 2. Divisions per class of parallel edges: the count whose element length is nearest the + // size, in ratio, over all the class's edges (the geometric mean of their projected + // lengths). A class runs through every block its edges are parallel in, so its edges can + // differ a lot in length (a cup's wall, 4 mm thick, and its floor's radial edges, twice + // that); the longest edge alone would cut the thin blocks into flat slivers. No element is + // longer than [`MAX_STRETCH`] sizes, so a coarse size still follows the geometry. let (classes, nclass) = edge_classes(t); - let mut class_len = vec![0.0f64; nclass]; + let mut class_log = vec![(0.0f64, 0usize); nclass]; + let mut class_max = vec![0.0f64; nclass]; + let mut seen: BTreeSet<[u32; 2]> = BTreeSet::new(); for (block, cls) in t.blocks.iter().zip(&classes) { for d in 0..3 { for (k, bits) in CORNERS.iter().enumerate() { @@ -372,6 +384,9 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o let mut e = *bits; e[d] = 1; let (a, b) = (block[k], block[corner_at(e)]); + if !seen.insert(edge_key(a, b)) { + continue; + } let (pa, pb) = (vpos.get(a as usize).copied().unwrap_or_default(), vpos.get(b as usize).copied().unwrap_or_default()); let c = edge_c(edge_key(a, b)); let mut len = 0.0; @@ -381,8 +396,12 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o len += dist(prev, x); prev = x; } - if let Some(l) = class_len.get_mut(cls[d]) { - *l = l.max(len); + if let Some(m) = class_max.get_mut(cls[d]) { + *m = m.max(len); + } + if let Some(l) = class_log.get_mut(cls[d]).filter(|_| len > tol) { + l.0 += (len / opts.size).ln(); + l.1 += 1; } } } @@ -441,8 +460,15 @@ pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, o } } } - let mut class_n: Vec = - class_len.iter().zip(&class_arc).map(|(l, &arc)| ((l / opts.size).round() as usize).max(arc).clamp(1, crate::block::MAX_DIVISIONS)).collect(); + let mut class_n: Vec = class_log + .iter() + .zip(&class_max) + .zip(&class_arc) + .map(|((&(sum, k), &longest), &arc)| { + let mean = if k == 0 { 1.0 } else { (sum / k as f64).exp().round() }; + (mean.max((longest / (MAX_STRETCH * opts.size)).ceil()) as usize).max(arc).clamp(1, crate::block::MAX_DIVISIONS) + }) + .collect(); for &(b, d, n) in &opts.divisions { let cls = classes.get(b).and_then(|c| c.get(d)).ok_or_else(|| Error::Invalid(format!("no direction {d} of block {b}")))?; if !(1..=crate::block::MAX_DIVISIONS).contains(&n) { diff --git a/crates/mesh/src/project/tests.rs b/crates/mesh/src/project/tests.rs index 86f527e..63f3bae 100644 --- a/crates/mesh/src/project/tests.rs +++ b/crates/mesh/src/project/tests.rs @@ -79,7 +79,7 @@ fn seven_block_sphere_converges_without_losing_quality() { assert_eq!(topo.blocks.len(), 7); let exact = 4.0 / 3.0 * PI * r * r * r; let mut last_error = f64::INFINITY; - for size in [4.0, 2.0, 1.0] { + for size in [3.0, 1.5, 0.75] { let m = project_mesh(&topo, &s, &t, &options(size)).unwrap(); let rep = &m.report; assert!(rep.min_scaled_jacobian >= 0.5, "size {size}: {rep:?}");