diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 1d79ff1..da83e56 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -18,6 +18,7 @@ pub mod octree; pub mod planar; pub mod quad; pub mod quality; +pub mod relax; pub mod seeds; pub mod size; pub mod surface; diff --git a/crates/mesh/src/relax.rs b/crates/mesh/src/relax.rs new file mode 100644 index 0000000..6f725fd --- /dev/null +++ b/crates/mesh/src/relax.rs @@ -0,0 +1,911 @@ +//! Node relaxation: moving a few nodes of a finished mesh to fix its worst elements, the nodes on +//! the geometry staying on it (`mesh.node.relax`). +//! +//! Each node to move carries a weight in `(0, 1]` (the strength, times a falloff that blends the +//! edited zone into the mesh around it). A sweep visits the nodes in label order and tries, for +//! each, two moves: +//! +//! 1. **Smart Laplacian.** The weighted step toward the centroid of the corner nodes it shares an +//! element edge with, kept only when the worst element around the node does not get worse and +//! no valid element becomes invalid; otherwise the step is halved (three times at most). +//! L. Freitag, *On combining Laplacian and optimization-based mesh smoothing techniques*, AMD +//! 220, 1997 (the quality-guarded Laplacian); S. Canann, J. Tristano, M. Staten, *An approach +//! to combined Laplacian and optimization-based smoothing for triangular, quadrilateral, and +//! quad-dominant meshes*, IMR 1998 (smart Laplacian: a move is accepted only if the local +//! quality does not drop). +//! 2. **Optimization.** A step up the gradient of the worst element's quality around the node +//! (central differences), with a backtracking line search, kept only when the worst quality +//! rises: the max–min formulation of L. Freitag, C. Ollivier-Gooch, *Tetrahedral mesh +//! improvement using swapping and smoothing*, IJNME 40, 1997, and P. Knupp, *Achieving finite +//! element mesh quality via optimization of the Jacobian matrix norm and associated +//! quantities*, IJNME 48, 2000, restricted to the active (worst) element. +//! +//! An element's quality is its scaled Jacobian (a triangle's: its smallest angle over 60°), from +//! [`crate::quality::element_quality`]; a solid is valid when its certified Jacobian check is, +//! a shell when its normal keeps its side. So no move inverts an element, and none makes the +//! worst element around a node worse. +//! +//! **Constraints** ([`Freedom`]): an interior node moves freely; a node on a geometry face moves +//! in the face's parameters and is evaluated back on it ([`Constraints::slide_face`]); a node on a +//! geometry edge slides along it ([`Constraints::slide_curve`]); a fixed node (a geometry vertex, +//! or a skin node whose geometry is unknown) does not move. Mid-side nodes of quadratic elements +//! follow their edge's corners, keeping their offset from the edge's midpoint (the curvature), +//! and are put back on their face or edge. + +use std::collections::{BTreeMap, BTreeSet}; + +use fieldcraft_model::{ElementKind, Model, RelaxationCell}; + +use crate::quality::{self, Metric, Validity}; +use crate::{Error, Result}; + +type P = [f64; 3]; + +/// Most sweeps over the nodes. +pub const MAX_ITERATIONS: usize = 100; +/// Most nodes one relaxation moves. +pub const MAX_NODES: usize = 200_000; +/// Elements listed as the worst of the region. +pub const WORST_LISTED: usize = 5; + +/// How a node may move. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum Freedom { + /// Anywhere (an interior node). + Free, + /// Not at all. + Fixed, + /// On the face with this id (the caller's). + Face(u32), + /// Along the curve with this id (the caller's). + Curve(u32), +} + +/// Where nodes may go: the caller's faces and curves. +pub trait Constraints { + /// How `node` may move; `None` when the caller does not know the node's geometry: it then + /// moves freely inside the mesh and is held on its skin. + fn freedom(&self, node: u32) -> Option; + /// `p` (on face `face`) moved by `d` in the face's parameters; `None` refuses the move. + fn slide_face(&self, face: u32, p: P, d: P) -> Option

; + /// `p` (on curve `curve`) moved by `d` along the curve; `None` refuses the move. + fn slide_curve(&self, curve: u32, p: P, d: P) -> Option

; +} + +/// No geometry: interior nodes move, every node on the skin is held. +#[derive(Debug, Clone, Copy, Default)] +pub struct NoGeometry; + +impl Constraints for NoGeometry { + fn freedom(&self, _: u32) -> Option { + None + } + fn slide_face(&self, _: u32, _: P, _: P) -> Option

{ + None + } + fn slide_curve(&self, _: u32, _: P, _: P) -> Option

{ + None + } +} + +/// How the weight falls from the centre of a sphere of influence to its radius. +#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)] +pub enum Falloff { + /// `(1 − t²)²` of the relative distance `t`: 1 at the centre, 0 with a zero slope at the + /// radius, so the edited zone blends into the mesh around it. + #[default] + Smooth, + /// `1 − t`. + Linear, + /// The full strength everywhere inside. + None, +} + +impl Falloff { + pub fn from_name(name: &str) -> Option { + match name { + "smooth" => Some(Falloff::Smooth), + "linear" => Some(Falloff::Linear), + "none" => Some(Falloff::None), + _ => None, + } + } + + /// The weight at relative distance `t` in `[0, 1]`. + pub fn weight(self, t: f64) -> f64 { + let t = t.clamp(0.0, 1.0); + match self { + Falloff::Smooth => (1.0 - t * t).powi(2), + Falloff::Linear => 1.0 - t, + Falloff::None => 1.0, + } + } +} + +/// The quality of the elements around the moved nodes. +#[derive(Debug, Clone, PartialEq, Default)] +pub struct RegionQuality { + pub elements: usize, + /// Smallest element quality (scaled Jacobian; a triangle's smallest angle over 60°). + pub min_quality: f64, + /// Largest ratio of an element's longest to its shortest edge. + pub max_edge_ratio: f64, + /// Elements whose Jacobian (or a shell's normal) is not valid. + pub invalid: usize, + /// The worst elements and their quality, worst first. + pub worst: Vec<(u32, f64)>, +} + +/// What [`relax`] did. +#[derive(Debug, Clone, PartialEq, Default)] +pub struct Relaxed { + /// The new position of every node that moved (corners and the mid-side nodes following them). + pub positions: BTreeMap, + /// Corner nodes that moved. + pub moved: usize, + /// Largest distance a node moved. + pub max_displacement: f64, + pub before: RegionQuality, + pub after: RegionQuality, + /// Sweeps made (fewer than asked when a sweep moved nothing). + pub iterations: usize, +} + +/// The corner nodes of measured elements (solids and shells) inside the sphere, weighted by +/// `strength` × the falloff; only nodes of `only` when given. +pub fn sphere_weights( + model: &Model, + center: P, + radius: f64, + strength: f64, + falloff: Falloff, + only: Option<&BTreeSet>, +) -> Result> { + if !(radius > 0.0 && radius.is_finite()) || !center.iter().all(|c| c.is_finite()) { + return Err(Error::Invalid("the sphere needs a finite centre and a positive radius".into())); + } + check_strength(strength)?; + let mut out = BTreeMap::new(); + for n in corner_nodes(model) { + if only.is_some_and(|o| !o.contains(&n)) { + continue; + } + let Some(node) = model.node(n) else { continue }; + let t = dist(node.xyz, center) / radius; + if t < 1.0 { + let w = strength * falloff.weight(t); + if w > 0.0 { + out.insert(n, w); + } + } + if out.len() > MAX_NODES { + return Err(Error::Limit(format!("more than {MAX_NODES} nodes inside the sphere; give a smaller radius"))); + } + } + Ok(out) +} + +/// The corner nodes among `nodes` (nodes of measured elements that are a corner of one), each +/// weighted `strength`. +pub fn uniform_weights(model: &Model, nodes: &BTreeSet, strength: f64) -> Result> { + check_strength(strength)?; + let corners = corner_nodes(model); + let out: BTreeMap = nodes.iter().filter(|n| corners.contains(n)).map(|&n| (n, strength)).collect(); + if out.len() > MAX_NODES { + return Err(Error::Limit(format!("more than {MAX_NODES} nodes to relax"))); + } + Ok(out) +} + +fn check_strength(strength: f64) -> Result<()> { + if strength > 0.0 && strength <= 1.0 { Ok(()) } else { Err(Error::Invalid(format!("strength must be in (0, 1], got {strength}"))) } +} + +/// Corner nodes of the solid and shell elements. +fn corner_nodes(model: &Model) -> BTreeSet { + let mut out = BTreeSet::new(); + for e in model.elements().values().filter(|e| quality::is_measured(e.kind)) { + out.extend(e.nodes.iter().take(e.kind.corner_count()).copied()); + } + out +} + +/// Nodes on the skin of the mesh: on a solid face no other solid shares, or on a shell. +fn skin_nodes(model: &Model, near: &BTreeSet) -> BTreeSet { + let mut faces: BTreeMap, (usize, Vec)> = BTreeMap::new(); + let mut skin = BTreeSet::new(); + for e in model.elements().values() { + if !e.nodes.iter().any(|n| near.contains(n)) { + continue; + } + match e.kind.dimension() { + 3 => { + for f in e.kind.faces() { + let nodes: Vec = f.iter().filter_map(|&i| e.nodes.get(i).copied()).collect(); + let mut key: Vec = nodes.iter().take(corner_len(f.len())).copied().collect(); + key.sort_unstable(); + faces.entry(key).or_insert((0, nodes)).0 += 1; + } + } + 2 => skin.extend(e.nodes.iter().copied()), + _ => {} + } + } + for (count, nodes) in faces.into_values() { + if count == 1 { + skin.extend(nodes); + } + } + skin +} + +/// Corners of a face with `n` nodes (3, 4 linear; 6, 8 quadratic). +fn corner_len(n: usize) -> usize { + match n { + 6 => 3, + 8 => 4, + n => n, + } +} + +fn sub(a: P, b: P) -> P { + [a[0] - b[0], a[1] - b[1], a[2] - b[2]] +} +fn add(a: P, b: P) -> P { + [a[0] + b[0], a[1] + b[1], a[2] + b[2]] +} +fn scale(a: P, s: f64) -> P { + [a[0] * s, a[1] * s, a[2] * s] +} +fn dot(a: P, b: P) -> f64 { + a[0] * b[0] + a[1] * b[1] + a[2] * b[2] +} +fn cross(a: P, b: P) -> P { + [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]] +} +fn dist(a: P, b: P) -> f64 { + dot(sub(a, b), sub(a, b)).sqrt() +} + +/// One element of the region. +struct Elem { + label: u32, + kind: ElementKind, + nodes: Vec, + /// A shell's normal as it was: its side, which a move must keep. + normal: Option

, +} + +/// The shell normal from its corners (`None` for solids). +fn shell_normal(kind: ElementKind, x: &[P]) -> Option

{ + if kind.dimension() != 2 { + return None; + } + let c = |i: usize| x.get(i).copied().unwrap_or([0.0; 3]); + Some(match kind.corner_count() { + 3 => cross(sub(c(1), c(0)), sub(c(2), c(0))), + _ => cross(sub(c(2), c(0)), sub(c(3), c(1))), + }) +} + +/// An element's quality and whether it is valid. +fn measure(e: &Elem, x: &[P]) -> (f64, bool) { + let Ok(q) = quality::element_quality(e.kind, x) else { return (-1.0, false) }; + let score = q.get(Metric::ScaledJacobian).or_else(|| q.get(Metric::MinAngle).map(|a| a / 60.0)).unwrap_or(-1.0); + let valid = match (q.validity(), e.normal) { + (Some(v), _) => v == Validity::Valid, + (None, Some(n0)) => shell_normal(e.kind, x).is_some_and(|n| dot(n, n0) > 0.0) && q.size > 0.0, + (None, None) => q.size > 0.0, + }; + (score, valid) +} + +/// Longest over shortest corner edge. +fn edge_ratio(kind: ElementKind, x: &[P]) -> f64 { + let (mut lo, mut hi) = (f64::INFINITY, 0.0f64); + for e in kind.edges() { + if let (Some(&a), Some(&b)) = (e.first().and_then(|&i| x.get(i)), e.get(1).and_then(|&i| x.get(i))) { + let l = dist(a, b); + lo = lo.min(l); + hi = hi.max(l); + } + } + if lo > 0.0 { hi / lo } else { f64::INFINITY } +} + +/// The working state: positions, the region's elements and who depends on whom. +struct State<'a> { + pos: BTreeMap, + elems: Vec, + /// Corner node → its elements (indices into `elems`). + around: BTreeMap>, + /// Corner node → corner nodes it shares an element edge with. + neighbours: BTreeMap>, + /// Mid-side node → its edge's corners and its offset from their midpoint. + mids: BTreeMap, + /// Corner node → the mid-side nodes on its edges. + mids_of: BTreeMap>, + skin: BTreeSet, + geometry: &'a dyn Constraints, +} + +impl State<'_> { + fn at(&self, n: u32) -> P { + self.pos.get(&n).copied().unwrap_or([0.0; 3]) + } + + fn coords(&self, e: &Elem) -> Vec

{ + e.nodes.iter().map(|&n| self.at(n)).collect() + } + + /// How `n` may move: a skin node the caller left free is held. + fn freedom(&self, n: u32) -> Freedom { + match self.geometry.freedom(n) { + Some(f) => f, + None if self.skin.contains(&n) => Freedom::Fixed, + None => Freedom::Free, + } + } + + /// `p` moved by `d` under `freedom`. + fn moved(&self, freedom: Freedom, p: P, d: P) -> Option

{ + let q = match freedom { + Freedom::Fixed => return None, + Freedom::Free => add(p, d), + Freedom::Face(f) => self.geometry.slide_face(f, p, d)?, + Freedom::Curve(c) => self.geometry.slide_curve(c, p, d)?, + }; + q.iter().all(|x| x.is_finite() && x.abs() <= fieldcraft_model::MAX_COORD).then_some(q) + } + + /// The worst quality around `n` and the validity of each element there. + fn local(&self, n: u32) -> (f64, Vec) { + let mut worst = f64::INFINITY; + let mut valid = Vec::new(); + for &i in self.around.get(&n).map(Vec::as_slice).unwrap_or(&[]) { + let Some(e) = self.elems.get(i) else { continue }; + let (q, v) = measure(e, &self.coords(e)); + worst = worst.min(q); + valid.push(v); + } + (worst, valid) + } + + /// The worst element around `n`: its index and quality. + fn worst_around(&self, n: u32) -> Option<(usize, f64)> { + let mut best: Option<(usize, f64)> = None; + for &i in self.around.get(&n)? { + let Some(e) = self.elems.get(i) else { continue }; + let (q, _) = measure(e, &self.coords(e)); + if best.is_none_or(|(_, b)| q < b) { + best = Some((i, q)); + } + } + best + } + + /// Put `n` at `p` with its mid-side nodes following; returns what to restore. + fn place(&mut self, n: u32, p: P) -> Vec<(u32, P)> { + let mut undo = vec![(n, self.at(n))]; + self.pos.insert(n, p); + for m in self.mids_of.get(&n).cloned().unwrap_or_default() { + let Some(&(a, b, bulge)) = self.mids.get(&m) else { continue }; + let old = self.at(m); + let guess = add(scale(add(self.at(a), self.at(b)), 0.5), bulge); + let q = match self.freedom(m) { + Freedom::Free => Some(guess), + Freedom::Fixed => None, + f => self.moved(f, old, sub(guess, old)), + }; + if let Some(q) = q { + undo.push((m, old)); + self.pos.insert(m, q); + } + } + undo + } + + fn restore(&mut self, undo: Vec<(u32, P)>) { + for (n, p) in undo.into_iter().rev() { + self.pos.insert(n, p); + } + } + + /// Try moving `n` by `d` (under its freedom): kept when no valid element around it becomes + /// invalid and its worst quality does not drop (rises by more than `gain`). + fn attempt(&mut self, n: u32, freedom: Freedom, d: P, before: &(f64, Vec), gain: f64) -> bool { + let p = self.at(n); + let Some(q) = self.moved(freedom, p, d) else { return false }; + if dist(p, q) == 0.0 { + return false; + } + let undo = self.place(n, q); + let (worst, valid) = self.local(n); + let ok = worst >= before.0 + gain && before.1.iter().zip(&valid).all(|(was, is)| !was || *is); + if !ok { + self.restore(undo); + } + ok + } + + /// The mean length of the element edges at `n`. + fn local_length(&self, n: u32) -> f64 { + let p = self.at(n); + let nb = self.neighbours.get(&n); + let (sum, count) = nb.into_iter().flatten().fold((0.0, 0usize), |(s, c), &m| (s + dist(p, self.at(m)), c + 1)); + if count > 0 { sum / count as f64 } else { 0.0 } + } + + /// One relaxation of `n` with weight `w`; whether it moved. + fn relax_node(&mut self, n: u32, w: f64) -> bool { + let freedom = self.freedom(n); + if freedom == Freedom::Fixed { + return false; + } + let p = self.at(n); + let mut moved = false; + // 1. Smart Laplacian. + let nb: Vec = self.neighbours.get(&n).map(|s| s.iter().copied().collect()).unwrap_or_default(); + if !nb.is_empty() { + let centroid = scale(nb.iter().fold([0.0; 3], |s, &m| add(s, self.at(m))), 1.0 / nb.len() as f64); + let before = self.local(n); + let mut d = scale(sub(centroid, p), w); + for _ in 0..4 { + if self.attempt(n, freedom, d, &before, 0.0) { + moved = true; + break; + } + d = scale(d, 0.5); + } + } + // 2. Optimization of the worst element around the node. + let h = self.local_length(n); + if !(h > 0.0) { + return moved; + } + let Some((worst, q0)) = self.worst_around(n) else { return moved }; + let Some(e) = self.elems.get(worst) else { return moved }; + let slot = e.nodes.iter().position(|&m| m == n); + let Some(slot) = slot else { return moved }; + let mut x = self.coords(e); + let here = self.at(n); + let eps = 1e-4 * h; + let mut grad = [0.0; 3]; + for k in 0..3 { + let mut f = [0.0; 2]; + for (j, sign) in [1.0, -1.0].into_iter().enumerate() { + if let Some(c) = x.get_mut(slot) { + c[k] = here[k] + sign * eps; + } + f[j] = measure(e, &x).0; + } + if let Some(c) = x.get_mut(slot) { + c[k] = here[k]; + } + grad[k] = (f[0] - f[1]) / (2.0 * eps); + } + let g = dot(grad, grad).sqrt(); + if !(g > 0.0 && g.is_finite()) || q0 >= 1.0 - 1e-9 { + return moved; + } + let before = self.local(n); + let mut step = 0.25 * w * h; + for _ in 0..5 { + if self.attempt(n, freedom, scale(grad, step / g), &before, 1e-9) { + return true; + } + step *= 0.5; + } + moved + } +} + +/// Quality over the region. +fn region_quality(s: &State<'_>) -> RegionQuality { + let mut r = RegionQuality { elements: s.elems.len(), min_quality: f64::INFINITY, max_edge_ratio: 0.0, invalid: 0, worst: Vec::new() }; + let mut all = Vec::with_capacity(s.elems.len()); + for e in &s.elems { + let x = s.coords(e); + let (q, valid) = measure(e, &x); + r.min_quality = r.min_quality.min(q); + r.max_edge_ratio = r.max_edge_ratio.max(edge_ratio(e.kind, &x)); + r.invalid += usize::from(!valid); + all.push((e.label, q)); + } + all.sort_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0))); + all.truncate(WORST_LISTED); + r.worst = all; + if r.elements == 0 { + r.min_quality = 0.0; + } + r +} + +/// Relax the nodes of `weights` (node → weight in `(0, 1]`) for `iterations` sweeps, nodes +/// moving as `geometry` allows. The model is not changed: the result has the new positions. +pub fn relax(model: &Model, weights: &BTreeMap, iterations: usize, geometry: &dyn Constraints) -> Result { + if iterations == 0 || iterations > MAX_ITERATIONS { + return Err(Error::Invalid(format!("iterations must be 1 to {MAX_ITERATIONS}, got {iterations}"))); + } + if weights.len() > MAX_NODES { + return Err(Error::Limit(format!("more than {MAX_NODES} nodes to relax"))); + } + if let Some((n, w)) = weights.iter().find(|(_, w)| !(**w > 0.0 && **w <= 1.0)) { + return Err(Error::Invalid(format!("node {n} has weight {w}; weights are in (0, 1]"))); + } + let movers: BTreeSet = weights.keys().copied().collect(); + let mut state = State { + pos: BTreeMap::new(), + elems: Vec::new(), + around: BTreeMap::new(), + neighbours: BTreeMap::new(), + mids: BTreeMap::new(), + mids_of: BTreeMap::new(), + skin: skin_nodes(model, &movers), + geometry, + }; + for (&label, e) in model.elements() { + if !quality::is_measured(e.kind) || e.nodes.len() != e.kind.node_count() { + continue; + } + let corners = e.kind.corner_count(); + if !e.nodes.iter().take(corners).any(|n| movers.contains(n)) { + continue; + } + let mut x = Vec::with_capacity(e.nodes.len()); + for &n in &e.nodes { + let p = model.node(n).map(|n| n.xyz).ok_or_else(|| Error::Invalid(format!("element {label} refers to missing node {n}")))?; + state.pos.insert(n, p); + x.push(p); + } + let index = state.elems.len(); + for &n in e.nodes.iter().take(corners) { + if movers.contains(&n) { + state.around.entry(n).or_default().push(index); + } + } + for edge in e.kind.edges() { + let (Some(&a), Some(&b)) = (edge.first().and_then(|&i| e.nodes.get(i)), edge.get(1).and_then(|&i| e.nodes.get(i))) else { continue }; + for (u, v) in [(a, b), (b, a)] { + if movers.contains(&u) { + state.neighbours.entry(u).or_default().insert(v); + } + } + if let Some(&m) = edge.get(2).and_then(|&i| e.nodes.get(i)) { + let (pa, pb, pm) = (x.get(edge[0]).copied(), x.get(edge[1]).copied(), x.get(edge[2]).copied()); + if let (Some(pa), Some(pb), Some(pm)) = (pa, pb, pm) { + state.mids.entry(m).or_insert((a, b, sub(pm, scale(add(pa, pb), 0.5)))); + for c in [a, b] { + let list = state.mids_of.entry(c).or_default(); + if !list.contains(&m) { + list.push(m); + } + } + } + } + } + state.elems.push(Elem { label, kind: e.kind, normal: shell_normal(e.kind, &x), nodes: e.nodes.clone() }); + } + let original = state.pos.clone(); + let before = region_quality(&state); + let mut done = 0; + for _ in 0..iterations { + done += 1; + let mut any = false; + for (&n, &w) in weights { + if state.around.contains_key(&n) { + any |= state.relax_node(n, w); + } + } + if !any { + break; + } + } + let after = region_quality(&state); + let mut out = Relaxed { before, after, iterations: done, ..Default::default() }; + for (n, p) in &state.pos { + let Some(o) = original.get(n) else { continue }; + let d = dist(*o, *p); + if d > 0.0 { + out.max_displacement = out.max_displacement.max(d); + out.moved += usize::from(movers.contains(n)); + out.positions.insert(*n, *p); + } + } + Ok(out) +} + +// ---- Relaxation fields ---- + +/// The linear kind of a quadratic element's corners. +fn linear(kind: ElementKind) -> Option { + Some(match kind { + ElementKind::Tri3 | ElementKind::Tri6 => ElementKind::Tri3, + ElementKind::Quad4 | ElementKind::Quad8 => ElementKind::Quad4, + ElementKind::Tet4 | ElementKind::Tet10 => ElementKind::Tet4, + ElementKind::Wedge6 => ElementKind::Wedge6, + ElementKind::Hex8 | ElementKind::Hex20 => ElementKind::Hex8, + _ => return None, + }) +} + +/// The relaxation as a field ([`fieldcraft_model::Relaxation`]'s cells): every element with a +/// moved node, its corners as they were in `model` and as `relaxed` left them. +pub fn field_cells(model: &Model, relaxed: &Relaxed) -> Vec { + let mut cells = Vec::new(); + for e in model.elements().values() { + let Some(kind) = linear(e.kind) else { continue }; + if !e.nodes.iter().any(|n| relaxed.positions.contains_key(n)) { + continue; + } + let corners: Vec = e.nodes.iter().take(kind.node_count()).copied().collect(); + let at: Vec

= corners.iter().filter_map(|n| model.node(*n).map(|x| x.xyz)).collect(); + if at.len() != corners.len() { + continue; + } + let moved = corners.iter().zip(&at).map(|(n, a)| relaxed.positions.get(n).copied().unwrap_or(*a)).collect(); + cells.push(RelaxationCell { kind, at, moved }); + } + cells +} + +/// Linear shape functions of a cell at parametric `xi`, and the domain's dimension. +fn shape(kind: ElementKind, xi: P) -> (Vec, usize) { + let [u, v, w] = xi; + match kind { + ElementKind::Tri3 => (vec![1.0 - u - v, u, v], 2), + ElementKind::Quad4 => (vec![(1.0 - u) * (1.0 - v), u * (1.0 - v), u * v, (1.0 - u) * v], 2), + ElementKind::Tet4 => (vec![1.0 - u - v - w, u, v, w], 3), + ElementKind::Wedge6 => { + let t = [1.0 - u - v, u, v]; + (t.iter().map(|a| a * (1.0 - w)).chain(t.iter().map(|a| a * w)).collect(), 3) + } + _ => { + let b = [(1.0 - u) * (1.0 - v), u * (1.0 - v), u * v, (1.0 - u) * v]; + (b.iter().map(|a| a * (1.0 - w)).chain(b.iter().map(|a| a * w)).collect(), 3) + } + } +} + +/// How far `xi` lies outside the cell's parametric domain (0 inside). +fn outside(kind: ElementKind, xi: P) -> f64 { + let [u, v, w] = xi; + let box01 = |x: f64| (-x).max(x - 1.0).max(0.0); + match kind { + ElementKind::Tri3 => (-u).max(-v).max(u + v - 1.0).max(0.0), + ElementKind::Tet4 => (-u).max(-v).max(-w).max(u + v + w - 1.0).max(0.0), + ElementKind::Quad4 => box01(u).max(box01(v)), + ElementKind::Wedge6 => (-u).max(-v).max(u + v - 1.0).max(0.0).max(box01(w)), + _ => box01(u).max(box01(v)).max(box01(w)), + } +} + +/// The nearest point of the parametric domain to `xi`. +fn clamp_domain(kind: ElementKind, xi: P) -> P { + let c = |x: f64| x.clamp(0.0, 1.0); + let simplex = |u: f64, v: f64, w: f64| { + let (u, v, w) = (u.max(0.0), v.max(0.0), w.max(0.0)); + let s = u + v + w; + if s > 1.0 { (u / s, v / s, w / s) } else { (u, v, w) } + }; + match kind { + ElementKind::Tri3 => { + let (u, v, _) = simplex(xi[0], xi[1], 0.0); + [u, v, 0.0] + } + ElementKind::Tet4 => { + let (u, v, w) = simplex(xi[0], xi[1], xi[2]); + [u, v, w] + } + ElementKind::Wedge6 => { + let (u, v, _) = simplex(xi[0], xi[1], 0.0); + [u, v, c(xi[2])] + } + ElementKind::Quad4 => [c(xi[0]), c(xi[1]), 0.0], + _ => [c(xi[0]), c(xi[1]), c(xi[2])], + } +} + +fn combine(n: &[f64], x: &[P]) -> P { + n.iter().zip(x).fold([0.0; 3], |s, (w, p)| add(s, scale(*p, *w))) +} + +/// The parametric coordinates of `x` in the cell (Gauss–Newton on the linear map; for a shell, +/// of the nearest point of its surface), and the distance left. +fn locate(c: &RelaxationCell, x: P) -> (P, f64) { + let (_, dim) = shape(c.kind, [0.0; 3]); + let mut xi = match c.kind { + ElementKind::Tri3 | ElementKind::Tet4 => [0.25; 3], + ElementKind::Wedge6 => [1.0 / 3.0, 1.0 / 3.0, 0.5], + _ => [0.5; 3], + }; + if dim == 2 { + xi[2] = 0.0; + } + let map = |xi: P| combine(&shape(c.kind, xi).0, &c.at); + for _ in 0..20 { + let r = sub(x, map(xi)); + // The Jacobian's columns (the map is at most trilinear: a unit-step difference of the + // two sides of `xi` is exact along each axis). + let h = 1e-3; + let mut j = [[0.0; 3]; 3]; + for (k, col) in j.iter_mut().enumerate().take(dim) { + let (mut a, mut b) = (xi, xi); + a[k] += h; + b[k] -= h; + *col = scale(sub(map(a), map(b)), 0.5 / h); + } + // Normal equations JᵀJ δ = Jᵀ r. + let mut m = [[0.0; 3]; 3]; + let mut rhs = [0.0; 3]; + for a in 0..dim { + rhs[a] = dot(j[a], r); + for b in 0..dim { + m[a][b] = dot(j[a], j[b]); + } + } + if dim == 2 { + m[2][2] = 1.0; + } + let det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) + + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]); + if !(det.abs() > 1e-300) || !det.is_finite() { + break; + } + let solve = |col: usize| { + let mut t = m; + for row in 0..3 { + t[row][col] = rhs[row]; + } + (t[0][0] * (t[1][1] * t[2][2] - t[1][2] * t[2][1]) - t[0][1] * (t[1][0] * t[2][2] - t[1][2] * t[2][0]) + + t[0][2] * (t[1][0] * t[2][1] - t[1][1] * t[2][0])) + / det + }; + let delta = [solve(0), solve(1), if dim == 3 { solve(2) } else { 0.0 }]; + xi = add(xi, delta); + if !xi.iter().all(|v| v.is_finite()) || dot(delta, delta) < 1e-28 { + break; + } + } + (xi, dist(x, map(clamp_domain(c.kind, xi)))) +} + +/// The cell's size: its longest corner-to-corner span. +fn cell_size(c: &RelaxationCell) -> f64 { + let mut s = 0.0f64; + for a in &c.at { + for b in &c.at { + s = s.max(dist(*a, *b)); + } + } + s +} + +/// Parametric slack accepted around a cell (a node of a finer mesh on a curved face lies off the +/// coarse cell's flat side by its chord error). +const SLACK: f64 = 0.25; + +/// The field's displacement at `x`: from the cell holding it (or the nearest within the slack), +/// with that cell's shape functions; `None` outside the field. +pub fn field_displacement(cells: &[RelaxationCell], x: P) -> Option

{ + let mut best: Option<(f64, &RelaxationCell, P)> = None; + for c in cells { + let size = cell_size(c); + let lo = c.at.iter().fold([f64::INFINITY; 3], |m, p| [m[0].min(p[0]), m[1].min(p[1]), m[2].min(p[2])]); + let hi = c.at.iter().fold([f64::NEG_INFINITY; 3], |m, p| [m[0].max(p[0]), m[1].max(p[1]), m[2].max(p[2])]); + let pad = SLACK * size; + if (0..3).any(|k| x[k] < lo[k] - pad || x[k] > hi[k] + pad) { + continue; + } + let (xi, gap) = locate(c, x); + let badness = outside(c.kind, xi).max(gap / size.max(f64::MIN_POSITIVE)); + if badness <= SLACK && best.is_none_or(|(b, _, _)| badness < b) { + best = Some((badness, c, xi)); + } + } + let (_, c, xi) = best?; + let n = shape(c.kind, clamp_domain(c.kind, xi)).0; + let d: Vec

= c.moved.iter().zip(&c.at).map(|(m, a)| sub(*m, *a)).collect(); + Some(combine(&n, &d)) +} + +/// Apply a relaxation field to the mesh: each node of `nodes` (default every node of a measured +/// element) inside the field moves by the interpolated displacement as `geometry` allows (on its +/// face or edge); a move that would make a valid element invalid is undone; then one guarded +/// smoothing sweep (strength `polish`, 0 for none) evens out what the interpolation left. The +/// model is not changed. +pub fn apply_field( + model: &Model, + cells: &[RelaxationCell], + nodes: Option<&BTreeSet>, + geometry: &dyn Constraints, + polish: f64, +) -> Result { + let candidates: BTreeSet = match nodes { + Some(n) => n.clone(), + None => model.elements().values().filter(|e| quality::is_measured(e.kind)).flat_map(|e| e.nodes.iter().copied()).collect(), + }; + if candidates.len() > fieldcraft_model::MAX_NODES { + return Err(Error::Limit("too many nodes".into())); + } + let skin = skin_nodes(model, &candidates); + let mut work = model.clone(); + let mut moved: BTreeMap = BTreeMap::new(); + for &n in &candidates { + let Some(x) = model.node(n).map(|p| p.xyz) else { continue }; + let Some(d) = field_displacement(cells, x) else { continue }; + if !(dot(d, d) > 0.0) { + continue; + } + let freedom = match geometry.freedom(n) { + Some(f) => f, + None if skin.contains(&n) => Freedom::Fixed, + None => Freedom::Free, + }; + let q = match freedom { + Freedom::Fixed => None, + Freedom::Free => Some(add(x, d)), + Freedom::Face(f) => geometry.slide_face(f, x, d), + Freedom::Curve(c) => geometry.slide_curve(c, x, d), + }; + if let Some(q) = q.filter(|q| q.iter().all(|v| v.is_finite() && v.abs() <= fieldcraft_model::MAX_COORD)) { + work.set_node(n, q)?; + moved.insert(n, q); + } + } + // Undo the moves around any element they would invalidate. + let touched: Vec<(u32, ElementKind, Vec)> = model + .elements() + .iter() + .filter(|(_, e)| quality::is_measured(e.kind) && e.nodes.iter().any(|n| moved.contains_key(n))) + .map(|(l, e)| (*l, e.kind, e.nodes.clone())) + .collect(); + let at = |m: &Model, nodes: &[u32]| -> Vec

{ nodes.iter().map(|n| m.node(*n).map_or([0.0; 3], |p| p.xyz)).collect() }; + for _ in 0..10 { + let mut reverted = false; + for (label, kind, nodes) in &touched { + let x0 = at(model, nodes); + let e = Elem { label: *label, kind: *kind, nodes: nodes.clone(), normal: shell_normal(*kind, &x0) }; + if measure(&e, &x0).1 && !measure(&e, &at(&work, nodes)).1 { + for n in nodes { + if moved.remove(n).is_some() + && let Some(p) = model.node(*n) + { + work.set_node(*n, p.xyz)?; + reverted = true; + } + } + } + } + if !reverted { + break; + } + } + let mut out = Relaxed::default(); + let corners = corner_nodes(model); + if polish > 0.0 { + let weights = uniform_weights(&work, &moved.keys().copied().filter(|n| corners.contains(n)).collect(), polish.min(1.0))?; + if !weights.is_empty() { + let r = relax(&work, &weights, 1, geometry)?; + for (n, p) in &r.positions { + work.set_node(*n, *p)?; + moved.insert(*n, *p); + } + out.after = r.after; + } + } + for (n, p) in moved { + let Some(o) = model.node(n) else { continue }; + let d = dist(o.xyz, p); + if d > 0.0 { + out.max_displacement = out.max_displacement.max(d); + out.moved += usize::from(corners.contains(&n)); + out.positions.insert(n, p); + } + } + out.iterations = 1; + Ok(out) +} + +#[cfg(test)] +mod tests; diff --git a/crates/mesh/src/relax/tests.rs b/crates/mesh/src/relax/tests.rs new file mode 100644 index 0000000..26accc5 --- /dev/null +++ b/crates/mesh/src/relax/tests.rs @@ -0,0 +1,344 @@ +use super::*; +use crate::tet::geometry::AnalyticSurface; + +const R_IN: f64 = 1.0; +const R_OUT: f64 = 2.0; +const HEIGHT: f64 = 1.0; + +/// A quarter annulus (r 1..2, θ 0..90°, z 0..1) on its analytic faces: 0 the outer cylinder, 1 +/// the inner, 2 and 3 the planes z = 0 and z = 1, 4 and 5 the planes θ = 0 and θ = 90°. +struct Annulus; + +impl Annulus { + fn surface(f: u32) -> AnalyticSurface { + match f { + 0 => AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: R_OUT }, + 1 => AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: R_IN }, + 2 => AnalyticSurface::Plane { normal: [0.0, 0.0, -1.0], point: [0.0; 3] }, + 3 => AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, HEIGHT] }, + 4 => AnalyticSurface::Plane { normal: [0.0, -1.0, 0.0], point: [0.0; 3] }, + _ => AnalyticSurface::Plane { normal: [-1.0, 0.0, 0.0], point: [0.0; 3] }, + } + } + + /// The faces a point lies on. + fn faces(p: P) -> Vec { + let r = (p[0] * p[0] + p[1] * p[1]).sqrt(); + let tol = 1e-9; + let mut f = Vec::new(); + for (i, on) in + [(r - R_OUT).abs() < tol, (r - R_IN).abs() < tol, p[2].abs() < tol, (p[2] - HEIGHT).abs() < tol, p[1].abs() < tol, p[0].abs() < tol] + .into_iter() + .enumerate() + { + if on { + f.push(i as u32); + } + } + f + } + + fn on(p: P, faces: &[u32]) -> bool { + faces.iter().all(|&f| Annulus::surface(f).project(p).is_some_and(|q| dist(p, q) < 1e-9)) + } +} + +/// Curves are encoded as `a * 8 + b`. +struct Geo { + model: Model, +} + +impl Constraints for Geo { + fn freedom(&self, node: u32) -> Option { + let p = self.model.node(node)?.xyz; + let f = Annulus::faces(p); + Some(match f.as_slice() { + [] => Freedom::Free, + [a] => Freedom::Face(*a), + [a, b] => Freedom::Curve(a * 8 + b), + _ => Freedom::Fixed, + }) + } + fn slide_face(&self, face: u32, p: P, d: P) -> Option

{ + Annulus::surface(face).slide(p, d) + } + fn slide_curve(&self, curve: u32, p: P, d: P) -> Option

{ + let (a, b) = (Annulus::surface(curve / 8), Annulus::surface(curve % 8)); + let mut q = add(p, d); + for _ in 0..100 { + q = b.project(a.project(q)?)?; + } + Some(q) + } +} + +fn polar(r: f64, t: f64, z: f64) -> P { + let (s, c) = crate::tet::geometry::sin_cos(t); + [r * c, r * s, z] +} + +/// A hex8 mesh of the annulus with `n` divisions each way; node label `1 + i + (n+1)(j + (n+1)k)`. +fn hex_annulus(n: usize) -> Model { + let mut m = Model::default(); + let id = |i: usize, j: usize, k: usize| (1 + i + (n + 1) * (j + (n + 1) * k)) as u32; + for k in 0..=n { + for j in 0..=n { + for i in 0..=n { + let p = polar( + R_IN + (R_OUT - R_IN) * i as f64 / n as f64, + std::f64::consts::FRAC_PI_2 * j as f64 / n as f64, + HEIGHT * k as f64 / n as f64, + ); + m.set_node(id(i, j, k), p).unwrap(); + } + } + } + let mut label = 1; + for k in 0..n { + for j in 0..n { + for i in 0..n { + let nodes = vec![ + id(i, j, k), + id(i + 1, j, k), + id(i + 1, j + 1, k), + id(i, j + 1, k), + id(i, j, k + 1), + id(i + 1, j, k + 1), + id(i + 1, j + 1, k + 1), + id(i, j + 1, k + 1), + ]; + m.set_element(label, ElementKind::Hex8, nodes).unwrap(); + label += 1; + } + } + } + m +} + +/// Move some nodes near `centre`, keeping each on its faces. +fn perturb(m: &mut Model, centre: P, radius: f64, amount: f64) -> BTreeSet { + let geo = Geo { model: m.clone() }; + let mut moved = BTreeSet::new(); + let labels: Vec = m.nodes().keys().copied().collect(); + for (k, n) in labels.into_iter().enumerate() { + let p = m.node(n).unwrap().xyz; + if dist(p, centre) > radius { + continue; + } + let d = [amount * ((k * 7 % 5) as f64 - 2.0) / 2.0, amount * ((k * 3 % 5) as f64 - 2.0) / 2.0, amount * ((k * 11 % 5) as f64 - 2.0) / 2.0]; + let q = match geo.freedom(n).unwrap() { + Freedom::Free => Some(add(p, d)), + Freedom::Face(f) => geo.slide_face(f, p, d), + _ => None, + }; + if let Some(q) = q { + m.set_node(n, q).unwrap(); + moved.insert(n); + } + } + moved +} + +fn region_min(m: &Model, labels: &[u32]) -> (f64, usize) { + let mut worst = f64::INFINITY; + let mut invalid = 0; + for l in labels { + let e = m.element(*l).unwrap(); + let x: Vec

= e.nodes.iter().map(|n| m.node(*n).unwrap().xyz).collect(); + let el = Elem { label: *l, kind: e.kind, nodes: e.nodes.clone(), normal: None }; + let (q, v) = measure(&el, &x); + worst = worst.min(q); + invalid += usize::from(!v); + } + (worst, invalid) +} + +fn applied(m: &Model, r: &Relaxed) -> Model { + let mut out = m.clone(); + for (n, p) in &r.positions { + out.set_node(*n, *p).unwrap(); + } + out +} + +const CENTRE: P = [1.414, 1.414, 0.5]; + +/// The issue's test: a perturbed hex mesh of a cylindrical part relaxes inside the sphere — every +/// surface node on its face, quality up, nothing inverted, nothing outside the sphere moved. +#[test] +fn hex_spikes_relax_with_surface_nodes_on_the_cylinder() { + let mut m = hex_annulus(6); + let original = m.clone(); + perturb(&mut m, CENTRE, 0.5, 0.08); + let elements: Vec = m.elements().keys().copied().collect(); + let (before, _) = region_min(&m, &elements); + let weights = sphere_weights(&m, CENTRE, 0.7, 1.0, Falloff::Smooth, None).unwrap(); + let geo = Geo { model: original.clone() }; + let r = relax(&m, &weights, 10, &geo).unwrap(); + assert!(r.moved > 0 && r.after.min_quality > r.before.min_quality, "{:?} → {:?}", r.before, r.after); + assert_eq!(r.after.invalid, 0); + let out = applied(&m, &r); + let (after, invalid) = region_min(&out, &elements); + assert!(after > before + 0.02 && invalid == 0, "{before} → {after}"); + for (n, p) in &r.positions { + assert!(dist(m.node(*n).unwrap().xyz, CENTRE) < 0.7, "node {n} outside the sphere moved"); + let faces = Annulus::faces(original.node(*n).unwrap().xyz); + assert!(Annulus::on(*p, &faces), "node {n} left its faces {faces:?}: {p:?}"); + } +} + +/// tet10 with curved mid-side nodes on the cylinder: after relaxing, corner and mid-side surface +/// nodes are still on the surface and no element is invalid. +#[test] +fn tet10_mid_side_nodes_follow_on_the_geometry() { + let hex = hex_annulus(4); + let mut m = Model::default(); + for (l, n) in hex.nodes() { + m.set_node(*l, n.xyz).unwrap(); + } + // Six tetrahedra per hexahedron around its 1-7 diagonal (Kuhn), consistent across faces. + const KUHN: [[usize; 4]; 6] = [[0, 1, 2, 6], [0, 2, 3, 6], [0, 3, 7, 6], [0, 7, 4, 6], [0, 4, 5, 6], [0, 5, 1, 6]]; + let mut mids: BTreeMap<(u32, u32), u32> = BTreeMap::new(); + let mut next = 100_000; + let mut label = 1; + for e in hex.elements().values() { + for t in KUHN { + let mut c: Vec = t.iter().map(|&i| e.nodes[i]).collect(); + let x: Vec

= c.iter().map(|n| m.node(*n).unwrap().xyz).collect(); + if dot(cross(sub(x[1], x[0]), sub(x[2], x[0])), sub(x[3], x[0])) < 0.0 { + c.swap(1, 2); + } + let mut nodes = c.clone(); + for [a, b] in [[0, 1], [1, 2], [2, 0], [0, 3], [1, 3], [2, 3]] { + let key = (c[a].min(c[b]), c[a].max(c[b])); + let mid = *mids.entry(key).or_insert_with(|| { + let (pa, pb) = (m.node(key.0).unwrap().xyz, m.node(key.1).unwrap().xyz); + let mut p = scale(add(pa, pb), 0.5); + let shared: Vec = Annulus::faces(pa).into_iter().filter(|f| Annulus::faces(pb).contains(f)).collect(); + if let [f] = shared.as_slice() { + p = Annulus::surface(*f).project(p).unwrap(); + } + next += 1; + m.set_node(next, p).unwrap(); + next + }); + nodes.push(mid); + } + m.set_element(label, ElementKind::Tet10, nodes).unwrap(); + label += 1; + } + } + let original = m.clone(); + // Perturb the corners only, then put the mid-side nodes at their edges' (curved) middles. + let corners = corner_nodes(&m); + perturb(&mut m, CENTRE, 0.6, 0.06); + for (&(a, b), &mid) in &mids { + let p = scale(add(m.node(a).unwrap().xyz, m.node(b).unwrap().xyz), 0.5); + let shared: Vec = + Annulus::faces(original.node(a).unwrap().xyz).into_iter().filter(|f| Annulus::faces(original.node(b).unwrap().xyz).contains(f)).collect(); + let p = if let [f] = shared.as_slice() { Annulus::surface(*f).project(p).unwrap() } else { p }; + m.set_node(mid, p).unwrap(); + } + let _ = corners; + let weights = sphere_weights(&m, CENTRE, 0.8, 1.0, Falloff::Smooth, None).unwrap(); + let r = relax(&m, &weights, 8, &Geo { model: original.clone() }).unwrap(); + assert!(r.after.min_quality >= r.before.min_quality && r.after.invalid == 0, "{:?} → {:?}", r.before, r.after); + assert!(r.after.min_quality > r.before.min_quality + 0.01, "{:?} → {:?}", r.before, r.after); + let mut mid_on_surface = 0; + for (n, p) in &r.positions { + let faces = Annulus::faces(original.node(*n).unwrap().xyz); + assert!(Annulus::on(*p, &faces), "node {n} left {faces:?}"); + mid_on_surface += usize::from(*n > 100_000 && faces.contains(&0)); + } + assert!(mid_on_surface > 0, "some curved mid-side nodes moved on the cylinder"); +} + +/// A cut between two cells is not geometry: its shared nodes move freely, and both cells keep +/// them (conforming). +#[test] +fn partition_cut_nodes_move_freely_and_stay_shared() { + let mut m = hex_annulus(6); + let original = m.clone(); + perturb(&mut m, CENTRE, 0.5, 0.08); + // The cut: the plane θ = 45° through the middle row of nodes (j = 3) — reported free. + struct Cut(Geo); + impl Constraints for Cut { + fn freedom(&self, node: u32) -> Option { + let f = self.0.freedom(node)?; + Some(if f == Freedom::Free { Freedom::Free } else { f }) + } + fn slide_face(&self, f: u32, p: P, d: P) -> Option

{ + self.0.slide_face(f, p, d) + } + fn slide_curve(&self, c: u32, p: P, d: P) -> Option

{ + self.0.slide_curve(c, p, d) + } + } + let weights = sphere_weights(&m, CENTRE, 0.7, 1.0, Falloff::Smooth, None).unwrap(); + let r = relax(&m, &weights, 10, &Cut(Geo { model: original.clone() })).unwrap(); + let on_cut = |p: P| (p[0] - p[1]).abs() < 1e-9; + let cut_moved = r.positions.iter().filter(|(n, p)| on_cut(original.node(**n).unwrap().xyz) && !on_cut(**p)).count(); + assert!(cut_moved > 0, "cut nodes are free to leave the cut plane"); + assert_eq!(r.after.invalid, 0); +} + +/// The issue's field test: relax a coarse mesh, mesh 2× finer, apply the field: the fine mesh +/// takes the relaxed shape (its nodes where the relaxed coarse mesh interpolates them) and +/// nothing inverts. +#[test] +fn a_relaxation_field_carries_over_to_a_finer_mesh() { + let mut coarse = hex_annulus(4); + let original = coarse.clone(); + // A deliberate interior shape change: pull the interior nodes toward the centre. + for (n, p) in original.nodes() { + if Annulus::faces(p.xyz).is_empty() && dist(p.xyz, CENTRE) < 0.6 { + coarse.set_node(*n, add(p.xyz, [0.03, -0.02, 0.04])).unwrap(); + } + } + let weights = sphere_weights(&coarse, CENTRE, 0.9, 1.0, Falloff::Smooth, None).unwrap(); + let r = relax(&coarse, &weights, 5, &Geo { model: coarse.clone() }).unwrap(); + // The field is the change from the original coarse mesh to the relaxed one. + let relaxed = applied(&coarse, &r); + let total = Relaxed { + positions: relaxed.nodes().iter().filter(|(n, p)| original.node(**n).unwrap().xyz != p.xyz).map(|(n, p)| (*n, p.xyz)).collect(), + ..Default::default() + }; + let cells = field_cells(&original, &total); + assert!(!cells.is_empty()); + let fine = hex_annulus(8); + let out = apply_field(&fine, &cells, None, &Geo { model: fine.clone() }, 0.0).unwrap(); + assert!(out.moved > 0); + let result = applied(&fine, &out); + let elements: Vec = fine.elements().keys().copied().collect(); + assert_eq!(region_min(&result, &elements).1, 0, "nothing inverts"); + for (n, p) in fine.nodes() { + let d = field_displacement(&cells, p.xyz).unwrap_or([0.0; 3]); + let got = result.node(*n).unwrap().xyz; + if Annulus::faces(p.xyz).is_empty() { + assert!(dist(got, add(p.xyz, d)) < 1e-9, "node {n}: {got:?} vs {:?}", add(p.xyz, d)); + } else { + assert!(Annulus::on(got, &Annulus::faces(p.xyz)), "node {n} stays on its faces"); + } + } + // A fine node at a coarse node takes exactly that node's displacement. + let coarse_node = original.nodes().iter().find(|(n, p)| r.positions.contains_key(n) && Annulus::faces(p.xyz).is_empty()).unwrap(); + let fine_twin = fine.nodes().values().find(|p| dist(p.xyz, coarse_node.1.xyz) < 1e-9).unwrap(); + let d = field_displacement(&cells, fine_twin.xyz).unwrap(); + assert!(dist(add(fine_twin.xyz, d), relaxed.node(*coarse_node.0).unwrap().xyz) < 1e-8); +} + +#[test] +fn hostile_parameters_are_refused() { + let m = hex_annulus(2); + assert!(sphere_weights(&m, [f64::NAN, 0.0, 0.0], 1.0, 1.0, Falloff::Smooth, None).is_err()); + assert!(sphere_weights(&m, CENTRE, -1.0, 1.0, Falloff::Smooth, None).is_err()); + assert!(sphere_weights(&m, CENTRE, 1.0, 1.5, Falloff::Smooth, None).is_err()); + let w = sphere_weights(&m, CENTRE, 1.0, 1.0, Falloff::Smooth, None).unwrap(); + assert!(relax(&m, &w, 0, &NoGeometry).is_err()); + assert!(relax(&m, &w, MAX_ITERATIONS + 1, &NoGeometry).is_err()); + let bad: BTreeMap = [(1, f64::NAN)].into(); + assert!(relax(&m, &bad, 1, &NoGeometry).is_err()); + // Without geometry, skin nodes stay put. + let r = relax(&m, &w, 3, &NoGeometry).unwrap(); + assert!(r.positions.is_empty() || r.positions.keys().all(|n| Annulus::faces(m.node(*n).unwrap().xyz).is_empty())); +} diff --git a/crates/mesh/src/tet/geometry.rs b/crates/mesh/src/tet/geometry.rs index aa4340f..d05bc21 100644 --- a/crates/mesh/src/tet/geometry.rs +++ b/crates/mesh/src/tet/geometry.rs @@ -117,6 +117,76 @@ impl AnalyticSurface { out.iter().all(|x| x.is_finite()).then_some(out) } + /// Move a point of the surface by the displacement `d` in the surface's own parameters (UV): + /// `p` is first put on the surface, `d` is resolved into the parameter steps along the + /// surface's coordinate directions there (the tangential part of `d`; the normal part is + /// dropped), and the surface is evaluated at the new parameters, so the result lies exactly + /// on the surface. Parameters are local to `p` (angles measured from it), so the chart has no + /// seam. A plane translates in its own plane; a cylinder, cone and torus turn about their + /// axis (and a torus round its tube) by the arc length over the radius; a sphere follows the + /// great circle in the direction of `d`. Angular steps are kept within half a turn. `None` + /// where the parameters are singular (on an axis, at an apex) or the data is not finite. + pub fn slide(&self, p: P, d: P) -> Option

{ + let q = self.project(p)?; + let rotate = |r: P, t: P, angle: f64| { + let (s, c) = sin_cos(angle.clamp(-std::f64::consts::PI, std::f64::consts::PI)); + add(scale(r, c), scale(t, s)) + }; + let out = match *self { + AnalyticSurface::Plane { normal, .. } => { + let n = unit(normal)?; + add(q, sub(d, scale(n, dot(d, n)))) + } + AnalyticSurface::Cylinder { origin, axis, radius } => { + let a = unit(axis)?; + let foot = add(origin, scale(a, dot(sub(q, origin), a))); + let r = unit(sub(q, foot))?; + let t = cross(a, r); + let radial = rotate(r, t, dot(d, t) / radius); + add(add(foot, scale(a, dot(d, a))), scale(radial, radius)) + } + AnalyticSurface::Sphere { center, radius } => { + let r = unit(sub(q, center))?; + let tangent = sub(d, scale(r, dot(d, r))); + let Some(t) = unit(tangent) else { return Some(q) }; + add(center, scale(rotate(r, t, norm(tangent) / radius), radius)) + } + AnalyticSurface::Cone { apex, axis, half_angle } => { + let a = unit(axis)?; + let g0 = sub(q, apex); + let along = norm(g0); + let g = unit(g0)?; + let r = unit(sub(g0, scale(a, dot(g0, a))))?; + let t = cross(a, r); + let rho = norm(sub(g0, scale(a, dot(g0, a)))); + let s = along + dot(d, g); + if !(s > 0.0 && rho > 0.0) { + return None; + } + let r2 = rotate(r, t, dot(d, t) / rho); + let (sa, ca) = sin_cos(half_angle); + add(apex, scale(add(scale(a, ca), scale(r2, sa)), s)) + } + AnalyticSurface::Torus { center, axis, major, minor } => { + let a = unit(axis)?; + let c = sub(q, center); + let radial = sub(c, scale(a, dot(c, a))); + let rho = norm(radial); + let r = unit(radial)?; + let ring = add(center, scale(r, major)); + let m = unit(sub(q, ring))?; + let around = cross(a, r); + let tube = cross(around, m); + // Round the tube, then about the axis. + let m2 = rotate(m, tube, dot(d, tube) / minor); + let r2 = rotate(r, around, dot(d, around) / rho.max(f64::MIN_POSITIVE)); + let m3 = add(scale(r2, dot(m2, r)), scale(a, dot(m2, a))); + add(add(center, scale(r2, major)), scale(m3, minor)) + } + }; + out.iter().all(|x| x.is_finite()).then_some(out) + } + /// The unit normal (up to sign) at a point on the surface, if defined there. pub fn normal_at(&self, q: P) -> Option

{ let n = match *self { @@ -473,6 +543,31 @@ impl BoundaryGeometry for Tessellation { } } +impl Tessellation { + /// Move a point of face `face` by the displacement `d` and keep it on the face: through the + /// face's own parameters ([`AnalyticSurface::slide`]) when its exact surface is known and the + /// result stays near the face's triangles (not on another part of the surface), else by + /// projecting `p + d` onto the face's triangles (free-form faces: accurate to the + /// tessellation's chordal tolerance, as [`BoundaryGeometry::project_face`]). + pub fn slide_face(&self, face: u32, p: P, d: P) -> Option

{ + let target = add(p, d); + let (near, scale_len) = self.on_tessellation(face, target)?; + if let Some(q) = self.analytic(face).and_then(|s| s.slide(p, d)) { + let tol = (0.5 * scale_len).max(2.0 * measure::dist(target, near)); + if measure::dist(q, near) <= tol { + return Some(q); + } + } + self.project_face(face, target) + } + + /// Move a point of the feature curve between faces `faces` by `d` and keep it on the curve + /// (the curve's closest point to `p + d`: the move along the curve's parameter). + pub fn slide_curve(&self, faces: [u32; 2], p: P, d: P) -> Option

{ + self.project_curve(faces, add(p, d)) + } +} + impl Tessellation { /// The local thickness at each vertex: the distance to the nearest triangle facing it from /// across the solid (a triangle that does not touch the vertex and whose outward normal @@ -737,4 +832,36 @@ mod tests { assert_eq!(far, 10.0); assert!(f.size([f64::NAN, 0.0, 0.0]) > 0.0); } + + /// A slide in UV keeps the point exactly on every kind of surface, moves it by the tangential + /// part of the displacement (to first order), and turns a cylinder point by arc / radius. + #[test] + fn slides_stay_on_the_surface() { + let surfaces = [ + AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, 2.0] }, + AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: 3.0 }, + AnalyticSurface::Sphere { center: [1.0, 0.0, 0.0], radius: 2.0 }, + AnalyticSurface::Cone { apex: [0.0; 3], axis: [0.0, 0.0, 1.0], half_angle: 0.4 }, + AnalyticSurface::Torus { center: [0.0; 3], axis: [0.0, 0.0, 1.0], major: 5.0, minor: 1.0 }, + ]; + let probe = [2.5, 1.5, 2.2]; + for s in surfaces { + let q = s.project(probe).unwrap(); + for d in [[0.01, 0.0, 0.0], [0.0, 0.02, 0.0], [0.0, 0.0, -0.015], [0.3, -0.2, 0.1]] { + let moved = s.slide(q, d).unwrap(); + let back = s.project(moved).unwrap(); + assert!(measure::dist(moved, back) < 1e-12, "{s:?} {d:?}: {moved:?} off by {}", measure::dist(moved, back)); + let n = s.normal_at(q).unwrap(); + let tangential = sub(d, scale(n, dot(d, n))); + if norm(d) < 0.05 { + assert!(measure::dist(moved, add(q, tangential)) < 0.05 * norm(d) + 1e-12, "{s:?} {d:?}"); + } + } + } + let cyl = AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: 2.0 }; + let q = cyl.slide([2.0, 0.0, 1.0], [0.0, std::f64::consts::PI, 0.5]).unwrap(); + assert!(measure::dist(q, [0.0, 2.0, 1.5]) < 1e-12, "{q:?}"); + // On the axis the parameters are singular. + assert!(cyl.slide([0.0, 0.0, 1.0], [1.0, 0.0, 0.0]).is_none()); + } } diff --git a/crates/model/src/lib.rs b/crates/model/src/lib.rs index 8669e81..3693a0e 100644 --- a/crates/model/src/lib.rs +++ b/crates/model/src/lib.rs @@ -44,6 +44,7 @@ pub mod materials; mod mesh; pub mod nodal; pub mod progress; +pub mod relaxations; pub mod rename; pub mod sections; pub mod seeds; @@ -70,6 +71,7 @@ pub use mesh::{ Element, ElementEdge, ElementFace, ElementKind, MIN_NORMALIZED_JACOBIAN, Node, admissible_jacobian, check_coordinates, normalized_jacobian, tet4_jacobian, }; +pub use relaxations::{Relaxation, RelaxationCell}; pub use sections::{BeamProfile, Extent, Formulation, Integration, Section, SectionKind}; pub use seeds::{Bias, BiasToward, Divisions, Seed, SeedConstraint}; pub use sets::Set; @@ -162,6 +164,9 @@ pub struct Model { /// Edge seeds, sorted by edge. #[serde(skip_serializing_if = "Vec::is_empty")] seeds: Vec, + /// Relaxation fields (`mesh.node.relax`), by name. + #[serde(skip_serializing_if = "BTreeMap::is_empty")] + relaxations: BTreeMap, } /// A [`Model`] as stored, before it is checked. @@ -210,6 +215,8 @@ struct ModelData { suppressed_datums: BTreeSet, #[serde(default)] seeds: Vec, + #[serde(default)] + relaxations: BTreeMap, } impl TryFrom for Model { @@ -293,6 +300,9 @@ impl TryFrom for Model { for s in d.seeds { m.set_seed(s)?; } + for (n, v) in d.relaxations { + m.set_relaxation(&n, v)?; + } Ok(m) } } @@ -423,9 +433,33 @@ impl Model { partitions: self.partitions.clone(), suppressed_datums: self.suppressed_datums.clone(), seeds: self.seeds.clone(), + relaxations: self.relaxations.clone(), } } + /// The relaxation fields, by name. + pub fn relaxations(&self) -> &BTreeMap { + &self.relaxations + } + + pub fn relaxation(&self, name: &str) -> Option<&Relaxation> { + self.relaxations.get(name) + } + + /// Add a relaxation field, or replace the one of that name. + pub fn set_relaxation(&mut self, name: &str, field: Relaxation) -> Result<()> { + field.check().map_err(|e| context("relaxation", name, e))?; + if self.relaxations.contains_key(name) { + replace_existing(&mut self.relaxations, "relaxation", name, field) + } else { + insert_new(&mut self.relaxations, "relaxation", name, field) + } + } + + pub fn remove_relaxation(&mut self, name: &str) -> Result { + remove_existing(&mut self.relaxations, "relaxation", name) + } + /// Counts of each kind of definition, for summaries. pub fn definition_counts(&self) -> [(&'static str, usize); 16] { [ diff --git a/crates/model/src/relaxations.rs b/crates/model/src/relaxations.rs new file mode 100644 index 0000000..35bc5a1 --- /dev/null +++ b/crates/model/src/relaxations.rs @@ -0,0 +1,60 @@ +//! Relaxation fields: what a node relaxation (`mesh.node.relax`) did, kept so it can be done again +//! on a new mesh of the same geometry. +//! +//! A field is the relaxed region of the mesh it was made on: each element there as its corners +//! were (`at`) and where they went (`moved`). A node of a later mesh that lies in one of these +//! cells takes the displacement interpolated with the cell's own shape functions at its +//! parametric coordinates (barycentric in a tetrahedron or triangle, trilinear `(u, v, w)` in a +//! hexahedron, …), so the relaxed shape carries over to a coarser or finer mesh. Fields are +//! definitions of the model: journalled with the command that made them and saved in the +//! `.fieldcraft` file. + +use serde::{Deserialize, Serialize}; + +use crate::{ElementKind, Error, Result, check_coordinates}; + +/// Most cells one field holds. +pub const MAX_RELAXATION_CELLS: usize = 1_000_000; + +/// One cell of a field: a linear element's corners before and after. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct RelaxationCell { + /// tri3, quad4, tet4, wedge6 or hex8 (a quadratic element's corners). + pub kind: ElementKind, + pub at: Vec<[f64; 3]>, + pub moved: Vec<[f64; 3]>, +} + +/// A relaxation field. +#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)] +#[serde(deny_unknown_fields)] +pub struct Relaxation { + /// The bodies it was made on (empty: any). + #[serde(default, skip_serializing_if = "Vec::is_empty")] + pub bodies: Vec, + pub cells: Vec, + /// Applied again by itself after each remesh of its bodies. + #[serde(default)] + pub reapply: bool, +} + +impl Relaxation { + pub fn check(&self) -> Result<()> { + if self.cells.len() > MAX_RELAXATION_CELLS { + return Err(Error::Limit(format!("a relaxation field holds at most {MAX_RELAXATION_CELLS} cells"))); + } + for c in &self.cells { + if !matches!(c.kind, ElementKind::Tri3 | ElementKind::Quad4 | ElementKind::Tet4 | ElementKind::Wedge6 | ElementKind::Hex8) { + return Err(Error::Invalid(format!("a relaxation cell is linear, not {}", c.kind.name()))); + } + let n = c.kind.node_count(); + if c.at.len() != n || c.moved.len() != n { + return Err(Error::Invalid(format!("a {} relaxation cell has {n} corners", c.kind.name()))); + } + check_coordinates(c.at.as_flattened())?; + check_coordinates(c.moved.as_flattened())?; + } + Ok(()) + } +}