From 1af87f1f98c02e87db0a86d17d6b28eebd3ceb09 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 23:49:21 -0400 Subject: [PATCH] mesh: a disk's O-grid follows the region's elongation (#198) A smooth loop (straight sides running tangent into arcs, as a slot or a rounded pocket) is meshed as a disk. Its middle block was a square halfway to the loop's quarter points, which on the crank's 6:1 pocket floor squashed elements to a scaled Jacobian of 0.16. The block is now a rectangle on the region's principal axes, inset by the same distance all round (the old square on a circle), its corners paired with the loop nodes nearest their diagonals, its sides spaced as the nodes they pair with, and the O-grid is smoothed (untangled first) before the validity check. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/engine/src/sample.rs | 17 +++ crates/mesh/src/hex/region.rs | 213 ++++++++++++++++++++++++++++++---- 2 files changed, 207 insertions(+), 23 deletions(-) diff --git a/crates/engine/src/sample.rs b/crates/engine/src/sample.rs index 8ee1250..08dce3b 100644 --- a/crates/engine/src/sample.rs +++ b/crates/engine/src/sample.rs @@ -769,6 +769,23 @@ mod tests { } } + /// The whole crank extruded all-hex (#198): its back pocket's floor is a rounded rectangle, + /// straight sides running tangent into corner arcs, so the region is a disk whose O-grid's + /// middle block follows its 6:1 elongation; a square middle block squashed its elements to a + /// scaled Jacobian of 0.16 (now 0.39 at 2 mm, the worst elements elsewhere at 3 and 5 mm). + #[test] + fn the_crank_pocket_floor_keeps_its_elements_square() { + let bytes = fieldcraft_geometry::script_design(CRANK_ARM).unwrap(); + let (_dir, path) = crate::tests::design_file(&bytes); + let mut s = Session::default(); + s.execute("geometry.open", &json!({"path": path})).unwrap(); + for size in [2.0, 3.0, 5.0] { + let r = s.execute("mesh.hex.generate", &json!({"body": "CrankArm", "size": size, "method": "extrude"})).unwrap(); + let min = r["bodies"][0]["quality"]["min_scaled_jacobian"].as_f64().unwrap(); + assert!(min > 0.35, "{size}: {r}"); + } + } + /// `file.sample NAME`, checked clean: ready, no errors or warnings. fn open(name: &str) -> Session { let mut s = Session::default(); diff --git a/crates/mesh/src/hex/region.rs b/crates/mesh/src/hex/region.rs index decea9d..df026b3 100644 --- a/crates/mesh/src/hex/region.rs +++ b/crates/mesh/src/hex/region.rs @@ -19,9 +19,12 @@ //! `m_{k-1}` and `m_{k+1}` divisions, the counts of the lines that end on the sides before //! and after it, so `m_k = (n_{k+1} + n_{k+2} − n_k) / 2`: the loop's count is even and each //! side shorter than the other two together. -//! - **Disk**: one smooth loop (no corners): an O-grid ("butterfly"): a mapped square in the -//! middle, its corners halfway from the centroid to the loop's quarter points, and rings of -//! quadrilaterals from the square out to the loop. The loop's count is a multiple of four. +//! - **Disk**: one smooth loop (no corners; a straight side running tangent into an arc is no +//! corner): an O-grid ("butterfly"): a mapped rectangle in the middle on the region's principal +//! axes, inset from the boundary by the same distance all round (on a circle, the square whose +//! corners are halfway to the boundary; on a slot, a long thin block), and rings of +//! quadrilaterals from it out to the loop, then smoothed. The loop's count is a multiple of +//! four. //! - **Free**: anything else: even loops for the free all-quad mesher //! ([`crate::quad::free_quad_mesh`]). //! @@ -363,7 +366,7 @@ fn annulus(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result let inner_ids: Vec = (0..n).map(|k| (n + ii.get((k + shift) % n).copied().unwrap_or(0)) as u32).collect(); let mut nodes: Vec = outer.iter().chain(inner.iter()).copied().collect(); let mut quads = Vec::with_capacity(quads_count); - if !rings(&mut nodes, &outer_ids, &inner_ids, m, &mut quads)? { + if !rings(&mut nodes, &outer_ids, &inner_ids, m, true, &mut quads)? { return Ok(None); } Ok(Some(QuadMesh { nodes, quads, boundary_nodes: 2 * n, stats: QuadStats::default() })) @@ -371,8 +374,9 @@ fn annulus(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result /// Rings of quadrilaterals between a closed outer node sequence (counter-clockwise) and an inner /// one of the same length, node `k` of each joined by a straight radial line divided into `m`. -/// `false` when a quadrilateral would not be strictly convex. -fn rings(nodes: &mut Vec, outer: &[u32], inner: &[u32], m: usize, quads: &mut Vec<[u32; 4]>) -> Result { +/// `false` when a quadrilateral would not be strictly convex, unless `check` is off (a caller +/// that smooths and checks the mesh itself). +fn rings(nodes: &mut Vec, outer: &[u32], inner: &[u32], m: usize, check: bool, quads: &mut Vec<[u32; 4]>) -> Result { let n = outer.len(); let at = |nodes: &[Q], i: u32| nodes.get(i as usize).copied().unwrap_or([0.0; 2]); let mut ring_ids: Vec> = vec![outer.to_vec()]; @@ -395,7 +399,7 @@ fn rings(nodes: &mut Vec, outer: &[u32], inner: &[u32], m: usize, quads: &mut }; // Outer ring counter-clockwise, inner ring inside it: (p0, p1, q1, q0) turns left. let q = [p0, p1, q1, q0]; - if planar_quad_validity(&q.map(|i| at(nodes, i)))? != Validity::Valid { + if check && planar_quad_validity(&q.map(|i| at(nodes, i)))? != Validity::Valid { return Ok(false); } quads.push(q); @@ -412,12 +416,10 @@ fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result MAX_ANNULUS_NODES { return Ok(None); } - let q = n / 4; let ccw = signed_area(outer) > 0.0; - let order: Vec = if ccw { (0..n as u32).collect() } else { (0..n as u32).rev().collect() }; let at = |nodes: &[Q], i: u32| nodes.get(i as usize).copied().unwrap_or([0.0; 2]); - // The area centroid. - let (mut cx, mut cy, mut a2) = (0.0, 0.0, 0.0); + // The area centroid and second moments about it (Green's theorem over the polygon). + let (mut a2, mut cx, mut cy) = (0.0, 0.0, 0.0); for k in 0..n { let (p, r) = (outer[k], outer[(k + 1) % n]); let c = p[0] * r[1] - r[0] * p[1]; @@ -429,17 +431,90 @@ fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result = outer.clone(); - // The middle square's corners, halfway to the loop's quarter points, and its perimeter. - let corner = |k: usize| { - let b = at(outer, order.get((k * q) % n).copied().unwrap_or(0)); - [centre[0] + 0.5 * (b[0] - centre[0]), centre[1] + 0.5 * (b[1] - centre[1])] + let (mut ixx, mut iyy, mut ixy) = (0.0, 0.0, 0.0); + for k in 0..n { + let (p, r) = (outer[k], outer[(k + 1) % n]); + let (p, r) = ([p[0] - centre[0], p[1] - centre[1]], [r[0] - centre[0], r[1] - centre[1]]); + let c = p[0] * r[1] - r[0] * p[1]; + ixx += c * (p[1] * p[1] + p[1] * r[1] + r[1] * r[1]); + iyy += c * (p[0] * p[0] + p[0] * r[0] + r[0] * r[0]); + ixy += c * (p[0] * r[1] + 2.0 * p[0] * p[1] + 2.0 * r[0] * r[1] + r[0] * p[1]); + } + // ∫y², ∫x² (both × 12) and ∫xy (× 24), signed as the loop: the region's long direction is + // the principal axis of the largest ∫ along², the angle of 0.5 · atan2(2 ∫xy, ∫x² − ∫y²). + let sign = a2.signum(); + let theta = 0.5 * (sign * ixy).atan2(sign * (iyy - ixx)); + let (u, v) = ([theta.cos(), theta.sin()], [-theta.sin(), theta.cos()]); + let local = |p: Q| { + let d = [p[0] - centre[0], p[1] - centre[1]]; + [d[0] * u[0] + d[1] * u[1], d[0] * v[0] + d[1] * v[1]] }; + // Half-widths of the region along its axes, from the centroid to the farther boundary. + let (mut a, mut b) = (0.0f64, 0.0f64); + for &p in outer.iter() { + let l = local(p); + a = a.max(l[0].abs()); + b = b.max(l[1].abs()); + } + if !(a > 0.0 && b > 0.0) { + return Ok(None); + } + // The middle block: a rectangle on the axes inset from the boundary by the same distance all + // round, its half-widths (a − b + c, c) with c = b / (2√2) (on a circle, the square whose + // corners are halfway to the boundary). A square middle block on an elongated region (a slot, + // a rounded pocket) squashes the elements at its long sides (#198). + let half_short = b / (2.0 * std::f64::consts::SQRT_2); + let half_long = (a - b + half_short).max(half_short); + let signs: [Q; 4] = [[-1.0, -1.0], [1.0, -1.0], [1.0, 1.0], [-1.0, 1.0]]; + let world = |l: Q| [centre[0] + l[0] * u[0] + l[1] * v[0], centre[1] + l[0] * u[1] + l[1] * v[1]]; + let corners: [Q; 4] = signs.map(|s| world([s[0] * half_long, s[1] * half_short])); + // Each corner's partner on the loop is a node near the point the inset carries the corner + // out to, along its diagonal; the loop nodes between two partners pair with the side between + // the corners, so the sides' divisions follow: q1 along the long axis, q2 across, their sum + // half the loop's. + let inset = b - half_short; + let ccw_order: Vec = if ccw { (0..n as u32).collect() } else { (0..n as u32).rev().collect() }; + let targets: [Q; 4] = std::array::from_fn(|k| world([signs[k][0] * (half_long + inset), signs[k][1] * (half_short + inset)])); + let dist = |i: usize, t: Q| { + let p = at(outer, ccw_order.get(i % n).copied().unwrap_or(0)); + (p[0] - t[0]).powi(2) + (p[1] - t[1]).powi(2) + }; + let nearest = |t: Q| (0..n).min_by(|&x, &y| dist(x, t).total_cmp(&dist(y, t))).unwrap_or(0); + let p: [usize; 4] = targets.map(nearest); + let run = |k: usize| (p[(k + 1) % 4] + n - p[k]) % n; + let half = n / 2; + // The four partners sit at `start`, `start + q1`, `start + half` and `start + half + q1`: + // the start and q1 (near the mean of the two long runs) that put them nearest their targets. + let mean = (run(0) + run(2)) / 2; + let mut best = (f64::INFINITY, p[0], mean.clamp(1, half - 1)); + for q1 in mean.saturating_sub(2).max(1)..=(mean + 2).min(half - 1) { + for shift in 0..5 { + let start = (p[0] + n + shift - 2) % n; + let offsets = [0, q1, half, half + q1]; + let cost: f64 = (0..4).map(|k| dist(start + offsets[k], targets[k])).sum(); + if cost < best.0 { + best = (cost, start, q1); + } + } + } + let (_, start, q1) = best; + let q2 = half - q1; + let offsets = [0, q1, q1 + q2, 2 * q1 + q2]; + let lengths = [q1, q2, q1, q2]; + let order: Vec = (0..n).map(|k| ccw_order.get((start + k) % n).copied().unwrap_or(0)).collect(); + let mut nodes: Vec = outer.clone(); let mut square: Vec = Vec::with_capacity(n); + // Each side's nodes spaced as the loop nodes they pair with (by chord length along the + // run), so the lines between the two run straight across and do not shear. for k in 0..4 { - let (c0, c1) = (corner(k), corner(k + 1)); - for j in 0..q { - let t = j as f64 / q as f64; + let (c0, c1) = (corners[k], corners[(k + 1) % 4]); + let run: Vec = (0..=lengths[k]).map(|j| at(outer, order.get((offsets[k] + j) % n).copied().unwrap_or(0))).collect(); + let along: Vec = run.windows(2).map(|w| ((w[1][0] - w[0][0]).powi(2) + (w[1][1] - w[0][1]).powi(2)).sqrt()).collect(); + let total: f64 = along.iter().sum(); + let mut acc = 0.0; + for j in 0..lengths[k] { + let t = if total > 0.0 { acc / total } else { j as f64 / lengths[k] as f64 }; + acc += along.get(j).copied().unwrap_or(0.0); nodes.push([c0[0] + t * (c1[0] - c0[0]), c0[1] + t * (c1[1] - c0[1])]); square.push(nodes.len() as u32 - 1); } @@ -451,17 +526,18 @@ fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result() / n as f64; let h = crate::size::checked_size(size, [centre[0], centre[1], 0.0])?; let m = (mean_gap / h).round().max(1.0) as usize; - let count = n.saturating_mul(m).saturating_add(q * q); + let count = n.saturating_mul(m).saturating_add(q1 * q2); if count > opts.max_elements { return Err(Error::Limit(format!("the mesh would have {count} quadrilaterals; the limit is {}", opts.max_elements))); } let mut quads = Vec::with_capacity(count); - if !rings(&mut nodes, &order, &square, m, &mut quads)? { + if !rings(&mut nodes, &order, &square, m, false, &mut quads)? { return Ok(None); } - // The square: a q × q transfinite grid on its perimeter (counter-clockwise from corner 0). + // The middle block: a q1 × q2 transfinite grid on its perimeter (counter-clockwise from + // corner 0). let side = |k: usize, rev: bool| -> Vec { - let mut v: Vec = (0..=q).map(|j| square.get((k * q + j) % n).copied().unwrap_or(0)).collect(); + let mut v: Vec = (0..=lengths[k]).map(|j| square.get((offsets[k] + j) % n).copied().unwrap_or(0)).collect(); if rev { v.reverse(); } @@ -471,6 +547,16 @@ fn disk(loops: &[Vec], size: &dyn SizeField, opts: &QuadOptions) -> Result (Vec, Vec>, Vec) { + let (a, b) = (0.5 * w - r, 0.5 * h - r); + let count = |len: f64| ((len / size).round() as usize).max(1); + let mut curves: Vec> = Vec::new(); + for k in 0..4 { + let (half, other) = if k % 2 == 0 { (a, b) } else { (b, a) }; + // Side k runs along direction d at distance other + r from the centre, then the arc + // round the corner after it. + let t = 0.5 * PI * k as f64; + let (d, nrm) = ([t.cos(), t.sin()], [t.sin(), -t.cos()]); + let at = |s: f64, off: f64| [d[0] * s + nrm[0] * off, d[1] * s + nrm[1] * off]; + let n = count(2.0 * half); + if half > 0.0 { + curves.push((0..=n).map(|j| at(-half + 2.0 * half * j as f64 / n as f64, other + r)).collect()); + } + let centre = at(half, other); + let m = count(0.5 * PI * r); + curves.push( + (0..=m) + .map(|j| { + let u = t - 0.5 * PI + 0.5 * PI * j as f64 / m as f64; + [centre[0] + r * u.cos(), centre[1] + r * u.sin()] + }) + .collect(), + ); + } + let mut pts = Vec::new(); + let mut starts = Vec::new(); + for c in &curves { + starts.push(pts.len()); + pts.extend(c.iter().take(c.len() - 1).copied()); + } + (pts, curves, starts) + } + + /// The smallest corner scaled Jacobian of a planar quad mesh. + fn min_scaled_jacobian(m: &QuadMesh) -> f64 { + m.quads + .iter() + .flat_map(|q| { + let p = q.map(|i| m.nodes[i as usize]); + (0..4).map(move |i| crate::quality::planar_corner(&p, i).1.unwrap_or(-1.0)) + }) + .fold(f64::INFINITY, f64::min) + } + + /// A slot-like region whose straight sides run tangent into its corner arcs (#198): its + /// tangent junctions are no corners, so it is a disk, and its O-grid's middle block follows + /// its elongation, so no element is squashed (a square middle block on a 6:1 region gave a + /// scaled Jacobian of 0.16). + #[test] + fn elongated_disks_keep_their_elements_square() { + for (w, h, r, size) in [ + (82.0, 14.0, 3.0, 1.0), + (82.0, 14.0, 5.0, 3.0), + (82.0, 14.0, 5.0, 2.0), + (82.0, 14.0, 7.0, 1.5), + (40.0, 30.0, 5.0, 2.0), + (10.0, 10.0, 5.0, 0.5), + ] { + let (mut pts, curves, starts) = rounded_rectangle(w, h, r, size); + assert_eq!(classify(std::slice::from_ref(&curves)), Scheme::Disk); + // Pad to a multiple of four nodes, as the interval assignment would. + while pts.len() % 4 != 0 { + let (p, q) = (pts[0], pts[1]); + pts.insert(1, [0.5 * (p[0] + q[0]), 0.5 * (p[1] + q[1])]); + } + let (m, s) = mesh(std::slice::from_ref(&pts), &[starts], Scheme::Disk, &UniformSize(size), &QuadOptions::default()).unwrap(); + assert_eq!(s, Scheme::Disk, "{w} × {h}"); + assert_eq!(&m.nodes[..pts.len()], &pts[..]); + let area: f64 = m.quads.iter().map(|q| signed_area(&q.map(|i| m.nodes[i as usize]))).sum(); + assert!((area - signed_area(&pts)).abs() < 1e-9 * w * h, "{w} × {h}: {area}"); + let min = min_scaled_jacobian(&m); + assert!(min > 0.5, "{w} × {h} r {r}: min scaled Jacobian {min}"); + } + } }