From 600b2ff633dad617844dee4a67e0b3c3ff97c72c Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Thu, 15 Oct 2026 16:02:11 -0400 Subject: [PATCH] mesh: hex extrusion of a quad section up to a target surface (#400) Each source node follows its own path (a direction or the section normal) to the target, found by Newton on the target's tangent plane; layers are graded by a bias; side faces, when given, hold the boundary nodes of every layer (alternating projections at the target), and interior nodes take the boundary-residual correction. Source nodes are kept, so a section taken from a neighbouring mesh conforms with it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/mesh/src/hex/extrude_to.rs | 580 ++++++++++++++++++++++++ crates/mesh/src/hex/extrude_to/tests.rs | 178 ++++++++ crates/mesh/src/hex/mod.rs | 2 + crates/mesh/src/hex/sweep.rs | 2 +- 4 files changed, 761 insertions(+), 1 deletion(-) create mode 100644 crates/mesh/src/hex/extrude_to.rs create mode 100644 crates/mesh/src/hex/extrude_to/tests.rs diff --git a/crates/mesh/src/hex/extrude_to.rs b/crates/mesh/src/hex/extrude_to.rs new file mode 100644 index 0000000..9b062e7 --- /dev/null +++ b/crates/mesh/src/hex/extrude_to.rs @@ -0,0 +1,580 @@ +//! Extruding a quadrilateral section up to a target surface. +//! +//! [`extrude_to_surface`] stacks layers of hexahedra on a quad mesh of a section (a face of a +//! body, a partition face or a 2D mesh) and ends them on one or more B-rep faces, the *target*, +//! instead of on a copy of the section: a boss meeting a curved wall, an air column up to a +//! domed lid, a fillet region capped by a curved face. The region need not be sweepable +//! between two parallel faces. +//! +//! **Paths.** Every source node has its own path: a straight line along the given direction, or +//! along the section's normal at the node. The path's end is where the line meets the target, +//! found by Newton's method on the tangent plane of the target's closest point (each step +//! intersects the line with the plane through the closest point, normal to the gap; it is exact +//! after one step on a plane and converges quadratically on smooth faces). The layers divide each +//! path by the same parameters, uniform or graded geometrically by `bias` (the last layer's +//! thickness over the first's). +//! +//! **Side faces.** When the section's sides are given, the boundary nodes of every layer are +//! projected onto them (onto the curve between two sides at a corner), and the boundary's last +//! layer is found by alternating projections onto the sides and the target (von Neumann's +//! alternating projections converge to the intersection of surfaces meeting transversally), so +//! the layers follow side faces that are not straight along the path. Interior nodes follow +//! their straight paths plus the boundary's displacement from its straight paths interpolated by +//! inverse-distance weights (the boundary-residual correction of Knupp, *Next-generation sweep +//! tool: a method for generating all-hex meshes on two-and-one-half dimensional geometries*, +//! Proc. 7th International Meshing Roundtable, 1998), and their last layer is placed back on the +//! target along the path direction. +//! +//! The source nodes are kept as given (the first nodes of the result, in order), so a section +//! taken from a neighbouring mesh conforms with it. Faces are projected onto through +//! [`crate::project::Geometry`], the trimmed B-rep projection of projection blocking: a point +//! beyond a face's boundary goes to the face, never to its surface's extension, so a path that +//! misses a target face does not converge on it. The result passes the certified validity +//! gate (every element's Jacobian determinant proven positive) and the tiling check (the +//! elements' volume equals the volume their boundary encloses). + +use std::collections::BTreeMap; + +use fieldcraft_model::ElementKind; + +use super::face_flux; +use super::geom::{add, dist, scale, unit}; +use super::sweep::shepard_weights; +use crate::project::Geometry; +use crate::quality::{GateOptions, GateReport, Metric, element_quality, validity_gate_checks}; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +/// The most layers. +pub const MAX_LAYERS: usize = 10_000; +/// Newton steps per path end, and alternating-projection rounds per node. +const MAX_NEWTON: usize = 60; +const MAX_ALTERNATING: usize = 100; +/// The most interior-node × boundary-node × layer products the residual correction may spend. +const MAX_CORRECTION_WORK: u64 = 2_000_000_000; +/// The largest layer bias (last over first thickness) and its inverse. +pub const MAX_BIAS: f64 = 1000.0; + +/// The path of each source node. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum ExtrudePath { + /// Every node along this direction (need not be unit length). + Direction(P), + /// Every node along the section's normal there, turned toward the target. + Normal, +} + +/// How many layers. +#[derive(Debug, Clone, Copy, PartialEq)] +pub enum ExtrudeLayers { + Count(usize), + /// Layers of about this thickness (along the mean path length). + Size(f64), +} + +/// What [`extrude_to_surface`] extrudes and where to. +#[derive(Debug, Clone)] +pub struct ExtrudeSpec<'a> { + /// The section's nodes and quadrilaterals (indices into `nodes`). + pub nodes: &'a [P], + pub quads: &'a [[u32; 4]], + pub path: ExtrudePath, + /// The target faces and the side faces (face ids of the geometry). + pub targets: &'a [u32], + pub sides: &'a [u32], + pub layers: ExtrudeLayers, + /// The last layer's thickness over the first's (1: uniform). + pub bias: f64, + /// How far a node may be from a face and still be on it (the geometry's accuracy). + pub tolerance: f64, + /// The most hexahedra ([`Error::Limit`] beyond). + pub max_elements: usize, +} + +/// What a boundary face of the extruded mesh lies on. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum ExtrudeRole { + /// The section. + Source, + /// A target face (its id). + Target(u32), + /// The side face it lies on, when sides were given. + Side(Option), +} + +/// A boundary face of the extruded mesh. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub struct ExtrudeBoundary { + /// Element index (0-based). + pub element: u32, + /// Face number (1-based, as in [`ElementKind::faces`]). + pub face: u8, + pub role: ExtrudeRole, +} + +/// The extruded hex8 mesh. +#[derive(Debug, Clone, PartialEq)] +pub struct Extruded { + /// The source nodes first (unchanged, in order), then each layer's nodes. + pub nodes: Vec

, + pub hexes: Vec<[u32; 8]>, + pub boundary: Vec, + pub layers: usize, + pub min_scaled_jacobian: f64, + pub mean_scaled_jacobian: f64, + /// The largest distance of a last-layer node from the target. + pub max_target_gap: f64, + pub gate: GateReport, +} + +/// Extrude the section `spec` up to its target surface. See the module docs. +pub fn extrude_to_surface(spec: &ExtrudeSpec, geometry: &dyn Geometry) -> Result { + let n_src = spec.nodes.len(); + check(spec)?; + let tol = spec.tolerance; + let ex = Extrusion { g: geometry, targets: spec.targets, tol }; + let boundary = section_boundary(spec.quads)?; + + // Path directions. + let dirs: Vec

= match spec.path { + ExtrudePath::Direction(d) => { + let d = unit(d).ok_or_else(|| Error::Invalid("the direction has no length".into()))?; + vec![d; n_src] + } + ExtrudePath::Normal => ex.normals(spec)?, + }; + + // Each side node's side faces. + let mut side_faces: BTreeMap> = BTreeMap::new(); + if !spec.sides.is_empty() { + for &v in &boundary.nodes { + let p = at(spec.nodes, v)?; + let on: Vec = spec.sides.iter().copied().filter(|&f| geometry.project(f, p).is_some_and(|q| dist(p, q) <= 10.0 * tol)).collect(); + if on.is_empty() { + return Err(Error::Invalid(format!("the section's boundary node at {} is on none of the side faces", fmt(p)))); + } + side_faces.insert(v, on); + } + } + + // Straight ends, and the boundary's ends on the sides. + let mut straight = Vec::with_capacity(n_src); + for (i, (&s, &d)) in spec.nodes.iter().zip(&dirs).enumerate() { + let (e, _) = ex.hit(s, d, true).ok_or_else(|| { + Error::Failed(format!("the path of section node {i} from {} along {} does not reach the target surface", fmt(s), fmt(d))) + })?; + straight.push(e); + } + let mut ends: BTreeMap = BTreeMap::new(); + for (&v, faces) in &side_faces { + let e = at(&straight, v)?; + ends.insert( + v, + ex.corner(faces, e).ok_or_else(|| Error::Failed(format!("the side faces and the target surface do not meet near {}", fmt(e))))?, + ); + } + + // Layer parameters. + let mean_length = spec.nodes.iter().zip(&straight).map(|(s, e)| dist(*s, *e)).sum::() / n_src as f64; + let n = match spec.layers { + ExtrudeLayers::Count(n) => n, + ExtrudeLayers::Size(h) => ((mean_length / h).round() as usize).max(1), + }; + if n == 0 || n > MAX_LAYERS { + return Err(Error::Invalid(format!("{n} layers; an extrusion has 1 to {MAX_LAYERS}"))); + } + let count = spec.quads.len().saturating_mul(n); + if count > spec.max_elements { + return Err(Error::Limit(format!("the extrusion would have {count} hexahedra; the limit is {}", spec.max_elements))); + } + let ts = parameters(n, spec.bias); + + // Interior correction weights, in the plane across the mean direction. + let correct = !side_faces.is_empty(); + let ring: Vec = side_faces.keys().copied().collect(); + let interior: Vec = (0..n_src as u32).filter(|v| !side_faces.contains_key(v)).collect(); + let work = (interior.len() as u64).saturating_mul(ring.len() as u64).saturating_mul(n as u64); + if correct && work > MAX_CORRECTION_WORK { + return Err(Error::Limit(format!("placing the interior nodes would take {work} steps; the limit is {MAX_CORRECTION_WORK}"))); + } + let weights: Vec> = if correct { + let mean_dir = unit(dirs.iter().fold([0.0; 3], |a, d| add(a, *d))).unwrap_or([0.0, 0.0, 1.0]); + let (u, w) = across(mean_dir); + let flat = |p: P| [dot(p, u), dot(p, w)]; + let ring2: Vec<[f64; 2]> = ring.iter().map(|&v| at(spec.nodes, v).map(flat)).collect::>()?; + interior.iter().map(|&v| at(spec.nodes, v).map(|p| shepard_weights(flat(p), &ring2))).collect::>()? + } else { + Vec::new() + }; + + let mut nodes: Vec

= Vec::with_capacity(n_src * (n + 1)); + nodes.extend_from_slice(spec.nodes); + let lerp = |s: P, e: P, t: f64| add(s, scale(sub(e, s), t)); + for (k, &t) in ts.iter().enumerate().skip(1) { + let last = k == n; + let mut layer: Vec

= spec.nodes.iter().zip(&straight).map(|(s, e)| lerp(*s, *e, t)).collect(); + // Boundary: onto the sides; their displacement from the straight paths. + let mut residual = Vec::with_capacity(ring.len()); + for &v in &ring { + let i = v as usize; + let straight_at = layer.get(i).copied().unwrap_or([0.0; 3]); + let x = if last { + ends.get(&v).copied().unwrap_or(straight_at) + } else { + let faces = side_faces.get(&v).map_or(&[][..], Vec::as_slice); + let s = at(spec.nodes, v)?; + let target = lerp(s, ends.get(&v).copied().unwrap_or(straight_at), t); + ex.onto(faces, target).ok_or_else(|| Error::Failed(format!("cannot place a layer node on the side faces near {}", fmt(target))))? + }; + residual.push(sub(x, straight_at)); + if let Some(slot) = layer.get_mut(i) { + *slot = x; + } + } + // Interior: the straight path plus the interpolated displacement; the last layer back on + // the target along the path. + for (j, &v) in interior.iter().enumerate() { + let i = v as usize; + let Some(&y0) = layer.get(i) else { continue }; + let correction = weights + .get(j) + .map_or([0.0; 3], |w| w.iter().fold([0.0; 3], |acc, &(r, wr)| add(acc, scale(residual.get(r).copied().unwrap_or([0.0; 3]), wr)))); + let mut y = add(y0, correction); + if last && correct { + let d = dirs.get(i).copied().unwrap_or([0.0, 0.0, 1.0]); + y = ex + .hit(y, d, false) + .map(|(p, _)| p) + .ok_or_else(|| Error::Failed(format!("the path near {} does not reach the target surface", fmt(y))))?; + } + if let Some(slot) = layer.get_mut(i) { + *slot = y; + } + } + nodes.extend(layer); + } + fieldcraft_model::check_coordinates(nodes.as_flattened())?; + + // Hexahedra, each quad turned so its normal runs along the path. + let mut hexes = Vec::with_capacity(count); + let mut quads = Vec::with_capacity(spec.quads.len()); + for q in spec.quads { + let [a, b, c, d] = q.map(|v| at(spec.nodes, v).unwrap_or([0.0; 3])); + let normal = cross(sub(c, a), sub(d, b)); + let along = q.iter().fold([0.0; 3], |acc, &v| { + let s = at(spec.nodes, v).unwrap_or([0.0; 3]); + add(acc, sub(nodes.get(n_src + v as usize).copied().unwrap_or(s), s)) + }); + quads.push(if dot(normal, along) < 0.0 { [q[0], q[3], q[2], q[1]] } else { *q }); + } + let stride = n_src as u32; + for q in &quads { + for k in 0..n as u32 { + let lo = q.map(|v| v + k * stride); + let hi = q.map(|v| v + (k + 1) * stride); + hexes.push([lo[0], lo[1], lo[2], lo[3], hi[0], hi[1], hi[2], hi[3]]); + } + } + + // Boundary faces: the section, the last layer and the sides. + let centre = |e: usize, face: usize| -> P { + let h = hexes.get(e).copied().unwrap_or([0; 8]); + let local = ElementKind::Hex8.faces().get(face).copied().unwrap_or(&[]); + let pts: Vec

= local.iter().filter_map(|&i| h.get(i)).filter_map(|&v| nodes.get(v as usize).copied()).collect(); + scale(pts.iter().fold([0.0; 3], |a, p| add(a, *p)), 1.0 / pts.len().max(1) as f64) + }; + let nearest = |faces: &[u32], p: P| -> Option { + faces + .iter() + .filter_map(|&f| geometry.project(f, p).map(|q| (f, dist(p, q)))) + .min_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0))) + .map(|(f, _)| f) + }; + let mut bfaces = Vec::new(); + for (qi, q) in quads.iter().enumerate() { + let first = qi * n; + let last = first + n - 1; + bfaces.push(ExtrudeBoundary { element: first as u32, face: 1, role: ExtrudeRole::Source }); + let f = nearest(spec.targets, centre(last, 1)).ok_or_else(|| Error::Failed("a last-layer face is on no target face".into()))?; + bfaces.push(ExtrudeBoundary { element: last as u32, face: 2, role: ExtrudeRole::Target(f) }); + for m in 0..4 { + let (a, b) = (q[m], q[(m + 1) % 4]); + if !boundary.edges.contains(&(a.min(b), a.max(b))) { + continue; + } + for e in first..=last { + let side = if spec.sides.is_empty() { None } else { nearest(spec.sides, centre(e, 2 + m)) }; + bfaces.push(ExtrudeBoundary { element: e as u32, face: 3 + m as u8, role: ExtrudeRole::Side(side) }); + } + } + } + + // Quality and the certified gate. + let kind = ElementKind::Hex8; + let coords = |e: usize| -> Result> { hexes.get(e).map_or(Ok(Vec::new()), |h| h.iter().map(|&v| at(&nodes, v)).collect()) }; + let mut enclosed = 0.0; + for b in &bfaces { + let x = coords(b.element as usize)?; + let local = kind.faces().get(b.face as usize - 1).copied().unwrap_or(&[]); + let face: Vec

= local.iter().filter_map(|&i| x.get(i).copied()).collect(); + enclosed += face_flux(&face); + } + let (mut min_sj, mut sum_sj, mut worst) = (f64::INFINITY, 0.0, 0usize); + let mut checks = Vec::with_capacity(hexes.len()); + for e in 0..hexes.len() { + let q = element_quality(kind, &coords(e)?)?; + let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0); + if sj < min_sj { + (min_sj, worst) = (sj, e); + } + sum_sj += sj; + checks.push((e as u32 + 1, q.jacobian.ok_or_else(|| Error::Failed("an element was not certified".into()))?)); + } + let gate = validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })? + .into_result() + .map_err(|e| { + Error::Failed(format!( + "{e} (worst element {} near {}: a path crosses another or a side; try other layers, a smaller bias or the side faces)", + worst + 1, + fmt(centre(worst, 0)) + )) + })?; + let max_target_gap = nodes + .iter() + .skip(n * n_src) + .map(|&p| spec.targets.iter().filter_map(|&f| geometry.project(f, p)).map(|q| dist(p, q)).fold(f64::INFINITY, f64::min)) + .fold(0.0, f64::max); + Ok(Extruded { + nodes, + hexes, + boundary: bfaces, + layers: n, + min_scaled_jacobian: min_sj, + mean_scaled_jacobian: sum_sj / count.max(1) as f64, + max_target_gap, + gate, + }) +} + +fn check(spec: &ExtrudeSpec) -> Result<()> { + if spec.nodes.is_empty() || spec.quads.is_empty() { + return Err(Error::Invalid("the section has no quadrilaterals".into())); + } + fieldcraft_model::check_coordinates(spec.nodes.as_flattened())?; + if let Some(v) = spec.quads.iter().flatten().find(|&&v| v as usize >= spec.nodes.len()) { + return Err(Error::Invalid(format!("section node {v} is out of range"))); + } + if spec.quads.iter().any(|q| (0..4).any(|i| (i + 1..4).any(|j| q[i] == q[j]))) { + return Err(Error::Invalid("a section quadrilateral repeats a node".into())); + } + if spec.targets.is_empty() { + return Err(Error::Invalid("no target surface".into())); + } + if !(spec.tolerance > 0.0 && spec.tolerance.is_finite()) { + return Err(Error::Invalid("the tolerance must be a positive number".into())); + } + if !(spec.bias >= 1.0 / MAX_BIAS && spec.bias <= MAX_BIAS) { + return Err(Error::Invalid(format!("the bias must be from {} to {MAX_BIAS}", 1.0 / MAX_BIAS))); + } + match spec.layers { + ExtrudeLayers::Count(n) if n == 0 || n > MAX_LAYERS => Err(Error::Invalid(format!("an extrusion has 1 to {MAX_LAYERS} layers"))), + ExtrudeLayers::Size(h) if !(h > 0.0 && h.is_finite()) => Err(Error::Invalid("the layer size must be a positive number".into())), + _ => Ok(()), + } +} + +/// The section's boundary: edges on one quadrilateral, and their nodes. +struct SectionBoundary { + edges: Vec<(u32, u32)>, + nodes: Vec, +} + +fn section_boundary(quads: &[[u32; 4]]) -> Result { + let mut uses: BTreeMap<(u32, u32), u32> = BTreeMap::new(); + for q in quads { + for m in 0..4 { + let (a, b) = (q[m], q[(m + 1) % 4]); + *uses.entry((a.min(b), a.max(b))).or_insert(0) += 1; + } + } + if let Some((&(a, b), _)) = uses.iter().find(|(_, n)| **n > 2) { + return Err(Error::Invalid(format!("the section is not a surface: the edge between nodes {a} and {b} is on more than two quadrilaterals"))); + } + let edges: Vec<(u32, u32)> = uses.into_iter().filter(|(_, n)| *n == 1).map(|(e, _)| e).collect(); + let mut nodes: Vec = edges.iter().flat_map(|&(a, b)| [a, b]).collect(); + nodes.sort_unstable(); + nodes.dedup(); + Ok(SectionBoundary { edges, nodes }) +} + +/// The layer parameters `0 = t_0 < … < t_n = 1`, thicknesses in geometric progression with the +/// last over the first equal to `bias`. +pub(crate) fn parameters(n: usize, bias: f64) -> Vec { + let r = if n > 1 { bias.powf(1.0 / (n - 1) as f64) } else { 1.0 }; + let mut out = Vec::with_capacity(n + 1); + let (mut sum, mut h) = (0.0, 1.0); + out.push(0.0); + for _ in 0..n { + sum += h; + out.push(sum); + h *= r; + } + let total = sum; + for t in &mut out { + *t /= total; + } + if let Some(t) = out.last_mut() { + *t = 1.0; + } + out +} + +/// The geometry queries of one extrusion. +struct Extrusion<'a> { + g: &'a dyn Geometry, + targets: &'a [u32], + tol: f64, +} + +impl Extrusion<'_> { + /// Where the line `o + t d` meets the target: the hit with the smallest `t` ahead of `o` + /// (`ahead`), or nearest to `o` either way. Returns the point and its target face. + fn hit(&self, o: P, d: P, ahead: bool) -> Option<(P, u32)> { + let mut best: Option<(P, u32, f64)> = None; + for &f in self.targets { + let Some((p, t)) = self.hit_face(f, o, d) else { continue }; + if ahead && t <= self.tol { + continue; + } + let key = if ahead { t } else { t.abs() }; + if best.is_none_or(|b| key < b.2) { + best = Some((p, f, key)); + } + } + best.map(|(p, f, _)| (p, f)) + } + + /// Newton's method for the line `o + t d` on face `f`: each step meets the tangent plane at + /// the closest point. `None` when it does not converge onto the face (the line misses it, + /// or grazes it). + fn hit_face(&self, f: u32, o: P, d: P) -> Option<(P, f64)> { + let q0 = self.g.project(f, o)?; + let mut t = dot(sub(q0, o), d); + let eps = 1e-10 * (1.0 + norm(o).max(dist(o, q0))); + for _ in 0..MAX_NEWTON { + let x = add(o, scale(d, t)); + let q = self.g.project(f, x)?; + let gap = sub(x, q); + let g2 = dot(gap, gap); + if g2.sqrt() <= eps { + return Some((q, t)); + } + let den = dot(d, gap); + if den.abs() <= 1e-6 * g2.sqrt() { + return None; + } + t -= g2 / den; + if !t.is_finite() { + return None; + } + } + // A free-form face's facets can make the last steps cycle within its accuracy. + let x = add(o, scale(d, t)); + let q = self.g.project(f, x)?; + (dist(x, q) <= self.tol).then_some((q, t)) + } + + /// `p` on all of `faces` (one face: its closest point; several: alternating projections). + fn onto(&self, faces: &[u32], p: P) -> Option

{ + let mut x = p; + if let [f] = faces { + return self.g.project(*f, x); + } + for _ in 0..MAX_ALTERNATING { + let before = x; + for &f in faces { + x = self.g.project(f, x)?; + } + if dist(before, x) <= 1e-12 * (1.0 + norm(x)) { + break; + } + } + faces.iter().all(|&f| self.g.project(f, x).is_some_and(|q| dist(q, x) <= self.tol)).then_some(x) + } + + /// The point near `p` on the side `faces` and the target (alternating projections, ending on + /// the target exactly). + fn corner(&self, faces: &[u32], p: P) -> Option

{ + let mut x = p; + let mut face = self.nearest_target(x)?; + for _ in 0..MAX_ALTERNATING { + let before = x; + for &f in faces { + x = self.g.project(f, x)?; + } + face = self.nearest_target(x)?; + x = self.g.project(face, x)?; + if dist(before, x) <= 1e-12 * (1.0 + norm(x)) { + break; + } + } + let on_sides = faces.iter().all(|&f| self.g.project(f, x).is_some_and(|q| dist(q, x) <= self.tol)); + let on_target = self.g.project(face, x).is_some_and(|q| dist(q, x) <= self.tol); + (on_sides && on_target).then_some(x) + } + + fn nearest_target(&self, p: P) -> Option { + self.targets + .iter() + .filter_map(|&f| self.g.project(f, p).map(|q| (f, dist(p, q)))) + .min_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0))) + .map(|(f, _)| f) + } + + /// Unit normals of the section at its nodes, turned toward the target. + fn normals(&self, spec: &ExtrudeSpec) -> Result> { + let mut sum = vec![[0.0; 3]; spec.nodes.len()]; + for q in spec.quads { + let [a, b, c, d] = q.map(|v| at(spec.nodes, v).unwrap_or([0.0; 3])); + let n = cross(sub(c, a), sub(d, b)); + for &v in q { + if let Some(s) = sum.get_mut(v as usize) { + *s = add(*s, n); + } + } + } + let mut out = Vec::with_capacity(sum.len()); + for (i, s) in sum.iter().enumerate() { + out.push(unit(*s).ok_or_else(|| Error::Invalid(format!("the section has no normal at node {i} (it folds there)")))?); + } + // Toward the target: the side of the section its nearest point is on. + let centre = scale(spec.nodes.iter().fold([0.0; 3], |a, p| add(a, *p)), 1.0 / spec.nodes.len() as f64); + let mean = out.iter().fold([0.0; 3], |a, n| add(a, *n)); + let f = self.nearest_target(centre).ok_or_else(|| Error::Failed("no target face to extrude toward".into()))?; + let q = self.g.project(f, centre).ok_or_else(|| Error::Failed("no target face to extrude toward".into()))?; + if dot(sub(q, centre), mean) < 0.0 { + for n in &mut out { + *n = scale(*n, -1.0); + } + } + Ok(out) + } +} + +/// Two unit vectors across `d`. +fn across(d: P) -> (P, P) { + let t = if d[0].abs() < 0.9 { [1.0, 0.0, 0.0] } else { [0.0, 1.0, 0.0] }; + let u = unit(cross(d, t)).unwrap_or([1.0, 0.0, 0.0]); + (u, cross(d, u)) +} + +fn at(nodes: &[P], v: u32) -> Result

{ + nodes.get(v as usize).copied().ok_or_else(|| Error::Failed(format!("node {v} out of range"))) +} + +fn fmt(p: P) -> String { + format!("({:.4}, {:.4}, {:.4})", p[0], p[1], p[2]) +} + +#[cfg(test)] +mod tests; diff --git a/crates/mesh/src/hex/extrude_to/tests.rs b/crates/mesh/src/hex/extrude_to/tests.rs new file mode 100644 index 0000000..d97552e --- /dev/null +++ b/crates/mesh/src/hex/extrude_to/tests.rs @@ -0,0 +1,178 @@ +use super::*; +use crate::tet::geometry::AnalyticSurface; + +/// Untrimmed analytic faces, by index. +struct Faces(Vec); + +impl Geometry for Faces { + fn project(&self, face: u32, p: P) -> Option

{ + self.0.get(face as usize)?.project(p) + } + fn normal(&self, face: u32, p: P) -> Option

{ + self.0.get(face as usize)?.normal_at(p) + } +} + +/// An `nx` × `ny` grid of quadrilaterals over `[x0, x1] × [y0, y1]` at `z`. +fn grid(x: [f64; 2], y: [f64; 2], z: f64, nx: u32, ny: u32) -> (Vec

, Vec<[u32; 4]>) { + let mut nodes = Vec::new(); + for j in 0..=ny { + for i in 0..=nx { + let (u, v) = (f64::from(i) / f64::from(nx), f64::from(j) / f64::from(ny)); + nodes.push([x[0] + u * (x[1] - x[0]), y[0] + v * (y[1] - y[0]), z]); + } + } + let id = |i: u32, j: u32| j * (nx + 1) + i; + let quads = (0..ny).flat_map(|j| (0..nx).map(move |i| [id(i, j), id(i + 1, j), id(i + 1, j + 1), id(i, j + 1)])).collect(); + (nodes, quads) +} + +fn spec<'a>(nodes: &'a [P], quads: &'a [[u32; 4]], targets: &'a [u32], sides: &'a [u32]) -> ExtrudeSpec<'a> { + ExtrudeSpec { + nodes, + quads, + path: ExtrudePath::Direction([0.0, 0.0, 1.0]), + targets, + sides, + layers: ExtrudeLayers::Count(8), + bias: 1.0, + tolerance: 1e-6, + max_elements: 1_000_000, + } +} + +fn last_layer(m: &Extruded, n_src: usize) -> &[P] { + &m.nodes[m.layers * n_src..] +} + +#[test] +fn a_square_section_reaches_a_cylinder() { + let g = Faces(vec![AnalyticSurface::Cylinder { origin: [0.0, 0.0, 30.0], axis: [1.0, 0.0, 0.0], radius: 20.0 }]); + let (nodes, quads) = grid([-4.0, 4.0], [-4.0, 4.0], 0.0, 4, 4); + let m = extrude_to_surface(&spec(&nodes, &quads, &[0], &[]), &g).unwrap(); + assert_eq!(m.hexes.len(), 16 * 8); + assert_eq!(&m.nodes[..nodes.len()], &nodes[..], "the source nodes are kept"); + for p in last_layer(&m, nodes.len()) { + let r = (p[1].powi(2) + (p[2] - 30.0).powi(2)).sqrt(); + assert!((r - 20.0).abs() < 1e-9, "{p:?} is {r} from the axis"); + } + assert!(m.max_target_gap < 1e-9); + assert!(m.min_scaled_jacobian >= 0.5, "{}", m.min_scaled_jacobian); + assert!(m.gate.passed()); + // Volume of the column under the arc: ∫∫ (30 − √(400 − y²)) dx dy over the square. + let exact = 8.0 * (30.0 * 8.0 - (4.0 * (400.0f64 - 16.0).sqrt() + 400.0 * (0.2f64).asin())); + assert!((m.gate.volume - exact).abs() < 1e-2 * exact, "{} vs {exact}", m.gate.volume); + // Every boundary face is tagged: 16 source, 16 target, 4 × 4 × 8 sides. + assert_eq!(m.boundary.iter().filter(|b| b.role == ExtrudeRole::Source).count(), 16); + assert_eq!(m.boundary.iter().filter(|b| b.role == ExtrudeRole::Target(0)).count(), 16); + assert_eq!(m.boundary.iter().filter(|b| b.role == ExtrudeRole::Side(None)).count(), 128); +} + +#[test] +fn a_square_section_reaches_a_sphere_cap_with_graded_layers() { + let g = Faces(vec![AnalyticSurface::Sphere { center: [0.0, 0.0, -5.0], radius: 25.0 }]); + let (nodes, quads) = grid([-6.0, 6.0], [-6.0, 6.0], 0.0, 6, 6); + let s = ExtrudeSpec { bias: 0.25, layers: ExtrudeLayers::Size(2.0), ..spec(&nodes, &quads, &[0], &[]) }; + let m = extrude_to_surface(&s, &g).unwrap(); + assert_eq!(m.layers, 10, "about 20 mm at 2 mm"); + for p in last_layer(&m, nodes.len()) { + let r = (p[0].powi(2) + p[1].powi(2) + (p[2] + 5.0).powi(2)).sqrt(); + assert!((r - 25.0).abs() < 1e-9, "{p:?}"); + } + // Graded: the last layer is a quarter of the first at the centre node. + let c = 3 * 7 + 3; + let z = |k: usize| m.nodes[k * nodes.len() + c][2]; + let (first, last) = (z(1) - z(0), z(10) - z(9)); + assert!((last / first - 0.25).abs() < 1e-9, "{first} {last}"); + assert!(m.min_scaled_jacobian >= 0.5, "{}", m.min_scaled_jacobian); +} + +#[test] +fn layers_follow_side_faces_that_are_not_straight() { + // Waisted sides x = ±(24 − √(400 − (z − 5)²)), flat sides y = ±4, a tilted target plane. + let x0 = 24.0 - 375.0f64.sqrt(); + let g = Faces(vec![ + AnalyticSurface::Plane { normal: [0.0, -0.1, 1.0], point: [0.0, 0.0, 10.0] }, + AnalyticSurface::Cylinder { origin: [24.0, 0.0, 5.0], axis: [0.0, 1.0, 0.0], radius: 20.0 }, + AnalyticSurface::Cylinder { origin: [-24.0, 0.0, 5.0], axis: [0.0, 1.0, 0.0], radius: 20.0 }, + AnalyticSurface::Plane { normal: [0.0, 1.0, 0.0], point: [0.0, 4.0, 0.0] }, + AnalyticSurface::Plane { normal: [0.0, 1.0, 0.0], point: [0.0, -4.0, 0.0] }, + ]); + let (nodes, quads) = grid([-x0, x0], [-4.0, 4.0], 0.0, 6, 6); + let m = extrude_to_surface(&spec(&nodes, &quads, &[0], &[1, 2, 3, 4]), &g).unwrap(); + let n = nodes.len(); + // Every boundary node of every layer is on its side face. + let on_side = |p: P| { + let waist = 24.0 - (400.0 - (p[2] - 5.0).powi(2)).sqrt(); + (p[0].abs() - waist).abs() < 1e-6 || (p[1].abs() - 4.0).abs() < 1e-9 + }; + for k in 0..=m.layers { + for j in 0..=6u32 { + for i in 0..=6u32 { + if i == 0 || j == 0 || i == 6 || j == 6 { + let p = m.nodes[k * n + (j * 7 + i) as usize]; + assert!(on_side(p), "layer {k} node {p:?}"); + } + } + } + } + for p in last_layer(&m, n) { + assert!((p[2] - (10.0 + 0.1 * p[1])).abs() < 1e-9, "{p:?}"); + } + assert!(m.boundary.iter().any(|b| b.role == ExtrudeRole::Side(Some(1)))); + assert!(m.boundary.iter().all(|b| b.role != ExtrudeRole::Side(None))); + assert!(m.min_scaled_jacobian >= 0.5, "{}", m.min_scaled_jacobian); +} + +#[test] +fn normal_paths_turn_toward_the_target() { + let g = Faces(vec![AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, -10.0] }]); + // Counterclockwise from above: the normal is +Z, the target below. + let (nodes, quads) = grid([0.0, 4.0], [0.0, 4.0], 0.0, 2, 2); + let s = ExtrudeSpec { path: ExtrudePath::Normal, layers: ExtrudeLayers::Count(5), ..spec(&nodes, &quads, &[0], &[]) }; + let m = extrude_to_surface(&s, &g).unwrap(); + assert!(last_layer(&m, nodes.len()).iter().all(|p| (p[2] + 10.0).abs() < 1e-12)); + assert!((m.gate.volume - 160.0).abs() < 1e-9 && m.min_scaled_jacobian > 0.99); +} + +#[test] +fn hostile_input_is_an_error() { + let g = Faces(vec![AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, 10.0] }]); + let (nodes, quads) = grid([0.0, 1.0], [0.0, 1.0], 0.0, 1, 1); + let ok = spec(&nodes, &quads, &[0], &[]); + assert!(extrude_to_surface(&ok, &g).is_ok()); + let bad = [ + ExtrudeSpec { bias: f64::NAN, ..ok.clone() }, + ExtrudeSpec { bias: 0.0, ..ok.clone() }, + ExtrudeSpec { layers: ExtrudeLayers::Count(0), ..ok.clone() }, + ExtrudeSpec { layers: ExtrudeLayers::Count(MAX_LAYERS + 1), ..ok.clone() }, + ExtrudeSpec { layers: ExtrudeLayers::Size(f64::INFINITY), ..ok.clone() }, + ExtrudeSpec { layers: ExtrudeLayers::Size(1e-12), ..ok.clone() }, + ExtrudeSpec { tolerance: -1.0, ..ok.clone() }, + ExtrudeSpec { targets: &[], ..ok.clone() }, + ExtrudeSpec { targets: &[7], ..ok.clone() }, + ExtrudeSpec { path: ExtrudePath::Direction([0.0; 3]), ..ok.clone() }, + // Away from the target. + ExtrudeSpec { path: ExtrudePath::Direction([0.0, 0.0, -1.0]), ..ok.clone() }, + // Parallel to it. + ExtrudeSpec { path: ExtrudePath::Direction([1.0, 0.0, 0.0]), ..ok.clone() }, + ExtrudeSpec { quads: &[[0, 1, 2, 9]], ..ok.clone() }, + ExtrudeSpec { quads: &[[0, 1, 1, 2]], ..ok.clone() }, + ExtrudeSpec { max_elements: 1, layers: ExtrudeLayers::Count(2), ..ok.clone() }, + // A side face the section is not on. + ExtrudeSpec { sides: &[0], ..ok.clone() }, + ]; + for (i, s) in bad.iter().enumerate() { + assert!(extrude_to_surface(s, &g).is_err(), "case {i}"); + } +} + +#[test] +fn layer_parameters_grade_geometrically() { + assert_eq!(parameters(1, 5.0), vec![0.0, 1.0]); + let t = parameters(4, 1.0); + assert!(t.iter().zip([0.0, 0.25, 0.5, 0.75, 1.0]).all(|(a, b)| (a - b).abs() < 1e-15)); + let t = parameters(3, 4.0); + let h: Vec = t.windows(2).map(|w| w[1] - w[0]).collect(); + assert!((h[2] / h[0] - 4.0).abs() < 1e-12 && (h[1] / h[0] - 2.0).abs() < 1e-12, "{h:?}"); +} diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index a9e5e2f..87e7564 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -40,6 +40,7 @@ mod arrangement; mod brep; mod cells; mod decompose; +mod extrude_to; pub(crate) mod geom; mod interval; mod layered; @@ -54,6 +55,7 @@ mod tests; use std::collections::HashMap; pub use cells::{Cell, CellsMesh, GROWTH, MAX_CELLS, hex_mesh_cells, hex_mesh_swept_cells}; +pub use extrude_to::{ExtrudeBoundary, ExtrudeLayers, ExtrudePath, ExtrudeRole, ExtrudeSpec, Extruded, MAX_BIAS, MAX_LAYERS, extrude_to_surface}; pub use geom::Surface; use crate::seeds::EdgeSeed; diff --git a/crates/mesh/src/hex/sweep.rs b/crates/mesh/src/hex/sweep.rs index a8568db..0800c76 100644 --- a/crates/mesh/src/hex/sweep.rs +++ b/crates/mesh/src/hex/sweep.rs @@ -518,7 +518,7 @@ impl Affine { /// Inverse-distance (Shepard, power 4) weights of the ring nodes at `p`, normalized: the /// residual correction at `p` is their weighted sum of the ring's residuals. -fn shepard_weights(p: [f64; 2], ring: &[[f64; 2]]) -> Vec<(usize, f64)> { +pub(super) fn shepard_weights(p: [f64; 2], ring: &[[f64; 2]]) -> Vec<(usize, f64)> { let mut w: Vec<(usize, f64)> = Vec::with_capacity(ring.len()); for (k, s) in ring.iter().enumerate() { let d2 = (p[0] - s[0]).powi(2) + (p[1] - s[1]).powi(2);