diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index da83e56..4cd6626 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -16,6 +16,7 @@ pub mod mirror; pub mod nodes; pub mod octree; pub mod planar; +pub mod project; pub mod quad; pub mod quality; pub mod relax; diff --git a/crates/mesh/src/project/lattice.rs b/crates/mesh/src/project/lattice.rs new file mode 100644 index 0000000..cb7d791 --- /dev/null +++ b/crates/mesh/src/project/lattice.rs @@ -0,0 +1,224 @@ +//! The node lattice of a block topology: each block's structured lattice, with the nodes on +//! shared vertices, edges and faces shared between the blocks. +//! +//! Division counts are per class of parallel block edges (opposite edges of a block have the +//! same count, and so every edge chained to them through other blocks): the classes are the +//! connected components of "parallel in some block". A node on a block's boundary is named by +//! the topology vertices it is blended from and its exact (rational) bilinear weights, which do +//! not depend on which block or orientation sees it, so the blocks around a face, edge or vertex +//! find the same node. + +use std::collections::{BTreeMap, HashMap}; + +use super::topology::{CORNERS, FACE_SIDES, Topology, corner_at}; +use crate::{Error, Result}; + +/// The topological entity a node lies inside. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Entity { + Vertex(u32), + /// Ascending vertex pair. + Edge([u32; 2]), + /// Ascending corner vertices. + Face([u32; 4]), + Inside, +} + +/// The ascending corner vertices of block face `(d, s)`. +pub(crate) fn face_key(block: &[u32; 8], d: usize, s: usize) -> [u32; 4] { + let mut k: Vec = CORNERS.iter().zip(block).filter(|(c, _)| c[d] == s).map(|(_, v)| *v).collect(); + k.sort_unstable(); + [0, 1, 2, 3].map(|i| k.get(i).copied().unwrap_or(u32::MAX)) +} + +/// The two corners (block-local) of each block edge along direction `d`. +fn edges_along(d: usize) -> Vec<[usize; 2]> { + CORNERS + .iter() + .enumerate() + .filter(|(_, c)| c[d] == 0) + .map(|(k, c)| { + let mut e = *c; + e[d] = 1; + [k, corner_at(e)] + }) + .collect() +} + +fn find(parent: &mut [usize], mut x: usize) -> usize { + while parent[x] != x { + parent[x] = parent[parent[x]]; + x = parent[x]; + } + x +} + +/// The class of parallel edges of each block direction, and the number of classes. +pub(crate) fn edge_classes(t: &Topology) -> (Vec<[usize; 3]>, usize) { + let mut ids: HashMap<[u32; 2], usize> = HashMap::new(); + let mut id_of = |a: u32, b: u32| { + let k = [a.min(b), a.max(b)]; + let n = ids.len(); + *ids.entry(k).or_insert(n) + }; + let mut block_edges: Vec<[[usize; 4]; 3]> = Vec::with_capacity(t.blocks.len()); + for block in &t.blocks { + let mut per = [[0; 4]; 3]; + for (d, row) in per.iter_mut().enumerate() { + for (slot, [a, b]) in row.iter_mut().zip(edges_along(d)) { + *slot = id_of(block[a], block[b]); + } + } + block_edges.push(per); + } + let mut parent: Vec = (0..ids.len()).collect(); + for per in &block_edges { + for row in per { + for &e in &row[1..] { + let (a, b) = (find(&mut parent, row[0]), find(&mut parent, e)); + if a != b { + parent[a.max(b)] = a.min(b); + } + } + } + } + let mut class: BTreeMap = BTreeMap::new(); + let out = block_edges + .iter() + .map(|per| { + [0, 1, 2].map(|d| { + let r = find(&mut parent, per[d][0]); + let n = class.len(); + *class.entry(r).or_insert(n) + }) + }) + .collect(); + (out, class.len()) +} + +/// The structured lattices of the blocks, joined. +#[derive(Debug, Clone)] +pub(crate) struct Lattice { + pub divisions: Vec<[usize; 3]>, + /// Per block, the global node of each lattice node ([`crate::quad::lattice_index`]). + pub block_nodes: Vec>, + pub entity: Vec, + pub hexes: Vec<[u32; 8]>, + /// The block and cell of each hexahedron. + pub hex_cell: Vec<(u32, [usize; 3])>, +} + +fn gcd(mut a: u64, mut b: u64) -> u64 { + while b != 0 { + (a, b) = (b, a % b); + } + a.max(1) +} + +/// Build the joined lattices; `divisions[b]` are block `b`'s counts per direction (equal across +/// shared edges, [`edge_classes`]). At most `max_elements` hexahedra. +pub(crate) fn build(t: &Topology, divisions: &[[usize; 3]], max_elements: usize) -> Result { + if divisions.len() != t.blocks.len() { + return Err(Error::Invalid("one set of divisions per block".into())); + } + let mut total = 0usize; + for n in divisions { + if n.iter().any(|&x| x == 0 || x > crate::block::MAX_DIVISIONS) { + return Err(Error::Invalid(format!("divisions {n:?} out of range"))); + } + total = total.saturating_add(n[0].saturating_mul(n[1]).saturating_mul(n[2])); + } + if total > max_elements { + return Err(Error::Limit(format!("{total} elements; at most {max_elements} (give a larger size)"))); + } + let mut keys: HashMap, u32> = HashMap::new(); + let mut entity: Vec = Vec::new(); + let mut block_nodes = Vec::with_capacity(t.blocks.len()); + let mut hexes = Vec::with_capacity(total); + let mut hex_cell = Vec::with_capacity(total); + for (b, (block, n)) in t.blocks.iter().zip(divisions).enumerate() { + let size = (n[0] + 1) * (n[1] + 1) * (n[2] + 1); + let mut nodes = Vec::with_capacity(size); + for k in 0..=n[2] { + for j in 0..=n[1] { + for i in 0..=n[0] { + let at = [i, j, k]; + let fixed: Vec<(usize, usize)> = (0..3) + .filter_map(|d| { + if at[d] == 0 { + Some((d, 0)) + } else if at[d] == n[d] { + Some((d, 1)) + } else { + None + } + }) + .collect(); + if fixed.is_empty() { + nodes.push(entity.len() as u32); + entity.push(Entity::Inside); + continue; + } + let free: Vec = (0..3).filter(|d| !fixed.iter().any(|f| f.0 == *d)).collect(); + let den: u64 = free.iter().map(|&d| n[d] as u64).product(); + let mut key: Vec<(u32, u64, u64)> = Vec::new(); + for (c, bits) in CORNERS.iter().enumerate() { + if fixed.iter().any(|&(d, s)| bits[d] != s) { + continue; + } + let num: u64 = free.iter().map(|&d| if bits[d] == 1 { at[d] as u64 } else { (n[d] - at[d]) as u64 }).product(); + if num == 0 { + continue; + } + let g = gcd(num, den); + key.push((block[c], num / g, den / g)); + } + key.sort_unstable(); + let next = entity.len() as u32; + let id = match keys.get(&key) { + Some(&id) => id, + None => { + let mut vs: Vec = key.iter().map(|x| x.0).collect(); + vs.sort_unstable(); + let e = match vs.len() { + 1 => Entity::Vertex(vs[0]), + 2 => Entity::Edge([vs[0], vs[1]]), + _ => Entity::Face([0, 1, 2, 3].map(|q| vs.get(q).copied().unwrap_or(u32::MAX))), + }; + entity.push(e); + keys.insert(key, next); + next + } + }; + nodes.push(id); + } + } + } + let idx = |c: [usize; 3]| crate::quad::lattice_index(*n, c); + for k in 0..n[2] { + for j in 0..n[1] { + for i in 0..n[0] { + let h = CORNERS.map(|c| nodes.get(idx([i + c[0], j + c[1], k + c[2]])).copied().unwrap_or(0)); + hexes.push(h); + hex_cell.push((b as u32, [i, j, k])); + } + } + } + block_nodes.push(nodes); + } + Ok(Lattice { divisions: divisions.to_vec(), block_nodes, entity, hexes, hex_cell }) +} + +/// How many blocks use each block face. +pub(crate) fn face_uses(t: &Topology) -> Result> { + let mut uses: HashMap<[u32; 4], usize> = HashMap::new(); + for block in &t.blocks { + for (d, s) in FACE_SIDES { + *uses.entry(face_key(block, d, s)).or_insert(0) += 1; + } + } + if uses.values().any(|&u| u > 2) { + return Err(Error::Invalid("a block face is shared by more than two blocks".into())); + } + Ok(uses) +} diff --git a/crates/mesh/src/project/mod.rs b/crates/mesh/src/project/mod.rs new file mode 100644 index 0000000..1d89e2e --- /dev/null +++ b/crates/mesh/src/project/mod.rs @@ -0,0 +1,621 @@ +//! Projection blocking: a block topology laid over the part and projected onto its B-rep. +//! +//! Instead of cutting the geometry into sweepable cells, the analyst lays down a clean block +//! topology — a grid with blocks carved away, butterfly (O-grid) templates for round sections, +//! the seven-block sphere — whose outer faces stand for the part's surface. Every block face on +//! the outside of the topology is projected onto a target: one B-rep face, or a union of faces +//! (by default the smooth group of faces nearest it). Where two outer block faces with different +//! targets meet, their edge goes onto the curve where the targets meet, and a vertex on three +//! targets onto their common point ([`target`]). Nodes are then placed outside in: block +//! vertices, edge nodes (evenly by arc length along their curve), face nodes by transfinite +//! interpolation of their edges projected onto the face's target, and interior nodes by the +//! transfinite interpolation of the block's faces ([`crate::quad::transfinite_volume`], Gordon +//! and Hall 1973, Gordon and Thiel 1982). Elliptic smoothing ([`smooth`]) finishes the grid with +//! the boundary held on the geometry. +//! +//! The method follows the public descriptions of block-structured projection meshing (blocks in +//! index space, carved, projected face → surface, edge → curve, vertex → point, then smoothed); +//! the algorithms are the papers cited in each module. + +mod lattice; +pub mod smooth; +pub mod target; +pub mod topology; + +#[cfg(test)] +mod tests; + +use std::collections::{BTreeMap, BTreeSet, HashMap}; + +use fieldcraft_model::ElementKind; + +pub use smooth::{QUALITY_FLOOR, SmoothReport, Smoothing}; +pub use target::Geometry; +pub use topology::Topology; + +use crate::hex::BoundaryFace; +use crate::quality::{Metric, element_quality}; +use crate::surface::{BoundarySurface, FeatureKind}; +use crate::vec3::{P, add, dist, lerp, scale}; +use crate::{Error, Result}; +use lattice::{Entity, build, edge_classes, face_key, face_uses}; +use target::{Constraint, on_target, place}; +use topology::{CORNERS, FACE_SIDES, corner_at, trilinear}; + +/// Faces meeting at less than this angle (degrees) along their shared edge belong to one smooth +/// group, the default target of a block face. +pub const SMOOTH_ANGLE: f64 = 20.0; + +/// The default [`ProjectOptions::min_quality`]: a corner flattened toward 180° (a block corner +/// projected onto a smooth surface) scores far below it. +pub const DEFAULT_MIN_QUALITY: f64 = 0.2; + +/// Sweeps of the inner vertices' displacement blend. +const VERTEX_SWEEPS: usize = 500; + +/// Rounds of arc-length redistribution of an edge's nodes. +const EDGE_ROUNDS: usize = 4; + +/// How to project. +#[derive(Debug, Clone, PartialEq)] +pub struct ProjectOptions { + /// The element size: each class of parallel block edges gets `round(length / size)` + /// divisions (at least one), from its longest projected edge. + pub size: f64, + /// Division counts set by hand: `(block, direction 0..3, count)`; the count applies to the + /// whole class of parallel edges. + pub divisions: Vec<(usize, usize, usize)>, + /// Targets set by hand: `(block, hex8 face, B-rep faces)`. Other outer faces take the smooth + /// group of B-rep faces nearest them. + pub faces: Vec<(usize, u8, Vec)>, + pub smoothing: Smoothing, + pub iterations: usize, + pub max_elements: usize, + /// The smallest scaled Jacobian accepted: a mesh with a worse element is an error that + /// says where a template is missing (0 accepts any valid mesh). + pub min_quality: f64, +} + +impl Default for ProjectOptions { + fn default() -> Self { + ProjectOptions { + size: 1.0, + divisions: Vec::new(), + faces: Vec::new(), + smoothing: Smoothing::Elliptic, + iterations: 100, + max_elements: 2_000_000, + min_quality: DEFAULT_MIN_QUALITY, + } + } +} + +/// How it went. +#[derive(Debug, Clone, PartialEq)] +pub struct ProjectReport { + pub blocks: usize, + pub elements: usize, + /// The smallest scaled Jacobian (corners) of the transfinite grid, before smoothing. + pub initial_min_scaled_jacobian: f64, + /// Certified smallest and mean scaled Jacobian of the result. + pub min_scaled_jacobian: f64, + pub mean_scaled_jacobian: f64, + pub smoothing: SmoothReport, + /// Farthest a boundary node lies from its target (surface, curve or point). + pub max_gap: f64, + pub volume: f64, + pub surface_volume: f64, + /// The B-rep faces each target stands for. + pub targets: Vec>, + /// Divisions per block. + pub divisions: Vec<[usize; 3]>, +} + +/// An all-hex mesh of a projected topology. +#[derive(Debug, Clone, PartialEq)] +pub struct ProjectedMesh { + pub nodes: Vec

, + pub hexes: Vec<[u32; 8]>, + pub boundary: Vec, + pub report: ProjectReport, +} + +/// The face groups joined across smooth B-rep edges. +fn smooth_groups(surface: &BoundarySurface, faces: usize) -> Vec { + let mut parent: Vec = (0..faces).collect(); + fn root(p: &mut [usize], mut x: usize) -> usize { + while p[x] != x { + p[x] = p[p[x]]; + x = p[x]; + } + x + } + let tf = surface.tri_face(); + for e in surface.feature_edges() { + if e.kind != FeatureKind::FaceBoundary || e.angle >= SMOOTH_ANGLE { + continue; + } + let (Some(&a), Some(&b)) = (tf.get(e.triangles[0] as usize), tf.get(e.triangles[1] as usize)) else { continue }; + let (a, b) = (a as usize, b as usize); + if a < faces && b < faces { + let (ra, rb) = (root(&mut parent, a), root(&mut parent, b)); + if ra != rb { + parent[ra.max(rb)] = ra.min(rb); + } + } + } + (0..faces).map(|f| root(&mut parent, f)).collect() +} + +/// The B-rep face a sample `u` of block face `(d, s)` stands for: the nearest face that faces the +/// same way (normals within 60°), else the nearest face. +fn facing(g: &dyn Geometry, faces: &[u32], c: &[P; 8], d: usize, s: usize, u: [f64; 3]) -> Option { + let p = trilinear(c, u); + let (e1, e2) = other(d); + let h = 1e-3; + let diff = |k: usize| { + let (mut a, mut b) = (u, u); + a[k] -= h; + b[k] += h; + crate::vec3::sub(trilinear(c, b), trilinear(c, a)) + }; + let mut n = crate::vec3::cross(diff(e1), diff(e2)); + // Outward: along +d on side 1, against it on side 0. + let inward = { + let mut q = u; + q[d] = if s == 1 { 1.0 - 0.1 } else { 0.1 }; + crate::vec3::sub(trilinear(c, q), p) + }; + if crate::vec3::dot(n, inward) > 0.0 { + n = scale(n, -1.0); + } + let ln = crate::vec3::norm(n); + let mut best: Option<(f64, u32)> = None; + let mut nearest: Option<(f64, u32)> = None; + for &f in faces { + let Some(q) = g.project(f, p) else { continue }; + let dq = dist(p, q); + if nearest.is_none_or(|b| dq < b.0) { + nearest = Some((dq, f)); + } + let agrees = ln > 0.0 && g.normal(f, q).is_some_and(|m| crate::vec3::dot(m, n) > 0.5 * ln); + if agrees && best.is_none_or(|b| dq < b.0) { + best = Some((dq, f)); + } + } + best.or(nearest).map(|b| b.1) +} + +/// The corners of block face `(d, s)` in cyclic order. +fn face_cycle(block: &[u32; 8], d: usize, s: usize) -> [u32; 4] { + let (e1, e2) = other(d); + [[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, b]| { + let mut bits = [0; 3]; + bits[d] = s; + bits[e1] = a; + bits[e2] = b; + block[corner_at(bits)] + }) +} + +pub(crate) fn other(d: usize) -> (usize, usize) { + match d { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + } +} + +fn edge_key(a: u32, b: u32) -> [u32; 2] { + [a.min(b), a.max(b)] +} + +/// Mesh the topology `t` projected onto the body bounded by `surface`, through `g`. +pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, opts: &ProjectOptions) -> Result { + t.check()?; + if !(opts.size > 0.0 && opts.size.is_finite()) { + return Err(Error::Invalid(format!("the element size must be a positive number, got {}", opts.size))); + } + let uses = face_uses(t)?; + let face_count = surface.tri_face().iter().max().map_or(1, |m| *m as usize + 1); + let groups = smooth_groups(surface, face_count); + let diag = t.diagonal().max(opts.size); + let tol = 1e-9 * diag; + + // Targets of the outer block faces. + let explicit: BTreeMap<(usize, u8), Vec> = opts.faces.iter().map(|(b, f, s)| ((*b, *f), s.clone())).collect(); + for ((b, f), set) in &explicit { + let (d, s) = FACE_SIDES + .get(usize::from(*f).wrapping_sub(1)) + .copied() + .ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face (1 to 6)")))?; + let block = t.blocks.get(*b).ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?; + if uses.get(&face_key(block, d, s)) != Some(&1) { + return Err(Error::Invalid(format!("face {f} of block {b} is inside the topology; only outer faces are projected"))); + } + if set.is_empty() || set.iter().any(|&x| x as usize >= face_count) { + return Err(Error::Invalid(format!("face {f} of block {b}: the target must name faces of the body"))); + } + } + let mut targets: Vec> = Vec::new(); + let mut target_of = |set: Vec| -> usize { + match targets.iter().position(|t| *t == set) { + Some(i) => i, + None => { + targets.push(set); + targets.len() - 1 + } + } + }; + let all_faces: Vec = (0..face_count as u32).collect(); + let mut face_target: HashMap<[u32; 4], usize> = HashMap::new(); + let mut outer: Vec<(usize, usize, usize)> = Vec::new(); + for (b, block) in t.blocks.iter().enumerate() { + let c = t.corners(b)?; + for (fi, (d, s)) in FACE_SIDES.iter().copied().enumerate() { + let key = face_key(block, d, s); + if uses.get(&key) != Some(&1) { + continue; + } + outer.push((b, d, s)); + let set = match explicit.get(&(b, fi as u8 + 1)) { + Some(set) => { + let mut set = set.clone(); + set.sort_unstable(); + set.dedup(); + set + } + None => { + // The smooth group most samples of the face lie nearest. + let (e1, e2) = other(d); + let mut votes: BTreeMap = BTreeMap::new(); + for (a, bb) in [(0.5, 0.5), (0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75)] { + let mut u = [0.0; 3]; + u[d] = s as f64; + u[e1] = a; + u[e2] = bb; + if let Some(f) = facing(g, &all_faces, &c, d, s, u) { + *votes.entry(groups.get(f as usize).copied().unwrap_or(0)).or_insert(0) += 1; + } + } + let group = votes + .iter() + .max_by_key(|(gr, n)| (**n, std::cmp::Reverse(**gr))) + .map(|(gr, _)| *gr) + .ok_or_else(|| Error::Failed(format!("face {} of block {b} has no B-rep face to project onto", fi + 1)))?; + (0..face_count as u32).filter(|&f| groups.get(f as usize) == Some(&group)).collect() + } + }; + face_target.insert(key, target_of(set)); + } + } + // Targets of the vertices and edges on the outside. + let mut vertex_on: BTreeMap> = BTreeMap::new(); + let mut edge_on: BTreeMap<[u32; 2], BTreeSet> = BTreeMap::new(); + for &(b, d, s) in &outer { + let block = t.blocks.get(b).copied().unwrap_or([0; 8]); + let tg = face_target.get(&face_key(&block, d, s)).copied().unwrap_or(0); + let cyc = face_cycle(&block, d, s); + for i in 0..4 { + vertex_on.entry(cyc[i]).or_default().insert(tg); + edge_on.entry(edge_key(cyc[i], cyc[(i + 1) % 4])).or_default().insert(tg); + } + } + let as_vec = |s: Option<&BTreeSet>| s.map(|s| s.iter().copied().collect::>()).unwrap_or_default(); + let vertex_c = |v: u32| Constraint::of(&as_vec(vertex_on.get(&v))); + let edge_c = |e: [u32; 2]| Constraint::of(&as_vec(edge_on.get(&e))); + + // 1. Vertices. + let mut vpos: Vec

= t.vertices.clone(); + for (v, p) in vpos.iter_mut().enumerate() { + *p = place(g, &targets, vertex_c(v as u32), *p, tol); + } + // Inner vertices follow: the projected vertices' moves blended through the block edges + // (a bounded Jacobi solve of Laplace's equation for the displacement), so a block the + // projection stretches does not leave its inner corners behind. + let mut nbr: Vec> = vec![BTreeSet::new(); vpos.len()]; + for block in &t.blocks { + for (k, bits) in CORNERS.iter().enumerate() { + for d in 0..3 { + let mut e = *bits; + e[d] = 1 - e[d]; + if let Some(set) = nbr.get_mut(block[k] as usize) { + set.insert(block[corner_at(e)]); + } + } + } + } + let free: Vec = (0..vpos.len()).map(|v| vertex_c(v as u32) == Constraint::Free).collect(); + if free.iter().any(|f| *f) { + let mut disp: Vec

= vpos.iter().zip(&t.vertices).map(|(a, b)| crate::vec3::sub(*a, *b)).collect(); + for _ in 0..VERTEX_SWEEPS { + let mut change = 0.0f64; + for v in 0..disp.len() { + if !free[v] { + continue; + } + let Some(ns) = nbr.get(v).filter(|n| !n.is_empty()) else { continue }; + let mean = scale(ns.iter().fold([0.0; 3], |a, &w| add(a, disp[w as usize])), 1.0 / ns.len() as f64); + change = change.max(dist(mean, disp[v])); + disp[v] = mean; + } + if change <= tol { + break; + } + } + for v in 0..vpos.len() { + if free[v] { + vpos[v] = add(t.vertices[v], disp[v]); + } + } + } + // 2. Divisions per class of parallel edges, from the longest projected edge. + let (classes, nclass) = edge_classes(t); + let mut class_len = vec![0.0f64; nclass]; + for (block, cls) in t.blocks.iter().zip(&classes) { + for d in 0..3 { + for (k, bits) in CORNERS.iter().enumerate() { + if bits[d] != 0 { + continue; + } + let mut e = *bits; + e[d] = 1; + let (a, b) = (block[k], block[corner_at(e)]); + let (pa, pb) = (vpos.get(a as usize).copied().unwrap_or_default(), vpos.get(b as usize).copied().unwrap_or_default()); + let c = edge_c(edge_key(a, b)); + let mut len = 0.0; + let mut prev = pa; + for q in 1..=16 { + let x = if q == 16 { pb } else { place(g, &targets, c, lerp(pa, pb, q as f64 / 16.0), tol) }; + len += dist(prev, x); + prev = x; + } + if let Some(l) = class_len.get_mut(cls[d]) { + *l = l.max(len); + } + } + } + } + let mut class_n: Vec = class_len.iter().map(|l| ((l / opts.size).round() as usize).clamp(1, crate::block::MAX_DIVISIONS)).collect(); + for &(b, d, n) in &opts.divisions { + let cls = classes.get(b).and_then(|c| c.get(d)).ok_or_else(|| Error::Invalid(format!("no direction {d} of block {b}")))?; + if !(1..=crate::block::MAX_DIVISIONS).contains(&n) { + return Err(Error::Invalid(format!("{n} divisions are out of range"))); + } + if let Some(slot) = class_n.get_mut(*cls) { + *slot = n; + } + } + let divisions: Vec<[usize; 3]> = classes.iter().map(|c| c.map(|k| class_n.get(k).copied().unwrap_or(1))).collect(); + let lat = build(t, &divisions, opts.max_elements)?; + + // Each node's constraint. + let constraint: Vec = lat + .entity + .iter() + .map(|e| match e { + Entity::Vertex(v) => vertex_c(*v), + Entity::Edge(k) => edge_c(*k), + Entity::Face(k) => face_target.get(k).map_or(Constraint::Free, |&tg| Constraint::Surface(tg)), + Entity::Inside => Constraint::Free, + }) + .collect(); + let n = lat.entity.len(); + let mut nodes = vec![[0.0; 3]; n]; + let mut placed = vec![false; n]; + for (i, e) in lat.entity.iter().enumerate() { + if let Entity::Vertex(v) = e { + nodes[i] = vpos.get(*v as usize).copied().unwrap_or_default(); + placed[i] = true; + } + } + // 3. Edge nodes: even by arc length along their curve or surface. + for (b, block) in t.blocks.iter().enumerate() { + let nd = lat.divisions[b]; + let bn = &lat.block_nodes[b]; + for d in 0..3 { + for (k, bits) in CORNERS.iter().enumerate() { + if bits[d] != 0 { + continue; + } + let mut e = *bits; + e[d] = 1; + let c = edge_c(edge_key(block[k], block[corner_at(e)])); + let ids: Vec = (0..=nd[d]) + .map(|q| { + let mut at = [bits[0] * nd[0], bits[1] * nd[1], bits[2] * nd[2]]; + at[d] = q; + bn.get(crate::quad::lattice_index(nd, at)).copied().unwrap_or(0) + }) + .collect(); + if ids.iter().skip(1).take(nd[d].saturating_sub(1)).all(|&i| placed[i as usize]) { + continue; + } + let (pa, pb) = (nodes[ids[0] as usize], nodes[ids[nd[d]] as usize]); + let m = nd[d]; + let mut line: Vec

= (0..=m).map(|q| place(g, &targets, c, lerp(pa, pb, q as f64 / m as f64), tol)).collect(); + line[0] = pa; + line[m] = pb; + if c != Constraint::Free { + for _ in 0..EDGE_ROUNDS { + line = even(&line).into_iter().map(|x| place(g, &targets, c, x, tol)).collect(); + line[0] = pa; + line[m] = pb; + } + } + for (q, &id) in ids.iter().enumerate() { + if !placed[id as usize] { + nodes[id as usize] = line[q]; + placed[id as usize] = true; + } + } + } + } + } + // 4. Faces (projected when outside), then 5. the interior. + for (b, block) in t.blocks.iter().enumerate() { + let nd = lat.divisions[b]; + let bn = &lat.block_nodes[b]; + let mut x: Vec> = bn.iter().map(|&i| placed[i as usize].then(|| nodes[i as usize])).collect(); + crate::quad::transfinite_faces(&mut x, nd, None, |d, side, grid| { + if let Some(&tg) = face_target.get(&face_key(block, d, side)) { + for p in grid.iter_mut() { + *p = place(g, &targets, Constraint::Surface(tg), *p, tol); + } + } + Ok(()) + })?; + for (slot, &i) in x.iter_mut().zip(bn) { + if placed[i as usize] { + *slot = Some(nodes[i as usize]); + } else if let Some(p) = *slot { + nodes[i as usize] = p; + placed[i as usize] = true; + } + } + let blend = nd.map(|m| (0..=m).map(|q| q as f64 / m as f64).collect::>()); + crate::quad::transfinite_volume(&mut x, nd, [&blend[0], &blend[1], &blend[2]])?; + for (slot, &i) in x.iter().zip(bn) { + if let (false, Some(p)) = (placed[i as usize], slot) { + nodes[i as usize] = *p; + placed[i as usize] = true; + } + } + } + if placed.iter().any(|p| !p) { + return Err(Error::Failed("a node of the block lattice was not placed".into())); + } + let corners_of = |nodes: &[P], h: &[u32; 8]| h.map(|v| nodes.get(v as usize).copied().unwrap_or_default()); + let initial_min = lat.hexes.iter().map(|h| smooth::hex_scaled_jacobian(&corners_of(&nodes, h))).fold(f64::INFINITY, f64::min); + + // The outer quads (outward, hex8 face order) and the target each lies on. + let hex8_faces = ElementKind::Hex8.faces(); + let mut quads: Vec<([u32; 4], u32, u8, usize)> = Vec::new(); + for (e, (hex, (b, cell))) in lat.hexes.iter().zip(&lat.hex_cell).enumerate() { + let block = t.blocks.get(*b as usize).copied().unwrap_or([0; 8]); + let nd = lat.divisions.get(*b as usize).copied().unwrap_or([1; 3]); + for (fi, (d, s)) in FACE_SIDES.iter().copied().enumerate() { + if cell[d] != if s == 0 { 0 } else { nd[d] - 1 } { + continue; + } + let Some(&tg) = face_target.get(&face_key(&block, d, s)) else { continue }; + let local = hex8_faces.get(fi).copied().unwrap_or(&[]); + let q = [0, 1, 2, 3].map(|i| local.get(i).and_then(|&l| hex.get(l)).copied().unwrap_or(0)); + quads.push((q, e as u32, fi as u8 + 1, tg)); + } + } + let boundary_quads: Vec<[u32; 4]> = quads.iter().map(|q| q.0).collect(); + let smoothing = smooth::smooth( + smooth::Mesh { nodes: &mut nodes, hexes: &lat.hexes, constraint: &constraint, boundary: &boundary_quads }, + g, + &targets, + opts.smoothing, + opts.iterations, + 1e-4 * opts.size, + ); + + // Certified quality, the gap to the geometry and the tags. + let (mut min_sj, mut sum_sj, mut volume) = (f64::INFINITY, 0.0, 0.0); + let mut invalid = 0usize; + for h in &lat.hexes { + let q = element_quality(ElementKind::Hex8, &corners_of(&nodes, h))?; + let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0); + if q.validity().is_none_or(|v| v != crate::quality::Validity::Valid) { + invalid += 1; + } + min_sj = min_sj.min(sj); + sum_sj += sj; + volume += q.size; + } + let mut max_gap = 0.0f64; + for (i, c) in constraint.iter().enumerate() { + if *c != Constraint::Free { + max_gap = max_gap.max(dist(nodes[i], place(g, &targets, *c, nodes[i], tol))); + } + } + let mut boundary = Vec::with_capacity(quads.len()); + let mut enclosed = 0.0; + for (q, e, f, tg) in &quads { + let pts = q.map(|v| nodes.get(v as usize).copied().unwrap_or_default()); + enclosed += crate::hex::face_flux(&pts); + let centre = scale(pts.iter().fold([0.0; 3], |a, p| add(a, *p)), 0.25); + let set = targets.get(*tg).map_or(&[][..], Vec::as_slice); + let brep_face = on_target(g, set, centre).map_or_else(|| set.first().copied().unwrap_or(0), |h| h.2); + boundary.push(BoundaryFace { element: *e, face: *f, brep_face }); + } + let surface_volume = surface.volume(); + let report = ProjectReport { + blocks: t.blocks.len(), + elements: lat.hexes.len(), + initial_min_scaled_jacobian: initial_min, + min_scaled_jacobian: min_sj, + mean_scaled_jacobian: sum_sj / lat.hexes.len().max(1) as f64, + smoothing, + max_gap, + volume, + surface_volume, + targets: targets.clone(), + divisions: lat.divisions.clone(), + }; + // Where the worst element is, for the error messages. + let worst = || { + lat.hexes + .iter() + .zip(&lat.hex_cell) + .min_by(|a, b| smooth::hex_scaled_jacobian(&corners_of(&nodes, a.0)).total_cmp(&smooth::hex_scaled_jacobian(&corners_of(&nodes, b.0)))) + .map_or_else(String::new, |(h, c)| { + let x = scale(corners_of(&nodes, h).iter().fold([0.0; 3], |a, p| add(a, *p)), 0.125); + format!("worst in block {} near ({:.3}, {:.3}, {:.3})", c.0, x[0], x[1], x[2]) + }) + }; + if invalid > 0 { + return Err(Error::Failed(format!( + "{invalid} of {} elements are inverted or degenerate after projection (smallest scaled Jacobian {min_sj:.3}, {}): the blocks do not fit the part — add an O-grid or a template where a block corner lands on a curved face", + lat.hexes.len(), + worst() + ))); + } + if min_sj < opts.min_quality { + return Err(Error::Failed(format!( + "the smallest scaled Jacobian is {min_sj:.3}, below {} ({}): a block corner flattens on a curved face — use an O-grid or a butterfly/sphere template there, or lower `min_quality`", + opts.min_quality, + worst() + ))); + } + if (enclosed - volume).abs() > 1e-6 * volume.abs().max(1e-300) { + return Err(Error::Failed(format!( + "the elements overlap or leave a gap inside (their volume {volume} differs from the {enclosed} their outer faces enclose): check that block faces meet face to face" + ))); + } + if (volume - surface_volume).abs() > 0.2 * surface_volume.abs() { + return Err(Error::Failed(format!( + "the mesh volume {volume} is far from the body's {surface_volume}: the blocks do not cover the part (every outer block face must stand for part of its surface)" + ))); + } + Ok(ProjectedMesh { nodes, hexes: lat.hexes, boundary, report }) +} + +/// The polyline's points redistributed evenly by arc length (ends kept). +fn even(line: &[P]) -> Vec

{ + let m = line.len().saturating_sub(1); + let mut acc = vec![0.0; line.len()]; + for q in 1..line.len() { + acc[q] = acc[q - 1] + dist(line[q - 1], line[q]); + } + let total = acc.last().copied().unwrap_or(0.0); + if !(total > 0.0) || m == 0 { + return line.to_vec(); + } + let mut out = Vec::with_capacity(line.len()); + let mut seg = 0; + for q in 0..=m { + let s = total * q as f64 / m as f64; + while seg + 1 < m && acc[seg + 1] < s { + seg += 1; + } + let l = acc[seg + 1] - acc[seg]; + let t = if l > 0.0 { ((s - acc[seg]) / l).clamp(0.0, 1.0) } else { 0.0 }; + out.push(lerp(line[seg], line[seg + 1], t)); + } + out +} diff --git a/crates/mesh/src/project/smooth.rs b/crates/mesh/src/project/smooth.rs new file mode 100644 index 0000000..dc2ebe5 --- /dev/null +++ b/crates/mesh/src/project/smooth.rs @@ -0,0 +1,401 @@ +//! Smoothing of a projected block mesh with its boundary held on the geometry. +//! +//! Interior nodes with the regular neighbourhood of a structured grid (eight hexahedra, six +//! edge neighbours in three opposite pairs) move by the elliptic (Winslow) grid equations, +//! `Σ g^{ij} x_{ξ^i ξ^j} + Σ g^{ii} φ_i x_{ξ^i} = 0` (A. M. Winslow, *Numerical solution of the +//! quasilinear Poisson equation in a nonuniform triangle mesh*, J. Comput. Phys. 2:149–172, 1967; +//! the three-dimensional form and its discretisation in J. F. Thompson, Z. U. A. Warsi, +//! C. W. Mastin, *Numerical Grid Generation*, North-Holland 1985, ch. 6), solved point by point +//! with central differences in the node's own stencil. The optional control functions `φ_i` +//! follow P. D. Thomas and J. F. Middlecoff, *Direct control of the grid point distribution in +//! meshes generated by elliptic equations*, AIAA J. 18(6):652–656, 1980, in their +//! one-dimensional form `φ = −x_ξ·x_ξξ / |x_ξ|²`, evaluated on the algebraic (transfinite) grid +//! at every node, so the smoothed grid keeps that grid's spacing along each line. Irregular nodes +//! (where three, five or more blocks meet) take the mean of their neighbours (Laplace). +//! +//! Boundary nodes move along their surface or curve: the mean of their boundary neighbours, +//! projected back. Block vertices on three targets stay put. +//! +//! Every move is checked: it is kept only when the smallest scaled Jacobian of the elements +//! around the node stays at or above the lower of its value before and [`QUALITY_FLOOR`], so +//! smoothing never makes a good element bad nor a bad one worse. + +use std::collections::{BTreeMap, BTreeSet}; + +use super::target::{Constraint, Geometry, place}; +use super::topology::{CORNERS, corner_at}; +use crate::vec3::{P, add, cross, dist, dot, norm, scale, sub}; + +/// Smoothing never takes an element's scaled Jacobian below this (unless it already was). +pub const QUALITY_FLOOR: f64 = 0.5; + +/// Rounds of the compass search on the worst elements. +const OPTIMIZE_ROUNDS: usize = 20; + +/// The compass search's directions: the axes and the cube diagonals. +const COMPASS: [P; 14] = { + const S: f64 = 0.577_350_269_189_625_8; + [ + [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], + [S, S, S], + [S, S, -S], + [S, -S, S], + [S, -S, -S], + [-S, S, S], + [-S, S, -S], + [-S, -S, S], + [-S, -S, -S], + ] +}; + +/// Largest magnitude of a control function. +const MAX_CONTROL: f64 = 1.0; + +/// How the interior is smoothed. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Smoothing { + /// The transfinite grid as it is. + None, + /// Each node to the mean of its neighbours. + Laplace, + /// Winslow's elliptic equations. + Elliptic, + /// Winslow's equations with Thomas–Middlecoff control functions from the algebraic grid. + EllipticControlled, +} + +/// The smallest scaled Jacobian over the eight corners of a hexahedron. +pub(crate) fn hex_scaled_jacobian(x: &[P; 8]) -> f64 { + let mut m = f64::INFINITY; + for (k, bits) in CORNERS.iter().enumerate() { + let along = |d: usize| { + let mut b = *bits; + b[d] = 1 - b[d]; + let e = sub(x[corner_at(b)], x[k]); + if bits[d] == 1 { scale(e, -1.0) } else { e } + }; + let (a, b, c) = (along(0), along(1), along(2)); + let l = norm(a) * norm(b) * norm(c); + let v = if l > 0.0 && l.is_finite() { dot(a, cross(b, c)) / l } else { -1.0 }; + m = m.min(v); + } + // Inside: the centre and the six face centres, where a twisted element folds first. + for u in [[0.5, 0.5, 0.5], [0.5, 0.5, 0.0], [0.5, 0.5, 1.0], [0.5, 0.0, 0.5], [0.5, 1.0, 0.5], [0.0, 0.5, 0.5], [1.0, 0.5, 0.5]] { + let der = |d: usize| { + let mut s3 = [0.0; 3]; + for (k, bits) in CORNERS.iter().enumerate() { + let w: f64 = (0..3) + .map(|e| { + if e == d { + if bits[e] == 1 { 1.0 } else { -1.0 } + } else if bits[e] == 1 { + u[e] + } else { + 1.0 - u[e] + } + }) + .product(); + s3 = add(s3, scale(x[k], w)); + } + s3 + }; + let (a, b, c) = (der(0), der(1), der(2)); + let l = norm(a) * norm(b) * norm(c); + m = m.min(if l > 0.0 && l.is_finite() { dot(a, cross(b, c)) / l } else { -1.0 }); + } + m +} + +/// The structured stencil of a regular interior node: `nb[d] = [minus, plus]` and the face +/// diagonals `diag[pair][sa][sb]` for the direction pairs (0, 1), (0, 2), (1, 2). +#[derive(Debug, Clone, Copy)] +struct Stencil { + nb: [[u32; 2]; 3], + diag: [[[u32; 2]; 2]; 3], + control: [f64; 3], +} + +const PAIRS: [(usize, usize); 3] = [(0, 1), (0, 2), (1, 2)]; + +fn stencil(node: u32, around: &[u32], hexes: &[[u32; 8]]) -> Option { + if around.len() != 8 { + return None; + } + // Per hexahedron: the node's three edge neighbours and three face diagonals. + let mut triples = Vec::with_capacity(8); + for &h in around { + let hex = hexes.get(h as usize)?; + let k = hex.iter().position(|&v| v == node)?; + let bits = CORNERS[k]; + let flip = |ds: &[usize]| { + let mut b = bits; + for &d in ds { + b[d] = 1 - b[d]; + } + hex[corner_at(b)] + }; + triples.push(([flip(&[0]), flip(&[1]), flip(&[2])], [flip(&[0, 1]), flip(&[0, 2]), flip(&[1, 2])])); + } + let neighbours: BTreeSet = triples.iter().flat_map(|t| t.0).collect(); + if neighbours.len() != 6 { + return None; + } + let together = |a: u32, b: u32| triples.iter().any(|t| t.0.contains(&a) && t.0.contains(&b)); + let mut axes: Vec<[u32; 2]> = Vec::new(); + let mut seen = BTreeSet::new(); + for &a in &neighbours { + if seen.contains(&a) { + continue; + } + let opposite: Vec = neighbours.iter().copied().filter(|&b| b != a && !together(a, b)).collect(); + let [b] = opposite.as_slice() else { return None }; + seen.insert(a); + seen.insert(*b); + axes.push([a, *b]); + } + let [a0, a1, a2] = axes.as_slice() else { return None }; + let nb = [*a0, *a1, *a2]; + let axis_of = |v: u32| (0..3).find_map(|d| nb[d].iter().position(|&x| x == v).map(|s| (d, s))); + let mut diag = [[[u32::MAX; 2]; 2]; 3]; + for (edge, faces) in &triples { + let dirs: Vec<(usize, usize)> = edge.iter().map(|&v| axis_of(v)).collect::>()?; + // Hex-local direction l went along global axis dirs[l].0, side dirs[l].1. + for (pi, (la, lb)) in PAIRS.iter().enumerate() { + let ((da, sa), (db, sb)) = (dirs[*la], dirs[*lb]); + let (pair, s1, s2) = match PAIRS.iter().position(|p| *p == (da, db)) { + Some(q) => (q, sa, sb), + None => (PAIRS.iter().position(|p| *p == (db, da))?, sb, sa), + }; + diag[pair][s1][s2] = faces[pi]; + } + } + if diag.iter().flatten().flatten().any(|&v| v == u32::MAX) { + return None; + } + Some(Stencil { nb, diag, control: [0.0; 3] }) +} + +fn at(x: &[P], v: u32) -> P { + x.get(v as usize).copied().unwrap_or([0.0; 3]) +} + +/// The Winslow update of a regular node. +fn winslow(x: &[P], s: &Stencil) -> Option

{ + let xi: [P; 3] = [0, 1, 2].map(|d| scale(sub(at(x, s.nb[d][1]), at(x, s.nb[d][0])), 0.5)); + let g = |i: usize, j: usize| dot(xi[i], xi[j]); + let (g11, g22, g33, g12, g13, g23) = (g(0, 0), g(1, 1), g(2, 2), g(0, 1), g(0, 2), g(1, 2)); + let cof_diag = [g22 * g33 - g23 * g23, g11 * g33 - g13 * g13, g11 * g22 - g12 * g12]; + let cof_off = [g13 * g23 - g12 * g33, g12 * g23 - g13 * g22, g13 * g12 - g11 * g23]; + let den = 2.0 * cof_diag.iter().sum::(); + if !(den > 0.0 && den.is_finite()) { + return None; + } + let mut num = [0.0; 3]; + for (d, cd) in cof_diag.iter().enumerate() { + let (m, pl) = (at(x, s.nb[d][0]), at(x, s.nb[d][1])); + num = add(num, scale(add(m, pl), *cd)); + num = add(num, scale(sub(pl, m), 0.5 * cd * s.control[d])); + } + for (q, c) in cof_off.iter().enumerate() { + let dg = &s.diag[q]; + let cross_d = scale(sub(add(at(x, dg[1][1]), at(x, dg[0][0])), add(at(x, dg[1][0]), at(x, dg[0][1]))), 0.25); + num = add(num, scale(cross_d, 2.0 * c)); + } + let out = scale(num, 1.0 / den); + out.iter().all(|v| v.is_finite()).then_some(out) +} + +/// What the smoother reports. +#[derive(Debug, Clone, Copy, PartialEq)] +pub struct SmoothReport { + pub iterations: usize, + /// Moves refused because they would have lowered quality. + pub refused: usize, +} + +/// Everything smoothing needs about the mesh. +pub(crate) struct Mesh<'a> { + pub nodes: &'a mut Vec

, + pub hexes: &'a [[u32; 8]], + pub constraint: &'a [Constraint], + /// Boundary quads (outward), for the neighbours of boundary nodes. + pub boundary: &'a [[u32; 4]], +} + +/// Smooth in place: at most `iterations` sweeps, stopping when no node moves more than `tol`. +pub(crate) fn smooth(m: Mesh, g: &dyn Geometry, targets: &[Vec], method: Smoothing, iterations: usize, tol: f64) -> SmoothReport { + let n = m.nodes.len(); + let mut around: Vec> = vec![Vec::new(); n]; + for (h, hex) in m.hexes.iter().enumerate() { + for &v in hex { + if let Some(a) = around.get_mut(v as usize) { + a.push(h as u32); + } + } + } + // Edge neighbours (all) and boundary neighbours (along boundary quads). + let mut edges: Vec> = vec![BTreeSet::new(); n]; + for hex in m.hexes { + for (k, bits) in CORNERS.iter().enumerate() { + for d in 0..3 { + let mut b = *bits; + b[d] = 1 - b[d]; + if let Some(e) = edges.get_mut(hex[k] as usize) { + e.insert(hex[corner_at(b)]); + } + } + } + } + let mut along_boundary: BTreeMap> = BTreeMap::new(); + for q in m.boundary { + for i in 0..4 { + let (a, b) = (q[i], q[(i + 1) % 4]); + along_boundary.entry(a).or_default().insert(b); + along_boundary.entry(b).or_default().insert(a); + } + } + let cons = |v: u32| m.constraint.get(v as usize).copied().unwrap_or(Constraint::Free); + let mut stencils: Vec> = vec![None; n]; + if matches!(method, Smoothing::Elliptic | Smoothing::EllipticControlled) { + for v in 0..n as u32 { + if cons(v) == Constraint::Free { + let mut s = around.get(v as usize).and_then(|a| stencil(v, a, m.hexes)); + if let (Some(st), Smoothing::EllipticControlled) = (s.as_mut(), method) { + let p = at(m.nodes, v); + for d in 0..3 { + let (mi, pl) = (at(m.nodes, st.nb[d][0]), at(m.nodes, st.nb[d][1])); + let xi = scale(sub(pl, mi), 0.5); + let xii = sub(add(pl, mi), scale(p, 2.0)); + let l2 = dot(xi, xi); + st.control[d] = if l2 > 0.0 { (-dot(xi, xii) / l2).clamp(-MAX_CONTROL, MAX_CONTROL) } else { 0.0 }; + } + } + stencils[v as usize] = s; + } + } + } + let quality = |x: &[P], v: u32| -> f64 { + around.get(v as usize).map_or(1.0, |hs| { + hs.iter().filter_map(|&h| m.hexes.get(h as usize)).map(|hex| hex_scaled_jacobian(&hex.map(|q| at(x, q)))).fold(f64::INFINITY, f64::min) + }) + }; + let mut report = SmoothReport { iterations: 0, refused: 0 }; + if method == Smoothing::None { + return report; + } + let ptol = 1e-3 * tol; + for it in 0..iterations { + report.iterations = it + 1; + let mut moved = 0.0f64; + for v in 0..n as u32 { + let c = cons(v); + let p = at(m.nodes, v); + let mean = |nbrs: &mut dyn Iterator| -> Option

{ + let (mut s, mut k) = ([0.0; 3], 0usize); + for w in nbrs { + s = add(s, at(m.nodes, w)); + k += 1; + } + (k > 0).then(|| scale(s, 1.0 / k as f64)) + }; + let target = match c { + Constraint::Corner(..) => continue, + Constraint::Free => match stencils.get(v as usize).copied().flatten() { + Some(s) => winslow(m.nodes, &s), + None => edges.get(v as usize).and_then(|e| mean(&mut e.iter().copied())), + }, + Constraint::Surface(_) => along_boundary.get(&v).and_then(|e| mean(&mut e.iter().copied())).map(|q| place(g, targets, c, q, ptol)), + Constraint::Curve(a, b) => { + let on_curve = |w: &u32| match cons(*w) { + Constraint::Curve(x, y) => (x, y) == (a, b), + Constraint::Corner(x, y, z) => [x, y, z].contains(&a) && [x, y, z].contains(&b), + _ => false, + }; + along_boundary.get(&v).and_then(|e| mean(&mut e.iter().copied().filter(on_curve))).map(|q| place(g, targets, c, q, ptol)) + } + }; + let Some(q) = target else { continue }; + if dist(q, p) <= 1e-3 * tol { + continue; + } + // The full move, else half, else a quarter of it, whichever keeps the quality. + let before = quality(m.nodes, v); + let mut accepted = false; + for step in [1.0, 0.5, 0.25] { + let x = if step == 1.0 { q } else { place(g, targets, c, add(p, scale(sub(q, p), step)), ptol) }; + if let Some(slot) = m.nodes.get_mut(v as usize) { + *slot = x; + } + if quality(m.nodes, v) >= before.min(QUALITY_FLOOR) { + moved = moved.max(dist(x, p)); + accepted = true; + break; + } + } + if !accepted { + report.refused += 1; + if let Some(slot) = m.nodes.get_mut(v as usize) { + *slot = p; + } + } + } + if moved <= tol { + break; + } + } + // The worst elements: each of their nodes searched for the place that raises the smallest + // scaled Jacobian around it (a compass search, the optimization-based smoothing of L. A. + // Freitag and C. Ollivier-Gooch, *Tetrahedral mesh improvement using swapping and smoothing*, + // IJNME 40:3979–4002, 1997, on Knupp's hexahedral metric), boundary nodes kept on their + // targets. + for _ in 0..OPTIMIZE_ROUNDS { + let poor: BTreeSet = + m.hexes.iter().filter(|h| hex_scaled_jacobian(&h.map(|q| at(m.nodes, q))) < QUALITY_FLOOR).flat_map(|h| h.iter().copied()).collect(); + if poor.is_empty() { + break; + } + let mut improved = false; + for &v in &poor { + let c = cons(v); + if matches!(c, Constraint::Corner(..)) { + continue; + } + let p = at(m.nodes, v); + let reach = edges.get(v as usize).map_or(0.0, |e| e.iter().map(|&w| dist(at(m.nodes, w), p)).fold(f64::INFINITY, f64::min)); + if !(reach > 0.0 && reach.is_finite()) { + continue; + } + let mut best = (quality(m.nodes, v), p); + let mut h = 0.25 * reach; + for _ in 0..6 { + let start = best.1; + for dir in COMPASS { + let x = place(g, targets, c, add(start, scale(dir, h)), ptol); + if let Some(slot) = m.nodes.get_mut(v as usize) { + *slot = x; + } + let qx = quality(m.nodes, v); + if qx > best.0 + 1e-9 { + best = (qx, x); + } + } + if best.1 == start { + h *= 0.5; + } + } + if let Some(slot) = m.nodes.get_mut(v as usize) { + *slot = best.1; + } + improved |= best.1 != p; + } + if !improved { + break; + } + } + report +} diff --git a/crates/mesh/src/project/target.rs b/crates/mesh/src/project/target.rs new file mode 100644 index 0000000..2ad3e5d --- /dev/null +++ b/crates/mesh/src/project/target.rs @@ -0,0 +1,143 @@ +//! Where projected nodes go: onto a target (one B-rep face or a union of faces), onto the curve +//! where two targets meet, or onto the point where three meet. +//! +//! A curve or corner node converges by tangent planes: project the node onto each target, +//! intersect the targets' tangent planes there, move to the point of that intersection nearest +//! the node, repeat. For two transversal surfaces this is Newton's method on the curve; it needs +//! no B-rep edge between the faces, so it tolerates small gaps and overlaps in the CAD (the nodes +//! end on the curve where the two surfaces, extended, meet). Faces meeting at a shallow angle +//! fall back to alternating projections. + +use crate::tet::geometry::Tessellation; +use crate::vec3::{P, add, cross, dist, dot, norm, scale, sub}; + +/// Rounds of the curve and corner iterations. +const ROUNDS: usize = 60; + +/// Below this `|n_a × n_b|` two targets are treated as tangent (sin 3°). +const TANGENT_SIN: f64 = 0.052; + +/// What the projection needs of the geometry. +pub trait Geometry { + /// The closest point on B-rep face `face`. + fn project(&self, face: u32, p: P) -> Option

; + /// The unit normal of B-rep face `face` near `p`. + fn normal(&self, face: u32, p: P) -> Option

; +} + +impl Geometry for Tessellation { + fn project(&self, face: u32, p: P) -> Option

{ + self.project_trimmed(face, p) + } + fn normal(&self, face: u32, p: P) -> Option

{ + self.face_normal(face, p) + } +} + +/// Where a node may go: targets by index into the projection's target list. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub(crate) enum Constraint { + Free, + Surface(usize), + Curve(usize, usize), + Corner(usize, usize, usize), +} + +impl Constraint { + /// The constraint of a node on the targets `on` (distinct, ascending). + pub(crate) fn of(on: &[usize]) -> Constraint { + match *on { + [] => Constraint::Free, + [a] => Constraint::Surface(a), + [a, b] => Constraint::Curve(a, b), + [a, b, c, ..] => Constraint::Corner(a, b, c), + } + } +} + +/// The closest point to `p` on any face of `faces`, with its normal and face. +pub(crate) fn on_target(g: &dyn Geometry, faces: &[u32], p: P) -> Option<(P, P, u32)> { + let mut best: Option<(f64, P, u32)> = None; + for &f in faces { + if let Some(q) = g.project(f, p) { + let d = dist(p, q); + if best.is_none_or(|b| d < b.0) { + best = Some((d, q, f)); + } + } + } + let (_, q, f) = best?; + let n = g.normal(f, q)?; + Some((q, n, f)) +} + +/// Place `p` by its constraint over `targets` (face lists); `tol` is the convergence distance. +/// A point the geometry cannot place stays where it is. +pub(crate) fn place(g: &dyn Geometry, targets: &[Vec], c: Constraint, p: P, tol: f64) -> P { + let set = |i: usize| targets.get(i).map_or(&[][..], Vec::as_slice); + match c { + Constraint::Free => p, + Constraint::Surface(a) => on_target(g, set(a), p).map_or(p, |h| h.0), + Constraint::Curve(a, b) => planes(g, &[set(a), set(b)], p, tol), + Constraint::Corner(a, b, c) => planes(g, &[set(a), set(b), set(c)], p, tol), + } +} + +/// The tangent-plane iteration onto the intersection of two or three targets. +fn planes(g: &dyn Geometry, sets: &[&[u32]], p: P, tol: f64) -> P { + let mut x = p; + for _ in 0..ROUNDS { + let mut hits = Vec::with_capacity(sets.len()); + for s in sets { + match on_target(g, s, x) { + Some(h) => hits.push(h), + None => return x, + } + } + let next = match hits.as_slice() { + [(qa, na, _), (qb, nb, _)] => { + if norm(cross(*na, *nb)) < TANGENT_SIN { + alternate(g, sets, x) + } else { + // The point of both planes nearest x: x − λa·na − λb·nb. + let m = dot(*na, *nb); + let (ra, rb) = (dot(*na, sub(x, *qa)), dot(*nb, sub(x, *qb))); + let den = 1.0 - m * m; + let (la, lb) = ((ra - m * rb) / den, (rb - m * ra) / den); + sub(x, add(scale(*na, la), scale(*nb, lb))) + } + } + [(qa, na, _), (qb, nb, _), (qc, nc, _)] => { + // The three planes' common point (Cramer's rule). + let d = dot(*na, cross(*nb, *nc)); + if d.abs() < TANGENT_SIN * TANGENT_SIN { + alternate(g, sets, x) + } else { + let (ha, hb, hc) = (dot(*na, *qa), dot(*nb, *qb), dot(*nc, *qc)); + scale(add(add(scale(cross(*nb, *nc), ha), scale(cross(*nc, *na), hb)), scale(cross(*na, *nb), hc)), 1.0 / d) + } + } + _ => return x, + }; + if !next.iter().all(|v| v.is_finite()) { + return x; + } + let step = dist(next, x); + x = next; + if step <= tol { + break; + } + } + x +} + +/// One round of alternating projections onto each target. +fn alternate(g: &dyn Geometry, sets: &[&[u32]], p: P) -> P { + let mut x = p; + for s in sets { + if let Some((q, _, _)) = on_target(g, s, x) { + x = q; + } + } + x +} diff --git a/crates/mesh/src/project/tests.rs b/crates/mesh/src/project/tests.rs new file mode 100644 index 0000000..3a9b0b5 --- /dev/null +++ b/crates/mesh/src/project/tests.rs @@ -0,0 +1,228 @@ +use std::f64::consts::{PI, TAU}; + +use super::target::{Constraint, Geometry, place}; +use super::*; +use crate::surface::{BoundarySurface, SurfaceOptions}; +use crate::tet::geometry::{AnalyticSurface, Tessellation}; + +/// A UV sphere of radius `r` about the origin, one face, with its exact surface. +fn sphere(r: f64, n: usize) -> (BoundarySurface, Tessellation) { + 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 = PI * j as f64 / rings as f64 - PI / 2.0; + for i in 0..n { + let a = 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()]; + let s = BoundarySurface::new(&p, &tris, &face, &SurfaceOptions::default()).unwrap(); + let t = Tessellation::new(&s, vec![Some(AnalyticSurface::Sphere { center: [0.0; 3], radius: r })]); + (s, t) +} + +/// A closed cylinder of radius `r` from z = 0 to `h`: faces 0 (bottom), 1 (top), 2 (side). +fn cylinder(r: f64, h: f64, n: usize) -> (BoundarySurface, Tessellation) { + 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 = 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]); + } + let s = BoundarySurface::new(&p, &tris, &face, &SurfaceOptions::default()).unwrap(); + let analytic = vec![ + Some(AnalyticSurface::Plane { normal: [0.0, 0.0, -1.0], point: [0.0; 3] }), + Some(AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, h] }), + Some(AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: r }), + ]; + let t = Tessellation::new(&s, analytic); + (s, t) +} + +fn options(size: f64) -> ProjectOptions { + ProjectOptions { size, ..ProjectOptions::default() } +} + +/// The 7-block sphere: all hex, every boundary node on the sphere, min scaled Jacobian at least +/// 0.5 at every refinement, and the volume converging to the ball's. +#[test] +fn seven_block_sphere_converges_without_losing_quality() { + let r = 10.0; + let (s, t) = sphere(r, 96); + let topo = Topology::sphere([0.0; 3], r, 0.55).unwrap(); + assert_eq!(topo.blocks.len(), 7); + let exact = 4.0 / 3.0 * PI * r * r * r; + let mut last_error = f64::INFINITY; + for size in [4.0, 2.0, 1.0] { + let m = project_mesh(&topo, &s, &t, &options(size)).unwrap(); + let rep = &m.report; + assert!(rep.min_scaled_jacobian >= 0.5, "size {size}: {rep:?}"); + assert!(rep.max_gap < 1e-9 * r, "size {size}: {rep:?}"); + for b in &m.boundary { + let hex = m.hexes[b.element as usize]; + for &l in ElementKind::Hex8.faces()[b.face as usize - 1] { + let p = m.nodes[hex[l] as usize]; + assert!((crate::vec3::norm(p) - r).abs() < 1e-9, "{p:?}"); + } + } + let error = (rep.volume - exact).abs() / exact; + assert!(error < last_error, "size {size}: {error} after {last_error}"); + last_error = error; + } + assert!(last_error < 0.01, "{last_error}"); +} + +/// A disc: the five-block butterfly, its rim on the circle where the side meets the caps. +#[test] +fn butterfly_disc_meshes_with_its_rim_on_the_circle() { + let (r, h) = (10.0, 2.0); + let (s, t) = cylinder(r, h, 128); + let topo = Topology::butterfly([0.0; 3], [0.0, 0.0, 1.0], r, &[0.0, h], 0.55).unwrap(); + assert_eq!(topo.blocks.len(), 5); + let mut last = f64::INFINITY; + for size in [2.0, 1.0, 0.5] { + let m = project_mesh(&topo, &s, &t, &options(size)).unwrap(); + let rep = &m.report; + assert!(rep.min_scaled_jacobian >= 0.5, "size {size}: {rep:?}"); + // Every node on the side lies on the cylinder; the rims on both circles. + let side: Vec<&BoundaryFace> = m.boundary.iter().filter(|b| b.brep_face == 2).collect(); + assert!(!side.is_empty()); + for b in side { + let hex = m.hexes[b.element as usize]; + for &l in ElementKind::Hex8.faces()[b.face as usize - 1] { + let p = m.nodes[hex[l] as usize]; + assert!((p[0].hypot(p[1]) - r).abs() < 1e-9, "{p:?}"); + } + } + let tags: std::collections::BTreeSet = m.boundary.iter().map(|b| b.brep_face).collect(); + assert_eq!(tags.len(), 3); + let error = (rep.volume - PI * r * r * h).abs(); + assert!(error < last, "size {size}"); + last = error; + } +} + +/// A plain grid on a box: corners land on the box's corners, edges on its edges. +#[test] +fn a_grid_on_a_cylinder_is_refused_where_corners_would_flatten() { + // A single block projected onto a cylinder puts its four side edges on the circle: the + // block's corners flatten to 180°, which the mesher reports instead of meshing. + let (s, t) = cylinder(5.0, 4.0, 64); + let topo = Topology::grid([&[-3.0, 3.0], &[-3.0, 3.0], &[0.0, 4.0]]).unwrap(); + let e = project_mesh(&topo, &s, &t, &options(1.0)).unwrap_err().to_string(); + assert!(e.contains("O-grid") || e.contains("inverted"), "{e}"); +} + +/// A geometry of exact surfaces, for the curve and corner iterations. +struct Exact(Vec); + +impl Geometry for Exact { + fn project(&self, face: u32, p: P) -> Option

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

{ + self.0.get(face as usize)?.normal_at(p) + } +} + +/// Nodes on two targets converge to their intersection curve; on three, to the common point. +#[test] +fn curve_and_corner_nodes_converge_by_tangent_planes() { + let g = Exact(vec![ + AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [1.0, 0.0, 0.0], radius: 10.0 }, + AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: 6.0 }, + AnalyticSurface::Plane { normal: [0.0, 1.0, 0.0], point: [0.0; 3] }, + ]); + let targets = vec![vec![0], vec![1], vec![2]]; + for start in [[3.0, 3.0, 9.0], [-5.0, 1.0, 7.0], [1.0, -4.0, 11.0]] { + let p = place(&g, &targets, Constraint::Curve(0, 1), start, 1e-13); + assert!((p[1].hypot(p[2]) - 10.0).abs() < 1e-9 && (p[0].hypot(p[1]) - 6.0).abs() < 1e-9, "{p:?}"); + } + let p = place(&g, &targets, Constraint::Corner(0, 1, 2), [5.0, 1.0, 9.0], 1e-13); + assert!((p[0] - 6.0).abs() < 1e-9 && p[1].abs() < 1e-9 && (p[2] - 10.0).abs() < 1e-9, "{p:?}"); +} + +#[test] +fn o_grid_keeps_neighbours_and_opens_faces() { + let mut t = Topology::grid([&[0.0, 1.0, 2.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap(); + let on = t.o_grid(0, &[2], [0.5; 3]).unwrap(); + assert_eq!(on[0].len(), 5); + // Inner block + five side blocks (top open), the neighbour untouched. + assert_eq!(t.blocks.len(), 7); + t.check().unwrap(); + let uses = lattice::face_uses(&t).unwrap(); + // The top of the old block is now five faces, each used once. + // Inner block to its five side blocks, the four ring blocks to the bottom one and to each + // other, and the untouched neighbour. + assert_eq!(uses.values().filter(|&&u| u == 2).count(), 5 + 4 + 4 + 1, "{uses:?}"); + let made = t.extrude(&on[0], [0.0, 0.0, 1.0], 2).unwrap(); + assert_eq!(made.len(), 10); + t.check().unwrap(); + let face = t.face_towards(made[8], [0.0, 0.0, 1.0]).unwrap(); + assert_eq!(face, 2); + let up = t.extrude(&[(made[8], 2)], [0.0, 0.0, 1.0], 1).unwrap(); + let side = t.face_towards(1, [1.0, 0.0, 0.0]).unwrap(); + t.connect((up[0], 2), (1, side), 2).unwrap(); + t.remove(&[1]).unwrap(); + t.check().unwrap(); +} + +#[test] +fn hostile_topologies_are_errors() { + let (s, g) = sphere(1.0, 32); + let cube = Topology::grid([&[0.0, 1.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap(); + let mut inside_out = cube.clone(); + inside_out.blocks[0].swap(1, 3); + let mut repeated = cube.clone(); + repeated.blocks[0][1] = repeated.blocks[0][0]; + let mut missing = cube.clone(); + missing.blocks[0][7] = 99; + let mut nan = cube.clone(); + nan.vertices[0] = [f64::NAN, 0.0, 0.0]; + for t in [Topology::default(), inside_out, repeated, missing, nan] { + assert!(project_mesh(&t, &s, &g, &options(0.5)).is_err(), "{t:?}"); + } + for size in [0.0, -1.0, f64::NAN, f64::INFINITY] { + assert!(project_mesh(&cube, &s, &g, &options(size)).is_err()); + } + assert!(project_mesh(&cube, &s, &g, &ProjectOptions { size: 1e-6, max_elements: 1000, ..ProjectOptions::default() }).is_err()); + let bad_face = ProjectOptions { faces: vec![(0, 9, vec![0])], ..options(0.5) }; + assert!(project_mesh(&cube, &s, &g, &bad_face).is_err()); + let bad_target = ProjectOptions { faces: vec![(0, 1, vec![7])], ..options(0.5) }; + assert!(project_mesh(&cube, &s, &g, &bad_target).is_err()); + let mut t = cube.clone(); + assert!(t.o_grid(0, &[], [1.5; 3]).is_err() && t.o_grid(0, &[7], [0.5; 3]).is_err() && t.o_grid(3, &[], [0.5; 3]).is_err()); + assert!(t.extrude(&[(0, 0)], [0.0, 0.0, 1.0], 1).is_err() && t.extrude(&[(0, 1)], [f64::NAN; 3], 1).is_err()); + assert!(Topology::grid([&[0.0], &[0.0, 1.0], &[0.0, 1.0]]).is_err()); + assert!(Topology::grid([&[1.0, 0.0], &[0.0, 1.0], &[0.0, 1.0]]).is_err()); + assert!(Topology::sphere([0.0; 3], -1.0, 0.5).is_err() && Topology::butterfly([0.0; 3], [0.0; 3], 1.0, &[0.0, 1.0], 0.5).is_err()); + let mut g3 = Topology::grid([&[0.0, 1.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap(); + assert!(g3.carve([1, 1, 1], &[[[0, 0], [0, 0], [0, 0]]]).is_err()); +} + diff --git a/crates/mesh/src/project/topology.rs b/crates/mesh/src/project/topology.rs new file mode 100644 index 0000000..e24a3b3 --- /dev/null +++ b/crates/mesh/src/project/topology.rs @@ -0,0 +1,514 @@ +//! The block topology of the projection method: hexahedral blocks over shared vertices, built +//! from a carved grid, from the round-feature templates (butterfly, seven-block sphere) and by the +//! O-grid and extrude operations. +//! +//! A block's corners are in hex8 order (1-2-3-4 counter-clockwise seen from 5-6-7-8's side) and +//! its parameter directions run 1→2, 1→4 and 1→5. Faces are numbered as hex8 faces: 1 = 1-2-3-4, +//! 2 = 5-6-7-8, 3 = 1-2-6-5, 4 = 2-3-7-6, 5 = 3-4-8-7, 6 = 4-1-5-8. Positions are only a first +//! guess: the projection moves every boundary vertex onto the geometry. + +use std::collections::BTreeMap; + +use crate::vec3::{P, add, det, dist, scale, sub}; +use crate::{Error, Result}; + +/// Most blocks one topology holds. +pub const MAX_BLOCKS: usize = 20_000; + +/// Corner `c` of a block in parameter space (hex8 order). +pub(crate) const CORNERS: [[usize; 3]; 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]]; + +/// The direction and side of each hex8 face (1-based face `f` at index `f - 1`). +pub(crate) const FACE_SIDES: [(usize, usize); 6] = [(2, 0), (2, 1), (1, 0), (0, 1), (1, 1), (0, 0)]; + +/// The corner at parameter `c`. +pub(crate) fn corner_at(c: [usize; 3]) -> usize { + CORNERS.iter().position(|x| *x == c).unwrap_or(0) +} + +/// Hex8 face number (1-based) of direction `d`, side `s`. +pub(crate) fn face_number(d: usize, s: usize) -> u8 { + FACE_SIDES.iter().position(|x| *x == (d, s)).map_or(1, |i| i as u8 + 1) +} + +/// Blocks over shared vertices. +#[derive(Debug, Clone, Default, PartialEq)] +pub struct Topology { + pub vertices: Vec

, + pub blocks: Vec<[u32; 8]>, +} + +/// The trilinear map of a block's corners at parameter `u`. +pub(crate) fn trilinear(c: &[P; 8], u: [f64; 3]) -> P { + let mut x = [0.0; 3]; + for (k, bits) in CORNERS.iter().enumerate() { + let w: f64 = (0..3).map(|d| if bits[d] == 1 { u[d] } else { 1.0 - u[d] }).product(); + x = add(x, scale(c[k], w)); + } + x +} + +/// The smallest corner Jacobian determinant of a trilinear block (positive for a valid one). +pub(crate) fn min_corner_det(c: &[P; 8]) -> f64 { + let mut m = f64::INFINITY; + for (k, bits) in CORNERS.iter().enumerate() { + // The three edges leaving corner k along increasing parameters, signed to keep the + // right-handed order. + let along = |d: usize| { + let mut b = *bits; + b[d] = 1 - b[d]; + let e = sub(c[corner_at(b)], c[k]); + if bits[d] == 1 { scale(e, -1.0) } else { e } + }; + m = m.min(det(along(0), along(1), along(2))); + } + m +} + +impl Topology { + /// The corner positions of block `b`. + pub fn corners(&self, b: usize) -> Result<[P; 8]> { + let block = self.blocks.get(b).ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?; + let mut out = [[0.0; 3]; 8]; + for (o, &v) in out.iter_mut().zip(block) { + *o = *self.vertices.get(v as usize).ok_or_else(|| Error::Invalid(format!("block {b} uses vertex {v}, which does not exist")))?; + } + Ok(out) + } + + /// Every block uses existing, distinct vertices, is right-handed and caps are kept. + pub fn check(&self) -> Result<()> { + if self.blocks.is_empty() { + return Err(Error::Invalid("the topology has no blocks".into())); + } + if self.blocks.len() > MAX_BLOCKS { + return Err(Error::Limit(format!("{} blocks; at most {MAX_BLOCKS}", self.blocks.len()))); + } + if let Some(v) = self.vertices.iter().position(|v| !v.iter().all(|x| x.is_finite() && x.abs() <= fieldcraft_model::MAX_COORD)) { + return Err(Error::Invalid(format!("vertex {v} is not a finite point"))); + } + for (b, block) in self.blocks.iter().enumerate() { + let mut sorted = *block; + sorted.sort_unstable(); + if sorted.windows(2).any(|w| w[0] == w[1]) { + return Err(Error::Invalid(format!("block {b} uses one vertex twice"))); + } + if !(min_corner_det(&self.corners(b)?) > 0.0) { + return Err(Error::Invalid(format!( + "block {b} is inside out or folded at a corner: give its corners in hex8 order (1-2-3-4 counter-clockwise seen from 5-6-7-8)" + ))); + } + } + Ok(()) + } + + /// The bounding-box diagonal of the vertices. + pub fn diagonal(&self) -> f64 { + let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]); + for v in &self.vertices { + for d in 0..3 { + lo[d] = lo[d].min(v[d]); + hi[d] = hi[d].max(v[d]); + } + } + if self.vertices.is_empty() { 0.0 } else { dist(lo, hi) } + } + + /// A grid of blocks between the planes `planes[d]` (increasing coordinates along X, Y and Z), + /// block `(i, j, k)` at index `(k·nj + j)·ni + i`. + pub fn grid(planes: [&[f64]; 3]) -> Result { + for (d, p) in planes.iter().enumerate() { + if p.len() < 2 || !p.iter().all(|x| x.is_finite()) || p.windows(2).any(|w| !(w[1] > w[0])) { + return Err(Error::Invalid(format!("the grid planes along {} must be at least two increasing numbers", ["X", "Y", "Z"][d]))); + } + } + let n = planes.map(|p| p.len() - 1); + if n.iter().try_fold(1usize, |a, &b| a.checked_mul(b)).is_none_or(|c| c > MAX_BLOCKS) { + return Err(Error::Limit(format!("a grid of {n:?} blocks; at most {MAX_BLOCKS}"))); + } + let vid = |i: usize, j: usize, k: usize| ((k * (n[1] + 1) + j) * (n[0] + 1) + i) as u32; + let mut vertices = Vec::new(); + for z in planes[2] { + for y in planes[1] { + for x in planes[0] { + vertices.push([*x, *y, *z]); + } + } + } + let mut blocks = Vec::new(); + for k in 0..n[2] { + for j in 0..n[1] { + for i in 0..n[0] { + blocks.push(CORNERS.map(|c| vid(i + c[0], j + c[1], k + c[2]))); + } + } + } + Ok(Topology { vertices, blocks }) + } + + /// Remove the blocks of a grid of `n` blocks per direction whose index lies in every range + /// of `ranges` (inclusive, per direction), keeping the others in order; vertices no block + /// uses are dropped. + pub fn carve(&mut self, n: [usize; 3], ranges: &[[[usize; 2]; 3]]) -> Result<()> { + if self.blocks.len() != n.iter().product::() { + return Err(Error::Invalid("carving needs the grid's own block numbering".into())); + } + let mut keep = Vec::new(); + for (b, block) in self.blocks.iter().enumerate() { + let at = [b % n[0], (b / n[0]) % n[1], b / (n[0] * n[1])]; + let gone = ranges.iter().any(|r| (0..3).all(|d| r[d][0] <= at[d] && at[d] <= r[d][1])); + if !gone { + keep.push(*block); + } + } + self.blocks = keep; + self.compact(); + if self.blocks.is_empty() { Err(Error::Invalid("carving removed every block".into())) } else { Ok(()) } + } + + /// Drop unused vertices and renumber. + fn compact(&mut self) { + let mut map: BTreeMap = BTreeMap::new(); + let mut vertices = Vec::new(); + for block in &mut self.blocks { + for v in block.iter_mut() { + let next = map.len() as u32; + let new = *map.entry(*v).or_insert_with(|| { + vertices.push(self.vertices.get(*v as usize).copied().unwrap_or([0.0; 3])); + next + }); + *v = new; + } + } + self.vertices = vertices; + } + + /// Join vertices closer than `tol` (the first one stays). + pub fn merge(&mut self, tol: f64) { + let mut map: Vec = (0..self.vertices.len() as u32).collect(); + // Sorted by X, so only a window is compared. + let mut order: Vec = (0..self.vertices.len()).collect(); + order.sort_by(|&a, &b| self.vertices[a][0].total_cmp(&self.vertices[b][0]).then(a.cmp(&b))); + for (oi, &a) in order.iter().enumerate() { + if map[a] as usize != a { + continue; + } + for &b in order.iter().skip(oi + 1) { + if self.vertices[b][0] - self.vertices[a][0] > tol { + break; + } + if map[b] as usize == b && dist(self.vertices[a], self.vertices[b]) <= tol { + map[b] = a as u32; + } + } + } + for block in &mut self.blocks { + for v in block.iter_mut() { + *v = map.get(*v as usize).copied().unwrap_or(*v); + } + } + self.compact(); + } + + /// Replace block `b` by an O-grid: an inner block, scaled by `inset[d]` (0 to 1) about the + /// block's centre along each of its directions (from the open face, along a direction with one + /// open face), and one block between each face of the inner + /// block and the block's own face, except across the `open` faces (hex8 face numbers), where + /// the inner block reaches the face. Faces of `b` stay faces of the result, so neighbours + /// still fit, except the open ones, which become the inner block's face ringed by the side + /// blocks' faces (the butterfly pattern on that face). The inner block takes index `b`; the + /// others are appended. Returns, per open face, the faces now on it: the inner block's first, + /// then the ring's. + pub fn o_grid(&mut self, b: usize, open: &[u8], inset: [f64; 3]) -> Result>> { + if !inset.iter().all(|i| *i > 0.0 && *i < 1.0) { + return Err(Error::Invalid(format!("the O-grid inset must be between 0 and 1, got {inset:?}"))); + } + if let Some(f) = open.iter().find(|f| !(1..=6).contains(*f)) { + return Err(Error::Invalid(format!("face {f} is not a hex8 face (1 to 6)"))); + } + let c = self.corners(b)?; + let outer = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?; + let is_open = |d: usize, s: usize| open.contains(&face_number(d, s)); + // Each corner's inner copy: on an open face it stays on the face. + let mut inner = [0u32; 8]; + for (k, bits) in CORNERS.iter().enumerate() { + // Along a direction with one open face the inner block reaches `inset` of the way + // across from it; with none it is centred. + let u = [0, 1, 2].map(|d| { + let b = bits[d] as f64; + match (is_open(d, 0), is_open(d, 1)) { + (true, true) => b, + (true, false) => { + if bits[d] == 0 { + 0.0 + } else { + inset[d] + } + } + (false, true) => { + if bits[d] == 1 { + 1.0 + } else { + 1.0 - inset[d] + } + } + (false, false) => 0.5 + inset[d] * (b - 0.5), + } + }); + inner[k] = self.vertices.len() as u32; + self.vertices.push(trilinear(&c, u)); + } + if let Some(slot) = self.blocks.get_mut(b) { + *slot = inner; + } + let mut on_open: Vec> = open.iter().map(|&f| vec![(b, f)]).collect(); + for (d, s) in FACE_SIDES { + if is_open(d, s) { + continue; + } + // Local directions: the face's two directions, then from the outer face inwards. + let (e1, e2) = match d { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + let pick = |a: usize, bb: usize, w: usize| { + let mut bits = [0; 3]; + bits[d] = s; + bits[e1] = a; + bits[e2] = bb; + let k = corner_at(bits); + if w == 0 { outer[k] } else { inner[k] } + }; + let mut block = CORNERS.map(|q| pick(q[0], q[1], q[2])); + let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3])); + let flipped = min_corner_det(&pos(&block)) < 0.0; + if flipped { + block = CORNERS.map(|q| pick(q[1], q[0], q[2])); + } + let index = self.blocks.len(); + self.blocks.push(block); + // This block's face on each open face of the old block. + for (list, &f) in on_open.iter_mut().zip(open) { + let Some(&(od, os)) = FACE_SIDES.get(usize::from(f) - 1) else { continue }; + let local = if od == e1 { + 0 + } else if od == e2 { + 1 + } else { + continue; + }; + let local = if flipped { 1 - local } else { local }; + list.push((index, face_number(local, os))); + } + } + Ok(on_open) + } + + /// The hex8 face of block `b` whose outward direction is nearest `dir`. + pub fn face_towards(&self, b: usize, dir: P) -> Result { + let c = self.corners(b)?; + let centre = scale(c.iter().fold([0.0; 3], |a, p| add(a, *p)), 0.125); + let mut best = (f64::NEG_INFINITY, 1u8); + for (fi, (d, s)) in FACE_SIDES.iter().enumerate() { + let fc = CORNERS.iter().zip(&c).filter(|(bits, _)| bits[*d] == *s).fold([0.0; 3], |a, (_, p)| add(a, *p)); + let out = sub(scale(fc, 0.25), centre); + let l = crate::vec3::norm(out); + let score = if l > 0.0 { crate::vec3::dot(out, dir) / l } else { f64::NEG_INFINITY }; + if score > best.0 { + best = (score, fi as u8 + 1); + } + } + Ok(best.1) + } + + /// The corners of a block face in cyclic order. + fn face_quad(&self, b: usize, f: u8) -> Result<[u32; 4]> { + let (d, s) = FACE_SIDES.get(usize::from(f).wrapping_sub(1)).copied().ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face")))?; + let block = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?; + let (e1, e2) = match d { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + Ok([[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, bb]| { + let mut bits = [0; 3]; + bits[d] = s; + bits[e1] = a; + bits[e2] = bb; + block[corner_at(bits)] + })) + } + + /// Fill the space between two block faces with `layers` blocks, each face's corner matched + /// to the nearest of the other's (over the eight ways of laying one quad on the other). + /// Returns the new blocks. + pub fn connect(&mut self, from: (usize, u8), to: (usize, u8), layers: usize) -> Result> { + if !(1..=1000).contains(&layers) { + return Err(Error::Invalid(format!("a connection has 1 to 1000 layers, got {layers}"))); + } + let a = self.face_quad(from.0, from.1)?; + let b = self.face_quad(to.0, to.1)?; + let verts = self.vertices.clone(); + let at = |v: u32| verts.get(v as usize).copied().unwrap_or([0.0; 3]); + let mut best: Option<(f64, [u32; 4])> = None; + for rev in [false, true] { + for shift in 0..4 { + let m: [u32; 4] = [0, 1, 2, 3].map(|i| { + let k = if rev { (4 + shift - i) % 4 } else { (shift + i) % 4 }; + b[k] + }); + let cost: f64 = (0..4).map(|i| dist(at(a[i]), at(m[i]))).sum(); + if best.is_none_or(|x| cost < x.0) { + best = Some((cost, m)); + } + } + } + let (_, m) = best.ok_or_else(|| Error::Invalid("no way to match the faces".into()))?; + // Intermediate layers of vertices, then the blocks. + let mut rows: Vec<[u32; 4]> = vec![a]; + for l in 1..layers { + let t = l as f64 / layers as f64; + let pts = [0, 1, 2, 3].map(|i| crate::vec3::lerp(at(a[i]), at(m[i]), t)); + let row = pts.map(|p| { + self.vertices.push(p); + self.vertices.len() as u32 - 1 + }); + rows.push(row); + } + rows.push(m); + let mut made = Vec::new(); + for w in rows.windows(2) { + let (q, t) = (w[0], w[1]); + let mut block = [q[0], q[1], q[2], q[3], t[0], t[1], t[2], t[3]]; + let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3])); + if min_corner_det(&pos(&block)) < 0.0 { + block = [q[0], q[3], q[2], q[1], t[0], t[3], t[2], t[1]]; + } + made.push(self.blocks.len()); + self.blocks.push(block); + } + Ok(made) + } + + /// Remove blocks (the others keep their order and are renumbered). + pub fn remove(&mut self, blocks: &[usize]) -> Result<()> { + if let Some(b) = blocks.iter().find(|b| **b >= self.blocks.len()) { + return Err(Error::Invalid(format!("there is no block {b}"))); + } + let mut k = 0; + self.blocks.retain(|_| { + k += 1; + !blocks.contains(&(k - 1)) + }); + self.compact(); + if self.blocks.is_empty() { Err(Error::Invalid("every block was removed".into())) } else { Ok(()) } + } + + /// Extrude block faces `(block, hex8 face)` by `offset` in `layers` blocks each; the faces' + /// shared vertices stay shared. Returns the new blocks' indices, the last layer's outer face + /// is face 2 of the blocks made from it. + pub fn extrude(&mut self, faces: &[(usize, u8)], offset: P, layers: usize) -> Result> { + if !(1..=1000).contains(&layers) { + return Err(Error::Invalid(format!("an extrusion has 1 to 1000 layers, got {layers}"))); + } + if !offset.iter().all(|x| x.is_finite()) { + return Err(Error::Invalid("the extrusion offset must be finite".into())); + } + let mut quads = Vec::new(); + for &(b, f) in faces { + let (d, s) = + FACE_SIDES.get(usize::from(f).wrapping_sub(1)).copied().ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face")))?; + let block = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?; + let (e1, e2) = match d { + 0 => (1, 2), + 1 => (0, 2), + _ => (0, 1), + }; + let q = [[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, bb]| { + let mut bits = [0; 3]; + bits[d] = s; + bits[e1] = a; + bits[e2] = bb; + block[corner_at(bits)] + }); + quads.push(q); + } + let mut made = Vec::new(); + let mut layer_of: BTreeMap = BTreeMap::new(); + let mut base: Vec<[u32; 4]> = quads; + for _ in 0..layers { + layer_of.clear(); + let mut next = Vec::new(); + for q in &base { + let mut top = [0u32; 4]; + for (t, &v) in top.iter_mut().zip(q) { + *t = match layer_of.get(&v) { + Some(&w) => w, + None => { + let p = self.vertices.get(v as usize).copied().ok_or_else(|| Error::Invalid(format!("vertex {v} does not exist")))?; + let w = self.vertices.len() as u32; + self.vertices.push(add(p, scale(offset, 1.0 / layers as f64))); + layer_of.insert(v, w); + w + } + }; + } + let mut block = [q[0], q[1], q[2], q[3], top[0], top[1], top[2], top[3]]; + let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3])); + if min_corner_det(&pos(&block)) < 0.0 { + block = [q[0], q[3], q[2], q[1], top[0], top[3], top[2], top[1]]; + } + made.push(self.blocks.len()); + self.blocks.push(block); + next.push(top); + } + base = next; + } + Ok(made) + } + + /// The seven-block sphere: an inner cube with its corners at `inset · radius` from the + /// centre and six blocks between the cube's faces and the sphere (the cube-to-sphere O-grid). + pub fn sphere(center: P, radius: f64, inset: f64) -> Result { + if !(radius > 0.0 && radius.is_finite()) || !center.iter().all(|x| x.is_finite()) { + return Err(Error::Invalid(format!("a sphere needs a finite centre and a positive radius, got {radius}"))); + } + let h = radius / 3f64.sqrt(); + let mut t = Topology::grid([&[center[0] - h, center[0] + h], &[center[1] - h, center[1] + h], &[center[2] - h, center[2] + h]])?; + t.o_grid(0, &[], [inset; 3])?; + Ok(t) + } + + /// The butterfly (five-block O-grid of a disc) of radius `radius` about `center` in the plane + /// normal to `axis`, stacked through `levels` (distances along `axis` from `center`, + /// increasing): a centre square with its corners at `inset · radius` and four blocks around + /// it out to the circle, per level. + pub fn butterfly(center: P, axis: P, radius: f64, levels: &[f64], inset: f64) -> Result { + if !(radius > 0.0 && radius.is_finite()) || !center.iter().all(|x| x.is_finite()) { + return Err(Error::Invalid(format!("a butterfly needs a finite centre and a positive radius, got {radius}"))); + } + let z = crate::vec3::norm(axis); + if !(z > 0.0 && z.is_finite()) { + return Err(Error::Invalid("the butterfly's axis must be a nonzero direction".into())); + } + let w = scale(axis, 1.0 / z); + let t0 = if w[0].abs() < 0.9 { [1.0, 0.0, 0.0] } else { [0.0, 1.0, 0.0] }; + let u = { + let c = crate::vec3::cross(w, t0); + scale(c, 1.0 / crate::vec3::norm(c)) + }; + let v = crate::vec3::cross(w, u); + let h = radius / 2f64.sqrt(); + let grid = Topology::grid([&[-h, h], &[-h, h], levels])?; + let mut t = Topology { + vertices: grid.vertices.iter().map(|p| add(center, add(add(scale(u, p[0]), scale(v, p[1])), scale(w, p[2])))).collect(), + blocks: grid.blocks, + }; + for b in 0..levels.len() - 1 { + t.o_grid(b, &[1, 2], [inset; 3])?; + } + t.merge(1e-9 * (radius + levels.iter().fold(0.0f64, |m, x| m.max(x.abs())))); + Ok(t) + } +} diff --git a/crates/mesh/src/tet/geometry.rs b/crates/mesh/src/tet/geometry.rs index d05bc21..e2a54e7 100644 --- a/crates/mesh/src/tet/geometry.rs +++ b/crates/mesh/src/tet/geometry.rs @@ -496,6 +496,53 @@ impl Tessellation { } } +impl Tessellation { + /// The closest point to `p` on the trimmed face `face`: the exact surface's closest point + /// when it lies on the face (within a fifth of the local triangle size of the face's + /// triangles), else the closest point of the face's triangles. Unlike + /// [`BoundaryGeometry::project_face`] it never returns a point of the surface's extension + /// beyond the face's boundary, whatever the distance from `p`. + pub fn project_trimmed(&self, face: u32, p: P) -> Option

{ + let accept = |t: u32| self.tri_face.get(t as usize) == Some(&face); + let nearest = |x: P| { + let dist = |t: u32| { + let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]); + let (a, b, c) = (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2])); + let q = closest_on_triangle(x, a, b, c); + let d = sub(q, x); + let mut e = [measure::dist(a, b), measure::dist(b, c), measure::dist(c, a)]; + e.sort_by(f64::total_cmp); + (dot(d, d), [q[0], q[1], q[2], e[1]]) + }; + self.tri_tree.nearest4(x, &accept, &dist).map(|r| ([r[0], r[1], r[2]], r[3])) + }; + let (near, _) = nearest(p)?; + let Some(q) = self.analytic(face).and_then(|s| s.project(p)) else { return Some(near) }; + let (back, size) = nearest(q)?; + Some(if measure::dist(back, q) <= 0.2 * size { q } else { near }) + } + + /// The outward unit normal of B-rep face `face` near `p`: the exact surface's (oriented like + /// the nearest triangle) when known, else the nearest triangle's. + pub fn face_normal(&self, face: u32, p: P) -> Option

{ + let accept = |t: u32| self.tri_face.get(t as usize) == Some(&face); + let dist = |t: u32| { + let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]); + let (a, b, c) = (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2])); + let d = sub(closest_on_triangle(p, a, b, c), p); + let n = unit(cross(sub(b, a), sub(c, a))).unwrap_or([0.0; 3]); + (dot(d, d), [n[0], n[1], n[2], 0.0]) + }; + let r = self.tri_tree.nearest4(p, &accept, &dist)?; + let flat = unit([r[0], r[1], r[2]])?; + match self.analytic(face).and_then(|s| s.normal_at(p)) { + Some(n) if dot(n, flat) < 0.0 => Some(scale(n, -1.0)), + Some(n) => Some(n), + None => Some(flat), + } + } +} + impl BoundaryGeometry for Tessellation { fn project_face(&self, face: u32, p: P) -> Option

{ let (near, scale_len) = self.on_tessellation(face, p)?;