diff --git a/crates/mesh/src/hex/cells.rs b/crates/mesh/src/hex/cells.rs index 2da773f..5c69914 100644 --- a/crates/mesh/src/hex/cells.rs +++ b/crates/mesh/src/hex/cells.rs @@ -36,8 +36,9 @@ use std::collections::BTreeMap; use super::brep::Brep; use super::geom::{TriangleTree, add, dist, scale}; use super::mapped::{self, Cube, face_ring, other}; -use super::{Cells, Ctx, HexMesh, HexMethod, HexOptions, RELATIVE_TOLERANCE, Raw, Surface, finish, project_to_face, project_to_faces}; +use super::{Cells, Ctx, HexMesh, HexMethod, HexOptions, Raw, Surface, finish, project_to_face, project_to_faces}; use crate::quad::{lattice_index, transfinite_grid, transfinite_volume}; +use crate::seeds::EdgeSeed; use crate::size::{SizeField, divide_polyline, ideal_divisions}; use crate::surface::BoundarySurface; use crate::vec3::{P, cross, dot, norm, sub}; @@ -112,8 +113,7 @@ pub fn hex_mesh_cells(cells: &[Cell], field: &dyn SizeField, opts: &HexOptions) super::check_estimate(&surfaces, field, opts)?; let breps: Vec = cells.iter().map(|c| Brep::composite(c.surface, c.faces)).collect::>()?; let trees: Vec = cells.iter().map(|c| TriangleTree::new(c.surface)).collect::>()?; - let ctxs: Vec = - breps.iter().zip(&trees).map(|(brep, tree)| Ctx { brep, tree, field, opts, tol: RELATIVE_TOLERANCE * brep.diagonal }).collect(); + let ctxs: Vec = breps.iter().zip(&trees).map(|(brep, tree)| Ctx::new(brep, tree, field, opts)).collect::>()?; let names: Vec<&str> = cells.iter().map(|c| c.name).collect(); let cubes: Vec = breps .iter() @@ -317,12 +317,24 @@ impl Topology { fn intervals(&self, ctxs: &[Ctx], cubes: &[Cube], field: &dyn SizeField) -> Result> { let polylines: Vec> = self.edges.iter().map(|e| self.polyline(ctxs, e)).collect::>()?; let goals: Vec = polylines.iter().map(|p| ideal_divisions(field, p)).collect::>()?; - let first = self.solve(&goals)?; + // The seed on each edge (a shared edge is one edge, so its seed holds in every cell it + // bounds), oriented from `a` to `b`. + let seeds: Vec> = self + .edges + .iter() + .map(|e| { + let (c, be, fwd) = e.owner; + ctxs.get(c).and_then(|x| x.seed_of(be)).map(|(seed, rev)| (seed, rev == fwd)) + }) + .collect(); + let first = self.solve(&goals, &seeds)?; // Ends of edges that leave a curved face: the spacing along that face there. let ends = self.graded_ends(ctxs, cubes, &first); let mut graded = goals; let mut divisions = Vec::with_capacity(polylines.len()); - for ((poly, ends), goal) in polylines.into_iter().zip(ends).zip(graded.iter_mut()) { + for (((poly, ends), goal), seed) in polylines.into_iter().zip(ends).zip(graded.iter_mut()).zip(&seeds) { + // A seeded edge keeps the seed's spacing, not a grading. + let ends = if seed.is_some() { [None, None] } else { ends }; let poly = if ends.iter().any(Option::is_some) { let dense = densify(&poly); *goal = goal.max(ideal_divisions(&Grading { field, ends }, &dense)?); @@ -330,18 +342,21 @@ impl Topology { } else { poly }; - divisions.push(Division { count: 1, ends, polyline: poly }); + divisions.push(Division { count: 1, ends, polyline: poly, seed: seed.map(|(s, rev)| (s.clone(), rev)) }); } - for (d, n) in divisions.iter_mut().zip(self.solve(&graded)?) { + for (d, n) in divisions.iter_mut().zip(self.solve(&graded, &seeds)?) { d.count = n; } Ok(divisions) } - fn solve(&self, goals: &[f64]) -> Result> { + fn solve(&self, goals: &[f64], seeds: &[Option<(&EdgeSeed, bool)>]) -> Result> { let mut ia = super::interval::Intervals::default(); - for &g in goals { - ia.var(g, 1); + for (&g, seed) in goals.iter().zip(seeds) { + let v = ia.var(g, 1); + if let Some((seed, _)) = seed { + ia.seed(v, seed); + } } for groups in &self.cell_edges { for g in groups { @@ -499,6 +514,8 @@ struct Division { count: usize, ends: GradedEnds, polyline: Vec

, + /// The edge's seed, and whether it runs from `b` to `a`. + seed: Option<(EdgeSeed, bool)>, } /// The size along a graded edge: the field's, capped near each graded end by the spacing there @@ -566,7 +583,9 @@ impl Builder<'_, '_, '_> { // Edges: divided once, on the geometry. for (edge, division) in topo.edges.iter().zip(self.divisions) { let n = division.count; - let pts = if division.ends.iter().any(Option::is_some) { + let pts = if let Some((seed, rev)) = &division.seed { + crate::seeds::place(&division.polyline, &seed.fractions(n, *rev))? + } else if division.ends.iter().any(Option::is_some) { divide_polyline(&Grading { field, ends: division.ends }, &division.polyline, n)? } else { divide_polyline(field, &division.polyline, n)? diff --git a/crates/mesh/src/hex/interval.rs b/crates/mesh/src/hex/interval.rs index cb604b7..1cb369c 100644 --- a/crates/mesh/src/hex/interval.rs +++ b/crates/mesh/src/hex/interval.rs @@ -37,6 +37,9 @@ pub(crate) struct Intervals { pub multiple: Vec<(Vec, usize)>, /// `Σ vars ≥ min`. pub at_least: Vec<(Vec, usize)>, + /// Variables held at their goal (an exact edge seed), with the seed in words; empty or + /// `None`: free. + pub fixed: Vec>, } /// The largest count a variable may reach. @@ -50,6 +53,43 @@ impl Intervals { self.goals.len() - 1 } + /// Apply an edge seed to variable `v`: an exact seed holds it at the seeded count, a minimum + /// raises its lower bound (and its goal), a free seed sets its goal. + pub(crate) fn seed(&mut self, v: usize, seed: &crate::seeds::EdgeSeed) { + use fieldcraft_model::SeedConstraint; + let n = seed.count.max(1); + let (Some(goal), Some(min)) = (self.goals.get_mut(v), self.mins.get_mut(v)) else { return }; + match seed.constraint { + SeedConstraint::Exact => { + *goal = n as f64; + *min = n; + if self.fixed.len() <= v { + self.fixed.resize(v + 1, None); + } + if let Some(slot) = self.fixed.get_mut(v) { + *slot = Some(seed.describe()); + } + } + SeedConstraint::Min => { + *min = (*min).max(n); + *goal = goal.max(n as f64); + } + SeedConstraint::Free => *goal = n as f64, + } + } + + fn is_fixed(&self, v: usize) -> bool { + self.fixed.get(v).is_some_and(Option::is_some) + } + + /// The conflict error: the exact seeds among `vars` cannot all keep their counts. + fn conflict(&self, vars: &[usize], what: &str) -> Error { + let mut seeds: Vec<&str> = vars.iter().filter_map(|&v| self.fixed.get(v).and_then(Option::as_deref)).collect(); + seeds.sort_unstable(); + seeds.dedup(); + Error::Invalid(format!("the exact edge seeds {} conflict: {what}; change a count or make a seed `min` or `free`", seeds.join(", "))) + } + pub(crate) fn solve(&self) -> Result> { let n = self.goals.len(); if self.mins.len() != n { @@ -110,13 +150,22 @@ impl Intervals { if cands.is_empty() { return Err(Error::Failed("interval assignment: a side constraint cannot be met (the same curves on both sides)".into())); } + let cands: Vec = cands.into_iter().filter(|&i| !self.is_fixed(i)).collect(); + if cands.is_empty() { + let all: Vec = a.iter().chain(b).copied().collect(); + return Err(self.conflict(&all, &format!("opposite sides of a mapped or swept face need equal counts ({sa} against {sb})"))); + } raise = Some(cands); break; } if raise.is_none() { for (l, m) in &self.at_least { if sum(l, &x) < *m { - raise = Some(l.clone()); + let free: Vec = l.iter().copied().filter(|&i| !self.is_fixed(i)).collect(); + if free.is_empty() { + return Err(self.conflict(l, &format!("a face's loop needs at least {m} divisions"))); + } + raise = Some(free); break; } } @@ -130,6 +179,10 @@ impl Intervals { if moving.is_empty() { return Err(Error::Failed("interval assignment: a loop's count cannot reach the multiple it needs".into())); } + let moving: Vec = moving.into_iter().filter(|&i| !self.is_fixed(i)).collect(); + if moving.is_empty() { + return Err(self.conflict(l, &format!("a face's loop needs a multiple of {m} divisions"))); + } raise = Some(moving); parity = true; break; @@ -205,6 +258,38 @@ mod tests { assert_eq!((x[a] + x[b] + x[c]) % 2, 0); } + #[test] + fn exact_seeds_hold_and_conflicts_name_them() { + use crate::seeds::EdgeSeed; + use fieldcraft_model::SeedConstraint; + let seed = + |name: &str, count: usize, constraint| EdgeSeed { name: name.into(), points: Vec::new(), tolerance: 0.0, count, bias: None, constraint }; + // The other side follows an exact seed. + let mut p = Intervals::default(); + let a = p.var(3.0, 1); + let b = p.var(3.0, 1); + p.equal.push((vec![a], vec![b])); + p.seed(a, &seed("A", 7, SeedConstraint::Exact)); + assert_eq!(p.solve().unwrap(), vec![7, 7]); + // Two exact seeds on opposite sides that differ. + p.seed(b, &seed("B", 5, SeedConstraint::Exact)); + let e = p.solve().unwrap_err().to_string(); + assert!(e.contains("A (7 exact)") && e.contains("B (5 exact)"), "{e}"); + // A minimum and a free goal. + let mut q = Intervals::default(); + let a = q.var(3.0, 1); + let b = q.var(3.0, 1); + q.seed(a, &seed("A", 6, SeedConstraint::Min)); + q.seed(b, &seed("B", 9, SeedConstraint::Free)); + assert_eq!(q.solve().unwrap(), vec![6, 9]); + // An exact odd count on a loop that must be even, with nothing else to move. + let mut r = Intervals::default(); + let a = r.var(3.0, 1); + r.multiple.push((vec![a], 2)); + r.seed(a, &seed("A", 3, SeedConstraint::Exact)); + assert!(r.solve().unwrap_err().to_string().contains("multiple of 2")); + } + #[test] fn infeasible_systems_are_errors_not_loops() { let mut p = Intervals::default(); diff --git a/crates/mesh/src/hex/layered.rs b/crates/mesh/src/hex/layered.rs index 7f81b40..955195c 100644 --- a/crates/mesh/src/hex/layered.rs +++ b/crates/mesh/src/hex/layered.rs @@ -34,9 +34,11 @@ use super::interval::Intervals; use super::region::{self, Scheme}; use super::{Ctx, Raw}; use crate::predicates::Sign; +use crate::seeds::EdgeSeed; use crate::size::{SizeField, divide_polyline, edge_divisions, ideal_divisions}; use crate::vec3::{P, cross, dot, norm, sub}; use crate::{Error, Result}; +use fieldcraft_model::SeedConstraint; /// Angular tolerance (radians) of "parallel" and "perpendicular" for the recognition. const ANGLE_TOL: f64 = 1e-6; @@ -262,6 +264,8 @@ pub(crate) struct Setup { roles: Vec, /// Sorted levels (heights or angles in `[0, 2π)`). levels: Vec, + /// Sweep coordinates closer than this are one level. + h_tol: f64, } /// Classify faces and edges and find the levels; an `Unsupported` error says why the body does @@ -322,7 +326,7 @@ pub(crate) fn setup(ctx: &Ctx, frame: &Frame) -> Result { if levels.len() < need { return Err(Error::Unsupported("no source and target faces across the sweep".into())); } - Ok(Setup { frame: *frame, roles, levels }) + Ok(Setup { frame: *frame, roles, levels, h_tol }) } impl Setup { @@ -452,6 +456,12 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { *slot = scheme; } } + let seeds = Seeds::of(ctx, frame, h_ref, intervals, setup.h_tol, &(0..arr.curves.len()).map(curve3d).collect::>())?; + for (c, seed) in seeds.curves.iter().enumerate() { + if let Some((seed, _)) = seed { + ia.seed(c, seed); + } + } let counts = ia.solve()?; let r_max = arr.points.iter().map(|q| q[1].abs()).fold(0.0, f64::max); let mut layers = Vec::with_capacity(intervals.len()); @@ -464,6 +474,13 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { if frame.periodic() { n = n.max(((hi - lo) / FRAC_PI_4).ceil() as usize); } + if let Some(Some(l)) = seeds.layers.get(layers.len()) { + let seeded = l.fractions.len().saturating_sub(1).max(1); + n = match l.seed.constraint { + SeedConstraint::Exact | SeedConstraint::Free => seeded, + SeedConstraint::Min => n.max(seeded), + }; + } layers.push(n); } // 2D nodes: curve ends, then curve interiors (projected onto analytic walls). @@ -472,7 +489,10 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { let mut curve_nodes: Vec> = Vec::with_capacity(arr.curves.len()); for (c, curve) in arr.curves.iter().enumerate() { let n = counts.get(c).copied().unwrap_or(1); - let pts3 = divide_polyline(ctx.field, &curve3d(c), n)?; + let pts3 = match seeds.curves.get(c).copied().flatten() { + Some((seed, rev)) => crate::seeds::place(&curve3d(c), &seed.fractions(n, rev))?, + None => divide_polyline(ctx.field, &curve3d(c), n)?, + }; let wall = curve.tags.iter().find_map(|&f| brep.faces.get(f).and_then(|x| x.surface)).filter(|s| !matches!(s, Surface::Plane { .. })); let mut ids = Vec::with_capacity(pts3.len()); let last = pts3.len().saturating_sub(1); @@ -548,8 +568,12 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { for (k, &(lo, hi)) in intervals.iter().enumerate() { first_layer.push(heights.len()); let n = layers.get(k).copied().unwrap_or(1); - for j in 0..n { - heights.push(lo + (hi - lo) * j as f64 / n as f64); + let fractions = match seeds.layers.get(k).and_then(Option::as_ref) { + Some(l) if l.fractions.len() == n + 1 => l.fractions.clone(), + _ => (0..=n).map(|j| j as f64 / n as f64).collect(), + }; + for f in fractions.iter().take(n) { + heights.push(lo + (hi - lo) * f); } } if !frame.periodic() @@ -614,6 +638,145 @@ pub(crate) fn mesh(ctx: &Ctx, prepared: &Prepared, refine: u32) -> Result { Ok(Raw { nodes, hexes, description }) } +/// The edge seeds as a layered mesh takes them: on the cross-section's curves (seeds of edges +/// across the sweep, projected into the cross-section) and on the level intervals (seeds of +/// edges along the sweep, which set the interval's layers). Each with whether it runs against +/// the curve, or from the interval's top down. +struct Seeds<'s> { + curves: Vec>, + layers: Vec>>, +} + +/// A seed's layers in one level interval: its nodes there, as fractions from the interval's +/// bottom (0) to its top (1). +#[derive(Clone)] +struct LayerSeed<'s> { + seed: &'s EdgeSeed, + fractions: Vec, +} + +/// Split a seed's node positions (fractions from the bottom of the span it seeds, 0 to 1) over +/// the level intervals in it (`bounds`: the levels from the bottom of the span to its top, as +/// fractions): each inner level takes the nearest node still free, and the nodes between two +/// levels are stretched to fit between them. Every interval gets at least one division; `None` +/// when there are fewer nodes than that needs. +fn split_over(nodes: &[f64], bounds: &[f64]) -> Option>> { + let (n, m) = (nodes.len().checked_sub(1)?, bounds.len().checked_sub(1)?); + if m == 0 || n < m { + return None; + } + // The node index at each level. + let mut at = vec![0usize; m + 1]; + at[m] = n; + for i in 1..m { + let b = bounds.get(i).copied()?; + let (lo, hi) = (at[i - 1] + 1, n - (m - i)); + at[i] = (lo..=hi).min_by(|&a, &c| (nodes[a] - b).abs().total_cmp(&(nodes[c] - b).abs()))?; + } + (0..m) + .map(|i| { + let (j0, j1) = (at[i], at[i + 1]); + let (z0, z1) = (nodes.get(j0).copied()?, nodes.get(j1).copied()?); + let span = z1 - z0; + Some( + (j0..=j1) + .map(|j| { + if j == j0 { + 0.0 + } else if j == j1 { + 1.0 + } else { + nodes.get(j).map_or(0.0, |z| (z - z0) / span) + } + }) + .collect(), + ) + }) + .collect() +} + +impl<'s> Seeds<'s> { + fn of(ctx: &'s Ctx, frame: &Frame, h_ref: f64, intervals: &[(f64, f64)], h_tol: f64, curves: &[Vec

]) -> Result> { + let brep = ctx.brep; + let tol = ctx.tol; + let mut out = Seeds { curves: vec![None; curves.len()], layers: vec![None; intervals.len()] }; + let stol = super::seed_tolerance(ctx.opts, tol); + let conflict = |a: &EdgeSeed, b: &EdgeSeed, why: &str| { + Error::Invalid(format!( + "the exact edge seeds {} and {} conflict: {why}; change a count or make a seed `min` or `free`", + a.describe(), + b.describe() + )) + }; + let unsupported = |seed: &EdgeSeed, why: &str| Error::Unsupported(format!("the exact edge seed {} {why}", seed.describe())); + for e in 0..brep.edges.len() { + let Some((seed, rev)) = ctx.seed_of(e) else { continue }; + let exact = seed.constraint == SeedConstraint::Exact; + let pts: Vec<(Q, f64)> = brep.edge_points(e, true).into_iter().map(|p| frame.coords(p)).collect(); + let (Some(&(_, h0)), Some(&(_, h1))) = (pts.first(), pts.last()) else { continue }; + let across = pts.iter().all(|&(q, h)| frame.span(q, frame.dh(h, h0)) <= tol); + if across { + // The seed moved into the cross-section plane at the reference level. + let projected = EdgeSeed { points: seed.points.iter().map(|&p| frame.point(frame.coords(p).0, h_ref)).collect(), ..seed.clone() }; + match curves.iter().position(|c| crate::seeds::matches(&projected, c, stol).is_some()) { + Some(c) => { + let r = crate::seeds::matches(&projected, curves.get(c).map(Vec::as_slice).unwrap_or(&[]), stol).unwrap_or(false); + let slot = out.curves.get_mut(c).ok_or_else(|| Error::Failed("curve index".into()))?; + if let Some((other, _)) = slot + && exact + && other.constraint == SeedConstraint::Exact + && other.count != seed.count + { + return Err(conflict(other, seed, "both lie on one curve of the swept cross-section")); + } + *slot = Some((seed, r)); + } + None if exact => return Err(unsupported(seed, "is not one curve of the swept cross-section")), + None => {} + } + continue; + } + // Along the sweep: the interval from its lower to its upper end. + let span = frame.dh(h1, h0); + let (lo, hi, down) = if span >= 0.0 { (h0, h0 + span, false) } else { (h0 + span, h0, true) }; + let same = |a: f64, b: f64| frame.dh(a, b).abs() <= h_tol; + let first = intervals.iter().position(|&(l, _)| same(l, lo)); + let last = intervals.iter().position(|&(_, h)| same(h, hi)); + let (Some(k0), Some(k1)) = (first, last) else { + if exact { + return Err(unsupported(seed, "does not run from one level of the sweep to another")); + } + continue; + }; + let spanned = intervals.get(k0..=k1).unwrap_or(&[]); + let length = frame.dh(hi, lo).abs().max(f64::MIN_POSITIVE); + let mut bounds: Vec = spanned.iter().map(|&(l, _)| frame.dh(l, lo) / length).collect(); + bounds.push(1.0); + // The seed's nodes from the bottom: the seed runs down when the edge does and the + // seed runs with it, or up against a downward edge. + let nodes = seed.fractions(seed.count, rev != down); + let Some(parts) = split_over(&nodes, &bounds) else { + if exact { + return Err(unsupported(seed, &format!("has fewer divisions than the {} level intervals it spans", spanned.len()))); + } + continue; + }; + for (k, fractions) in (k0..=k1).zip(parts) { + let slot = out.layers.get_mut(k).ok_or_else(|| Error::Failed("interval index".into()))?; + if let Some(other) = slot + && exact + && other.seed.constraint == SeedConstraint::Exact + && other.fractions.len() != fractions.len() + { + return Err(conflict(other.seed, seed, "both set the layers between the same two levels")); + } + *slot = Some(LayerSeed { seed, fractions }); + } + } + Ok(out) + } +} + /// Which level intervals each region of the arrangement fills (inside test at the region's /// centre, mid-interval), and which regions fill any. fn occupancy(arr: &Arrangement, frame: &Frame, intervals: &[(f64, f64)], inside: &mut Inside) -> Result<(Vec>, Vec)> { @@ -791,3 +954,22 @@ fn interior_point(arr: &Arrangement, region: &arrangement::Region) -> Result } best.filter(|(c, _)| *c > 0.0).map(|(_, p)| p).ok_or_else(|| Error::Failed("no point found inside a cross-section region".into())) } + +#[cfg(test)] +mod seed_tests { + use super::split_over; + + #[test] + fn a_seed_over_several_levels_is_split_at_them() { + // Six even divisions over levels at 0, 0.4 and 1: the level takes the node at 1/3. + let nodes: Vec = (0..=6).map(|j| j as f64 / 6.0).collect(); + let parts = split_over(&nodes, &[0.0, 0.4, 1.0]).unwrap(); + assert_eq!(parts.iter().map(|p| p.len() - 1).collect::>(), vec![2, 4]); + assert!(parts.iter().all(|p| p[0] == 0.0 && *p.last().unwrap() == 1.0 && p.windows(2).all(|w| w[0] < w[1]))); + // Every interval keeps at least one division; too few nodes is None. + let parts = split_over(&[0.0, 0.5, 1.0], &[0.0, 0.1, 0.2, 1.0]); + assert!(parts.is_none()); + let parts = split_over(&[0.0, 0.9, 0.95, 1.0], &[0.0, 0.1, 0.2, 1.0]).unwrap(); + assert_eq!(parts.iter().map(|p| p.len() - 1).collect::>(), vec![1, 1, 1]); + } +} diff --git a/crates/mesh/src/hex/mapped.rs b/crates/mesh/src/hex/mapped.rs index de79600..d7fbcca 100644 --- a/crates/mesh/src/hex/mapped.rs +++ b/crates/mesh/src/hex/mapped.rs @@ -16,7 +16,7 @@ use super::brep::Brep; use super::interval::Intervals; use super::{Ctx, Raw, project_to_face, project_to_faces}; use crate::quad::{lattice_index, transfinite_faces, transfinite_volume}; -use crate::size::{divide_polyline, ideal_divisions}; +use crate::size::ideal_divisions; use crate::vec3::{P, det, sub}; use crate::{Error, Result}; @@ -155,7 +155,7 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result { let groups = cube.groups; let mut ia = Intervals::default(); for e in 0..brep.edges.len() { - ia.var(ideal_divisions(ctx.field, &brep.edge_points(e, true))?, 1); + ctx.edge_var(&mut ia, e, ideal_divisions(ctx.field, &brep.edge_points(e, true))?); } for g in &groups { let first = g.first().map(|x| x.0).unwrap_or(0); @@ -173,7 +173,7 @@ pub(crate) fn mesh(ctx: &Ctx) -> Result { } // Edge node positions, low corner to high. let edge_pts = |(e, fwd): (usize, bool), n: usize| -> Result> { - let pts = divide_polyline(ctx.field, &brep.edge_points(e, fwd), n)?; + let pts = ctx.divide_edge(e, fwd, n)?; let faces = brep.edges.get(e).map(|x| x.faces.to_vec()).unwrap_or_default(); let last = pts.len().saturating_sub(1); Ok(pts.iter().enumerate().map(|(k, q)| if k == 0 || k == last { *q } else { project_to_faces(ctx, &faces, *q) }).collect()) diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index 5aa8db9..f278f6f 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -52,6 +52,8 @@ use std::collections::HashMap; pub use cells::{Cell, CellsMesh, GROWTH, MAX_CELLS, hex_mesh_cells}; pub use geom::Surface; +use crate::seeds::EdgeSeed; + use brep::Brep; use geom::{TriangleTree, add, dist, scale, unit}; @@ -109,11 +111,21 @@ pub struct HexOptions { pub max_elements: usize, /// Options of the face meshers. pub quad: QuadOptions, + /// Edge seeds: the B-rep edges they lie on get their counts (exact seeds are held, and a + /// method that cannot keep them names them) and their spacing. + pub seeds: Vec, } impl Default for HexOptions { fn default() -> Self { - HexOptions { method: HexMethod::Auto, quadratic: false, chordal_tolerance: 0.0, max_elements: 5_000_000, quad: QuadOptions::default() } + HexOptions { + method: HexMethod::Auto, + quadratic: false, + chordal_tolerance: 0.0, + max_elements: 5_000_000, + quad: QuadOptions::default(), + seeds: Vec::new(), + } } } @@ -176,6 +188,49 @@ pub(crate) struct Ctx<'a> { pub opts: &'a HexOptions, /// Classification tolerance (model units). pub tol: f64, + /// The seed on each B-rep edge (index into `opts.seeds`, and whether the edge runs from the + /// seed's end to its start). + pub edge_seeds: Vec>, +} + +impl<'a> Ctx<'a> { + pub(crate) fn new(brep: &'a Brep<'a>, tree: &'a TriangleTree<'a>, field: &'a dyn SizeField, opts: &'a HexOptions) -> Result> { + let tol = RELATIVE_TOLERANCE * brep.diagonal; + let curves: Vec> = (0..brep.edges.len()).map(|e| brep.edge_points(e, true)).collect(); + let edge_seeds = crate::seeds::match_curves(&opts.seeds, &curves, seed_tolerance(opts, tol))?; + Ok(Ctx { brep, tree, field, opts, tol, edge_seeds }) + } + + /// The seed on B-rep edge `e`, and whether it runs against the edge. + pub(crate) fn seed_of(&self, e: usize) -> Option<(&EdgeSeed, bool)> { + let (i, rev) = self.edge_seeds.get(e).copied().flatten()?; + Some((self.opts.seeds.get(i)?, rev)) + } + + /// A division-count variable for B-rep edge `e` with goal `goal`, its seed applied. + pub(crate) fn edge_var(&self, ia: &mut interval::Intervals, e: usize, goal: f64) -> usize { + let v = ia.var(goal, 1); + if let Some((seed, _)) = self.seed_of(e) { + ia.seed(v, seed); + } + v + } + + /// `n + 1` points along B-rep edge `e` (traversed forwards or not): at the seed's spacing + /// when it has one, else evenly in the size field. + pub(crate) fn divide_edge(&self, e: usize, forward: bool, n: usize) -> Result> { + let pts = self.brep.edge_points(e, forward); + match self.seed_of(e) { + Some((seed, rev)) => crate::seeds::place(&pts, &seed.fractions(n, rev == forward)), + None => crate::size::divide_polyline(self.field, &pts, n), + } + } +} + +/// How far apart a seed's polyline and a mesher's curve of the same edge may lie: both are +/// tessellations of it. +pub(crate) fn seed_tolerance(opts: &HexOptions, tol: f64) -> f64 { + 4.0 * opts.chordal_tolerance.max(0.0) + 16.0 * tol } /// What a method produces: linear hexahedra in the model's node order. @@ -208,7 +263,7 @@ pub fn hex_mesh(surface: &BoundarySurface, faces: &[Option], field: &dy check_estimate(&[surface], field, opts)?; let brep = Brep::new(surface, faces)?; let tree = TriangleTree::new(surface)?; - let ctx = Ctx { brep: &brep, tree: &tree, field, opts, tol: RELATIVE_TOLERANCE * brep.diagonal }; + let ctx = Ctx::new(&brep, &tree, field, opts)?; let methods: &[HexMethod] = match opts.method { HexMethod::Auto => &[HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve, HexMethod::Mapped], HexMethod::Swept => &[HexMethod::Extrude, HexMethod::Sweep, HexMethod::Revolve], diff --git a/crates/mesh/src/hex/sweep.rs b/crates/mesh/src/hex/sweep.rs index 1eaa38a..864dfe5 100644 --- a/crates/mesh/src/hex/sweep.rs +++ b/crates/mesh/src/hex/sweep.rs @@ -32,7 +32,7 @@ use super::layered::{Frame, PlaneField}; use super::region; use super::{Ctx, Raw, project_to_face, project_to_faces}; use crate::quad::transfinite_grid; -use crate::size::{divide_polyline, ideal_divisions}; +use crate::size::ideal_divisions; use crate::vec3::{P, sub}; use crate::{Error, Result}; @@ -215,7 +215,7 @@ fn build(ctx: &Ctx, topo: &Topology, refine: u32) -> Result { let fine = f64::from(1u32 << refine.min(8)); for e in 0..brep.edges.len() { let goal = ideal_divisions(ctx.field, &brep.edge_points(e, true))?; - ia.var(if sides.contains(&e) { goal } else { goal * fine }, 1); + ctx.edge_var(&mut ia, e, if sides.contains(&e) { goal } else { goal * fine }); } let vars = |uses: &[(usize, bool)]| uses.iter().map(|x| x.0).collect::>(); for stages in &topo.rings { @@ -257,7 +257,7 @@ fn build(ctx: &Ctx, topo: &Topology, refine: u32) -> Result { let mut edge_nodes: Vec> = Vec::with_capacity(brep.edges.len()); for (e, edge) in brep.edges.iter().enumerate() { let n = counts.get(e).copied().unwrap_or(1); - let pts = divide_polyline(ctx.field, &brep.edge_points(e, true), n)?; + let pts = ctx.divide_edge(e, true, n)?; let last = pts.len().saturating_sub(1); let mut ids = Vec::with_capacity(pts.len()); for (k, p) in pts.iter().enumerate() { diff --git a/crates/mesh/src/hex/tests.rs b/crates/mesh/src/hex/tests.rs index 877780f..28d34ff 100644 --- a/crates/mesh/src/hex/tests.rs +++ b/crates/mesh/src/hex/tests.rs @@ -480,3 +480,105 @@ fn an_o_grid_round_a_hole_conforms_with_graded_radial_lines() { assert!(r.windows(3).all(|w| (w[2] - w[1]) <= 1.35 * (w[1] - w[0])), "{r:?}"); } } + +fn seed( + name: &str, + points: Vec<[f64; 3]>, + count: usize, + bias: Option, + constraint: fieldcraft_model::SeedConstraint, +) -> crate::seeds::EdgeSeed { + crate::seeds::EdgeSeed { name: name.into(), points, tolerance: 0.0, count, bias, constraint } +} + +/// The distinct values of coordinate `k` of the nodes `on` accepts, sorted. +fn along(m: &HexMesh, k: usize, on: impl Fn([f64; 3]) -> bool) -> Vec { + let mut v: Vec = m.nodes.iter().filter(|p| on(**p)).map(|p| p[k]).collect(); + v.sort_by(f64::total_cmp); + v.dedup_by(|a, b| (*a - *b).abs() < 1e-9); + v +} + +#[test] +fn seeds_give_a_swept_block_exact_counts_and_geometric_spacing() { + use fieldcraft_model::{Bias, BiasToward, SeedConstraint}; + let s = prism(&[[0.0, 0.0], [4.0, 0.0], [4.0, 2.0], [0.0, 2.0]], 0.0, 3.0); + let bias = Some(Bias { ratio: 4.0, toward: BiasToward::Start }); + let opts = HexOptions { + seeds: vec![ + // Through the thickness, 7 layers small at z = 0 (the seed runs up). + seed("thickness", vec![[0.0, 0.0, 0.0], [0.0, 0.0, 3.0]], 7, bias, SeedConstraint::Exact), + // Along x on the top: 5 divisions, the seed runs from x = 4 to 0. + seed("length", vec![[4.0, 2.0, 3.0], [0.0, 2.0, 3.0]], 5, None, SeedConstraint::Exact), + ], + ..HexOptions::default() + }; + let m = hex_mesh(&s, &[], &UniformSize(1.0), &opts).unwrap(); + check(&m, 24.0); + assert_eq!(m.element_count(), 5 * 2 * 7, "{}", m.report.description); + // The nodes on the seeded edges. + assert_eq!(along(&m, 0, |p| (p[1] - 2.0).abs() < 1e-9 && (p[2] - 3.0).abs() < 1e-9).len(), 6); + let z = along(&m, 2, |p| p[0].abs() < 1e-9 && p[1].abs() < 1e-9); + assert_eq!(z.len(), 8, "{z:?}"); + // Element heights grow geometrically from the bottom, the last 4 times the first. + let dz: Vec = z.windows(2).map(|w| w[1] - w[0]).collect(); + let q = 4f64.powf(1.0 / 6.0); + assert!(dz.windows(2).all(|w| (w[1] / w[0] - q).abs() < 1e-9), "{dz:?}"); + // Seeded the other way round: small at the top. + let down = HexOptions { seeds: vec![seed("t", vec![[0.0, 0.0, 3.0], [0.0, 0.0, 0.0]], 7, bias, SeedConstraint::Exact)], ..opts.clone() }; + let m = hex_mesh(&s, &[], &UniformSize(1.0), &down).unwrap(); + let z = along(&m, 2, |p| p[0].abs() < 1e-9 && p[1].abs() < 1e-9); + assert!(z[7] - z[6] < z[1] - z[0], "{z:?}"); + // The one-to-one sweep and the mapped mesher keep the same seeds. + for method in [HexMethod::Sweep, HexMethod::Mapped] { + let m = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { method, ..opts.clone() }).unwrap(); + let z = along(&m, 2, |p| p[0].abs() < 1e-9 && p[1].abs() < 1e-9); + let dz: Vec = z.windows(2).map(|w| w[1] - w[0]).collect(); + assert!(dz.len() == 7 && dz.windows(2).all(|w| (w[1] / w[0] - q).abs() < 1e-9), "{method:?} {dz:?}"); + assert_eq!(m.element_count(), 5 * 2 * 7, "{method:?}"); + } + // A minimum raises a count; a free seed is a goal the other constraints may move (its + // opposite edges want 4, so 6 holds and 2 would not). + let min = HexOptions { seeds: vec![seed("x", vec![[0.0, 0.0, 0.0], [4.0, 0.0, 0.0]], 2, None, SeedConstraint::Min)], ..opts.clone() }; + let front = |p: [f64; 3]| p[1].abs() < 1e-9 && p[2].abs() < 1e-9; + assert_eq!(along(&hex_mesh(&s, &[], &UniformSize(1.0), &min).unwrap(), 0, front).len(), 5); + let free = HexOptions { seeds: vec![seed("x", vec![[0.0, 0.0, 0.0], [4.0, 0.0, 0.0]], 6, None, SeedConstraint::Free)], ..opts.clone() }; + assert_eq!(along(&hex_mesh(&s, &[], &UniformSize(1.0), &free).unwrap(), 0, front).len(), 7); +} + +#[test] +fn conflicting_exact_seeds_are_refused_with_their_names() { + use fieldcraft_model::SeedConstraint; + let s = prism(&[[0.0, 0.0], [4.0, 0.0], [4.0, 2.0], [0.0, 2.0]], 0.0, 3.0); + // Opposite edges of the mapped block (and of the swept cross-section) with different counts. + let opts = HexOptions { + method: HexMethod::Mapped, + seeds: vec![ + seed("front", vec![[0.0, 0.0, 0.0], [4.0, 0.0, 0.0]], 5, None, SeedConstraint::Exact), + seed("back", vec![[0.0, 2.0, 3.0], [4.0, 2.0, 3.0]], 3, None, SeedConstraint::Exact), + ], + ..HexOptions::default() + }; + let e = hex_mesh(&s, &[], &UniformSize(1.0), &opts).unwrap_err().to_string(); + assert!(e.contains("front (5 exact)") && e.contains("back (3 exact)"), "{e}"); + let e = hex_mesh(&s, &[], &UniformSize(1.0), &HexOptions { method: HexMethod::Auto, ..opts }).unwrap_err().to_string(); + assert!(e.contains("front (5 exact)") && e.contains("back (3 exact)"), "{e}"); +} + +#[test] +fn a_seed_on_a_shared_edge_holds_in_every_cell() { + use fieldcraft_model::{Bias, BiasToward, SeedConstraint}; + let pairs = boxes(&[([0.0, 2.0], [0.0, 1.0]), ([2.0, 3.0], [0.0, 1.0])], 1.0); + let none = vec![Vec::new(); 2]; + // The edge both cells share, at x = 2, along z. + let bias = Some(Bias { ratio: 3.0, toward: BiasToward::Both }); + let opts = + HexOptions { seeds: vec![seed("shared", vec![[2.0, 0.0, 0.0], [2.0, 0.0, 1.0]], 6, bias, SeedConstraint::Exact)], ..HexOptions::default() }; + let m = hex_mesh_cells(&cells_of(&pairs, &none), &UniformSize(0.5), &opts).unwrap(); + assert!(m.blocks.iter().all(|b| b[2] == 6), "{:?}", m.blocks); + assert!(distinct_positions(&m.mesh)); + let z = along(&m.mesh, 2, |p| (p[0] - 2.0).abs() < 1e-9 && p[1].abs() < 1e-9); + assert_eq!(z.len(), 7, "{z:?}"); + // Small at both ends, symmetric. + assert!((z[1] - z[0] - (z[6] - z[5])).abs() < 1e-9 && z[3] - z[2] > z[1] - z[0], "{z:?}"); +} diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index a4cdc72..685dc7b 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -17,6 +17,7 @@ pub mod octree; pub mod planar; pub mod quad; pub mod quality; +pub mod seeds; pub mod size; pub mod surface; pub mod tet; diff --git a/crates/mesh/src/seeds.rs b/crates/mesh/src/seeds.rs new file mode 100644 index 0000000..6af654a --- /dev/null +++ b/crates/mesh/src/seeds.rs @@ -0,0 +1,209 @@ +//! Edge seeds as the meshers take them: a design edge's polyline with the count, bias and +//! constraint the model's seed gives it (`fieldcraft_model::seeds`). +//! +//! The meshers recover their own curves from the face-tagged boundary surface, so a seed is +//! matched to a mesher's curve by position: the curve's end points lie on the seed's end points +//! (either way round, which gives the orientation the bias is measured in) and each one's middle +//! lies on the other, within the tolerance of the two tessellations. The seed's node positions +//! are placed by arc length along the mesher's curve ([`place`]). + +use fieldcraft_model::{Bias, SeedConstraint, seeds::spacing}; + +use crate::vec3::{P, dist, lerp, sub}; +use crate::{Error, Result}; + +/// A seed on one design edge. +#[derive(Debug, Clone, PartialEq)] +pub struct EdgeSeed { + /// The edge's name, for messages. + pub name: String, + /// The edge's polyline, start to end. + pub points: Vec<[f64; 3]>, + /// How far the polyline may lie from the exact edge (its tessellation's chordal tolerance). + pub tolerance: f64, + /// The divisions the seed asks for. + pub count: usize, + pub bias: Option, + pub constraint: SeedConstraint, +} + +impl EdgeSeed { + /// `n + 1` node positions as fractions of the arc length from the seed's start (reversed: + /// from its end). + pub fn fractions(&self, n: usize, reversed: bool) -> Vec { + let f = spacing(u32::try_from(n).unwrap_or(u32::MAX), self.bias); + if reversed { f.iter().rev().map(|x| 1.0 - x).collect() } else { f } + } + + /// The seed in words: `F2:a|F2:b (8 exact)`. + pub fn describe(&self) -> String { + format!("{} ({} {})", self.name, self.count, self.constraint.name()) + } +} + +fn length(poly: &[P]) -> f64 { + poly.windows(2).map(|w| dist(w[0], w[1])).sum() +} + +/// The point at fraction `t` of a polyline's arc length. +fn at_fraction(poly: &[P], t: f64) -> Option

{ + let total = length(poly); + let target = t.clamp(0.0, 1.0) * total; + let mut walked = 0.0; + for w in poly.windows(2) { + let l = dist(w[0], w[1]); + if walked + l >= target && l > 0.0 { + return Some(lerp(w[0], w[1], ((target - walked) / l).clamp(0.0, 1.0))); + } + walked += l; + } + poly.last().copied() +} + +fn point_segment(p: P, a: P, b: P) -> f64 { + let ab = sub(b, a); + let l2 = crate::vec3::dot(ab, ab); + let t = if l2 > 0.0 { (crate::vec3::dot(sub(p, a), ab) / l2).clamp(0.0, 1.0) } else { 0.0 }; + dist(p, lerp(a, b, t)) +} + +fn point_polyline(p: P, poly: &[P]) -> f64 { + match poly { + [] => f64::INFINITY, + [a] => dist(p, *a), + _ => poly.windows(2).map(|w| point_segment(p, w[0], w[1])).fold(f64::INFINITY, f64::min), + } +} + +/// Whether polyline `curve` is the seed's edge within `tol`: `Some(reversed)` when it runs from +/// the seed's end to its start. +pub(crate) fn matches(seed: &EdgeSeed, curve: &[P], tol: f64) -> Option { + let tol = tol + seed.tolerance; + let (Some(&a), Some(&b), Some(&s0), Some(&s1)) = (curve.first(), curve.last(), seed.points.first(), seed.points.last()) else { + return None; + }; + let forward = dist(a, s0) <= tol && dist(b, s1) <= tol; + let backward = dist(a, s1) <= tol && dist(b, s0) <= tol; + if !(forward || backward) { + return None; + } + // A closed edge (both ends at one point) is oriented by its first quarter. + let reversed = if forward && backward { + let q = at_fraction(curve, 0.25)?; + let (f, r) = (at_fraction(&seed.points, 0.25)?, at_fraction(&seed.points, 0.75)?); + dist(q, r) < dist(q, f) + } else { + backward + }; + let (cm, sm) = (at_fraction(curve, 0.5)?, at_fraction(&seed.points, 0.5)?); + (point_polyline(cm, &seed.points) <= tol && point_polyline(sm, curve) <= tol).then_some(reversed) +} + +/// The seed (index, reversed) of each curve, the first that matches; two seeds of different +/// counts on one curve are an error naming both. +pub(crate) fn match_curves(seeds: &[EdgeSeed], curves: &[Vec

], tol: f64) -> Result>> { + let mut out = vec![None; curves.len()]; + if seeds.is_empty() { + return Ok(out); + } + for (c, curve) in curves.iter().enumerate() { + for (i, s) in seeds.iter().enumerate() { + let Some(rev) = matches(s, curve, tol) else { continue }; + match out.get(c).copied().flatten() { + None => { + if let Some(slot) = out.get_mut(c) { + *slot = Some((i, rev)); + } + } + Some((j, _)) => { + if let Some(o) = seeds.get(j) + && (o.count, o.bias, o.constraint) != (s.count, s.bias, s.constraint) + { + return Err(Error::Invalid(format!("the edge seeds {} and {} are on the same edge", o.describe(), s.describe()))); + } + } + } + } + } + Ok(out) +} + +/// The points at the given arc-length fractions along `poly` (its end points exact). +pub(crate) fn place(poly: &[P], fractions: &[f64]) -> Result> { + let (Some(&first), Some(&last)) = (poly.first(), poly.last()) else { + return Err(Error::Invalid("empty polyline".into())); + }; + let n = fractions.len(); + fractions + .iter() + .enumerate() + .map(|(k, &t)| { + if k == 0 { + Ok(first) + } else if k + 1 == n { + Ok(last) + } else { + at_fraction(poly, t).ok_or_else(|| Error::Failed("seed spacing walked off the curve".into())) + } + }) + .collect() +} + +/// The element size a seed asks for at each of its nodes (the mean length of the elements on +/// either side), for size fields. +pub fn node_sizes(seed: &EdgeSeed) -> Vec<(P, f64)> { + let f = seed.fractions(seed.count, false); + let total = length(&seed.points); + let Ok(pts) = place(&seed.points, &f) else { return Vec::new() }; + pts.iter() + .enumerate() + .filter_map(|(k, p)| { + let before = k.checked_sub(1).and_then(|j| Some(f.get(k)? - f.get(j)?)); + let after = f.get(k + 1).zip(f.get(k)).map(|(b, a)| b - a); + let h = match (before, after) { + (Some(x), Some(y)) => 0.5 * (x + y), + (Some(x), None) | (None, Some(x)) => x, + (None, None) => return None, + } * total; + (h > 0.0 && h.is_finite()).then_some((*p, h)) + }) + .collect() +} + +#[cfg(test)] +mod tests { + use fieldcraft_model::BiasToward; + + use super::*; + + fn seed(points: Vec

) -> EdgeSeed { + EdgeSeed { + name: "e".into(), + points, + tolerance: 0.0, + count: 3, + bias: Some(Bias { ratio: 4.0, toward: BiasToward::Start }), + constraint: SeedConstraint::Exact, + } + } + + #[test] + fn curves_match_either_way_round_and_carry_the_orientation() { + let s = seed(vec![[0.0; 3], [5.0, 0.0, 0.0], [10.0, 0.0, 0.0]]); + assert_eq!(matches(&s, &[[0.0; 3], [10.0, 0.0, 0.0]], 1e-6), Some(false)); + assert_eq!(matches(&s, &[[10.0, 0.0, 0.0], [0.0; 3]], 1e-6), Some(true)); + // Same ends, another path. + assert_eq!(matches(&s, &[[0.0; 3], [5.0, 3.0, 0.0], [10.0, 0.0, 0.0]], 1e-6), None); + assert_eq!(matches(&s, &[[0.0; 3], [9.0, 0.0, 0.0]], 1e-6), None); + // Reversed fractions mirror the spacing. + let f = s.fractions(3, false); + let r = s.fractions(3, true); + assert!(f.iter().zip(r.iter().rev()).all(|(a, b)| (a - (1.0 - b)).abs() < 1e-12)); + let pts = place(&[[0.0; 3], [10.0, 0.0, 0.0]], &f).unwrap(); + assert!((pts[1][0] - 10.0 * f[1]).abs() < 1e-12); + // Sizes: the mean of the neighbouring elements. + let h = node_sizes(&s); + assert_eq!(h.len(), 4); + assert!((h[0].1 - 10.0 * f[1]).abs() < 1e-9); + } +}