diff --git a/crates/mesh/src/size.rs b/crates/mesh/src/size.rs index 1448cde..907e35d 100644 --- a/crates/mesh/src/size.rs +++ b/crates/mesh/src/size.rs @@ -206,6 +206,150 @@ pub fn divide_polyline(field: &dyn SizeField, polyline: &[[f64; 3]], n: usize) - Ok(out) } +/// The most triangles the local sizes of one field hold. +pub const MAX_REGION_TRIANGLES: usize = 1 << 20; + +/// Where a [`LocalSize`] applies. +#[derive(Debug, Clone, PartialEq)] +pub enum Region { + /// An axis-aligned box. + Box { lo: [f64; 3], hi: [f64; 3] }, + /// A ball. + Sphere { center: [f64; 3], radius: f64 }, + /// A surface, as triangles (a B-rep face's tessellation). + Triangles(Vec<[[f64; 3]; 3]>), +} + +impl Region { + /// The distance from `p` to the region (0 inside a box or a ball, and on a surface). + pub fn distance(&self, p: [f64; 3]) -> f64 { + match self { + Region::Box { lo, hi } => { + let d = [0, 1, 2].map(|k| (lo[k] - p[k]).max(p[k] - hi[k]).max(0.0)); + (d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt() + } + Region::Sphere { center, radius } => (dist(p, *center) - radius).max(0.0), + Region::Triangles(tris) => { + tris.iter().map(|[a, b, c]| dist(p, crate::tet::geometry::closest_on_triangle(p, *a, *b, *c))).fold(f64::INFINITY, f64::min) + } + } + } + + /// The region's bounding box. + pub fn bounds(&self) -> ([f64; 3], [f64; 3]) { + match self { + Region::Box { lo, hi } => (*lo, *hi), + Region::Sphere { center, radius } => (center.map(|c| c - radius), center.map(|c| c + radius)), + Region::Triangles(tris) => { + let mut lo = [f64::INFINITY; 3]; + let mut hi = [f64::NEG_INFINITY; 3]; + for p in tris.iter().flatten() { + for k in 0..3 { + lo[k] = lo[k].min(p[k]); + hi[k] = hi[k].max(p[k]); + } + } + (lo, hi) + } + } + } +} + +/// A local size: elements of about `size` in `region`. +#[derive(Debug, Clone, PartialEq)] +pub struct LocalSize { + pub region: Region, + pub size: f64, +} + +/// A global size refined by local sizes and graded: the size at `p` is +/// `min(global, min_i (h_i + g · d_i(p)))` over the local sizes `h_i` with `d_i` the distance to +/// their regions, so it changes by at most `g` per unit length. Exact (no background grid): +/// meant for the structured meshers, which query it along curves; the tetrahedral refinement +/// bakes the same local sizes into its graded grid ([`crate::tet::geometry::GradedSize`]). +#[derive(Debug, Clone, PartialEq)] +pub struct RegionSize { + global: f64, + growth: f64, + /// Each local size with its region's bounds grown by the reach of its refinement. + locals: Vec<(LocalSize, [f64; 3], [f64; 3])>, +} + +impl RegionSize { + /// A checked field: a finite, positive `global` and local sizes, `growth` finite and + /// positive, regions with finite coordinates and at most [`MAX_REGION_TRIANGLES`] triangles + /// in all. Local sizes not below `global` change nothing and are dropped. + pub fn new(global: f64, growth: f64, locals: Vec) -> Result { + let positive = |what: &str, v: f64| { + if v > 0.0 && v.is_finite() { Ok(v) } else { Err(Error::Invalid(format!("{what} must be finite and positive, got {v}"))) } + }; + positive("the global size", global)?; + positive("the growth rate", growth)?; + let triangles: usize = locals.iter().map(|l| if let Region::Triangles(t) = &l.region { t.len() } else { 0 }).sum(); + if triangles > MAX_REGION_TRIANGLES { + return Err(Error::Limit(format!("size regions hold {triangles} triangles; the limit is {MAX_REGION_TRIANGLES}"))); + } + let mut out = Vec::with_capacity(locals.len()); + for l in locals { + positive("a local size", l.size)?; + match &l.region { + Region::Box { lo, hi } => fieldcraft_model::check_coordinates(&[lo.as_slice(), hi.as_slice()].concat())?, + Region::Sphere { center, radius } => { + fieldcraft_model::check_coordinates(center)?; + positive("a sphere's radius", *radius)?; + } + Region::Triangles(t) => fieldcraft_model::check_coordinates(t.as_flattened().as_flattened())?, + } + if l.size >= global { + continue; + } + let (lo, hi) = l.region.bounds(); + if (0..3).any(|k| !(lo[k] <= hi[k])) { + continue; // no triangles + } + let reach = (global - l.size) / growth; + out.push((l, lo.map(|c| c - reach), hi.map(|c| c + reach))); + } + Ok(RegionSize { global, growth, locals: out }) + } + + pub fn global(&self) -> f64 { + self.global + } + + pub fn growth(&self) -> f64 { + self.growth + } + + /// Whether any local size refines the global one. + pub fn is_uniform(&self) -> bool { + self.locals.is_empty() + } + + /// The local sizes that refine the global one. + pub fn locals(&self) -> impl Iterator { + self.locals.iter().map(|l| &l.0) + } + + /// The smallest size anywhere. + pub fn finest(&self) -> f64 { + self.locals().map(|l| l.size).fold(self.global, f64::min) + } +} + +impl SizeField for RegionSize { + fn size(&self, p: [f64; 3]) -> f64 { + let mut h = self.global; + for (l, lo, hi) in &self.locals { + if (0..3).any(|k| p[k] < lo[k] || p[k] > hi[k]) { + continue; + } + h = h.min(l.size + self.growth * l.region.distance(p)); + } + h + } +} + #[cfg(test)] mod tests { use super::*; @@ -260,4 +404,45 @@ mod tests { assert!(loop_divisions(&[f64::INFINITY]).is_err()); assert!(even_divisions(f64::NAN).is_err()); } + + #[test] + fn region_size_refines_inside_and_grows_away() { + let boxed = LocalSize { region: Region::Box { lo: [0.0; 3], hi: [2.0; 3] }, size: 0.5 }; + let ball = LocalSize { region: Region::Sphere { center: [10.0, 0.0, 0.0], radius: 1.0 }, size: 0.25 }; + let face = LocalSize { region: Region::Triangles(vec![[[0.0, 10.0, 0.0], [1.0, 10.0, 0.0], [0.0, 11.0, 0.0]]]), size: 1.0 }; + let coarse = LocalSize { region: Region::Box { lo: [0.0; 3], hi: [1.0; 3] }, size: 5.0 }; + let f = RegionSize::new(2.0, 0.5, vec![boxed, ball, face, coarse]).unwrap(); + assert_eq!(f.locals().count(), 3, "the coarser control is dropped"); + assert_eq!(f.size([1.0, 1.0, 1.0]), 0.5); + assert_eq!(f.size([3.0, 1.0, 1.0]), 1.0, "one unit out: 0.5 + 0.5"); + assert_eq!(f.size([20.0, 20.0, 20.0]), 2.0); + assert_eq!(f.size([10.5, 0.0, 0.0]), 0.25); + assert_eq!(f.size([0.2, 10.2, 0.0]), 1.0); + assert!((f.size([0.2, 10.2, 1.0]) - 1.5).abs() < 1e-12); + assert_eq!(f.finest(), 0.25); + // Graded: the size never changes faster than the growth rate. + let pts: Vec<[f64; 3]> = (0..400).map(|i| [-5.0 + 0.05 * i as f64, 1.0, 1.0]).collect(); + for w in pts.windows(2) { + assert!((f.size(w[1]) - f.size(w[0])).abs() <= 0.5 * 0.05 + 1e-12); + } + // Divisions follow it: a line through the box is cut finer inside. + let line = [[-4.0, 1.0, 1.0], [6.0, 1.0, 1.0]]; + let n = edge_divisions(&f, &line).unwrap(); + let pts = divide_polyline(&f, &line, n).unwrap(); + let inside = pts.windows(2).filter(|w| w[0][0] >= 0.0 && w[1][0] <= 2.0).map(|w| w[1][0] - w[0][0]).fold(0.0, f64::max); + assert!(inside < 0.6, "{inside}"); + } + + #[test] + fn hostile_region_sizes_are_refused() { + let ok = || LocalSize { region: Region::Sphere { center: [0.0; 3], radius: 1.0 }, size: 0.5 }; + assert!(RegionSize::new(0.0, 0.3, vec![]).is_err()); + assert!(RegionSize::new(1.0, f64::NAN, vec![]).is_err()); + assert!(RegionSize::new(1.0, 0.3, vec![LocalSize { size: -1.0, ..ok() }]).is_err()); + assert!(RegionSize::new(1.0, 0.3, vec![LocalSize { region: Region::Sphere { center: [f64::NAN; 3], radius: 1.0 }, size: 0.5 }]).is_err()); + assert!(RegionSize::new(1.0, 0.3, vec![LocalSize { region: Region::Triangles(vec![[[f64::INFINITY; 3]; 3]]), size: 0.5 }]).is_err()); + let empty = RegionSize::new(1.0, 0.3, vec![LocalSize { region: Region::Triangles(vec![]), size: 0.5 }]).unwrap(); + assert!(empty.is_uniform()); + assert!(RegionSize::new(1.0, 0.3, vec![ok()]).unwrap().size([f64::NAN; 3]) > 0.0); + } } diff --git a/crates/mesh/src/tet/geometry.rs b/crates/mesh/src/tet/geometry.rs index 122a2bb..34df939 100644 --- a/crates/mesh/src/tet/geometry.rs +++ b/crates/mesh/src/tet/geometry.rs @@ -13,7 +13,7 @@ //! analytic (von Neumann's alternating projections converge to the intersection for surfaces //! meeting transversally); the result is kept only when it stays close to the polyline. -use crate::size::SizeField; +use crate::size::{LocalSize, Region, SizeField}; use crate::surface::BoundarySurface; use super::measure::{self, P, add, cross, dot, norm, scale, sub}; @@ -842,9 +842,13 @@ pub const MAX_GRID_NODES: usize = 1 << 21; impl GradedSize { /// `samples`: points with their local size. `grading`: the largest change of size per unit /// length. - pub fn new(global: f64, samples: &[(P, f64)], lo: P, hi: P, grading: f64) -> GradedSize { + /// + /// `locals`: regions with their size ([`crate::size::LocalSize`]), seeded at every grid node + /// within a cell of the region (boxes and balls: inside too) at the local size plus the + /// grading times the distance, then graded with the rest. + pub fn new(global: f64, samples: &[(P, f64)], locals: &[LocalSize], lo: P, hi: P, grading: f64) -> GradedSize { let diag = measure::dist(lo, hi).max(f64::MIN_POSITIVE); - let finest = samples.iter().map(|s| s.1).fold(global, f64::min).max(diag * 1e-4); + let finest = samples.iter().map(|s| s.1).chain(locals.iter().map(|l| l.size)).filter(|h| *h > 0.0).fold(global, f64::min).max(diag * 1e-4); let mut cell = finest.max(diag / 64.0) * 0.5; let dims_for = |c: f64| [0, 1, 2].map(|k| (((hi[k] - lo[k]) / c).ceil().max(1.0) as usize) + 2); let mut dims = dims_for(cell); @@ -879,6 +883,47 @@ impl GradedSize { } } } + // Seed: the grid nodes near each local region. + let reach = cell * 3f64.sqrt(); + let range = |a: f64, b: f64, k: usize| -> (usize, usize) { + let first = ((a - reach - lo[k]) / cell).ceil().max(0.0); + let last = ((b + reach - lo[k]) / cell).floor().min((dims[k] - 1) as f64); + if first.is_finite() && last.is_finite() && first <= last { (first as usize, last as usize) } else { (1, 0) } + }; + let mut seed_near = |a: P, b: P, h: f64, distance: &dyn Fn(P) -> f64| { + let r = [0, 1, 2].map(|k| range(a[k], b[k], k)); + for k in r[2].0..=r[2].1 { + for j in r[1].0..=r[1].1 { + for i in r[0].0..=r[0].1 { + let node = [lo[0] + i as f64 * cell, lo[1] + j as f64 * cell, lo[2] + k as f64 * cell]; + let d = distance(node); + if d <= reach + && let Some(s) = values.get_mut(idx(i, j, k)) + { + *s = s.min(h + g * d); + } + } + } + } + }; + for l in locals { + if !(l.size > 0.0 && l.size.is_finite() && l.size < global) { + continue; + } + match &l.region { + Region::Triangles(tris) => { + for t in tris { + let a = [0, 1, 2].map(|k| t[0][k].min(t[1][k]).min(t[2][k])); + let b = [0, 1, 2].map(|k| t[0][k].max(t[1][k]).max(t[2][k])); + seed_near(a, b, l.size, &|p| measure::dist(p, closest_on_triangle(p, t[0], t[1], t[2]))); + } + } + region => { + let (a, b) = region.bounds(); + seed_near(a, b, l.size, &|p| region.distance(p)); + } + } + } // Fast sweeping: eight sweep orders, each relaxing against the three axis neighbours // behind the sweep direction (and diagonals through them over passes). if g > 0.0 { @@ -983,9 +1028,24 @@ mod tests { assert_eq!(closest_on_segment([0.5, 2.0, 0.0], a, b), [0.5, 0.0, 0.0]); } + #[test] + fn graded_size_takes_local_regions() { + let boxed = LocalSize { region: Region::Box { lo: [8.0; 3], hi: [12.0; 3] }, size: 2.0 }; + let face = LocalSize { region: Region::Triangles(vec![[[0.0, 0.0, 20.0], [20.0, 0.0, 20.0], [0.0, 20.0, 20.0]]]), size: 1.0 }; + let f = GradedSize::new(8.0, &[], &[boxed, face], [0.0; 3], [20.0, 20.0, 20.0], 0.5); + let inside = f.size([10.0, 10.0, 10.0]); + assert!((inside - 2.0).abs() < 0.3, "{inside}"); + let on_face = f.size([5.0, 5.0, 20.0]); + assert!((on_face - 1.0).abs() < 0.3, "{on_face}"); + // Four units off the box: about 2 + 0.5 × 4. + let off = f.size([10.0, 10.0, 4.0]); + assert!((off - 4.0).abs() < 0.6, "{off}"); + assert_eq!(f.size([20.0, 20.0, 0.0]), 8.0); + } + #[test] fn graded_size_grows_away_from_samples() { - let f = GradedSize::new(10.0, &[([0.0; 3], 1.0)], [0.0; 3], [20.0, 20.0, 20.0], 0.5); + let f = GradedSize::new(10.0, &[([0.0; 3], 1.0)], &[], [0.0; 3], [20.0, 20.0, 20.0], 0.5); let near = f.size([0.0; 3]); let mid = f.size([4.0, 0.0, 0.0]); let far = f.size([20.0, 20.0, 20.0]); diff --git a/crates/mesh/src/tet/mod.rs b/crates/mesh/src/tet/mod.rs index 72e1a8f..6dacbe7 100644 --- a/crates/mesh/src/tet/mod.rs +++ b/crates/mesh/src/tet/mod.rs @@ -97,6 +97,9 @@ pub struct QualityOptions { /// 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, + /// Local sizes (size controls): the size field is at most each one's size in its region and + /// grows from it at the grading rate. + pub locals: Vec, /// B-rep face tags whose boundary triangles are kept exactly: their vertices stay where they /// are and no triangle of theirs is split, flipped or removed. The faces a pyramid /// transition shares with the tetrahedra ([`crate::transition`]) are frozen, so the two @@ -122,6 +125,7 @@ impl Default for QualityOptions { max_tets: 5_000_000, max_points: 1_000_000, seeds: Vec::new(), + locals: Vec::new(), frozen_faces: Vec::new(), } } @@ -237,7 +241,8 @@ fn bounds(points: &[[f64; 3]]) -> ([f64; 3], [f64; 3]) { (lo, hi) } -/// The size field of a boundary: `global` limited by curvature, graded. +/// The size field of a boundary: `global` limited by curvature, thin walls, edge seeds and local +/// sizes, graded. pub fn size_field(surface: &BoundarySurface, geometry: &Tessellation, global: f64, o: &QualityOptions) -> GradedSize { let (lo, hi) = bounds(surface.points()); let floor = o.min_size * global; @@ -268,7 +273,7 @@ pub fn size_field(surface: &BoundarySurface, geometry: &Tessellation, global: f6 } } } - GradedSize::new(global, &samples, lo, hi, o.grading) + GradedSize::new(global, &samples, &o.locals, lo, hi, o.grading) } /// Put the seeds' nodes on their feature curves, pinned (corner vertices, which neither