mesh: local sizes in boxes, balls and on surfaces, graded (#435)

RegionSize is the exact graded field the structured meshers query;
the tetrahedral refinement bakes the same local sizes into its graded
background grid.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo

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<LocalSize>) -> Result<RegionSize> {
+        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<Item = &LocalSize> {
+        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);
+    }
 }