From 3587b0de7498e18ace7b96570726b21610b02235 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Wed, 14 Oct 2026 05:34:57 -0400 Subject: [PATCH] mesh: unstructured all-hex octree mesher `fieldcraft_mesh::octree::octree_hex_mesh` meshes any closed part with hex8 only: - an adaptive 3-refinement grid (27-tree) with strong balance and closure, filled with conforming node-marking transition templates (Schneiders; checked over all 48 symmetries); - a core of grid cells kept clear of the boundary (the HybridOctree buffer rule), with non-manifold cleanup and a genus check against the part; - a feature graph from the B-rep (corners, curves, patches) and a partition-first map of the core's boundary onto it, then two boundary layers: a padding layer and a pillow per patch, so every face boundary and sharp edge is captured by element edges; - a Garanzha-type untangling and fitting energy (shape term |J|^3/chi, a volume term and a soft fitting term, deterministic across platforms), a compass search on the worst elements and a certified repair; failures refine the grid near them and rebuild, within bounded attempts and growth. The output passes the certified Jacobian check of every element and a gate that the boundary encloses the elements' volume; boundary quads carry their B-rep faces. Quality is reported with the same scaled Jacobian as mesh.quality. Literature notes for each step are in docs/meshing/literature.md (B1.6-B1.10). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm --- crates/mesh/src/hex/mod.rs | 4 +- crates/mesh/src/lib.rs | 1 + crates/mesh/src/octree/geom.rs | 636 ++++++++++++++++++++++ crates/mesh/src/octree/grid.rs | 376 +++++++++++++ crates/mesh/src/octree/map.rs | 282 ++++++++++ crates/mesh/src/octree/mod.rs | 801 ++++++++++++++++++++++++++++ crates/mesh/src/octree/smooth.rs | 641 ++++++++++++++++++++++ crates/mesh/src/octree/templates.rs | 374 +++++++++++++ crates/mesh/src/octree/tests.rs | 98 ++++ crates/mesh/src/octree/topo.rs | 404 ++++++++++++++ docs/architecture.md | 19 +- docs/meshing/literature.md | 47 ++ 12 files changed, 3676 insertions(+), 7 deletions(-) create mode 100644 crates/mesh/src/octree/geom.rs create mode 100644 crates/mesh/src/octree/grid.rs create mode 100644 crates/mesh/src/octree/map.rs create mode 100644 crates/mesh/src/octree/mod.rs create mode 100644 crates/mesh/src/octree/smooth.rs create mode 100644 crates/mesh/src/octree/templates.rs create mode 100644 crates/mesh/src/octree/tests.rs create mode 100644 crates/mesh/src/octree/topo.rs diff --git a/crates/mesh/src/hex/mod.rs b/crates/mesh/src/hex/mod.rs index 83dae45..5aa8db9 100644 --- a/crates/mesh/src/hex/mod.rs +++ b/crates/mesh/src/hex/mod.rs @@ -37,7 +37,7 @@ mod arrangement; mod brep; mod cells; -mod geom; +pub(crate) mod geom; mod interval; mod layered; mod mapped; @@ -583,7 +583,7 @@ fn finish(cells: &Cells, raw: Raw, method: HexMethod) -> Result { /// `∫ x n_x dA` over a quad4 or quad8 face (outward normal): its contribution to the enclosed /// volume. The 3 × 3 Gauss rule integrates it exactly (degree at most 5 per direction). -fn face_flux(face: &[P]) -> f64 { +pub(crate) fn face_flux(face: &[P]) -> f64 { const G: [(f64, f64); 3] = [(-0.774_596_669_241_483_4, 5.0 / 9.0), (0.0, 8.0 / 9.0), (0.774_596_669_241_483_4, 5.0 / 9.0)]; let mut sum = 0.0; for (r, wr) in G { diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index dab51ad..a5626aa 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -12,6 +12,7 @@ pub mod cdt; pub mod delaunay; pub mod hex; pub mod join; +pub mod octree; pub mod planar; pub mod quad; pub mod quality; diff --git a/crates/mesh/src/octree/geom.rs b/crates/mesh/src/octree/geom.rs new file mode 100644 index 0000000..83f10ec --- /dev/null +++ b/crates/mesh/src/octree/geom.rs @@ -0,0 +1,636 @@ +//! The boundary as the hex mesher sees it: a triangle tree for nearest-point, ray and box +//! queries, the exact inside test, and the feature graph (corners, curves, patches). +//! +//! # Feature graph +//! Every feature edge of the [`BoundarySurface`] (between two B-rep faces, or a crease inside a +//! face) whose dihedral angle reaches the sharp angle is kept; smooth face boundaries (a fillet +//! meeting its faces tangentially) are not features, the faces on either side share a patch. +//! Neither are blades and slits, edges whose normals turn by more than the blade angle (a 3° +//! wedge, a sphere touching a plane): no valid hexahedra can follow them, so the mesh rounds +//! them within the boundary tolerance. *Corners* are vertices where the number of feature edges is not two, or where a +//! curve turns by more than [`CORNER_TURN`]; *curves* are the chains of feature edges between +//! corners, or closed loops; *patches* are the connected sets of triangles between feature edges, +//! so every patch lies in one B-rep face and is bounded by curves. This is the partition the hex +//! mesh's boundary must reproduce (Gao, Shen, Panozzo 2019, §3.3). + +use std::collections::BTreeMap; + +pub(crate) use crate::hex::geom::closest_on_triangle; +use crate::predicates::{Axis, Point3, Sign, orient2d_projected, orient3d, snap_to_domain}; +use crate::surface::BoundarySurface; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +/// A degree-2 feature vertex whose curve turns by more than this (degrees) is a corner. +pub(crate) const CORNER_TURN: f64 = 30.0; +/// Triangles per tree leaf. +const LEAF: usize = 4; + +#[derive(Debug, Clone, Copy, PartialEq)] +pub(crate) struct Aabb { + pub lo: P, + pub hi: P, +} + +impl Aabb { + pub const EMPTY: Aabb = Aabb { lo: [f64::INFINITY; 3], hi: [f64::NEG_INFINITY; 3] }; + pub fn add(&mut self, p: P) { + for ((lo, hi), x) in self.lo.iter_mut().zip(&mut self.hi).zip(p) { + *lo = lo.min(x); + *hi = hi.max(x); + } + } + fn join(&mut self, o: &Aabb) { + self.add(o.lo); + self.add(o.hi); + } + pub fn overlaps(&self, o: &Aabb) -> bool { + (0..3).all(|k| self.lo[k] <= o.hi[k] && o.lo[k] <= self.hi[k]) + } + /// Squared distance from a point to the box (0 inside). + fn dist2(&self, p: P) -> f64 { + (0..3).map(|k| (self.lo[k] - p[k]).max(0.0).max(p[k] - self.hi[k])).map(|d| d * d).sum() + } + /// Whether the segment `a + t (b − a)`, `t ∈ [0, 1]`, may meet the box (slab test, widened by + /// a relative margin so rounding never drops a hit). + fn meets_segment(&self, a: P, b: P) -> bool { + let (mut t0, mut t1) = (0.0f64, 1.0f64); + for k in 0..3 { + let pad = 1e-9 * (self.hi[k] - self.lo[k]).abs().max(1e-300) + 1e-12 * (a[k].abs() + b[k].abs()); + let (lo, hi) = (self.lo[k] - pad, self.hi[k] + pad); + let d = b[k] - a[k]; + if d == 0.0 { + if a[k] < lo || a[k] > hi { + return false; + } + continue; + } + let (mut u, mut v) = ((lo - a[k]) / d, (hi - a[k]) / d); + if u > v { + std::mem::swap(&mut u, &mut v); + } + t0 = t0.max(u); + t1 = t1.min(v); + if t0 > t1 { + return false; + } + } + true + } +} + +#[derive(Debug, Clone, Copy)] +struct Node { + bounds: Aabb, + children: Option, + start: usize, + end: usize, +} + +/// A median-split tree over triangle boxes. +#[derive(Debug, Clone)] +struct Tree { + nodes: Vec, + order: Vec, +} + +impl Tree { + fn build(boxes: &[Aabb]) -> Tree { + let mut order: Vec = (0..boxes.len() as u32).collect(); + let empty = Node { bounds: Aabb::EMPTY, children: None, start: 0, end: 0 }; + let mut nodes = vec![empty]; + let mut work = vec![(0usize, 0usize, boxes.len())]; + while let Some((ni, start, end)) = work.pop() { + let Some(slice) = order.get_mut(start..end) else { continue }; + let (mut bounds, mut centres) = (Aabb::EMPTY, Aabb::EMPTY); + for &t in slice.iter() { + if let Some(b) = boxes.get(t as usize) { + bounds.join(b); + centres.add([0, 1, 2].map(|k| b.lo[k] + b.hi[k])); + } + } + let mut node = Node { bounds, children: None, start, end }; + if end - start > LEAF { + let ext = [0, 1, 2].map(|k| centres.hi[k] - centres.lo[k]); + let axis = if ext[0] >= ext[1] && ext[0] >= ext[2] { + 0 + } else if ext[1] >= ext[2] { + 1 + } else { + 2 + }; + let key = |t: &u32| boxes.get(*t as usize).map_or(0.0, |b| b.lo[axis] + b.hi[axis]); + let mid = (end - start) / 2; + slice.select_nth_unstable_by(mid, |a, b| key(a).total_cmp(&key(b)).then(a.cmp(b))); + let l = nodes.len(); + nodes.push(empty); + nodes.push(empty); + node.children = Some(l); + work.push((l, start, start + mid)); + work.push((l + 1, start + mid, end)); + } + if let Some(n) = nodes.get_mut(ni) { + *n = node; + } + } + Tree { nodes, order } + } + + fn leaf(&self, n: &Node) -> &[u32] { + self.order.get(n.start..n.end).unwrap_or(&[]) + } +} + +pub(crate) fn lerp(a: P, b: P, t: f64) -> P { + [0, 1, 2].map(|k| a[k] + (b[k] - a[k]) * t) +} + +/// The closest point of segment `ab` to `p` and its parameter. +pub(crate) fn closest_on_segment(p: P, a: P, b: P) -> (P, f64) { + let ab = sub(b, a); + let l2 = dot(ab, ab); + let t = if l2 > 0.0 { (dot(sub(p, a), ab) / l2).clamp(0.0, 1.0) } else { 0.0 }; + (lerp(a, b, t), t) +} + +pub(crate) fn dist2(a: P, b: P) -> f64 { + let d = sub(a, b); + dot(d, d) +} + +/// A chain of feature edges. +#[derive(Debug, Clone, PartialEq)] +pub(crate) struct Curve { + /// Surface vertices in order; a closed loop repeats its first vertex at the end. + pub vertices: Vec, + pub closed: bool, + /// The patches on either side (sorted, deduplicated: one entry for a crease seam). + pub patches: Vec, +} + +#[derive(Debug, Clone, PartialEq)] +pub(crate) struct Patch { + pub triangles: Vec, + /// The B-rep faces of its triangles (sorted, deduplicated; one for a well-formed boundary). + pub faces: Vec, + pub area: f64, + /// Twice the area over the length of the bounding curves: the width of a strip, the + /// diameter of a disc, `∞` for a closed surface without curves. + pub width: f64, +} + +/// The boundary with its queries and feature graph. +pub(crate) struct Geometry<'a> { + pub surface: &'a BoundarySurface, + tree: Tree, + pub boxes: Vec, + pub normals: Vec

, + pub tri_patch: Vec, + pub patches: Vec, + pub curves: Vec, + /// Corner vertices, ascending. + pub corners: Vec, + pub bounds: Aabb, + /// Far enough to leave the boundary from any point of its box. + far: f64, +} + +impl<'a> Geometry<'a> { + pub fn new(surface: &'a BoundarySurface, sharp_angle: f64, blade_angle: f64) -> Result> { + let pts = surface.points(); + let tris = surface.triangles(); + let xyz = |v: u32| pts.get(v as usize).copied().unwrap_or([0.0; 3]); + let mut bounds = Aabb::EMPTY; + let mut boxes = Vec::with_capacity(tris.len()); + let mut normals = Vec::with_capacity(tris.len()); + for t in tris { + let mut b = Aabb::EMPTY; + for &v in t { + b.add(xyz(v)); + } + bounds.join(&b); + boxes.push(b); + let n = cross(sub(xyz(t[1]), xyz(t[0])), sub(xyz(t[2]), xyz(t[0]))); + let l = norm(n); + normals.push(if l > 0.0 { n.map(|c| c / l) } else { [0.0; 3] }); + } + if tris.is_empty() || !(0..3).all(|k| bounds.hi[k] >= bounds.lo[k]) { + return Err(Error::Invalid("the boundary has no triangles".into())); + } + let tree = Tree::build(&boxes); + let diag = norm(sub(bounds.hi, bounds.lo)); + let far = 4.0 * diag.max(1e-9) + (0..3).map(|k| bounds.lo[k].abs().max(bounds.hi[k].abs())).fold(0.0, f64::max); + + // Feature edges and corners. + // A boundary between two B-rep faces is kept or dropped as a whole, by the median angle + // along it, so tessellation noise never breaks a curve into pieces; creases inside a + // face go edge by edge. + let face_of = |t: u32| surface.tri_face().get(t as usize).copied().unwrap_or(0); + let mut pair_angles: BTreeMap<(u32, u32), Vec> = BTreeMap::new(); + for e in surface.feature_edges() { + let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1])); + if a != b { + pair_angles.entry((a.min(b), a.max(b))).or_default().push(e.angle); + } + } + let pair_kept: BTreeMap<(u32, u32), bool> = pair_angles + .into_iter() + .map(|(k, mut v)| { + v.sort_by(f64::total_cmp); + let m = v.get(v.len() / 2).copied().unwrap_or(0.0); + (k, m >= sharp_angle && m <= blade_angle) + }) + .collect(); + let kept = |e: &&crate::surface::FeatureEdge| { + let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1])); + if a != b { pair_kept.get(&(a.min(b), a.max(b))).copied().unwrap_or(false) } else { e.angle >= sharp_angle && e.angle <= blade_angle } + }; + let mut edge_set: BTreeMap<(u32, u32), [u32; 2]> = BTreeMap::new(); + for e in surface.feature_edges().iter().filter(kept) { + edge_set.insert((e.vertices[0], e.vertices[1]), e.triangles); + } + let mut adj: BTreeMap> = BTreeMap::new(); + for &(a, b) in edge_set.keys() { + adj.entry(a).or_default().push(b); + adj.entry(b).or_default().push(a); + } + let mut corners: Vec = adj + .iter() + .filter(|(v, n)| { + if n.len() != 2 { + return true; + } + let (p, a, b) = (xyz(**v), xyz(n[0]), xyz(n[1])); + let (u, w) = (sub(p, a), sub(b, p)); + crate::vec3::angle_deg(u, w) > CORNER_TURN + }) + .map(|(v, _)| *v) + .collect(); + corners.sort_unstable(); + + // Patches: flood fill across non-feature edges. + let neighbors = surface.neighbors(); + let mut tri_patch = vec![u32::MAX; tris.len()]; + let mut patches: Vec = Vec::new(); + for seed in 0..tris.len() { + if tri_patch[seed] != u32::MAX { + continue; + } + let id = patches.len() as u32; + let mut patch = Patch { triangles: Vec::new(), faces: Vec::new(), area: 0.0, width: f64::INFINITY }; + tri_patch[seed] = id; + let mut stack = vec![seed as u32]; + while let Some(t) = stack.pop() { + patch.triangles.push(t); + let tri = tris.get(t as usize).copied().unwrap_or([0; 3]); + let n = neighbors.get(t as usize).copied().unwrap_or([u32::MAX; 3]); + for k in 0..3 { + let (a, b) = (tri[k], tri[(k + 1) % 3]); + if edge_set.contains_key(&(a.min(b), a.max(b))) { + continue; + } + if let Some(slot) = tri_patch.get_mut(n[k] as usize) + && *slot == u32::MAX + { + *slot = id; + stack.push(n[k]); + } + } + } + patch.triangles.sort_unstable(); + for &t in &patch.triangles { + let tri = tris.get(t as usize).copied().unwrap_or([0; 3]); + patch.area += 0.5 * norm(cross(sub(xyz(tri[1]), xyz(tri[0])), sub(xyz(tri[2]), xyz(tri[0])))); + patch.faces.push(surface.tri_face().get(t as usize).copied().unwrap_or(0)); + } + patch.faces.sort_unstable(); + patch.faces.dedup(); + let mut perimeter = 0.0; + for &t in &patch.triangles { + let tri = tris.get(t as usize).copied().unwrap_or([0; 3]); + for k in 0..3 { + let (a, b) = (tri[k], tri[(k + 1) % 3]); + if edge_set.contains_key(&(a.min(b), a.max(b))) { + perimeter += norm(sub(xyz(b), xyz(a))); + } + } + } + if perimeter > 0.0 { + patch.width = 2.0 * patch.area / perimeter; + } + patches.push(patch); + } + + // Curves: chains from corners, then closed loops. + let mut used: BTreeMap<(u32, u32), bool> = edge_set.keys().map(|k| (*k, false)).collect(); + let is_corner = |v: u32| corners.binary_search(&v).is_ok(); + let mut curves: Vec = Vec::new(); + let trace = |start: u32, next: u32, used: &mut BTreeMap<(u32, u32), bool>, curves: &mut Vec| { + let mut verts = vec![start]; + let (mut prev, mut cur) = (start, next); + loop { + let key = (prev.min(cur), prev.max(cur)); + if let Some(u) = used.get_mut(&key) { + *u = true; + } + verts.push(cur); + if cur == start || is_corner(cur) || verts.len() > edge_set.len() + 2 { + break; + } + let nxt = adj.get(&cur).and_then(|n| n.iter().copied().find(|&w| w != prev && used.get(&(cur.min(w), cur.max(w))) == Some(&false))); + match nxt { + Some(w) => { + prev = cur; + cur = w; + } + None => break, + } + } + let closed = verts.first() == verts.last() && verts.len() > 2; + let mut ps: Vec = Vec::new(); + for w in verts.windows(2) { + if let Some(t) = edge_set.get(&(w[0].min(w[1]), w[0].max(w[1]))) { + ps.extend(t.iter().filter_map(|&t| tri_patch.get(t as usize).copied())); + } + } + ps.sort_unstable(); + ps.dedup(); + curves.push(Curve { vertices: verts, closed, patches: ps }); + }; + for &c in &corners { + let ns = adj.get(&c).cloned().unwrap_or_default(); + for n in ns { + if used.get(&(c.min(n), c.max(n))) == Some(&false) { + trace(c, n, &mut used, &mut curves); + } + } + } + let keys: Vec<(u32, u32)> = used.keys().copied().collect(); + for k in keys { + if used.get(&k) == Some(&false) { + trace(k.0, k.1, &mut used, &mut curves); + } + } + Ok(Geometry { surface, tree, boxes, normals, tri_patch, patches, curves, corners, bounds, far }) + } + + pub fn xyz(&self, v: u32) -> P { + self.surface.points().get(v as usize).copied().unwrap_or([0.0; 3]) + } + + /// The centroid of triangle `t`. + pub fn triangle_centre(&self, t: u32) -> P { + let [a, b, c] = self.triangle(t); + [0, 1, 2].map(|k| (a[k] + b[k] + c[k]) / 3.0) + } + + pub fn triangle(&self, t: u32) -> [P; 3] { + let tri = self.surface.triangles().get(t as usize).copied().unwrap_or([0; 3]); + tri.map(|v| self.xyz(v)) + } + + /// The closest boundary point to `p` among triangles that pass `keep`: (triangle, point, + /// squared distance). + pub fn nearest(&self, p: P, keep: impl Fn(u32) -> bool) -> Option<(u32, P, f64)> { + let mut best: Option<(u32, P, f64)> = None; + let mut stack = vec![0usize]; + while let Some(n) = stack.pop() { + let Some(node) = self.tree.nodes.get(n) else { continue }; + let bd = node.bounds.dist2(p); + if best.is_some_and(|b| bd >= b.2) { + continue; + } + match node.children { + Some(l) => { + let (dl, dr) = ( + self.tree.nodes.get(l).map_or(f64::INFINITY, |x| x.bounds.dist2(p)), + self.tree.nodes.get(l + 1).map_or(f64::INFINITY, |x| x.bounds.dist2(p)), + ); + if dl <= dr { + stack.extend([l + 1, l]); + } else { + stack.extend([l, l + 1]); + } + } + None => { + for &t in self.tree.leaf(node) { + if !keep(t) { + continue; + } + let [a, b, c] = self.triangle(t); + let q = closest_on_triangle(p, a, b, c); + let d = dist2(p, q); + if best.is_none_or(|x| d < x.2 || (d == x.2 && t < x.0)) { + best = Some((t, q, d)); + } + } + } + } + } + best + } + + /// The triangles whose boxes overlap `b`, ascending. + pub fn overlapping(&self, b: &Aabb) -> Vec { + let mut out = Vec::new(); + let mut stack = vec![0usize]; + while let Some(n) = stack.pop() { + let Some(node) = self.tree.nodes.get(n) else { continue }; + if !node.bounds.overlaps(b) { + continue; + } + match node.children { + Some(l) => stack.extend([l, l + 1]), + None => out.extend(self.tree.leaf(node).iter().copied().filter(|&t| self.boxes.get(t as usize).is_some_and(|x| x.overlaps(b)))), + } + } + out.sort_unstable(); + out + } + + /// The first triangle hit by the ray from `p` along `d` (unit) within `max`, with its + /// distance (floating point; used to measure, never to decide). + pub fn ray(&self, p: P, d: P, max: f64, skip: u32) -> Option<(u32, f64)> { + let q = [0, 1, 2].map(|k| p[k] + d[k] * max); + let mut best: Option<(u32, f64)> = None; + let mut stack = vec![0usize]; + while let Some(n) = stack.pop() { + let Some(node) = self.tree.nodes.get(n) else { continue }; + if !node.bounds.meets_segment(p, q) { + continue; + } + match node.children { + Some(l) => stack.extend([l, l + 1]), + None => { + for &t in self.tree.leaf(node) { + if t == skip { + continue; + } + let [a, b, c] = self.triangle(t); + if let Some(s) = ray_triangle(p, d, a, b, c) + && s > 0.0 + && s <= max + && best.is_none_or(|x| s < x.1) + { + best = Some((t, s)); + } + } + } + } + } + best + } + + /// Whether `p` is strictly inside the solid, decided exactly: the parity of the boundary + /// crossings of a segment from `p` to a far point, tried along fixed directions until one + /// meets no edge or vertex and does not start on the boundary. `None` when every direction is + /// degenerate (p lies on the boundary) or `p` is not admissible. + pub fn inside(&self, p: P) -> Option { + const DIRS: [P; 7] = [ + [0.5773, 0.6123, 0.5401], + [-0.6301, 0.4372, 0.6417], + [0.3989, -0.7068, 0.5839], + [0.7211, 0.4109, -0.5577], + [-0.4402, -0.6005, -0.6675], + [0.1351, 0.9013, -0.4116], + [-0.8803, 0.1237, 0.4583], + ]; + let pp = Point3::new(snap_to_domain(p[0]), snap_to_domain(p[1]), snap_to_domain(p[2])).ok()?; + 'dirs: for d in DIRS { + let q = [0, 1, 2].map(|k| snap_to_domain(p[k] + d[k] * self.far)); + let Ok(qq) = Point3::new(q[0], q[1], q[2]) else { continue }; + let mut parity = false; + let mut stack = vec![0usize]; + while let Some(n) = stack.pop() { + let Some(node) = self.tree.nodes.get(n) else { continue }; + if !node.bounds.meets_segment(p, q) { + continue; + } + match node.children { + Some(l) => stack.extend([l, l + 1]), + None => { + for &t in self.tree.leaf(node) { + match self.crossing(&pp, &qq, t) { + Some(true) => parity = !parity, + Some(false) => {} + None => continue 'dirs, + } + } + } + } + } + return Some(parity); + } + None + } + + /// Whether segment `pq` crosses triangle `t` properly (`Some(true)`), misses it + /// (`Some(false)`), or touches it degenerately (`None`): exact. + fn crossing(&self, p: &Point3, q: &Point3, t: u32) -> Option { + let tri = self.surface.triangles().get(t as usize)?; + let a = self.surface.point3(tri[0])?; + let b = self.surface.point3(tri[1])?; + let c = self.surface.point3(tri[2])?; + let (sp, sq) = (orient3d(&a, &b, &c, p), orient3d(&a, &b, &c, q)); + if sq.is_zero() { + // The far end is outside the boundary's box, so never on a triangle; a segment in + // the plane is ambiguous. + return None; + } + if sp.is_zero() { + // p in the triangle's plane: the segment leaves the plane at p, so it touches the + // triangle only if p lies in it. + let axis = [Axis::Z, Axis::X, Axis::Y].into_iter().find(|&ax| !orient2d_projected(&a, &b, &c, ax).is_zero())?; + let o = orient2d_projected(&a, &b, &c, axis); + let inside = [(&a, &b), (&b, &c), (&c, &a)].into_iter().all(|(u, v)| orient2d_projected(u, v, p, axis) != -o); + return if inside { None } else { Some(false) }; + } + if sp == sq { + return Some(false); + } + let s = [orient3d(p, q, &a, &b), orient3d(p, q, &b, &c), orient3d(p, q, &c, &a)]; + if s.iter().any(|x| x.is_zero()) { + let nz: Vec = s.iter().copied().filter(|x| !x.is_zero()).collect(); + return if nz.windows(2).all(|w| w[0] == w[1]) { None } else { Some(false) }; + } + Some(s[0] == s[1] && s[1] == s[2]) + } + + /// The closest point of curve `c` to `p`: (point, segment index, parameter on it). + pub fn nearest_on_curve(&self, c: u32, p: P) -> Option<(P, usize, f64)> { + let curve = self.curves.get(c as usize)?; + let mut best: Option<(P, usize, f64, f64)> = None; + for (i, w) in curve.vertices.windows(2).enumerate() { + let (q, t) = closest_on_segment(p, self.xyz(w[0]), self.xyz(w[1])); + let d = dist2(p, q); + if best.is_none_or(|b| d < b.3) { + best = Some((q, i, t, d)); + } + } + best.map(|(q, i, t, _)| (q, i, t)) + } +} + +/// Möller–Trumbore ray–triangle distance (floating point). +fn ray_triangle(o: P, d: P, a: P, b: P, c: P) -> Option { + let (e1, e2) = (sub(b, a), sub(c, a)); + let h = cross(d, e2); + let det = dot(e1, h); + if det.abs() < 1e-300 { + return None; + } + let inv = 1.0 / det; + let s = sub(o, a); + let u = inv * dot(s, h); + if !(-1e-12..=1.0 + 1e-12).contains(&u) { + return None; + } + let q = cross(s, e1); + let v = inv * dot(d, q); + if v < -1e-12 || u + v > 1.0 + 1e-12 { + return None; + } + Some(inv * dot(e2, q)) +} + +#[cfg(test)] +pub(crate) mod tests { + use super::*; + use crate::surface::SurfaceOptions; + + /// A closed axis-aligned box surface, one B-rep face per side. + pub(crate) fn box_surface(lo: P, hi: P) -> BoundarySurface { + let c = |i: usize| [if i & 1 == 0 { lo[0] } else { hi[0] }, if i & 2 == 0 { lo[1] } else { hi[1] }, if i & 4 == 0 { lo[2] } else { hi[2] }]; + let pts: Vec

= (0..8).map(c).collect(); + let quads = [[0, 2, 3, 1], [4, 5, 7, 6], [0, 1, 5, 4], [2, 6, 7, 3], [0, 4, 6, 2], [1, 3, 7, 5]]; + let mut tris = Vec::new(); + let mut face = Vec::new(); + for (f, q) in quads.iter().enumerate() { + tris.push([q[0], q[1], q[2]]); + tris.push([q[0], q[2], q[3]]); + face.extend([f as u32; 2]); + } + BoundarySurface::new(&pts, &tris, &face, &SurfaceOptions::default()).unwrap() + } + + #[test] + fn box_feature_graph_inside_and_nearest() { + let s = box_surface([0.0; 3], [2.0, 3.0, 4.0]); + let g = Geometry::new(&s, 30.0, 160.0).unwrap(); + assert_eq!(g.corners.len(), 8); + assert_eq!(g.curves.len(), 12); + assert_eq!(g.patches.len(), 6); + assert!(g.curves.iter().all(|c| c.vertices.len() == 2 && c.patches.len() == 2 && !c.closed)); + assert_eq!(g.inside([1.0, 1.5, 2.0]), Some(true)); + assert_eq!(g.inside([3.0, 1.5, 2.0]), Some(false)); + assert_eq!(g.inside([1.0, 1.5, 0.0]), None, "on the boundary"); + assert_eq!(g.inside([0.0, 0.0, 0.0]), None, "a corner"); + let (_, q, d) = g.nearest([1.0, 1.0, -1.0], |_| true).unwrap(); + assert!((d - 1.0).abs() < 1e-12 && q[2] == 0.0); + let top = g.tri_patch[2]; + let (_, q, _) = g.nearest([1.0, 1.0, -1.0], |t| g.tri_patch[t as usize] == top).unwrap(); + assert_eq!(q[2], 4.0, "restricted to the top patch"); + let (t, s) = g.ray([1.0, 1.5, 2.0], [0.0, 0.0, 1.0], 10.0, u32::MAX).unwrap(); + assert!((s - 2.0).abs() < 1e-12 && g.tri_patch[t as usize] == top); + } +} diff --git a/crates/mesh/src/octree/grid.rs b/crates/mesh/src/octree/grid.rs new file mode 100644 index 0000000..33b6995 --- /dev/null +++ b/crates/mesh/src/octree/grid.rs @@ -0,0 +1,376 @@ +//! The adaptive grid: a box of root cubes, each refined as a 27-tree (a leaf splits into +//! 3 × 3 × 3), strongly balanced and closed under the template patterns, then cut into a +//! conforming hexahedral mesh by [`super::templates`]. +//! +//! # Balance and closure +//! *Strong balance*: leaves that share a vertex differ by at most one level, so a leaf's finer +//! neighbours are exactly one level finer and a leaf next to a coarser one never has a finer one +//! on that side. A vertex of a leaf at level `L` is *marked* when a finer leaf touches it, that +//! is, when one of the eight level-`L` cells around the vertex is split. *Closure*: a leaf whose +//! marks fit no template (two diagonal corners of a face, three of a face, …) is split too. +//! Balancing and closing repeat until neither splits anything; both only split, so this ends, +//! and every round runs over the leaves in sorted order, so the result is the same on every run. +//! +//! # Coordinates +//! Everything is integer. A cell at level `L` has index `i` (per axis, in units of its own edge) +//! and spans `6 · 3^(D − L)` lattice units, `D` the deepest level allowed; template vertices are +//! in sixths of a leaf, so every vertex is a lattice point and vertices shared by neighbours are +//! found by exact lookup. + +use std::collections::HashMap; + +use super::templates::{self, Pattern, Sym}; +use crate::{Error, Result}; + +/// A cell: its level and index (in units of its edge at that level). +pub(crate) type Key = (u8, [i64; 3]); + +/// The deepest level a tree may have (`3^12` lattice steps per root cell edge, × 6). +pub(crate) const MAX_LEVEL: u8 = 12; + +fn pow3(n: u8) -> i64 { + 3i64.pow(u32::from(n)) +} + +/// What covers a cell-sized region at some level. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +enum Cover { + /// A leaf (at that level or coarser). + Leaf(Key), + /// The region is split at that level (finer leaves inside). + Split, + /// Outside the root box. + Outside, +} + +/// The tree. `nodes` maps every cell that exists to whether it is a leaf. +#[derive(Debug, Clone)] +pub(crate) struct Tree { + /// Root cubes per axis. + pub roots: [i64; 3], + /// The deepest level allowed. + pub depth: u8, + nodes: HashMap, + leaves: usize, + max_leaves: usize, +} + +/// The sort order of [`Tree::sorted_leaves`]: level, then z, y, x. +fn order(k: &Key) -> (u8, i64, i64, i64) { + (k.0, k.1[2], k.1[1], k.1[0]) +} + +/// The 26 neighbour offsets, in a fixed order. +fn offsets() -> impl Iterator { + (0..27).map(|k| [k % 3 - 1, (k / 3) % 3 - 1, k / 9 - 1]).filter(|d| *d != [0, 0, 0]) +} + +impl Tree { + /// One leaf per root cube. + pub fn new(roots: [i64; 3], depth: u8, max_leaves: usize) -> Result { + if roots.iter().any(|&n| n < 1) || depth > MAX_LEVEL { + return Err(Error::Invalid(format!("grid of {roots:?} root cells to depth {depth}"))); + } + let count = roots.iter().try_fold(1usize, |a, &n| usize::try_from(n).ok().and_then(|n| a.checked_mul(n))); + let count = count.filter(|&c| c <= max_leaves).ok_or_else(|| Error::Limit(format!("{roots:?} root cells exceed {max_leaves} cells")))?; + let mut nodes = HashMap::with_capacity(count); + for z in 0..roots[2] { + for y in 0..roots[1] { + for x in 0..roots[0] { + nodes.insert((0, [x, y, z]), true); + } + } + } + Ok(Tree { roots, depth, nodes, leaves: count, max_leaves }) + } + + pub fn leaf_count(&self) -> usize { + self.leaves + } + + /// Whether `key` is a leaf. + pub fn is_leaf(&self, key: &Key) -> bool { + self.nodes.get(key).copied().unwrap_or(false) + } + + fn inside(&self, (level, i): Key) -> bool { + let n = pow3(level); + (0..3).all(|k| i[k] >= 0 && i[k] < self.roots[k] * n) + } + + /// What covers the level-`L` cell `key`. + fn cover(&self, key: Key) -> Cover { + if !self.inside(key) { + return Cover::Outside; + } + let (level, i) = key; + for l in (0..=level).rev() { + let s = pow3(level - l); + let k = (l, [i[0].div_euclid(s), i[1].div_euclid(s), i[2].div_euclid(s)]); + match self.nodes.get(&k) { + Some(true) => return Cover::Leaf(k), + Some(false) => return Cover::Split, + None => {} + } + } + Cover::Outside + } + + /// Split a leaf into its 27 children. + pub fn split(&mut self, key: Key) -> Result<()> { + let level = key.0; + if !self.is_leaf(&key) { + return Ok(()); + } + if level >= self.depth { + return Err(Error::Limit(format!("refinement needs more than {} levels", self.depth))); + } + if self.leaves + 26 > self.max_leaves { + return Err(Error::Limit(format!("the grid needs more than {} cells", self.max_leaves))); + } + self.nodes.insert(key, false); + for child in Self::children(key) { + self.nodes.insert(child, true); + } + self.leaves += 26; + Ok(()) + } + + /// The leaves in sorted order (level, then z, y, x). + pub fn sorted_leaves(&self) -> Vec { + let mut v: Vec = self.nodes.iter().filter(|(_, leaf)| **leaf).map(|(k, _)| *k).collect(); + v.sort_by_key(order); + v + } + + /// The children of a split cell. + pub fn children(key: Key) -> impl Iterator { + let (level, i) = key; + (0..27).map(move |k| (level + 1, [3 * i[0] + k % 3, 3 * i[1] + (k / 3) % 3, 3 * i[2] + k / 9])) + } + + /// Split until leaves sharing a vertex differ by at most one level, checking only the leaves + /// `work` (in sorted order) and the children of the cells it splits: a split only makes + /// neighbours finer, so only new leaves can be too fine for a neighbour, and the result is + /// the unique coarsest balanced refinement whatever the order. Every split cell goes to `log`. + fn balance_from(&mut self, mut work: Vec, log: &mut Vec) -> Result<()> { + work.reverse(); + while let Some(key) = work.pop() { + if !self.is_leaf(&key) || key.0 < 2 { + continue; + } + let (level, i) = key; + for d in offsets() { + if let Cover::Leaf(k) = self.cover((level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]])) + && k.0 + 1 < level + { + self.split(k)?; + log.push(k); + let mut kids: Vec = Self::children(k).collect(); + kids.reverse(); + work.extend(kids); + work.push(key); + break; + } + } + } + Ok(()) + } + + /// The leaves whose marks a split of `key` may change: its children and the leaves beside it + /// on its level. + fn touched_by(&self, key: Key, out: &mut Vec) { + out.extend(Self::children(key)); + let (level, i) = key; + out.extend(offsets().map(|d| (level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]])).filter(|k| self.is_leaf(k))); + } + + /// The marked corners of a leaf (bit `c` for corner `c`: bit 0 x, bit 1 y, bit 2 z). + pub fn marks(&self, (level, i): Key) -> u8 { + let mut mask = 0u8; + for c in 0..8usize { + let o = [(c & 1) as i64, ((c >> 1) & 1) as i64, ((c >> 2) & 1) as i64]; + let split = (0..8usize).any(|e| { + let n = [i[0] + o[0] - 1 + (e & 1) as i64, i[1] + o[1] - 1 + ((e >> 1) & 1) as i64, i[2] + o[2] - 1 + ((e >> 2) & 1) as i64]; + self.nodes.get(&(level, n)) == Some(&false) + }); + if split { + mask |= 1 << c; + } + } + mask + } + + /// Balance, then split every leaf whose marks fit no template (and every leaf marked at all + /// eight corners, whose template is the split), until nothing changes. `max_rounds` bounds the + /// rounds. Each round re-checks only the leaves the previous splits can have changed, so a + /// round costs what it splits, not the size of the tree. + pub fn close(&mut self, max_rounds: usize, skip: &dyn Fn(Key) -> bool) -> Result<()> { + let mut fresh = self.sorted_leaves(); + let mut touched = fresh.clone(); + for _ in 0..max_rounds { + let mut log = Vec::new(); + self.balance_from(std::mem::take(&mut fresh), &mut log)?; + for &k in &log { + self.touched_by(k, &mut touched); + } + touched.retain(|k| self.is_leaf(k)); + touched.sort_by_key(order); + touched.dedup(); + let bad: Vec = touched + .drain(..) + .filter(|&k| !skip(k)) + .filter(|&k| !matches!(templates::classify(self.marks(k)), Some(Pattern::None | Pattern::Corner | Pattern::Edge | Pattern::Face))) + .collect(); + if bad.is_empty() { + return Ok(()); + } + for &k in &bad { + self.split(k)?; + } + for &k in &bad { + fresh.extend(Self::children(k)); + self.touched_by(k, &mut touched); + } + fresh.sort_by_key(order); + } + Err(Error::Limit(format!("grid closure did not settle in {max_rounds} rounds"))) + } + + /// Lattice units per edge of a level-`L` cell. + pub fn span(&self, level: u8) -> i64 { + 6 * pow3(self.depth - level) + } + + /// Cut every leaf by its template into a conforming mesh on the integer lattice. Fails if a + /// leaf's marks fit no template (call [`Tree::close`] first). + pub fn mesh(&self, skip: &dyn Fn(Key) -> bool) -> Result { + let syms: Vec = templates::symmetries(); + let canon: Vec<(Pattern, Vec)> = + [Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face].into_iter().map(|p| (p, templates::template(p))).collect(); + let mut out = LatticeMesh::default(); + let mut index: HashMap<[i64; 3], u32> = HashMap::new(); + for key in self.sorted_leaves() { + if skip(key) { + continue; + } + let mask = self.marks(key); + let pattern = templates::classify(mask).filter(|p| *p != Pattern::Full); + let (pattern, sym) = pattern + .and_then(|p| templates::orientation(&syms, p, mask).map(|s| (p, s))) + .ok_or_else(|| Error::Failed(format!("grid cell {key:?} has marks {mask:#010b} that fit no template")))?; + let hexes = canon.iter().find(|(p, _)| *p == pattern).map(|(_, h)| h.as_slice()).unwrap_or(&[]); + let (level, i) = key; + let unit = pow3(self.depth - level); + for h in hexes { + let mut g = h.map(|p| { + let q = sym.apply(p); + [0, 1, 2].map(|k| (6 * i[k] + q[k]) * unit) + }); + if templates::corner0_orientation(&g) < 0 { + g = templates::reversed(g); + } + let mut ids = [0u32; 8]; + for (id, p) in ids.iter_mut().zip(g) { + *id = match index.get(&p) { + Some(&v) => v, + None => { + let v = u32::try_from(out.points.len()).map_err(|_| Error::Limit("too many grid points".into()))?; + index.insert(p, v); + out.points.push(p); + v + } + }; + } + out.hexes.push(ids); + out.level.push(level); + } + } + Ok(out) + } +} + +/// A conforming hexahedral mesh on the integer lattice. +#[derive(Debug, Clone, Default)] +pub(crate) struct LatticeMesh { + pub points: Vec<[i64; 3]>, + pub hexes: Vec<[u32; 8]>, + /// The tree level of the leaf each hexahedron came from. + pub level: Vec, +} + +#[cfg(test)] +mod tests { + use std::collections::HashMap; + + use super::*; + use crate::quality::{Validity, element_validity}; + use fieldcraft_model::ElementKind; + + /// A tiny deterministic generator for test trees. + struct Lcg(u64); + impl Lcg { + fn next(&mut self) -> u64 { + self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + self.0 >> 33 + } + } + + fn check_conforming(tree: &Tree) { + let m = tree.mesh(&|_| false).unwrap(); + let mut faces: HashMap<[u32; 4], usize> = HashMap::new(); + let mut volume = 0.0; + for h in &m.hexes { + let xyz: Vec<[f64; 3]> = h.iter().map(|&v| m.points[v as usize].map(|c| c as f64)).collect(); + let c = element_validity(ElementKind::Hex8, &xyz).unwrap(); + assert_eq!(c.validity, Validity::Valid); + volume += c.volume; + for f in ElementKind::Hex8.faces() { + let mut q = [h[f[0]], h[f[1]], h[f[2]], h[f[3]]]; + q.sort(); + *faces.entry(q).or_default() += 1; + } + } + let side = tree.span(0) as f64; + let expected = side.powi(3) * (tree.roots[0] * tree.roots[1] * tree.roots[2]) as f64; + assert!((volume - expected).abs() < 1e-6 * expected, "volume {volume} vs {expected}"); + let hi = tree.roots.map(|n| n * tree.span(0)); + for (q, n) in faces { + let on_box = (0..3).any(|k| q.iter().all(|&v| m.points[v as usize][k] == 0) || q.iter().all(|&v| m.points[v as usize][k] == hi[k])); + assert_eq!(n, if on_box { 1 } else { 2 }, "face {q:?}"); + } + } + + #[test] + fn random_trees_close_into_conforming_valid_meshes() { + let mut rng = Lcg(7); + for case in 0..12 { + let roots = [1 + (rng.next() % 2) as i64, 1 + (rng.next() % 2) as i64, 1]; + let mut tree = Tree::new(roots, 3, 1_000_000).unwrap(); + for _ in 0..(1 + rng.next() % 6) { + let leaves = tree.sorted_leaves(); + let k = leaves[(rng.next() as usize) % leaves.len()]; + if k.0 < 3 { + tree.split(k).unwrap(); + } + } + tree.close(100, &|_| false).unwrap(); + for k in tree.sorted_leaves() { + assert!(templates::classify(tree.marks(k)).is_some(), "case {case}"); + } + check_conforming(&tree); + } + } + + #[test] + fn caps_are_errors() { + assert!(Tree::new([0, 1, 1], 2, 10).is_err()); + assert!(Tree::new([10, 10, 10], 2, 999).is_err()); + let mut t = Tree::new([1, 1, 1], 1, 100).unwrap(); + t.split((0, [0, 0, 0])).unwrap(); + assert!(t.split((1, [0, 0, 0])).is_err(), "deeper than the depth"); + let mut t = Tree::new([1, 1, 1], 3, 30).unwrap(); + t.split((0, [0, 0, 0])).unwrap(); + assert!(t.split((1, [0, 0, 0])).is_err(), "more cells than allowed"); + } +} diff --git a/crates/mesh/src/octree/map.rs b/crates/mesh/src/octree/map.rs new file mode 100644 index 0000000..ac176c9 --- /dev/null +++ b/crates/mesh/src/octree/map.rs @@ -0,0 +1,282 @@ +//! The topological map from the boundary's feature graph (corners, curves, patches) onto the +//! core's boundary quads, after X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based +//! hexahedral meshing*, CGF 38(5), 2019, §3.3 (paper only; the authors' code has no licence and +//! is not read). Gao et al. trace the curves first and read the patches off them; here the +//! patches come first and the curves are their borders, which keeps every region on the side of +//! the part it faces: +//! +//! 1. **Label.** Each quad takes the patch its outward normal meets first (within two quad +//! sizes), or else the patch nearest its centre. +//! 2. **Clean.** A quad whose label differs from all its edge neighbours takes their majority +//! label; then every patch keeps only its largest connected group of quads, the other +//! groups joining the label they border most. Both repeat a bounded number of times. +//! 3. **Check.** The map must be a bijection: every patch has a region, and two regions border +//! each other exactly where the part has a curve between their patches. +//! 4. **Features.** Region borders are the curves (each border edge takes the curve between its +//! two patches nearest to it); a vertex where three or more regions meet is a corner and takes +//! the nearest corner of the part around those patches. Corners that are not region junctions +//! (where one curve turns sharply) go to the nearest vertex on their curve's border. +//! +//! Any failure names where the grid is too coarse for the features; the caller refines there and +//! tries again (§3.3, "Bijective map"). + +use std::collections::{BTreeMap, BTreeSet}; + +use super::geom::{Geometry, closest_on_segment, dist2}; +use super::smooth::Constraint; +use super::topo::quad_edges; +use crate::vec3::{P, cross, norm, sub}; + +/// Rounds of label cleaning. +const CLEAN_ROUNDS: usize = 16; + +/// Where and why the map failed. +#[derive(Debug, Clone)] +pub(crate) struct MapFailure { + pub at: Vec

, + pub why: String, +} + +/// The map: a constraint per quad vertex and a patch per quad. +#[derive(Debug, Clone)] +pub(crate) struct FeatureMap { + pub role: BTreeMap, + pub quad_patch: Vec, +} + +/// Map the feature graph onto `quads` (the core's boundary, outward) over `xyz`. +pub(crate) fn map_features(g: &Geometry, quads: &[[u32; 4]], xyz: &dyn Fn(u32) -> P) -> Result { + let n = quads.len(); + let centre = |q: usize| { + let mut c = [0.0; 3]; + for &v in &quads[q] { + let p = xyz(v); + for k in 0..3 { + c[k] += 0.25 * p[k]; + } + } + c + }; + let edges = quad_edges(quads); + let mut nbr: Vec> = vec![Vec::new(); n]; + for qs in edges.values() { + if let [a, b] = qs.as_slice() { + nbr[*a as usize].push(*b); + nbr[*b as usize].push(*a); + } + } + + // 1. Label: the patch the quad faces (a ray along its outward normal), else the nearest. + let mut label: Vec = (0..n) + .map(|q| { + let c = centre(q); + let p = quads[q].map(xyz); + let nrm = cross(sub(p[2], p[0]), sub(p[3], p[1])); + let l = norm(nrm); + let reach = 2.0 * (norm(sub(p[2], p[0])) + norm(sub(p[3], p[1]))); + let hit = if l > 0.0 { g.ray(c, nrm.map(|x| x / l), reach, u32::MAX).map(|(t, _)| t) } else { None }; + hit.or_else(|| g.nearest(c, |_| true).map(|(t, _, _)| t)).and_then(|t| g.tri_patch.get(t as usize).copied()).unwrap_or(0) + }) + .collect(); + + // A patch no quad faces (a narrow face the core's boundary passes by) takes the quad + // nearest to it, which the cleaning keeps. + let mut seeded: BTreeSet = BTreeSet::new(); + let have: BTreeSet = label.iter().copied().collect(); + for p in 0..g.patches.len() as u32 { + if have.contains(&p) { + continue; + } + let best = (0..n) + .filter(|q| !seeded.contains(q)) + .filter_map(|q| g.nearest(centre(q), |t| g.tri_patch.get(t as usize) == Some(&p)).map(|x| (q, x.2))) + .min_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0))); + if let Some((q, _)) = best { + label[q] = p; + seeded.insert(q); + } + } + // 2. Clean. + for _ in 0..CLEAN_ROUNDS { + let mut changed = false; + for q in 0..n { + let mut count: BTreeMap = BTreeMap::new(); + for &m in &nbr[q] { + *count.entry(label[m as usize]).or_default() += 1; + } + let same = count.get(&label[q]).copied().unwrap_or(0); + if same == 0 + && !nbr[q].is_empty() + && !seeded.contains(&q) + && let Some((&best, _)) = count.iter().filter(|(l, _)| **l != label[q]).max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0))) + { + label[q] = best; + changed = true; + } + } + let mut comp = vec![usize::MAX; n]; + let mut comps: Vec<(u32, Vec)> = Vec::new(); + for s in 0..n { + if comp[s] != usize::MAX { + continue; + } + let id = comps.len(); + let mut members = vec![s]; + comp[s] = id; + let mut i = 0; + while i < members.len() { + let q = members[i]; + i += 1; + for &m in &nbr[q] { + let m = m as usize; + if comp[m] == usize::MAX && label[m] == label[s] { + comp[m] = id; + members.push(m); + } + } + } + comps.push((label[s], members)); + } + let mut largest: BTreeMap = BTreeMap::new(); + for (i, (l, m)) in comps.iter().enumerate() { + let e = largest.entry(*l).or_insert(i); + if comps[*e].1.len() < m.len() { + *e = i; + } + } + for (i, (l, members)) in comps.iter().enumerate() { + if largest.get(l) == Some(&i) { + continue; + } + let mut count: BTreeMap = BTreeMap::new(); + for &q in members { + for &m in &nbr[q] { + if comp[m as usize] != i { + *count.entry(label[m as usize]).or_default() += 1; + } + } + } + if let Some((&best, _)) = count.iter().max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0))) { + for &q in members { + label[q] = best; + } + changed = true; + } + } + if !changed { + break; + } + } + + // 3. Check: every patch has a region; borders match curves. + let present: BTreeSet = label.iter().copied().collect(); + for (p, patch) in g.patches.iter().enumerate() { + if !present.contains(&(p as u32)) { + let at: Vec

= patch.triangles.iter().take(64).map(|&t| g.triangle_centre(t)).collect(); + return Err(MapFailure { at, why: format!("patch {p} (B-rep faces {:?}) has no region", patch.faces) }); + } + } + let mut curves_between: BTreeMap<(u32, u32), Vec> = BTreeMap::new(); + for (i, c) in g.curves.iter().enumerate() { + match c.patches.as_slice() { + [a, b] => curves_between.entry((*a, *b)).or_default().push(i as u32), + [a] => curves_between.entry((*a, *a)).or_default().push(i as u32), + _ => {} + } + } + let mut role: BTreeMap = BTreeMap::new(); + let mut border_pairs: BTreeSet<(u32, u32)> = BTreeSet::new(); + let mut vertex_labels: BTreeMap> = BTreeMap::new(); + for (q, quad) in quads.iter().enumerate() { + for &v in quad { + vertex_labels.entry(v).or_default().insert(label[q]); + } + } + for ((a, b), qs) in &edges { + let (q0, q1) = match qs.as_slice() { + [x, y] => (*x as usize, *y as usize), + _ => continue, + }; + let (la, lb) = (label[q0].min(label[q1]), label[q0].max(label[q1])); + if la == lb { + continue; + } + let Some(cs) = curves_between.get(&(la, lb)) else { + return Err(MapFailure { at: vec![xyz(*a), xyz(*b)], why: format!("regions of patches {la} and {lb} meet where the part has no curve") }); + }; + border_pairs.insert((la, lb)); + let mid = [0, 1, 2].map(|k| 0.5 * (xyz(*a)[k] + xyz(*b)[k])); + let c = nearest_curve(g, cs, mid); + for v in [*a, *b] { + role.entry(v).or_insert(Constraint::Curve(c)); + } + } + for (pair, cs) in &curves_between { + if pair.0 != pair.1 && !border_pairs.contains(pair) { + let at: Vec

= cs.iter().flat_map(|&c| g.curves[c as usize].vertices.iter().map(|&v| g.xyz(v)).collect::>()).collect(); + return Err(MapFailure { at, why: format!("the regions of patches {} and {} do not meet along their curve", pair.0, pair.1) }); + } + } + + // 4. Corners. + let mut taken: BTreeSet = BTreeSet::new(); + let junctions: Vec = vertex_labels.iter().filter(|(_, l)| l.len() >= 3).map(|(v, _)| *v).collect(); + for &v in &junctions { + let p = xyz(v); + let labels = &vertex_labels[&v]; + let best = g + .corners + .iter() + .copied() + .filter(|c| !taken.contains(c) && corner_patches(g, *c).is_superset(labels)) + .min_by(|&a, &b| dist2(g.xyz(a), p).total_cmp(&dist2(g.xyz(b), p)).then(a.cmp(&b))); + match best { + Some(c) => { + taken.insert(c); + role.insert(v, Constraint::Corner(c)); + } + None => return Err(MapFailure { at: vec![p], why: format!("regions {labels:?} meet where the part has no corner") }), + } + } + for &c in &g.corners { + if taken.contains(&c) { + continue; + } + let p = g.xyz(c); + let cp = corner_patches(g, c); + let best = role + .iter() + .filter(|(_, r)| matches!(r, Constraint::Curve(_))) + .map(|(v, _)| *v) + .filter(|v| vertex_labels.get(v).is_some_and(|l| l.is_subset(&cp))) + .min_by(|&a, &b| dist2(xyz(a), p).total_cmp(&dist2(xyz(b), p)).then(a.cmp(&b))); + match best { + Some(v) => { + role.insert(v, Constraint::Corner(c)); + } + None => return Err(MapFailure { at: vec![p], why: format!("the feature corner at {p:?} has no vertex") }), + } + } + for (q, quad) in quads.iter().enumerate() { + for &v in quad { + role.entry(v).or_insert(Constraint::Patch(label[q])); + } + } + Ok(FeatureMap { role, quad_patch: label }) +} + +/// The patches around a corner of the part (on either side of the curves that end there). +fn corner_patches(g: &Geometry, c: u32) -> BTreeSet { + g.curves.iter().filter(|cv| cv.vertices.first() == Some(&c) || cv.vertices.last() == Some(&c)).flat_map(|cv| cv.patches.iter().copied()).collect() +} + +/// The curve of `cs` nearest to `p`. +fn nearest_curve(g: &Geometry, cs: &[u32], p: P) -> u32 { + cs.iter().copied().min_by(|&a, &b| curve_distance(g, a, p).total_cmp(&curve_distance(g, b, p)).then(a.cmp(&b))).unwrap_or(0) +} + +fn curve_distance(g: &Geometry, c: u32, p: P) -> f64 { + g.curves.get(c as usize).map_or(f64::INFINITY, |cv| { + cv.vertices.windows(2).map(|w| dist2(p, closest_on_segment(p, g.xyz(w[0]), g.xyz(w[1])).0)).fold(f64::INFINITY, f64::min) + }) +} diff --git a/crates/mesh/src/octree/mod.rs b/crates/mesh/src/octree/mod.rs new file mode 100644 index 0000000..b4ecd32 --- /dev/null +++ b/crates/mesh/src/octree/mod.rs @@ -0,0 +1,801 @@ +//! Unstructured all-hexahedral meshing of any closed boundary: an adaptive grid cut into +//! conforming hexahedra, a core of grid elements inside the part, two boundary layers that carry +//! the part's faces, sharp curves and corners, and an optimisation that fits them to the +//! geometry while keeping every element valid. +//! +//! # Pipeline +//! 1. **Frame.** The grid is aligned with the part: the normals of its two largest planar +//! patches when they are orthogonal and cover more area than the coordinate axes do, the +//! coordinate axes otherwise (octree methods are orientation dependent; N. Pietroni et al., +//! *Hex-mesh generation and processing: a survey*, ACM TOG 42(2), 2022, §4.6). +//! 2. **Grid** ([`grid`]). Root cubes of the target size over the part's box, refined as a +//! 27-tree where a B-rep face curves (its normals turn by more than +//! [`OctreeOptions::curvature_angle`] inside a cell), where the part is thin (fewer than +//! [`OctreeOptions::thickness_cells`] cells across, measured along the inward normal of each +//! boundary triangle) and where a patch is narrow, then strongly balanced and closed under the +//! transition templates. +//! 3. **Conforming hexahedra** ([`templates`]): the 3-refinement node-marking templates +//! (R. Schneiders, 1996; Y. Ito, A. M. Shih, B. K. Soni, IJNME 77(13), 2009), so every grid +//! leaf becomes 1, 4, 11, 22 or 27 hexahedra and the grid is conforming by construction. +//! 4. **Core.** Elements whose eight corners are inside the part (exact ray parity on the +//! boundary triangles) and at least [`OctreeOptions::core_offset`] of their cell size from it. +//! Elements are dropped until the core's boundary is a 2-manifold and every boundary vertex can +//! move outward (its quads' normals lie in one half-space: the buffer rule of H. Tong, +//! E. Halilaj, Y. J. Zhang, *HybridOctree_Hex*, arXiv 2401.05984, 2024); its components and +//! genus must match the part's shells. +//! 5. **Feature map** ([`map`]): the part's patches, sharp curves and corners mapped onto the +//! core's boundary quads (X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based +//! hexahedral meshing*, CGF 38(5), 2019, §3.3), patches first. +//! 6. **Padding and pillowing** ([`topo`]): one layer of hexahedra on the core's boundary (so no +//! element keeps two faces on the boundary), then each patch's layer is pillowed, which inserts +//! a ring of quads along every sharp curve on both sides: a curve then runs between quads that +//! each have one edge on it, so straight and curved features are captured without degenerate +//! quads (Gao et al. 2019, §3.4, Fig. 9). +//! 7. **Fit** ([`smooth`]): the foldover-free energy of V. Garanzha et al., *Foldover-free maps +//! in 50 lines of code*, ACM TOG 40(4), 2021, plus a fitting term that draws every boundary +//! vertex onto its corner, curve or patch with a growing weight (Gao et al. 2019, §3.5); the +//! barrier keeps every element valid while the boundary moves. Each boundary vertex then goes +//! exactly onto its feature where its elements stay valid; where they would not (a blade, a +//! near-tangent contact) the mesh stays off by [`OctreeReport::boundary_distance`]. +//! 8. **Gate.** The worst elements are improved by compass search, then every element is proven +//! valid by the certified Jacobian check (Bernstein bounds, [`crate::quality`]; corner +//! Jacobians alone do not prove a trilinear hexahedron valid, J. Marschner et al., CGF 39(5), +//! 2020); elements that are not are repaired locally first. A failure in 4–8 names where the +//! grid is too coarse; the leaves there are split once more and the mesh is rebuilt, at most +//! [`OctreeOptions::attempts`] times, and a rebuild may not grow the grid more than +//! [`RETRY_GROWTH`] times (Gao et al. 2019, §3.3). +//! +//! # Why 3-refinement templates +//! The dual schemes of M. Livesu, L. Pitzalis, G. Cherchi, *Optimal dual schemes for adaptive +//! grid based hexmeshing* (ACM TOG 41(2), 2022) give a simpler singular structure, but their +//! installation depends on a paired grid whose transition vertices the reference pipeline +//! (L. Pitzalis et al., *Generalized adaptive refinement for grid-based hexahedral meshing*, ACM +//! TOG 40(6), 2021) selects with an integer linear program, and on heuristics that are retried +//! with other parameters when they fail. The node-marking templates need only the balance and +//! closure rules and are proven conforming by an exhaustive test. Their limit is the one +//! Schneiders states: they close convex refined regions only, so a leaf next to a concave corner +//! of a finer region is split, and the split spreads to the region's orthogonal hull (H. Tong, +//! Y. J. Zhang, *MCHex*, arXiv 2511.02064, 2025, meet the same limit with the same templates). +//! Curvature therefore refines only [`OctreeOptions::curvature_levels`] levels; deeper levels +//! follow thickness and narrow faces, whose refined regions are mostly convex. Schneiders' +//! modified edge rule (every edge with a marked end split in three) or the dual schemes are the +//! upgrade path. +//! +//! # Determinism and bounds +//! The same boundary and options give the same mesh: grid cells are processed in sorted order, +//! maps are only looked up, the optimiser sweeps vertices in index order. Grid cells, elements, +//! refinement depth, closure rounds, optimiser sweeps and attempts are capped. + +mod geom; +mod grid; +mod map; +mod smooth; +mod templates; +mod topo; + +#[cfg(test)] +mod tests; + +use std::cell::{Cell, RefCell}; +use std::collections::{BTreeMap, BTreeSet, HashMap, HashSet}; + +use fieldcraft_model::ElementKind; + +pub use crate::hex::BoundaryFace; +use crate::hex::face_flux; +use crate::quality::{GateOptions, Metric, Validity, element_quality, validity_gate_checks}; +use crate::size::{SizeField, checked_size}; +use crate::surface::BoundarySurface; +use crate::vec3::{P, cross, dot, norm, sub}; +use crate::{Error, Result}; + +use geom::{Aabb, Geometry}; +use grid::Tree; +use smooth::{Constraint, Optimizer}; + +/// Options and caps of [`octree_hex_mesh`]. +#[derive(Debug, Clone, PartialEq)] +pub struct OctreeOptions { + /// Refinement levels below the root cells (each one a third of the size), at most 8. + pub max_depth: u8, + /// The most grid cells and the most output elements. + pub max_elements: usize, + /// A cell is split where the normals of one B-rep face inside it turn by more than this + /// (degrees). + pub curvature_angle: f64, + /// Levels below the root that curvature may refine (see the module docs). + pub curvature_levels: u8, + /// A cell is split where the part is less than this many cells thick. + pub thickness_cells: f64, + /// Core elements keep this fraction of their cell size from the boundary. + pub core_offset: f64, + /// Optimiser sweeps over all vertices. + pub max_sweeps: usize, + /// Builds, counting the rebuilds after refining where a stage failed (1 = no retry). + pub attempts: usize, + /// Boundary edges whose faces meet at more than this angle (degrees between normals) are + /// sharp features the mesh captures with its own edges; smoother face boundaries are not. + pub sharp_angle: f64, + /// Edges whose normals turn by more than this (degrees: thin blades, near-tangent contacts) + /// are not captured; the mesh rounds them. + pub blade_angle: f64, +} + +impl Default for OctreeOptions { + fn default() -> Self { + OctreeOptions { + max_depth: 2, + max_elements: 5_000_000, + curvature_angle: 25.0, + curvature_levels: 1, + thickness_cells: 3.0, + core_offset: 0.3, + max_sweeps: 300, + attempts: 4, + sharp_angle: 30.0, + blade_angle: 160.0, + } + } +} + +/// What [`octree_hex_mesh`] did. +#[derive(Debug, Clone, PartialEq)] +pub struct OctreeReport { + /// Builds run (1 + rebuilds after a failure). + pub attempts: usize, + /// Grid leaves of the last build. + pub grid_cells: usize, + /// Core elements (from the grid) and boundary-layer elements. + pub core_elements: usize, + pub layer_elements: usize, + /// The grid's axes. + pub frame: [[f64; 3]; 3], + /// Sharp corners, curves and patches of the boundary. + pub corners: usize, + pub curves: usize, + pub patches: usize, + /// Optimiser sweeps of the last build. + pub sweeps: usize, + /// Smallest and mean corner scaled Jacobian. + pub min_scaled_jacobian: f64, + pub mean_scaled_jacobian: f64, + /// Exact volume of the elements and the boundary's enclosed volume. + pub volume: f64, + pub surface_volume: f64, + /// The largest distance of a boundary vertex from its corner, curve or patch: zero when + /// every boundary vertex lies on the geometry; where the geometry is sharper than valid + /// elements can follow (blades, slits), the mesh stays this far off. + pub boundary_distance: f64, + /// Why earlier builds were refined. + pub retries: Vec, +} + +/// An all-hex8 mesh. +#[derive(Debug, Clone, PartialEq)] +pub struct OctreeHexMesh { + pub nodes: Vec<[f64; 3]>, + /// Model node order; every element is proven valid. + pub hexes: Vec<[u32; 8]>, + pub boundary: Vec, + pub report: OctreeReport, +} + +/// A rebuild after a failure may have at most this many times the grid cells of the build +/// before… +pub const RETRY_GROWTH: usize = 4; +/// … or this many cells, whichever is more. +const RETRY_MIN_CELLS: usize = 250_000; +/// The fitting weights of the successive optimisation stages. +const FIT_STAGES: &[f64] = &[1e-3, 1e-2, 1e-1, 1.0, 10.0, 100.0]; +/// A boundary vertex goes exactly onto its feature only if its elements keep at least this +/// corner scaled Jacobian. +const MIN_SNAP_QUALITY: f64 = 0.05; +/// Elements below this corner scaled Jacobian are improved by compass search. +const IMPROVE_BELOW: f64 = 0.3; +/// Closure rounds of the grid (each splits at least one leaf; the cell cap bounds them too). +const CLOSURE_ROUNDS: usize = 1000; +/// Rounds of the core clean-up (each drops at least one element). +const CORE_ROUNDS: usize = 1000; + +/// A failure that refining at `at` may cure. +struct Retry { + at: Vec

, + why: String, +} + +enum Attempt { + Done(Box), + Retry(Retry), +} + +/// Mesh the solid bounded by `surface` with hexahedra of about `field`'s size (the root cell +/// size is the field at the centre of the part's box; finer cells follow the part's curvature, +/// thickness and narrow faces). See the module docs. +pub fn octree_hex_mesh(surface: &BoundarySurface, field: &dyn SizeField, options: &OctreeOptions) -> Result { + check_options(options)?; + let g = Geometry::new(surface, options.sharp_angle, options.blade_angle)?; + let frame = choose_frame(&g); + let mut lo = [f64::INFINITY; 3]; + let mut hi = [f64::NEG_INFINITY; 3]; + for &p in surface.points() { + let u = to_frame(&frame, p); + for k in 0..3 { + lo[k] = lo[k].min(u[k]); + hi[k] = hi[k].max(u[k]); + } + } + let centre_u = [0, 1, 2].map(|k| 0.5 * (lo[k] + hi[k])); + let h = checked_size(field, from_frame(&frame, centre_u))?; + let mut roots = [0i64; 3]; + for k in 0..3 { + let n = ((hi[k] - lo[k]) / h + 1.0).ceil(); + if !(n.is_finite() && n >= 1.0 && n <= options.max_elements as f64) { + return Err(Error::Limit(format!("a size of {h} needs {n} root cells along an axis"))); + } + roots[k] = n as i64; + } + let origin = [0, 1, 2].map(|k| centre_u[k] - 0.5 * roots[k] as f64 * h); + let ctx = Context { + g: &g, + field, + options, + frame, + origin, + h, + roots, + thickness: triangle_thickness(&g), + outside_memo: RefCell::new(HashMap::new()), + last_leaves: Cell::new(0), + too_thin: Cell::new(f64::INFINITY), + }; + let mut hints: Vec

= Vec::new(); + let mut retries: Vec = Vec::new(); + for attempt in 0..options.attempts { + let cap = if attempt == 0 { + options.max_elements + } else { + options.max_elements.min(RETRY_GROWTH.saturating_mul(ctx.last_leaves.get()).max(RETRY_MIN_CELLS)) + }; + let built = match ctx.build(&hints, cap) { + Err(Error::Limit(why)) if attempt > 0 => { + retries.push(format!("refining for another build needs more than {cap} grid cells ({why})")); + break; + } + other => other?, + }; + match built { + Attempt::Done(mut m) => { + m.report.attempts = attempt + 1; + m.report.retries = retries; + return Ok(*m); + } + Attempt::Retry(r) => { + retries.push(r.why); + hints.extend(r.at); + } + } + } + Err(Error::Failed(format!("octree hex meshing failed after {} builds: {}", retries.len(), retries.join("; ")))) +} + +fn check_options(o: &OctreeOptions) -> Result<()> { + let ok = o.max_depth <= 8 + && o.curvature_levels <= o.max_depth + && o.max_elements >= 1 + && o.curvature_angle.is_finite() + && o.curvature_angle > 0.0 + && o.thickness_cells.is_finite() + && o.thickness_cells >= 0.0 + && o.core_offset.is_finite() + && (0.0..1.0).contains(&o.core_offset) + && (1..=16).contains(&o.attempts) + && o.sharp_angle.is_finite() + && o.sharp_angle > 0.0 + && o.blade_angle.is_finite() + && o.blade_angle > o.sharp_angle + && o.blade_angle <= 180.0 + && o.max_sweeps <= 100_000; + if ok { Ok(()) } else { Err(Error::Invalid(format!("octree options out of range: {o:?}"))) } +} + +fn from_frame(f: &[[f64; 3]; 3], u: P) -> P { + [0, 1, 2].map(|k| u[0] * f[0][k] + u[1] * f[1][k] + u[2] * f[2][k]) +} + +fn to_frame(f: &[[f64; 3]; 3], p: P) -> P { + [dot(p, f[0]), dot(p, f[1]), dot(p, f[2])] +} + +/// The grid axes (rows): the coordinate axes, or the normals of the two largest orthogonal planar +/// patches when more boundary area is aligned with them. +fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] { + let tol = 1f64.to_radians().sin(); + let normal = |t: u32| g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); + let mut planar: Vec<(f64, usize, P)> = Vec::new(); + for (i, p) in g.patches.iter().enumerate() { + let Some(&t0) = p.triangles.first() else { continue }; + let n0 = normal(t0); + if p.triangles.iter().all(|&t| norm(cross(normal(t), n0)) < tol) { + planar.push((p.area, i, n0)); + } + } + planar.sort_by(|a, b| b.0.total_cmp(&a.0).then(a.1.cmp(&b.1))); + let identity = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]; + let Some(&(_, _, a0)) = planar.first() else { return identity }; + let Some(&(_, _, a1)) = planar.iter().find(|x| dot(x.2, a0).abs() < tol) else { return identity }; + let a1 = { + let v = sub(a1, a0.map(|x| x * dot(a1, a0))); + let l = norm(v); + v.map(|x| x / l) + }; + let candidate = [a0, a1, cross(a0, a1)]; + let score = |f: &[[f64; 3]; 3]| -> f64 { + (0..g.surface.triangles().len() as u32) + .filter(|&t| f.iter().any(|a| norm(cross(normal(t), *a)) < tol)) + .map(|t| { + let [a, b, c] = g.triangle(t); + 0.5 * norm(cross(sub(b, a), sub(c, a))) + }) + .sum() + }; + if score(&candidate) > score(&identity) * (1.0 + 1e-9) { candidate } else { identity } +} + +/// The distance through the part from each boundary triangle: a ray from its centre along the +/// inward normal to the next triangle (`∞` when none). +fn triangle_thickness(g: &Geometry) -> Vec { + let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0; + (0..g.surface.triangles().len() as u32) + .map(|t| { + let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); + g.ray(g.triangle_centre(t), n.map(|x| -x), far, t).map_or(f64::INFINITY, |(_, s)| s) + }) + .collect() +} + +fn centroid(points: &[P], hx: &[u32; 8]) -> P { + [0, 1, 2].map(|k| hx.iter().map(|&v| points.get(v as usize).map_or(0.0, |p| p[k])).sum::() / 8.0) +} + +fn edge_mean(points: &[P], hx: &[u32; 8]) -> f64 { + let x = hx.map(|v| points.get(v as usize).copied().unwrap_or([0.0; 3])); + let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)]; + let m = edges.iter().map(|&(a, b)| norm(sub(x[a], x[b]))).sum::() / 12.0; + if m > 0.0 { m } else { 1.0 } +} + +/// The part being meshed and the state carried from one build to the next. +struct Context<'a> { + g: &'a Geometry<'a>, + field: &'a dyn SizeField, + options: &'a OctreeOptions, + frame: [[f64; 3]; 3], + origin: P, + h: f64, + roots: [i64; 3], + thickness: Vec, + outside_memo: RefCell>, + /// Grid leaves of the last build. + last_leaves: Cell, + /// The thinnest wall the deepest cells cannot cross in `thickness_cells` (∞ when none). + too_thin: Cell, +} + +/// The mesh while its boundary layers are built: points, elements, their target sizes and how +/// each point may move. +struct Work { + points: Vec

, + hexes: Vec<[u32; 8]>, + size: Vec, + role: Vec, + core_count: usize, +} + +impl Context<'_> { + fn edge(&self, level: u8) -> f64 { + self.h / 3f64.powi(i32::from(level)) + } + + /// The world box (inflated by `pad` cell edges) and the centre of a cell. + fn cell_box(&self, (level, i): grid::Key, pad: f64) -> (Aabb, P) { + let e = self.edge(level); + let mut b = Aabb::EMPTY; + for c in 0..8 { + b.add(from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + ((c >> k) & 1) as f64) * e))); + } + let centre = from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + 0.5) * e)); + let p = pad * e; + (Aabb { lo: b.lo.map(|x| x - p), hi: b.hi.map(|x| x + p) }, centre) + } + + fn needs_split(&self, key: grid::Key) -> Result { + let e = self.edge(key.0); + let (bx, centre) = self.cell_box(key, 0.05); + if e > checked_size(self.field, centre)? * (1.0 + 1e-9) { + return Ok(true); + } + let reach = 0.5 * 3f64.sqrt() * e * 1.05; + let cos = self.options.curvature_angle.to_radians().cos(); + let mut first: BTreeMap = BTreeMap::new(); + for t in self.g.overlapping(&bx) { + let [a, b, c] = self.g.triangle(t); + if geom::dist2(centre, geom::closest_on_triangle(centre, a, b, c)) > reach * reach { + continue; + } + let n = self.g.normals.get(t as usize).copied().unwrap_or([0.0; 3]); + let face = self.g.surface.tri_face().get(t as usize).copied().unwrap_or(0); + let f = *first.entry(face).or_insert(n); + if key.0 < self.options.curvature_levels && dot(f, n) < cos { + return Ok(true); + } + if let Some(&th) = self.thickness.get(t as usize) + && th < self.options.thickness_cells * e + { + if key.0 + 1 >= self.options.max_depth && th < self.options.thickness_cells * e / 3.0 { + self.too_thin.set(self.too_thin.get().min(th)); + } + return Ok(true); + } + // Narrow patches (a hole's wall, a thin rim) need cells across them. + if self.g.tri_patch.get(t as usize).and_then(|&p| self.g.patches.get(p as usize)).is_some_and(|pt| pt.width < 2.0 * e) { + return Ok(true); + } + } + Ok(false) + } + + /// Whether a cell lies wholly outside the part (its centre outside and farther from the + /// boundary than its half diagonal): no core element can come from it, so its transition + /// pattern does not matter and it is not meshed. Memoised. + fn outside(&self, key: grid::Key) -> bool { + if let Some(&v) = self.outside_memo.borrow().get(&key) { + return v; + } + let (_, centre) = self.cell_box(key, 0.0); + let half = 0.5 * 3f64.sqrt() * self.edge(key.0); + let far = self.g.nearest(centre, |_| true).is_some_and(|x| x.2 > half * half * (1.0 + 1e-9)); + let v = far && self.g.inside(centre) == Some(false); + self.outside_memo.borrow_mut().insert(key, v); + v + } + + /// The level-`level` cells within one cell of a hint point. + fn hint_cells(&self, hints: &[P], level: u8) -> HashSet { + let e = self.edge(level); + let mut out = HashSet::new(); + for &p in hints { + let u = to_frame(&self.frame, p); + let i = [0, 1, 2].map(|k| ((u[k] - self.origin[k]) / e).floor()); + if !i.iter().all(|x| x.is_finite() && x.abs() < 1e15) { + continue; + } + for d in 0..27i64 { + out.insert((level, [i[0] as i64 + d % 3 - 1, i[1] as i64 + (d / 3) % 3 - 1, i[2] as i64 + d / 9 - 1])); + } + } + out + } + + /// The refined, balanced and closed grid; leaves near `hints` (where an earlier build + /// failed) are split once more. + fn tree(&self, hints: &[P], max_leaves: usize) -> Result { + let mut tree = Tree::new(self.roots, self.options.max_depth, max_leaves)?; + for level in 0..self.options.max_depth { + for key in tree.sorted_leaves() { + if key.0 == level && self.needs_split(key)? { + tree.split(key)?; + } + } + } + if !hints.is_empty() { + let near: Vec> = (0..self.options.max_depth).map(|l| self.hint_cells(hints, l)).collect(); + let split: Vec = + tree.sorted_leaves().into_iter().filter(|k| near.get(usize::from(k.0)).is_some_and(|s| s.contains(k))).collect(); + for key in split { + tree.split(key)?; + } + } + tree.close(CLOSURE_ROUNDS, &|k| self.outside(k))?; + Ok(tree) + } + + fn build(&self, hints: &[P], max_leaves: usize) -> Result { + let g = self.g; + let tree = self.tree(hints, max_leaves)?; + self.last_leaves.set(tree.leaf_count()); + let lattice = tree.mesh(&|k| self.outside(k))?; + if lattice.hexes.len() > self.options.max_elements { + return Err(Error::Limit(format!("the grid has {} elements; the limit is {}", lattice.hexes.len(), self.options.max_elements))); + } + let unit = self.h / tree.span(0) as f64; + let points: Vec

= lattice.points.iter().map(|q| from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + q[k] as f64 * unit))).collect(); + + // Core. + let dist: Vec = points.iter().map(|&p| g.nearest(p, |_| true).map_or(0.0, |x| x.2.sqrt())).collect(); + let keep = self.core(&lattice, &points, &dist); + let core: Vec = (0..lattice.hexes.len()).filter(|&i| keep[i]).collect(); + if core.is_empty() { + let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect(); + return Ok(Attempt::Retry(Retry { at, why: "no grid element lies inside the part".into() })); + } + let mut w = Work { + hexes: core.iter().map(|&i| lattice.hexes[i]).collect(), + size: core.iter().map(|&i| self.edge(lattice.level[i])).collect(), + role: Vec::new(), + points, + core_count: core.len(), + }; + let qb = topo::boundary(&w.hexes); + let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect(); + if let Some(r) = check_topology(g, &q) { + let thin = self.too_thin.get(); + if thin.is_finite() { + return Err(Error::Failed(format!( + "{}: the part has walls {thin:.3} thick, thinner than {} cells of the deepest level ({:.3}); use a smaller size, more levels, or a structured (swept) hex mesh", + r.why, + self.options.thickness_cells, + self.edge(self.options.max_depth) + ))); + } + return Ok(Attempt::Retry(r)); + } + + // Feature map, padding and pillows. + let fmap = { + let pts = &w.points; + match map::map_features(g, &q, &|v| pts.get(v as usize).copied().unwrap_or([0.0; 3])) { + Ok(m) => m, + Err(f) => return Ok(Attempt::Retry(Retry { at: f.at, why: f.why })), + } + }; + self.layers(&mut w, &qb, &fmap, &dist)?; + if w.hexes.len() > self.options.max_elements { + return Err(Error::Limit(format!("the mesh needs {} elements; the limit is {}", w.hexes.len(), self.options.max_elements))); + } + + // Fit, then improve and repair. + let mut opt = Optimizer::new(g, std::mem::take(&mut w.points), &w.hexes, &w.size, w.role.clone()); + let mut ran = 0; + let sweeps = (self.options.max_sweeps / FIT_STAGES.len()).max(1); + for &weight in FIT_STAGES { + opt.fit = weight; + ran = opt.run(sweeps, 1e-4); + } + let (_, boundary_distance) = opt.snap(MIN_SNAP_QUALITY); + let inverted: Vec

= w + .hexes + .iter() + .filter(|hx| smooth::min_scaled_jacobian(&hx.map(|v| opt.points[v as usize])) <= 0.0) + .map(|hx| centroid(&opt.points, hx)) + .collect(); + if !inverted.is_empty() { + let why = format!("{} elements stay inverted after untangling", inverted.len()); + return Ok(Attempt::Retry(Retry { at: inverted, why })); + } + opt.improve_worst(IMPROVE_BELOW, 4); + let uncertified = opt.repair_certified(6); + if !uncertified.is_empty() { + let at: Vec

= uncertified.iter().filter_map(|&e| w.hexes.get(e as usize)).map(|hx| centroid(&opt.points, hx)).collect(); + let why = format!("{} elements fail the certified validity check", at.len()); + return Ok(Attempt::Retry(Retry { at, why })); + } + w.points = opt.points; + self.finish(&w, &tree, ran, boundary_distance) + } + + /// The core: elements wholly inside the part, cleaned to a 2-manifold boundary whose + /// vertices can all move outward. + fn core(&self, lattice: &grid::LatticeMesh, points: &[P], dist: &[f64]) -> Vec { + let inside: Vec = points.iter().map(|&p| self.g.inside(p) == Some(true)).collect(); + let mut keep: Vec = lattice + .hexes + .iter() + .zip(&lattice.level) + .map(|(hx, &l)| { + let floor = self.options.core_offset * self.edge(l); + hx.iter().all(|&v| inside.get(v as usize) == Some(&true) && dist.get(v as usize).is_some_and(|&d| d >= floor)) + }) + .collect(); + for _ in 0..CORE_ROUNDS { + let mut drop = topo::nonmanifold(&lattice.hexes, &keep); + if drop.is_empty() { + drop = topo::hidden_vertices(&lattice.hexes, &keep, &|v| points.get(v as usize).copied().unwrap_or([0.0; 3]), 0.1); + } + if drop.is_empty() { + break; + } + for d in drop { + if let Some(k) = keep.get_mut(d as usize) { + *k = false; + } + } + } + keep + } + + /// The boundary layers: padding on the core's boundary quads `qb`, whose outer copies take the + /// map's roles and start halfway to the boundary along the core's normals, then a pillow + /// under each patch, whose copies start a quarter of the way into their elements. Every + /// element starts valid or nearly so; the fit untangles the rest. + fn layers(&self, w: &mut Work, qb: &[(topo::FaceRef, [u32; 4])], fmap: &map::FeatureMap, dist: &[f64]) -> Result<()> { + let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect(); + let core_count = w.core_count; + let copy = topo::pad(&mut w.points, &mut w.hexes, &q)?; + w.role = vec![Constraint::Free; w.points.len()]; + let mut normal: BTreeMap = BTreeMap::new(); + for f in &q { + let p = f.map(|v| w.points[v as usize]); + let n = cross(sub(p[2], p[0]), sub(p[3], p[1])); + let l = norm(n); + if l > 0.0 { + for &v in f { + let e = normal.entry(v).or_insert([0.0; 3]); + for k in 0..3 { + e[k] += n[k] / l; + } + } + } + } + for (&v, &c) in © { + if let Some(r) = fmap.role.get(&v) { + w.role[c as usize] = *r; + } + let n = normal.get(&v).copied().unwrap_or([0.0; 3]); + let l = norm(n); + let d = 0.5 * dist.get(v as usize).copied().unwrap_or(0.0); + if l > 0.0 { + w.points[c as usize] = [0, 1, 2].map(|k| w.points[v as usize][k] + d * n[k] / l); + } + } + for ((owner, _), _) in qb { + w.size.push(0.5 * w.size.get(*owner as usize).copied().unwrap_or(self.h)); + } + let patches: BTreeSet = fmap.quad_patch.iter().copied().collect(); + for p in patches { + let inside: Vec = + (0..w.hexes.len()).map(|i| i >= core_count && i < core_count + q.len() && fmap.quad_patch.get(i - core_count) == Some(&p)).collect(); + let before = w.hexes.len(); + let centroids: Vec

= w.hexes.iter().map(|hx| centroid(&w.points, hx)).collect(); + let copies = topo::pillow(&mut w.points, &mut w.hexes, &inside)?; + w.role.resize(w.points.len(), Constraint::Free); + for pc in &copies { + let orig = w.points[pc.original as usize]; + let n = pc.elements.len().max(1) as f64; + let mut c = [0.0; 3]; + for &e in &pc.elements { + let ce = centroids.get(e as usize).copied().unwrap_or(orig); + for k in 0..3 { + c[k] += ce[k] / n; + } + } + w.points[pc.copy as usize] = [0, 1, 2].map(|k| orig[k] + 0.25 * (c[k] - orig[k])); + w.role[pc.copy as usize] = match w.role.get(pc.original as usize).copied().unwrap_or(Constraint::Free) { + Constraint::Free => Constraint::Free, + _ => Constraint::Patch(p), + }; + } + for hx in &w.hexes[before..] { + w.size.push(edge_mean(&w.points, hx)); + } + } + Ok(()) + } + + /// The gate and the output: every element proven valid, nodes renumbered, boundary faces with + /// their B-rep faces. + fn finish(&self, w: &Work, tree: &Tree, sweeps: usize, boundary_distance: f64) -> Result { + let g = self.g; + // One certified pass: every element's quality (the same metric as mesh.quality and + // mesh.hex.generate report) and its Jacobian check for the gate. + let mut bad: Vec

= Vec::new(); + let mut checks = Vec::with_capacity(w.hexes.len()); + let (mut sj_min, mut sj_sum) = (f64::INFINITY, 0.0); + for (e, hx) in w.hexes.iter().enumerate() { + let q = element_quality(ElementKind::Hex8, &hx.map(|v| w.points[v as usize]))?; + let check = q.jacobian.ok_or_else(|| Error::Failed("an element was not certified".into()))?; + if check.validity != Validity::Valid { + bad.push(centroid(&w.points, hx)); + } + let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0); + sj_min = sj_min.min(sj); + sj_sum += sj; + checks.push((e as u32 + 1, check)); + } + if !bad.is_empty() { + let why = format!("{} elements fail the certified validity check", bad.len()); + return Ok(Attempt::Retry(Retry { at: bad, why })); + } + // The volume the mesh's own boundary encloses equals the elements' total exactly when they + // tile it without overlapping. + let faces = topo::boundary(&w.hexes); + let enclosed: f64 = faces.iter().map(|(_, quad)| face_flux(&quad.map(|v| w.points[v as usize]))).sum(); + let gate = + validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })? + .into_result()?; + let volume = gate.volume; + let mut new_id = vec![u32::MAX; w.points.len()]; + let mut nodes = Vec::new(); + let mut hexes = Vec::with_capacity(w.hexes.len()); + for hx in &w.hexes { + let mut o = [0u32; 8]; + for (k, &v) in hx.iter().enumerate() { + if new_id[v as usize] == u32::MAX { + new_id[v as usize] = u32::try_from(nodes.len()).map_err(|_| Error::Limit("too many nodes".into()))?; + nodes.push(w.points[v as usize]); + } + o[k] = new_id[v as usize]; + } + hexes.push(o); + } + let mut boundary = Vec::new(); + for ((e, f), quad) in faces { + let patch = quad.iter().find_map(|&v| match w.role.get(v as usize) { + Some(Constraint::Patch(p)) => Some(*p), + _ => None, + }); + let c = [0, 1, 2].map(|k| quad.iter().map(|&v| w.points[v as usize][k]).sum::() / 4.0); + let brep_face = match patch.and_then(|p| g.patches.get(p as usize)) { + Some(pt) if pt.faces.len() == 1 => pt.faces[0], + _ => { + let keep = |t: u32| patch.is_none_or(|p| g.tri_patch.get(t as usize) == Some(&p)); + g.nearest(c, keep).and_then(|(t, _, _)| g.surface.tri_face().get(t as usize).copied()).unwrap_or(0) + } + }; + boundary.push(BoundaryFace { element: e, face: f + 1, brep_face }); + } + let report = OctreeReport { + attempts: 0, + grid_cells: tree.leaf_count(), + core_elements: w.core_count, + layer_elements: w.hexes.len() - w.core_count, + frame: self.frame, + corners: g.corners.len(), + curves: g.curves.len(), + patches: g.patches.len(), + sweeps, + min_scaled_jacobian: sj_min, + mean_scaled_jacobian: sj_sum / w.hexes.len().max(1) as f64, + volume, + surface_volume: g.surface.volume(), + boundary_distance, + retries: Vec::new(), + }; + Ok(Attempt::Done(Box::new(OctreeHexMesh { nodes, hexes, boundary, report }))) + } +} + +/// The core boundary's components and genus against the part's shells; on a mismatch, a retry +/// that refines near the whole boundary. +fn check_topology(g: &Geometry, quads: &[[u32; 4]]) -> Option { + let edges = topo::quad_edges(quads); + let mut comp: Vec = (0..quads.len()).collect(); + for qs in edges.values() { + for w in qs.windows(2) { + let (a, b) = (topo::find(&mut comp, w[0] as usize), topo::find(&mut comp, w[1] as usize)); + comp[a] = b; + } + } + // Per component: vertices, edges, faces. + let mut stats: BTreeMap, i64, i64)> = BTreeMap::new(); + for (i, q) in quads.iter().enumerate() { + let s = stats.entry(topo::find(&mut comp, i)).or_default(); + s.2 += 1; + s.0.extend(q.iter().copied()); + } + for qs in edges.values() { + if let Some(&q0) = qs.first() + && let Some(s) = stats.get_mut(&topo::find(&mut comp, q0 as usize)) + { + s.1 += 1; + } + } + let mut got: Vec = stats.values().map(|s| (2 - (s.0.len() as i64 - s.1 + s.2)) / 2).collect(); + let mut want: Vec = g.surface.shells().iter().map(|s| i64::from(s.genus)).collect(); + got.sort_unstable(); + want.sort_unstable(); + if got == want { + return None; + } + let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect(); + Some(Retry { at, why: format!("the core's boundary has components of genus {got:?}, the part {want:?}") }) +} diff --git a/crates/mesh/src/octree/smooth.rs b/crates/mesh/src/octree/smooth.rs new file mode 100644 index 0000000..b271a5e --- /dev/null +++ b/crates/mesh/src/octree/smooth.rs @@ -0,0 +1,641 @@ +//! Untangling, fitting and shape optimisation of the hexahedral mesh. +//! +//! # Energy +//! The foldover-free mapping energy of V. Garanzha, I. Kaporin, L. Kudryavtseva, F. Protais, +//! N. Ray, D. Sokolov, *Foldover-free maps in 50 lines of code*, ACM TOG 40(4):102, 2021 (paper +//! only; the authors' code is AGPL and is not read). Each hexahedron is measured at nine points, +//! its eight corners (the trapezoidal rule the paper uses for quads) and its centre, where a +//! twisted element folds first: `J` is the sample's three edge directions against the ideal cube +//! of the element's target edge `s`, `D = det J`, and the sample costs +//! +//! `f_ε + λ g_ε`, `f_ε = |J|³_F / χ(D, ε)`, `g_ε = (D² + 1) / χ(D, ε)`, +//! `χ(D, ε) = (D + √(ε² + D²)) / 2`, +//! +//! weighted by `s³`. `f` is the shape (conformal) term, `g` the volume term. The paper's shape +//! term is `|J|²/D^{2/3}`; its 3/2 power `|J|³/D`, the classical polyconvex distortion in three +//! dimensions (J. M. Ball, 1977), is minimal at the same ideal shapes and needs only square +//! roots, so the optimiser gives the same bits on every platform (no `pow` of the system maths +//! library). `χ` keeps both finite for inverted samples, and `ε` follows the paper's rule (7), +//! `ε = √(10⁻¹² + 0.04 min(0, min D)²)`, so it vanishes once every sample is positive and the +//! energy becomes a barrier against inversion. +//! +//! A vertex on a corner, curve or patch adds `fit · 8 s̄ · dist²` (`s̄` the mean target edge of its +//! elements) to its distance from that feature: with `fit` raised stage by stage the boundary is +//! drawn onto the geometry while the barrier keeps every element valid (Gao, Shen, Panozzo 2019, +//! §3.5). With `fit = 0` constrained vertices are held on their features instead (each step is +//! restricted to the tangent and projected back). +//! +//! # Minimisation +//! Nonlinear Gauss–Seidel: each vertex in turn takes a Newton step on its own energy with the +//! paper's positive approximation of the Hessian (the second derivatives of `χ` and `det` are +//! dropped, §2.3, App. A), accepted only if the energy decreased (halving line search). The paper +//! minimises the same energy globally with L-BFGS; the vertex sweep is the local form of the same +//! descent. Then a compass search raises the worst corner scaled Jacobians, and another repairs +//! elements whose trilinear map the certified check does not prove positive. Sweeps, line-search +//! steps, compass steps and rounds are capped. + +use fieldcraft_model::ElementKind; + +use crate::quality::{Validity, element_validity}; +use crate::vec3::{P, cross, dot, norm, sub}; + +use super::geom::dist2; + +use super::geom::Geometry; + +/// How a vertex may move. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Constraint { + Free, + /// On a corner of the part (the boundary surface's vertex). + Corner(u32), + /// On a feature curve. + Curve(u32), + /// On a patch of the boundary. + Patch(u32), +} + +/// Model-order corner `c` of a hex8: its three neighbours along x, y, z and the signs that make +/// the reference corner frame right-handed. +const CORNER_FRAMES: [([usize; 3], [f64; 3]); 8] = [ + ([1, 3, 4], [1.0, 1.0, 1.0]), + ([0, 2, 5], [-1.0, 1.0, 1.0]), + ([3, 1, 6], [-1.0, -1.0, 1.0]), + ([2, 0, 7], [1.0, -1.0, 1.0]), + ([5, 7, 0], [1.0, 1.0, -1.0]), + ([4, 6, 1], [-1.0, 1.0, -1.0]), + ([7, 5, 2], [-1.0, -1.0, -1.0]), + ([6, 4, 3], [1.0, -1.0, -1.0]), +]; + +/// Weight of the volume term. +const LAMBDA: f64 = 0.1; + +fn chi(d: f64, eps: f64) -> f64 { + 0.5 * (d + (eps * eps + d * d).sqrt()) +} + +fn chi_d(d: f64, eps: f64) -> f64 { + 0.5 * (1.0 + d / (eps * eps + d * d).sqrt()) +} + +/// The corner Jacobian (columns) of corner `c` of `x` against an ideal cube of edge `s`. +fn corner_jacobian(x: &[P; 8], c: usize, s: f64) -> [P; 3] { + let (n, sg) = CORNER_FRAMES[c]; + [0, 1, 2].map(|a| sub(x[n[a]], x[c]).map(|v| v * sg[a] / s)) +} + +/// The sample points of the energy: the eight corners and the centre. Column `a` of sample `s`'s +/// Jacobian is `Σ_k W[s][a][k] x_k` (before dividing by the target edge). The centre keeps +/// twisted elements (positive corners, folded inside) out of reach, as the nine-point measure of +/// ParaView and of Zhang et al.'s meshers does. +const SAMPLES: usize = 9; + +fn sample_weights(s: usize) -> [[f64; 8]; 3] { + let mut w = [[0.0; 8]; 3]; + if s < 8 { + let (n, sg) = CORNER_FRAMES[s]; + for a in 0..3 { + w[a][n[a]] += sg[a]; + w[a][s] -= sg[a]; + } + } else { + w = [ + [-0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25, -0.25], + [-0.25, -0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25], + [-0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25], + ]; + } + w +} + +/// The Jacobian (columns) of sample `smp` of `x` against an ideal cube of edge `s`. +fn sample_jacobian(x: &[P; 8], smp: usize, s: f64) -> [P; 3] { + let w = sample_weights(smp); + [0, 1, 2].map(|a| { + let mut c = [0.0; 3]; + for (k, p) in x.iter().enumerate() { + for r in 0..3 { + c[r] += w[a][k] * p[r]; + } + } + c.map(|v| v / s) + }) +} + +fn det3(j: &[P; 3]) -> f64 { + dot(j[0], cross(j[1], j[2])) +} + +/// The cofactor columns: `∂ det / ∂ column a`. +fn cofactors(j: &[P; 3]) -> [P; 3] { + [cross(j[1], j[2]), cross(j[2], j[0]), cross(j[0], j[1])] +} + +/// The energy of one corner. +fn corner_energy(j: &[P; 3], eps: f64) -> f64 { + let d = det3(j); + let c = chi(d, eps); + let fro: f64 = j.iter().map(|v| dot(*v, *v)).sum(); + fro * fro.sqrt() / c + LAMBDA * (d * d + 1.0) / c +} + +/// The scaled Jacobian of each corner, smallest (`min over corners of det / product of edge +/// lengths`), the usual hex quality measure. +pub(crate) fn min_scaled_jacobian(x: &[P; 8]) -> f64 { + (0..8) + .map(|c| { + let j = corner_jacobian(x, c, 1.0); + let l = norm(j[0]) * norm(j[1]) * norm(j[2]); + if l > 0.0 { det3(&j) / l } else { -1.0 } + }) + .fold(f64::INFINITY, f64::min) +} + +/// The optimiser's state. +pub(crate) struct Optimizer<'a> { + pub geom: &'a Geometry<'a>, + pub points: Vec

, + pub hexes: &'a [[u32; 8]], + /// Target edge per element. + pub size: &'a [f64], + pub constraint: Vec, + /// Elements around each vertex with the vertex's local position. + around: Vec>, + pub eps: f64, + /// Weight of the fitting term; 0 holds boundary vertices on their features exactly. + pub fit: f64, +} + +impl<'a> Optimizer<'a> { + pub fn new(geom: &'a Geometry<'a>, points: Vec

, hexes: &'a [[u32; 8]], size: &'a [f64], constraint: Vec) -> Optimizer<'a> { + let mut around: Vec> = vec![Vec::new(); points.len()]; + for (e, h) in hexes.iter().enumerate() { + for (k, &v) in h.iter().enumerate() { + if let Some(a) = around.get_mut(v as usize) { + a.push((e as u32, k as u8)); + } + } + } + Optimizer { geom, points, hexes, size, constraint, around, eps: 1.0, fit: 0.0 } + } + + fn element(&self, e: u32) -> ([P; 8], f64) { + let h = self.hexes.get(e as usize).copied().unwrap_or([0; 8]); + (h.map(|v| self.points.get(v as usize).copied().unwrap_or([0.0; 3])), self.size.get(e as usize).copied().unwrap_or(1.0)) + } + + /// The smallest normalised corner determinant over the mesh, and the total energy. + pub fn measure(&self) -> (f64, f64) { + let (mut m, mut total) = (f64::INFINITY, 0.0); + for e in 0..self.hexes.len() as u32 { + let (x, s) = self.element(e); + for c in 0..SAMPLES { + let j = sample_jacobian(&x, c, s); + m = m.min(det3(&j)); + total += s * s * s * corner_energy(&j, self.eps); + } + } + (m, total) + } + + /// Where vertex `v`'s feature is closest to `p` (`None` for a free vertex). + pub fn target(&self, v: u32, p: P) -> Option

{ + match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) { + Constraint::Free => None, + Constraint::Corner(c) => Some(self.geom.xyz(c)), + Constraint::Curve(c) => self.geom.nearest_on_curve(c, p).map(|x| x.0), + Constraint::Patch(pt) => self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)).map(|x| x.1), + } + } + + /// The fitting weight of vertex `v`: the fit factor times its elements' size (so a distance + /// of one element weighs like the shape of its corners), ten times more for a corner. + fn fit_weight(&self, v: u32) -> f64 { + let a = self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()); + if a.is_empty() || self.fit == 0.0 { + return 0.0; + } + let s = a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(1.0)).sum::() / a.len() as f64; + let corner = if matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_))) { 10.0 } else { 1.0 }; + 8.0 * self.fit * corner * s + } + + /// The energy of the corners that involve vertex `v`, with its fitting term. + fn local_energy(&self, v: u32) -> f64 { + let mut total = 0.0; + if self.fit > 0.0 { + let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + if let Some(t) = self.target(v, p) { + total += self.fit_weight(v) * dist2(p, t); + } + } + for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) { + let (x, s) = self.element(e); + let (n, _) = CORNER_FRAMES[k as usize]; + for c in [k as usize, n[0], n[1], n[2], 8] { + total += s * s * s * corner_energy(&sample_jacobian(&x, c, s), self.eps); + } + } + total + } + + /// The gradient and approximate Hessian of the local energy at vertex `v`. + fn derivatives(&self, v: u32) -> (P, [[f64; 3]; 3]) { + let mut g = [0.0; 3]; + let mut h = [[0.0; 3]; 3]; + for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) { + let (x, s) = self.element(e); + let k = k as usize; + let w = s * s * s; + let (nk, _) = CORNER_FRAMES[k]; + for c in [k, nk[0], nk[1], nk[2], 8] { + // ∂ column a / ∂ x_v = coef[a] I. + let wt = sample_weights(c); + let coef: [f64; 3] = [0, 1, 2].map(|a| wt[a][k] / s); + let j = sample_jacobian(&x, c, s); + let d = det3(&j); + let ch = chi(d, self.eps); + let chd = chi_d(d, self.eps); + let fro: f64 = j.iter().map(|q| dot(*q, *q)).sum(); + let cof = cofactors(&j); + let root = fro.sqrt(); + // dE/dJ = A J + B cof(J). + let a_coef = 3.0 * root / ch; + let f_d = -fro * root * chd / (ch * ch); + let g_d = 2.0 * d / ch - (d * d + 1.0) * chd / (ch * ch); + let b_coef = f_d + LAMBDA * g_d; + let mut dd = [0.0; 3]; + let mut sum_c2 = 0.0; + for a in 0..3 { + if coef[a] == 0.0 { + continue; + } + sum_c2 += coef[a] * coef[a]; + for r in 0..3 { + g[r] += w * coef[a] * (a_coef * j[a][r] + b_coef * cof[a][r]); + dd[r] += coef[a] * cof[a][r]; + } + } + // Positive part of the second derivative along D, χ'' and det'' dropped. + let f_dd = 2.0 * fro * root * chd * chd / (ch * ch * ch); + let g_dd = 2.0 / ch - 4.0 * d * chd / (ch * ch) + 2.0 * (d * d + 1.0) * chd * chd / (ch * ch * ch); + let k_dd = (f_dd + LAMBDA * g_dd).max(0.0); + for r in 0..3 { + h[r][r] += w * a_coef * sum_c2; + for q in 0..3 { + h[r][q] += w * k_dd * dd[r] * dd[q]; + } + } + } + } + if self.fit > 0.0 { + let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + if let Some(t) = self.target(v, p) { + let w = self.fit_weight(v); + for r in 0..3 { + g[r] += 2.0 * w * (p[r] - t[r]); + h[r][r] += 2.0 * w; + } + } + } + (g, h) + } + + /// Project a point onto vertex `v`'s constraint. + fn project(&self, v: u32, p: P) -> P { + self.target(v, p).unwrap_or(p) + } + + /// The direction a constrained vertex may move along: the step restricted to its tangent. + fn restrict(&self, v: u32, d: P) -> P { + let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) { + Constraint::Free => d, + Constraint::Corner(_) => [0.0; 3], + Constraint::Curve(c) => { + let Some((_, i, _)) = self.geom.nearest_on_curve(c, p) else { return [0.0; 3] }; + let Some(cv) = self.geom.curves.get(c as usize) else { return [0.0; 3] }; + let (a, b) = (cv.vertices.get(i).copied().unwrap_or(0), cv.vertices.get(i + 1).copied().unwrap_or(0)); + let t = sub(self.geom.xyz(b), self.geom.xyz(a)); + let l2 = dot(t, t); + if l2 > 0.0 { t.map(|x| x * dot(d, t) / l2) } else { [0.0; 3] } + } + Constraint::Patch(pt) => { + let Some((t, _, _)) = self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)) else { return [0.0; 3] }; + let n = self.geom.normals.get(t as usize).copied().unwrap_or([0.0; 3]); + let dn = dot(d, n); + [0, 1, 2].map(|k| d[k] - dn * n[k]) + } + } + } + + /// One Newton step of vertex `v`; whether it moved. + fn relax(&mut self, v: u32) -> bool { + let soft = self.fit > 0.0; + if (!soft && matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_)))) + || self.around.get(v as usize).is_none_or(|a| a.is_empty()) + { + return false; + } + let (g, h) = self.derivatives(v); + let Some(step) = solve3(h, g.map(|x| -x)) else { return false }; + let step = if soft { step } else { self.restrict(v, step) }; + if !step.iter().all(|x| x.is_finite()) || norm(step) == 0.0 { + return false; + } + let old = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + let e0 = self.local_energy(v); + let mut tau = 1.0; + for _ in 0..12 { + let moved = [0, 1, 2].map(|k| old[k] + tau * step[k]); + let trial = if soft { moved } else { self.project(v, moved) }; + if let Some(p) = self.points.get_mut(v as usize) { + *p = trial; + } + let e1 = self.local_energy(v); + if e1 < e0 && e1.is_finite() { + return true; + } + tau *= 0.5; + } + if let Some(p) = self.points.get_mut(v as usize) { + *p = old; + } + false + } + + /// Set the regularisation from the smallest corner determinant (Garanzha et al. 2021: ε + /// shrinks as the mesh untangles); that determinant and the energy at the new ε. + fn update_eps(&mut self) -> (f64, f64) { + let (min_d, energy) = self.measure(); + let old = self.eps; + self.eps = (1e-12 + 0.04 * min_d.min(0.0).powi(2)).sqrt(); + (min_d, if self.eps == old { energy } else { self.measure().1 }) + } + + /// Sweep all vertices until every corner is positive and the energy settles, or `max_sweeps`; + /// the sweeps run. + pub fn run(&mut self, max_sweeps: usize, tolerance: f64) -> usize { + let n = self.points.len() as u32; + let (_, mut energy) = self.update_eps(); + let mut sweeps = 0; + for sweep in 0..max_sweeps { + let mut moved = 0usize; + for v in 0..n { + if self.relax(v) { + moved += 1; + } + } + let old_eps = self.eps; + let (min_d, e) = self.update_eps(); + sweeps = sweep + 1; + let settled = old_eps == self.eps && (energy - e).abs() <= tolerance * energy.abs(); + energy = e; + if moved == 0 || (min_d > 0.0 && settled) { + break; + } + } + sweeps + } +} + +impl Optimizer<'_> { + /// The smallest corner scaled Jacobian of the elements around vertex `v`. + fn local_quality(&self, v: u32) -> f64 { + self.around + .get(v as usize) + .map_or(f64::INFINITY, |a| a.iter().map(|&(e, _)| min_scaled_jacobian(&self.element(e).0)).fold(f64::INFINITY, f64::min)) + } + + /// Move every constrained vertex exactly onto its feature where its elements keep a corner + /// scaled Jacobian of at least `floor`; the others stay where the fit left them. Returns + /// how many stayed off and the largest distance of a constrained vertex from its feature. + pub fn snap(&mut self, floor: f64) -> (usize, f64) { + let (mut off, mut worst) = (0usize, 0.0f64); + for v in 0..self.points.len() as u32 { + let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + let Some(t) = self.target(v, p) else { continue }; + if p == t { + continue; + } + if let Some(x) = self.points.get_mut(v as usize) { + *x = t; + } + if self.local_quality(v) < floor { + if let Some(x) = self.points.get_mut(v as usize) { + *x = p; + } + off += 1; + worst = worst.max(dist2(p, t).sqrt()); + } + } + (off, worst) + } + + /// Raise the worst elements: every vertex of an element whose corner scaled Jacobian is + /// below `threshold` moves, within its constraint, by a compass search that maximises the + /// smallest scaled Jacobian around it (steps shrink by halves; at most `passes` rounds over + /// the worst elements, 12 step sizes per vertex). + pub fn improve_worst(&mut self, threshold: f64, passes: usize) { + for _ in 0..passes { + let mut verts: Vec = Vec::new(); + for (e, h) in self.hexes.iter().enumerate() { + if min_scaled_jacobian(&self.element(e as u32).0) < threshold { + verts.extend(h.iter().copied()); + } + } + verts.sort_unstable(); + verts.dedup(); + if verts.is_empty() { + break; + } + let mut moved = false; + for &v in &verts { + if matches!(self.constraint.get(v as usize), Some(Constraint::Free)) { + moved |= self.compass(v, &|o: &Self, v| o.local_quality(v)); + } + } + if !moved { + break; + } + } + } + + /// Compass search on vertex `v` (within its constraint) for the largest `quality`; steps + /// start at 0.15 of its largest element size and halve when no direction improves, 12 times + /// at most. Whether it moved. + fn compass(&mut self, v: u32, quality: &dyn Fn(&Self, u32) -> f64) -> bool { + let p0 = self.points.get(v as usize).copied().unwrap_or([0.0; 3]); + let reach = self + .around + .get(v as usize) + .map_or(0.0, |a| a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(0.0)).fold(0.0, f64::max)); + let mut step = 0.15 * reach; + let mut best = quality(self, v); + let mut cur = p0; + let mut moved = false; + for _ in 0..12 { + if step <= 1e-9 * reach.max(1e-300) { + break; + } + let mut improved = false; + for d in [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0]] { + let dir = self.restrict(v, d.map(|x| x * step)); + if norm(dir) == 0.0 { + continue; + } + let trial = self.project(v, [0, 1, 2].map(|k| cur[k] + dir[k])); + if let Some(p) = self.points.get_mut(v as usize) { + *p = trial; + } + let q = quality(self, v); + if q > best + 1e-12 { + best = q; + cur = trial; + improved = true; + moved = true; + } else if let Some(p) = self.points.get_mut(v as usize) { + *p = cur; + } + } + if !improved { + step *= 0.5; + } + } + if let Some(p) = self.points.get_mut(v as usize) { + *p = cur; + } + moved + } + + /// The smallest certified Jacobian ratio of the elements around `v` (negative where the + /// trilinear map folds somewhere, even with positive corners). + fn certified_quality(&self, v: u32) -> f64 { + self.around.get(v as usize).map_or(f64::INFINITY, |a| { + a.iter() + .map(|&(e, _)| { + let x = self.element(e).0; + if min_scaled_jacobian(&x) <= 0.0 { + return -1e9; + } + element_validity(ElementKind::Hex8, &x).map_or(-1.0, |c| match c.validity { + Validity::Valid => c.ratio.max(1e-300), + _ => c.ratio.min(0.0) - 1e-12, + }) + }) + .fold(f64::INFINITY, f64::min) + }) + } + + /// Repair elements whose trilinear map is not proven positive although their corners are: + /// the vertices of each such element (free ones, and boundary ones along their features) + /// move by compass search to raise the smallest certified bound around them; `rounds` rounds + /// at most. Returns the elements still not proven valid. + pub fn repair_certified(&mut self, rounds: usize) -> Vec { + let check = |o: &Self, e: u32| element_validity(ElementKind::Hex8, &o.element(e).0).is_ok_and(|c| c.validity == Validity::Valid); + let mut bad: Vec = (0..self.hexes.len() as u32).filter(|&e| !check(self, e)).collect(); + for _ in 0..rounds { + if bad.is_empty() { + break; + } + let mut verts: Vec = bad.iter().flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8])).collect(); + verts.sort_unstable(); + verts.dedup(); + for v in verts { + if !matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_)) | None) { + self.compass(v, &|o: &Self, v| o.certified_quality(v)); + } + } + let mut around: Vec = bad + .iter() + .flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8])) + .flat_map(|v| self.around.get(v as usize).map_or(Vec::new(), |a| a.iter().map(|x| x.0).collect())) + .collect(); + around.sort_unstable(); + around.dedup(); + bad = around.into_iter().filter(|&e| !check(self, e)).collect(); + } + bad + } +} + +/// Solve the symmetric 3 × 3 system `h x = b` (Cramer), `None` when singular. +fn solve3(h: [[f64; 3]; 3], b: P) -> Option

{ + let d = dot(h[0], cross(h[1], h[2])); + if !(d.abs() > 1e-300) || !d.is_finite() { + return None; + } + let col = |k: usize| [h[0][k], h[1][k], h[2][k]]; + let (c0, c1, c2) = (col(0), col(1), col(2)); + Some([dot(b, cross(c1, c2)) / d, dot(c0, cross(b, c2)) / d, dot(c0, cross(c1, b)) / d]) +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::octree::geom::tests::box_surface; + + #[test] + fn corner_frames_are_right_handed_on_the_unit_cube() { + let cube: [P; 8] = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]]; + for c in 0..8 { + let j = corner_jacobian(&cube, c, 1.0); + assert!((det3(&j) - 1.0).abs() < 1e-15, "corner {c}"); + } + assert!((min_scaled_jacobian(&cube) - 1.0).abs() < 1e-15); + } + + #[test] + fn a_tangled_interior_vertex_is_untangled() { + // A 2 × 2 × 2 block of cubes with the centre vertex pushed through a face; the boundary + // fixed. + let s = box_surface([0.0; 3], [2.0; 3]); + let g = Geometry::new(&s, 30.0, 160.0).unwrap(); + let id = |x: u32, y: u32, z: u32| x + 3 * (y + 3 * z); + let mut p = Vec::new(); + for z in 0..3 { + for y in 0..3 { + for x in 0..3 { + p.push([x as f64, y as f64, z as f64]); + } + } + } + let mut h = Vec::new(); + for z in 0..2 { + for y in 0..2 { + for x in 0..2 { + h.push([ + id(x, y, z), + id(x + 1, y, z), + id(x + 1, y + 1, z), + id(x, y + 1, z), + id(x, y, z + 1), + id(x + 1, y, z + 1), + id(x + 1, y + 1, z + 1), + id(x, y + 1, z + 1), + ]); + } + } + } + let centre = id(1, 1, 1) as usize; + p[centre] = [1.7, 1.9, 0.2]; + let size = vec![1.0; h.len()]; + // Boundary vertices slide on the box's faces (hard constraints, patch of the nearest + // face), the centre is free. + let cons: Vec = (0..p.len()) + .map(|i| { + if i == centre { + return Constraint::Free; + } + let t = g.nearest(p[i], |_| true).unwrap().0; + Constraint::Patch(g.tri_patch[t as usize]) + }) + .collect(); + let mut o = Optimizer::new(&g, p, &h, &size, cons); + assert!(o.measure().0 < 0.0); + let st = o.run(200, 1e-9); + assert!(o.measure().0 > 0.9, "{st:?}"); + assert!(sub(o.points[centre], [1.0; 3]).iter().all(|c| c.abs() < 1e-3), "{:?}", o.points[centre]); + } +} diff --git a/crates/mesh/src/octree/templates.rs b/crates/mesh/src/octree/templates.rs new file mode 100644 index 0000000..78cdedd --- /dev/null +++ b/crates/mesh/src/octree/templates.rs @@ -0,0 +1,374 @@ +//! The 3-refinement transition templates: how a leaf of the balanced tree is cut into +//! hexahedra so that it conforms to finer neighbours. +//! +//! Every vertex of a leaf is *marked* when a finer leaf touches it. Under strong balance the +//! marked set of a leaf is one of five patterns (after the closure in [`super::grid`] has split +//! every leaf that had another one), each with one template: +//! +//! | marks | template | hexahedra | +//! |---|---|---| +//! | none | the leaf itself | 1 | +//! | one corner | corner cube and three blocks fanning to the far corner | 4 | +//! | the two ends of an edge | three cubes along the edge and two four-hex slabs | 11 | +//! | the four corners of a face | nine cubes on the face, nine above them, four on top | 22 | +//! | all eight | the 3 × 3 × 3 split | 27 | +//! +//! A face of a leaf is cut according to its own marks only (none, one corner, one edge, all +//! four), and an edge into thirds at its marked ends, so two leaves sharing a face cut it the same +//! way: the mesh is conforming by construction. This is the node-marking form of Schneiders' +//! 3-refinement (R. Schneiders, *Refining quadrilateral and hexahedral element meshes*, NUMGRID +//! 1996; Y. Ito, A. M. Shih, B. K. Soni, *Octree-based reasonable-quality hexahedral mesh +//! generation using a new set of refinement templates*, IJNME 77(13):1809–1833, 2009). The +//! templates are our own construction from those papers' description: the face template's top +//! block is the 3 × 3 → 1 transition traced from its dual sheets, the edge template is two +//! copies of the planar 3 → 1 transition swept across the leaf, and the tests prove every +//! template's faces, volume and validity. +//! +//! Coordinates are integers in sixths of the leaf's edge, so every vertex shared with a neighbour +//! lies on the global integer lattice of [`super::grid`] and coincident vertices merge exactly. + +/// A point of a leaf in sixths of its edge, `0..=6` per axis. +pub(crate) type Local = [i64; 3]; +/// A hexahedron in the model's node order (corners 1-2-3-4 then 5-6-7-8 above them). +pub(crate) type Hex = [Local; 8]; + +/// The marked-corner pattern of a leaf, up to symmetry. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Pattern { + None, + Corner, + Edge, + Face, + Full, +} + +/// A symmetry of the cube: coordinate `k` of the image is coordinate `perm[k]` of the point, +/// reflected (`6 − v`) when `flip[k]`. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) struct Sym { + perm: [usize; 3], + flip: [bool; 3], +} + +impl Sym { + pub(crate) fn apply(&self, p: Local) -> Local { + let mut q = [0; 3]; + for (k, out) in q.iter_mut().enumerate() { + let v = p[self.perm[k]]; + *out = if self.flip[k] { 6 - v } else { v }; + } + q + } +} + +/// The 48 symmetries in a fixed order. +pub(crate) fn symmetries() -> Vec { + let perms = [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]]; + let mut out = Vec::with_capacity(48); + for perm in perms { + for bits in 0..8u8 { + out.push(Sym { perm, flip: [bits & 1 != 0, bits & 2 != 0, bits & 4 != 0] }); + } + } + out +} + +/// The leaf corner `c` (bit 0: x, bit 1: y, bit 2: z) in sixths. +pub(crate) fn corner(c: usize) -> Local { + [(c & 1) as i64 * 6, ((c >> 1) & 1) as i64 * 6, ((c >> 2) & 1) as i64 * 6] +} + +/// The corner index of a corner point. +fn corner_index(p: Local) -> usize { + usize::from(p[0] == 6) | (usize::from(p[1] == 6) << 1) | (usize::from(p[2] == 6) << 2) +} + +/// The image of a corner mask under a symmetry. +fn map_mask(sym: &Sym, mask: u8) -> u8 { + (0..8).filter(|c| mask & (1 << c) != 0).fold(0u8, |m, c| m | (1 << corner_index(sym.apply(corner(c))))) +} + +/// The pattern of a corner mask, or `None` when no template fits it. +pub(crate) fn classify(mask: u8) -> Option { + let corners: Vec = (0..8).filter(|c| mask & (1 << c) != 0).map(corner).collect(); + let shared = |n: usize| (0..3).filter(|&k| corners.iter().all(|p| p[k] == corners[0][k])).count() == n; + match corners.len() { + 0 => Some(Pattern::None), + 1 => Some(Pattern::Corner), + 2 if shared(2) => Some(Pattern::Edge), + 4 if shared(1) => Some(Pattern::Face), + 8 => Some(Pattern::Full), + _ => None, + } +} + +/// The canonical mask of a pattern: corner 0; edge 0–1 (along x); face z = 0. +fn canonical_mask(p: Pattern) -> u8 { + match p { + Pattern::None => 0, + Pattern::Corner => 0b1, + Pattern::Edge => 0b11, + Pattern::Face => 0b1111, + Pattern::Full => 0xff, + } +} + +/// The first symmetry (in [`symmetries`] order) taking the pattern's canonical mask to `mask`. +pub(crate) fn orientation(syms: &[Sym], pattern: Pattern, mask: u8) -> Option { + let c = canonical_mask(pattern); + syms.iter().find(|s| map_mask(s, c) == mask).copied() +} + +/// Thirds to sixths. +fn t(x: i64, y: i64, z: i64) -> Local { + [2 * x, 2 * y, 2 * z] +} + +/// A thirds-aligned cube with lower corner `(x, y, z)` (in thirds) and edge `e` thirds. +fn cube(x: i64, y: i64, z: i64, e: i64) -> Hex { + let (a, b) = (t(x, y, z), t(x + e, y + e, z + e)); + [ + [a[0], a[1], a[2]], + [b[0], a[1], a[2]], + [b[0], b[1], a[2]], + [a[0], b[1], a[2]], + [a[0], a[1], b[2]], + [b[0], a[1], b[2]], + [b[0], b[1], b[2]], + [a[0], b[1], b[2]], + ] +} + +/// The planar 3 → 1 transition in `(x, s)`: x in thirds, `s` in halves of the transition height, +/// three segments at `s = 0` and one at `s = 2`. +const THREE_TO_ONE: [[(i64, i64); 4]; 4] = + [[(0, 0), (1, 0), (1, 1), (0, 2)], [(1, 0), (2, 0), (2, 1), (1, 1)], [(2, 0), (3, 0), (3, 2), (2, 1)], [(1, 1), (2, 1), (3, 2), (0, 2)]]; + +/// The canonical template of a pattern, in sixths of the leaf. +pub(crate) fn template(p: Pattern) -> Vec { + match p { + Pattern::None => vec![cube(0, 0, 0, 3)], + Pattern::Full => (0..27).map(|k| cube(k % 3, (k / 3) % 3, k / 9, 1)).collect(), + Pattern::Corner => vec![ + cube(0, 0, 0, 1), + [t(1, 0, 0), t(3, 0, 0), t(3, 3, 0), t(1, 1, 0), t(1, 0, 1), t(3, 0, 3), t(3, 3, 3), t(1, 1, 1)], + [t(0, 1, 0), t(1, 1, 0), t(3, 3, 0), t(0, 3, 0), t(0, 1, 1), t(1, 1, 1), t(3, 3, 3), t(0, 3, 3)], + [t(0, 0, 1), t(1, 0, 1), t(1, 1, 1), t(0, 1, 1), t(0, 0, 3), t(3, 0, 3), t(3, 3, 3), t(0, 3, 3)], + ], + Pattern::Edge => { + let mut out: Vec = (0..3).map(|i| cube(i, 0, 0, 1)).collect(); + // Two slabs over the edge column's faces y = 1 and z = 1, each the 3 → 1 transition + // swept from the leaf face to the diagonal plane y = z. + for q in THREE_TO_ONE { + let b = |(x, s): (i64, i64), w: i64| t(x, 1 + s, w * (1 + s)); + let c = |(x, s): (i64, i64), w: i64| t(x, w * (1 + s), 1 + s); + out.push([b(q[0], 0), b(q[1], 0), b(q[2], 0), b(q[3], 0), b(q[0], 1), b(q[1], 1), b(q[2], 1), b(q[3], 1)]); + out.push([c(q[0], 0), c(q[1], 0), c(q[2], 0), c(q[3], 0), c(q[0], 1), c(q[1], 1), c(q[2], 1), c(q[3], 1)]); + } + out + } + Pattern::Face => { + let mut out: Vec = (0..9).map(|k| cube(k % 3, k / 3, 0, 1)).collect(); + // The layer above: its top is the corners raised to the leaf's top, the side nodes of + // the 3 → 1 transitions at two thirds, the inside nodes at half height. + let top = |a: i64, b: i64| -> Local { + let edge = |v: i64| v == 0 || v == 3; + match (edge(a), edge(b)) { + (true, true) => [2 * a, 2 * b, 6], + (false, false) => [2 * a, 2 * b, 3], + _ => [2 * a, 2 * b, 4], + } + }; + for k in 0..9 { + let (i, j) = (k % 3, k / 3); + out.push([ + t(i, j, 1), + t(i + 1, j, 1), + t(i + 1, j + 1, 1), + t(i, j + 1, 1), + top(i, j), + top(i + 1, j), + top(i + 1, j + 1), + top(i, j + 1), + ]); + } + // The four on top: the y-arch (low), the x-arch (high) and the two side blocks. + out.push([[2, 2, 3], [4, 2, 3], [4, 4, 3], [2, 4, 3], [2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4]]); + out.push([[2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4], [0, 0, 6], [6, 0, 6], [6, 6, 6], [0, 6, 6]]); + out.push([[2, 2, 3], [2, 4, 3], [2, 6, 4], [2, 0, 4], [0, 2, 4], [0, 4, 4], [0, 6, 6], [0, 0, 6]]); + out.push([[4, 2, 3], [4, 4, 3], [4, 6, 4], [4, 0, 4], [6, 2, 4], [6, 4, 4], [6, 6, 6], [6, 0, 6]]); + out + } + } +} + +/// Six times the signed volume of the corner tetrahedron at node 0 (exact on the lattice). +pub(crate) fn corner0_orientation(h: &[[i64; 3]; 8]) -> i128 { + let e = |j: usize| [0, 1, 2].map(|k| i128::from(h[j][k] - h[0][k])); + let (a, b, c) = (e(1), e(3), e(4)); + a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0]) +} + +/// The hexahedron with its orientation reversed (the same cell, mirrored node order). +pub(crate) fn reversed(h: [T; 8]) -> [T; 8] { + [h[0], h[3], h[2], h[1], h[4], h[7], h[6], h[5]] +} + +#[cfg(test)] +mod tests { + use std::collections::HashMap; + + use super::*; + use crate::quality::{Validity, element_validity}; + use fieldcraft_model::ElementKind; + + /// The face of the leaf a quad lies on (axis, value), if any. + fn on_leaf_face(q: &[Local; 4]) -> Option<(usize, i64)> { + (0..3).find_map(|k| [0, 6].into_iter().find(|&v| q.iter().all(|p| p[k] == v)).map(|v| (k, v))) + } + + fn sorted(mut q: [Local; 4]) -> [Local; 4] { + q.sort(); + q + } + + /// The quads a leaf face must carry for its four corner marks, as sorted vertex sets, on the + /// face `(axis, value)`. + fn face_template(axis: usize, value: i64, marks: [bool; 4]) -> Vec<[Local; 4]> { + // Face coordinates (u, v) are the other two axes in increasing order; corners in order + // (0,0), (6,0), (0,6), (6,6). + let (ua, va) = match axis { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + let lift = |u: i64, v: i64| { + let mut p = [0; 3]; + p[axis] = value; + p[ua] = u; + p[va] = v; + p + }; + // Build in a canonical frame, then map the 2D symmetry that takes the canonical marks to + // these marks. + let count = marks.iter().filter(|m| **m).count(); + let canon: Vec<[(i64, i64); 4]> = match count { + 0 => vec![[(0, 0), (6, 0), (6, 6), (0, 6)]], + 1 => vec![[(0, 0), (2, 0), (2, 2), (0, 2)], [(2, 0), (6, 0), (6, 6), (2, 2)], [(0, 2), (2, 2), (6, 6), (0, 6)]], + 2 => { + let mut v: Vec<[(i64, i64); 4]> = (0..3).map(|i| [(2 * i, 0), (2 * i + 2, 0), (2 * i + 2, 2), (2 * i, 2)]).collect(); + v.extend([ + [(0, 2), (2, 2), (2, 4), (0, 6)], + [(2, 2), (4, 2), (4, 4), (2, 4)], + [(4, 2), (6, 2), (6, 6), (4, 4)], + [(2, 4), (4, 4), (6, 6), (0, 6)], + ]); + v + } + _ => (0..9) + .map(|k| { + [(2 * (k % 3), 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3) + 2), (2 * (k % 3), 2 * (k / 3) + 2)] + }) + .collect(), + }; + let canon_marks: [bool; 4] = match count { + 0 => [false; 4], + 1 => [true, false, false, false], + 2 => [true, true, false, false], + _ => [true; 4], + }; + let maps: [fn(i64, i64) -> (i64, i64); 8] = [ + |u, v| (u, v), + |u, v| (6 - u, v), + |u, v| (u, 6 - v), + |u, v| (6 - u, 6 - v), + |u, v| (v, u), + |u, v| (6 - v, u), + |u, v| (v, 6 - u), + |u, v| (6 - v, 6 - u), + ]; + let corners = [(0, 0), (6, 0), (0, 6), (6, 6)]; + let idx = |c: (i64, i64)| corners.iter().position(|&x| x == c).unwrap(); + let m = + maps.iter().find(|m| (0..4).all(|c| canon_marks[c] == marks[idx(m(corners[c].0, corners[c].1))])).expect("a face symmetry for the marks"); + canon.iter().map(|q| sorted(q.map(|(u, v)| m(u, v)).map(|(u, v)| lift(u, v)))).collect() + } + + fn face_corner_marks(axis: usize, value: i64, mask: u8) -> [bool; 4] { + let (ua, va) = match axis { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + [(0, 0), (6, 0), (0, 6), (6, 6)].map(|(u, v)| { + let mut p = [0; 3]; + p[axis] = value; + p[ua] = u; + p[va] = v; + mask & (1 << corner_index(p)) != 0 + }) + } + + /// Every template under every symmetry: valid elements, the leaf's volume, closed (each inner + /// face shared by two elements) and, on each leaf face, exactly the face template of that + /// face's marks. + #[test] + fn templates_are_valid_closed_and_conform_to_their_face_marks() { + let syms = symmetries(); + for pattern in [Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face, Pattern::Full] { + for sym in &syms { + let mask = map_mask(sym, canonical_mask(pattern)); + assert_eq!(classify(mask), Some(pattern)); + let hexes: Vec = template(pattern) + .into_iter() + .map(|h| { + let h = h.map(|p| sym.apply(p)); + if corner0_orientation(&h) < 0 { reversed(h) } else { h } + }) + .collect(); + let mut volume = 0.0; + let mut faces: HashMap<[Local; 4], usize> = HashMap::new(); + for h in &hexes { + let xyz: Vec<[f64; 3]> = h.iter().map(|p| p.map(|v| v as f64)).collect(); + let check = element_validity(ElementKind::Hex8, &xyz).unwrap(); + assert_eq!(check.validity, Validity::Valid, "{pattern:?} {sym:?} {h:?}"); + volume += check.volume; + for f in ElementKind::Hex8.faces() { + *faces.entry(sorted([h[f[0]], h[f[1]], h[f[2]], h[f[3]]])).or_default() += 1; + } + } + assert!((volume - 216.0).abs() < 1e-9, "{pattern:?}: volume {volume}"); + let mut outer: HashMap<(usize, i64), Vec<[Local; 4]>> = HashMap::new(); + for (q, n) in &faces { + match on_leaf_face(q) { + Some(f) => { + assert_eq!(*n, 1); + outer.entry(f).or_default().push(*q); + } + None => assert_eq!(*n, 2, "{pattern:?} {sym:?}: inner face {q:?} used {n} times"), + } + } + for axis in 0..3 { + for value in [0, 6] { + let mut got = outer.remove(&(axis, value)).unwrap_or_default(); + let mut want = face_template(axis, value, face_corner_marks(axis, value, mask)); + got.sort(); + want.sort(); + assert_eq!(got, want, "{pattern:?} {sym:?} face {axis}={value}"); + } + } + } + } + } + + #[test] + fn patterns_without_a_template_are_recognised() { + assert_eq!(classify(0b1001), None, "diagonal of a face"); + assert_eq!(classify(0b1000_0001), None, "opposite corners"); + assert_eq!(classify(0b0111), None, "three of a face"); + assert_eq!(classify(0b0110_1001), None, "four not on a face"); + assert_eq!(classify(0b0111_1111), None); + assert_eq!(classify(0b0011_0011), Some(Pattern::Face)); + assert_eq!(classify(0b0101_0000), Some(Pattern::Edge)); + } +} diff --git a/crates/mesh/src/octree/tests.rs b/crates/mesh/src/octree/tests.rs new file mode 100644 index 0000000..7fa388b --- /dev/null +++ b/crates/mesh/src/octree/tests.rs @@ -0,0 +1,98 @@ +use super::*; +use crate::octree::geom::tests::box_surface; +use crate::size::UniformSize; + +#[test] +fn a_box_meshes_with_exact_volume_and_tagged_faces() { + let s = box_surface([0.0; 3], [10.0, 8.0, 6.0]); + let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap(); + let r = &m.report; + assert!((r.volume - 480.0).abs() < 1e-9 * 480.0, "{r:?}"); + assert!(r.min_scaled_jacobian > 0.3, "{r:?}"); + let faces: BTreeSet = m.boundary.iter().map(|b| b.brep_face).collect(); + assert_eq!(faces.len(), 6); + // Each boundary quad lies in the plane of its face. + for b in &m.boundary { + let hx = m.hexes[b.element as usize]; + let q = topo::FACES[(b.face - 1) as usize].map(|i| m.nodes[hx[i] as usize]); + let (axis, value) = match b.brep_face { + 0 => (2, 0.0), + 1 => (2, 6.0), + 2 => (1, 0.0), + 3 => (1, 8.0), + 4 => (0, 0.0), + _ => (0, 10.0), + }; + assert!(q.iter().all(|p| (p[axis] - value).abs() < 1e-9), "{b:?} {q:?}"); + } +} + +/// A closed cylinder of radius `r` and height `h` on the z axis, `n` segments: faces 0 (bottom +/// cap), 1 (top cap), 2 (side). +pub(crate) fn cylinder_surface(r: f64, h: f64, n: usize) -> BoundarySurface { + let mut p = vec![[0.0, 0.0, 0.0], [0.0, 0.0, h]]; + for z in [0.0, h] { + for i in 0..n { + let a = std::f64::consts::TAU * i as f64 / n as f64; + p.push([r * a.cos(), r * a.sin(), z]); + } + } + let (b, t) = (|i: usize| (2 + i % n) as u32, |i: usize| (2 + n + i % n) as u32); + let mut tris = Vec::new(); + let mut face = Vec::new(); + for i in 0..n { + tris.push([0, b(i + 1), b(i)]); + face.push(0); + tris.push([1, t(i), t(i + 1)]); + face.push(1); + tris.push([b(i), b(i + 1), t(i + 1)]); + tris.push([b(i), t(i + 1), t(i)]); + face.extend([2, 2]); + } + BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap() +} + +/// A UV sphere of radius `r`, one face. +pub(crate) fn sphere_surface(r: f64, n: usize) -> BoundarySurface { + let rings = n / 2; + let mut p = vec![[0.0, 0.0, -r], [0.0, 0.0, r]]; + for j in 1..rings { + let th = std::f64::consts::PI * j as f64 / rings as f64 - std::f64::consts::FRAC_PI_2; + for i in 0..n { + let a = std::f64::consts::TAU * i as f64 / n as f64; + p.push([r * th.cos() * a.cos(), r * th.cos() * a.sin(), r * th.sin()]); + } + } + let id = |j: usize, i: usize| (2 + (j - 1) * n + i % n) as u32; + let mut tris = Vec::new(); + for i in 0..n { + tris.push([0, id(1, i + 1), id(1, i)]); + tris.push([1, id(rings - 1, i), id(rings - 1, i + 1)]); + } + for j in 1..rings - 1 { + for i in 0..n { + tris.push([id(j, i), id(j, i + 1), id(j + 1, i + 1)]); + tris.push([id(j, i), id(j + 1, i + 1), id(j + 1, i)]); + } + } + let face = vec![0; tris.len()]; + BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap() +} + +#[test] +fn a_cylinder_meshes_with_its_rims_captured() { + let s = cylinder_surface(5.0, 8.0, 64); + let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap(); + let r = &m.report; + assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}"); + let faces: BTreeSet = m.boundary.iter().map(|b| b.brep_face).collect(); + assert_eq!(faces.len(), 3); +} + +#[test] +fn a_sphere_meshes() { + let s = sphere_surface(5.0, 48); + let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap(); + let r = &m.report; + assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}"); +} diff --git a/crates/mesh/src/octree/topo.rs b/crates/mesh/src/octree/topo.rs new file mode 100644 index 0000000..1877717 --- /dev/null +++ b/crates/mesh/src/octree/topo.rs @@ -0,0 +1,404 @@ +//! Hexahedral mesh topology for the boundary stages: boundary faces, the manifold clean-up of +//! the core, padding a boundary with a layer of hexahedra, and pillowing a shrink set. +//! +//! *Padding* puts one hexahedron on every boundary quad, between the quad and a copy of it, so +//! no hexahedron keeps more than one face on the boundary. *Pillowing* a shrink set `S` inserts +//! a layer of hexahedra along the faces between `S` and the rest of the mesh; the layer reaches +//! the boundary along the loop where `S`'s boundary faces meet the others' (S. A. Mitchell, +//! T. J. Tautges, *Pillowing doublets: refining a mesh to ensure that faces share at most one +//! edge*, IMR 1995; Gao, Shen, Panozzo 2019, §3.4 use the same insertion around features). + +use std::collections::{BTreeMap, HashMap}; + +use crate::{Error, Result}; + +/// The six faces of a hex8 as local node positions, right-hand normals outward (the model's +/// `ElementKind::Hex8.faces()` order). +pub(crate) const FACES: [[usize; 4]; 6] = [[0, 3, 2, 1], [4, 5, 6, 7], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7]]; + +/// A face key: the vertices sorted. +pub(crate) fn key(q: [u32; 4]) -> [u32; 4] { + let mut k = q; + k.sort_unstable(); + k +} + +/// The root of `i` in the union-find forest `c`, compressing the path from `i`. +pub(crate) fn find(c: &mut [usize], i: usize) -> usize { + let mut r = i; + while let Some(&p) = c.get(r).filter(|&&p| p != r) { + r = p; + } + if let Some(x) = c.get_mut(i) { + *x = r; + } + r +} + +/// A face of an element: (hexahedron, local face 0..6). +pub(crate) type FaceRef = (u32, u8); + +/// The vertices of a local face, outward. +pub(crate) fn face(h: &[u32; 8], f: u8) -> [u32; 4] { + FACES.get(f as usize).map_or([0; 4], |l| l.map(|i| h[i])) +} + +/// Every face's elements, keyed by sorted vertices (deterministic map). +pub(crate) fn face_map(hexes: &[[u32; 8]]) -> BTreeMap<[u32; 4], Vec> { + let mut m: BTreeMap<[u32; 4], Vec> = BTreeMap::new(); + for (i, h) in hexes.iter().enumerate() { + for f in 0..6u8 { + m.entry(key(face(h, f))).or_default().push((i as u32, f)); + } + } + m +} + +/// The boundary faces (on one element only), outward, in face-map order. +pub(crate) fn boundary(hexes: &[[u32; 8]]) -> Vec<(FaceRef, [u32; 4])> { + face_map(hexes) + .into_values() + .filter(|v| v.len() == 1) + .filter_map(|v| v.first().copied()) + .map(|(h, f)| ((h, f), hexes.get(h as usize).map_or([0; 4], |x| face(x, f)))) + .collect() +} + +/// Edges of a quad list with the quads using each (ascending vertex pair → quad indices). +pub(crate) fn quad_edges(quads: &[[u32; 4]]) -> BTreeMap<(u32, u32), Vec> { + let mut m: BTreeMap<(u32, u32), Vec> = BTreeMap::new(); + for (i, q) in quads.iter().enumerate() { + for k in 0..4 { + let (a, b) = (q[k], q[(k + 1) % 4]); + m.entry((a.min(b), a.max(b))).or_default().push(i as u32); + } + } + m +} + +/// The hexahedra to drop so that the boundary of `keep` is a 2-manifold: elements at edges with +/// four boundary faces (elements touching along an edge only) and at vertices whose boundary +/// faces form more than one fan. Returns the indices to drop (empty when manifold). +pub(crate) fn nonmanifold(hexes: &[[u32; 8]], keep: &[bool]) -> Vec { + let kept: Vec = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect(); + let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect(); + let b = boundary(&sub); + let quads: Vec<[u32; 4]> = b.iter().map(|(_, q)| *q).collect(); + let mut bad_vertices: Vec = Vec::new(); + for ((a, c), qs) in quad_edges(&quads) { + if qs.len() != 2 { + bad_vertices.extend([a, c]); + } + } + // Fans: around each vertex, boundary quads connected through shared edges at the vertex. + let mut at: BTreeMap> = BTreeMap::new(); + for (i, q) in quads.iter().enumerate() { + for &v in q { + at.entry(v).or_default().push(i as u32); + } + } + for (v, qs) in &at { + let mut comp: Vec = (0..qs.len()).collect(); + for i in 0..qs.len() { + for j in i + 1..qs.len() { + let (qa, qb) = (quads[qs[i] as usize], quads[qs[j] as usize]); + // Share an edge through v: another common vertex adjacent to v in both. + let nb = |q: [u32; 4]| { + let p = q.iter().position(|&x| x == *v).unwrap_or(0); + [q[(p + 1) % 4], q[(p + 3) % 4]] + }; + if nb(qa).iter().any(|x| nb(qb).contains(x)) { + let (ri, rj) = (find(&mut comp, i), find(&mut comp, j)); + comp[ri] = rj; + } + } + } + let roots = (0..qs.len()).filter(|&i| find(&mut comp, i) == i).count(); + if roots > 1 { + bad_vertices.push(*v); + } + } + bad_vertices.sort_unstable(); + bad_vertices.dedup(); + if bad_vertices.is_empty() { + return Vec::new(); + } + kept.into_iter().filter(|&i| hexes.get(i as usize).is_some_and(|h| h.iter().any(|v| bad_vertices.binary_search(v).is_ok()))).collect() +} + +/// Pad: a hexahedron on each listed boundary face (outward), between it and a copy; the copies +/// of the faces' vertices are new points (appended, initially at the original positions). +/// Returns the copy of each original vertex. +pub(crate) fn pad(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, faces: &[[u32; 4]]) -> Result> { + let mut copy: BTreeMap = BTreeMap::new(); + for q in faces { + for &v in q { + if let std::collections::btree_map::Entry::Vacant(slot) = copy.entry(v) { + let n = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?; + points.push(points.get(v as usize).copied().unwrap_or([0.0; 3])); + slot.insert(n); + } + } + } + for q in faces { + let c = q.map(|v| copy.get(&v).copied().unwrap_or(v)); + hexes.push([q[0], q[1], q[2], q[3], c[0], c[1], c[2], c[3]]); + } + Ok(copy) +} + +/// One vertex copy made by [`pillow`]: the original, the copy, and the shrink-set elements +/// around the original that now use the copy. +#[derive(Debug, Clone)] +pub(crate) struct PillowCopy { + pub original: u32, + pub copy: u32, + pub elements: Vec, +} + +/// Pillow the shrink set `inside` (one flag per element): a new layer of hexahedra on every face +/// between an element of the set and one outside it. Vertices of those faces get one copy per +/// group of set elements around them (connected through faces inside the set), which the set's +/// elements now use; the new layer joins the originals to the copies. Faces on the boundary are +/// not pillowed: where the layer meets the boundary, the set's boundary faces move to the copies, +/// so the boundary gains a ring of quads along the loop. New points are appended at the +/// originals' positions. Returns the copies made. +pub(crate) fn pillow(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, inside: &[bool]) -> Result> { + let fm = face_map(hexes); + let in_set = |h: u32| inside.get(h as usize).copied().unwrap_or(false); + // Faces to pillow, oriented from the outside element into the set. + let mut sheet: Vec<(u32, [u32; 4])> = Vec::new(); + for refs in fm.values() { + if let [(a, fa), (b, fb)] = refs.as_slice() { + match (in_set(*a), in_set(*b)) { + (true, false) => sheet.push((*a, hexes.get(*a as usize).map_or([0; 4], |h| face(h, *fa)))), + (false, true) => sheet.push((*b, hexes.get(*b as usize).map_or([0; 4], |h| face(h, *fb)))), + _ => {} + } + } + } + if sheet.is_empty() { + return Ok(Vec::new()); + } + let mut on_sheet: Vec = sheet.iter().flat_map(|(_, q)| *q).collect(); + on_sheet.sort_unstable(); + on_sheet.dedup(); + // Set elements around each sheet vertex. + let mut around: HashMap> = HashMap::new(); + for (i, h) in hexes.iter().enumerate() { + if !in_set(i as u32) { + continue; + } + for &v in h { + if on_sheet.binary_search(&v).is_ok() { + let e = around.entry(v).or_default(); + if e.last() != Some(&(i as u32)) { + e.push(i as u32); + } + } + } + } + let mut copies: Vec = Vec::new(); + // (vertex, element) → copy. + let mut copy_of: HashMap<(u32, u32), u32> = HashMap::new(); + for &v in &on_sheet { + let mut els = around.get(&v).cloned().unwrap_or_default(); + els.sort_unstable(); + els.dedup(); + // Group by faces inside the set that contain v. + let mut group: Vec = (0..els.len()).collect(); + for i in 0..els.len() { + for j in i + 1..els.len() { + let (Some(hi), Some(hj)) = (hexes.get(els[i] as usize), hexes.get(els[j] as usize)) else { continue }; + let shares = (0..6u8).any(|f| { + let q = face(hi, f); + q.contains(&v) && (0..6u8).any(|g| key(face(hj, g)) == key(q)) + }); + if shares { + let (a, b) = (find(&mut group, i), find(&mut group, j)); + group[a] = b; + } + } + } + // Per group: its copy and that copy's index in `copies`. + let mut made: BTreeMap = BTreeMap::new(); + for (i, &e) in els.iter().enumerate() { + let r = find(&mut group, i); + let (c, at) = match made.get(&r) { + Some(&m) => m, + None => { + let c = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?; + points.push(points.get(v as usize).copied().unwrap_or([0.0; 3])); + copies.push(PillowCopy { original: v, copy: c, elements: Vec::new() }); + made.insert(r, (c, copies.len() - 1)); + (c, copies.len() - 1) + } + }; + copy_of.insert((v, e), c); + if let Some(pc) = copies.get_mut(at) { + pc.elements.push(e); + } + } + } + // The new layer: original face (seen from outside, normal into the set) below, copies above. + let mut layer: Vec<[u32; 8]> = Vec::with_capacity(sheet.len()); + for (h, q) in &sheet { + let r = [q[0], q[3], q[2], q[1]]; + let c = r.map(|v| copy_of.get(&(v, *h)).copied().unwrap_or(v)); + layer.push([r[0], r[1], r[2], r[3], c[0], c[1], c[2], c[3]]); + } + for (i, h) in hexes.iter_mut().enumerate() { + if in_set(i as u32) { + for v in h.iter_mut() { + if let Some(&c) = copy_of.get(&(*v, i as u32)) { + *v = c; + } + } + } + } + hexes.extend(layer); + Ok(copies) +} + +/// Hexahedra to drop so that every vertex of the boundary of `keep` can be pushed outward: the +/// normals of its boundary quads must lie strictly inside one half-space (their mean has a +/// positive product with each, at least `margin`). Where this fails, padding from the vertex +/// to the part's surface inverts elements, so for each such vertex the kept element at it with +/// the most boundary faces is dropped (the buffer-zone rule of H. Tong, E. Halilaj, Y. J. Zhang, +/// *HybridOctree_Hex*, arXiv 2401.05984, 2024, §2.3, in a mean-normal form). +pub(crate) fn hidden_vertices(hexes: &[[u32; 8]], keep: &[bool], xyz: &dyn Fn(u32) -> [f64; 3], margin: f64) -> Vec { + let kept: Vec = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect(); + let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect(); + let b = boundary(&sub); + let mut normals: BTreeMap> = BTreeMap::new(); + let mut faces_of: BTreeMap = BTreeMap::new(); + for ((h, _), q) in &b { + let p = q.map(xyz); + let d1 = [0, 1, 2].map(|k| p[2][k] - p[0][k]); + let d2 = [0, 1, 2].map(|k| p[3][k] - p[1][k]); + let n = [d1[1] * d2[2] - d1[2] * d2[1], d1[2] * d2[0] - d1[0] * d2[2], d1[0] * d2[1] - d1[1] * d2[0]]; + let l = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt(); + if l > 0.0 { + for &v in q { + normals.entry(v).or_default().push(n.map(|c| c / l)); + } + } + *faces_of.entry(*h).or_default() += 1; + } + let mut drop: Vec = Vec::new(); + for (v, ns) in &normals { + let mut m = [0.0; 3]; + for n in ns { + for k in 0..3 { + m[k] += n[k]; + } + } + let l = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]).sqrt(); + let ok = l > 0.0 && ns.iter().all(|n| (n[0] * m[0] + n[1] * m[1] + n[2] * m[2]) / l > margin); + if ok { + continue; + } + // The kept element at v with the most boundary faces (lowest index on ties). + let best = sub + .iter() + .enumerate() + .filter(|(_, h)| h.contains(v)) + .max_by(|a, b| faces_of.get(&(a.0 as u32)).unwrap_or(&0).cmp(faces_of.get(&(b.0 as u32)).unwrap_or(&0)).then(b.0.cmp(&a.0))) + .and_then(|(i, _)| kept.get(i).copied()); + if let Some(e) = best { + drop.push(e); + } + } + drop.sort_unstable(); + drop.dedup(); + drop +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn faces_follow_the_model_hex8_face_order() { + let model = fieldcraft_model::ElementKind::Hex8.faces(); + assert_eq!(model.len(), FACES.len()); + for (m, f) in model.iter().zip(FACES) { + assert_eq!(*m, &f[..]); + } + } + + /// An `n × n × n` block of unit cubes. + pub(crate) fn block(n: u32) -> (Vec<[f64; 3]>, Vec<[u32; 8]>) { + let id = |x: u32, y: u32, z: u32| x + (n + 1) * (y + (n + 1) * z); + let mut p = Vec::new(); + for z in 0..=n { + for y in 0..=n { + for x in 0..=n { + p.push([x as f64, y as f64, z as f64]); + } + } + } + let mut h = Vec::new(); + for z in 0..n { + for y in 0..n { + for x in 0..n { + h.push([ + id(x, y, z), + id(x + 1, y, z), + id(x + 1, y + 1, z), + id(x, y + 1, z), + id(x, y, z + 1), + id(x + 1, y, z + 1), + id(x + 1, y + 1, z + 1), + id(x, y + 1, z + 1), + ]); + } + } + } + (p, h) + } + + fn closed(hexes: &[[u32; 8]]) { + for (k, v) in face_map(hexes) { + assert!(v.len() <= 2, "face {k:?} on {} elements", v.len()); + } + let b = boundary(hexes); + for (_, qs) in quad_edges(&b.iter().map(|x| x.1).collect::>()) { + assert_eq!(qs.len(), 2); + } + } + + #[test] + fn padding_and_pillowing_keep_a_closed_conforming_mesh() { + let (mut p, mut h) = block(3); + let b: Vec<[u32; 4]> = boundary(&h).into_iter().map(|x| x.1).collect(); + assert_eq!(b.len(), 54); + let copy = pad(&mut p, &mut h, &b).unwrap(); + assert_eq!(copy.len(), 56); + assert_eq!(h.len(), 27 + 54); + closed(&h); + // Pillow the padding hexes on the face z = 0 (outer faces down). + let pad_start = 27; + let inside: Vec = (0..h.len()).map(|i| i >= pad_start && (0..4).all(|k| p[h[i][k] as usize][2] == 0.0)).collect(); + assert_eq!(inside.iter().filter(|x| **x).count(), 9); + let before = h.len(); + let copies = pillow(&mut p, &mut h, &inside).unwrap(); + closed(&h); + // Sheet: the 9 inner faces plus the 12 side faces to the neighbouring padding. + assert_eq!(h.len(), before + 9 + 12); + // The boundary gains the ring: 54 + 12 quads. + assert_eq!(boundary(&h).len(), 54 + 12); + assert!(!copies.is_empty()); + } + + #[test] + fn nonmanifold_edges_and_vertices_are_found() { + let (_, h) = block(2); + // Two cubes touching along an edge. + let keep = |ids: &[usize]| (0..8).map(|i| ids.contains(&i)).collect::>(); + assert!(!nonmanifold(&h, &keep(&[0, 3])).is_empty(), "edge contact"); + assert!(!nonmanifold(&h, &keep(&[0, 7])).is_empty(), "vertex contact"); + assert!(nonmanifold(&h, &keep(&[0, 1, 3])).is_empty()); + assert!(nonmanifold(&h, &keep(&[0, 1, 2, 3, 4, 5, 6, 7])).is_empty()); + } +} diff --git a/docs/architecture.md b/docs/architecture.md index 4df04e0..7818c96 100644 --- a/docs/architecture.md +++ b/docs/architecture.md @@ -68,12 +68,21 @@ Rules beyond the layer order: its B-rep face tag survives. Then constrained Delaunay refinement to a size field, ODT smoothing interleaved as in Tournois et al. 2009, and the Klingner–Shewchuk 2007 improvement schedule against slivers. -4. **Hexahedral meshing (the core).** Octree/grid dual all-hex first, the one family that meshes - any input validly (Livesu et al. 2021 dual schemes, Pitzalis et al. 2021 refinement, Gao et al. - 2019 feature projection and padding, Garanzha et al. 2021 untangling). Then hex-dominant by tet +4. **Hexahedral meshing (the core).** Grid-based all-hex first, the one family that meshes any + input validly. `fieldcraft-mesh::octree` (command `mesh.hex.octree`): a 27-tree over the part, + refined by curvature, thickness and narrow faces, balanced and cut by the 3-refinement + node-marking templates (Schneiders 1996, Ito et al. 2009) into a conforming grid; the core of + grid elements inside the part (exact inside test) is cleaned to a manifold boundary; the part's + patches, sharp curves and corners are mapped onto it (Gao et al. 2019); one padding layer and a + pillow under each patch carry the boundary and capture every sharp curve; the foldover-free + energy (Garanzha et al. 2021) with a growing fitting term moves the boundary onto the geometry + while every element stays valid; the worst elements are improved, and every element is proven + valid by the certified Bernstein Jacobian check (Johnen et al. 2013; Marschner et al. 2020) + before the boundary quads are tagged with their B-rep faces. Where a stage fails the grid is + refined there and the mesh rebuilt. The dual schemes of Livesu et al. 2021 (fewer elements, + concave transitions) are the next upgrade of the transitions. Then hex-dominant by tet recombination (Pellerin et al. 2018), then polycube (DualCube labelling) for blocky parts; - frame-field/integer-grid-map methods stay research. Validity uses a certified trilinear Jacobian - check (Marschner et al. 2020), not corner Jacobians alone. + frame-field/integer-grid-map methods stay research. 5. **Structured hex and face meshing.** Faces get all-quad meshes on shared, prescribed edge nodes (`fieldcraft-mesh::quad`): transfinite interpolation for four-sided faces (Gordon and Hall 1973), and for free faces a half-resolution constrained Delaunay triangulation, diff --git a/docs/meshing/literature.md b/docs/meshing/literature.md index 6d60515..8dfd24b 100644 --- a/docs/meshing/literature.md +++ b/docs/meshing/literature.md @@ -531,6 +531,53 @@ Delaunay methods with Voronoi/ODT improvement, then (2) all-hex and hex-dominant with conforming templates and no dualization. - **Use.** Lets us refine a coarse aligned hex mesh, for example polycube output, as well as an octree. - **Code.** `cg3hci/Gen-Adapt-Ref-for-Hexmeshing` is **MIT**: READ OK. +- **Finding (2026-10).** The reference pipeline chooses transition vertices with an integer linear + program (Gurobi) and installs CinoLib's schemes by heuristics that are retried with other + parameters when they fail; neither was adopted (see B1.10). + +**B1.7 Schneiders, Algorithms for quadrilateral and hexahedral mesh generation (VKI lecture notes, +2000), §5.** **Verdict: built on (the transition templates of `fieldcraft-mesh::octree`).** +- **Link.** https://www.robertschneiders.de/papers/vki.pdf (author's copy). +- **Method.** 3-refinement driven by nodal marks: an element is split by the template of its + marked corners, and faces split by their own marks only, so the refinement is conforming. +- **Limit (stated there).** The 2D templates of the original set make some 3D mark patterns + impossible: a hexahedron with seven marked corners has 51 boundary quads, an odd number, so no + hexahedral mesh of it exists (Mitchell's parity condition). Concave refined regions therefore + need their concave cells refined. The fix is a modified edge rule (an edge with a marked end + splits in three) with buffer layers in the face templates. + +**B1.8 Ito, Shih, Soni 2009, Octree-based reasonable-quality hexahedral mesh generation using a +new set of refinement templates.** **Verdict: reference.** +- **Citation.** IJNME 77(13):1809–1833, 2009. Templates for concave edges of the refined region + (paper not obtained; cited after Tong and Zhang 2025). + +**B1.9 Tong, Halilaj, Zhang 2024, HybridOctree_Hex.** **Verdict: reference (buffer-zone rule adopted).** +- **Citation.** H. Tong, E. Halilaj, Y. J. Zhang, *HybridOctree_Hex: Hybrid octree-based adaptive + all-hexahedral mesh generation with Jacobian control*, arXiv 2401.05984, 2024. +- **Method.** Dual of a strongly balanced, paired octree with five templates (a 13-hexahedron + frustum per transition face plus four transition-edge templates); a core cleared at the + boundary so that the normals of the quads around each core boundary vertex lie in one + half-space, then a buffer layer projected onto the surface and optimised. +- **Code.** `CMU-CBML/HybridOctree_Hex` is **BSD-3-Clause**: READ OK (hard-coded template + geometry, not adopted). + +**B1.10 Tong, Zhang 2025, MCHex.** **Verdict: reference (next boundary stage to evaluate).** +- **Citation.** H. Tong, Y. J. Zhang, *MCHex: Marching cubes based adaptive hexahedral mesh + generation with guaranteed positive Jacobian*, arXiv 2511.02064, 2025. +- **Method.** The same 3-refinement templates (B1.7); the boundary cells are cut by marching + cubes on a signed distance field, every cut cell's pieces are 3-regular polyhedra, and their + midpoint subdivision gives hexahedra with a proven positive Jacobian (Bernstein coefficients) and + a boundary that is pillowed by construction; 17 of the 31 cube cases are refined locally first. +- **Limits (stated there).** Sharp edges and corners are lost (marching cubes smooths them), and + the convex-only transitions over-refine. ModelCraft keeps the core-padding-pillow boundary with + feature capture for CAD parts; MCHex's cut cells are the candidate for organic shapes. + +**Implemented (`fieldcraft-mesh::octree`, `mesh.hex.octree`, 2026-10).** B1.7 templates on a +27-tree refined by curvature, thickness and narrow faces; a core cleared by the B1.9 rule; the +feature map of B1.4 (patches first); padding plus a pillow per patch (B1.4 §3.4); the B6.4 energy +with a fitting term (B1.4 §3.5); the certified check of B6.6 with local repair. The dual schemes +of B1.5 were not adopted (B1.6 finding); concave refined regions are closed by refinement (B1.7 +limit), so curvature refines one level only. ---