forked from 1b-agent/fieldcraft-upstream
mesh: projection blocking — block topology, projection onto the B-rep, elliptic smoothing (#399)
A block topology (carved grid, butterfly and seven-block sphere templates, O-grid, extrude, connect) projected onto a body: outer block faces onto their smooth face groups, edges onto the curve of two targets and vertices onto the point of three by tangent-plane iteration, nodes by transfinite interpolation, then Winslow smoothing with the boundary on the geometry and a compass search on the worst elements. 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
8092c39892
commit
6efad8feeb
8 files changed
+2179
No files matched your search
@@ -16,6 +16,7 @@ pub mod mirror;
|
||||
pub mod nodes;
|
||||
pub mod octree;
|
||||
pub mod planar;
|
||||
pub mod project;
|
||||
pub mod quad;
|
||||
pub mod quality;
|
||||
pub mod relax;
|
||||
|
||||
@@ -0,0 +1,224 @@
|
||||
//! The node lattice of a block topology: each block's structured lattice, with the nodes on
|
||||
//! shared vertices, edges and faces shared between the blocks.
|
||||
//!
|
||||
//! Division counts are per class of parallel block edges (opposite edges of a block have the
|
||||
//! same count, and so every edge chained to them through other blocks): the classes are the
|
||||
//! connected components of "parallel in some block". A node on a block's boundary is named by
|
||||
//! the topology vertices it is blended from and its exact (rational) bilinear weights, which do
|
||||
//! not depend on which block or orientation sees it, so the blocks around a face, edge or vertex
|
||||
//! find the same node.
|
||||
|
||||
use std::collections::{BTreeMap, HashMap};
|
||||
|
||||
use super::topology::{CORNERS, FACE_SIDES, Topology, corner_at};
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// The topological entity a node lies inside.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub(crate) enum Entity {
|
||||
Vertex(u32),
|
||||
/// Ascending vertex pair.
|
||||
Edge([u32; 2]),
|
||||
/// Ascending corner vertices.
|
||||
Face([u32; 4]),
|
||||
Inside,
|
||||
}
|
||||
|
||||
/// The ascending corner vertices of block face `(d, s)`.
|
||||
pub(crate) fn face_key(block: &[u32; 8], d: usize, s: usize) -> [u32; 4] {
|
||||
let mut k: Vec<u32> = CORNERS.iter().zip(block).filter(|(c, _)| c[d] == s).map(|(_, v)| *v).collect();
|
||||
k.sort_unstable();
|
||||
[0, 1, 2, 3].map(|i| k.get(i).copied().unwrap_or(u32::MAX))
|
||||
}
|
||||
|
||||
/// The two corners (block-local) of each block edge along direction `d`.
|
||||
fn edges_along(d: usize) -> Vec<[usize; 2]> {
|
||||
CORNERS
|
||||
.iter()
|
||||
.enumerate()
|
||||
.filter(|(_, c)| c[d] == 0)
|
||||
.map(|(k, c)| {
|
||||
let mut e = *c;
|
||||
e[d] = 1;
|
||||
[k, corner_at(e)]
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn find(parent: &mut [usize], mut x: usize) -> usize {
|
||||
while parent[x] != x {
|
||||
parent[x] = parent[parent[x]];
|
||||
x = parent[x];
|
||||
}
|
||||
x
|
||||
}
|
||||
|
||||
/// The class of parallel edges of each block direction, and the number of classes.
|
||||
pub(crate) fn edge_classes(t: &Topology) -> (Vec<[usize; 3]>, usize) {
|
||||
let mut ids: HashMap<[u32; 2], usize> = HashMap::new();
|
||||
let mut id_of = |a: u32, b: u32| {
|
||||
let k = [a.min(b), a.max(b)];
|
||||
let n = ids.len();
|
||||
*ids.entry(k).or_insert(n)
|
||||
};
|
||||
let mut block_edges: Vec<[[usize; 4]; 3]> = Vec::with_capacity(t.blocks.len());
|
||||
for block in &t.blocks {
|
||||
let mut per = [[0; 4]; 3];
|
||||
for (d, row) in per.iter_mut().enumerate() {
|
||||
for (slot, [a, b]) in row.iter_mut().zip(edges_along(d)) {
|
||||
*slot = id_of(block[a], block[b]);
|
||||
}
|
||||
}
|
||||
block_edges.push(per);
|
||||
}
|
||||
let mut parent: Vec<usize> = (0..ids.len()).collect();
|
||||
for per in &block_edges {
|
||||
for row in per {
|
||||
for &e in &row[1..] {
|
||||
let (a, b) = (find(&mut parent, row[0]), find(&mut parent, e));
|
||||
if a != b {
|
||||
parent[a.max(b)] = a.min(b);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut class: BTreeMap<usize, usize> = BTreeMap::new();
|
||||
let out = block_edges
|
||||
.iter()
|
||||
.map(|per| {
|
||||
[0, 1, 2].map(|d| {
|
||||
let r = find(&mut parent, per[d][0]);
|
||||
let n = class.len();
|
||||
*class.entry(r).or_insert(n)
|
||||
})
|
||||
})
|
||||
.collect();
|
||||
(out, class.len())
|
||||
}
|
||||
|
||||
/// The structured lattices of the blocks, joined.
|
||||
#[derive(Debug, Clone)]
|
||||
pub(crate) struct Lattice {
|
||||
pub divisions: Vec<[usize; 3]>,
|
||||
/// Per block, the global node of each lattice node ([`crate::quad::lattice_index`]).
|
||||
pub block_nodes: Vec<Vec<u32>>,
|
||||
pub entity: Vec<Entity>,
|
||||
pub hexes: Vec<[u32; 8]>,
|
||||
/// The block and cell of each hexahedron.
|
||||
pub hex_cell: Vec<(u32, [usize; 3])>,
|
||||
}
|
||||
|
||||
fn gcd(mut a: u64, mut b: u64) -> u64 {
|
||||
while b != 0 {
|
||||
(a, b) = (b, a % b);
|
||||
}
|
||||
a.max(1)
|
||||
}
|
||||
|
||||
/// Build the joined lattices; `divisions[b]` are block `b`'s counts per direction (equal across
|
||||
/// shared edges, [`edge_classes`]). At most `max_elements` hexahedra.
|
||||
pub(crate) fn build(t: &Topology, divisions: &[[usize; 3]], max_elements: usize) -> Result<Lattice> {
|
||||
if divisions.len() != t.blocks.len() {
|
||||
return Err(Error::Invalid("one set of divisions per block".into()));
|
||||
}
|
||||
let mut total = 0usize;
|
||||
for n in divisions {
|
||||
if n.iter().any(|&x| x == 0 || x > crate::block::MAX_DIVISIONS) {
|
||||
return Err(Error::Invalid(format!("divisions {n:?} out of range")));
|
||||
}
|
||||
total = total.saturating_add(n[0].saturating_mul(n[1]).saturating_mul(n[2]));
|
||||
}
|
||||
if total > max_elements {
|
||||
return Err(Error::Limit(format!("{total} elements; at most {max_elements} (give a larger size)")));
|
||||
}
|
||||
let mut keys: HashMap<Vec<(u32, u64, u64)>, u32> = HashMap::new();
|
||||
let mut entity: Vec<Entity> = Vec::new();
|
||||
let mut block_nodes = Vec::with_capacity(t.blocks.len());
|
||||
let mut hexes = Vec::with_capacity(total);
|
||||
let mut hex_cell = Vec::with_capacity(total);
|
||||
for (b, (block, n)) in t.blocks.iter().zip(divisions).enumerate() {
|
||||
let size = (n[0] + 1) * (n[1] + 1) * (n[2] + 1);
|
||||
let mut nodes = Vec::with_capacity(size);
|
||||
for k in 0..=n[2] {
|
||||
for j in 0..=n[1] {
|
||||
for i in 0..=n[0] {
|
||||
let at = [i, j, k];
|
||||
let fixed: Vec<(usize, usize)> = (0..3)
|
||||
.filter_map(|d| {
|
||||
if at[d] == 0 {
|
||||
Some((d, 0))
|
||||
} else if at[d] == n[d] {
|
||||
Some((d, 1))
|
||||
} else {
|
||||
None
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
if fixed.is_empty() {
|
||||
nodes.push(entity.len() as u32);
|
||||
entity.push(Entity::Inside);
|
||||
continue;
|
||||
}
|
||||
let free: Vec<usize> = (0..3).filter(|d| !fixed.iter().any(|f| f.0 == *d)).collect();
|
||||
let den: u64 = free.iter().map(|&d| n[d] as u64).product();
|
||||
let mut key: Vec<(u32, u64, u64)> = Vec::new();
|
||||
for (c, bits) in CORNERS.iter().enumerate() {
|
||||
if fixed.iter().any(|&(d, s)| bits[d] != s) {
|
||||
continue;
|
||||
}
|
||||
let num: u64 = free.iter().map(|&d| if bits[d] == 1 { at[d] as u64 } else { (n[d] - at[d]) as u64 }).product();
|
||||
if num == 0 {
|
||||
continue;
|
||||
}
|
||||
let g = gcd(num, den);
|
||||
key.push((block[c], num / g, den / g));
|
||||
}
|
||||
key.sort_unstable();
|
||||
let next = entity.len() as u32;
|
||||
let id = match keys.get(&key) {
|
||||
Some(&id) => id,
|
||||
None => {
|
||||
let mut vs: Vec<u32> = key.iter().map(|x| x.0).collect();
|
||||
vs.sort_unstable();
|
||||
let e = match vs.len() {
|
||||
1 => Entity::Vertex(vs[0]),
|
||||
2 => Entity::Edge([vs[0], vs[1]]),
|
||||
_ => Entity::Face([0, 1, 2, 3].map(|q| vs.get(q).copied().unwrap_or(u32::MAX))),
|
||||
};
|
||||
entity.push(e);
|
||||
keys.insert(key, next);
|
||||
next
|
||||
}
|
||||
};
|
||||
nodes.push(id);
|
||||
}
|
||||
}
|
||||
}
|
||||
let idx = |c: [usize; 3]| crate::quad::lattice_index(*n, c);
|
||||
for k in 0..n[2] {
|
||||
for j in 0..n[1] {
|
||||
for i in 0..n[0] {
|
||||
let h = CORNERS.map(|c| nodes.get(idx([i + c[0], j + c[1], k + c[2]])).copied().unwrap_or(0));
|
||||
hexes.push(h);
|
||||
hex_cell.push((b as u32, [i, j, k]));
|
||||
}
|
||||
}
|
||||
}
|
||||
block_nodes.push(nodes);
|
||||
}
|
||||
Ok(Lattice { divisions: divisions.to_vec(), block_nodes, entity, hexes, hex_cell })
|
||||
}
|
||||
|
||||
/// How many blocks use each block face.
|
||||
pub(crate) fn face_uses(t: &Topology) -> Result<HashMap<[u32; 4], usize>> {
|
||||
let mut uses: HashMap<[u32; 4], usize> = HashMap::new();
|
||||
for block in &t.blocks {
|
||||
for (d, s) in FACE_SIDES {
|
||||
*uses.entry(face_key(block, d, s)).or_insert(0) += 1;
|
||||
}
|
||||
}
|
||||
if uses.values().any(|&u| u > 2) {
|
||||
return Err(Error::Invalid("a block face is shared by more than two blocks".into()));
|
||||
}
|
||||
Ok(uses)
|
||||
}
|
||||
@@ -0,0 +1,621 @@
|
||||
//! Projection blocking: a block topology laid over the part and projected onto its B-rep.
|
||||
//!
|
||||
//! Instead of cutting the geometry into sweepable cells, the analyst lays down a clean block
|
||||
//! topology — a grid with blocks carved away, butterfly (O-grid) templates for round sections,
|
||||
//! the seven-block sphere — whose outer faces stand for the part's surface. Every block face on
|
||||
//! the outside of the topology is projected onto a target: one B-rep face, or a union of faces
|
||||
//! (by default the smooth group of faces nearest it). Where two outer block faces with different
|
||||
//! targets meet, their edge goes onto the curve where the targets meet, and a vertex on three
|
||||
//! targets onto their common point ([`target`]). Nodes are then placed outside in: block
|
||||
//! vertices, edge nodes (evenly by arc length along their curve), face nodes by transfinite
|
||||
//! interpolation of their edges projected onto the face's target, and interior nodes by the
|
||||
//! transfinite interpolation of the block's faces ([`crate::quad::transfinite_volume`], Gordon
|
||||
//! and Hall 1973, Gordon and Thiel 1982). Elliptic smoothing ([`smooth`]) finishes the grid with
|
||||
//! the boundary held on the geometry.
|
||||
//!
|
||||
//! The method follows the public descriptions of block-structured projection meshing (blocks in
|
||||
//! index space, carved, projected face → surface, edge → curve, vertex → point, then smoothed);
|
||||
//! the algorithms are the papers cited in each module.
|
||||
|
||||
mod lattice;
|
||||
pub mod smooth;
|
||||
pub mod target;
|
||||
pub mod topology;
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet, HashMap};
|
||||
|
||||
use fieldcraft_model::ElementKind;
|
||||
|
||||
pub use smooth::{QUALITY_FLOOR, SmoothReport, Smoothing};
|
||||
pub use target::Geometry;
|
||||
pub use topology::Topology;
|
||||
|
||||
use crate::hex::BoundaryFace;
|
||||
use crate::quality::{Metric, element_quality};
|
||||
use crate::surface::{BoundarySurface, FeatureKind};
|
||||
use crate::vec3::{P, add, dist, lerp, scale};
|
||||
use crate::{Error, Result};
|
||||
use lattice::{Entity, build, edge_classes, face_key, face_uses};
|
||||
use target::{Constraint, on_target, place};
|
||||
use topology::{CORNERS, FACE_SIDES, corner_at, trilinear};
|
||||
|
||||
/// Faces meeting at less than this angle (degrees) along their shared edge belong to one smooth
|
||||
/// group, the default target of a block face.
|
||||
pub const SMOOTH_ANGLE: f64 = 20.0;
|
||||
|
||||
/// The default [`ProjectOptions::min_quality`]: a corner flattened toward 180° (a block corner
|
||||
/// projected onto a smooth surface) scores far below it.
|
||||
pub const DEFAULT_MIN_QUALITY: f64 = 0.2;
|
||||
|
||||
/// Sweeps of the inner vertices' displacement blend.
|
||||
const VERTEX_SWEEPS: usize = 500;
|
||||
|
||||
/// Rounds of arc-length redistribution of an edge's nodes.
|
||||
const EDGE_ROUNDS: usize = 4;
|
||||
|
||||
/// How to project.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct ProjectOptions {
|
||||
/// The element size: each class of parallel block edges gets `round(length / size)`
|
||||
/// divisions (at least one), from its longest projected edge.
|
||||
pub size: f64,
|
||||
/// Division counts set by hand: `(block, direction 0..3, count)`; the count applies to the
|
||||
/// whole class of parallel edges.
|
||||
pub divisions: Vec<(usize, usize, usize)>,
|
||||
/// Targets set by hand: `(block, hex8 face, B-rep faces)`. Other outer faces take the smooth
|
||||
/// group of B-rep faces nearest them.
|
||||
pub faces: Vec<(usize, u8, Vec<u32>)>,
|
||||
pub smoothing: Smoothing,
|
||||
pub iterations: usize,
|
||||
pub max_elements: usize,
|
||||
/// The smallest scaled Jacobian accepted: a mesh with a worse element is an error that
|
||||
/// says where a template is missing (0 accepts any valid mesh).
|
||||
pub min_quality: f64,
|
||||
}
|
||||
|
||||
impl Default for ProjectOptions {
|
||||
fn default() -> Self {
|
||||
ProjectOptions {
|
||||
size: 1.0,
|
||||
divisions: Vec::new(),
|
||||
faces: Vec::new(),
|
||||
smoothing: Smoothing::Elliptic,
|
||||
iterations: 100,
|
||||
max_elements: 2_000_000,
|
||||
min_quality: DEFAULT_MIN_QUALITY,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// How it went.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct ProjectReport {
|
||||
pub blocks: usize,
|
||||
pub elements: usize,
|
||||
/// The smallest scaled Jacobian (corners) of the transfinite grid, before smoothing.
|
||||
pub initial_min_scaled_jacobian: f64,
|
||||
/// Certified smallest and mean scaled Jacobian of the result.
|
||||
pub min_scaled_jacobian: f64,
|
||||
pub mean_scaled_jacobian: f64,
|
||||
pub smoothing: SmoothReport,
|
||||
/// Farthest a boundary node lies from its target (surface, curve or point).
|
||||
pub max_gap: f64,
|
||||
pub volume: f64,
|
||||
pub surface_volume: f64,
|
||||
/// The B-rep faces each target stands for.
|
||||
pub targets: Vec<Vec<u32>>,
|
||||
/// Divisions per block.
|
||||
pub divisions: Vec<[usize; 3]>,
|
||||
}
|
||||
|
||||
/// An all-hex mesh of a projected topology.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct ProjectedMesh {
|
||||
pub nodes: Vec<P>,
|
||||
pub hexes: Vec<[u32; 8]>,
|
||||
pub boundary: Vec<BoundaryFace>,
|
||||
pub report: ProjectReport,
|
||||
}
|
||||
|
||||
/// The face groups joined across smooth B-rep edges.
|
||||
fn smooth_groups(surface: &BoundarySurface, faces: usize) -> Vec<usize> {
|
||||
let mut parent: Vec<usize> = (0..faces).collect();
|
||||
fn root(p: &mut [usize], mut x: usize) -> usize {
|
||||
while p[x] != x {
|
||||
p[x] = p[p[x]];
|
||||
x = p[x];
|
||||
}
|
||||
x
|
||||
}
|
||||
let tf = surface.tri_face();
|
||||
for e in surface.feature_edges() {
|
||||
if e.kind != FeatureKind::FaceBoundary || e.angle >= SMOOTH_ANGLE {
|
||||
continue;
|
||||
}
|
||||
let (Some(&a), Some(&b)) = (tf.get(e.triangles[0] as usize), tf.get(e.triangles[1] as usize)) else { continue };
|
||||
let (a, b) = (a as usize, b as usize);
|
||||
if a < faces && b < faces {
|
||||
let (ra, rb) = (root(&mut parent, a), root(&mut parent, b));
|
||||
if ra != rb {
|
||||
parent[ra.max(rb)] = ra.min(rb);
|
||||
}
|
||||
}
|
||||
}
|
||||
(0..faces).map(|f| root(&mut parent, f)).collect()
|
||||
}
|
||||
|
||||
/// The B-rep face a sample `u` of block face `(d, s)` stands for: the nearest face that faces the
|
||||
/// same way (normals within 60°), else the nearest face.
|
||||
fn facing(g: &dyn Geometry, faces: &[u32], c: &[P; 8], d: usize, s: usize, u: [f64; 3]) -> Option<u32> {
|
||||
let p = trilinear(c, u);
|
||||
let (e1, e2) = other(d);
|
||||
let h = 1e-3;
|
||||
let diff = |k: usize| {
|
||||
let (mut a, mut b) = (u, u);
|
||||
a[k] -= h;
|
||||
b[k] += h;
|
||||
crate::vec3::sub(trilinear(c, b), trilinear(c, a))
|
||||
};
|
||||
let mut n = crate::vec3::cross(diff(e1), diff(e2));
|
||||
// Outward: along +d on side 1, against it on side 0.
|
||||
let inward = {
|
||||
let mut q = u;
|
||||
q[d] = if s == 1 { 1.0 - 0.1 } else { 0.1 };
|
||||
crate::vec3::sub(trilinear(c, q), p)
|
||||
};
|
||||
if crate::vec3::dot(n, inward) > 0.0 {
|
||||
n = scale(n, -1.0);
|
||||
}
|
||||
let ln = crate::vec3::norm(n);
|
||||
let mut best: Option<(f64, u32)> = None;
|
||||
let mut nearest: Option<(f64, u32)> = None;
|
||||
for &f in faces {
|
||||
let Some(q) = g.project(f, p) else { continue };
|
||||
let dq = dist(p, q);
|
||||
if nearest.is_none_or(|b| dq < b.0) {
|
||||
nearest = Some((dq, f));
|
||||
}
|
||||
let agrees = ln > 0.0 && g.normal(f, q).is_some_and(|m| crate::vec3::dot(m, n) > 0.5 * ln);
|
||||
if agrees && best.is_none_or(|b| dq < b.0) {
|
||||
best = Some((dq, f));
|
||||
}
|
||||
}
|
||||
best.or(nearest).map(|b| b.1)
|
||||
}
|
||||
|
||||
/// The corners of block face `(d, s)` in cyclic order.
|
||||
fn face_cycle(block: &[u32; 8], d: usize, s: usize) -> [u32; 4] {
|
||||
let (e1, e2) = other(d);
|
||||
[[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, b]| {
|
||||
let mut bits = [0; 3];
|
||||
bits[d] = s;
|
||||
bits[e1] = a;
|
||||
bits[e2] = b;
|
||||
block[corner_at(bits)]
|
||||
})
|
||||
}
|
||||
|
||||
pub(crate) fn other(d: usize) -> (usize, usize) {
|
||||
match d {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
}
|
||||
}
|
||||
|
||||
fn edge_key(a: u32, b: u32) -> [u32; 2] {
|
||||
[a.min(b), a.max(b)]
|
||||
}
|
||||
|
||||
/// Mesh the topology `t` projected onto the body bounded by `surface`, through `g`.
|
||||
pub fn project_mesh(t: &Topology, surface: &BoundarySurface, g: &dyn Geometry, opts: &ProjectOptions) -> Result<ProjectedMesh> {
|
||||
t.check()?;
|
||||
if !(opts.size > 0.0 && opts.size.is_finite()) {
|
||||
return Err(Error::Invalid(format!("the element size must be a positive number, got {}", opts.size)));
|
||||
}
|
||||
let uses = face_uses(t)?;
|
||||
let face_count = surface.tri_face().iter().max().map_or(1, |m| *m as usize + 1);
|
||||
let groups = smooth_groups(surface, face_count);
|
||||
let diag = t.diagonal().max(opts.size);
|
||||
let tol = 1e-9 * diag;
|
||||
|
||||
// Targets of the outer block faces.
|
||||
let explicit: BTreeMap<(usize, u8), Vec<u32>> = opts.faces.iter().map(|(b, f, s)| ((*b, *f), s.clone())).collect();
|
||||
for ((b, f), set) in &explicit {
|
||||
let (d, s) = FACE_SIDES
|
||||
.get(usize::from(*f).wrapping_sub(1))
|
||||
.copied()
|
||||
.ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face (1 to 6)")))?;
|
||||
let block = t.blocks.get(*b).ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
|
||||
if uses.get(&face_key(block, d, s)) != Some(&1) {
|
||||
return Err(Error::Invalid(format!("face {f} of block {b} is inside the topology; only outer faces are projected")));
|
||||
}
|
||||
if set.is_empty() || set.iter().any(|&x| x as usize >= face_count) {
|
||||
return Err(Error::Invalid(format!("face {f} of block {b}: the target must name faces of the body")));
|
||||
}
|
||||
}
|
||||
let mut targets: Vec<Vec<u32>> = Vec::new();
|
||||
let mut target_of = |set: Vec<u32>| -> usize {
|
||||
match targets.iter().position(|t| *t == set) {
|
||||
Some(i) => i,
|
||||
None => {
|
||||
targets.push(set);
|
||||
targets.len() - 1
|
||||
}
|
||||
}
|
||||
};
|
||||
let all_faces: Vec<u32> = (0..face_count as u32).collect();
|
||||
let mut face_target: HashMap<[u32; 4], usize> = HashMap::new();
|
||||
let mut outer: Vec<(usize, usize, usize)> = Vec::new();
|
||||
for (b, block) in t.blocks.iter().enumerate() {
|
||||
let c = t.corners(b)?;
|
||||
for (fi, (d, s)) in FACE_SIDES.iter().copied().enumerate() {
|
||||
let key = face_key(block, d, s);
|
||||
if uses.get(&key) != Some(&1) {
|
||||
continue;
|
||||
}
|
||||
outer.push((b, d, s));
|
||||
let set = match explicit.get(&(b, fi as u8 + 1)) {
|
||||
Some(set) => {
|
||||
let mut set = set.clone();
|
||||
set.sort_unstable();
|
||||
set.dedup();
|
||||
set
|
||||
}
|
||||
None => {
|
||||
// The smooth group most samples of the face lie nearest.
|
||||
let (e1, e2) = other(d);
|
||||
let mut votes: BTreeMap<usize, usize> = BTreeMap::new();
|
||||
for (a, bb) in [(0.5, 0.5), (0.25, 0.25), (0.75, 0.25), (0.25, 0.75), (0.75, 0.75)] {
|
||||
let mut u = [0.0; 3];
|
||||
u[d] = s as f64;
|
||||
u[e1] = a;
|
||||
u[e2] = bb;
|
||||
if let Some(f) = facing(g, &all_faces, &c, d, s, u) {
|
||||
*votes.entry(groups.get(f as usize).copied().unwrap_or(0)).or_insert(0) += 1;
|
||||
}
|
||||
}
|
||||
let group = votes
|
||||
.iter()
|
||||
.max_by_key(|(gr, n)| (**n, std::cmp::Reverse(**gr)))
|
||||
.map(|(gr, _)| *gr)
|
||||
.ok_or_else(|| Error::Failed(format!("face {} of block {b} has no B-rep face to project onto", fi + 1)))?;
|
||||
(0..face_count as u32).filter(|&f| groups.get(f as usize) == Some(&group)).collect()
|
||||
}
|
||||
};
|
||||
face_target.insert(key, target_of(set));
|
||||
}
|
||||
}
|
||||
// Targets of the vertices and edges on the outside.
|
||||
let mut vertex_on: BTreeMap<u32, BTreeSet<usize>> = BTreeMap::new();
|
||||
let mut edge_on: BTreeMap<[u32; 2], BTreeSet<usize>> = BTreeMap::new();
|
||||
for &(b, d, s) in &outer {
|
||||
let block = t.blocks.get(b).copied().unwrap_or([0; 8]);
|
||||
let tg = face_target.get(&face_key(&block, d, s)).copied().unwrap_or(0);
|
||||
let cyc = face_cycle(&block, d, s);
|
||||
for i in 0..4 {
|
||||
vertex_on.entry(cyc[i]).or_default().insert(tg);
|
||||
edge_on.entry(edge_key(cyc[i], cyc[(i + 1) % 4])).or_default().insert(tg);
|
||||
}
|
||||
}
|
||||
let as_vec = |s: Option<&BTreeSet<usize>>| s.map(|s| s.iter().copied().collect::<Vec<_>>()).unwrap_or_default();
|
||||
let vertex_c = |v: u32| Constraint::of(&as_vec(vertex_on.get(&v)));
|
||||
let edge_c = |e: [u32; 2]| Constraint::of(&as_vec(edge_on.get(&e)));
|
||||
|
||||
// 1. Vertices.
|
||||
let mut vpos: Vec<P> = t.vertices.clone();
|
||||
for (v, p) in vpos.iter_mut().enumerate() {
|
||||
*p = place(g, &targets, vertex_c(v as u32), *p, tol);
|
||||
}
|
||||
// Inner vertices follow: the projected vertices' moves blended through the block edges
|
||||
// (a bounded Jacobi solve of Laplace's equation for the displacement), so a block the
|
||||
// projection stretches does not leave its inner corners behind.
|
||||
let mut nbr: Vec<BTreeSet<u32>> = vec![BTreeSet::new(); vpos.len()];
|
||||
for block in &t.blocks {
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
for d in 0..3 {
|
||||
let mut e = *bits;
|
||||
e[d] = 1 - e[d];
|
||||
if let Some(set) = nbr.get_mut(block[k] as usize) {
|
||||
set.insert(block[corner_at(e)]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let free: Vec<bool> = (0..vpos.len()).map(|v| vertex_c(v as u32) == Constraint::Free).collect();
|
||||
if free.iter().any(|f| *f) {
|
||||
let mut disp: Vec<P> = vpos.iter().zip(&t.vertices).map(|(a, b)| crate::vec3::sub(*a, *b)).collect();
|
||||
for _ in 0..VERTEX_SWEEPS {
|
||||
let mut change = 0.0f64;
|
||||
for v in 0..disp.len() {
|
||||
if !free[v] {
|
||||
continue;
|
||||
}
|
||||
let Some(ns) = nbr.get(v).filter(|n| !n.is_empty()) else { continue };
|
||||
let mean = scale(ns.iter().fold([0.0; 3], |a, &w| add(a, disp[w as usize])), 1.0 / ns.len() as f64);
|
||||
change = change.max(dist(mean, disp[v]));
|
||||
disp[v] = mean;
|
||||
}
|
||||
if change <= tol {
|
||||
break;
|
||||
}
|
||||
}
|
||||
for v in 0..vpos.len() {
|
||||
if free[v] {
|
||||
vpos[v] = add(t.vertices[v], disp[v]);
|
||||
}
|
||||
}
|
||||
}
|
||||
// 2. Divisions per class of parallel edges, from the longest projected edge.
|
||||
let (classes, nclass) = edge_classes(t);
|
||||
let mut class_len = vec![0.0f64; nclass];
|
||||
for (block, cls) in t.blocks.iter().zip(&classes) {
|
||||
for d in 0..3 {
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
if bits[d] != 0 {
|
||||
continue;
|
||||
}
|
||||
let mut e = *bits;
|
||||
e[d] = 1;
|
||||
let (a, b) = (block[k], block[corner_at(e)]);
|
||||
let (pa, pb) = (vpos.get(a as usize).copied().unwrap_or_default(), vpos.get(b as usize).copied().unwrap_or_default());
|
||||
let c = edge_c(edge_key(a, b));
|
||||
let mut len = 0.0;
|
||||
let mut prev = pa;
|
||||
for q in 1..=16 {
|
||||
let x = if q == 16 { pb } else { place(g, &targets, c, lerp(pa, pb, q as f64 / 16.0), tol) };
|
||||
len += dist(prev, x);
|
||||
prev = x;
|
||||
}
|
||||
if let Some(l) = class_len.get_mut(cls[d]) {
|
||||
*l = l.max(len);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut class_n: Vec<usize> = class_len.iter().map(|l| ((l / opts.size).round() as usize).clamp(1, crate::block::MAX_DIVISIONS)).collect();
|
||||
for &(b, d, n) in &opts.divisions {
|
||||
let cls = classes.get(b).and_then(|c| c.get(d)).ok_or_else(|| Error::Invalid(format!("no direction {d} of block {b}")))?;
|
||||
if !(1..=crate::block::MAX_DIVISIONS).contains(&n) {
|
||||
return Err(Error::Invalid(format!("{n} divisions are out of range")));
|
||||
}
|
||||
if let Some(slot) = class_n.get_mut(*cls) {
|
||||
*slot = n;
|
||||
}
|
||||
}
|
||||
let divisions: Vec<[usize; 3]> = classes.iter().map(|c| c.map(|k| class_n.get(k).copied().unwrap_or(1))).collect();
|
||||
let lat = build(t, &divisions, opts.max_elements)?;
|
||||
|
||||
// Each node's constraint.
|
||||
let constraint: Vec<Constraint> = lat
|
||||
.entity
|
||||
.iter()
|
||||
.map(|e| match e {
|
||||
Entity::Vertex(v) => vertex_c(*v),
|
||||
Entity::Edge(k) => edge_c(*k),
|
||||
Entity::Face(k) => face_target.get(k).map_or(Constraint::Free, |&tg| Constraint::Surface(tg)),
|
||||
Entity::Inside => Constraint::Free,
|
||||
})
|
||||
.collect();
|
||||
let n = lat.entity.len();
|
||||
let mut nodes = vec![[0.0; 3]; n];
|
||||
let mut placed = vec![false; n];
|
||||
for (i, e) in lat.entity.iter().enumerate() {
|
||||
if let Entity::Vertex(v) = e {
|
||||
nodes[i] = vpos.get(*v as usize).copied().unwrap_or_default();
|
||||
placed[i] = true;
|
||||
}
|
||||
}
|
||||
// 3. Edge nodes: even by arc length along their curve or surface.
|
||||
for (b, block) in t.blocks.iter().enumerate() {
|
||||
let nd = lat.divisions[b];
|
||||
let bn = &lat.block_nodes[b];
|
||||
for d in 0..3 {
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
if bits[d] != 0 {
|
||||
continue;
|
||||
}
|
||||
let mut e = *bits;
|
||||
e[d] = 1;
|
||||
let c = edge_c(edge_key(block[k], block[corner_at(e)]));
|
||||
let ids: Vec<u32> = (0..=nd[d])
|
||||
.map(|q| {
|
||||
let mut at = [bits[0] * nd[0], bits[1] * nd[1], bits[2] * nd[2]];
|
||||
at[d] = q;
|
||||
bn.get(crate::quad::lattice_index(nd, at)).copied().unwrap_or(0)
|
||||
})
|
||||
.collect();
|
||||
if ids.iter().skip(1).take(nd[d].saturating_sub(1)).all(|&i| placed[i as usize]) {
|
||||
continue;
|
||||
}
|
||||
let (pa, pb) = (nodes[ids[0] as usize], nodes[ids[nd[d]] as usize]);
|
||||
let m = nd[d];
|
||||
let mut line: Vec<P> = (0..=m).map(|q| place(g, &targets, c, lerp(pa, pb, q as f64 / m as f64), tol)).collect();
|
||||
line[0] = pa;
|
||||
line[m] = pb;
|
||||
if c != Constraint::Free {
|
||||
for _ in 0..EDGE_ROUNDS {
|
||||
line = even(&line).into_iter().map(|x| place(g, &targets, c, x, tol)).collect();
|
||||
line[0] = pa;
|
||||
line[m] = pb;
|
||||
}
|
||||
}
|
||||
for (q, &id) in ids.iter().enumerate() {
|
||||
if !placed[id as usize] {
|
||||
nodes[id as usize] = line[q];
|
||||
placed[id as usize] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// 4. Faces (projected when outside), then 5. the interior.
|
||||
for (b, block) in t.blocks.iter().enumerate() {
|
||||
let nd = lat.divisions[b];
|
||||
let bn = &lat.block_nodes[b];
|
||||
let mut x: Vec<Option<P>> = bn.iter().map(|&i| placed[i as usize].then(|| nodes[i as usize])).collect();
|
||||
crate::quad::transfinite_faces(&mut x, nd, None, |d, side, grid| {
|
||||
if let Some(&tg) = face_target.get(&face_key(block, d, side)) {
|
||||
for p in grid.iter_mut() {
|
||||
*p = place(g, &targets, Constraint::Surface(tg), *p, tol);
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
})?;
|
||||
for (slot, &i) in x.iter_mut().zip(bn) {
|
||||
if placed[i as usize] {
|
||||
*slot = Some(nodes[i as usize]);
|
||||
} else if let Some(p) = *slot {
|
||||
nodes[i as usize] = p;
|
||||
placed[i as usize] = true;
|
||||
}
|
||||
}
|
||||
let blend = nd.map(|m| (0..=m).map(|q| q as f64 / m as f64).collect::<Vec<f64>>());
|
||||
crate::quad::transfinite_volume(&mut x, nd, [&blend[0], &blend[1], &blend[2]])?;
|
||||
for (slot, &i) in x.iter().zip(bn) {
|
||||
if let (false, Some(p)) = (placed[i as usize], slot) {
|
||||
nodes[i as usize] = *p;
|
||||
placed[i as usize] = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
if placed.iter().any(|p| !p) {
|
||||
return Err(Error::Failed("a node of the block lattice was not placed".into()));
|
||||
}
|
||||
let corners_of = |nodes: &[P], h: &[u32; 8]| h.map(|v| nodes.get(v as usize).copied().unwrap_or_default());
|
||||
let initial_min = lat.hexes.iter().map(|h| smooth::hex_scaled_jacobian(&corners_of(&nodes, h))).fold(f64::INFINITY, f64::min);
|
||||
|
||||
// The outer quads (outward, hex8 face order) and the target each lies on.
|
||||
let hex8_faces = ElementKind::Hex8.faces();
|
||||
let mut quads: Vec<([u32; 4], u32, u8, usize)> = Vec::new();
|
||||
for (e, (hex, (b, cell))) in lat.hexes.iter().zip(&lat.hex_cell).enumerate() {
|
||||
let block = t.blocks.get(*b as usize).copied().unwrap_or([0; 8]);
|
||||
let nd = lat.divisions.get(*b as usize).copied().unwrap_or([1; 3]);
|
||||
for (fi, (d, s)) in FACE_SIDES.iter().copied().enumerate() {
|
||||
if cell[d] != if s == 0 { 0 } else { nd[d] - 1 } {
|
||||
continue;
|
||||
}
|
||||
let Some(&tg) = face_target.get(&face_key(&block, d, s)) else { continue };
|
||||
let local = hex8_faces.get(fi).copied().unwrap_or(&[]);
|
||||
let q = [0, 1, 2, 3].map(|i| local.get(i).and_then(|&l| hex.get(l)).copied().unwrap_or(0));
|
||||
quads.push((q, e as u32, fi as u8 + 1, tg));
|
||||
}
|
||||
}
|
||||
let boundary_quads: Vec<[u32; 4]> = quads.iter().map(|q| q.0).collect();
|
||||
let smoothing = smooth::smooth(
|
||||
smooth::Mesh { nodes: &mut nodes, hexes: &lat.hexes, constraint: &constraint, boundary: &boundary_quads },
|
||||
g,
|
||||
&targets,
|
||||
opts.smoothing,
|
||||
opts.iterations,
|
||||
1e-4 * opts.size,
|
||||
);
|
||||
|
||||
// Certified quality, the gap to the geometry and the tags.
|
||||
let (mut min_sj, mut sum_sj, mut volume) = (f64::INFINITY, 0.0, 0.0);
|
||||
let mut invalid = 0usize;
|
||||
for h in &lat.hexes {
|
||||
let q = element_quality(ElementKind::Hex8, &corners_of(&nodes, h))?;
|
||||
let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0);
|
||||
if q.validity().is_none_or(|v| v != crate::quality::Validity::Valid) {
|
||||
invalid += 1;
|
||||
}
|
||||
min_sj = min_sj.min(sj);
|
||||
sum_sj += sj;
|
||||
volume += q.size;
|
||||
}
|
||||
let mut max_gap = 0.0f64;
|
||||
for (i, c) in constraint.iter().enumerate() {
|
||||
if *c != Constraint::Free {
|
||||
max_gap = max_gap.max(dist(nodes[i], place(g, &targets, *c, nodes[i], tol)));
|
||||
}
|
||||
}
|
||||
let mut boundary = Vec::with_capacity(quads.len());
|
||||
let mut enclosed = 0.0;
|
||||
for (q, e, f, tg) in &quads {
|
||||
let pts = q.map(|v| nodes.get(v as usize).copied().unwrap_or_default());
|
||||
enclosed += crate::hex::face_flux(&pts);
|
||||
let centre = scale(pts.iter().fold([0.0; 3], |a, p| add(a, *p)), 0.25);
|
||||
let set = targets.get(*tg).map_or(&[][..], Vec::as_slice);
|
||||
let brep_face = on_target(g, set, centre).map_or_else(|| set.first().copied().unwrap_or(0), |h| h.2);
|
||||
boundary.push(BoundaryFace { element: *e, face: *f, brep_face });
|
||||
}
|
||||
let surface_volume = surface.volume();
|
||||
let report = ProjectReport {
|
||||
blocks: t.blocks.len(),
|
||||
elements: lat.hexes.len(),
|
||||
initial_min_scaled_jacobian: initial_min,
|
||||
min_scaled_jacobian: min_sj,
|
||||
mean_scaled_jacobian: sum_sj / lat.hexes.len().max(1) as f64,
|
||||
smoothing,
|
||||
max_gap,
|
||||
volume,
|
||||
surface_volume,
|
||||
targets: targets.clone(),
|
||||
divisions: lat.divisions.clone(),
|
||||
};
|
||||
// Where the worst element is, for the error messages.
|
||||
let worst = || {
|
||||
lat.hexes
|
||||
.iter()
|
||||
.zip(&lat.hex_cell)
|
||||
.min_by(|a, b| smooth::hex_scaled_jacobian(&corners_of(&nodes, a.0)).total_cmp(&smooth::hex_scaled_jacobian(&corners_of(&nodes, b.0))))
|
||||
.map_or_else(String::new, |(h, c)| {
|
||||
let x = scale(corners_of(&nodes, h).iter().fold([0.0; 3], |a, p| add(a, *p)), 0.125);
|
||||
format!("worst in block {} near ({:.3}, {:.3}, {:.3})", c.0, x[0], x[1], x[2])
|
||||
})
|
||||
};
|
||||
if invalid > 0 {
|
||||
return Err(Error::Failed(format!(
|
||||
"{invalid} of {} elements are inverted or degenerate after projection (smallest scaled Jacobian {min_sj:.3}, {}): the blocks do not fit the part — add an O-grid or a template where a block corner lands on a curved face",
|
||||
lat.hexes.len(),
|
||||
worst()
|
||||
)));
|
||||
}
|
||||
if min_sj < opts.min_quality {
|
||||
return Err(Error::Failed(format!(
|
||||
"the smallest scaled Jacobian is {min_sj:.3}, below {} ({}): a block corner flattens on a curved face — use an O-grid or a butterfly/sphere template there, or lower `min_quality`",
|
||||
opts.min_quality,
|
||||
worst()
|
||||
)));
|
||||
}
|
||||
if (enclosed - volume).abs() > 1e-6 * volume.abs().max(1e-300) {
|
||||
return Err(Error::Failed(format!(
|
||||
"the elements overlap or leave a gap inside (their volume {volume} differs from the {enclosed} their outer faces enclose): check that block faces meet face to face"
|
||||
)));
|
||||
}
|
||||
if (volume - surface_volume).abs() > 0.2 * surface_volume.abs() {
|
||||
return Err(Error::Failed(format!(
|
||||
"the mesh volume {volume} is far from the body's {surface_volume}: the blocks do not cover the part (every outer block face must stand for part of its surface)"
|
||||
)));
|
||||
}
|
||||
Ok(ProjectedMesh { nodes, hexes: lat.hexes, boundary, report })
|
||||
}
|
||||
|
||||
/// The polyline's points redistributed evenly by arc length (ends kept).
|
||||
fn even(line: &[P]) -> Vec<P> {
|
||||
let m = line.len().saturating_sub(1);
|
||||
let mut acc = vec![0.0; line.len()];
|
||||
for q in 1..line.len() {
|
||||
acc[q] = acc[q - 1] + dist(line[q - 1], line[q]);
|
||||
}
|
||||
let total = acc.last().copied().unwrap_or(0.0);
|
||||
if !(total > 0.0) || m == 0 {
|
||||
return line.to_vec();
|
||||
}
|
||||
let mut out = Vec::with_capacity(line.len());
|
||||
let mut seg = 0;
|
||||
for q in 0..=m {
|
||||
let s = total * q as f64 / m as f64;
|
||||
while seg + 1 < m && acc[seg + 1] < s {
|
||||
seg += 1;
|
||||
}
|
||||
let l = acc[seg + 1] - acc[seg];
|
||||
let t = if l > 0.0 { ((s - acc[seg]) / l).clamp(0.0, 1.0) } else { 0.0 };
|
||||
out.push(lerp(line[seg], line[seg + 1], t));
|
||||
}
|
||||
out
|
||||
}
|
||||
@@ -0,0 +1,401 @@
|
||||
//! Smoothing of a projected block mesh with its boundary held on the geometry.
|
||||
//!
|
||||
//! Interior nodes with the regular neighbourhood of a structured grid (eight hexahedra, six
|
||||
//! edge neighbours in three opposite pairs) move by the elliptic (Winslow) grid equations,
|
||||
//! `Σ g^{ij} x_{ξ^i ξ^j} + Σ g^{ii} φ_i x_{ξ^i} = 0` (A. M. Winslow, *Numerical solution of the
|
||||
//! quasilinear Poisson equation in a nonuniform triangle mesh*, J. Comput. Phys. 2:149–172, 1967;
|
||||
//! the three-dimensional form and its discretisation in J. F. Thompson, Z. U. A. Warsi,
|
||||
//! C. W. Mastin, *Numerical Grid Generation*, North-Holland 1985, ch. 6), solved point by point
|
||||
//! with central differences in the node's own stencil. The optional control functions `φ_i`
|
||||
//! follow P. D. Thomas and J. F. Middlecoff, *Direct control of the grid point distribution in
|
||||
//! meshes generated by elliptic equations*, AIAA J. 18(6):652–656, 1980, in their
|
||||
//! one-dimensional form `φ = −x_ξ·x_ξξ / |x_ξ|²`, evaluated on the algebraic (transfinite) grid
|
||||
//! at every node, so the smoothed grid keeps that grid's spacing along each line. Irregular nodes
|
||||
//! (where three, five or more blocks meet) take the mean of their neighbours (Laplace).
|
||||
//!
|
||||
//! Boundary nodes move along their surface or curve: the mean of their boundary neighbours,
|
||||
//! projected back. Block vertices on three targets stay put.
|
||||
//!
|
||||
//! Every move is checked: it is kept only when the smallest scaled Jacobian of the elements
|
||||
//! around the node stays at or above the lower of its value before and [`QUALITY_FLOOR`], so
|
||||
//! smoothing never makes a good element bad nor a bad one worse.
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet};
|
||||
|
||||
use super::target::{Constraint, Geometry, place};
|
||||
use super::topology::{CORNERS, corner_at};
|
||||
use crate::vec3::{P, add, cross, dist, dot, norm, scale, sub};
|
||||
|
||||
/// Smoothing never takes an element's scaled Jacobian below this (unless it already was).
|
||||
pub const QUALITY_FLOOR: f64 = 0.5;
|
||||
|
||||
/// Rounds of the compass search on the worst elements.
|
||||
const OPTIMIZE_ROUNDS: usize = 20;
|
||||
|
||||
/// The compass search's directions: the axes and the cube diagonals.
|
||||
const COMPASS: [P; 14] = {
|
||||
const S: f64 = 0.577_350_269_189_625_8;
|
||||
[
|
||||
[1.0, 0.0, 0.0],
|
||||
[-1.0, 0.0, 0.0],
|
||||
[0.0, 1.0, 0.0],
|
||||
[0.0, -1.0, 0.0],
|
||||
[0.0, 0.0, 1.0],
|
||||
[0.0, 0.0, -1.0],
|
||||
[S, S, S],
|
||||
[S, S, -S],
|
||||
[S, -S, S],
|
||||
[S, -S, -S],
|
||||
[-S, S, S],
|
||||
[-S, S, -S],
|
||||
[-S, -S, S],
|
||||
[-S, -S, -S],
|
||||
]
|
||||
};
|
||||
|
||||
/// Largest magnitude of a control function.
|
||||
const MAX_CONTROL: f64 = 1.0;
|
||||
|
||||
/// How the interior is smoothed.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub enum Smoothing {
|
||||
/// The transfinite grid as it is.
|
||||
None,
|
||||
/// Each node to the mean of its neighbours.
|
||||
Laplace,
|
||||
/// Winslow's elliptic equations.
|
||||
Elliptic,
|
||||
/// Winslow's equations with Thomas–Middlecoff control functions from the algebraic grid.
|
||||
EllipticControlled,
|
||||
}
|
||||
|
||||
/// The smallest scaled Jacobian over the eight corners of a hexahedron.
|
||||
pub(crate) fn hex_scaled_jacobian(x: &[P; 8]) -> f64 {
|
||||
let mut m = f64::INFINITY;
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
let along = |d: usize| {
|
||||
let mut b = *bits;
|
||||
b[d] = 1 - b[d];
|
||||
let e = sub(x[corner_at(b)], x[k]);
|
||||
if bits[d] == 1 { scale(e, -1.0) } else { e }
|
||||
};
|
||||
let (a, b, c) = (along(0), along(1), along(2));
|
||||
let l = norm(a) * norm(b) * norm(c);
|
||||
let v = if l > 0.0 && l.is_finite() { dot(a, cross(b, c)) / l } else { -1.0 };
|
||||
m = m.min(v);
|
||||
}
|
||||
// Inside: the centre and the six face centres, where a twisted element folds first.
|
||||
for u in [[0.5, 0.5, 0.5], [0.5, 0.5, 0.0], [0.5, 0.5, 1.0], [0.5, 0.0, 0.5], [0.5, 1.0, 0.5], [0.0, 0.5, 0.5], [1.0, 0.5, 0.5]] {
|
||||
let der = |d: usize| {
|
||||
let mut s3 = [0.0; 3];
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
let w: f64 = (0..3)
|
||||
.map(|e| {
|
||||
if e == d {
|
||||
if bits[e] == 1 { 1.0 } else { -1.0 }
|
||||
} else if bits[e] == 1 {
|
||||
u[e]
|
||||
} else {
|
||||
1.0 - u[e]
|
||||
}
|
||||
})
|
||||
.product();
|
||||
s3 = add(s3, scale(x[k], w));
|
||||
}
|
||||
s3
|
||||
};
|
||||
let (a, b, c) = (der(0), der(1), der(2));
|
||||
let l = norm(a) * norm(b) * norm(c);
|
||||
m = m.min(if l > 0.0 && l.is_finite() { dot(a, cross(b, c)) / l } else { -1.0 });
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
/// The structured stencil of a regular interior node: `nb[d] = [minus, plus]` and the face
|
||||
/// diagonals `diag[pair][sa][sb]` for the direction pairs (0, 1), (0, 2), (1, 2).
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
struct Stencil {
|
||||
nb: [[u32; 2]; 3],
|
||||
diag: [[[u32; 2]; 2]; 3],
|
||||
control: [f64; 3],
|
||||
}
|
||||
|
||||
const PAIRS: [(usize, usize); 3] = [(0, 1), (0, 2), (1, 2)];
|
||||
|
||||
fn stencil(node: u32, around: &[u32], hexes: &[[u32; 8]]) -> Option<Stencil> {
|
||||
if around.len() != 8 {
|
||||
return None;
|
||||
}
|
||||
// Per hexahedron: the node's three edge neighbours and three face diagonals.
|
||||
let mut triples = Vec::with_capacity(8);
|
||||
for &h in around {
|
||||
let hex = hexes.get(h as usize)?;
|
||||
let k = hex.iter().position(|&v| v == node)?;
|
||||
let bits = CORNERS[k];
|
||||
let flip = |ds: &[usize]| {
|
||||
let mut b = bits;
|
||||
for &d in ds {
|
||||
b[d] = 1 - b[d];
|
||||
}
|
||||
hex[corner_at(b)]
|
||||
};
|
||||
triples.push(([flip(&[0]), flip(&[1]), flip(&[2])], [flip(&[0, 1]), flip(&[0, 2]), flip(&[1, 2])]));
|
||||
}
|
||||
let neighbours: BTreeSet<u32> = triples.iter().flat_map(|t| t.0).collect();
|
||||
if neighbours.len() != 6 {
|
||||
return None;
|
||||
}
|
||||
let together = |a: u32, b: u32| triples.iter().any(|t| t.0.contains(&a) && t.0.contains(&b));
|
||||
let mut axes: Vec<[u32; 2]> = Vec::new();
|
||||
let mut seen = BTreeSet::new();
|
||||
for &a in &neighbours {
|
||||
if seen.contains(&a) {
|
||||
continue;
|
||||
}
|
||||
let opposite: Vec<u32> = neighbours.iter().copied().filter(|&b| b != a && !together(a, b)).collect();
|
||||
let [b] = opposite.as_slice() else { return None };
|
||||
seen.insert(a);
|
||||
seen.insert(*b);
|
||||
axes.push([a, *b]);
|
||||
}
|
||||
let [a0, a1, a2] = axes.as_slice() else { return None };
|
||||
let nb = [*a0, *a1, *a2];
|
||||
let axis_of = |v: u32| (0..3).find_map(|d| nb[d].iter().position(|&x| x == v).map(|s| (d, s)));
|
||||
let mut diag = [[[u32::MAX; 2]; 2]; 3];
|
||||
for (edge, faces) in &triples {
|
||||
let dirs: Vec<(usize, usize)> = edge.iter().map(|&v| axis_of(v)).collect::<Option<_>>()?;
|
||||
// Hex-local direction l went along global axis dirs[l].0, side dirs[l].1.
|
||||
for (pi, (la, lb)) in PAIRS.iter().enumerate() {
|
||||
let ((da, sa), (db, sb)) = (dirs[*la], dirs[*lb]);
|
||||
let (pair, s1, s2) = match PAIRS.iter().position(|p| *p == (da, db)) {
|
||||
Some(q) => (q, sa, sb),
|
||||
None => (PAIRS.iter().position(|p| *p == (db, da))?, sb, sa),
|
||||
};
|
||||
diag[pair][s1][s2] = faces[pi];
|
||||
}
|
||||
}
|
||||
if diag.iter().flatten().flatten().any(|&v| v == u32::MAX) {
|
||||
return None;
|
||||
}
|
||||
Some(Stencil { nb, diag, control: [0.0; 3] })
|
||||
}
|
||||
|
||||
fn at(x: &[P], v: u32) -> P {
|
||||
x.get(v as usize).copied().unwrap_or([0.0; 3])
|
||||
}
|
||||
|
||||
/// The Winslow update of a regular node.
|
||||
fn winslow(x: &[P], s: &Stencil) -> Option<P> {
|
||||
let xi: [P; 3] = [0, 1, 2].map(|d| scale(sub(at(x, s.nb[d][1]), at(x, s.nb[d][0])), 0.5));
|
||||
let g = |i: usize, j: usize| dot(xi[i], xi[j]);
|
||||
let (g11, g22, g33, g12, g13, g23) = (g(0, 0), g(1, 1), g(2, 2), g(0, 1), g(0, 2), g(1, 2));
|
||||
let cof_diag = [g22 * g33 - g23 * g23, g11 * g33 - g13 * g13, g11 * g22 - g12 * g12];
|
||||
let cof_off = [g13 * g23 - g12 * g33, g12 * g23 - g13 * g22, g13 * g12 - g11 * g23];
|
||||
let den = 2.0 * cof_diag.iter().sum::<f64>();
|
||||
if !(den > 0.0 && den.is_finite()) {
|
||||
return None;
|
||||
}
|
||||
let mut num = [0.0; 3];
|
||||
for (d, cd) in cof_diag.iter().enumerate() {
|
||||
let (m, pl) = (at(x, s.nb[d][0]), at(x, s.nb[d][1]));
|
||||
num = add(num, scale(add(m, pl), *cd));
|
||||
num = add(num, scale(sub(pl, m), 0.5 * cd * s.control[d]));
|
||||
}
|
||||
for (q, c) in cof_off.iter().enumerate() {
|
||||
let dg = &s.diag[q];
|
||||
let cross_d = scale(sub(add(at(x, dg[1][1]), at(x, dg[0][0])), add(at(x, dg[1][0]), at(x, dg[0][1]))), 0.25);
|
||||
num = add(num, scale(cross_d, 2.0 * c));
|
||||
}
|
||||
let out = scale(num, 1.0 / den);
|
||||
out.iter().all(|v| v.is_finite()).then_some(out)
|
||||
}
|
||||
|
||||
/// What the smoother reports.
|
||||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||||
pub struct SmoothReport {
|
||||
pub iterations: usize,
|
||||
/// Moves refused because they would have lowered quality.
|
||||
pub refused: usize,
|
||||
}
|
||||
|
||||
/// Everything smoothing needs about the mesh.
|
||||
pub(crate) struct Mesh<'a> {
|
||||
pub nodes: &'a mut Vec<P>,
|
||||
pub hexes: &'a [[u32; 8]],
|
||||
pub constraint: &'a [Constraint],
|
||||
/// Boundary quads (outward), for the neighbours of boundary nodes.
|
||||
pub boundary: &'a [[u32; 4]],
|
||||
}
|
||||
|
||||
/// Smooth in place: at most `iterations` sweeps, stopping when no node moves more than `tol`.
|
||||
pub(crate) fn smooth(m: Mesh, g: &dyn Geometry, targets: &[Vec<u32>], method: Smoothing, iterations: usize, tol: f64) -> SmoothReport {
|
||||
let n = m.nodes.len();
|
||||
let mut around: Vec<Vec<u32>> = vec![Vec::new(); n];
|
||||
for (h, hex) in m.hexes.iter().enumerate() {
|
||||
for &v in hex {
|
||||
if let Some(a) = around.get_mut(v as usize) {
|
||||
a.push(h as u32);
|
||||
}
|
||||
}
|
||||
}
|
||||
// Edge neighbours (all) and boundary neighbours (along boundary quads).
|
||||
let mut edges: Vec<BTreeSet<u32>> = vec![BTreeSet::new(); n];
|
||||
for hex in m.hexes {
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
for d in 0..3 {
|
||||
let mut b = *bits;
|
||||
b[d] = 1 - b[d];
|
||||
if let Some(e) = edges.get_mut(hex[k] as usize) {
|
||||
e.insert(hex[corner_at(b)]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut along_boundary: BTreeMap<u32, BTreeSet<u32>> = BTreeMap::new();
|
||||
for q in m.boundary {
|
||||
for i in 0..4 {
|
||||
let (a, b) = (q[i], q[(i + 1) % 4]);
|
||||
along_boundary.entry(a).or_default().insert(b);
|
||||
along_boundary.entry(b).or_default().insert(a);
|
||||
}
|
||||
}
|
||||
let cons = |v: u32| m.constraint.get(v as usize).copied().unwrap_or(Constraint::Free);
|
||||
let mut stencils: Vec<Option<Stencil>> = vec![None; n];
|
||||
if matches!(method, Smoothing::Elliptic | Smoothing::EllipticControlled) {
|
||||
for v in 0..n as u32 {
|
||||
if cons(v) == Constraint::Free {
|
||||
let mut s = around.get(v as usize).and_then(|a| stencil(v, a, m.hexes));
|
||||
if let (Some(st), Smoothing::EllipticControlled) = (s.as_mut(), method) {
|
||||
let p = at(m.nodes, v);
|
||||
for d in 0..3 {
|
||||
let (mi, pl) = (at(m.nodes, st.nb[d][0]), at(m.nodes, st.nb[d][1]));
|
||||
let xi = scale(sub(pl, mi), 0.5);
|
||||
let xii = sub(add(pl, mi), scale(p, 2.0));
|
||||
let l2 = dot(xi, xi);
|
||||
st.control[d] = if l2 > 0.0 { (-dot(xi, xii) / l2).clamp(-MAX_CONTROL, MAX_CONTROL) } else { 0.0 };
|
||||
}
|
||||
}
|
||||
stencils[v as usize] = s;
|
||||
}
|
||||
}
|
||||
}
|
||||
let quality = |x: &[P], v: u32| -> f64 {
|
||||
around.get(v as usize).map_or(1.0, |hs| {
|
||||
hs.iter().filter_map(|&h| m.hexes.get(h as usize)).map(|hex| hex_scaled_jacobian(&hex.map(|q| at(x, q)))).fold(f64::INFINITY, f64::min)
|
||||
})
|
||||
};
|
||||
let mut report = SmoothReport { iterations: 0, refused: 0 };
|
||||
if method == Smoothing::None {
|
||||
return report;
|
||||
}
|
||||
let ptol = 1e-3 * tol;
|
||||
for it in 0..iterations {
|
||||
report.iterations = it + 1;
|
||||
let mut moved = 0.0f64;
|
||||
for v in 0..n as u32 {
|
||||
let c = cons(v);
|
||||
let p = at(m.nodes, v);
|
||||
let mean = |nbrs: &mut dyn Iterator<Item = u32>| -> Option<P> {
|
||||
let (mut s, mut k) = ([0.0; 3], 0usize);
|
||||
for w in nbrs {
|
||||
s = add(s, at(m.nodes, w));
|
||||
k += 1;
|
||||
}
|
||||
(k > 0).then(|| scale(s, 1.0 / k as f64))
|
||||
};
|
||||
let target = match c {
|
||||
Constraint::Corner(..) => continue,
|
||||
Constraint::Free => match stencils.get(v as usize).copied().flatten() {
|
||||
Some(s) => winslow(m.nodes, &s),
|
||||
None => edges.get(v as usize).and_then(|e| mean(&mut e.iter().copied())),
|
||||
},
|
||||
Constraint::Surface(_) => along_boundary.get(&v).and_then(|e| mean(&mut e.iter().copied())).map(|q| place(g, targets, c, q, ptol)),
|
||||
Constraint::Curve(a, b) => {
|
||||
let on_curve = |w: &u32| match cons(*w) {
|
||||
Constraint::Curve(x, y) => (x, y) == (a, b),
|
||||
Constraint::Corner(x, y, z) => [x, y, z].contains(&a) && [x, y, z].contains(&b),
|
||||
_ => false,
|
||||
};
|
||||
along_boundary.get(&v).and_then(|e| mean(&mut e.iter().copied().filter(on_curve))).map(|q| place(g, targets, c, q, ptol))
|
||||
}
|
||||
};
|
||||
let Some(q) = target else { continue };
|
||||
if dist(q, p) <= 1e-3 * tol {
|
||||
continue;
|
||||
}
|
||||
// The full move, else half, else a quarter of it, whichever keeps the quality.
|
||||
let before = quality(m.nodes, v);
|
||||
let mut accepted = false;
|
||||
for step in [1.0, 0.5, 0.25] {
|
||||
let x = if step == 1.0 { q } else { place(g, targets, c, add(p, scale(sub(q, p), step)), ptol) };
|
||||
if let Some(slot) = m.nodes.get_mut(v as usize) {
|
||||
*slot = x;
|
||||
}
|
||||
if quality(m.nodes, v) >= before.min(QUALITY_FLOOR) {
|
||||
moved = moved.max(dist(x, p));
|
||||
accepted = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
if !accepted {
|
||||
report.refused += 1;
|
||||
if let Some(slot) = m.nodes.get_mut(v as usize) {
|
||||
*slot = p;
|
||||
}
|
||||
}
|
||||
}
|
||||
if moved <= tol {
|
||||
break;
|
||||
}
|
||||
}
|
||||
// The worst elements: each of their nodes searched for the place that raises the smallest
|
||||
// scaled Jacobian around it (a compass search, the optimization-based smoothing of L. A.
|
||||
// Freitag and C. Ollivier-Gooch, *Tetrahedral mesh improvement using swapping and smoothing*,
|
||||
// IJNME 40:3979–4002, 1997, on Knupp's hexahedral metric), boundary nodes kept on their
|
||||
// targets.
|
||||
for _ in 0..OPTIMIZE_ROUNDS {
|
||||
let poor: BTreeSet<u32> =
|
||||
m.hexes.iter().filter(|h| hex_scaled_jacobian(&h.map(|q| at(m.nodes, q))) < QUALITY_FLOOR).flat_map(|h| h.iter().copied()).collect();
|
||||
if poor.is_empty() {
|
||||
break;
|
||||
}
|
||||
let mut improved = false;
|
||||
for &v in &poor {
|
||||
let c = cons(v);
|
||||
if matches!(c, Constraint::Corner(..)) {
|
||||
continue;
|
||||
}
|
||||
let p = at(m.nodes, v);
|
||||
let reach = edges.get(v as usize).map_or(0.0, |e| e.iter().map(|&w| dist(at(m.nodes, w), p)).fold(f64::INFINITY, f64::min));
|
||||
if !(reach > 0.0 && reach.is_finite()) {
|
||||
continue;
|
||||
}
|
||||
let mut best = (quality(m.nodes, v), p);
|
||||
let mut h = 0.25 * reach;
|
||||
for _ in 0..6 {
|
||||
let start = best.1;
|
||||
for dir in COMPASS {
|
||||
let x = place(g, targets, c, add(start, scale(dir, h)), ptol);
|
||||
if let Some(slot) = m.nodes.get_mut(v as usize) {
|
||||
*slot = x;
|
||||
}
|
||||
let qx = quality(m.nodes, v);
|
||||
if qx > best.0 + 1e-9 {
|
||||
best = (qx, x);
|
||||
}
|
||||
}
|
||||
if best.1 == start {
|
||||
h *= 0.5;
|
||||
}
|
||||
}
|
||||
if let Some(slot) = m.nodes.get_mut(v as usize) {
|
||||
*slot = best.1;
|
||||
}
|
||||
improved |= best.1 != p;
|
||||
}
|
||||
if !improved {
|
||||
break;
|
||||
}
|
||||
}
|
||||
report
|
||||
}
|
||||
@@ -0,0 +1,143 @@
|
||||
//! Where projected nodes go: onto a target (one B-rep face or a union of faces), onto the curve
|
||||
//! where two targets meet, or onto the point where three meet.
|
||||
//!
|
||||
//! A curve or corner node converges by tangent planes: project the node onto each target,
|
||||
//! intersect the targets' tangent planes there, move to the point of that intersection nearest
|
||||
//! the node, repeat. For two transversal surfaces this is Newton's method on the curve; it needs
|
||||
//! no B-rep edge between the faces, so it tolerates small gaps and overlaps in the CAD (the nodes
|
||||
//! end on the curve where the two surfaces, extended, meet). Faces meeting at a shallow angle
|
||||
//! fall back to alternating projections.
|
||||
|
||||
use crate::tet::geometry::Tessellation;
|
||||
use crate::vec3::{P, add, cross, dist, dot, norm, scale, sub};
|
||||
|
||||
/// Rounds of the curve and corner iterations.
|
||||
const ROUNDS: usize = 60;
|
||||
|
||||
/// Below this `|n_a × n_b|` two targets are treated as tangent (sin 3°).
|
||||
const TANGENT_SIN: f64 = 0.052;
|
||||
|
||||
/// What the projection needs of the geometry.
|
||||
pub trait Geometry {
|
||||
/// The closest point on B-rep face `face`.
|
||||
fn project(&self, face: u32, p: P) -> Option<P>;
|
||||
/// The unit normal of B-rep face `face` near `p`.
|
||||
fn normal(&self, face: u32, p: P) -> Option<P>;
|
||||
}
|
||||
|
||||
impl Geometry for Tessellation {
|
||||
fn project(&self, face: u32, p: P) -> Option<P> {
|
||||
self.project_trimmed(face, p)
|
||||
}
|
||||
fn normal(&self, face: u32, p: P) -> Option<P> {
|
||||
self.face_normal(face, p)
|
||||
}
|
||||
}
|
||||
|
||||
/// Where a node may go: targets by index into the projection's target list.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub(crate) enum Constraint {
|
||||
Free,
|
||||
Surface(usize),
|
||||
Curve(usize, usize),
|
||||
Corner(usize, usize, usize),
|
||||
}
|
||||
|
||||
impl Constraint {
|
||||
/// The constraint of a node on the targets `on` (distinct, ascending).
|
||||
pub(crate) fn of(on: &[usize]) -> Constraint {
|
||||
match *on {
|
||||
[] => Constraint::Free,
|
||||
[a] => Constraint::Surface(a),
|
||||
[a, b] => Constraint::Curve(a, b),
|
||||
[a, b, c, ..] => Constraint::Corner(a, b, c),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The closest point to `p` on any face of `faces`, with its normal and face.
|
||||
pub(crate) fn on_target(g: &dyn Geometry, faces: &[u32], p: P) -> Option<(P, P, u32)> {
|
||||
let mut best: Option<(f64, P, u32)> = None;
|
||||
for &f in faces {
|
||||
if let Some(q) = g.project(f, p) {
|
||||
let d = dist(p, q);
|
||||
if best.is_none_or(|b| d < b.0) {
|
||||
best = Some((d, q, f));
|
||||
}
|
||||
}
|
||||
}
|
||||
let (_, q, f) = best?;
|
||||
let n = g.normal(f, q)?;
|
||||
Some((q, n, f))
|
||||
}
|
||||
|
||||
/// Place `p` by its constraint over `targets` (face lists); `tol` is the convergence distance.
|
||||
/// A point the geometry cannot place stays where it is.
|
||||
pub(crate) fn place(g: &dyn Geometry, targets: &[Vec<u32>], c: Constraint, p: P, tol: f64) -> P {
|
||||
let set = |i: usize| targets.get(i).map_or(&[][..], Vec::as_slice);
|
||||
match c {
|
||||
Constraint::Free => p,
|
||||
Constraint::Surface(a) => on_target(g, set(a), p).map_or(p, |h| h.0),
|
||||
Constraint::Curve(a, b) => planes(g, &[set(a), set(b)], p, tol),
|
||||
Constraint::Corner(a, b, c) => planes(g, &[set(a), set(b), set(c)], p, tol),
|
||||
}
|
||||
}
|
||||
|
||||
/// The tangent-plane iteration onto the intersection of two or three targets.
|
||||
fn planes(g: &dyn Geometry, sets: &[&[u32]], p: P, tol: f64) -> P {
|
||||
let mut x = p;
|
||||
for _ in 0..ROUNDS {
|
||||
let mut hits = Vec::with_capacity(sets.len());
|
||||
for s in sets {
|
||||
match on_target(g, s, x) {
|
||||
Some(h) => hits.push(h),
|
||||
None => return x,
|
||||
}
|
||||
}
|
||||
let next = match hits.as_slice() {
|
||||
[(qa, na, _), (qb, nb, _)] => {
|
||||
if norm(cross(*na, *nb)) < TANGENT_SIN {
|
||||
alternate(g, sets, x)
|
||||
} else {
|
||||
// The point of both planes nearest x: x − λa·na − λb·nb.
|
||||
let m = dot(*na, *nb);
|
||||
let (ra, rb) = (dot(*na, sub(x, *qa)), dot(*nb, sub(x, *qb)));
|
||||
let den = 1.0 - m * m;
|
||||
let (la, lb) = ((ra - m * rb) / den, (rb - m * ra) / den);
|
||||
sub(x, add(scale(*na, la), scale(*nb, lb)))
|
||||
}
|
||||
}
|
||||
[(qa, na, _), (qb, nb, _), (qc, nc, _)] => {
|
||||
// The three planes' common point (Cramer's rule).
|
||||
let d = dot(*na, cross(*nb, *nc));
|
||||
if d.abs() < TANGENT_SIN * TANGENT_SIN {
|
||||
alternate(g, sets, x)
|
||||
} else {
|
||||
let (ha, hb, hc) = (dot(*na, *qa), dot(*nb, *qb), dot(*nc, *qc));
|
||||
scale(add(add(scale(cross(*nb, *nc), ha), scale(cross(*nc, *na), hb)), scale(cross(*na, *nb), hc)), 1.0 / d)
|
||||
}
|
||||
}
|
||||
_ => return x,
|
||||
};
|
||||
if !next.iter().all(|v| v.is_finite()) {
|
||||
return x;
|
||||
}
|
||||
let step = dist(next, x);
|
||||
x = next;
|
||||
if step <= tol {
|
||||
break;
|
||||
}
|
||||
}
|
||||
x
|
||||
}
|
||||
|
||||
/// One round of alternating projections onto each target.
|
||||
fn alternate(g: &dyn Geometry, sets: &[&[u32]], p: P) -> P {
|
||||
let mut x = p;
|
||||
for s in sets {
|
||||
if let Some((q, _, _)) = on_target(g, s, x) {
|
||||
x = q;
|
||||
}
|
||||
}
|
||||
x
|
||||
}
|
||||
@@ -0,0 +1,228 @@
|
||||
use std::f64::consts::{PI, TAU};
|
||||
|
||||
use super::target::{Constraint, Geometry, place};
|
||||
use super::*;
|
||||
use crate::surface::{BoundarySurface, SurfaceOptions};
|
||||
use crate::tet::geometry::{AnalyticSurface, Tessellation};
|
||||
|
||||
/// A UV sphere of radius `r` about the origin, one face, with its exact surface.
|
||||
fn sphere(r: f64, n: usize) -> (BoundarySurface, Tessellation) {
|
||||
let rings = n / 2;
|
||||
let mut p = vec![[0.0, 0.0, -r], [0.0, 0.0, r]];
|
||||
for j in 1..rings {
|
||||
let th = PI * j as f64 / rings as f64 - PI / 2.0;
|
||||
for i in 0..n {
|
||||
let a = TAU * i as f64 / n as f64;
|
||||
p.push([r * th.cos() * a.cos(), r * th.cos() * a.sin(), r * th.sin()]);
|
||||
}
|
||||
}
|
||||
let id = |j: usize, i: usize| (2 + (j - 1) * n + i % n) as u32;
|
||||
let mut tris = Vec::new();
|
||||
for i in 0..n {
|
||||
tris.push([0, id(1, i + 1), id(1, i)]);
|
||||
tris.push([1, id(rings - 1, i), id(rings - 1, i + 1)]);
|
||||
}
|
||||
for j in 1..rings - 1 {
|
||||
for i in 0..n {
|
||||
tris.push([id(j, i), id(j, i + 1), id(j + 1, i + 1)]);
|
||||
tris.push([id(j, i), id(j + 1, i + 1), id(j + 1, i)]);
|
||||
}
|
||||
}
|
||||
let face = vec![0; tris.len()];
|
||||
let s = BoundarySurface::new(&p, &tris, &face, &SurfaceOptions::default()).unwrap();
|
||||
let t = Tessellation::new(&s, vec![Some(AnalyticSurface::Sphere { center: [0.0; 3], radius: r })]);
|
||||
(s, t)
|
||||
}
|
||||
|
||||
/// A closed cylinder of radius `r` from z = 0 to `h`: faces 0 (bottom), 1 (top), 2 (side).
|
||||
fn cylinder(r: f64, h: f64, n: usize) -> (BoundarySurface, Tessellation) {
|
||||
let mut p = vec![[0.0, 0.0, 0.0], [0.0, 0.0, h]];
|
||||
for z in [0.0, h] {
|
||||
for i in 0..n {
|
||||
let a = TAU * i as f64 / n as f64;
|
||||
p.push([r * a.cos(), r * a.sin(), z]);
|
||||
}
|
||||
}
|
||||
let (b, t) = (|i: usize| (2 + i % n) as u32, |i: usize| (2 + n + i % n) as u32);
|
||||
let mut tris = Vec::new();
|
||||
let mut face = Vec::new();
|
||||
for i in 0..n {
|
||||
tris.push([0, b(i + 1), b(i)]);
|
||||
face.push(0);
|
||||
tris.push([1, t(i), t(i + 1)]);
|
||||
face.push(1);
|
||||
tris.push([b(i), b(i + 1), t(i + 1)]);
|
||||
tris.push([b(i), t(i + 1), t(i)]);
|
||||
face.extend([2, 2]);
|
||||
}
|
||||
let s = BoundarySurface::new(&p, &tris, &face, &SurfaceOptions::default()).unwrap();
|
||||
let analytic = vec![
|
||||
Some(AnalyticSurface::Plane { normal: [0.0, 0.0, -1.0], point: [0.0; 3] }),
|
||||
Some(AnalyticSurface::Plane { normal: [0.0, 0.0, 1.0], point: [0.0, 0.0, h] }),
|
||||
Some(AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: r }),
|
||||
];
|
||||
let t = Tessellation::new(&s, analytic);
|
||||
(s, t)
|
||||
}
|
||||
|
||||
fn options(size: f64) -> ProjectOptions {
|
||||
ProjectOptions { size, ..ProjectOptions::default() }
|
||||
}
|
||||
|
||||
/// The 7-block sphere: all hex, every boundary node on the sphere, min scaled Jacobian at least
|
||||
/// 0.5 at every refinement, and the volume converging to the ball's.
|
||||
#[test]
|
||||
fn seven_block_sphere_converges_without_losing_quality() {
|
||||
let r = 10.0;
|
||||
let (s, t) = sphere(r, 96);
|
||||
let topo = Topology::sphere([0.0; 3], r, 0.55).unwrap();
|
||||
assert_eq!(topo.blocks.len(), 7);
|
||||
let exact = 4.0 / 3.0 * PI * r * r * r;
|
||||
let mut last_error = f64::INFINITY;
|
||||
for size in [4.0, 2.0, 1.0] {
|
||||
let m = project_mesh(&topo, &s, &t, &options(size)).unwrap();
|
||||
let rep = &m.report;
|
||||
assert!(rep.min_scaled_jacobian >= 0.5, "size {size}: {rep:?}");
|
||||
assert!(rep.max_gap < 1e-9 * r, "size {size}: {rep:?}");
|
||||
for b in &m.boundary {
|
||||
let hex = m.hexes[b.element as usize];
|
||||
for &l in ElementKind::Hex8.faces()[b.face as usize - 1] {
|
||||
let p = m.nodes[hex[l] as usize];
|
||||
assert!((crate::vec3::norm(p) - r).abs() < 1e-9, "{p:?}");
|
||||
}
|
||||
}
|
||||
let error = (rep.volume - exact).abs() / exact;
|
||||
assert!(error < last_error, "size {size}: {error} after {last_error}");
|
||||
last_error = error;
|
||||
}
|
||||
assert!(last_error < 0.01, "{last_error}");
|
||||
}
|
||||
|
||||
/// A disc: the five-block butterfly, its rim on the circle where the side meets the caps.
|
||||
#[test]
|
||||
fn butterfly_disc_meshes_with_its_rim_on_the_circle() {
|
||||
let (r, h) = (10.0, 2.0);
|
||||
let (s, t) = cylinder(r, h, 128);
|
||||
let topo = Topology::butterfly([0.0; 3], [0.0, 0.0, 1.0], r, &[0.0, h], 0.55).unwrap();
|
||||
assert_eq!(topo.blocks.len(), 5);
|
||||
let mut last = f64::INFINITY;
|
||||
for size in [2.0, 1.0, 0.5] {
|
||||
let m = project_mesh(&topo, &s, &t, &options(size)).unwrap();
|
||||
let rep = &m.report;
|
||||
assert!(rep.min_scaled_jacobian >= 0.5, "size {size}: {rep:?}");
|
||||
// Every node on the side lies on the cylinder; the rims on both circles.
|
||||
let side: Vec<&BoundaryFace> = m.boundary.iter().filter(|b| b.brep_face == 2).collect();
|
||||
assert!(!side.is_empty());
|
||||
for b in side {
|
||||
let hex = m.hexes[b.element as usize];
|
||||
for &l in ElementKind::Hex8.faces()[b.face as usize - 1] {
|
||||
let p = m.nodes[hex[l] as usize];
|
||||
assert!((p[0].hypot(p[1]) - r).abs() < 1e-9, "{p:?}");
|
||||
}
|
||||
}
|
||||
let tags: std::collections::BTreeSet<u32> = m.boundary.iter().map(|b| b.brep_face).collect();
|
||||
assert_eq!(tags.len(), 3);
|
||||
let error = (rep.volume - PI * r * r * h).abs();
|
||||
assert!(error < last, "size {size}");
|
||||
last = error;
|
||||
}
|
||||
}
|
||||
|
||||
/// A plain grid on a box: corners land on the box's corners, edges on its edges.
|
||||
#[test]
|
||||
fn a_grid_on_a_cylinder_is_refused_where_corners_would_flatten() {
|
||||
// A single block projected onto a cylinder puts its four side edges on the circle: the
|
||||
// block's corners flatten to 180°, which the mesher reports instead of meshing.
|
||||
let (s, t) = cylinder(5.0, 4.0, 64);
|
||||
let topo = Topology::grid([&[-3.0, 3.0], &[-3.0, 3.0], &[0.0, 4.0]]).unwrap();
|
||||
let e = project_mesh(&topo, &s, &t, &options(1.0)).unwrap_err().to_string();
|
||||
assert!(e.contains("O-grid") || e.contains("inverted"), "{e}");
|
||||
}
|
||||
|
||||
/// A geometry of exact surfaces, for the curve and corner iterations.
|
||||
struct Exact(Vec<AnalyticSurface>);
|
||||
|
||||
impl Geometry for Exact {
|
||||
fn project(&self, face: u32, p: P) -> Option<P> {
|
||||
self.0.get(face as usize)?.project(p)
|
||||
}
|
||||
fn normal(&self, face: u32, p: P) -> Option<P> {
|
||||
self.0.get(face as usize)?.normal_at(p)
|
||||
}
|
||||
}
|
||||
|
||||
/// Nodes on two targets converge to their intersection curve; on three, to the common point.
|
||||
#[test]
|
||||
fn curve_and_corner_nodes_converge_by_tangent_planes() {
|
||||
let g = Exact(vec![
|
||||
AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [1.0, 0.0, 0.0], radius: 10.0 },
|
||||
AnalyticSurface::Cylinder { origin: [0.0; 3], axis: [0.0, 0.0, 1.0], radius: 6.0 },
|
||||
AnalyticSurface::Plane { normal: [0.0, 1.0, 0.0], point: [0.0; 3] },
|
||||
]);
|
||||
let targets = vec![vec![0], vec![1], vec![2]];
|
||||
for start in [[3.0, 3.0, 9.0], [-5.0, 1.0, 7.0], [1.0, -4.0, 11.0]] {
|
||||
let p = place(&g, &targets, Constraint::Curve(0, 1), start, 1e-13);
|
||||
assert!((p[1].hypot(p[2]) - 10.0).abs() < 1e-9 && (p[0].hypot(p[1]) - 6.0).abs() < 1e-9, "{p:?}");
|
||||
}
|
||||
let p = place(&g, &targets, Constraint::Corner(0, 1, 2), [5.0, 1.0, 9.0], 1e-13);
|
||||
assert!((p[0] - 6.0).abs() < 1e-9 && p[1].abs() < 1e-9 && (p[2] - 10.0).abs() < 1e-9, "{p:?}");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn o_grid_keeps_neighbours_and_opens_faces() {
|
||||
let mut t = Topology::grid([&[0.0, 1.0, 2.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap();
|
||||
let on = t.o_grid(0, &[2], [0.5; 3]).unwrap();
|
||||
assert_eq!(on[0].len(), 5);
|
||||
// Inner block + five side blocks (top open), the neighbour untouched.
|
||||
assert_eq!(t.blocks.len(), 7);
|
||||
t.check().unwrap();
|
||||
let uses = lattice::face_uses(&t).unwrap();
|
||||
// The top of the old block is now five faces, each used once.
|
||||
// Inner block to its five side blocks, the four ring blocks to the bottom one and to each
|
||||
// other, and the untouched neighbour.
|
||||
assert_eq!(uses.values().filter(|&&u| u == 2).count(), 5 + 4 + 4 + 1, "{uses:?}");
|
||||
let made = t.extrude(&on[0], [0.0, 0.0, 1.0], 2).unwrap();
|
||||
assert_eq!(made.len(), 10);
|
||||
t.check().unwrap();
|
||||
let face = t.face_towards(made[8], [0.0, 0.0, 1.0]).unwrap();
|
||||
assert_eq!(face, 2);
|
||||
let up = t.extrude(&[(made[8], 2)], [0.0, 0.0, 1.0], 1).unwrap();
|
||||
let side = t.face_towards(1, [1.0, 0.0, 0.0]).unwrap();
|
||||
t.connect((up[0], 2), (1, side), 2).unwrap();
|
||||
t.remove(&[1]).unwrap();
|
||||
t.check().unwrap();
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn hostile_topologies_are_errors() {
|
||||
let (s, g) = sphere(1.0, 32);
|
||||
let cube = Topology::grid([&[0.0, 1.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap();
|
||||
let mut inside_out = cube.clone();
|
||||
inside_out.blocks[0].swap(1, 3);
|
||||
let mut repeated = cube.clone();
|
||||
repeated.blocks[0][1] = repeated.blocks[0][0];
|
||||
let mut missing = cube.clone();
|
||||
missing.blocks[0][7] = 99;
|
||||
let mut nan = cube.clone();
|
||||
nan.vertices[0] = [f64::NAN, 0.0, 0.0];
|
||||
for t in [Topology::default(), inside_out, repeated, missing, nan] {
|
||||
assert!(project_mesh(&t, &s, &g, &options(0.5)).is_err(), "{t:?}");
|
||||
}
|
||||
for size in [0.0, -1.0, f64::NAN, f64::INFINITY] {
|
||||
assert!(project_mesh(&cube, &s, &g, &options(size)).is_err());
|
||||
}
|
||||
assert!(project_mesh(&cube, &s, &g, &ProjectOptions { size: 1e-6, max_elements: 1000, ..ProjectOptions::default() }).is_err());
|
||||
let bad_face = ProjectOptions { faces: vec![(0, 9, vec![0])], ..options(0.5) };
|
||||
assert!(project_mesh(&cube, &s, &g, &bad_face).is_err());
|
||||
let bad_target = ProjectOptions { faces: vec![(0, 1, vec![7])], ..options(0.5) };
|
||||
assert!(project_mesh(&cube, &s, &g, &bad_target).is_err());
|
||||
let mut t = cube.clone();
|
||||
assert!(t.o_grid(0, &[], [1.5; 3]).is_err() && t.o_grid(0, &[7], [0.5; 3]).is_err() && t.o_grid(3, &[], [0.5; 3]).is_err());
|
||||
assert!(t.extrude(&[(0, 0)], [0.0, 0.0, 1.0], 1).is_err() && t.extrude(&[(0, 1)], [f64::NAN; 3], 1).is_err());
|
||||
assert!(Topology::grid([&[0.0], &[0.0, 1.0], &[0.0, 1.0]]).is_err());
|
||||
assert!(Topology::grid([&[1.0, 0.0], &[0.0, 1.0], &[0.0, 1.0]]).is_err());
|
||||
assert!(Topology::sphere([0.0; 3], -1.0, 0.5).is_err() && Topology::butterfly([0.0; 3], [0.0; 3], 1.0, &[0.0, 1.0], 0.5).is_err());
|
||||
let mut g3 = Topology::grid([&[0.0, 1.0], &[0.0, 1.0], &[0.0, 1.0]]).unwrap();
|
||||
assert!(g3.carve([1, 1, 1], &[[[0, 0], [0, 0], [0, 0]]]).is_err());
|
||||
}
|
||||
|
||||
@@ -0,0 +1,514 @@
|
||||
//! The block topology of the projection method: hexahedral blocks over shared vertices, built
|
||||
//! from a carved grid, from the round-feature templates (butterfly, seven-block sphere) and by the
|
||||
//! O-grid and extrude operations.
|
||||
//!
|
||||
//! A block's corners are in hex8 order (1-2-3-4 counter-clockwise seen from 5-6-7-8's side) and
|
||||
//! its parameter directions run 1→2, 1→4 and 1→5. Faces are numbered as hex8 faces: 1 = 1-2-3-4,
|
||||
//! 2 = 5-6-7-8, 3 = 1-2-6-5, 4 = 2-3-7-6, 5 = 3-4-8-7, 6 = 4-1-5-8. Positions are only a first
|
||||
//! guess: the projection moves every boundary vertex onto the geometry.
|
||||
|
||||
use std::collections::BTreeMap;
|
||||
|
||||
use crate::vec3::{P, add, det, dist, scale, sub};
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// Most blocks one topology holds.
|
||||
pub const MAX_BLOCKS: usize = 20_000;
|
||||
|
||||
/// Corner `c` of a block in parameter space (hex8 order).
|
||||
pub(crate) const CORNERS: [[usize; 3]; 8] = [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]];
|
||||
|
||||
/// The direction and side of each hex8 face (1-based face `f` at index `f - 1`).
|
||||
pub(crate) const FACE_SIDES: [(usize, usize); 6] = [(2, 0), (2, 1), (1, 0), (0, 1), (1, 1), (0, 0)];
|
||||
|
||||
/// The corner at parameter `c`.
|
||||
pub(crate) fn corner_at(c: [usize; 3]) -> usize {
|
||||
CORNERS.iter().position(|x| *x == c).unwrap_or(0)
|
||||
}
|
||||
|
||||
/// Hex8 face number (1-based) of direction `d`, side `s`.
|
||||
pub(crate) fn face_number(d: usize, s: usize) -> u8 {
|
||||
FACE_SIDES.iter().position(|x| *x == (d, s)).map_or(1, |i| i as u8 + 1)
|
||||
}
|
||||
|
||||
/// Blocks over shared vertices.
|
||||
#[derive(Debug, Clone, Default, PartialEq)]
|
||||
pub struct Topology {
|
||||
pub vertices: Vec<P>,
|
||||
pub blocks: Vec<[u32; 8]>,
|
||||
}
|
||||
|
||||
/// The trilinear map of a block's corners at parameter `u`.
|
||||
pub(crate) fn trilinear(c: &[P; 8], u: [f64; 3]) -> P {
|
||||
let mut x = [0.0; 3];
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
let w: f64 = (0..3).map(|d| if bits[d] == 1 { u[d] } else { 1.0 - u[d] }).product();
|
||||
x = add(x, scale(c[k], w));
|
||||
}
|
||||
x
|
||||
}
|
||||
|
||||
/// The smallest corner Jacobian determinant of a trilinear block (positive for a valid one).
|
||||
pub(crate) fn min_corner_det(c: &[P; 8]) -> f64 {
|
||||
let mut m = f64::INFINITY;
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
// The three edges leaving corner k along increasing parameters, signed to keep the
|
||||
// right-handed order.
|
||||
let along = |d: usize| {
|
||||
let mut b = *bits;
|
||||
b[d] = 1 - b[d];
|
||||
let e = sub(c[corner_at(b)], c[k]);
|
||||
if bits[d] == 1 { scale(e, -1.0) } else { e }
|
||||
};
|
||||
m = m.min(det(along(0), along(1), along(2)));
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
impl Topology {
|
||||
/// The corner positions of block `b`.
|
||||
pub fn corners(&self, b: usize) -> Result<[P; 8]> {
|
||||
let block = self.blocks.get(b).ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
|
||||
let mut out = [[0.0; 3]; 8];
|
||||
for (o, &v) in out.iter_mut().zip(block) {
|
||||
*o = *self.vertices.get(v as usize).ok_or_else(|| Error::Invalid(format!("block {b} uses vertex {v}, which does not exist")))?;
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Every block uses existing, distinct vertices, is right-handed and caps are kept.
|
||||
pub fn check(&self) -> Result<()> {
|
||||
if self.blocks.is_empty() {
|
||||
return Err(Error::Invalid("the topology has no blocks".into()));
|
||||
}
|
||||
if self.blocks.len() > MAX_BLOCKS {
|
||||
return Err(Error::Limit(format!("{} blocks; at most {MAX_BLOCKS}", self.blocks.len())));
|
||||
}
|
||||
if let Some(v) = self.vertices.iter().position(|v| !v.iter().all(|x| x.is_finite() && x.abs() <= fieldcraft_model::MAX_COORD)) {
|
||||
return Err(Error::Invalid(format!("vertex {v} is not a finite point")));
|
||||
}
|
||||
for (b, block) in self.blocks.iter().enumerate() {
|
||||
let mut sorted = *block;
|
||||
sorted.sort_unstable();
|
||||
if sorted.windows(2).any(|w| w[0] == w[1]) {
|
||||
return Err(Error::Invalid(format!("block {b} uses one vertex twice")));
|
||||
}
|
||||
if !(min_corner_det(&self.corners(b)?) > 0.0) {
|
||||
return Err(Error::Invalid(format!(
|
||||
"block {b} is inside out or folded at a corner: give its corners in hex8 order (1-2-3-4 counter-clockwise seen from 5-6-7-8)"
|
||||
)));
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// The bounding-box diagonal of the vertices.
|
||||
pub fn diagonal(&self) -> f64 {
|
||||
let (mut lo, mut hi) = ([f64::INFINITY; 3], [f64::NEG_INFINITY; 3]);
|
||||
for v in &self.vertices {
|
||||
for d in 0..3 {
|
||||
lo[d] = lo[d].min(v[d]);
|
||||
hi[d] = hi[d].max(v[d]);
|
||||
}
|
||||
}
|
||||
if self.vertices.is_empty() { 0.0 } else { dist(lo, hi) }
|
||||
}
|
||||
|
||||
/// A grid of blocks between the planes `planes[d]` (increasing coordinates along X, Y and Z),
|
||||
/// block `(i, j, k)` at index `(k·nj + j)·ni + i`.
|
||||
pub fn grid(planes: [&[f64]; 3]) -> Result<Topology> {
|
||||
for (d, p) in planes.iter().enumerate() {
|
||||
if p.len() < 2 || !p.iter().all(|x| x.is_finite()) || p.windows(2).any(|w| !(w[1] > w[0])) {
|
||||
return Err(Error::Invalid(format!("the grid planes along {} must be at least two increasing numbers", ["X", "Y", "Z"][d])));
|
||||
}
|
||||
}
|
||||
let n = planes.map(|p| p.len() - 1);
|
||||
if n.iter().try_fold(1usize, |a, &b| a.checked_mul(b)).is_none_or(|c| c > MAX_BLOCKS) {
|
||||
return Err(Error::Limit(format!("a grid of {n:?} blocks; at most {MAX_BLOCKS}")));
|
||||
}
|
||||
let vid = |i: usize, j: usize, k: usize| ((k * (n[1] + 1) + j) * (n[0] + 1) + i) as u32;
|
||||
let mut vertices = Vec::new();
|
||||
for z in planes[2] {
|
||||
for y in planes[1] {
|
||||
for x in planes[0] {
|
||||
vertices.push([*x, *y, *z]);
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut blocks = Vec::new();
|
||||
for k in 0..n[2] {
|
||||
for j in 0..n[1] {
|
||||
for i in 0..n[0] {
|
||||
blocks.push(CORNERS.map(|c| vid(i + c[0], j + c[1], k + c[2])));
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(Topology { vertices, blocks })
|
||||
}
|
||||
|
||||
/// Remove the blocks of a grid of `n` blocks per direction whose index lies in every range
|
||||
/// of `ranges` (inclusive, per direction), keeping the others in order; vertices no block
|
||||
/// uses are dropped.
|
||||
pub fn carve(&mut self, n: [usize; 3], ranges: &[[[usize; 2]; 3]]) -> Result<()> {
|
||||
if self.blocks.len() != n.iter().product::<usize>() {
|
||||
return Err(Error::Invalid("carving needs the grid's own block numbering".into()));
|
||||
}
|
||||
let mut keep = Vec::new();
|
||||
for (b, block) in self.blocks.iter().enumerate() {
|
||||
let at = [b % n[0], (b / n[0]) % n[1], b / (n[0] * n[1])];
|
||||
let gone = ranges.iter().any(|r| (0..3).all(|d| r[d][0] <= at[d] && at[d] <= r[d][1]));
|
||||
if !gone {
|
||||
keep.push(*block);
|
||||
}
|
||||
}
|
||||
self.blocks = keep;
|
||||
self.compact();
|
||||
if self.blocks.is_empty() { Err(Error::Invalid("carving removed every block".into())) } else { Ok(()) }
|
||||
}
|
||||
|
||||
/// Drop unused vertices and renumber.
|
||||
fn compact(&mut self) {
|
||||
let mut map: BTreeMap<u32, u32> = BTreeMap::new();
|
||||
let mut vertices = Vec::new();
|
||||
for block in &mut self.blocks {
|
||||
for v in block.iter_mut() {
|
||||
let next = map.len() as u32;
|
||||
let new = *map.entry(*v).or_insert_with(|| {
|
||||
vertices.push(self.vertices.get(*v as usize).copied().unwrap_or([0.0; 3]));
|
||||
next
|
||||
});
|
||||
*v = new;
|
||||
}
|
||||
}
|
||||
self.vertices = vertices;
|
||||
}
|
||||
|
||||
/// Join vertices closer than `tol` (the first one stays).
|
||||
pub fn merge(&mut self, tol: f64) {
|
||||
let mut map: Vec<u32> = (0..self.vertices.len() as u32).collect();
|
||||
// Sorted by X, so only a window is compared.
|
||||
let mut order: Vec<usize> = (0..self.vertices.len()).collect();
|
||||
order.sort_by(|&a, &b| self.vertices[a][0].total_cmp(&self.vertices[b][0]).then(a.cmp(&b)));
|
||||
for (oi, &a) in order.iter().enumerate() {
|
||||
if map[a] as usize != a {
|
||||
continue;
|
||||
}
|
||||
for &b in order.iter().skip(oi + 1) {
|
||||
if self.vertices[b][0] - self.vertices[a][0] > tol {
|
||||
break;
|
||||
}
|
||||
if map[b] as usize == b && dist(self.vertices[a], self.vertices[b]) <= tol {
|
||||
map[b] = a as u32;
|
||||
}
|
||||
}
|
||||
}
|
||||
for block in &mut self.blocks {
|
||||
for v in block.iter_mut() {
|
||||
*v = map.get(*v as usize).copied().unwrap_or(*v);
|
||||
}
|
||||
}
|
||||
self.compact();
|
||||
}
|
||||
|
||||
/// Replace block `b` by an O-grid: an inner block, scaled by `inset[d]` (0 to 1) about the
|
||||
/// block's centre along each of its directions (from the open face, along a direction with one
|
||||
/// open face), and one block between each face of the inner
|
||||
/// block and the block's own face, except across the `open` faces (hex8 face numbers), where
|
||||
/// the inner block reaches the face. Faces of `b` stay faces of the result, so neighbours
|
||||
/// still fit, except the open ones, which become the inner block's face ringed by the side
|
||||
/// blocks' faces (the butterfly pattern on that face). The inner block takes index `b`; the
|
||||
/// others are appended. Returns, per open face, the faces now on it: the inner block's first,
|
||||
/// then the ring's.
|
||||
pub fn o_grid(&mut self, b: usize, open: &[u8], inset: [f64; 3]) -> Result<Vec<Vec<(usize, u8)>>> {
|
||||
if !inset.iter().all(|i| *i > 0.0 && *i < 1.0) {
|
||||
return Err(Error::Invalid(format!("the O-grid inset must be between 0 and 1, got {inset:?}")));
|
||||
}
|
||||
if let Some(f) = open.iter().find(|f| !(1..=6).contains(*f)) {
|
||||
return Err(Error::Invalid(format!("face {f} is not a hex8 face (1 to 6)")));
|
||||
}
|
||||
let c = self.corners(b)?;
|
||||
let outer = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
|
||||
let is_open = |d: usize, s: usize| open.contains(&face_number(d, s));
|
||||
// Each corner's inner copy: on an open face it stays on the face.
|
||||
let mut inner = [0u32; 8];
|
||||
for (k, bits) in CORNERS.iter().enumerate() {
|
||||
// Along a direction with one open face the inner block reaches `inset` of the way
|
||||
// across from it; with none it is centred.
|
||||
let u = [0, 1, 2].map(|d| {
|
||||
let b = bits[d] as f64;
|
||||
match (is_open(d, 0), is_open(d, 1)) {
|
||||
(true, true) => b,
|
||||
(true, false) => {
|
||||
if bits[d] == 0 {
|
||||
0.0
|
||||
} else {
|
||||
inset[d]
|
||||
}
|
||||
}
|
||||
(false, true) => {
|
||||
if bits[d] == 1 {
|
||||
1.0
|
||||
} else {
|
||||
1.0 - inset[d]
|
||||
}
|
||||
}
|
||||
(false, false) => 0.5 + inset[d] * (b - 0.5),
|
||||
}
|
||||
});
|
||||
inner[k] = self.vertices.len() as u32;
|
||||
self.vertices.push(trilinear(&c, u));
|
||||
}
|
||||
if let Some(slot) = self.blocks.get_mut(b) {
|
||||
*slot = inner;
|
||||
}
|
||||
let mut on_open: Vec<Vec<(usize, u8)>> = open.iter().map(|&f| vec![(b, f)]).collect();
|
||||
for (d, s) in FACE_SIDES {
|
||||
if is_open(d, s) {
|
||||
continue;
|
||||
}
|
||||
// Local directions: the face's two directions, then from the outer face inwards.
|
||||
let (e1, e2) = match d {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
let pick = |a: usize, bb: usize, w: usize| {
|
||||
let mut bits = [0; 3];
|
||||
bits[d] = s;
|
||||
bits[e1] = a;
|
||||
bits[e2] = bb;
|
||||
let k = corner_at(bits);
|
||||
if w == 0 { outer[k] } else { inner[k] }
|
||||
};
|
||||
let mut block = CORNERS.map(|q| pick(q[0], q[1], q[2]));
|
||||
let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
let flipped = min_corner_det(&pos(&block)) < 0.0;
|
||||
if flipped {
|
||||
block = CORNERS.map(|q| pick(q[1], q[0], q[2]));
|
||||
}
|
||||
let index = self.blocks.len();
|
||||
self.blocks.push(block);
|
||||
// This block's face on each open face of the old block.
|
||||
for (list, &f) in on_open.iter_mut().zip(open) {
|
||||
let Some(&(od, os)) = FACE_SIDES.get(usize::from(f) - 1) else { continue };
|
||||
let local = if od == e1 {
|
||||
0
|
||||
} else if od == e2 {
|
||||
1
|
||||
} else {
|
||||
continue;
|
||||
};
|
||||
let local = if flipped { 1 - local } else { local };
|
||||
list.push((index, face_number(local, os)));
|
||||
}
|
||||
}
|
||||
Ok(on_open)
|
||||
}
|
||||
|
||||
/// The hex8 face of block `b` whose outward direction is nearest `dir`.
|
||||
pub fn face_towards(&self, b: usize, dir: P) -> Result<u8> {
|
||||
let c = self.corners(b)?;
|
||||
let centre = scale(c.iter().fold([0.0; 3], |a, p| add(a, *p)), 0.125);
|
||||
let mut best = (f64::NEG_INFINITY, 1u8);
|
||||
for (fi, (d, s)) in FACE_SIDES.iter().enumerate() {
|
||||
let fc = CORNERS.iter().zip(&c).filter(|(bits, _)| bits[*d] == *s).fold([0.0; 3], |a, (_, p)| add(a, *p));
|
||||
let out = sub(scale(fc, 0.25), centre);
|
||||
let l = crate::vec3::norm(out);
|
||||
let score = if l > 0.0 { crate::vec3::dot(out, dir) / l } else { f64::NEG_INFINITY };
|
||||
if score > best.0 {
|
||||
best = (score, fi as u8 + 1);
|
||||
}
|
||||
}
|
||||
Ok(best.1)
|
||||
}
|
||||
|
||||
/// The corners of a block face in cyclic order.
|
||||
fn face_quad(&self, b: usize, f: u8) -> Result<[u32; 4]> {
|
||||
let (d, s) = FACE_SIDES.get(usize::from(f).wrapping_sub(1)).copied().ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face")))?;
|
||||
let block = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
|
||||
let (e1, e2) = match d {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
Ok([[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, bb]| {
|
||||
let mut bits = [0; 3];
|
||||
bits[d] = s;
|
||||
bits[e1] = a;
|
||||
bits[e2] = bb;
|
||||
block[corner_at(bits)]
|
||||
}))
|
||||
}
|
||||
|
||||
/// Fill the space between two block faces with `layers` blocks, each face's corner matched
|
||||
/// to the nearest of the other's (over the eight ways of laying one quad on the other).
|
||||
/// Returns the new blocks.
|
||||
pub fn connect(&mut self, from: (usize, u8), to: (usize, u8), layers: usize) -> Result<Vec<usize>> {
|
||||
if !(1..=1000).contains(&layers) {
|
||||
return Err(Error::Invalid(format!("a connection has 1 to 1000 layers, got {layers}")));
|
||||
}
|
||||
let a = self.face_quad(from.0, from.1)?;
|
||||
let b = self.face_quad(to.0, to.1)?;
|
||||
let verts = self.vertices.clone();
|
||||
let at = |v: u32| verts.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
let mut best: Option<(f64, [u32; 4])> = None;
|
||||
for rev in [false, true] {
|
||||
for shift in 0..4 {
|
||||
let m: [u32; 4] = [0, 1, 2, 3].map(|i| {
|
||||
let k = if rev { (4 + shift - i) % 4 } else { (shift + i) % 4 };
|
||||
b[k]
|
||||
});
|
||||
let cost: f64 = (0..4).map(|i| dist(at(a[i]), at(m[i]))).sum();
|
||||
if best.is_none_or(|x| cost < x.0) {
|
||||
best = Some((cost, m));
|
||||
}
|
||||
}
|
||||
}
|
||||
let (_, m) = best.ok_or_else(|| Error::Invalid("no way to match the faces".into()))?;
|
||||
// Intermediate layers of vertices, then the blocks.
|
||||
let mut rows: Vec<[u32; 4]> = vec![a];
|
||||
for l in 1..layers {
|
||||
let t = l as f64 / layers as f64;
|
||||
let pts = [0, 1, 2, 3].map(|i| crate::vec3::lerp(at(a[i]), at(m[i]), t));
|
||||
let row = pts.map(|p| {
|
||||
self.vertices.push(p);
|
||||
self.vertices.len() as u32 - 1
|
||||
});
|
||||
rows.push(row);
|
||||
}
|
||||
rows.push(m);
|
||||
let mut made = Vec::new();
|
||||
for w in rows.windows(2) {
|
||||
let (q, t) = (w[0], w[1]);
|
||||
let mut block = [q[0], q[1], q[2], q[3], t[0], t[1], t[2], t[3]];
|
||||
let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
if min_corner_det(&pos(&block)) < 0.0 {
|
||||
block = [q[0], q[3], q[2], q[1], t[0], t[3], t[2], t[1]];
|
||||
}
|
||||
made.push(self.blocks.len());
|
||||
self.blocks.push(block);
|
||||
}
|
||||
Ok(made)
|
||||
}
|
||||
|
||||
/// Remove blocks (the others keep their order and are renumbered).
|
||||
pub fn remove(&mut self, blocks: &[usize]) -> Result<()> {
|
||||
if let Some(b) = blocks.iter().find(|b| **b >= self.blocks.len()) {
|
||||
return Err(Error::Invalid(format!("there is no block {b}")));
|
||||
}
|
||||
let mut k = 0;
|
||||
self.blocks.retain(|_| {
|
||||
k += 1;
|
||||
!blocks.contains(&(k - 1))
|
||||
});
|
||||
self.compact();
|
||||
if self.blocks.is_empty() { Err(Error::Invalid("every block was removed".into())) } else { Ok(()) }
|
||||
}
|
||||
|
||||
/// Extrude block faces `(block, hex8 face)` by `offset` in `layers` blocks each; the faces'
|
||||
/// shared vertices stay shared. Returns the new blocks' indices, the last layer's outer face
|
||||
/// is face 2 of the blocks made from it.
|
||||
pub fn extrude(&mut self, faces: &[(usize, u8)], offset: P, layers: usize) -> Result<Vec<usize>> {
|
||||
if !(1..=1000).contains(&layers) {
|
||||
return Err(Error::Invalid(format!("an extrusion has 1 to 1000 layers, got {layers}")));
|
||||
}
|
||||
if !offset.iter().all(|x| x.is_finite()) {
|
||||
return Err(Error::Invalid("the extrusion offset must be finite".into()));
|
||||
}
|
||||
let mut quads = Vec::new();
|
||||
for &(b, f) in faces {
|
||||
let (d, s) =
|
||||
FACE_SIDES.get(usize::from(f).wrapping_sub(1)).copied().ok_or_else(|| Error::Invalid(format!("face {f} is not a hex8 face")))?;
|
||||
let block = self.blocks.get(b).copied().ok_or_else(|| Error::Invalid(format!("there is no block {b}")))?;
|
||||
let (e1, e2) = match d {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
let q = [[0, 0], [1, 0], [1, 1], [0, 1]].map(|[a, bb]| {
|
||||
let mut bits = [0; 3];
|
||||
bits[d] = s;
|
||||
bits[e1] = a;
|
||||
bits[e2] = bb;
|
||||
block[corner_at(bits)]
|
||||
});
|
||||
quads.push(q);
|
||||
}
|
||||
let mut made = Vec::new();
|
||||
let mut layer_of: BTreeMap<u32, u32> = BTreeMap::new();
|
||||
let mut base: Vec<[u32; 4]> = quads;
|
||||
for _ in 0..layers {
|
||||
layer_of.clear();
|
||||
let mut next = Vec::new();
|
||||
for q in &base {
|
||||
let mut top = [0u32; 4];
|
||||
for (t, &v) in top.iter_mut().zip(q) {
|
||||
*t = match layer_of.get(&v) {
|
||||
Some(&w) => w,
|
||||
None => {
|
||||
let p = self.vertices.get(v as usize).copied().ok_or_else(|| Error::Invalid(format!("vertex {v} does not exist")))?;
|
||||
let w = self.vertices.len() as u32;
|
||||
self.vertices.push(add(p, scale(offset, 1.0 / layers as f64)));
|
||||
layer_of.insert(v, w);
|
||||
w
|
||||
}
|
||||
};
|
||||
}
|
||||
let mut block = [q[0], q[1], q[2], q[3], top[0], top[1], top[2], top[3]];
|
||||
let pos = |blk: &[u32; 8]| blk.map(|v| self.vertices.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
if min_corner_det(&pos(&block)) < 0.0 {
|
||||
block = [q[0], q[3], q[2], q[1], top[0], top[3], top[2], top[1]];
|
||||
}
|
||||
made.push(self.blocks.len());
|
||||
self.blocks.push(block);
|
||||
next.push(top);
|
||||
}
|
||||
base = next;
|
||||
}
|
||||
Ok(made)
|
||||
}
|
||||
|
||||
/// The seven-block sphere: an inner cube with its corners at `inset · radius` from the
|
||||
/// centre and six blocks between the cube's faces and the sphere (the cube-to-sphere O-grid).
|
||||
pub fn sphere(center: P, radius: f64, inset: f64) -> Result<Topology> {
|
||||
if !(radius > 0.0 && radius.is_finite()) || !center.iter().all(|x| x.is_finite()) {
|
||||
return Err(Error::Invalid(format!("a sphere needs a finite centre and a positive radius, got {radius}")));
|
||||
}
|
||||
let h = radius / 3f64.sqrt();
|
||||
let mut t = Topology::grid([&[center[0] - h, center[0] + h], &[center[1] - h, center[1] + h], &[center[2] - h, center[2] + h]])?;
|
||||
t.o_grid(0, &[], [inset; 3])?;
|
||||
Ok(t)
|
||||
}
|
||||
|
||||
/// The butterfly (five-block O-grid of a disc) of radius `radius` about `center` in the plane
|
||||
/// normal to `axis`, stacked through `levels` (distances along `axis` from `center`,
|
||||
/// increasing): a centre square with its corners at `inset · radius` and four blocks around
|
||||
/// it out to the circle, per level.
|
||||
pub fn butterfly(center: P, axis: P, radius: f64, levels: &[f64], inset: f64) -> Result<Topology> {
|
||||
if !(radius > 0.0 && radius.is_finite()) || !center.iter().all(|x| x.is_finite()) {
|
||||
return Err(Error::Invalid(format!("a butterfly needs a finite centre and a positive radius, got {radius}")));
|
||||
}
|
||||
let z = crate::vec3::norm(axis);
|
||||
if !(z > 0.0 && z.is_finite()) {
|
||||
return Err(Error::Invalid("the butterfly's axis must be a nonzero direction".into()));
|
||||
}
|
||||
let w = scale(axis, 1.0 / z);
|
||||
let t0 = if w[0].abs() < 0.9 { [1.0, 0.0, 0.0] } else { [0.0, 1.0, 0.0] };
|
||||
let u = {
|
||||
let c = crate::vec3::cross(w, t0);
|
||||
scale(c, 1.0 / crate::vec3::norm(c))
|
||||
};
|
||||
let v = crate::vec3::cross(w, u);
|
||||
let h = radius / 2f64.sqrt();
|
||||
let grid = Topology::grid([&[-h, h], &[-h, h], levels])?;
|
||||
let mut t = Topology {
|
||||
vertices: grid.vertices.iter().map(|p| add(center, add(add(scale(u, p[0]), scale(v, p[1])), scale(w, p[2])))).collect(),
|
||||
blocks: grid.blocks,
|
||||
};
|
||||
for b in 0..levels.len() - 1 {
|
||||
t.o_grid(b, &[1, 2], [inset; 3])?;
|
||||
}
|
||||
t.merge(1e-9 * (radius + levels.iter().fold(0.0f64, |m, x| m.max(x.abs()))));
|
||||
Ok(t)
|
||||
}
|
||||
}
|
||||
@@ -496,6 +496,53 @@ impl Tessellation {
|
||||
}
|
||||
}
|
||||
|
||||
impl Tessellation {
|
||||
/// The closest point to `p` on the trimmed face `face`: the exact surface's closest point
|
||||
/// when it lies on the face (within a fifth of the local triangle size of the face's
|
||||
/// triangles), else the closest point of the face's triangles. Unlike
|
||||
/// [`BoundaryGeometry::project_face`] it never returns a point of the surface's extension
|
||||
/// beyond the face's boundary, whatever the distance from `p`.
|
||||
pub fn project_trimmed(&self, face: u32, p: P) -> Option<P> {
|
||||
let accept = |t: u32| self.tri_face.get(t as usize) == Some(&face);
|
||||
let nearest = |x: P| {
|
||||
let dist = |t: u32| {
|
||||
let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]);
|
||||
let (a, b, c) = (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2]));
|
||||
let q = closest_on_triangle(x, a, b, c);
|
||||
let d = sub(q, x);
|
||||
let mut e = [measure::dist(a, b), measure::dist(b, c), measure::dist(c, a)];
|
||||
e.sort_by(f64::total_cmp);
|
||||
(dot(d, d), [q[0], q[1], q[2], e[1]])
|
||||
};
|
||||
self.tri_tree.nearest4(x, &accept, &dist).map(|r| ([r[0], r[1], r[2]], r[3]))
|
||||
};
|
||||
let (near, _) = nearest(p)?;
|
||||
let Some(q) = self.analytic(face).and_then(|s| s.project(p)) else { return Some(near) };
|
||||
let (back, size) = nearest(q)?;
|
||||
Some(if measure::dist(back, q) <= 0.2 * size { q } else { near })
|
||||
}
|
||||
|
||||
/// The outward unit normal of B-rep face `face` near `p`: the exact surface's (oriented like
|
||||
/// the nearest triangle) when known, else the nearest triangle's.
|
||||
pub fn face_normal(&self, face: u32, p: P) -> Option<P> {
|
||||
let accept = |t: u32| self.tri_face.get(t as usize) == Some(&face);
|
||||
let dist = |t: u32| {
|
||||
let tri = self.triangles.get(t as usize).copied().unwrap_or([0; 3]);
|
||||
let (a, b, c) = (self.pt(tri[0]), self.pt(tri[1]), self.pt(tri[2]));
|
||||
let d = sub(closest_on_triangle(p, a, b, c), p);
|
||||
let n = unit(cross(sub(b, a), sub(c, a))).unwrap_or([0.0; 3]);
|
||||
(dot(d, d), [n[0], n[1], n[2], 0.0])
|
||||
};
|
||||
let r = self.tri_tree.nearest4(p, &accept, &dist)?;
|
||||
let flat = unit([r[0], r[1], r[2]])?;
|
||||
match self.analytic(face).and_then(|s| s.normal_at(p)) {
|
||||
Some(n) if dot(n, flat) < 0.0 => Some(scale(n, -1.0)),
|
||||
Some(n) => Some(n),
|
||||
None => Some(flat),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl BoundaryGeometry for Tessellation {
|
||||
fn project_face(&self, face: u32, p: P) -> Option<P> {
|
||||
let (near, scale_len) = self.on_tessellation(face, p)?;
|
||||
|
||||
Reference in new issue
Block a user