diff --git a/crates/mesh/src/hex/decompose.rs b/crates/mesh/src/hex/decompose.rs new file mode 100644 index 0000000..7cc9b45 --- /dev/null +++ b/crates/mesh/src/hex/decompose.rs @@ -0,0 +1,1331 @@ +//! Automatic decomposition of a swept cross-section into mappable blocks. +//! +//! The regions of a cross-section ([`super::arrangement`]) are meshed by structured schemes +//! ([`super::region`]): mapped four-sided regions, the triangle primitive for three-cornered +//! ones, O-grids for rings and disks. A region none of them fits (an L, a U round a pocket, a +//! plate with several holes) is cut here by straight lines into pieces that fit, before interval +//! assignment sees the cross-section, so the cuts are curves like any other and their counts are +//! assigned with the rest. Nothing of this reaches the model: the cuts exist only inside the +//! mesher. +//! +//! The method follows medial-axis subdivision (Tam and Armstrong, *2D finite element mesh +//! generation by medial axis subdivision*, Advances in Engineering Software 13(5–6):313–324, +//! 1991): a region is split at its concave features, each by a line from the feature into the +//! region, and where it branches, until every piece is a simple primitive (Li, Armstrong, +//! McKeag, *Quad mesh generation for k-sided faces and hex mesh generation for trivalent +//! polyhedra*, Finite Elements in Analysis and Design 26(4):279–301, 1997). Corners are +//! classified by their interior angle as in submapping (Whiteley, White, Benzley, Blacker, *Two +//! and three-quarter dimensional meshing facilitators*, Engineering with Computers 12:144–154, +//! 1996). Each filled region is decomposed on its own, piece by piece: +//! 1. **Holes apart.** A piece with several holes is split by lines across the gap between +//! neighbouring holes (perpendicular to the line between their centres), so that the pieces +//! hold one hole each. +//! 2. **O-grids.** A piece with one hole that is not a ring of two smooth loops (the +//! [`Scheme::Annulus`]) is cut into pieces round the hole: from the outer loop's corners towards +//! the hole's centre or along their bisectors, from the hole's concave corners outwards, or +//! from the hole along four directions (a long plate's layout, whose corners the next step cuts +//! off). The layout whose pieces score best is taken. +//! 3. **Pieces without holes** are split by a search over candidate cuts: lines from concave +//! corners (along the bisector or the extensions of the corner's edges), from the ends and the +//! middle of strongly turning curves (fillets) along their normal and square to the curves +//! beside them, and from convex corners along the bisector, each to the first boundary it +//! meets (or a curve end near there). The cut taken is the best of a two-ply look-ahead: the +//! fewest free pieces (counting their concave features), then the least quality shortfall, +//! then the fewest pieces. Pieces are cut until they are good or no cut helps. +//! +//! **Scoring.** A piece is scored by meshing it alone with its scheme on the target size (coarser +//! for a fine size: the quality of a block depends on its shape): the smallest corner scaled +//! Jacobian of its quadrilaterals. A piece no structured scheme fits counts as free. The score +//! takes the place of the medial axis's geometric rules in Tam and Armstrong: it prefers cuts +//! that cross the region where it narrows, square to its sides, and that leave blocks that do +//! not shear. +//! +//! **Do no harm.** A region's blocks are kept when they are fair (smallest corner scaled Jacobian +//! at least 0.5) or no worse than the free all-quad mesh of the whole region; otherwise the +//! region is left to the free mesher. Cuts never end inside a seeded curve, so edge seeds keep +//! their curves. + +use super::arrangement::{Arrangement, Q}; +use super::interval::Intervals; +use super::region::{self, Scheme}; +use crate::quad::QuadOptions; +use crate::quality::planar_corner; +use crate::size::SizeField; +use crate::{Error, Result}; + +/// The tag of cut lines (no wall behind them). +pub(crate) const CUT: usize = usize::MAX; + +/// Steps (one cut, or one hole step) per region at most. +const MAX_STEPS: usize = 256; +/// Blocks this good (smallest corner scaled Jacobian) are taken whatever the free mesh would do. +const FAIR: f64 = 0.5; +/// How much worse than the region's free mesh (smallest corner scaled Jacobian) its blocks may be. +const HARM_MARGIN: f64 = 0.05; + +/// A piece whose score is at least this is not cut further. +const GOOD: f64 = 0.8; +/// The quality the shortfall is measured from. +const TARGET: f64 = 0.9; +/// The shortfall a cut must save to be worth its extra piece. +const MIN_GAIN: f64 = 2e-3; +/// Candidate cuts tried in one search at most, and over one region (keeps a pathological region +/// bounded). +const SEARCH_BUDGET: usize = 2_000; +const REGION_BUDGET: usize = 20_000; +/// Boundary divisions a piece is scored with at most. +const SCORE_RESOLUTION: f64 = 96.0; +/// Candidate cuts compared per piece at most. +const MAX_CANDIDATES: usize = 40; +/// Boundary points of a piece the search takes at most (larger regions are left as they are). +const MAX_PIECE_POINTS: usize = 4096; +/// The smallest angle (radians) a cut may make with the boundary at its ends (about 25°). +const MIN_CUT_ANGLE: f64 = 0.436; +/// Interior angles beyond this (radians, about 200°) are concave features. +const CONCAVE: f64 = 3.49; +/// Turning (radians) along one curve from which it is a rounded corner (about 40°). +const ROUNDED: f64 = 0.70; + +/// What the cutting needs to know: the target size in the cross-section plane, the merge +/// tolerance and the face mesher options. +pub(crate) struct Settings<'a> { + pub field: &'a dyn SizeField, + pub tol: f64, + pub opts: &'a QuadOptions, +} + +impl Settings<'_> { + fn size(&self, p: Q) -> f64 { + let h = self.field.size([p[0], p[1], 0.0]); + if h > 0.0 && h.is_finite() { h } else { f64::INFINITY } + } +} + +/// One loop of a region as a polygon: its points (counter-clockwise for the outer loop, clockwise +/// for holes, so the region is always on the left), which of them are curve ends, and which +/// edges (`k` → `k + 1`) lie on curves a cut may not end inside (seeded ones). +#[derive(Debug, Clone)] +struct Poly { + pts: Vec, + end: Vec, + locked: Vec, +} + +impl Poly { + fn len(&self) -> usize { + self.pts.len() + } + + fn at(&self, i: usize) -> Q { + self.pts.get(i % self.len().max(1)).copied().unwrap_or([0.0; 2]) + } + + fn is_end(&self, i: usize) -> bool { + self.end.get(i % self.len().max(1)).copied().unwrap_or(false) + } + + /// The interior angle at point `i` (the region on the left of the loop). + fn angle(&self, i: usize) -> f64 { + let n = self.len(); + interior(self.at(i + n - 1), self.at(i), self.at(i + 1)) + } + + /// The curves of the loop, each the points from one curve end to the next (closed), starting + /// at the first end; a loop without ends is one closed curve. With the position where each + /// curve starts. + fn curves(&self) -> (Vec>, Vec) { + let n = self.len(); + let first = (0..n).find(|&i| self.is_end(i)).unwrap_or(0); + let mut curves = Vec::new(); + let mut starts = Vec::new(); + let mut cur = vec![self.at(first)]; + starts.push(first); + for k in 1..=n { + let i = (first + k) % n; + cur.push(self.at(i)); + if self.is_end(i) || k == n { + curves.push(std::mem::take(&mut cur)); + if k < n { + cur.push(self.at(i)); + starts.push(i); + } + } + } + (curves, starts) + } +} + +fn sub(a: Q, b: Q) -> Q { + [a[0] - b[0], a[1] - b[1]] +} + +fn len(a: Q) -> f64 { + a[0].hypot(a[1]) +} + +fn cross(a: Q, b: Q) -> f64 { + a[0] * b[1] - a[1] * b[0] +} + +fn dot(a: Q, b: Q) -> f64 { + a[0] * b[0] + a[1] * b[1] +} + +fn unit(a: Q) -> Option { + let l = len(a); + (l > 0.0 && l.is_finite()).then(|| [a[0] / l, a[1] / l]) +} + +/// The interior angle at `o` between the incoming edge from `a` and the outgoing edge to `b`, the +/// region on the left: π on a straight line, less at a convex corner, more at a concave one. +fn interior(a: Q, o: Q, b: Q) -> f64 { + let (u, v) = (sub(o, a), sub(b, o)); + std::f64::consts::PI - cross(u, v).atan2(dot(u, v)) +} + +/// The inward bisector at a corner of interior angle `theta` whose outgoing edge runs along unit +/// `w`: `w` turned left by half the angle. +fn bisector(w: Q, theta: f64) -> Q { + let half = 0.5 * theta; + [w[0] * half.cos() - w[1] * half.sin(), w[0] * half.sin() + w[1] * half.cos()] +} + +fn centroid(pts: &[Q]) -> Q { + let n = pts.len(); + let (mut cx, mut cy, mut a2) = (0.0, 0.0, 0.0); + for k in 0..n { + let (p, q) = (pts[k], pts[(k + 1) % n]); + let c = cross(p, q); + a2 += c; + cx += (p[0] + q[0]) * c; + cy += (p[1] + q[1]) * c; + } + if a2.abs() > 0.0 { + [cx / (3.0 * a2), cy / (3.0 * a2)] + } else { + let s = pts.iter().fold([0.0, 0.0], |s, p| [s[0] + p[0], s[1] + p[1]]); + [s[0] / n.max(1) as f64, s[1] / n.max(1) as f64] + } +} + +/// The region `r`'s loops as polygons (the region on the left of each). +fn region_polys(arr: &Arrangement, r: usize, locked_curves: &[bool]) -> Option> { + let region = arr.regions.get(r)?; + let mut out = Vec::with_capacity(region.loops.len()); + for lp in ®ion.loops { + let mut poly = Poly { pts: Vec::new(), end: Vec::new(), locked: Vec::new() }; + for &(c, fwd) in lp { + let curve = arr.curves.get(c)?; + let mut ids = curve.points.clone(); + if !fwd { + ids.reverse(); + } + ids.pop(); + let locked = locked_curves.get(c).copied().unwrap_or(false); + for (k, id) in ids.into_iter().enumerate() { + poly.pts.push(*arr.points.get(id)?); + poly.end.push(k == 0); + poly.locked.push(locked); + } + } + if poly.len() < 2 { + return None; + } + out.push(poly); + } + Some(out) +} + +/// A point where a cut ends: a loop point, or a point inside edge `k` (`k` → `k + 1`). +#[derive(Debug, Clone, Copy, PartialEq)] +enum End { + Point(usize, usize), + Edge(usize, usize, Q), +} + +impl End { + fn loop_index(self) -> usize { + match self { + End::Point(l, _) | End::Edge(l, _, _) => l, + } + } + + fn at(self, loops: &[Poly]) -> Q { + match self { + End::Point(l, i) => loops.get(l).map_or([0.0; 2], |p| p.at(i)), + End::Edge(_, _, p) => p, + } + } +} + +/// The first boundary point the ray from `o` along unit `d` meets further than `min_t`, skipping +/// the edges at loop point `from` (where the ray starts), as an end snapped to a loop point +/// within `snap`, with its distance. +fn cast(loops: &[Poly], from: Option<(usize, usize)>, o: Q, d: Q, snap: f64, min_t: f64) -> Option<(f64, End)> { + let mut best: Option<(f64, usize, usize)> = None; + for (l, poly) in loops.iter().enumerate() { + let n = poly.len(); + for k in 0..n { + if let Some((fl, fi)) = from + && fl == l + && (k == fi || (k + 1) % n == fi) + { + continue; + } + let (a, b) = (poly.at(k), poly.at(k + 1)); + let e = sub(b, a); + let den = cross(d, e); + if den.abs() <= 1e-12 * len(e) { + continue; + } + let w = sub(a, o); + let t = cross(w, e) / den; + let u = cross(w, d) / den; + if t > min_t && (-1e-9..=1.0 + 1e-9).contains(&u) && best.is_none_or(|(bt, ..)| t < bt) { + best = Some((t, l, k)); + } + } + } + let (t, l, k) = best?; + let poly = loops.get(l)?; + let x = [o[0] + t * d[0], o[1] + t * d[1]]; + let n = poly.len(); + let (a, b) = (poly.at(k), poly.at(k + 1)); + let near = [(k, a), ((k + 1) % n, b)] + .into_iter() + .filter(|(_, p)| len(sub(*p, x)) <= snap) + .min_by(|x1, x2| len(sub(x1.1, x)).total_cmp(&len(sub(x2.1, x)))); + Some(( + t, + match near { + Some((i, _)) => End::Point(l, i), + None => End::Edge(l, k, x), + }, + )) +} + +/// The nearest curve end of loop `l` to `x` within `radius` (an alternative end for a cut). +fn nearby_end(loops: &[Poly], l: usize, x: Q, radius: f64) -> Option { + let poly = loops.get(l)?; + (0..poly.len()) + .filter(|&i| poly.is_end(i) && len(sub(poly.at(i), x)) <= radius) + .min_by(|&i, &j| len(sub(poly.at(i), x)).total_cmp(&len(sub(poly.at(j), x)))) + .map(|i| End::Point(l, i)) +} + +/// Do the open segments `ab` and `cd` cross or touch? +fn segments_meet(a: Q, b: Q, c: Q, d: Q) -> bool { + let o = |p: Q, q: Q, r: Q| cross(sub(q, p), sub(r, p)); + let (d1, d2, d3, d4) = (o(a, b, c), o(a, b, d), o(c, d, a), o(c, d, b)); + (d1 > 0.0) != (d2 > 0.0) && (d3 > 0.0) != (d4 > 0.0) && d1 != 0.0 && d2 != 0.0 && d3 != 0.0 && d4 != 0.0 +} + +/// Distance from `p` to segment `ab`. +fn seg_dist(p: Q, a: Q, b: Q) -> f64 { + super::arrangement::seg_dist2(p, a, b).0.sqrt() +} + +/// Is `p` inside the region (inside the outer loop, outside the holes)? +fn inside(loops: &[Poly], p: Q) -> bool { + let winding = |poly: &Poly| -> i32 { + let mut wn = 0; + for k in 0..poly.len() { + let (a, b) = (poly.at(k), poly.at(k + 1)); + let side = cross(sub(b, a), sub(p, a)); + if a[1] <= p[1] { + if b[1] > p[1] && side > 0.0 { + wn += 1; + } + } else if b[1] <= p[1] && side < 0.0 { + wn -= 1; + } + } + wn + }; + loops.iter().map(winding).sum::() != 0 +} + +/// Whether the cut from `a` to `b` is admissible: inside the region, clear of other boundary +/// points, square enough to the boundary at both ends, not ending inside a locked edge, and not +/// along the boundary. +fn admissible(loops: &[Poly], a: End, b: End, h: f64) -> bool { + let (pa, pb) = (a.at(loops), b.at(loops)); + let l = len(sub(pb, pa)); + if !(l > 0.25 * h) || a == b { + return false; + } + for e in [a, b] { + if let End::Edge(li, k, _) = e + && loops.get(li).and_then(|p| p.locked.get(k)).copied().unwrap_or(true) + { + return false; + } + } + // The angles the cut makes with the boundary at each end, on both sides. + let Some(d) = unit(sub(pb, pa)) else { return false }; + for (e, dir) in [(a, d), (b, [-d[0], -d[1]])] { + let (prev, next) = match e { + End::Point(li, i) => { + let Some(poly) = loops.get(li) else { return false }; + let n = poly.len(); + (poly.at(i + n - 1), poly.at(i + 1)) + } + End::Edge(li, k, _) => { + let Some(poly) = loops.get(li) else { return false }; + (poly.at(k), poly.at(k + 1)) + } + }; + let o = e.at(loops); + let tip = [o[0] + dir[0], o[1] + dir[1]]; + // The cut must leave into the region: between the outgoing and the incoming edge. + let (a1, a2) = (interior(prev, o, tip), interior(tip, o, next)); + let total = interior(prev, o, next); + if (a1 + a2 - total).abs() > 1e-6 || a1 < MIN_CUT_ANGLE || a2 < MIN_CUT_ANGLE { + return false; + } + } + // No crossing, and clear of every other boundary point. + let clear = 0.05 * h.min(l); + for (li, poly) in loops.iter().enumerate() { + for k in 0..poly.len() { + let (c, e) = (poly.at(k), poly.at(k + 1)); + if segments_meet(pa, pb, c, e) { + return false; + } + let is_end = |x: End| match x { + End::Point(l2, i) => l2 == li && i == k, + End::Edge(..) => false, + }; + if !is_end(a) && !is_end(b) && seg_dist(c, pa, pb) < clear { + return false; + } + } + } + inside(loops, [0.5 * (pa[0] + pb[0]), 0.5 * (pa[1] + pb[1])]) +} + +/// The region split by the cut from `a` to `b` (both on one loop) into two, each a list of loops +/// (holes go to the piece that contains them); `None` when the cut does not split it. +fn split(loops: &[Poly], a: End, b: End) -> Option>> { + // Rebuild the loops with the cut's ends inserted, and find their new positions. + let mut rebuilt: Vec = Vec::with_capacity(loops.len()); + let mut pos: [Option<(usize, usize)>; 2] = [None, None]; + for (l, poly) in loops.iter().enumerate() { + let mut p = Poly { pts: Vec::with_capacity(poly.len() + 2), end: Vec::new(), locked: Vec::new() }; + for k in 0..poly.len() { + let locked = poly.locked.get(k).copied().unwrap_or(false); + p.pts.push(poly.at(k)); + p.end.push(poly.is_end(k)); + p.locked.push(locked); + for (slot, e) in [a, b].into_iter().enumerate() { + if e == End::Point(l, k) + && let (Some(x), Some(last)) = (pos.get_mut(slot), p.end.last_mut()) + { + *last = true; + *x = Some((l, p.pts.len() - 1)); + } + } + // Points inside this edge, in order along it. + let mut inside: Vec<(f64, usize, Q)> = [a, b] + .into_iter() + .enumerate() + .filter_map(|(slot, e)| match e { + End::Edge(le, ke, x) if le == l && ke == k => Some((len(sub(x, poly.at(k))), slot, x)), + _ => None, + }) + .collect(); + inside.sort_by(|x, y| x.0.total_cmp(&y.0)); + for (_, slot, x) in inside { + p.pts.push(x); + p.end.push(true); + p.locked.push(locked); + if let Some(s) = pos.get_mut(slot) { + *s = Some((l, p.pts.len() - 1)); + } + } + } + rebuilt.push(p); + } + let [Some((la, ia)), Some((lb, ib))] = pos else { return None }; + if (la, ia) == (lb, ib) { + return None; + } + let loops = rebuilt; + // A cut between two loops would join them (the hole steps handle holes). + if la != lb { + return None; + } + let poly = loops.get(la)?; + let n = poly.len(); + let piece = |from: usize, to: usize| -> Poly { + let mut p = Poly { pts: Vec::new(), end: Vec::new(), locked: Vec::new() }; + let count = (to + n - from) % n + 1; + for k in 0..count { + let i = (from + k) % n; + p.pts.push(poly.at(i)); + p.end.push(poly.is_end(i)); + // The last edge of the piece is the cut. + p.locked.push(if k + 1 == count { false } else { poly.locked.get(i).copied().unwrap_or(false) }); + } + p + }; + let (pa, pb) = (piece(ia, ib), piece(ib, ia)); + if pa.len() < 3 || pb.len() < 3 { + return None; + } + let holes: Vec = loops.iter().enumerate().filter(|(l, _)| *l != la).map(|(_, p)| p.clone()).collect(); + let mut out = vec![vec![pa], vec![pb]]; + for h in holes { + let probe = h.at(0); + let first = out.first().and_then(|p| p.first()).is_some_and(|p| inside(std::slice::from_ref(p), probe)); + if let Some(piece) = out.get_mut(usize::from(!first)) { + piece.push(h); + } + } + Some(out) +} + +/// How good a piece is: `None` when no structured scheme fits it, else the smallest corner +/// scaled Jacobian of its mesh on the target size. +fn score(piece: &[Poly], s: &Settings) -> Option { + let [poly] = piece else { return ring_score(piece, s) }; + if poly.len() > MAX_PIECE_POINTS { + return None; + } + let (curves, _) = poly.curves(); + let scheme = region::classify(std::slice::from_ref(&curves)); + if scheme == Scheme::Free { + return None; + } + mesh_score(&[curves], scheme, s) +} + +/// A piece with one hole: an annulus when its two loops are smooth, else not structured. +fn ring_score(piece: &[Poly], s: &Settings) -> Option { + if piece.len() != 2 || piece.iter().map(Poly::len).sum::() > MAX_PIECE_POINTS { + return None; + } + let curves: Vec>> = piece.iter().map(|p| p.curves().0).collect(); + let scheme = region::classify(&curves); + if scheme != Scheme::Annulus { + return None; + } + mesh_score(&curves, scheme, s) +} + +/// Mesh loops of curves with a structured scheme on counts from their lengths, and measure. +fn mesh_score(loops: &[Vec>], scheme: Scheme, s: &Settings) -> Option { + let mut goals: Vec> = Vec::with_capacity(loops.len()); + for lp in loops { + let mut g = Vec::with_capacity(lp.len()); + for c in lp { + let l: f64 = c.windows(2).map(|w| len(sub(w[1], w[0]))).sum(); + let mid = c.get(c.len() / 2).copied().unwrap_or([0.0; 2]); + let goal = l / s.size(mid); + if !(goal.is_finite() && goal <= 1e5) { + return None; + } + g.push(goal); + } + goals.push(g); + } + // A fine mesh is scored on a coarser one of the same shape (the quality of a structured + // block depends on its shape, the cost on its size). + let total: f64 = goals.iter().flatten().sum(); + let coarsen = if total > SCORE_RESOLUTION { SCORE_RESOLUTION / total } else { 1.0 }; + let mut ia = Intervals::default(); + let vars: Vec> = goals.iter().map(|g| g.iter().map(|&x| ia.var(x * coarsen, 1)).collect()).collect(); + region::constrain(&mut ia, &vars, scheme); + let counts = ia.solve().ok()?; + let mut nodes: Vec> = Vec::with_capacity(loops.len()); + let mut starts: Vec> = Vec::with_capacity(loops.len()); + for (lp, v) in loops.iter().zip(&vars) { + let mut pts = Vec::new(); + let mut st = Vec::new(); + for (c, &var) in lp.iter().zip(v) { + st.push(pts.len()); + let n = counts.get(var).copied().unwrap_or(1); + let mut d = divide(c, n); + d.pop(); + pts.extend(d); + } + nodes.push(pts); + starts.push(st); + } + let m = region::mesh_structured(&nodes, &starts, scheme, s.field, s.opts).ok()??; + let mut worst = f64::INFINITY; + for q in &m.quads { + let c = q.map(|i| m.nodes.get(i as usize).copied().unwrap_or([0.0; 2])); + for i in 0..4 { + worst = worst.min(planar_corner(&c, i).1.unwrap_or(-1.0)); + } + } + worst.is_finite().then_some(worst) +} + +/// `n + 1` points evenly along a polyline by length. +fn divide(c: &[Q], n: usize) -> Vec { + let n = n.max(1); + let lengths: Vec = c.windows(2).map(|w| len(sub(w[1], w[0]))).collect(); + let total: f64 = lengths.iter().sum(); + let mut out = Vec::with_capacity(n + 1); + let (mut seg, mut acc) = (0usize, 0.0); + for j in 0..=n { + let target = total * j as f64 / n as f64; + while seg + 1 < lengths.len() && acc + lengths.get(seg).copied().unwrap_or(0.0) < target { + acc += lengths.get(seg).copied().unwrap_or(0.0); + seg += 1; + } + let (a, b) = (c.get(seg).copied().unwrap_or([0.0; 2]), c.get(seg + 1).copied().unwrap_or([0.0; 2])); + let l = lengths.get(seg).copied().unwrap_or(0.0); + let t = if l > 0.0 { ((target - acc) / l).clamp(0.0, 1.0) } else { 0.0 }; + out.push([a[0] + t * (b[0] - a[0]), a[1] + t * (b[1] - a[1])]); + } + if let (Some(first), Some(c0)) = (out.first_mut(), c.first()) { + *first = *c0; + } + if let (Some(last), Some(c1)) = (out.last_mut(), c.last()) { + *last = *c1; + } + out +} + +/// The concave corners and concave rounded curves of a piece's loops. +fn concave_features(piece: &[Poly]) -> usize { + let mut count = 0; + for poly in piece { + let n = poly.len(); + count += (0..n).filter(|&i| poly.is_end(i) && poly.angle(i) >= CONCAVE).count(); + let (curves, _) = poly.curves(); + for c in curves { + let turning: f64 = (1..c.len().saturating_sub(1)).map(|k| std::f64::consts::PI - interior(c[k - 1], c[k], c[k + 1])).sum(); + if turning <= -ROUNDED { + count += 1; + } + } + } + count +} + +/// The value of a set of pieces: free pieces, quality shortfall, pieces. +#[derive(Debug, Clone, Copy, PartialEq)] +struct Value { + free: usize, + shortfall: f64, + pieces: usize, +} + +impl Value { + /// A piece's value from its score, or (free) from its concave features: each cut that takes + /// one off counts as progress even while the piece is still free. + fn of(score: Option, piece: &[Poly]) -> Value { + match score { + Some(q) => Value { free: 0, shortfall: (TARGET - q).max(0.0).powi(2), pieces: 1 }, + None => Value { free: 1 + concave_features(piece), shortfall: 0.0, pieces: 1 }, + } + } + + fn add(self, o: Value) -> Value { + Value { free: self.free + o.free, shortfall: self.shortfall + o.shortfall, pieces: self.pieces + o.pieces } + } + + fn better(self, o: Value) -> bool { + if self.free != o.free { + return self.free < o.free; + } + if (self.shortfall - o.shortfall).abs() > MIN_GAIN * (self.pieces.max(o.pieces) - self.pieces.min(o.pieces)).max(1) as f64 { + return self.shortfall < o.shortfall; + } + self.pieces < o.pieces + } +} + +/// The candidate cuts of a piece without holes (see the module docs). +fn candidates(loops: &[Poly], s: &Settings) -> Vec<(End, End)> { + let mut out: Vec<(End, End, f64)> = Vec::new(); + let Some(poly) = loops.first() else { return Vec::new() }; + let n = poly.len(); + // Rays: (loop point, direction). + let mut rays: Vec<(usize, Q)> = Vec::new(); + for i in 0..n { + let (a, o, b) = (poly.at(i + n - 1), poly.at(i), poly.at(i + 1)); + let (Some(u), Some(w)) = (unit(sub(o, a)), unit(sub(b, o))) else { continue }; + let theta = poly.angle(i); + // The inward bisector: the incoming direction turned left by half the interior angle + // from the outgoing one. + let bis = bisector(w, theta); + if poly.is_end(i) { + if theta >= CONCAVE { + rays.push((i, u)); + rays.push((i, [-w[0], -w[1]])); + rays.push((i, bis)); + } else if theta <= 0.83 * std::f64::consts::PI { + rays.push((i, bis)); + } + } + } + // Strongly turning curves: their ends and the point of half their turning, along the normal. + let (curves, starts) = poly.curves(); + for (c, &st) in curves.iter().zip(&starts) { + let m = c.len(); + if m < 3 { + continue; + } + let turns: Vec = (1..m - 1).map(|k| std::f64::consts::PI - interior(c[k - 1], c[k], c[k + 1])).collect(); + let total: f64 = turns.iter().sum(); + if total.abs() < ROUNDED { + continue; + } + let mut acc = 0.0; + let mut half = 1; + for (k, t) in turns.iter().enumerate() { + acc += t; + if acc.abs() >= 0.5 * total.abs() { + half = k + 1; + break; + } + } + for k in [0, half, m - 1] { + let i = (st + k) % n; + let (a, o, b) = (poly.at(i + n - 1), poly.at(i), poly.at(i + 1)); + let (Some(u), Some(w)) = (unit(sub(o, a)), unit(sub(b, o))) else { continue }; + // The inward normal: left of the mean tangent. + if let Some(t) = unit([u[0] + w[0], u[1] + w[1]]) { + rays.push((i, [-t[1], t[0]])); + } + if k == 0 || k == m - 1 { + // At a curve end, also square to the neighbouring curve. + rays.push((i, [-u[1], u[0]])); + rays.push((i, [-w[1], w[0]])); + } + } + } + for (i, d) in rays { + let o = poly.at(i); + let h = s.size(o); + let Some((t, end)) = cast(loops, Some((0, i)), o, d, 0.3 * h, s.tol) else { continue }; + let a = End::Point(0, i); + let x = [o[0] + t * d[0], o[1] + t * d[1]]; + let mut push = |b: End| { + let same = |p: End, q: End| len(sub(p.at(loops), q.at(loops))) <= s.tol; + let known = out.iter().any(|&(p, q, _)| (same(p, a) && same(q, b)) || (same(p, b) && same(q, a))); + if !known && admissible(loops, a, b, h) { + out.push((a, b, t)); + } + }; + push(end); + if let Some(alt) = nearby_end(loops, end.loop_index(), x, (0.75 * h).min(0.3 * t)) { + push(alt); + } + } + out.sort_by(|x, y| x.2.total_cmp(&y.2)); + out.truncate(MAX_CANDIDATES); + out.into_iter().map(|(a, b, _)| (a, b)).collect() +} + +/// The best value reachable from a piece with at most `depth` more cuts, and the first cut. +fn search(piece: &[Poly], depth: usize, s: &Settings, budget: &mut usize) -> (Value, Option<(End, End)>) { + let here = score(piece, s); + let leaf = Value::of(here, piece); + if depth == 0 || here.is_some_and(|q| q >= GOOD) || piece.len() != 1 || *budget == 0 { + return (leaf, None); + } + let mut best = (leaf, None); + for (a, b) in candidates(piece, s) { + if *budget == 0 { + break; + } + *budget -= 1; + let Some(parts) = split(piece, a, b) else { continue }; + let mut total: Option = None; + for part in &parts { + let (v, _) = search(part, depth - 1, s, budget); + total = Some(total.map_or(v, |t| t.add(v))); + } + if let Some(v) = total + && v.better(best.0) + { + best = (v, Some((a, b))); + } + } + best +} + +/// Lines across the gaps between neighbouring holes of a region with several (step 1). +fn separate_holes(loops: &[Poly], s: &Settings) -> Vec> { + let centres: Vec = loops.iter().map(|p| centroid(&p.pts)).collect(); + let mut out: Vec> = Vec::new(); + let mut pairs: Vec<(usize, usize)> = Vec::new(); + for i in 1..loops.len() { + let centre = |k: usize| centres.get(k).copied().unwrap_or([0.0; 2]); + let nearest = (1..loops.len()).filter(|&j| j != i).min_by(|&a, &b| len(sub(centre(a), centre(i))).total_cmp(&len(sub(centre(b), centre(i))))); + let Some(j) = nearest else { continue }; + let key = (i.min(j), i.max(j)); + if pairs.contains(&key) { + continue; + } + pairs.push(key); + let (ci, cj) = (centre(i), centre(j)); + let Some(axis) = unit(sub(cj, ci)) else { continue }; + // The middle of the gap: halfway between where the line of centres leaves each hole. + let leave = |from: Q, d: Q, hole: usize| -> Option { + let poly = loops.get(hole)?; + let mut best: Option = None; + for k in 0..poly.len() { + let (a, b) = (poly.at(k), poly.at(k + 1)); + let e = sub(b, a); + let den = cross(d, e); + if den.abs() <= 1e-12 * len(e) { + continue; + } + let w = sub(a, from); + let (t, u) = (cross(w, e) / den, cross(w, d) / den); + if t > 0.0 && (0.0..=1.0).contains(&u) && best.is_none_or(|bt| t > bt) { + best = Some(t); + } + } + best.map(|t| [from[0] + t * d[0], from[1] + t * d[1]]) + }; + let (Some(pi), Some(pj)) = (leave(ci, axis, i), leave(cj, [-axis[0], -axis[1]], j)) else { continue }; + let m = [0.5 * (pi[0] + pj[0]), 0.5 * (pi[1] + pj[1])]; + let perp = [-axis[1], axis[0]]; + let h = s.size(m); + let hit = |d: Q| cast(loops, None, m, d, 0.3 * h, s.tol); + let (Some((_, e1)), Some((_, e2))) = (hit(perp), hit([-perp[0], -perp[1]])) else { continue }; + // A separator that ends on either of the two holes does not separate them. + if [e1, e2].iter().any(|e| e.loop_index() == i || e.loop_index() == j) { + continue; + } + let (p1, p2) = (e1.at(loops), e2.at(loops)); + if len(sub(p1, p2)) > 0.25 * h { + out.push(vec![p1, p2]); + } + } + out +} + +/// The pieces a region falls into with `cuts` added (a new arrangement of its loops and the +/// cuts), each as its loops; `None` when the arrangement fails. +fn pieces_after(loops: &[Poly], cuts: &[Vec], tol: f64) -> Option>> { + let mut polylines: Vec<(Vec, usize)> = Vec::new(); + let mut vertices: Vec = Vec::new(); + for poly in loops { + let mut pts = poly.pts.clone(); + pts.push(*poly.pts.first()?); + polylines.push((pts, 0)); + vertices.extend((0..poly.len()).filter(|&i| poly.is_end(i)).map(|i| poly.at(i))); + if !(0..poly.len()).any(|i| poly.is_end(i)) { + vertices.push(poly.at(0)); + } + } + for c in cuts { + vertices.extend([*c.first()?, *c.last()?]); + polylines.push((c.clone(), CUT)); + } + let arr = super::arrangement::build(&vertices, &polylines, tol, 0.0).ok()?; + let mut out = Vec::new(); + for (r, region) in arr.regions.iter().enumerate() { + let p = super::layered::interior_point(&arr, region).ok()?; + if inside(loops, p) { + let mut piece = region_polys(&arr, r, &[])?; + // An edge on a locked edge of the region stays locked. + for poly in &mut piece { + for k in 0..poly.len() { + let (a, b) = (poly.at(k), poly.at(k + 1)); + let m = [0.5 * (a[0] + b[0]), 0.5 * (a[1] + b[1])]; + let on_locked = loops + .iter() + .any(|l| (0..l.len()).any(|j| l.locked.get(j).copied().unwrap_or(false) && seg_dist(m, l.at(j), l.at(j + 1)) <= 4.0 * tol)); + if let Some(x) = poly.locked.get_mut(k) { + *x = on_locked; + } + } + } + out.push(piece); + } + } + Some(out) +} + +/// O-grid cuts round the one hole of a region (step 2): of the layouts below, the one whose +/// pieces (each allowed one more cut) score best. +/// - From the outer loop's corners towards the hole's centre. +/// - From the outer loop's corners along their bisectors (towards the centre where a bisector +/// misses the hole). +/// - From the hole's concave corners away from its centre. +/// - From the hole along four directions (the axes, then the diagonals): the outer loop's +/// corners are then cut off by the pieces' own search, the layout of a long plate. +fn o_grid(loops: &[Poly], s: &Settings, budget: &mut usize) -> Vec> { + let (Some(outer), Some(hole)) = (loops.first(), loops.get(1)) else { return Vec::new() }; + let c = centroid(&hole.pts); + let corner = |p: &Poly, i: usize| p.is_end(i) && p.angle(i) <= 0.75 * std::f64::consts::PI; + let concave = |p: &Poly, i: usize| p.is_end(i) && p.angle(i) >= 1.25 * std::f64::consts::PI; + // A cut from loop point (l, i) along d that must end on loop `target`. + let ray = |l: usize, i: usize, d: Q, target: usize| -> Option> { + let o = loops.get(l)?.at(i); + let h = s.size(o); + let (_, end) = cast(loops, Some((l, i)), o, d, 0.3 * h, s.tol)?; + (end.loop_index() == target && admissible(loops, End::Point(l, i), end, h)).then(|| vec![o, end.at(loops)]) + }; + let outer_corners: Vec = (0..outer.len()).filter(|&i| corner(outer, i)).collect(); + let mut layouts: Vec>> = Vec::new(); + layouts.push(outer_corners.iter().filter_map(|&i| ray(0, i, unit(sub(c, outer.at(i)))?, 1)).collect()); + layouts.push( + outer_corners + .iter() + .filter_map(|&i| { + let o = outer.at(i); + let w = unit(sub(outer.at(i + 1), o))?; + let bis = bisector(w, outer.angle(i)); + ray(0, i, bis, 1).or_else(|| ray(0, i, unit(sub(c, o))?, 1)) + }) + .collect(), + ); + layouts.push((0..hole.len()).filter(|&i| concave(hole, i)).filter_map(|i| ray(1, i, unit(sub(hole.at(i), c))?, 0)).collect()); + let s2 = std::f64::consts::FRAC_1_SQRT_2; + for dirs in [[[1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, -1.0]], [[s2, s2], [-s2, s2], [-s2, -s2], [s2, -s2]]] { + let mut cuts = Vec::new(); + for d in dirs { + let Some((_, start)) = cast(std::slice::from_ref(hole), None, c, d, 0.0, 0.0) else { continue }; + let start = match start { + End::Point(_, i) => End::Point(1, i), + End::Edge(_, k, x) => End::Edge(1, k, x), + }; + let o = start.at(loops); + let h = s.size(o); + // From just off the hole, so the ray starts inside the region. + let o2 = [o[0] + d[0] * s.tol * 4.0, o[1] + d[1] * s.tol * 4.0]; + let Some((_, end)) = cast(loops, None, o2, d, 0.3 * h, s.tol) else { continue }; + if end.loop_index() == 0 && admissible(loops, start, end, h) { + cuts.push(vec![o, end.at(loops)]); + } + } + layouts.push(cuts); + } + let mut best: Option<(Value, Vec>)> = None; + for cuts in layouts { + if cuts.len() < 2 { + continue; + } + let Some(pieces) = pieces_after(loops, &cuts, s.tol) else { continue }; + let mut value: Option = None; + for piece in &pieces { + let (v, _) = search(piece, 1, s, budget); + value = Some(value.map_or(v, |t| t.add(v))); + } + if let Some(v) = value + && best.as_ref().is_none_or(|(b, _)| v.better(*b)) + { + best = Some((v, cuts)); + } + } + best.map(|(_, c)| c).unwrap_or_default() +} + +/// The cuts that decompose one region (its loops), and the pieces it ends in. +fn decompose_region(loops: &[Poly], s: &Settings) -> (Vec>, Vec>) { + let mut cuts: Vec> = Vec::new(); + let mut leaves: Vec> = Vec::new(); + let mut work = vec![loops.to_vec()]; + let mut steps = 0; + let mut remaining = REGION_BUDGET; + while let Some(piece) = work.pop() { + steps += 1; + if steps > MAX_STEPS || remaining == 0 || piece.iter().map(Poly::len).sum::() > MAX_PIECE_POINTS { + leaves.push(piece); + continue; + } + let step: Vec> = match piece.len() { + 0 => continue, + 1 => { + let mut budget = SEARCH_BUDGET.min(remaining); + let (_, first) = search(&piece, 2, s, &mut budget); + remaining -= SEARCH_BUDGET.min(remaining) - budget; + if let Some((a, b)) = first + && let Some(parts) = split(&piece, a, b) + { + cuts.push(vec![a.at(&piece), b.at(&piece)]); + work.extend(parts); + continue; + } + Vec::new() + } + 2 if ring_score(&piece, s).is_some() => Vec::new(), + 2 => o_grid(&piece, s, &mut remaining), + _ => separate_holes(&piece, s), + }; + match (step.is_empty(), pieces_after(&piece, &step, s.tol)) { + (false, Some(parts)) if parts.len() > 1 => { + cuts.extend(step); + work.extend(parts); + } + _ => leaves.push(piece), + } + } + (cuts, leaves) +} + +/// The smallest corner scaled Jacobian of the free all-quad mesh of a region on the target +/// size (0 when it cannot be meshed): what the region gets without cuts. +fn free_score(loops: &[Poly], s: &Settings) -> f64 { + let min = if loops.len() > 1 { 6 } else { 4 }; + let mut nodes: Vec> = Vec::with_capacity(loops.len()); + for poly in loops { + let (curves, _) = poly.curves(); + let goals: Vec = + curves.iter().map(|c| c.windows(2).map(|w| len(sub(w[1], w[0]))).sum::() / s.size(c.first().copied().unwrap_or([0.0; 2]))).collect(); + if goals.iter().any(|g| !(g.is_finite() && *g <= 1e5)) { + return 0.0; + } + // Counts near the goals, the loop even and long enough for the free mesher (as the + // section's interval assignment makes them). + let mut counts: Vec = goals.iter().map(|g| (g.round() as usize).max(1)).collect(); + let coarsest = |counts: &[usize]| { + (0..counts.len()).max_by(|&i, &j| { + let r = |k: usize| goals.get(k).copied().unwrap_or(0.0) / counts.get(k).copied().unwrap_or(1) as f64; + r(i).total_cmp(&r(j)).then(j.cmp(&i)) + }) + }; + while counts.iter().sum::() < min || counts.iter().sum::() % 2 == 1 { + let Some(k) = coarsest(&counts) else { return 0.0 }; + if let Some(c) = counts.get_mut(k) { + *c += 1; + } + } + let mut pts = Vec::new(); + for (c, &n) in curves.iter().zip(&counts) { + let mut d = divide(c, n); + d.pop(); + pts.extend(d); + } + nodes.push(pts); + } + let Ok(m) = crate::quad::free_quad_mesh(&nodes, s.field, s.opts) else { return 0.0 }; + m.stats.min_scaled_jacobian.max(0.0) +} + +/// The cuts that decompose the filled regions of `arr` (`used`) into blocks the structured +/// schemes fit, as polylines to add to the arrangement. Each region is decomposed on its own +/// (holes apart, O-grids, then the pieces), and kept when its blocks are fair or not worse than +/// the region's free mesh. `locked[c]`: curve `c` carries an edge seed, so no cut may end inside it. +pub(crate) fn cuts(arr: &Arrangement, used: &[bool], locked: &[bool], s: &Settings) -> Result>> { + if !(s.tol > 0.0 && s.tol.is_finite()) { + return Err(Error::Invalid("the cut tolerance must be a positive number".into())); + } + let mut out = Vec::new(); + for r in 0..arr.regions.len() { + if !used.get(r).copied().unwrap_or(false) { + continue; + } + let Some(loops) = region_polys(arr, r, locked) else { continue }; + // A structured region with a seeded curve is meshed as its seeds set it up: cuts would + // add side equalities the exact seeds may not meet. + let seeded = loops.iter().any(|p| p.locked.iter().any(|&l| l)); + if seeded && score(&loops, s).is_some() { + continue; + } + let (cuts, leaves) = decompose_region(&loops, s); + if cuts.is_empty() { + continue; + } + // Do no harm: blocks of poor quality that are worse than the free mesh of the whole + // region are not taken. + let blocks = leaves.iter().filter_map(|l| score(l, s)).fold(f64::INFINITY, f64::min); + if blocks.is_finite() && blocks < FAIR && blocks < free_score(&loops, s) - HARM_MARGIN { + continue; + } + out.extend(cuts); + } + out.retain(|c| c.len() >= 2 && c.windows(2).all(|w| len(sub(w[1], w[0])) > s.tol) && c.iter().all(|p| p[0].is_finite() && p[1].is_finite())); + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::hex::arrangement; + use crate::hex::section; + use crate::size::UniformSize; + use std::f64::consts::PI; + + /// A profile: loops of curves (polylines), the outer loop counter-clockwise, holes clockwise. + type Profile = Vec>>; + + fn line(a: Q, b: Q) -> Vec { + vec![a, b] + } + + /// An arc about `c` of radius `r` from angle `a0` to `a1` (degrees), in steps of at most 5°. + fn arc(c: Q, r: f64, a0: f64, a1: f64) -> Vec { + let n = ((a1 - a0).abs() / 5.0).ceil().max(1.0) as usize; + (0..=n) + .map(|k| { + let t = (a0 + (a1 - a0) * k as f64 / n as f64).to_radians(); + [c[0] + r * t.cos(), c[1] + r * t.sin()] + }) + .collect() + } + + /// A closed loop through `pts` as straight curves. + fn polygon(pts: &[Q]) -> Vec> { + (0..pts.len()).map(|k| line(pts[k], pts[(k + 1) % pts.len()])).collect() + } + + /// A circle as two half arcs, counter-clockwise or not. + fn circle(c: Q, r: f64, ccw: bool) -> Vec> { + if ccw { vec![arc(c, r, 0.0, 180.0), arc(c, r, 180.0, 360.0)] } else { vec![arc(c, r, 360.0, 180.0), arc(c, r, 180.0, 0.0)] } + } + + /// The cross-section mesh of a profile and its extra (imprinted) curves at size `h`, with its + /// regions decomposed: the mesh, which regions are filled and the arrangement. + struct Meshed { + section: section::Section, + used: Vec, + cuts: usize, + } + + fn mesh_profile(profile: &Profile, extra: &[Vec], h: f64) -> Meshed { + let mut polylines: Vec<(Vec, usize)> = Vec::new(); + let mut vertices: Vec = Vec::new(); + for lp in profile { + for c in lp { + vertices.push(c[0]); + vertices.push(*c.last().unwrap()); + polylines.push((c.clone(), 0)); + } + } + for c in extra { + vertices.push(c[0]); + vertices.push(*c.last().unwrap()); + polylines.push((c.clone(), 1)); + } + let polys: Vec = profile + .iter() + .map(|lp| { + let pts: Vec = lp.iter().flat_map(|c| c[..c.len() - 1].to_vec()).collect(); + Poly { end: vec![false; pts.len()], locked: vec![false; pts.len()], pts } + }) + .collect(); + let tol = 1e-9; + let field = UniformSize(h); + let opts = QuadOptions::default(); + let settings = Settings { field: &field, tol, opts: &opts }; + let filled = |arr: &Arrangement| -> Vec { + arr.regions.iter().map(|r| inside(&polys, super::super::layered::interior_point(arr, r).unwrap())).collect() + }; + let mut arr = arrangement::build(&vertices, &polylines, tol, 0.0).unwrap(); + let mut used = filled(&arr); + let c = cuts(&arr, &used, &vec![false; arr.curves.len()], &settings).unwrap(); + let cuts_added = c.len(); + if !c.is_empty() { + polylines.extend(c.into_iter().map(|c| (c, CUT))); + arr = arrangement::build(&vertices, &polylines, tol, 0.0).unwrap(); + used = filled(&arr); + } + let mut ia = Intervals::default(); + let pts = |c: usize| -> Vec { arr.curves[c].points.iter().map(|&p| arr.points[p]).collect() }; + for c in 0..arr.curves.len() { + let l: f64 = pts(c).windows(2).map(|w| len(sub(w[1], w[0]))).sum(); + ia.var(l / h, 1); + } + let schemes = section::schemes(&arr, &used, &mut ia); + let counts = ia.solve().unwrap(); + let place = |c: usize, n: usize| -> Result> { Ok(divide(&pts(c), n)) }; + let section = section::mesh(&arr, &used, schemes, &counts, &place, &field, &opts).unwrap(); + Meshed { section, used, cuts: cuts_added } + } + + /// Every filled region structured, the quadrilaterals valid with corner scaled Jacobians of at + /// least `min_sj`, conforming (each edge on one or two quadrilaterals, the ones on one along + /// the profile, the others shared by two regions alike), and the area filled. + fn check(m: &Meshed, profile: &Profile, min_sj: f64) -> f64 { + let s = &m.section; + for (r, scheme) in s.schemes.iter().enumerate() { + if m.used[r] { + assert_ne!(*scheme, Scheme::Free, "region {r} is free"); + } + } + let mut worst = f64::INFINITY; + let mut area_sum = 0.0; + let mut edges: std::collections::HashMap<(u32, u32), usize> = std::collections::HashMap::new(); + for q in s.quads.iter().flatten() { + let c = q.map(|v| s.nodes[v as usize]); + for i in 0..4 { + worst = worst.min(planar_corner(&c, i).1.unwrap()); + let (a, b) = (q[i], q[(i + 1) % 4]); + *edges.entry((a.min(b), a.max(b))).or_insert(0) += 1; + } + area_sum += 0.5 * (0..4).map(|k| cross(c[k], c[(k + 1) % 4])).sum::(); + } + assert!(worst >= min_sj, "smallest corner scaled Jacobian {worst} below {min_sj}"); + let boundary: Vec = profile.iter().flatten().flatten().copied().collect(); + let on_profile = |p: Q| profile.iter().flatten().any(|c| c.windows(2).any(|w| seg_dist(p, w[0], w[1]) < 1e-6)); + for (&(a, b), &n) in &edges { + assert!(n <= 2, "edge used {n} times"); + if n == 1 { + let mid = [0.5 * (s.nodes[a as usize][0] + s.nodes[b as usize][0]), 0.5 * (s.nodes[a as usize][1] + s.nodes[b as usize][1])]; + assert!(on_profile(s.nodes[a as usize]) && on_profile(s.nodes[b as usize]), "an open edge inside the profile near {mid:?}"); + } + } + let _ = boundary; + let profile_area: f64 = profile + .iter() + .map(|lp| { + let pts: Vec = lp.iter().flat_map(|c| c[..c.len() - 1].to_vec()).collect(); + 0.5 * (0..pts.len()).map(|k| cross(pts[k], pts[(k + 1) % pts.len()])).sum::() + }) + .sum(); + // Curved boundaries are chords of the profile's polylines between the mesh nodes. + assert!((area_sum - profile_area).abs() < 2e-3 * profile_area.abs(), "area {area_sum} against {profile_area}"); + worst + } + + #[test] + fn an_l_shape_is_cut_at_its_concave_corner() { + let profile = vec![polygon(&[[0.0, 0.0], [40.0, 0.0], [40.0, 15.0], [15.0, 15.0], [15.0, 40.0], [0.0, 40.0]])]; + let m = mesh_profile(&profile, &[], 2.5); + assert_eq!(m.cuts, 1); + assert!(check(&m, &profile, 0.99) > 0.99); + } + + #[test] + fn a_plate_with_a_hole_gets_an_o_grid() { + let mut profile = vec![polygon(&[[0.0, 0.0], [60.0, 0.0], [60.0, 40.0], [0.0, 40.0]])]; + profile.push(circle([30.0, 20.0], 8.0, false)); + let m = mesh_profile(&profile, &[], 2.5); + assert_eq!(m.cuts, 4); + // The O-grid's lines split the plate's corners (atan(20 / 30) = 34° of 90°). + check(&m, &profile, 0.5); + } + + #[test] + fn a_plate_with_a_slot_gets_an_o_grid() { + let mut profile = vec![polygon(&[[0.0, 0.0], [80.0, 0.0], [80.0, 40.0], [0.0, 40.0]])]; + // An obround slot: two R6 ends at x = 30 and 50 (clockwise). + profile.push(vec![ + arc([30.0, 20.0], 6.0, 270.0, 90.0), + line([30.0, 26.0], [50.0, 26.0]), + arc([50.0, 20.0], 6.0, 90.0, -90.0), + line([50.0, 14.0], [30.0, 14.0]), + ]); + let m = mesh_profile(&profile, &[], 2.5); + check(&m, &profile, 0.5); + } + + #[test] + fn a_plate_with_several_holes_is_split_between_them() { + let mut profile = vec![polygon(&[[0.0, 0.0], [90.0, 0.0], [90.0, 40.0], [0.0, 40.0]])]; + for x in [20.0, 45.0, 70.0] { + profile.push(circle([x, 20.0], 6.0, false)); + } + let m = mesh_profile(&profile, &[], 2.0); + check(&m, &profile, 0.5); + } + + /// The half crank's cross-section (y ≤ 0) as the extrusion projects it: the outline, the + /// bosses' circles where they rise above the arm, and the pocket's outline. + fn crank_half() -> (Profile, Vec>) { + let deg = |p: Q, c: Q| (p[1] - c[1]).atan2(p[0] - c[0]).to_degrees(); + let spindle = [0.0, 0.0]; + let pedal = [170.0, 0.0]; + let j1 = [13.455063072993756, -18.65372021077189]; + let r1 = [15.471735240951828, -18.362928548901984]; + let r2 = [160.67227699330076, -12.384082712040556]; + let j2 = [161.58293237729853, -12.415835559271503]; + let c1 = [15.883151316213851, -28.354461805265395]; + let c2 = [160.91912663845798, -18.3790026658586]; + let a = |c: Q, r: f64, p: Q, q: Q, ccw: bool| { + let (mut a0, mut a1) = (deg(p, c), deg(q, c)); + if ccw && a1 < a0 { + a1 += 360.0; + } + if !ccw && a1 > a0 { + a0 += 360.0; + } + let mut v = arc(c, r, a0, a1); + v[0] = p; + *v.last_mut().unwrap() = q; + v + }; + let outer = vec![ + a(spindle, 23.0, [-23.0, 0.0], j1, true), + a(c1, (j1[0] - c1[0]).hypot(j1[1] - c1[1]), j1, r1, false), + line(r1, r2), + a(c2, 6.0, r2, j2, false), + a(pedal, 15.0, j2, [185.0, 0.0], true), + line([185.0, 0.0], [177.0, 0.0]), + a(pedal, 7.0, [177.0, 0.0], [163.0, 0.0], false), + line([163.0, 0.0], [6.35, 0.0]), + line([6.35, 0.0], [6.35, -4.85]), + a([4.85, -4.85], 1.5, [6.35, -4.85], [4.85, -6.35], false), + line([4.85, -6.35], [-4.85, -6.35]), + a([-4.85, -4.85], 1.5, [-4.85, -6.35], [-6.35, -4.85], false), + line([-6.35, -4.85], [-6.35, 0.0]), + line([-6.35, 0.0], [-23.0, 0.0]), + ]; + let extra = vec![ + a(spindle, 23.0, j1, [23.0, 0.0], true), + a(pedal, 15.0, [155.0, 0.0], j2, true), + line([48.0, 0.0], [48.0, -2.0]), + a([53.0, -2.0], 5.0, [48.0, -2.0], [53.0, -7.0], true), + line([53.0, -7.0], [125.0, -7.0]), + a([125.0, -2.0], 5.0, [125.0, -7.0], [130.0, -2.0], true), + line([130.0, -2.0], [130.0, 0.0]), + ]; + (vec![outer], extra) + } + + #[test] + fn the_crank_half_section_is_all_blocks() { + let (profile, extra) = crank_half(); + let m = mesh_profile(&profile, &extra, 3.0); + // The arm meets the bosses about 20° off tangent: the corner there bounds the quality. + let worst = check(&m, &profile, 0.37); + assert!(m.cuts >= 4, "{}", m.cuts); + let _ = worst; + } + + /// One cut joining the T's two concave corners leaves two rectangles. + #[test] + fn a_t_shape_is_cut_between_its_concave_corners() { + let profile = vec![polygon(&[[0.0, 0.0], [60.0, 0.0], [60.0, 10.0], [35.0, 10.0], [35.0, 40.0], [25.0, 40.0], [25.0, 10.0], [0.0, 10.0]])]; + let m = mesh_profile(&profile, &[], 2.0); + assert_eq!(m.cuts, 1); + check(&m, &profile, 0.99); + } + + /// A cut never ends inside a curve that carries a seed: the L's long bottom side is locked, + /// so its corner is cut the other way (along the bottom side's direction, ending on the + /// unlocked left side). + #[test] + fn cuts_do_not_end_inside_seeded_curves() { + let pts = [[0.0, 0.0], [40.0, 0.0], [40.0, 15.0], [15.0, 15.0], [15.0, 40.0], [0.0, 40.0]]; + let mut polylines: Vec<(Vec, usize)> = Vec::new(); + for k in 0..pts.len() { + polylines.push((vec![pts[k], pts[(k + 1) % pts.len()]], 0)); + } + let arr = arrangement::build(&pts, &polylines, 1e-9, 0.0).unwrap(); + let used: Vec = arr.regions.iter().map(|r| r.loops.len() == 1).collect(); + let field = UniformSize(2.5); + let opts = QuadOptions::default(); + let settings = Settings { field: &field, tol: 1e-9, opts: &opts }; + let lies_on = |c: &arrangement::Curve, a: Q, b: Q| c.points.iter().all(|&p| seg_dist(arr.points[p], a, b) < 1e-9); + for (a, b) in [([0.0, 0.0], [40.0, 0.0]), ([0.0, 40.0], [0.0, 0.0])] { + let locked: Vec = arr.curves.iter().map(|c| lies_on(c, a, b)).collect(); + let c = cuts(&arr, &used, &locked, &settings).unwrap(); + assert_eq!(c.len(), 1); + for p in [c[0][0], *c[0].last().unwrap()] { + assert!(!(seg_dist(p, a, b) < 1e-9 && len(sub(p, a)) > 1e-9 && len(sub(p, b)) > 1e-9), "{c:?} ends inside the seeded side"); + } + } + } + + #[test] + fn hostile_settings_are_an_error() { + let arr = arrangement::build(&[], &[(vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 0.0]], 0)], 1e-9, 0.0).unwrap(); + let field = UniformSize(1.0); + let opts = QuadOptions::default(); + assert!(cuts(&arr, &[true], &[], &Settings { field: &field, tol: f64::NAN, opts: &opts }).is_err()); + let _ = PI; + } +} diff --git a/crates/mesh/src/hex/layered.rs b/crates/mesh/src/hex/layered.rs index fa0314a..2f218ec 100644 --- a/crates/mesh/src/hex/layered.rs +++ b/crates/mesh/src/hex/layered.rs @@ -17,6 +17,10 @@ //! dimensional geometries*, Proc. 7th International Meshing Roundtable, 1998). Which intervals a //! region fills is decided by the inside test at its centre ([`super::geom::Inside`]). //! +//! **Blocks.** Regions no structured scheme fits are cut into blocks that fit +//! ([`super::decompose`]: holes apart, O-grids, cuts from concave corners and fillets), and the +//! arrangement is laid out again with the cuts, so they are curves like the others. +//! //! **Meshing.** Interval assignment ([`super::interval`]) gives every curve its count so that //! every region's scheme fits ([`super::section`]: four-sided regions mapped with opposite sides //! equal, three-sided ones the triangle primitive, rings and disks O-grids, the rest even loops @@ -30,6 +34,7 @@ use std::f64::consts::{FRAC_PI_4, PI, TAU}; use super::arrangement::{self, Arrangement, Q}; use super::brep::{Brep, Face}; +use super::decompose; use super::geom::{Inside, Surface, add, dist, scale, unit}; use super::interval::Intervals; use super::region::Scheme; @@ -351,7 +356,7 @@ impl SizeField for PlaneField<'_> { } } -/// The cross-section of a frame the body sweeps in: the arrangement (cut at reflex corners), +/// The cross-section of a frame the body sweeps in: the arrangement (cut into blocks), /// the level intervals and which of them each region fills. It does not depend on the element /// size, so refinement retries reuse it. pub(crate) struct Prepared { @@ -403,17 +408,20 @@ pub(crate) fn prepare(ctx: &Ctx, setup: Setup) -> Result { let ends: Vec = if frame.periodic() { rest.iter().copied().chain(std::iter::once(first_level + TAU)).collect() } else { rest.to_vec() }; let intervals: Vec<(f64, f64)> = levels.iter().copied().zip(ends).collect(); let mut inside = Inside::new(brep.surface); - let first = arrangement::build(&vertices, &polylines, tol, coincide)?; - let (_, used_first) = occupancy(&first, frame, &intervals, &mut inside)?; - // Cut the filled regions at their reflex corners, then lay the arrangement out again. - let cut_lines = cuts(&first, &used_first, tol); - let arr = if cut_lines.is_empty() { - first - } else { - polylines.extend(cut_lines); - arrangement::build(&vertices, &polylines, tol, coincide)? - }; - let (filled, used) = occupancy(&arr, frame, &intervals, &mut inside)?; + let mut arr = arrangement::build(&vertices, &polylines, tol, coincide)?; + let (mut filled, mut used) = occupancy(&arr, frame, &intervals, &mut inside)?; + // Cut the filled regions into blocks the structured schemes fit, and lay the arrangement out + // again with the cuts. + let h_ref = levels.first().copied().unwrap_or(0.0); + let plane = PlaneField { field: ctx.field, frame: *frame, h: h_ref }; + let settings = decompose::Settings { field: &plane, tol, opts: &ctx.opts.quad }; + let locked = seeded_curves(ctx, frame, h_ref, &intervals, setup.h_tol, &arr); + let cut_lines = decompose::cuts(&arr, &used, &locked, &settings)?; + if !cut_lines.is_empty() { + polylines.extend(cut_lines.into_iter().map(|c| (c, decompose::CUT))); + arr = arrangement::build(&vertices, &polylines, tol, coincide)?; + (filled, used) = occupancy(&arr, frame, &intervals, &mut inside)?; + } Ok(Prepared { setup, arr, intervals, filled, used }) } @@ -717,95 +725,19 @@ fn occupancy(arr: &Arrangement, frame: &Frame, intervals: &[(f64, f64)], inside: Ok((filled, used)) } -/// The tag of cut lines (no wall behind them). -const CUT: usize = usize::MAX; - -/// Straight cuts from the reflex corners of the filled regions into them, each along one of the -/// corner's edges (the shorter way) to the first boundary it meets. A rectilinear region falls -/// into rectangles this way; others lose their reflex corners. This is the explicit form of the -/// virtual decomposition behind submapping (Whiteley, White, Benzley, Blacker 1996), done in the -/// arrangement so that interval assignment sees every piece. -fn cuts(arr: &Arrangement, used: &[bool], tol: f64) -> Vec<(Vec, usize)> { - let at = |p: usize| arr.points.get(p).copied().unwrap_or([0.0; 2]); - let segments: Vec<(usize, usize)> = arr.curves.iter().flat_map(|c| c.points.windows(2).map(|w| (w[0], w[1]))).collect(); - // The nearest boundary hit of the ray from point v along unit d (ignoring segments at v). - let hit = |v: usize, d: Q| -> Option<(f64, Q, Option)> { - let o = at(v); - let mut best: Option<(f64, Q, Option)> = None; - for &(a, b) in &segments { - if a == v || b == v { - continue; - } - let (pa, pb) = (at(a), at(b)); - let e = [pb[0] - pa[0], pb[1] - pa[1]]; - let den = d[0] * e[1] - d[1] * e[0]; - if den.abs() <= f64::EPSILON * (e[0].abs() + e[1].abs()) { - continue; - } - let w = [pa[0] - o[0], pa[1] - o[1]]; - let t = (w[0] * e[1] - w[1] * e[0]) / den; - let u = (w[0] * d[1] - w[1] * d[0]) / den; - if t > tol && (-1e-12..=1.0 + 1e-12).contains(&u) && best.is_none_or(|(bt, _, _)| t < bt) { - let x = [o[0] + t * d[0], o[1] + t * d[1]]; - // End at a segment end that is this close to the hit. - let near = [(a, pa), (b, pb)].into_iter().find(|(_, p)| ((p[0] - x[0]).powi(2) + (p[1] - x[1]).powi(2)).sqrt() <= 0.1 * t); - best = Some((t, near.map_or(x, |n| n.1), near.map(|n| n.0))); - } - } - best - }; - let mut out: Vec<(Vec, usize)> = Vec::new(); - let mut seen: Vec<(usize, usize)> = Vec::new(); - for (r, region) in arr.regions.iter().enumerate() { - if !used.get(r).copied().unwrap_or(false) { - continue; - } - for lp in ®ion.loops { - let n = lp.len(); - for k in 0..n { - let (Some(&(pc, pf)), Some(&(nc, nf))) = (lp.get((k + n - 1) % n), lp.get(k)) else { continue }; - let pts = |c: usize, f: bool| { - let mut v = arr.curves.get(c).map(|x| x.points.clone()).unwrap_or_default(); - if !f { - v.reverse(); - } - v - }; - let (prev, next) = (pts(pc, pf), pts(nc, nf)); - let (Some(&a), Some(&v), Some(&b)) = (prev.get(prev.len().wrapping_sub(2)), next.first(), next.get(1)) else { continue }; - let (pa, pv, pb) = (at(a), at(v), at(b)); - let (u, w) = ([pv[0] - pa[0], pv[1] - pa[1]], [pb[0] - pv[0], pb[1] - pv[1]]); - let turn = (u[0] * w[1] - u[1] * w[0]).atan2(u[0] * w[0] + u[1] * w[1]); - // The region is on the left: a right turn of more than 45° is a reflex corner. - if PI - turn < 1.25 * PI { - continue; - } - let norm2 = |x: Q| (x[0] * x[0] + x[1] * x[1]).sqrt(); - let (lu, lw) = (norm2(u), norm2(w)); - if !(lu > 0.0 && lw > 0.0) { - continue; - } - let rays = [[u[0] / lu, u[1] / lu], [-w[0] / lw, -w[1] / lw]]; - let best = rays.iter().filter_map(|&d| hit(v, d)).min_by(|x, y| x.0.total_cmp(&y.0)); - if let Some((_, end, end_point)) = best { - if let Some(e) = end_point { - let key = (v.min(e), v.max(e)); - if seen.contains(&key) { - continue; - } - seen.push(key); - } - out.push((vec![pv, end], CUT)); - } - } - } +/// Which curves of `arr` carry an edge seed (a cut may not end inside them). +fn seeded_curves(ctx: &Ctx, frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, arr: &Arrangement) -> Vec { + let curves: Vec> = + arr.curves.iter().map(|cv| cv.points.iter().map(|&p| frame.point(arr.points.get(p).copied().unwrap_or([0.0; 2]), h_ref)).collect()).collect(); + match Seeds::of(ctx, frame, h_ref, intervals, h_tol, &curves) { + Ok(seeds) => seeds.curves.iter().map(Option::is_some).collect(), + Err(_) => vec![false; arr.curves.len()], } - out } /// A point well inside a region: off the middle of one of its outer edges, inside the outer /// loop, outside the holes and clear of every boundary segment. -fn interior_point(arr: &Arrangement, region: &arrangement::Region) -> Result { +pub(crate) fn interior_point(arr: &Arrangement, region: &arrangement::Region) -> Result { let at = |p: usize| arr.points.get(p).copied().unwrap_or([0.0; 2]); let polygon = |lp: &[(usize, bool)]| -> Vec { let mut out = Vec::new(); diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index 8ca17d7..4babfb6 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -37,6 +37,7 @@ mod arrangement; mod brep; mod cells; +mod decompose; pub(crate) mod geom; mod interval; mod layered; diff --git a/docs/architecture.md b/docs/architecture.md index 305f9ea..3b80fb3 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -95,8 +95,10 @@ Rules beyond the layer order: (`fieldcraft-mesh::size`) sets interior sizes and edge divisions, with even counts per loop. Mapped and swept hex build on these face meshes (`fieldcraft-mesh::hex`, `mesh.hex.generate`): the B-rep topology is recovered from the face-tagged surface; a body is extruded (straight - multi-level sweep: the levels' outlines imprinted into one cross-section, cut at reflex corners, - regions mapped, O-gridded or free-meshed, then stacked), revolved (the same in a rotating + multi-level sweep: the levels' outlines imprinted into one cross-section, cut automatically + into blocks (holes apart, O-grids, cuts from concave corners and fillets chosen by the quality + of the blocks they leave), regions mapped, triangle primitives, O-gridded or free-meshed, then + stacked), revolved (the same in a rotating frame), swept one-to-one through rings of linking faces (least-squares affine layers, Roca and Sarrate 2010), or mapped as a logical cube (Gordon and Thiel 1982); interval assignment couples the counts of shared curves. Nodes on analytic faces (`Shape::face_surfaces`) are diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index 8e605fd..634ee48 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -1001,8 +1001,8 @@ record or the publisher; the others are cited by venue only. - **Citation.** M. Whiteley, D. White, S. Benzley, T. Blacker. *Two and three-quarter dimensional meshing facilitators.* Engineering with Computers 12:144–154, 1996. - **Use.** Corners by interior angle (≤ 135° end, ≥ 225° reflex); cross-section regions are cut from - their reflex corners to the boundary, which decomposes rectilinear sections into rectangles (an explicit - form of submapping's virtual decomposition), then mapped. + their reflex corners to the boundary (B8.6b), which decomposes rectilinear sections into rectangles (an + explicit form of submapping's virtual decomposition), then mapped. **B8.6a Midpoint subdivision: Tam and Armstrong 1991; Li, Armstrong, McKeag 1997.** **Verdict: BUILD ON (triangle primitive).** @@ -1017,6 +1017,24 @@ record or the publisher; the others are cited by venue only. shorter than the other two together). The centre is the mean of the split nodes; the interior is then smoothed (`quad::smooth`) and every quadrilateral checked; a fold falls back to the free mesher. +**B8.6b Automatic cross-section decomposition (medial-axis subdivision).** **Verdict: BUILD ON (blocks +of 2.5D bodies).** +- **Citations.** Tam and Armstrong 1991 (B8.6a): a region is subdivided at its concave features and the + branches of its medial axis into simple primitives. Li, Armstrong, McKeag 1997 (B8.6a): the primitives. + Whiteley et al. 1996 (B8.6): vertex classification. +- **Use.** `hex::decompose` cuts the regions of a multi-level sweep's cross-section that no structured + scheme fits into blocks, before interval assignment (B8.5), so the cuts are curves of the arrangement + and their counts are assigned with the rest; nothing reaches the model. Per region: several holes are + separated by lines across the gaps between neighbouring holes; a piece with one hole gets an O-grid, + the layout (lines from the outer corners towards the hole or along their bisectors, from the hole's + corners, or along four directions from the hole) chosen by the score of its pieces; pieces without + holes are cut by a two-ply search over lines from concave corners and fillets (bisector, edge + extensions, normals of strongly turning curves) and convex corners' bisectors. In place of the + medial axis's geometric rules each candidate is scored by meshing its pieces with their schemes (the + smallest corner scaled Jacobian), preferring fewer free pieces, then less shortfall from 0.9, then + fewer pieces. A region's blocks are kept when fair (≥ 0.5) or no worse than its free mesh; cuts never + end inside a seeded curve. + **B8.7 Gordon and Thiel 1982, trivariate transfinite interpolation.** **Verdict: BUILD ON (mapped blocks).** - **Citation.** W. J. Gordon, L. C. Thiel. *Transfinite mappings and their application to grid generation.* Applied Mathematics and Computation 10–11:171–233, 1982.