forked from 1b-agent/fieldcraft-upstream
mesh: node relaxation with UV surface slides and relaxation fields (#398)
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
This commit is contained in:
jyvikoandClaude Opus 5.5 authored and Irakli Kourtis committed
1 parent
2030a81327
commit
75f9e5ae6e
6 files changed
+1477
No files matched your search
@@ -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;
|
||||
|
||||
@@ -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<Freedom>;
|
||||
/// `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>;
|
||||
/// `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<P>;
|
||||
}
|
||||
|
||||
/// 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<Freedom> {
|
||||
None
|
||||
}
|
||||
fn slide_face(&self, _: u32, _: P, _: P) -> Option<P> {
|
||||
None
|
||||
}
|
||||
fn slide_curve(&self, _: u32, _: P, _: P) -> Option<P> {
|
||||
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<Falloff> {
|
||||
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<u32, P>,
|
||||
/// 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<u32>>,
|
||||
) -> Result<BTreeMap<u32, f64>> {
|
||||
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<u32>, strength: f64) -> Result<BTreeMap<u32, f64>> {
|
||||
check_strength(strength)?;
|
||||
let corners = corner_nodes(model);
|
||||
let out: BTreeMap<u32, f64> = 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<u32> {
|
||||
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<u32>) -> BTreeSet<u32> {
|
||||
let mut faces: BTreeMap<Vec<u32>, (usize, Vec<u32>)> = 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<u32> = f.iter().filter_map(|&i| e.nodes.get(i).copied()).collect();
|
||||
let mut key: Vec<u32> = 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<u32>,
|
||||
/// A shell's normal as it was: its side, which a move must keep.
|
||||
normal: Option<P>,
|
||||
}
|
||||
|
||||
/// The shell normal from its corners (`None` for solids).
|
||||
fn shell_normal(kind: ElementKind, x: &[P]) -> Option<P> {
|
||||
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<u32, P>,
|
||||
elems: Vec<Elem>,
|
||||
/// Corner node → its elements (indices into `elems`).
|
||||
around: BTreeMap<u32, Vec<usize>>,
|
||||
/// Corner node → corner nodes it shares an element edge with.
|
||||
neighbours: BTreeMap<u32, BTreeSet<u32>>,
|
||||
/// Mid-side node → its edge's corners and its offset from their midpoint.
|
||||
mids: BTreeMap<u32, (u32, u32, P)>,
|
||||
/// Corner node → the mid-side nodes on its edges.
|
||||
mids_of: BTreeMap<u32, Vec<u32>>,
|
||||
skin: BTreeSet<u32>,
|
||||
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<P> {
|
||||
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<P> {
|
||||
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<bool>) {
|
||||
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<bool>), 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<u32> = 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<u32, f64>, iterations: usize, geometry: &dyn Constraints) -> Result<Relaxed> {
|
||||
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<u32> = 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<ElementKind> {
|
||||
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<RelaxationCell> {
|
||||
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<u32> = e.nodes.iter().take(kind.node_count()).copied().collect();
|
||||
let at: Vec<P> = 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<f64>, 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<P> {
|
||||
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<P> = 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<u32>>,
|
||||
geometry: &dyn Constraints,
|
||||
polish: f64,
|
||||
) -> Result<Relaxed> {
|
||||
let candidates: BTreeSet<u32> = 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<u32, P> = 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<u32>)> = 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<P> { 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;
|
||||
@@ -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<u32> {
|
||||
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<Freedom> {
|
||||
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<P> {
|
||||
Annulus::surface(face).slide(p, d)
|
||||
}
|
||||
fn slide_curve(&self, curve: u32, p: P, d: P) -> Option<P> {
|
||||
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<u32> {
|
||||
let geo = Geo { model: m.clone() };
|
||||
let mut moved = BTreeSet::new();
|
||||
let labels: Vec<u32> = 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<P> = 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<u32> = 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<u32> = t.iter().map(|&i| e.nodes[i]).collect();
|
||||
let x: Vec<P> = 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<u32> = 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<u32> =
|
||||
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<Freedom> {
|
||||
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<P> {
|
||||
self.0.slide_face(f, p, d)
|
||||
}
|
||||
fn slide_curve(&self, c: u32, p: P, d: P) -> Option<P> {
|
||||
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<u32> = 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<u32, f64> = [(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()));
|
||||
}
|
||||
@@ -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<P> {
|
||||
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<P> {
|
||||
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<P> {
|
||||
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<P> {
|
||||
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());
|
||||
}
|
||||
}
|
||||
@@ -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<Seed>,
|
||||
/// Relaxation fields (`mesh.node.relax`), by name.
|
||||
#[serde(skip_serializing_if = "BTreeMap::is_empty")]
|
||||
relaxations: BTreeMap<String, Relaxation>,
|
||||
}
|
||||
|
||||
/// A [`Model`] as stored, before it is checked.
|
||||
@@ -210,6 +215,8 @@ struct ModelData {
|
||||
suppressed_datums: BTreeSet<String>,
|
||||
#[serde(default)]
|
||||
seeds: Vec<Seed>,
|
||||
#[serde(default)]
|
||||
relaxations: BTreeMap<String, Relaxation>,
|
||||
}
|
||||
|
||||
impl TryFrom<ModelData> for Model {
|
||||
@@ -293,6 +300,9 @@ impl TryFrom<ModelData> 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<String, Relaxation> {
|
||||
&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<Relaxation> {
|
||||
remove_existing(&mut self.relaxations, "relaxation", name)
|
||||
}
|
||||
|
||||
/// Counts of each kind of definition, for summaries.
|
||||
pub fn definition_counts(&self) -> [(&'static str, usize); 16] {
|
||||
[
|
||||
|
||||
@@ -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<String>,
|
||||
pub cells: Vec<RelaxationCell>,
|
||||
/// 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(())
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user