diff --git a/crates/mesh/src/seeds.rs b/crates/mesh/src/seeds.rs index 6af654a..0101582 100644 --- a/crates/mesh/src/seeds.rs +++ b/crates/mesh/src/seeds.rs @@ -67,6 +67,11 @@ fn point_segment(p: P, a: P, b: P) -> f64 { dist(p, lerp(a, b, t)) } +/// The distance from `p` to a polyline. +pub(crate) fn distance_to(poly: &[P], p: P) -> f64 { + point_polyline(p, poly) +} + fn point_polyline(p: P, poly: &[P]) -> f64 { match poly { [] => f64::INFINITY, diff --git a/crates/mesh/src/tet/mod.rs b/crates/mesh/src/tet/mod.rs index 3d5614d..b96b6ed 100644 --- a/crates/mesh/src/tet/mod.rs +++ b/crates/mesh/src/tet/mod.rs @@ -93,6 +93,10 @@ pub struct QualityOptions { pub order: Order, pub max_tets: usize, pub max_points: usize, + /// Edge seeds: the size at each seeded node is the seed's spacing there, and the size field + /// grows from it at the grading rate. A seed finer than the global size refines along its + /// edge; one coarser leaves it to the global size. + pub seeds: Vec, } impl Default for QualityOptions { @@ -112,6 +116,7 @@ impl Default for QualityOptions { order: Order::Linear, max_tets: 5_000_000, max_points: 1_000_000, + seeds: Vec::new(), } } } @@ -197,6 +202,9 @@ pub struct QualityReport { pub straightened: usize, /// Why a stage stopped early, if one did. pub stopped: Vec, + /// Per edge seed: its name, the count it asks for and the boundary edges the mesh has along + /// it. + pub seeds: Vec<(String, usize, usize)>, } /// The quality mesh. @@ -246,9 +254,94 @@ pub fn size_field(surface: &BoundarySurface, geometry: &Tessellation, global: f6 } } } + for seed in &o.seeds { + for (p, h) in crate::seeds::node_sizes(seed) { + let h = h.max(floor); + if h < global { + samples.push((p, h)); + } + } + } GradedSize::new(global, &samples, lo, hi, o.grading) } +/// Put the seeds' nodes on their feature curves, pinned (corner vertices, which neither +/// smoothing nor coarsening moves), so the curve is divided where the seed says. A node within +/// [`refine::MIN_GAP`] of the local size of a vertex already on the curve pins that vertex +/// instead. +fn seed_curves(cx: &mut TetComplex, ctx: &Context, seeds: &[crate::seeds::EdgeSeed], global: f64) -> Result<()> { + for seed in seeds { + let tol = 2.0 * seed.tolerance + 1e-6 * global; + let on = |x: [f64; 3]| crate::seeds::distance_to(&seed.points, x) <= tol; + let nodes = crate::seeds::place(&seed.points, &seed.fractions(seed.count, false))?; + for p in nodes.iter().skip(1).take(nodes.len().saturating_sub(2)) { + // The curve edge along the seed nearest the node. + let mut best: Option<(f64, u32, u32, u32)> = None; + for ((a, b), c) in cx.curve_edges() { + let (xa, xb) = (cx.xyz(a)?, cx.xyz(b)?); + if !(on(xa) && on(xb)) { + continue; + } + let d = crate::seeds::distance_to(&[xa, xb], *p); + if best.is_none_or(|x| d < x.0) { + best = Some((d, a, b, c)); + } + } + let Some((_, a, b, c)) = best else { continue }; + let gap = refine::MIN_GAP * ctx.h(*p, global); + let near = [a, b].into_iter().map(|v| Ok((measure::dist(cx.xyz(v)?, *p), v))).collect::>>()?; + if let Some(&(_, v)) = near.iter().filter(|x| x.0 <= gap).min_by(|x, y| x.0.total_cmp(&y.0)) { + if cx.kind(v) != VertexKind::Corner { + cx.set_kind(v, VertexKind::Corner); + } + continue; + } + let Some(q) = ctx.place(VertexKind::Curve(c), *p) else { continue }; + if let Some(v) = ops::insert_delaunay_apart(cx, q, VertexKind::Curve(c), ops::Site::Edge(a, b), refine::MAX_CAVITY, gap)? { + cx.set_kind(v, VertexKind::Corner); + } + } + } + Ok(()) +} + +/// How many boundary edges of a mesh lie along each seed's edge: both ends on its polyline, and +/// the middle within a fifth of the edge's length of it (an element edge on a curved edge cuts +/// its chord). +fn seed_counts(points: &[[f64; 3]], tets: &[[u32; 4]], boundary: &[TaggedFace], seeds: &[crate::seeds::EdgeSeed]) -> Vec<(String, usize, usize)> { + if seeds.is_empty() { + return Vec::new(); + } + let mut edges: Vec<(u32, u32)> = Vec::with_capacity(boundary.len() * 3); + for b in boundary { + let Some(t) = tets.get(b.tet as usize) else { continue }; + let f = face_toward(*t, usize::from(b.opposite)); + for (a, c) in [(f[0], f[1]), (f[1], f[2]), (f[2], f[0])] { + edges.push(sorted2(a, c)); + } + } + edges.sort_unstable(); + edges.dedup(); + let (lo, hi) = bounds(points); + let eps = 1e-6 * measure::dist(lo, hi); + seeds + .iter() + .map(|seed| { + let tol = 2.0 * seed.tolerance + eps; + let on = |p: [f64; 3], slack: f64| crate::seeds::distance_to(&seed.points, p) <= tol + slack; + let n = edges + .iter() + .filter(|(a, c)| { + let (Some(&pa), Some(&pc)) = (points.get(*a as usize), points.get(*c as usize)) else { return false }; + let m = [(pa[0] + pc[0]) / 2.0, (pa[1] + pc[1]) / 2.0, (pa[2] + pc[2]) / 2.0]; + on(pa, 0.0) && on(pc, 0.0) && on(m, 0.2 * measure::dist(pa, pc)) + }) + .count(); + (seed.name.clone(), seed.count, n) + }) + .collect() +} + /// The CDT as a [`TetComplex`] with B-rep face tags, feature curves and vertex kinds; and the /// feature curves' face pairs. pub fn complex_from_cdt(surface: &BoundarySurface, cdt: &TetMesh, o: &QualityOptions) -> Result<(TetComplex, Vec<[u32; 2]>)> { @@ -391,6 +484,7 @@ pub fn quality_mesh(surface: &BoundarySurface, cdt: &TetMesh, analytic: Vec = c.boundary.iter().map(|&(tet, k, face)| TaggedFace { tet, opposite: k as u8, face }).collect(); report.output = tet_quality_stats(&c.points, &c.tets)?; + report.seeds = seed_counts(&c.points, &c.tets, &boundary, &o.seeds); report.boundary_volume = boundary_volume(&c.points, &c.tets, &boundary)?; let gate = linear_gate(&c.points, &c.tets, report.boundary_volume)?; report.volume = gate.volume; diff --git a/crates/mesh/src/tet/tests.rs b/crates/mesh/src/tet/tests.rs index ed2fc78..1986b17 100644 --- a/crates/mesh/src/tet/tests.rs +++ b/crates/mesh/src/tet/tests.rs @@ -195,6 +195,35 @@ mod pipeline { assert_eq!(a, b); } + #[test] + fn a_seeded_edge_gets_the_seeded_spacing_and_the_size_grows_from_it() { + use crate::seeds::EdgeSeed; + use fieldcraft_model::SeedConstraint; + let s = box_surface(); + let cdt = tetrahedralize(&s, &CdtOptions::default()).unwrap(); + let edge = vec![[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]; + let seed = EdgeSeed { name: "front".into(), points: edge, tolerance: 0.0, count: 16, bias: None, constraint: SeedConstraint::Exact }; + let plain = QualityOptions { size: Some(5.0), seeds: vec![EdgeSeed { count: 1, ..seed.clone() }], ..QualityOptions::default() }; + let seeded = QualityOptions { seeds: vec![seed], ..plain.clone() }; + let a = quality_mesh(&s, &cdt, Vec::new(), &plain).unwrap(); + let b = quality_mesh(&s, &cdt, Vec::new(), &seeded).unwrap(); + let (na, nb) = (a.report.seeds[0].2, b.report.seeds[0].2); + // Unseeded, the 10 mm edge gets about two 5 mm elements; seeded, about sixteen. + assert!(na <= 4, "{na}"); + assert_eq!(nb, 16); + assert_eq!(b.report.seeds[0].1, 16); + assert!(b.report.output.tets > a.report.output.tets && b.report.output.inadmissible == 0); + // Biased: the mesh's nodes on the edge are the seed's. + let bias = Some(fieldcraft_model::Bias { ratio: 4.0, toward: fieldcraft_model::BiasToward::Start }); + let biased = QualityOptions { seeds: vec![EdgeSeed { count: 8, bias, ..seeded.seeds[0].clone() }], ..plain.clone() }; + let c = quality_mesh(&s, &cdt, Vec::new(), &biased).unwrap(); + let mut x: Vec = c.points.iter().filter(|p| p[1].abs() < 1e-9 && p[2].abs() < 1e-9).map(|p| p[0]).collect(); + x.sort_by(f64::total_cmp); + let want: Vec = biased.seeds[0].fractions(8, false).iter().map(|f| 10.0 * f).collect(); + assert_eq!(x.len(), want.len(), "{x:?}"); + assert!(x.iter().zip(&want).all(|(a, b)| (a - b).abs() < 1e-9), "{x:?} {want:?}"); + } + #[test] fn a_sphere_keeps_its_boundary_on_the_sphere_to_second_order() { let (p, t) = uv_sphere(10, 20);