forked from 1b-agent/fieldcraft-upstream
mesh: unstructured all-hex octree mesher
`fieldcraft_mesh::octree::octree_hex_mesh` meshes any closed part with hex8 only: - an adaptive 3-refinement grid (27-tree) with strong balance and closure, filled with conforming node-marking transition templates (Schneiders; checked over all 48 symmetries); - a core of grid cells kept clear of the boundary (the HybridOctree buffer rule), with non-manifold cleanup and a genus check against the part; - a feature graph from the B-rep (corners, curves, patches) and a partition-first map of the core's boundary onto it, then two boundary layers: a padding layer and a pillow per patch, so every face boundary and sharp edge is captured by element edges; - a Garanzha-type untangling and fitting energy (shape term |J|^3/chi, a volume term and a soft fitting term, deterministic across platforms), a compass search on the worst elements and a certified repair; failures refine the grid near them and rebuild, within bounded attempts and growth. The output passes the certified Jacobian check of every element and a gate that the boundary encloses the elements' volume; boundary quads carry their B-rep faces. Quality is reported with the same scaled Jacobian as mesh.quality. Literature notes for each step are in docs/meshing/literature.md (B1.6-B1.10). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01NmAttHhUZTokgPeEJia8sm
This commit is contained in:
jyvikoandClaude Opus 5.5 authored and Irakli Kourtis committed
1 parent
741a803785
commit
3587b0de74
12 files changed
+3676
-7
No files matched your search
@@ -37,7 +37,7 @@
|
||||
mod arrangement;
|
||||
mod brep;
|
||||
mod cells;
|
||||
mod geom;
|
||||
pub(crate) mod geom;
|
||||
mod interval;
|
||||
mod layered;
|
||||
mod mapped;
|
||||
@@ -583,7 +583,7 @@ fn finish(cells: &Cells, raw: Raw, method: HexMethod) -> Result<HexMesh> {
|
||||
|
||||
/// `∫ x n_x dA` over a quad4 or quad8 face (outward normal): its contribution to the enclosed
|
||||
/// volume. The 3 × 3 Gauss rule integrates it exactly (degree at most 5 per direction).
|
||||
fn face_flux(face: &[P]) -> f64 {
|
||||
pub(crate) fn face_flux(face: &[P]) -> f64 {
|
||||
const G: [(f64, f64); 3] = [(-0.774_596_669_241_483_4, 5.0 / 9.0), (0.0, 8.0 / 9.0), (0.774_596_669_241_483_4, 5.0 / 9.0)];
|
||||
let mut sum = 0.0;
|
||||
for (r, wr) in G {
|
||||
|
||||
@@ -12,6 +12,7 @@ pub mod cdt;
|
||||
pub mod delaunay;
|
||||
pub mod hex;
|
||||
pub mod join;
|
||||
pub mod octree;
|
||||
pub mod planar;
|
||||
pub mod quad;
|
||||
pub mod quality;
|
||||
|
||||
@@ -0,0 +1,636 @@
|
||||
//! The boundary as the hex mesher sees it: a triangle tree for nearest-point, ray and box
|
||||
//! queries, the exact inside test, and the feature graph (corners, curves, patches).
|
||||
//!
|
||||
//! # Feature graph
|
||||
//! Every feature edge of the [`BoundarySurface`] (between two B-rep faces, or a crease inside a
|
||||
//! face) whose dihedral angle reaches the sharp angle is kept; smooth face boundaries (a fillet
|
||||
//! meeting its faces tangentially) are not features, the faces on either side share a patch.
|
||||
//! Neither are blades and slits, edges whose normals turn by more than the blade angle (a 3°
|
||||
//! wedge, a sphere touching a plane): no valid hexahedra can follow them, so the mesh rounds
|
||||
//! them within the boundary tolerance. *Corners* are vertices where the number of feature edges is not two, or where a
|
||||
//! curve turns by more than [`CORNER_TURN`]; *curves* are the chains of feature edges between
|
||||
//! corners, or closed loops; *patches* are the connected sets of triangles between feature edges,
|
||||
//! so every patch lies in one B-rep face and is bounded by curves. This is the partition the hex
|
||||
//! mesh's boundary must reproduce (Gao, Shen, Panozzo 2019, §3.3).
|
||||
|
||||
use std::collections::BTreeMap;
|
||||
|
||||
pub(crate) use crate::hex::geom::closest_on_triangle;
|
||||
use crate::predicates::{Axis, Point3, Sign, orient2d_projected, orient3d, snap_to_domain};
|
||||
use crate::surface::BoundarySurface;
|
||||
use crate::vec3::{P, cross, dot, norm, sub};
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// A degree-2 feature vertex whose curve turns by more than this (degrees) is a corner.
|
||||
pub(crate) const CORNER_TURN: f64 = 30.0;
|
||||
/// Triangles per tree leaf.
|
||||
const LEAF: usize = 4;
|
||||
|
||||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||||
pub(crate) struct Aabb {
|
||||
pub lo: P,
|
||||
pub hi: P,
|
||||
}
|
||||
|
||||
impl Aabb {
|
||||
pub const EMPTY: Aabb = Aabb { lo: [f64::INFINITY; 3], hi: [f64::NEG_INFINITY; 3] };
|
||||
pub fn add(&mut self, p: P) {
|
||||
for ((lo, hi), x) in self.lo.iter_mut().zip(&mut self.hi).zip(p) {
|
||||
*lo = lo.min(x);
|
||||
*hi = hi.max(x);
|
||||
}
|
||||
}
|
||||
fn join(&mut self, o: &Aabb) {
|
||||
self.add(o.lo);
|
||||
self.add(o.hi);
|
||||
}
|
||||
pub fn overlaps(&self, o: &Aabb) -> bool {
|
||||
(0..3).all(|k| self.lo[k] <= o.hi[k] && o.lo[k] <= self.hi[k])
|
||||
}
|
||||
/// Squared distance from a point to the box (0 inside).
|
||||
fn dist2(&self, p: P) -> f64 {
|
||||
(0..3).map(|k| (self.lo[k] - p[k]).max(0.0).max(p[k] - self.hi[k])).map(|d| d * d).sum()
|
||||
}
|
||||
/// Whether the segment `a + t (b − a)`, `t ∈ [0, 1]`, may meet the box (slab test, widened by
|
||||
/// a relative margin so rounding never drops a hit).
|
||||
fn meets_segment(&self, a: P, b: P) -> bool {
|
||||
let (mut t0, mut t1) = (0.0f64, 1.0f64);
|
||||
for k in 0..3 {
|
||||
let pad = 1e-9 * (self.hi[k] - self.lo[k]).abs().max(1e-300) + 1e-12 * (a[k].abs() + b[k].abs());
|
||||
let (lo, hi) = (self.lo[k] - pad, self.hi[k] + pad);
|
||||
let d = b[k] - a[k];
|
||||
if d == 0.0 {
|
||||
if a[k] < lo || a[k] > hi {
|
||||
return false;
|
||||
}
|
||||
continue;
|
||||
}
|
||||
let (mut u, mut v) = ((lo - a[k]) / d, (hi - a[k]) / d);
|
||||
if u > v {
|
||||
std::mem::swap(&mut u, &mut v);
|
||||
}
|
||||
t0 = t0.max(u);
|
||||
t1 = t1.min(v);
|
||||
if t0 > t1 {
|
||||
return false;
|
||||
}
|
||||
}
|
||||
true
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Debug, Clone, Copy)]
|
||||
struct Node {
|
||||
bounds: Aabb,
|
||||
children: Option<usize>,
|
||||
start: usize,
|
||||
end: usize,
|
||||
}
|
||||
|
||||
/// A median-split tree over triangle boxes.
|
||||
#[derive(Debug, Clone)]
|
||||
struct Tree {
|
||||
nodes: Vec<Node>,
|
||||
order: Vec<u32>,
|
||||
}
|
||||
|
||||
impl Tree {
|
||||
fn build(boxes: &[Aabb]) -> Tree {
|
||||
let mut order: Vec<u32> = (0..boxes.len() as u32).collect();
|
||||
let empty = Node { bounds: Aabb::EMPTY, children: None, start: 0, end: 0 };
|
||||
let mut nodes = vec![empty];
|
||||
let mut work = vec![(0usize, 0usize, boxes.len())];
|
||||
while let Some((ni, start, end)) = work.pop() {
|
||||
let Some(slice) = order.get_mut(start..end) else { continue };
|
||||
let (mut bounds, mut centres) = (Aabb::EMPTY, Aabb::EMPTY);
|
||||
for &t in slice.iter() {
|
||||
if let Some(b) = boxes.get(t as usize) {
|
||||
bounds.join(b);
|
||||
centres.add([0, 1, 2].map(|k| b.lo[k] + b.hi[k]));
|
||||
}
|
||||
}
|
||||
let mut node = Node { bounds, children: None, start, end };
|
||||
if end - start > LEAF {
|
||||
let ext = [0, 1, 2].map(|k| centres.hi[k] - centres.lo[k]);
|
||||
let axis = if ext[0] >= ext[1] && ext[0] >= ext[2] {
|
||||
0
|
||||
} else if ext[1] >= ext[2] {
|
||||
1
|
||||
} else {
|
||||
2
|
||||
};
|
||||
let key = |t: &u32| boxes.get(*t as usize).map_or(0.0, |b| b.lo[axis] + b.hi[axis]);
|
||||
let mid = (end - start) / 2;
|
||||
slice.select_nth_unstable_by(mid, |a, b| key(a).total_cmp(&key(b)).then(a.cmp(b)));
|
||||
let l = nodes.len();
|
||||
nodes.push(empty);
|
||||
nodes.push(empty);
|
||||
node.children = Some(l);
|
||||
work.push((l, start, start + mid));
|
||||
work.push((l + 1, start + mid, end));
|
||||
}
|
||||
if let Some(n) = nodes.get_mut(ni) {
|
||||
*n = node;
|
||||
}
|
||||
}
|
||||
Tree { nodes, order }
|
||||
}
|
||||
|
||||
fn leaf(&self, n: &Node) -> &[u32] {
|
||||
self.order.get(n.start..n.end).unwrap_or(&[])
|
||||
}
|
||||
}
|
||||
|
||||
pub(crate) fn lerp(a: P, b: P, t: f64) -> P {
|
||||
[0, 1, 2].map(|k| a[k] + (b[k] - a[k]) * t)
|
||||
}
|
||||
|
||||
/// The closest point of segment `ab` to `p` and its parameter.
|
||||
pub(crate) fn closest_on_segment(p: P, a: P, b: P) -> (P, f64) {
|
||||
let ab = sub(b, a);
|
||||
let l2 = dot(ab, ab);
|
||||
let t = if l2 > 0.0 { (dot(sub(p, a), ab) / l2).clamp(0.0, 1.0) } else { 0.0 };
|
||||
(lerp(a, b, t), t)
|
||||
}
|
||||
|
||||
pub(crate) fn dist2(a: P, b: P) -> f64 {
|
||||
let d = sub(a, b);
|
||||
dot(d, d)
|
||||
}
|
||||
|
||||
/// A chain of feature edges.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub(crate) struct Curve {
|
||||
/// Surface vertices in order; a closed loop repeats its first vertex at the end.
|
||||
pub vertices: Vec<u32>,
|
||||
pub closed: bool,
|
||||
/// The patches on either side (sorted, deduplicated: one entry for a crease seam).
|
||||
pub patches: Vec<u32>,
|
||||
}
|
||||
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub(crate) struct Patch {
|
||||
pub triangles: Vec<u32>,
|
||||
/// The B-rep faces of its triangles (sorted, deduplicated; one for a well-formed boundary).
|
||||
pub faces: Vec<u32>,
|
||||
pub area: f64,
|
||||
/// Twice the area over the length of the bounding curves: the width of a strip, the
|
||||
/// diameter of a disc, `∞` for a closed surface without curves.
|
||||
pub width: f64,
|
||||
}
|
||||
|
||||
/// The boundary with its queries and feature graph.
|
||||
pub(crate) struct Geometry<'a> {
|
||||
pub surface: &'a BoundarySurface,
|
||||
tree: Tree,
|
||||
pub boxes: Vec<Aabb>,
|
||||
pub normals: Vec<P>,
|
||||
pub tri_patch: Vec<u32>,
|
||||
pub patches: Vec<Patch>,
|
||||
pub curves: Vec<Curve>,
|
||||
/// Corner vertices, ascending.
|
||||
pub corners: Vec<u32>,
|
||||
pub bounds: Aabb,
|
||||
/// Far enough to leave the boundary from any point of its box.
|
||||
far: f64,
|
||||
}
|
||||
|
||||
impl<'a> Geometry<'a> {
|
||||
pub fn new(surface: &'a BoundarySurface, sharp_angle: f64, blade_angle: f64) -> Result<Geometry<'a>> {
|
||||
let pts = surface.points();
|
||||
let tris = surface.triangles();
|
||||
let xyz = |v: u32| pts.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
let mut bounds = Aabb::EMPTY;
|
||||
let mut boxes = Vec::with_capacity(tris.len());
|
||||
let mut normals = Vec::with_capacity(tris.len());
|
||||
for t in tris {
|
||||
let mut b = Aabb::EMPTY;
|
||||
for &v in t {
|
||||
b.add(xyz(v));
|
||||
}
|
||||
bounds.join(&b);
|
||||
boxes.push(b);
|
||||
let n = cross(sub(xyz(t[1]), xyz(t[0])), sub(xyz(t[2]), xyz(t[0])));
|
||||
let l = norm(n);
|
||||
normals.push(if l > 0.0 { n.map(|c| c / l) } else { [0.0; 3] });
|
||||
}
|
||||
if tris.is_empty() || !(0..3).all(|k| bounds.hi[k] >= bounds.lo[k]) {
|
||||
return Err(Error::Invalid("the boundary has no triangles".into()));
|
||||
}
|
||||
let tree = Tree::build(&boxes);
|
||||
let diag = norm(sub(bounds.hi, bounds.lo));
|
||||
let far = 4.0 * diag.max(1e-9) + (0..3).map(|k| bounds.lo[k].abs().max(bounds.hi[k].abs())).fold(0.0, f64::max);
|
||||
|
||||
// Feature edges and corners.
|
||||
// A boundary between two B-rep faces is kept or dropped as a whole, by the median angle
|
||||
// along it, so tessellation noise never breaks a curve into pieces; creases inside a
|
||||
// face go edge by edge.
|
||||
let face_of = |t: u32| surface.tri_face().get(t as usize).copied().unwrap_or(0);
|
||||
let mut pair_angles: BTreeMap<(u32, u32), Vec<f64>> = BTreeMap::new();
|
||||
for e in surface.feature_edges() {
|
||||
let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1]));
|
||||
if a != b {
|
||||
pair_angles.entry((a.min(b), a.max(b))).or_default().push(e.angle);
|
||||
}
|
||||
}
|
||||
let pair_kept: BTreeMap<(u32, u32), bool> = pair_angles
|
||||
.into_iter()
|
||||
.map(|(k, mut v)| {
|
||||
v.sort_by(f64::total_cmp);
|
||||
let m = v.get(v.len() / 2).copied().unwrap_or(0.0);
|
||||
(k, m >= sharp_angle && m <= blade_angle)
|
||||
})
|
||||
.collect();
|
||||
let kept = |e: &&crate::surface::FeatureEdge| {
|
||||
let (a, b) = (face_of(e.triangles[0]), face_of(e.triangles[1]));
|
||||
if a != b { pair_kept.get(&(a.min(b), a.max(b))).copied().unwrap_or(false) } else { e.angle >= sharp_angle && e.angle <= blade_angle }
|
||||
};
|
||||
let mut edge_set: BTreeMap<(u32, u32), [u32; 2]> = BTreeMap::new();
|
||||
for e in surface.feature_edges().iter().filter(kept) {
|
||||
edge_set.insert((e.vertices[0], e.vertices[1]), e.triangles);
|
||||
}
|
||||
let mut adj: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
|
||||
for &(a, b) in edge_set.keys() {
|
||||
adj.entry(a).or_default().push(b);
|
||||
adj.entry(b).or_default().push(a);
|
||||
}
|
||||
let mut corners: Vec<u32> = adj
|
||||
.iter()
|
||||
.filter(|(v, n)| {
|
||||
if n.len() != 2 {
|
||||
return true;
|
||||
}
|
||||
let (p, a, b) = (xyz(**v), xyz(n[0]), xyz(n[1]));
|
||||
let (u, w) = (sub(p, a), sub(b, p));
|
||||
crate::vec3::angle_deg(u, w) > CORNER_TURN
|
||||
})
|
||||
.map(|(v, _)| *v)
|
||||
.collect();
|
||||
corners.sort_unstable();
|
||||
|
||||
// Patches: flood fill across non-feature edges.
|
||||
let neighbors = surface.neighbors();
|
||||
let mut tri_patch = vec![u32::MAX; tris.len()];
|
||||
let mut patches: Vec<Patch> = Vec::new();
|
||||
for seed in 0..tris.len() {
|
||||
if tri_patch[seed] != u32::MAX {
|
||||
continue;
|
||||
}
|
||||
let id = patches.len() as u32;
|
||||
let mut patch = Patch { triangles: Vec::new(), faces: Vec::new(), area: 0.0, width: f64::INFINITY };
|
||||
tri_patch[seed] = id;
|
||||
let mut stack = vec![seed as u32];
|
||||
while let Some(t) = stack.pop() {
|
||||
patch.triangles.push(t);
|
||||
let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
|
||||
let n = neighbors.get(t as usize).copied().unwrap_or([u32::MAX; 3]);
|
||||
for k in 0..3 {
|
||||
let (a, b) = (tri[k], tri[(k + 1) % 3]);
|
||||
if edge_set.contains_key(&(a.min(b), a.max(b))) {
|
||||
continue;
|
||||
}
|
||||
if let Some(slot) = tri_patch.get_mut(n[k] as usize)
|
||||
&& *slot == u32::MAX
|
||||
{
|
||||
*slot = id;
|
||||
stack.push(n[k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
patch.triangles.sort_unstable();
|
||||
for &t in &patch.triangles {
|
||||
let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
|
||||
patch.area += 0.5 * norm(cross(sub(xyz(tri[1]), xyz(tri[0])), sub(xyz(tri[2]), xyz(tri[0]))));
|
||||
patch.faces.push(surface.tri_face().get(t as usize).copied().unwrap_or(0));
|
||||
}
|
||||
patch.faces.sort_unstable();
|
||||
patch.faces.dedup();
|
||||
let mut perimeter = 0.0;
|
||||
for &t in &patch.triangles {
|
||||
let tri = tris.get(t as usize).copied().unwrap_or([0; 3]);
|
||||
for k in 0..3 {
|
||||
let (a, b) = (tri[k], tri[(k + 1) % 3]);
|
||||
if edge_set.contains_key(&(a.min(b), a.max(b))) {
|
||||
perimeter += norm(sub(xyz(b), xyz(a)));
|
||||
}
|
||||
}
|
||||
}
|
||||
if perimeter > 0.0 {
|
||||
patch.width = 2.0 * patch.area / perimeter;
|
||||
}
|
||||
patches.push(patch);
|
||||
}
|
||||
|
||||
// Curves: chains from corners, then closed loops.
|
||||
let mut used: BTreeMap<(u32, u32), bool> = edge_set.keys().map(|k| (*k, false)).collect();
|
||||
let is_corner = |v: u32| corners.binary_search(&v).is_ok();
|
||||
let mut curves: Vec<Curve> = Vec::new();
|
||||
let trace = |start: u32, next: u32, used: &mut BTreeMap<(u32, u32), bool>, curves: &mut Vec<Curve>| {
|
||||
let mut verts = vec![start];
|
||||
let (mut prev, mut cur) = (start, next);
|
||||
loop {
|
||||
let key = (prev.min(cur), prev.max(cur));
|
||||
if let Some(u) = used.get_mut(&key) {
|
||||
*u = true;
|
||||
}
|
||||
verts.push(cur);
|
||||
if cur == start || is_corner(cur) || verts.len() > edge_set.len() + 2 {
|
||||
break;
|
||||
}
|
||||
let nxt = adj.get(&cur).and_then(|n| n.iter().copied().find(|&w| w != prev && used.get(&(cur.min(w), cur.max(w))) == Some(&false)));
|
||||
match nxt {
|
||||
Some(w) => {
|
||||
prev = cur;
|
||||
cur = w;
|
||||
}
|
||||
None => break,
|
||||
}
|
||||
}
|
||||
let closed = verts.first() == verts.last() && verts.len() > 2;
|
||||
let mut ps: Vec<u32> = Vec::new();
|
||||
for w in verts.windows(2) {
|
||||
if let Some(t) = edge_set.get(&(w[0].min(w[1]), w[0].max(w[1]))) {
|
||||
ps.extend(t.iter().filter_map(|&t| tri_patch.get(t as usize).copied()));
|
||||
}
|
||||
}
|
||||
ps.sort_unstable();
|
||||
ps.dedup();
|
||||
curves.push(Curve { vertices: verts, closed, patches: ps });
|
||||
};
|
||||
for &c in &corners {
|
||||
let ns = adj.get(&c).cloned().unwrap_or_default();
|
||||
for n in ns {
|
||||
if used.get(&(c.min(n), c.max(n))) == Some(&false) {
|
||||
trace(c, n, &mut used, &mut curves);
|
||||
}
|
||||
}
|
||||
}
|
||||
let keys: Vec<(u32, u32)> = used.keys().copied().collect();
|
||||
for k in keys {
|
||||
if used.get(&k) == Some(&false) {
|
||||
trace(k.0, k.1, &mut used, &mut curves);
|
||||
}
|
||||
}
|
||||
Ok(Geometry { surface, tree, boxes, normals, tri_patch, patches, curves, corners, bounds, far })
|
||||
}
|
||||
|
||||
pub fn xyz(&self, v: u32) -> P {
|
||||
self.surface.points().get(v as usize).copied().unwrap_or([0.0; 3])
|
||||
}
|
||||
|
||||
/// The centroid of triangle `t`.
|
||||
pub fn triangle_centre(&self, t: u32) -> P {
|
||||
let [a, b, c] = self.triangle(t);
|
||||
[0, 1, 2].map(|k| (a[k] + b[k] + c[k]) / 3.0)
|
||||
}
|
||||
|
||||
pub fn triangle(&self, t: u32) -> [P; 3] {
|
||||
let tri = self.surface.triangles().get(t as usize).copied().unwrap_or([0; 3]);
|
||||
tri.map(|v| self.xyz(v))
|
||||
}
|
||||
|
||||
/// The closest boundary point to `p` among triangles that pass `keep`: (triangle, point,
|
||||
/// squared distance).
|
||||
pub fn nearest(&self, p: P, keep: impl Fn(u32) -> bool) -> Option<(u32, P, f64)> {
|
||||
let mut best: Option<(u32, P, f64)> = None;
|
||||
let mut stack = vec![0usize];
|
||||
while let Some(n) = stack.pop() {
|
||||
let Some(node) = self.tree.nodes.get(n) else { continue };
|
||||
let bd = node.bounds.dist2(p);
|
||||
if best.is_some_and(|b| bd >= b.2) {
|
||||
continue;
|
||||
}
|
||||
match node.children {
|
||||
Some(l) => {
|
||||
let (dl, dr) = (
|
||||
self.tree.nodes.get(l).map_or(f64::INFINITY, |x| x.bounds.dist2(p)),
|
||||
self.tree.nodes.get(l + 1).map_or(f64::INFINITY, |x| x.bounds.dist2(p)),
|
||||
);
|
||||
if dl <= dr {
|
||||
stack.extend([l + 1, l]);
|
||||
} else {
|
||||
stack.extend([l, l + 1]);
|
||||
}
|
||||
}
|
||||
None => {
|
||||
for &t in self.tree.leaf(node) {
|
||||
if !keep(t) {
|
||||
continue;
|
||||
}
|
||||
let [a, b, c] = self.triangle(t);
|
||||
let q = closest_on_triangle(p, a, b, c);
|
||||
let d = dist2(p, q);
|
||||
if best.is_none_or(|x| d < x.2 || (d == x.2 && t < x.0)) {
|
||||
best = Some((t, q, d));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
best
|
||||
}
|
||||
|
||||
/// The triangles whose boxes overlap `b`, ascending.
|
||||
pub fn overlapping(&self, b: &Aabb) -> Vec<u32> {
|
||||
let mut out = Vec::new();
|
||||
let mut stack = vec![0usize];
|
||||
while let Some(n) = stack.pop() {
|
||||
let Some(node) = self.tree.nodes.get(n) else { continue };
|
||||
if !node.bounds.overlaps(b) {
|
||||
continue;
|
||||
}
|
||||
match node.children {
|
||||
Some(l) => stack.extend([l, l + 1]),
|
||||
None => out.extend(self.tree.leaf(node).iter().copied().filter(|&t| self.boxes.get(t as usize).is_some_and(|x| x.overlaps(b)))),
|
||||
}
|
||||
}
|
||||
out.sort_unstable();
|
||||
out
|
||||
}
|
||||
|
||||
/// The first triangle hit by the ray from `p` along `d` (unit) within `max`, with its
|
||||
/// distance (floating point; used to measure, never to decide).
|
||||
pub fn ray(&self, p: P, d: P, max: f64, skip: u32) -> Option<(u32, f64)> {
|
||||
let q = [0, 1, 2].map(|k| p[k] + d[k] * max);
|
||||
let mut best: Option<(u32, f64)> = None;
|
||||
let mut stack = vec![0usize];
|
||||
while let Some(n) = stack.pop() {
|
||||
let Some(node) = self.tree.nodes.get(n) else { continue };
|
||||
if !node.bounds.meets_segment(p, q) {
|
||||
continue;
|
||||
}
|
||||
match node.children {
|
||||
Some(l) => stack.extend([l, l + 1]),
|
||||
None => {
|
||||
for &t in self.tree.leaf(node) {
|
||||
if t == skip {
|
||||
continue;
|
||||
}
|
||||
let [a, b, c] = self.triangle(t);
|
||||
if let Some(s) = ray_triangle(p, d, a, b, c)
|
||||
&& s > 0.0
|
||||
&& s <= max
|
||||
&& best.is_none_or(|x| s < x.1)
|
||||
{
|
||||
best = Some((t, s));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
best
|
||||
}
|
||||
|
||||
/// Whether `p` is strictly inside the solid, decided exactly: the parity of the boundary
|
||||
/// crossings of a segment from `p` to a far point, tried along fixed directions until one
|
||||
/// meets no edge or vertex and does not start on the boundary. `None` when every direction is
|
||||
/// degenerate (p lies on the boundary) or `p` is not admissible.
|
||||
pub fn inside(&self, p: P) -> Option<bool> {
|
||||
const DIRS: [P; 7] = [
|
||||
[0.5773, 0.6123, 0.5401],
|
||||
[-0.6301, 0.4372, 0.6417],
|
||||
[0.3989, -0.7068, 0.5839],
|
||||
[0.7211, 0.4109, -0.5577],
|
||||
[-0.4402, -0.6005, -0.6675],
|
||||
[0.1351, 0.9013, -0.4116],
|
||||
[-0.8803, 0.1237, 0.4583],
|
||||
];
|
||||
let pp = Point3::new(snap_to_domain(p[0]), snap_to_domain(p[1]), snap_to_domain(p[2])).ok()?;
|
||||
'dirs: for d in DIRS {
|
||||
let q = [0, 1, 2].map(|k| snap_to_domain(p[k] + d[k] * self.far));
|
||||
let Ok(qq) = Point3::new(q[0], q[1], q[2]) else { continue };
|
||||
let mut parity = false;
|
||||
let mut stack = vec![0usize];
|
||||
while let Some(n) = stack.pop() {
|
||||
let Some(node) = self.tree.nodes.get(n) else { continue };
|
||||
if !node.bounds.meets_segment(p, q) {
|
||||
continue;
|
||||
}
|
||||
match node.children {
|
||||
Some(l) => stack.extend([l, l + 1]),
|
||||
None => {
|
||||
for &t in self.tree.leaf(node) {
|
||||
match self.crossing(&pp, &qq, t) {
|
||||
Some(true) => parity = !parity,
|
||||
Some(false) => {}
|
||||
None => continue 'dirs,
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
return Some(parity);
|
||||
}
|
||||
None
|
||||
}
|
||||
|
||||
/// Whether segment `pq` crosses triangle `t` properly (`Some(true)`), misses it
|
||||
/// (`Some(false)`), or touches it degenerately (`None`): exact.
|
||||
fn crossing(&self, p: &Point3, q: &Point3, t: u32) -> Option<bool> {
|
||||
let tri = self.surface.triangles().get(t as usize)?;
|
||||
let a = self.surface.point3(tri[0])?;
|
||||
let b = self.surface.point3(tri[1])?;
|
||||
let c = self.surface.point3(tri[2])?;
|
||||
let (sp, sq) = (orient3d(&a, &b, &c, p), orient3d(&a, &b, &c, q));
|
||||
if sq.is_zero() {
|
||||
// The far end is outside the boundary's box, so never on a triangle; a segment in
|
||||
// the plane is ambiguous.
|
||||
return None;
|
||||
}
|
||||
if sp.is_zero() {
|
||||
// p in the triangle's plane: the segment leaves the plane at p, so it touches the
|
||||
// triangle only if p lies in it.
|
||||
let axis = [Axis::Z, Axis::X, Axis::Y].into_iter().find(|&ax| !orient2d_projected(&a, &b, &c, ax).is_zero())?;
|
||||
let o = orient2d_projected(&a, &b, &c, axis);
|
||||
let inside = [(&a, &b), (&b, &c), (&c, &a)].into_iter().all(|(u, v)| orient2d_projected(u, v, p, axis) != -o);
|
||||
return if inside { None } else { Some(false) };
|
||||
}
|
||||
if sp == sq {
|
||||
return Some(false);
|
||||
}
|
||||
let s = [orient3d(p, q, &a, &b), orient3d(p, q, &b, &c), orient3d(p, q, &c, &a)];
|
||||
if s.iter().any(|x| x.is_zero()) {
|
||||
let nz: Vec<Sign> = s.iter().copied().filter(|x| !x.is_zero()).collect();
|
||||
return if nz.windows(2).all(|w| w[0] == w[1]) { None } else { Some(false) };
|
||||
}
|
||||
Some(s[0] == s[1] && s[1] == s[2])
|
||||
}
|
||||
|
||||
/// The closest point of curve `c` to `p`: (point, segment index, parameter on it).
|
||||
pub fn nearest_on_curve(&self, c: u32, p: P) -> Option<(P, usize, f64)> {
|
||||
let curve = self.curves.get(c as usize)?;
|
||||
let mut best: Option<(P, usize, f64, f64)> = None;
|
||||
for (i, w) in curve.vertices.windows(2).enumerate() {
|
||||
let (q, t) = closest_on_segment(p, self.xyz(w[0]), self.xyz(w[1]));
|
||||
let d = dist2(p, q);
|
||||
if best.is_none_or(|b| d < b.3) {
|
||||
best = Some((q, i, t, d));
|
||||
}
|
||||
}
|
||||
best.map(|(q, i, t, _)| (q, i, t))
|
||||
}
|
||||
}
|
||||
|
||||
/// Möller–Trumbore ray–triangle distance (floating point).
|
||||
fn ray_triangle(o: P, d: P, a: P, b: P, c: P) -> Option<f64> {
|
||||
let (e1, e2) = (sub(b, a), sub(c, a));
|
||||
let h = cross(d, e2);
|
||||
let det = dot(e1, h);
|
||||
if det.abs() < 1e-300 {
|
||||
return None;
|
||||
}
|
||||
let inv = 1.0 / det;
|
||||
let s = sub(o, a);
|
||||
let u = inv * dot(s, h);
|
||||
if !(-1e-12..=1.0 + 1e-12).contains(&u) {
|
||||
return None;
|
||||
}
|
||||
let q = cross(s, e1);
|
||||
let v = inv * dot(d, q);
|
||||
if v < -1e-12 || u + v > 1.0 + 1e-12 {
|
||||
return None;
|
||||
}
|
||||
Some(inv * dot(e2, q))
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
pub(crate) mod tests {
|
||||
use super::*;
|
||||
use crate::surface::SurfaceOptions;
|
||||
|
||||
/// A closed axis-aligned box surface, one B-rep face per side.
|
||||
pub(crate) fn box_surface(lo: P, hi: P) -> BoundarySurface {
|
||||
let c = |i: usize| [if i & 1 == 0 { lo[0] } else { hi[0] }, if i & 2 == 0 { lo[1] } else { hi[1] }, if i & 4 == 0 { lo[2] } else { hi[2] }];
|
||||
let pts: Vec<P> = (0..8).map(c).collect();
|
||||
let quads = [[0, 2, 3, 1], [4, 5, 7, 6], [0, 1, 5, 4], [2, 6, 7, 3], [0, 4, 6, 2], [1, 3, 7, 5]];
|
||||
let mut tris = Vec::new();
|
||||
let mut face = Vec::new();
|
||||
for (f, q) in quads.iter().enumerate() {
|
||||
tris.push([q[0], q[1], q[2]]);
|
||||
tris.push([q[0], q[2], q[3]]);
|
||||
face.extend([f as u32; 2]);
|
||||
}
|
||||
BoundarySurface::new(&pts, &tris, &face, &SurfaceOptions::default()).unwrap()
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn box_feature_graph_inside_and_nearest() {
|
||||
let s = box_surface([0.0; 3], [2.0, 3.0, 4.0]);
|
||||
let g = Geometry::new(&s, 30.0, 160.0).unwrap();
|
||||
assert_eq!(g.corners.len(), 8);
|
||||
assert_eq!(g.curves.len(), 12);
|
||||
assert_eq!(g.patches.len(), 6);
|
||||
assert!(g.curves.iter().all(|c| c.vertices.len() == 2 && c.patches.len() == 2 && !c.closed));
|
||||
assert_eq!(g.inside([1.0, 1.5, 2.0]), Some(true));
|
||||
assert_eq!(g.inside([3.0, 1.5, 2.0]), Some(false));
|
||||
assert_eq!(g.inside([1.0, 1.5, 0.0]), None, "on the boundary");
|
||||
assert_eq!(g.inside([0.0, 0.0, 0.0]), None, "a corner");
|
||||
let (_, q, d) = g.nearest([1.0, 1.0, -1.0], |_| true).unwrap();
|
||||
assert!((d - 1.0).abs() < 1e-12 && q[2] == 0.0);
|
||||
let top = g.tri_patch[2];
|
||||
let (_, q, _) = g.nearest([1.0, 1.0, -1.0], |t| g.tri_patch[t as usize] == top).unwrap();
|
||||
assert_eq!(q[2], 4.0, "restricted to the top patch");
|
||||
let (t, s) = g.ray([1.0, 1.5, 2.0], [0.0, 0.0, 1.0], 10.0, u32::MAX).unwrap();
|
||||
assert!((s - 2.0).abs() < 1e-12 && g.tri_patch[t as usize] == top);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,376 @@
|
||||
//! The adaptive grid: a box of root cubes, each refined as a 27-tree (a leaf splits into
|
||||
//! 3 × 3 × 3), strongly balanced and closed under the template patterns, then cut into a
|
||||
//! conforming hexahedral mesh by [`super::templates`].
|
||||
//!
|
||||
//! # Balance and closure
|
||||
//! *Strong balance*: leaves that share a vertex differ by at most one level, so a leaf's finer
|
||||
//! neighbours are exactly one level finer and a leaf next to a coarser one never has a finer one
|
||||
//! on that side. A vertex of a leaf at level `L` is *marked* when a finer leaf touches it, that
|
||||
//! is, when one of the eight level-`L` cells around the vertex is split. *Closure*: a leaf whose
|
||||
//! marks fit no template (two diagonal corners of a face, three of a face, …) is split too.
|
||||
//! Balancing and closing repeat until neither splits anything; both only split, so this ends,
|
||||
//! and every round runs over the leaves in sorted order, so the result is the same on every run.
|
||||
//!
|
||||
//! # Coordinates
|
||||
//! Everything is integer. A cell at level `L` has index `i` (per axis, in units of its own edge)
|
||||
//! and spans `6 · 3^(D − L)` lattice units, `D` the deepest level allowed; template vertices are
|
||||
//! in sixths of a leaf, so every vertex is a lattice point and vertices shared by neighbours are
|
||||
//! found by exact lookup.
|
||||
|
||||
use std::collections::HashMap;
|
||||
|
||||
use super::templates::{self, Pattern, Sym};
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// A cell: its level and index (in units of its edge at that level).
|
||||
pub(crate) type Key = (u8, [i64; 3]);
|
||||
|
||||
/// The deepest level a tree may have (`3^12` lattice steps per root cell edge, × 6).
|
||||
pub(crate) const MAX_LEVEL: u8 = 12;
|
||||
|
||||
fn pow3(n: u8) -> i64 {
|
||||
3i64.pow(u32::from(n))
|
||||
}
|
||||
|
||||
/// What covers a cell-sized region at some level.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
enum Cover {
|
||||
/// A leaf (at that level or coarser).
|
||||
Leaf(Key),
|
||||
/// The region is split at that level (finer leaves inside).
|
||||
Split,
|
||||
/// Outside the root box.
|
||||
Outside,
|
||||
}
|
||||
|
||||
/// The tree. `nodes` maps every cell that exists to whether it is a leaf.
|
||||
#[derive(Debug, Clone)]
|
||||
pub(crate) struct Tree {
|
||||
/// Root cubes per axis.
|
||||
pub roots: [i64; 3],
|
||||
/// The deepest level allowed.
|
||||
pub depth: u8,
|
||||
nodes: HashMap<Key, bool>,
|
||||
leaves: usize,
|
||||
max_leaves: usize,
|
||||
}
|
||||
|
||||
/// The sort order of [`Tree::sorted_leaves`]: level, then z, y, x.
|
||||
fn order(k: &Key) -> (u8, i64, i64, i64) {
|
||||
(k.0, k.1[2], k.1[1], k.1[0])
|
||||
}
|
||||
|
||||
/// The 26 neighbour offsets, in a fixed order.
|
||||
fn offsets() -> impl Iterator<Item = [i64; 3]> {
|
||||
(0..27).map(|k| [k % 3 - 1, (k / 3) % 3 - 1, k / 9 - 1]).filter(|d| *d != [0, 0, 0])
|
||||
}
|
||||
|
||||
impl Tree {
|
||||
/// One leaf per root cube.
|
||||
pub fn new(roots: [i64; 3], depth: u8, max_leaves: usize) -> Result<Tree> {
|
||||
if roots.iter().any(|&n| n < 1) || depth > MAX_LEVEL {
|
||||
return Err(Error::Invalid(format!("grid of {roots:?} root cells to depth {depth}")));
|
||||
}
|
||||
let count = roots.iter().try_fold(1usize, |a, &n| usize::try_from(n).ok().and_then(|n| a.checked_mul(n)));
|
||||
let count = count.filter(|&c| c <= max_leaves).ok_or_else(|| Error::Limit(format!("{roots:?} root cells exceed {max_leaves} cells")))?;
|
||||
let mut nodes = HashMap::with_capacity(count);
|
||||
for z in 0..roots[2] {
|
||||
for y in 0..roots[1] {
|
||||
for x in 0..roots[0] {
|
||||
nodes.insert((0, [x, y, z]), true);
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(Tree { roots, depth, nodes, leaves: count, max_leaves })
|
||||
}
|
||||
|
||||
pub fn leaf_count(&self) -> usize {
|
||||
self.leaves
|
||||
}
|
||||
|
||||
/// Whether `key` is a leaf.
|
||||
pub fn is_leaf(&self, key: &Key) -> bool {
|
||||
self.nodes.get(key).copied().unwrap_or(false)
|
||||
}
|
||||
|
||||
fn inside(&self, (level, i): Key) -> bool {
|
||||
let n = pow3(level);
|
||||
(0..3).all(|k| i[k] >= 0 && i[k] < self.roots[k] * n)
|
||||
}
|
||||
|
||||
/// What covers the level-`L` cell `key`.
|
||||
fn cover(&self, key: Key) -> Cover {
|
||||
if !self.inside(key) {
|
||||
return Cover::Outside;
|
||||
}
|
||||
let (level, i) = key;
|
||||
for l in (0..=level).rev() {
|
||||
let s = pow3(level - l);
|
||||
let k = (l, [i[0].div_euclid(s), i[1].div_euclid(s), i[2].div_euclid(s)]);
|
||||
match self.nodes.get(&k) {
|
||||
Some(true) => return Cover::Leaf(k),
|
||||
Some(false) => return Cover::Split,
|
||||
None => {}
|
||||
}
|
||||
}
|
||||
Cover::Outside
|
||||
}
|
||||
|
||||
/// Split a leaf into its 27 children.
|
||||
pub fn split(&mut self, key: Key) -> Result<()> {
|
||||
let level = key.0;
|
||||
if !self.is_leaf(&key) {
|
||||
return Ok(());
|
||||
}
|
||||
if level >= self.depth {
|
||||
return Err(Error::Limit(format!("refinement needs more than {} levels", self.depth)));
|
||||
}
|
||||
if self.leaves + 26 > self.max_leaves {
|
||||
return Err(Error::Limit(format!("the grid needs more than {} cells", self.max_leaves)));
|
||||
}
|
||||
self.nodes.insert(key, false);
|
||||
for child in Self::children(key) {
|
||||
self.nodes.insert(child, true);
|
||||
}
|
||||
self.leaves += 26;
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// The leaves in sorted order (level, then z, y, x).
|
||||
pub fn sorted_leaves(&self) -> Vec<Key> {
|
||||
let mut v: Vec<Key> = self.nodes.iter().filter(|(_, leaf)| **leaf).map(|(k, _)| *k).collect();
|
||||
v.sort_by_key(order);
|
||||
v
|
||||
}
|
||||
|
||||
/// The children of a split cell.
|
||||
pub fn children(key: Key) -> impl Iterator<Item = Key> {
|
||||
let (level, i) = key;
|
||||
(0..27).map(move |k| (level + 1, [3 * i[0] + k % 3, 3 * i[1] + (k / 3) % 3, 3 * i[2] + k / 9]))
|
||||
}
|
||||
|
||||
/// Split until leaves sharing a vertex differ by at most one level, checking only the leaves
|
||||
/// `work` (in sorted order) and the children of the cells it splits: a split only makes
|
||||
/// neighbours finer, so only new leaves can be too fine for a neighbour, and the result is
|
||||
/// the unique coarsest balanced refinement whatever the order. Every split cell goes to `log`.
|
||||
fn balance_from(&mut self, mut work: Vec<Key>, log: &mut Vec<Key>) -> Result<()> {
|
||||
work.reverse();
|
||||
while let Some(key) = work.pop() {
|
||||
if !self.is_leaf(&key) || key.0 < 2 {
|
||||
continue;
|
||||
}
|
||||
let (level, i) = key;
|
||||
for d in offsets() {
|
||||
if let Cover::Leaf(k) = self.cover((level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]]))
|
||||
&& k.0 + 1 < level
|
||||
{
|
||||
self.split(k)?;
|
||||
log.push(k);
|
||||
let mut kids: Vec<Key> = Self::children(k).collect();
|
||||
kids.reverse();
|
||||
work.extend(kids);
|
||||
work.push(key);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// The leaves whose marks a split of `key` may change: its children and the leaves beside it
|
||||
/// on its level.
|
||||
fn touched_by(&self, key: Key, out: &mut Vec<Key>) {
|
||||
out.extend(Self::children(key));
|
||||
let (level, i) = key;
|
||||
out.extend(offsets().map(|d| (level, [i[0] + d[0], i[1] + d[1], i[2] + d[2]])).filter(|k| self.is_leaf(k)));
|
||||
}
|
||||
|
||||
/// The marked corners of a leaf (bit `c` for corner `c`: bit 0 x, bit 1 y, bit 2 z).
|
||||
pub fn marks(&self, (level, i): Key) -> u8 {
|
||||
let mut mask = 0u8;
|
||||
for c in 0..8usize {
|
||||
let o = [(c & 1) as i64, ((c >> 1) & 1) as i64, ((c >> 2) & 1) as i64];
|
||||
let split = (0..8usize).any(|e| {
|
||||
let n = [i[0] + o[0] - 1 + (e & 1) as i64, i[1] + o[1] - 1 + ((e >> 1) & 1) as i64, i[2] + o[2] - 1 + ((e >> 2) & 1) as i64];
|
||||
self.nodes.get(&(level, n)) == Some(&false)
|
||||
});
|
||||
if split {
|
||||
mask |= 1 << c;
|
||||
}
|
||||
}
|
||||
mask
|
||||
}
|
||||
|
||||
/// Balance, then split every leaf whose marks fit no template (and every leaf marked at all
|
||||
/// eight corners, whose template is the split), until nothing changes. `max_rounds` bounds the
|
||||
/// rounds. Each round re-checks only the leaves the previous splits can have changed, so a
|
||||
/// round costs what it splits, not the size of the tree.
|
||||
pub fn close(&mut self, max_rounds: usize, skip: &dyn Fn(Key) -> bool) -> Result<()> {
|
||||
let mut fresh = self.sorted_leaves();
|
||||
let mut touched = fresh.clone();
|
||||
for _ in 0..max_rounds {
|
||||
let mut log = Vec::new();
|
||||
self.balance_from(std::mem::take(&mut fresh), &mut log)?;
|
||||
for &k in &log {
|
||||
self.touched_by(k, &mut touched);
|
||||
}
|
||||
touched.retain(|k| self.is_leaf(k));
|
||||
touched.sort_by_key(order);
|
||||
touched.dedup();
|
||||
let bad: Vec<Key> = touched
|
||||
.drain(..)
|
||||
.filter(|&k| !skip(k))
|
||||
.filter(|&k| !matches!(templates::classify(self.marks(k)), Some(Pattern::None | Pattern::Corner | Pattern::Edge | Pattern::Face)))
|
||||
.collect();
|
||||
if bad.is_empty() {
|
||||
return Ok(());
|
||||
}
|
||||
for &k in &bad {
|
||||
self.split(k)?;
|
||||
}
|
||||
for &k in &bad {
|
||||
fresh.extend(Self::children(k));
|
||||
self.touched_by(k, &mut touched);
|
||||
}
|
||||
fresh.sort_by_key(order);
|
||||
}
|
||||
Err(Error::Limit(format!("grid closure did not settle in {max_rounds} rounds")))
|
||||
}
|
||||
|
||||
/// Lattice units per edge of a level-`L` cell.
|
||||
pub fn span(&self, level: u8) -> i64 {
|
||||
6 * pow3(self.depth - level)
|
||||
}
|
||||
|
||||
/// Cut every leaf by its template into a conforming mesh on the integer lattice. Fails if a
|
||||
/// leaf's marks fit no template (call [`Tree::close`] first).
|
||||
pub fn mesh(&self, skip: &dyn Fn(Key) -> bool) -> Result<LatticeMesh> {
|
||||
let syms: Vec<Sym> = templates::symmetries();
|
||||
let canon: Vec<(Pattern, Vec<templates::Hex>)> =
|
||||
[Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face].into_iter().map(|p| (p, templates::template(p))).collect();
|
||||
let mut out = LatticeMesh::default();
|
||||
let mut index: HashMap<[i64; 3], u32> = HashMap::new();
|
||||
for key in self.sorted_leaves() {
|
||||
if skip(key) {
|
||||
continue;
|
||||
}
|
||||
let mask = self.marks(key);
|
||||
let pattern = templates::classify(mask).filter(|p| *p != Pattern::Full);
|
||||
let (pattern, sym) = pattern
|
||||
.and_then(|p| templates::orientation(&syms, p, mask).map(|s| (p, s)))
|
||||
.ok_or_else(|| Error::Failed(format!("grid cell {key:?} has marks {mask:#010b} that fit no template")))?;
|
||||
let hexes = canon.iter().find(|(p, _)| *p == pattern).map(|(_, h)| h.as_slice()).unwrap_or(&[]);
|
||||
let (level, i) = key;
|
||||
let unit = pow3(self.depth - level);
|
||||
for h in hexes {
|
||||
let mut g = h.map(|p| {
|
||||
let q = sym.apply(p);
|
||||
[0, 1, 2].map(|k| (6 * i[k] + q[k]) * unit)
|
||||
});
|
||||
if templates::corner0_orientation(&g) < 0 {
|
||||
g = templates::reversed(g);
|
||||
}
|
||||
let mut ids = [0u32; 8];
|
||||
for (id, p) in ids.iter_mut().zip(g) {
|
||||
*id = match index.get(&p) {
|
||||
Some(&v) => v,
|
||||
None => {
|
||||
let v = u32::try_from(out.points.len()).map_err(|_| Error::Limit("too many grid points".into()))?;
|
||||
index.insert(p, v);
|
||||
out.points.push(p);
|
||||
v
|
||||
}
|
||||
};
|
||||
}
|
||||
out.hexes.push(ids);
|
||||
out.level.push(level);
|
||||
}
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
}
|
||||
|
||||
/// A conforming hexahedral mesh on the integer lattice.
|
||||
#[derive(Debug, Clone, Default)]
|
||||
pub(crate) struct LatticeMesh {
|
||||
pub points: Vec<[i64; 3]>,
|
||||
pub hexes: Vec<[u32; 8]>,
|
||||
/// The tree level of the leaf each hexahedron came from.
|
||||
pub level: Vec<u8>,
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use std::collections::HashMap;
|
||||
|
||||
use super::*;
|
||||
use crate::quality::{Validity, element_validity};
|
||||
use fieldcraft_model::ElementKind;
|
||||
|
||||
/// A tiny deterministic generator for test trees.
|
||||
struct Lcg(u64);
|
||||
impl Lcg {
|
||||
fn next(&mut self) -> u64 {
|
||||
self.0 = self.0.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407);
|
||||
self.0 >> 33
|
||||
}
|
||||
}
|
||||
|
||||
fn check_conforming(tree: &Tree) {
|
||||
let m = tree.mesh(&|_| false).unwrap();
|
||||
let mut faces: HashMap<[u32; 4], usize> = HashMap::new();
|
||||
let mut volume = 0.0;
|
||||
for h in &m.hexes {
|
||||
let xyz: Vec<[f64; 3]> = h.iter().map(|&v| m.points[v as usize].map(|c| c as f64)).collect();
|
||||
let c = element_validity(ElementKind::Hex8, &xyz).unwrap();
|
||||
assert_eq!(c.validity, Validity::Valid);
|
||||
volume += c.volume;
|
||||
for f in ElementKind::Hex8.faces() {
|
||||
let mut q = [h[f[0]], h[f[1]], h[f[2]], h[f[3]]];
|
||||
q.sort();
|
||||
*faces.entry(q).or_default() += 1;
|
||||
}
|
||||
}
|
||||
let side = tree.span(0) as f64;
|
||||
let expected = side.powi(3) * (tree.roots[0] * tree.roots[1] * tree.roots[2]) as f64;
|
||||
assert!((volume - expected).abs() < 1e-6 * expected, "volume {volume} vs {expected}");
|
||||
let hi = tree.roots.map(|n| n * tree.span(0));
|
||||
for (q, n) in faces {
|
||||
let on_box = (0..3).any(|k| q.iter().all(|&v| m.points[v as usize][k] == 0) || q.iter().all(|&v| m.points[v as usize][k] == hi[k]));
|
||||
assert_eq!(n, if on_box { 1 } else { 2 }, "face {q:?}");
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn random_trees_close_into_conforming_valid_meshes() {
|
||||
let mut rng = Lcg(7);
|
||||
for case in 0..12 {
|
||||
let roots = [1 + (rng.next() % 2) as i64, 1 + (rng.next() % 2) as i64, 1];
|
||||
let mut tree = Tree::new(roots, 3, 1_000_000).unwrap();
|
||||
for _ in 0..(1 + rng.next() % 6) {
|
||||
let leaves = tree.sorted_leaves();
|
||||
let k = leaves[(rng.next() as usize) % leaves.len()];
|
||||
if k.0 < 3 {
|
||||
tree.split(k).unwrap();
|
||||
}
|
||||
}
|
||||
tree.close(100, &|_| false).unwrap();
|
||||
for k in tree.sorted_leaves() {
|
||||
assert!(templates::classify(tree.marks(k)).is_some(), "case {case}");
|
||||
}
|
||||
check_conforming(&tree);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn caps_are_errors() {
|
||||
assert!(Tree::new([0, 1, 1], 2, 10).is_err());
|
||||
assert!(Tree::new([10, 10, 10], 2, 999).is_err());
|
||||
let mut t = Tree::new([1, 1, 1], 1, 100).unwrap();
|
||||
t.split((0, [0, 0, 0])).unwrap();
|
||||
assert!(t.split((1, [0, 0, 0])).is_err(), "deeper than the depth");
|
||||
let mut t = Tree::new([1, 1, 1], 3, 30).unwrap();
|
||||
t.split((0, [0, 0, 0])).unwrap();
|
||||
assert!(t.split((1, [0, 0, 0])).is_err(), "more cells than allowed");
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,282 @@
|
||||
//! The topological map from the boundary's feature graph (corners, curves, patches) onto the
|
||||
//! core's boundary quads, after X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based
|
||||
//! hexahedral meshing*, CGF 38(5), 2019, §3.3 (paper only; the authors' code has no licence and
|
||||
//! is not read). Gao et al. trace the curves first and read the patches off them; here the
|
||||
//! patches come first and the curves are their borders, which keeps every region on the side of
|
||||
//! the part it faces:
|
||||
//!
|
||||
//! 1. **Label.** Each quad takes the patch its outward normal meets first (within two quad
|
||||
//! sizes), or else the patch nearest its centre.
|
||||
//! 2. **Clean.** A quad whose label differs from all its edge neighbours takes their majority
|
||||
//! label; then every patch keeps only its largest connected group of quads, the other
|
||||
//! groups joining the label they border most. Both repeat a bounded number of times.
|
||||
//! 3. **Check.** The map must be a bijection: every patch has a region, and two regions border
|
||||
//! each other exactly where the part has a curve between their patches.
|
||||
//! 4. **Features.** Region borders are the curves (each border edge takes the curve between its
|
||||
//! two patches nearest to it); a vertex where three or more regions meet is a corner and takes
|
||||
//! the nearest corner of the part around those patches. Corners that are not region junctions
|
||||
//! (where one curve turns sharply) go to the nearest vertex on their curve's border.
|
||||
//!
|
||||
//! Any failure names where the grid is too coarse for the features; the caller refines there and
|
||||
//! tries again (§3.3, "Bijective map").
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet};
|
||||
|
||||
use super::geom::{Geometry, closest_on_segment, dist2};
|
||||
use super::smooth::Constraint;
|
||||
use super::topo::quad_edges;
|
||||
use crate::vec3::{P, cross, norm, sub};
|
||||
|
||||
/// Rounds of label cleaning.
|
||||
const CLEAN_ROUNDS: usize = 16;
|
||||
|
||||
/// Where and why the map failed.
|
||||
#[derive(Debug, Clone)]
|
||||
pub(crate) struct MapFailure {
|
||||
pub at: Vec<P>,
|
||||
pub why: String,
|
||||
}
|
||||
|
||||
/// The map: a constraint per quad vertex and a patch per quad.
|
||||
#[derive(Debug, Clone)]
|
||||
pub(crate) struct FeatureMap {
|
||||
pub role: BTreeMap<u32, Constraint>,
|
||||
pub quad_patch: Vec<u32>,
|
||||
}
|
||||
|
||||
/// Map the feature graph onto `quads` (the core's boundary, outward) over `xyz`.
|
||||
pub(crate) fn map_features(g: &Geometry, quads: &[[u32; 4]], xyz: &dyn Fn(u32) -> P) -> Result<FeatureMap, MapFailure> {
|
||||
let n = quads.len();
|
||||
let centre = |q: usize| {
|
||||
let mut c = [0.0; 3];
|
||||
for &v in &quads[q] {
|
||||
let p = xyz(v);
|
||||
for k in 0..3 {
|
||||
c[k] += 0.25 * p[k];
|
||||
}
|
||||
}
|
||||
c
|
||||
};
|
||||
let edges = quad_edges(quads);
|
||||
let mut nbr: Vec<Vec<u32>> = vec![Vec::new(); n];
|
||||
for qs in edges.values() {
|
||||
if let [a, b] = qs.as_slice() {
|
||||
nbr[*a as usize].push(*b);
|
||||
nbr[*b as usize].push(*a);
|
||||
}
|
||||
}
|
||||
|
||||
// 1. Label: the patch the quad faces (a ray along its outward normal), else the nearest.
|
||||
let mut label: Vec<u32> = (0..n)
|
||||
.map(|q| {
|
||||
let c = centre(q);
|
||||
let p = quads[q].map(xyz);
|
||||
let nrm = cross(sub(p[2], p[0]), sub(p[3], p[1]));
|
||||
let l = norm(nrm);
|
||||
let reach = 2.0 * (norm(sub(p[2], p[0])) + norm(sub(p[3], p[1])));
|
||||
let hit = if l > 0.0 { g.ray(c, nrm.map(|x| x / l), reach, u32::MAX).map(|(t, _)| t) } else { None };
|
||||
hit.or_else(|| g.nearest(c, |_| true).map(|(t, _, _)| t)).and_then(|t| g.tri_patch.get(t as usize).copied()).unwrap_or(0)
|
||||
})
|
||||
.collect();
|
||||
|
||||
// A patch no quad faces (a narrow face the core's boundary passes by) takes the quad
|
||||
// nearest to it, which the cleaning keeps.
|
||||
let mut seeded: BTreeSet<usize> = BTreeSet::new();
|
||||
let have: BTreeSet<u32> = label.iter().copied().collect();
|
||||
for p in 0..g.patches.len() as u32 {
|
||||
if have.contains(&p) {
|
||||
continue;
|
||||
}
|
||||
let best = (0..n)
|
||||
.filter(|q| !seeded.contains(q))
|
||||
.filter_map(|q| g.nearest(centre(q), |t| g.tri_patch.get(t as usize) == Some(&p)).map(|x| (q, x.2)))
|
||||
.min_by(|a, b| a.1.total_cmp(&b.1).then(a.0.cmp(&b.0)));
|
||||
if let Some((q, _)) = best {
|
||||
label[q] = p;
|
||||
seeded.insert(q);
|
||||
}
|
||||
}
|
||||
// 2. Clean.
|
||||
for _ in 0..CLEAN_ROUNDS {
|
||||
let mut changed = false;
|
||||
for q in 0..n {
|
||||
let mut count: BTreeMap<u32, usize> = BTreeMap::new();
|
||||
for &m in &nbr[q] {
|
||||
*count.entry(label[m as usize]).or_default() += 1;
|
||||
}
|
||||
let same = count.get(&label[q]).copied().unwrap_or(0);
|
||||
if same == 0
|
||||
&& !nbr[q].is_empty()
|
||||
&& !seeded.contains(&q)
|
||||
&& let Some((&best, _)) = count.iter().filter(|(l, _)| **l != label[q]).max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0)))
|
||||
{
|
||||
label[q] = best;
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
let mut comp = vec![usize::MAX; n];
|
||||
let mut comps: Vec<(u32, Vec<usize>)> = Vec::new();
|
||||
for s in 0..n {
|
||||
if comp[s] != usize::MAX {
|
||||
continue;
|
||||
}
|
||||
let id = comps.len();
|
||||
let mut members = vec![s];
|
||||
comp[s] = id;
|
||||
let mut i = 0;
|
||||
while i < members.len() {
|
||||
let q = members[i];
|
||||
i += 1;
|
||||
for &m in &nbr[q] {
|
||||
let m = m as usize;
|
||||
if comp[m] == usize::MAX && label[m] == label[s] {
|
||||
comp[m] = id;
|
||||
members.push(m);
|
||||
}
|
||||
}
|
||||
}
|
||||
comps.push((label[s], members));
|
||||
}
|
||||
let mut largest: BTreeMap<u32, usize> = BTreeMap::new();
|
||||
for (i, (l, m)) in comps.iter().enumerate() {
|
||||
let e = largest.entry(*l).or_insert(i);
|
||||
if comps[*e].1.len() < m.len() {
|
||||
*e = i;
|
||||
}
|
||||
}
|
||||
for (i, (l, members)) in comps.iter().enumerate() {
|
||||
if largest.get(l) == Some(&i) {
|
||||
continue;
|
||||
}
|
||||
let mut count: BTreeMap<u32, usize> = BTreeMap::new();
|
||||
for &q in members {
|
||||
for &m in &nbr[q] {
|
||||
if comp[m as usize] != i {
|
||||
*count.entry(label[m as usize]).or_default() += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
if let Some((&best, _)) = count.iter().max_by(|a, b| a.1.cmp(b.1).then(b.0.cmp(a.0))) {
|
||||
for &q in members {
|
||||
label[q] = best;
|
||||
}
|
||||
changed = true;
|
||||
}
|
||||
}
|
||||
if !changed {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// 3. Check: every patch has a region; borders match curves.
|
||||
let present: BTreeSet<u32> = label.iter().copied().collect();
|
||||
for (p, patch) in g.patches.iter().enumerate() {
|
||||
if !present.contains(&(p as u32)) {
|
||||
let at: Vec<P> = patch.triangles.iter().take(64).map(|&t| g.triangle_centre(t)).collect();
|
||||
return Err(MapFailure { at, why: format!("patch {p} (B-rep faces {:?}) has no region", patch.faces) });
|
||||
}
|
||||
}
|
||||
let mut curves_between: BTreeMap<(u32, u32), Vec<u32>> = BTreeMap::new();
|
||||
for (i, c) in g.curves.iter().enumerate() {
|
||||
match c.patches.as_slice() {
|
||||
[a, b] => curves_between.entry((*a, *b)).or_default().push(i as u32),
|
||||
[a] => curves_between.entry((*a, *a)).or_default().push(i as u32),
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
let mut role: BTreeMap<u32, Constraint> = BTreeMap::new();
|
||||
let mut border_pairs: BTreeSet<(u32, u32)> = BTreeSet::new();
|
||||
let mut vertex_labels: BTreeMap<u32, BTreeSet<u32>> = BTreeMap::new();
|
||||
for (q, quad) in quads.iter().enumerate() {
|
||||
for &v in quad {
|
||||
vertex_labels.entry(v).or_default().insert(label[q]);
|
||||
}
|
||||
}
|
||||
for ((a, b), qs) in &edges {
|
||||
let (q0, q1) = match qs.as_slice() {
|
||||
[x, y] => (*x as usize, *y as usize),
|
||||
_ => continue,
|
||||
};
|
||||
let (la, lb) = (label[q0].min(label[q1]), label[q0].max(label[q1]));
|
||||
if la == lb {
|
||||
continue;
|
||||
}
|
||||
let Some(cs) = curves_between.get(&(la, lb)) else {
|
||||
return Err(MapFailure { at: vec![xyz(*a), xyz(*b)], why: format!("regions of patches {la} and {lb} meet where the part has no curve") });
|
||||
};
|
||||
border_pairs.insert((la, lb));
|
||||
let mid = [0, 1, 2].map(|k| 0.5 * (xyz(*a)[k] + xyz(*b)[k]));
|
||||
let c = nearest_curve(g, cs, mid);
|
||||
for v in [*a, *b] {
|
||||
role.entry(v).or_insert(Constraint::Curve(c));
|
||||
}
|
||||
}
|
||||
for (pair, cs) in &curves_between {
|
||||
if pair.0 != pair.1 && !border_pairs.contains(pair) {
|
||||
let at: Vec<P> = cs.iter().flat_map(|&c| g.curves[c as usize].vertices.iter().map(|&v| g.xyz(v)).collect::<Vec<_>>()).collect();
|
||||
return Err(MapFailure { at, why: format!("the regions of patches {} and {} do not meet along their curve", pair.0, pair.1) });
|
||||
}
|
||||
}
|
||||
|
||||
// 4. Corners.
|
||||
let mut taken: BTreeSet<u32> = BTreeSet::new();
|
||||
let junctions: Vec<u32> = vertex_labels.iter().filter(|(_, l)| l.len() >= 3).map(|(v, _)| *v).collect();
|
||||
for &v in &junctions {
|
||||
let p = xyz(v);
|
||||
let labels = &vertex_labels[&v];
|
||||
let best = g
|
||||
.corners
|
||||
.iter()
|
||||
.copied()
|
||||
.filter(|c| !taken.contains(c) && corner_patches(g, *c).is_superset(labels))
|
||||
.min_by(|&a, &b| dist2(g.xyz(a), p).total_cmp(&dist2(g.xyz(b), p)).then(a.cmp(&b)));
|
||||
match best {
|
||||
Some(c) => {
|
||||
taken.insert(c);
|
||||
role.insert(v, Constraint::Corner(c));
|
||||
}
|
||||
None => return Err(MapFailure { at: vec![p], why: format!("regions {labels:?} meet where the part has no corner") }),
|
||||
}
|
||||
}
|
||||
for &c in &g.corners {
|
||||
if taken.contains(&c) {
|
||||
continue;
|
||||
}
|
||||
let p = g.xyz(c);
|
||||
let cp = corner_patches(g, c);
|
||||
let best = role
|
||||
.iter()
|
||||
.filter(|(_, r)| matches!(r, Constraint::Curve(_)))
|
||||
.map(|(v, _)| *v)
|
||||
.filter(|v| vertex_labels.get(v).is_some_and(|l| l.is_subset(&cp)))
|
||||
.min_by(|&a, &b| dist2(xyz(a), p).total_cmp(&dist2(xyz(b), p)).then(a.cmp(&b)));
|
||||
match best {
|
||||
Some(v) => {
|
||||
role.insert(v, Constraint::Corner(c));
|
||||
}
|
||||
None => return Err(MapFailure { at: vec![p], why: format!("the feature corner at {p:?} has no vertex") }),
|
||||
}
|
||||
}
|
||||
for (q, quad) in quads.iter().enumerate() {
|
||||
for &v in quad {
|
||||
role.entry(v).or_insert(Constraint::Patch(label[q]));
|
||||
}
|
||||
}
|
||||
Ok(FeatureMap { role, quad_patch: label })
|
||||
}
|
||||
|
||||
/// The patches around a corner of the part (on either side of the curves that end there).
|
||||
fn corner_patches(g: &Geometry, c: u32) -> BTreeSet<u32> {
|
||||
g.curves.iter().filter(|cv| cv.vertices.first() == Some(&c) || cv.vertices.last() == Some(&c)).flat_map(|cv| cv.patches.iter().copied()).collect()
|
||||
}
|
||||
|
||||
/// The curve of `cs` nearest to `p`.
|
||||
fn nearest_curve(g: &Geometry, cs: &[u32], p: P) -> u32 {
|
||||
cs.iter().copied().min_by(|&a, &b| curve_distance(g, a, p).total_cmp(&curve_distance(g, b, p)).then(a.cmp(&b))).unwrap_or(0)
|
||||
}
|
||||
|
||||
fn curve_distance(g: &Geometry, c: u32, p: P) -> f64 {
|
||||
g.curves.get(c as usize).map_or(f64::INFINITY, |cv| {
|
||||
cv.vertices.windows(2).map(|w| dist2(p, closest_on_segment(p, g.xyz(w[0]), g.xyz(w[1])).0)).fold(f64::INFINITY, f64::min)
|
||||
})
|
||||
}
|
||||
@@ -0,0 +1,801 @@
|
||||
//! Unstructured all-hexahedral meshing of any closed boundary: an adaptive grid cut into
|
||||
//! conforming hexahedra, a core of grid elements inside the part, two boundary layers that carry
|
||||
//! the part's faces, sharp curves and corners, and an optimisation that fits them to the
|
||||
//! geometry while keeping every element valid.
|
||||
//!
|
||||
//! # Pipeline
|
||||
//! 1. **Frame.** The grid is aligned with the part: the normals of its two largest planar
|
||||
//! patches when they are orthogonal and cover more area than the coordinate axes do, the
|
||||
//! coordinate axes otherwise (octree methods are orientation dependent; N. Pietroni et al.,
|
||||
//! *Hex-mesh generation and processing: a survey*, ACM TOG 42(2), 2022, §4.6).
|
||||
//! 2. **Grid** ([`grid`]). Root cubes of the target size over the part's box, refined as a
|
||||
//! 27-tree where a B-rep face curves (its normals turn by more than
|
||||
//! [`OctreeOptions::curvature_angle`] inside a cell), where the part is thin (fewer than
|
||||
//! [`OctreeOptions::thickness_cells`] cells across, measured along the inward normal of each
|
||||
//! boundary triangle) and where a patch is narrow, then strongly balanced and closed under the
|
||||
//! transition templates.
|
||||
//! 3. **Conforming hexahedra** ([`templates`]): the 3-refinement node-marking templates
|
||||
//! (R. Schneiders, 1996; Y. Ito, A. M. Shih, B. K. Soni, IJNME 77(13), 2009), so every grid
|
||||
//! leaf becomes 1, 4, 11, 22 or 27 hexahedra and the grid is conforming by construction.
|
||||
//! 4. **Core.** Elements whose eight corners are inside the part (exact ray parity on the
|
||||
//! boundary triangles) and at least [`OctreeOptions::core_offset`] of their cell size from it.
|
||||
//! Elements are dropped until the core's boundary is a 2-manifold and every boundary vertex can
|
||||
//! move outward (its quads' normals lie in one half-space: the buffer rule of H. Tong,
|
||||
//! E. Halilaj, Y. J. Zhang, *HybridOctree_Hex*, arXiv 2401.05984, 2024); its components and
|
||||
//! genus must match the part's shells.
|
||||
//! 5. **Feature map** ([`map`]): the part's patches, sharp curves and corners mapped onto the
|
||||
//! core's boundary quads (X. Gao, H. Shen, D. Panozzo, *Feature preserving octree-based
|
||||
//! hexahedral meshing*, CGF 38(5), 2019, §3.3), patches first.
|
||||
//! 6. **Padding and pillowing** ([`topo`]): one layer of hexahedra on the core's boundary (so no
|
||||
//! element keeps two faces on the boundary), then each patch's layer is pillowed, which inserts
|
||||
//! a ring of quads along every sharp curve on both sides: a curve then runs between quads that
|
||||
//! each have one edge on it, so straight and curved features are captured without degenerate
|
||||
//! quads (Gao et al. 2019, §3.4, Fig. 9).
|
||||
//! 7. **Fit** ([`smooth`]): the foldover-free energy of V. Garanzha et al., *Foldover-free maps
|
||||
//! in 50 lines of code*, ACM TOG 40(4), 2021, plus a fitting term that draws every boundary
|
||||
//! vertex onto its corner, curve or patch with a growing weight (Gao et al. 2019, §3.5); the
|
||||
//! barrier keeps every element valid while the boundary moves. Each boundary vertex then goes
|
||||
//! exactly onto its feature where its elements stay valid; where they would not (a blade, a
|
||||
//! near-tangent contact) the mesh stays off by [`OctreeReport::boundary_distance`].
|
||||
//! 8. **Gate.** The worst elements are improved by compass search, then every element is proven
|
||||
//! valid by the certified Jacobian check (Bernstein bounds, [`crate::quality`]; corner
|
||||
//! Jacobians alone do not prove a trilinear hexahedron valid, J. Marschner et al., CGF 39(5),
|
||||
//! 2020); elements that are not are repaired locally first. A failure in 4–8 names where the
|
||||
//! grid is too coarse; the leaves there are split once more and the mesh is rebuilt, at most
|
||||
//! [`OctreeOptions::attempts`] times, and a rebuild may not grow the grid more than
|
||||
//! [`RETRY_GROWTH`] times (Gao et al. 2019, §3.3).
|
||||
//!
|
||||
//! # Why 3-refinement templates
|
||||
//! The dual schemes of M. Livesu, L. Pitzalis, G. Cherchi, *Optimal dual schemes for adaptive
|
||||
//! grid based hexmeshing* (ACM TOG 41(2), 2022) give a simpler singular structure, but their
|
||||
//! installation depends on a paired grid whose transition vertices the reference pipeline
|
||||
//! (L. Pitzalis et al., *Generalized adaptive refinement for grid-based hexahedral meshing*, ACM
|
||||
//! TOG 40(6), 2021) selects with an integer linear program, and on heuristics that are retried
|
||||
//! with other parameters when they fail. The node-marking templates need only the balance and
|
||||
//! closure rules and are proven conforming by an exhaustive test. Their limit is the one
|
||||
//! Schneiders states: they close convex refined regions only, so a leaf next to a concave corner
|
||||
//! of a finer region is split, and the split spreads to the region's orthogonal hull (H. Tong,
|
||||
//! Y. J. Zhang, *MCHex*, arXiv 2511.02064, 2025, meet the same limit with the same templates).
|
||||
//! Curvature therefore refines only [`OctreeOptions::curvature_levels`] levels; deeper levels
|
||||
//! follow thickness and narrow faces, whose refined regions are mostly convex. Schneiders'
|
||||
//! modified edge rule (every edge with a marked end split in three) or the dual schemes are the
|
||||
//! upgrade path.
|
||||
//!
|
||||
//! # Determinism and bounds
|
||||
//! The same boundary and options give the same mesh: grid cells are processed in sorted order,
|
||||
//! maps are only looked up, the optimiser sweeps vertices in index order. Grid cells, elements,
|
||||
//! refinement depth, closure rounds, optimiser sweeps and attempts are capped.
|
||||
|
||||
mod geom;
|
||||
mod grid;
|
||||
mod map;
|
||||
mod smooth;
|
||||
mod templates;
|
||||
mod topo;
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests;
|
||||
|
||||
use std::cell::{Cell, RefCell};
|
||||
use std::collections::{BTreeMap, BTreeSet, HashMap, HashSet};
|
||||
|
||||
use fieldcraft_model::ElementKind;
|
||||
|
||||
pub use crate::hex::BoundaryFace;
|
||||
use crate::hex::face_flux;
|
||||
use crate::quality::{GateOptions, Metric, Validity, element_quality, validity_gate_checks};
|
||||
use crate::size::{SizeField, checked_size};
|
||||
use crate::surface::BoundarySurface;
|
||||
use crate::vec3::{P, cross, dot, norm, sub};
|
||||
use crate::{Error, Result};
|
||||
|
||||
use geom::{Aabb, Geometry};
|
||||
use grid::Tree;
|
||||
use smooth::{Constraint, Optimizer};
|
||||
|
||||
/// Options and caps of [`octree_hex_mesh`].
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct OctreeOptions {
|
||||
/// Refinement levels below the root cells (each one a third of the size), at most 8.
|
||||
pub max_depth: u8,
|
||||
/// The most grid cells and the most output elements.
|
||||
pub max_elements: usize,
|
||||
/// A cell is split where the normals of one B-rep face inside it turn by more than this
|
||||
/// (degrees).
|
||||
pub curvature_angle: f64,
|
||||
/// Levels below the root that curvature may refine (see the module docs).
|
||||
pub curvature_levels: u8,
|
||||
/// A cell is split where the part is less than this many cells thick.
|
||||
pub thickness_cells: f64,
|
||||
/// Core elements keep this fraction of their cell size from the boundary.
|
||||
pub core_offset: f64,
|
||||
/// Optimiser sweeps over all vertices.
|
||||
pub max_sweeps: usize,
|
||||
/// Builds, counting the rebuilds after refining where a stage failed (1 = no retry).
|
||||
pub attempts: usize,
|
||||
/// Boundary edges whose faces meet at more than this angle (degrees between normals) are
|
||||
/// sharp features the mesh captures with its own edges; smoother face boundaries are not.
|
||||
pub sharp_angle: f64,
|
||||
/// Edges whose normals turn by more than this (degrees: thin blades, near-tangent contacts)
|
||||
/// are not captured; the mesh rounds them.
|
||||
pub blade_angle: f64,
|
||||
}
|
||||
|
||||
impl Default for OctreeOptions {
|
||||
fn default() -> Self {
|
||||
OctreeOptions {
|
||||
max_depth: 2,
|
||||
max_elements: 5_000_000,
|
||||
curvature_angle: 25.0,
|
||||
curvature_levels: 1,
|
||||
thickness_cells: 3.0,
|
||||
core_offset: 0.3,
|
||||
max_sweeps: 300,
|
||||
attempts: 4,
|
||||
sharp_angle: 30.0,
|
||||
blade_angle: 160.0,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// What [`octree_hex_mesh`] did.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct OctreeReport {
|
||||
/// Builds run (1 + rebuilds after a failure).
|
||||
pub attempts: usize,
|
||||
/// Grid leaves of the last build.
|
||||
pub grid_cells: usize,
|
||||
/// Core elements (from the grid) and boundary-layer elements.
|
||||
pub core_elements: usize,
|
||||
pub layer_elements: usize,
|
||||
/// The grid's axes.
|
||||
pub frame: [[f64; 3]; 3],
|
||||
/// Sharp corners, curves and patches of the boundary.
|
||||
pub corners: usize,
|
||||
pub curves: usize,
|
||||
pub patches: usize,
|
||||
/// Optimiser sweeps of the last build.
|
||||
pub sweeps: usize,
|
||||
/// Smallest and mean corner scaled Jacobian.
|
||||
pub min_scaled_jacobian: f64,
|
||||
pub mean_scaled_jacobian: f64,
|
||||
/// Exact volume of the elements and the boundary's enclosed volume.
|
||||
pub volume: f64,
|
||||
pub surface_volume: f64,
|
||||
/// The largest distance of a boundary vertex from its corner, curve or patch: zero when
|
||||
/// every boundary vertex lies on the geometry; where the geometry is sharper than valid
|
||||
/// elements can follow (blades, slits), the mesh stays this far off.
|
||||
pub boundary_distance: f64,
|
||||
/// Why earlier builds were refined.
|
||||
pub retries: Vec<String>,
|
||||
}
|
||||
|
||||
/// An all-hex8 mesh.
|
||||
#[derive(Debug, Clone, PartialEq)]
|
||||
pub struct OctreeHexMesh {
|
||||
pub nodes: Vec<[f64; 3]>,
|
||||
/// Model node order; every element is proven valid.
|
||||
pub hexes: Vec<[u32; 8]>,
|
||||
pub boundary: Vec<BoundaryFace>,
|
||||
pub report: OctreeReport,
|
||||
}
|
||||
|
||||
/// A rebuild after a failure may have at most this many times the grid cells of the build
|
||||
/// before…
|
||||
pub const RETRY_GROWTH: usize = 4;
|
||||
/// … or this many cells, whichever is more.
|
||||
const RETRY_MIN_CELLS: usize = 250_000;
|
||||
/// The fitting weights of the successive optimisation stages.
|
||||
const FIT_STAGES: &[f64] = &[1e-3, 1e-2, 1e-1, 1.0, 10.0, 100.0];
|
||||
/// A boundary vertex goes exactly onto its feature only if its elements keep at least this
|
||||
/// corner scaled Jacobian.
|
||||
const MIN_SNAP_QUALITY: f64 = 0.05;
|
||||
/// Elements below this corner scaled Jacobian are improved by compass search.
|
||||
const IMPROVE_BELOW: f64 = 0.3;
|
||||
/// Closure rounds of the grid (each splits at least one leaf; the cell cap bounds them too).
|
||||
const CLOSURE_ROUNDS: usize = 1000;
|
||||
/// Rounds of the core clean-up (each drops at least one element).
|
||||
const CORE_ROUNDS: usize = 1000;
|
||||
|
||||
/// A failure that refining at `at` may cure.
|
||||
struct Retry {
|
||||
at: Vec<P>,
|
||||
why: String,
|
||||
}
|
||||
|
||||
enum Attempt {
|
||||
Done(Box<OctreeHexMesh>),
|
||||
Retry(Retry),
|
||||
}
|
||||
|
||||
/// Mesh the solid bounded by `surface` with hexahedra of about `field`'s size (the root cell
|
||||
/// size is the field at the centre of the part's box; finer cells follow the part's curvature,
|
||||
/// thickness and narrow faces). See the module docs.
|
||||
pub fn octree_hex_mesh(surface: &BoundarySurface, field: &dyn SizeField, options: &OctreeOptions) -> Result<OctreeHexMesh> {
|
||||
check_options(options)?;
|
||||
let g = Geometry::new(surface, options.sharp_angle, options.blade_angle)?;
|
||||
let frame = choose_frame(&g);
|
||||
let mut lo = [f64::INFINITY; 3];
|
||||
let mut hi = [f64::NEG_INFINITY; 3];
|
||||
for &p in surface.points() {
|
||||
let u = to_frame(&frame, p);
|
||||
for k in 0..3 {
|
||||
lo[k] = lo[k].min(u[k]);
|
||||
hi[k] = hi[k].max(u[k]);
|
||||
}
|
||||
}
|
||||
let centre_u = [0, 1, 2].map(|k| 0.5 * (lo[k] + hi[k]));
|
||||
let h = checked_size(field, from_frame(&frame, centre_u))?;
|
||||
let mut roots = [0i64; 3];
|
||||
for k in 0..3 {
|
||||
let n = ((hi[k] - lo[k]) / h + 1.0).ceil();
|
||||
if !(n.is_finite() && n >= 1.0 && n <= options.max_elements as f64) {
|
||||
return Err(Error::Limit(format!("a size of {h} needs {n} root cells along an axis")));
|
||||
}
|
||||
roots[k] = n as i64;
|
||||
}
|
||||
let origin = [0, 1, 2].map(|k| centre_u[k] - 0.5 * roots[k] as f64 * h);
|
||||
let ctx = Context {
|
||||
g: &g,
|
||||
field,
|
||||
options,
|
||||
frame,
|
||||
origin,
|
||||
h,
|
||||
roots,
|
||||
thickness: triangle_thickness(&g),
|
||||
outside_memo: RefCell::new(HashMap::new()),
|
||||
last_leaves: Cell::new(0),
|
||||
too_thin: Cell::new(f64::INFINITY),
|
||||
};
|
||||
let mut hints: Vec<P> = Vec::new();
|
||||
let mut retries: Vec<String> = Vec::new();
|
||||
for attempt in 0..options.attempts {
|
||||
let cap = if attempt == 0 {
|
||||
options.max_elements
|
||||
} else {
|
||||
options.max_elements.min(RETRY_GROWTH.saturating_mul(ctx.last_leaves.get()).max(RETRY_MIN_CELLS))
|
||||
};
|
||||
let built = match ctx.build(&hints, cap) {
|
||||
Err(Error::Limit(why)) if attempt > 0 => {
|
||||
retries.push(format!("refining for another build needs more than {cap} grid cells ({why})"));
|
||||
break;
|
||||
}
|
||||
other => other?,
|
||||
};
|
||||
match built {
|
||||
Attempt::Done(mut m) => {
|
||||
m.report.attempts = attempt + 1;
|
||||
m.report.retries = retries;
|
||||
return Ok(*m);
|
||||
}
|
||||
Attempt::Retry(r) => {
|
||||
retries.push(r.why);
|
||||
hints.extend(r.at);
|
||||
}
|
||||
}
|
||||
}
|
||||
Err(Error::Failed(format!("octree hex meshing failed after {} builds: {}", retries.len(), retries.join("; "))))
|
||||
}
|
||||
|
||||
fn check_options(o: &OctreeOptions) -> Result<()> {
|
||||
let ok = o.max_depth <= 8
|
||||
&& o.curvature_levels <= o.max_depth
|
||||
&& o.max_elements >= 1
|
||||
&& o.curvature_angle.is_finite()
|
||||
&& o.curvature_angle > 0.0
|
||||
&& o.thickness_cells.is_finite()
|
||||
&& o.thickness_cells >= 0.0
|
||||
&& o.core_offset.is_finite()
|
||||
&& (0.0..1.0).contains(&o.core_offset)
|
||||
&& (1..=16).contains(&o.attempts)
|
||||
&& o.sharp_angle.is_finite()
|
||||
&& o.sharp_angle > 0.0
|
||||
&& o.blade_angle.is_finite()
|
||||
&& o.blade_angle > o.sharp_angle
|
||||
&& o.blade_angle <= 180.0
|
||||
&& o.max_sweeps <= 100_000;
|
||||
if ok { Ok(()) } else { Err(Error::Invalid(format!("octree options out of range: {o:?}"))) }
|
||||
}
|
||||
|
||||
fn from_frame(f: &[[f64; 3]; 3], u: P) -> P {
|
||||
[0, 1, 2].map(|k| u[0] * f[0][k] + u[1] * f[1][k] + u[2] * f[2][k])
|
||||
}
|
||||
|
||||
fn to_frame(f: &[[f64; 3]; 3], p: P) -> P {
|
||||
[dot(p, f[0]), dot(p, f[1]), dot(p, f[2])]
|
||||
}
|
||||
|
||||
/// The grid axes (rows): the coordinate axes, or the normals of the two largest orthogonal planar
|
||||
/// patches when more boundary area is aligned with them.
|
||||
fn choose_frame(g: &Geometry) -> [[f64; 3]; 3] {
|
||||
let tol = 1f64.to_radians().sin();
|
||||
let normal = |t: u32| g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
|
||||
let mut planar: Vec<(f64, usize, P)> = Vec::new();
|
||||
for (i, p) in g.patches.iter().enumerate() {
|
||||
let Some(&t0) = p.triangles.first() else { continue };
|
||||
let n0 = normal(t0);
|
||||
if p.triangles.iter().all(|&t| norm(cross(normal(t), n0)) < tol) {
|
||||
planar.push((p.area, i, n0));
|
||||
}
|
||||
}
|
||||
planar.sort_by(|a, b| b.0.total_cmp(&a.0).then(a.1.cmp(&b.1)));
|
||||
let identity = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
|
||||
let Some(&(_, _, a0)) = planar.first() else { return identity };
|
||||
let Some(&(_, _, a1)) = planar.iter().find(|x| dot(x.2, a0).abs() < tol) else { return identity };
|
||||
let a1 = {
|
||||
let v = sub(a1, a0.map(|x| x * dot(a1, a0)));
|
||||
let l = norm(v);
|
||||
v.map(|x| x / l)
|
||||
};
|
||||
let candidate = [a0, a1, cross(a0, a1)];
|
||||
let score = |f: &[[f64; 3]; 3]| -> f64 {
|
||||
(0..g.surface.triangles().len() as u32)
|
||||
.filter(|&t| f.iter().any(|a| norm(cross(normal(t), *a)) < tol))
|
||||
.map(|t| {
|
||||
let [a, b, c] = g.triangle(t);
|
||||
0.5 * norm(cross(sub(b, a), sub(c, a)))
|
||||
})
|
||||
.sum()
|
||||
};
|
||||
if score(&candidate) > score(&identity) * (1.0 + 1e-9) { candidate } else { identity }
|
||||
}
|
||||
|
||||
/// The distance through the part from each boundary triangle: a ray from its centre along the
|
||||
/// inward normal to the next triangle (`∞` when none).
|
||||
fn triangle_thickness(g: &Geometry) -> Vec<f64> {
|
||||
let far = 2.0 * norm(sub(g.bounds.hi, g.bounds.lo)) + 1.0;
|
||||
(0..g.surface.triangles().len() as u32)
|
||||
.map(|t| {
|
||||
let n = g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
|
||||
g.ray(g.triangle_centre(t), n.map(|x| -x), far, t).map_or(f64::INFINITY, |(_, s)| s)
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn centroid(points: &[P], hx: &[u32; 8]) -> P {
|
||||
[0, 1, 2].map(|k| hx.iter().map(|&v| points.get(v as usize).map_or(0.0, |p| p[k])).sum::<f64>() / 8.0)
|
||||
}
|
||||
|
||||
fn edge_mean(points: &[P], hx: &[u32; 8]) -> f64 {
|
||||
let x = hx.map(|v| points.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
let edges = [(0, 1), (1, 2), (2, 3), (3, 0), (4, 5), (5, 6), (6, 7), (7, 4), (0, 4), (1, 5), (2, 6), (3, 7)];
|
||||
let m = edges.iter().map(|&(a, b)| norm(sub(x[a], x[b]))).sum::<f64>() / 12.0;
|
||||
if m > 0.0 { m } else { 1.0 }
|
||||
}
|
||||
|
||||
/// The part being meshed and the state carried from one build to the next.
|
||||
struct Context<'a> {
|
||||
g: &'a Geometry<'a>,
|
||||
field: &'a dyn SizeField,
|
||||
options: &'a OctreeOptions,
|
||||
frame: [[f64; 3]; 3],
|
||||
origin: P,
|
||||
h: f64,
|
||||
roots: [i64; 3],
|
||||
thickness: Vec<f64>,
|
||||
outside_memo: RefCell<HashMap<grid::Key, bool>>,
|
||||
/// Grid leaves of the last build.
|
||||
last_leaves: Cell<usize>,
|
||||
/// The thinnest wall the deepest cells cannot cross in `thickness_cells` (∞ when none).
|
||||
too_thin: Cell<f64>,
|
||||
}
|
||||
|
||||
/// The mesh while its boundary layers are built: points, elements, their target sizes and how
|
||||
/// each point may move.
|
||||
struct Work {
|
||||
points: Vec<P>,
|
||||
hexes: Vec<[u32; 8]>,
|
||||
size: Vec<f64>,
|
||||
role: Vec<Constraint>,
|
||||
core_count: usize,
|
||||
}
|
||||
|
||||
impl Context<'_> {
|
||||
fn edge(&self, level: u8) -> f64 {
|
||||
self.h / 3f64.powi(i32::from(level))
|
||||
}
|
||||
|
||||
/// The world box (inflated by `pad` cell edges) and the centre of a cell.
|
||||
fn cell_box(&self, (level, i): grid::Key, pad: f64) -> (Aabb, P) {
|
||||
let e = self.edge(level);
|
||||
let mut b = Aabb::EMPTY;
|
||||
for c in 0..8 {
|
||||
b.add(from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + ((c >> k) & 1) as f64) * e)));
|
||||
}
|
||||
let centre = from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + (i[k] as f64 + 0.5) * e));
|
||||
let p = pad * e;
|
||||
(Aabb { lo: b.lo.map(|x| x - p), hi: b.hi.map(|x| x + p) }, centre)
|
||||
}
|
||||
|
||||
fn needs_split(&self, key: grid::Key) -> Result<bool> {
|
||||
let e = self.edge(key.0);
|
||||
let (bx, centre) = self.cell_box(key, 0.05);
|
||||
if e > checked_size(self.field, centre)? * (1.0 + 1e-9) {
|
||||
return Ok(true);
|
||||
}
|
||||
let reach = 0.5 * 3f64.sqrt() * e * 1.05;
|
||||
let cos = self.options.curvature_angle.to_radians().cos();
|
||||
let mut first: BTreeMap<u32, P> = BTreeMap::new();
|
||||
for t in self.g.overlapping(&bx) {
|
||||
let [a, b, c] = self.g.triangle(t);
|
||||
if geom::dist2(centre, geom::closest_on_triangle(centre, a, b, c)) > reach * reach {
|
||||
continue;
|
||||
}
|
||||
let n = self.g.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
|
||||
let face = self.g.surface.tri_face().get(t as usize).copied().unwrap_or(0);
|
||||
let f = *first.entry(face).or_insert(n);
|
||||
if key.0 < self.options.curvature_levels && dot(f, n) < cos {
|
||||
return Ok(true);
|
||||
}
|
||||
if let Some(&th) = self.thickness.get(t as usize)
|
||||
&& th < self.options.thickness_cells * e
|
||||
{
|
||||
if key.0 + 1 >= self.options.max_depth && th < self.options.thickness_cells * e / 3.0 {
|
||||
self.too_thin.set(self.too_thin.get().min(th));
|
||||
}
|
||||
return Ok(true);
|
||||
}
|
||||
// Narrow patches (a hole's wall, a thin rim) need cells across them.
|
||||
if self.g.tri_patch.get(t as usize).and_then(|&p| self.g.patches.get(p as usize)).is_some_and(|pt| pt.width < 2.0 * e) {
|
||||
return Ok(true);
|
||||
}
|
||||
}
|
||||
Ok(false)
|
||||
}
|
||||
|
||||
/// Whether a cell lies wholly outside the part (its centre outside and farther from the
|
||||
/// boundary than its half diagonal): no core element can come from it, so its transition
|
||||
/// pattern does not matter and it is not meshed. Memoised.
|
||||
fn outside(&self, key: grid::Key) -> bool {
|
||||
if let Some(&v) = self.outside_memo.borrow().get(&key) {
|
||||
return v;
|
||||
}
|
||||
let (_, centre) = self.cell_box(key, 0.0);
|
||||
let half = 0.5 * 3f64.sqrt() * self.edge(key.0);
|
||||
let far = self.g.nearest(centre, |_| true).is_some_and(|x| x.2 > half * half * (1.0 + 1e-9));
|
||||
let v = far && self.g.inside(centre) == Some(false);
|
||||
self.outside_memo.borrow_mut().insert(key, v);
|
||||
v
|
||||
}
|
||||
|
||||
/// The level-`level` cells within one cell of a hint point.
|
||||
fn hint_cells(&self, hints: &[P], level: u8) -> HashSet<grid::Key> {
|
||||
let e = self.edge(level);
|
||||
let mut out = HashSet::new();
|
||||
for &p in hints {
|
||||
let u = to_frame(&self.frame, p);
|
||||
let i = [0, 1, 2].map(|k| ((u[k] - self.origin[k]) / e).floor());
|
||||
if !i.iter().all(|x| x.is_finite() && x.abs() < 1e15) {
|
||||
continue;
|
||||
}
|
||||
for d in 0..27i64 {
|
||||
out.insert((level, [i[0] as i64 + d % 3 - 1, i[1] as i64 + (d / 3) % 3 - 1, i[2] as i64 + d / 9 - 1]));
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// The refined, balanced and closed grid; leaves near `hints` (where an earlier build
|
||||
/// failed) are split once more.
|
||||
fn tree(&self, hints: &[P], max_leaves: usize) -> Result<Tree> {
|
||||
let mut tree = Tree::new(self.roots, self.options.max_depth, max_leaves)?;
|
||||
for level in 0..self.options.max_depth {
|
||||
for key in tree.sorted_leaves() {
|
||||
if key.0 == level && self.needs_split(key)? {
|
||||
tree.split(key)?;
|
||||
}
|
||||
}
|
||||
}
|
||||
if !hints.is_empty() {
|
||||
let near: Vec<HashSet<grid::Key>> = (0..self.options.max_depth).map(|l| self.hint_cells(hints, l)).collect();
|
||||
let split: Vec<grid::Key> =
|
||||
tree.sorted_leaves().into_iter().filter(|k| near.get(usize::from(k.0)).is_some_and(|s| s.contains(k))).collect();
|
||||
for key in split {
|
||||
tree.split(key)?;
|
||||
}
|
||||
}
|
||||
tree.close(CLOSURE_ROUNDS, &|k| self.outside(k))?;
|
||||
Ok(tree)
|
||||
}
|
||||
|
||||
fn build(&self, hints: &[P], max_leaves: usize) -> Result<Attempt> {
|
||||
let g = self.g;
|
||||
let tree = self.tree(hints, max_leaves)?;
|
||||
self.last_leaves.set(tree.leaf_count());
|
||||
let lattice = tree.mesh(&|k| self.outside(k))?;
|
||||
if lattice.hexes.len() > self.options.max_elements {
|
||||
return Err(Error::Limit(format!("the grid has {} elements; the limit is {}", lattice.hexes.len(), self.options.max_elements)));
|
||||
}
|
||||
let unit = self.h / tree.span(0) as f64;
|
||||
let points: Vec<P> = lattice.points.iter().map(|q| from_frame(&self.frame, [0, 1, 2].map(|k| self.origin[k] + q[k] as f64 * unit))).collect();
|
||||
|
||||
// Core.
|
||||
let dist: Vec<f64> = points.iter().map(|&p| g.nearest(p, |_| true).map_or(0.0, |x| x.2.sqrt())).collect();
|
||||
let keep = self.core(&lattice, &points, &dist);
|
||||
let core: Vec<usize> = (0..lattice.hexes.len()).filter(|&i| keep[i]).collect();
|
||||
if core.is_empty() {
|
||||
let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect();
|
||||
return Ok(Attempt::Retry(Retry { at, why: "no grid element lies inside the part".into() }));
|
||||
}
|
||||
let mut w = Work {
|
||||
hexes: core.iter().map(|&i| lattice.hexes[i]).collect(),
|
||||
size: core.iter().map(|&i| self.edge(lattice.level[i])).collect(),
|
||||
role: Vec::new(),
|
||||
points,
|
||||
core_count: core.len(),
|
||||
};
|
||||
let qb = topo::boundary(&w.hexes);
|
||||
let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect();
|
||||
if let Some(r) = check_topology(g, &q) {
|
||||
let thin = self.too_thin.get();
|
||||
if thin.is_finite() {
|
||||
return Err(Error::Failed(format!(
|
||||
"{}: the part has walls {thin:.3} thick, thinner than {} cells of the deepest level ({:.3}); use a smaller size, more levels, or a structured (swept) hex mesh",
|
||||
r.why,
|
||||
self.options.thickness_cells,
|
||||
self.edge(self.options.max_depth)
|
||||
)));
|
||||
}
|
||||
return Ok(Attempt::Retry(r));
|
||||
}
|
||||
|
||||
// Feature map, padding and pillows.
|
||||
let fmap = {
|
||||
let pts = &w.points;
|
||||
match map::map_features(g, &q, &|v| pts.get(v as usize).copied().unwrap_or([0.0; 3])) {
|
||||
Ok(m) => m,
|
||||
Err(f) => return Ok(Attempt::Retry(Retry { at: f.at, why: f.why })),
|
||||
}
|
||||
};
|
||||
self.layers(&mut w, &qb, &fmap, &dist)?;
|
||||
if w.hexes.len() > self.options.max_elements {
|
||||
return Err(Error::Limit(format!("the mesh needs {} elements; the limit is {}", w.hexes.len(), self.options.max_elements)));
|
||||
}
|
||||
|
||||
// Fit, then improve and repair.
|
||||
let mut opt = Optimizer::new(g, std::mem::take(&mut w.points), &w.hexes, &w.size, w.role.clone());
|
||||
let mut ran = 0;
|
||||
let sweeps = (self.options.max_sweeps / FIT_STAGES.len()).max(1);
|
||||
for &weight in FIT_STAGES {
|
||||
opt.fit = weight;
|
||||
ran = opt.run(sweeps, 1e-4);
|
||||
}
|
||||
let (_, boundary_distance) = opt.snap(MIN_SNAP_QUALITY);
|
||||
let inverted: Vec<P> = w
|
||||
.hexes
|
||||
.iter()
|
||||
.filter(|hx| smooth::min_scaled_jacobian(&hx.map(|v| opt.points[v as usize])) <= 0.0)
|
||||
.map(|hx| centroid(&opt.points, hx))
|
||||
.collect();
|
||||
if !inverted.is_empty() {
|
||||
let why = format!("{} elements stay inverted after untangling", inverted.len());
|
||||
return Ok(Attempt::Retry(Retry { at: inverted, why }));
|
||||
}
|
||||
opt.improve_worst(IMPROVE_BELOW, 4);
|
||||
let uncertified = opt.repair_certified(6);
|
||||
if !uncertified.is_empty() {
|
||||
let at: Vec<P> = uncertified.iter().filter_map(|&e| w.hexes.get(e as usize)).map(|hx| centroid(&opt.points, hx)).collect();
|
||||
let why = format!("{} elements fail the certified validity check", at.len());
|
||||
return Ok(Attempt::Retry(Retry { at, why }));
|
||||
}
|
||||
w.points = opt.points;
|
||||
self.finish(&w, &tree, ran, boundary_distance)
|
||||
}
|
||||
|
||||
/// The core: elements wholly inside the part, cleaned to a 2-manifold boundary whose
|
||||
/// vertices can all move outward.
|
||||
fn core(&self, lattice: &grid::LatticeMesh, points: &[P], dist: &[f64]) -> Vec<bool> {
|
||||
let inside: Vec<bool> = points.iter().map(|&p| self.g.inside(p) == Some(true)).collect();
|
||||
let mut keep: Vec<bool> = lattice
|
||||
.hexes
|
||||
.iter()
|
||||
.zip(&lattice.level)
|
||||
.map(|(hx, &l)| {
|
||||
let floor = self.options.core_offset * self.edge(l);
|
||||
hx.iter().all(|&v| inside.get(v as usize) == Some(&true) && dist.get(v as usize).is_some_and(|&d| d >= floor))
|
||||
})
|
||||
.collect();
|
||||
for _ in 0..CORE_ROUNDS {
|
||||
let mut drop = topo::nonmanifold(&lattice.hexes, &keep);
|
||||
if drop.is_empty() {
|
||||
drop = topo::hidden_vertices(&lattice.hexes, &keep, &|v| points.get(v as usize).copied().unwrap_or([0.0; 3]), 0.1);
|
||||
}
|
||||
if drop.is_empty() {
|
||||
break;
|
||||
}
|
||||
for d in drop {
|
||||
if let Some(k) = keep.get_mut(d as usize) {
|
||||
*k = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
keep
|
||||
}
|
||||
|
||||
/// The boundary layers: padding on the core's boundary quads `qb`, whose outer copies take the
|
||||
/// map's roles and start halfway to the boundary along the core's normals, then a pillow
|
||||
/// under each patch, whose copies start a quarter of the way into their elements. Every
|
||||
/// element starts valid or nearly so; the fit untangles the rest.
|
||||
fn layers(&self, w: &mut Work, qb: &[(topo::FaceRef, [u32; 4])], fmap: &map::FeatureMap, dist: &[f64]) -> Result<()> {
|
||||
let q: Vec<[u32; 4]> = qb.iter().map(|x| x.1).collect();
|
||||
let core_count = w.core_count;
|
||||
let copy = topo::pad(&mut w.points, &mut w.hexes, &q)?;
|
||||
w.role = vec![Constraint::Free; w.points.len()];
|
||||
let mut normal: BTreeMap<u32, P> = BTreeMap::new();
|
||||
for f in &q {
|
||||
let p = f.map(|v| w.points[v as usize]);
|
||||
let n = cross(sub(p[2], p[0]), sub(p[3], p[1]));
|
||||
let l = norm(n);
|
||||
if l > 0.0 {
|
||||
for &v in f {
|
||||
let e = normal.entry(v).or_insert([0.0; 3]);
|
||||
for k in 0..3 {
|
||||
e[k] += n[k] / l;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
for (&v, &c) in © {
|
||||
if let Some(r) = fmap.role.get(&v) {
|
||||
w.role[c as usize] = *r;
|
||||
}
|
||||
let n = normal.get(&v).copied().unwrap_or([0.0; 3]);
|
||||
let l = norm(n);
|
||||
let d = 0.5 * dist.get(v as usize).copied().unwrap_or(0.0);
|
||||
if l > 0.0 {
|
||||
w.points[c as usize] = [0, 1, 2].map(|k| w.points[v as usize][k] + d * n[k] / l);
|
||||
}
|
||||
}
|
||||
for ((owner, _), _) in qb {
|
||||
w.size.push(0.5 * w.size.get(*owner as usize).copied().unwrap_or(self.h));
|
||||
}
|
||||
let patches: BTreeSet<u32> = fmap.quad_patch.iter().copied().collect();
|
||||
for p in patches {
|
||||
let inside: Vec<bool> =
|
||||
(0..w.hexes.len()).map(|i| i >= core_count && i < core_count + q.len() && fmap.quad_patch.get(i - core_count) == Some(&p)).collect();
|
||||
let before = w.hexes.len();
|
||||
let centroids: Vec<P> = w.hexes.iter().map(|hx| centroid(&w.points, hx)).collect();
|
||||
let copies = topo::pillow(&mut w.points, &mut w.hexes, &inside)?;
|
||||
w.role.resize(w.points.len(), Constraint::Free);
|
||||
for pc in &copies {
|
||||
let orig = w.points[pc.original as usize];
|
||||
let n = pc.elements.len().max(1) as f64;
|
||||
let mut c = [0.0; 3];
|
||||
for &e in &pc.elements {
|
||||
let ce = centroids.get(e as usize).copied().unwrap_or(orig);
|
||||
for k in 0..3 {
|
||||
c[k] += ce[k] / n;
|
||||
}
|
||||
}
|
||||
w.points[pc.copy as usize] = [0, 1, 2].map(|k| orig[k] + 0.25 * (c[k] - orig[k]));
|
||||
w.role[pc.copy as usize] = match w.role.get(pc.original as usize).copied().unwrap_or(Constraint::Free) {
|
||||
Constraint::Free => Constraint::Free,
|
||||
_ => Constraint::Patch(p),
|
||||
};
|
||||
}
|
||||
for hx in &w.hexes[before..] {
|
||||
w.size.push(edge_mean(&w.points, hx));
|
||||
}
|
||||
}
|
||||
Ok(())
|
||||
}
|
||||
|
||||
/// The gate and the output: every element proven valid, nodes renumbered, boundary faces with
|
||||
/// their B-rep faces.
|
||||
fn finish(&self, w: &Work, tree: &Tree, sweeps: usize, boundary_distance: f64) -> Result<Attempt> {
|
||||
let g = self.g;
|
||||
// One certified pass: every element's quality (the same metric as mesh.quality and
|
||||
// mesh.hex.generate report) and its Jacobian check for the gate.
|
||||
let mut bad: Vec<P> = Vec::new();
|
||||
let mut checks = Vec::with_capacity(w.hexes.len());
|
||||
let (mut sj_min, mut sj_sum) = (f64::INFINITY, 0.0);
|
||||
for (e, hx) in w.hexes.iter().enumerate() {
|
||||
let q = element_quality(ElementKind::Hex8, &hx.map(|v| w.points[v as usize]))?;
|
||||
let check = q.jacobian.ok_or_else(|| Error::Failed("an element was not certified".into()))?;
|
||||
if check.validity != Validity::Valid {
|
||||
bad.push(centroid(&w.points, hx));
|
||||
}
|
||||
let sj = q.get(Metric::ScaledJacobian).unwrap_or(-1.0);
|
||||
sj_min = sj_min.min(sj);
|
||||
sj_sum += sj;
|
||||
checks.push((e as u32 + 1, check));
|
||||
}
|
||||
if !bad.is_empty() {
|
||||
let why = format!("{} elements fail the certified validity check", bad.len());
|
||||
return Ok(Attempt::Retry(Retry { at: bad, why }));
|
||||
}
|
||||
// The volume the mesh's own boundary encloses equals the elements' total exactly when they
|
||||
// tile it without overlapping.
|
||||
let faces = topo::boundary(&w.hexes);
|
||||
let enclosed: f64 = faces.iter().map(|(_, quad)| face_flux(&quad.map(|v| w.points[v as usize]))).sum();
|
||||
let gate =
|
||||
validity_gate_checks(checks.iter().map(|(l, c)| (*l, c)), &GateOptions { expected_volume: Some(enclosed), volume_tolerance: 1e-9 })?
|
||||
.into_result()?;
|
||||
let volume = gate.volume;
|
||||
let mut new_id = vec![u32::MAX; w.points.len()];
|
||||
let mut nodes = Vec::new();
|
||||
let mut hexes = Vec::with_capacity(w.hexes.len());
|
||||
for hx in &w.hexes {
|
||||
let mut o = [0u32; 8];
|
||||
for (k, &v) in hx.iter().enumerate() {
|
||||
if new_id[v as usize] == u32::MAX {
|
||||
new_id[v as usize] = u32::try_from(nodes.len()).map_err(|_| Error::Limit("too many nodes".into()))?;
|
||||
nodes.push(w.points[v as usize]);
|
||||
}
|
||||
o[k] = new_id[v as usize];
|
||||
}
|
||||
hexes.push(o);
|
||||
}
|
||||
let mut boundary = Vec::new();
|
||||
for ((e, f), quad) in faces {
|
||||
let patch = quad.iter().find_map(|&v| match w.role.get(v as usize) {
|
||||
Some(Constraint::Patch(p)) => Some(*p),
|
||||
_ => None,
|
||||
});
|
||||
let c = [0, 1, 2].map(|k| quad.iter().map(|&v| w.points[v as usize][k]).sum::<f64>() / 4.0);
|
||||
let brep_face = match patch.and_then(|p| g.patches.get(p as usize)) {
|
||||
Some(pt) if pt.faces.len() == 1 => pt.faces[0],
|
||||
_ => {
|
||||
let keep = |t: u32| patch.is_none_or(|p| g.tri_patch.get(t as usize) == Some(&p));
|
||||
g.nearest(c, keep).and_then(|(t, _, _)| g.surface.tri_face().get(t as usize).copied()).unwrap_or(0)
|
||||
}
|
||||
};
|
||||
boundary.push(BoundaryFace { element: e, face: f + 1, brep_face });
|
||||
}
|
||||
let report = OctreeReport {
|
||||
attempts: 0,
|
||||
grid_cells: tree.leaf_count(),
|
||||
core_elements: w.core_count,
|
||||
layer_elements: w.hexes.len() - w.core_count,
|
||||
frame: self.frame,
|
||||
corners: g.corners.len(),
|
||||
curves: g.curves.len(),
|
||||
patches: g.patches.len(),
|
||||
sweeps,
|
||||
min_scaled_jacobian: sj_min,
|
||||
mean_scaled_jacobian: sj_sum / w.hexes.len().max(1) as f64,
|
||||
volume,
|
||||
surface_volume: g.surface.volume(),
|
||||
boundary_distance,
|
||||
retries: Vec::new(),
|
||||
};
|
||||
Ok(Attempt::Done(Box::new(OctreeHexMesh { nodes, hexes, boundary, report })))
|
||||
}
|
||||
}
|
||||
|
||||
/// The core boundary's components and genus against the part's shells; on a mismatch, a retry
|
||||
/// that refines near the whole boundary.
|
||||
fn check_topology(g: &Geometry, quads: &[[u32; 4]]) -> Option<Retry> {
|
||||
let edges = topo::quad_edges(quads);
|
||||
let mut comp: Vec<usize> = (0..quads.len()).collect();
|
||||
for qs in edges.values() {
|
||||
for w in qs.windows(2) {
|
||||
let (a, b) = (topo::find(&mut comp, w[0] as usize), topo::find(&mut comp, w[1] as usize));
|
||||
comp[a] = b;
|
||||
}
|
||||
}
|
||||
// Per component: vertices, edges, faces.
|
||||
let mut stats: BTreeMap<usize, (BTreeSet<u32>, i64, i64)> = BTreeMap::new();
|
||||
for (i, q) in quads.iter().enumerate() {
|
||||
let s = stats.entry(topo::find(&mut comp, i)).or_default();
|
||||
s.2 += 1;
|
||||
s.0.extend(q.iter().copied());
|
||||
}
|
||||
for qs in edges.values() {
|
||||
if let Some(&q0) = qs.first()
|
||||
&& let Some(s) = stats.get_mut(&topo::find(&mut comp, q0 as usize))
|
||||
{
|
||||
s.1 += 1;
|
||||
}
|
||||
}
|
||||
let mut got: Vec<i64> = stats.values().map(|s| (2 - (s.0.len() as i64 - s.1 + s.2)) / 2).collect();
|
||||
let mut want: Vec<i64> = g.surface.shells().iter().map(|s| i64::from(s.genus)).collect();
|
||||
got.sort_unstable();
|
||||
want.sort_unstable();
|
||||
if got == want {
|
||||
return None;
|
||||
}
|
||||
let at = (0..g.surface.triangles().len() as u32).map(|t| g.triangle_centre(t)).collect();
|
||||
Some(Retry { at, why: format!("the core's boundary has components of genus {got:?}, the part {want:?}") })
|
||||
}
|
||||
@@ -0,0 +1,641 @@
|
||||
//! Untangling, fitting and shape optimisation of the hexahedral mesh.
|
||||
//!
|
||||
//! # Energy
|
||||
//! The foldover-free mapping energy of V. Garanzha, I. Kaporin, L. Kudryavtseva, F. Protais,
|
||||
//! N. Ray, D. Sokolov, *Foldover-free maps in 50 lines of code*, ACM TOG 40(4):102, 2021 (paper
|
||||
//! only; the authors' code is AGPL and is not read). Each hexahedron is measured at nine points,
|
||||
//! its eight corners (the trapezoidal rule the paper uses for quads) and its centre, where a
|
||||
//! twisted element folds first: `J` is the sample's three edge directions against the ideal cube
|
||||
//! of the element's target edge `s`, `D = det J`, and the sample costs
|
||||
//!
|
||||
//! `f_ε + λ g_ε`, `f_ε = |J|³_F / χ(D, ε)`, `g_ε = (D² + 1) / χ(D, ε)`,
|
||||
//! `χ(D, ε) = (D + √(ε² + D²)) / 2`,
|
||||
//!
|
||||
//! weighted by `s³`. `f` is the shape (conformal) term, `g` the volume term. The paper's shape
|
||||
//! term is `|J|²/D^{2/3}`; its 3/2 power `|J|³/D`, the classical polyconvex distortion in three
|
||||
//! dimensions (J. M. Ball, 1977), is minimal at the same ideal shapes and needs only square
|
||||
//! roots, so the optimiser gives the same bits on every platform (no `pow` of the system maths
|
||||
//! library). `χ` keeps both finite for inverted samples, and `ε` follows the paper's rule (7),
|
||||
//! `ε = √(10⁻¹² + 0.04 min(0, min D)²)`, so it vanishes once every sample is positive and the
|
||||
//! energy becomes a barrier against inversion.
|
||||
//!
|
||||
//! A vertex on a corner, curve or patch adds `fit · 8 s̄ · dist²` (`s̄` the mean target edge of its
|
||||
//! elements) to its distance from that feature: with `fit` raised stage by stage the boundary is
|
||||
//! drawn onto the geometry while the barrier keeps every element valid (Gao, Shen, Panozzo 2019,
|
||||
//! §3.5). With `fit = 0` constrained vertices are held on their features instead (each step is
|
||||
//! restricted to the tangent and projected back).
|
||||
//!
|
||||
//! # Minimisation
|
||||
//! Nonlinear Gauss–Seidel: each vertex in turn takes a Newton step on its own energy with the
|
||||
//! paper's positive approximation of the Hessian (the second derivatives of `χ` and `det` are
|
||||
//! dropped, §2.3, App. A), accepted only if the energy decreased (halving line search). The paper
|
||||
//! minimises the same energy globally with L-BFGS; the vertex sweep is the local form of the same
|
||||
//! descent. Then a compass search raises the worst corner scaled Jacobians, and another repairs
|
||||
//! elements whose trilinear map the certified check does not prove positive. Sweeps, line-search
|
||||
//! steps, compass steps and rounds are capped.
|
||||
|
||||
use fieldcraft_model::ElementKind;
|
||||
|
||||
use crate::quality::{Validity, element_validity};
|
||||
use crate::vec3::{P, cross, dot, norm, sub};
|
||||
|
||||
use super::geom::dist2;
|
||||
|
||||
use super::geom::Geometry;
|
||||
|
||||
/// How a vertex may move.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub(crate) enum Constraint {
|
||||
Free,
|
||||
/// On a corner of the part (the boundary surface's vertex).
|
||||
Corner(u32),
|
||||
/// On a feature curve.
|
||||
Curve(u32),
|
||||
/// On a patch of the boundary.
|
||||
Patch(u32),
|
||||
}
|
||||
|
||||
/// Model-order corner `c` of a hex8: its three neighbours along x, y, z and the signs that make
|
||||
/// the reference corner frame right-handed.
|
||||
const CORNER_FRAMES: [([usize; 3], [f64; 3]); 8] = [
|
||||
([1, 3, 4], [1.0, 1.0, 1.0]),
|
||||
([0, 2, 5], [-1.0, 1.0, 1.0]),
|
||||
([3, 1, 6], [-1.0, -1.0, 1.0]),
|
||||
([2, 0, 7], [1.0, -1.0, 1.0]),
|
||||
([5, 7, 0], [1.0, 1.0, -1.0]),
|
||||
([4, 6, 1], [-1.0, 1.0, -1.0]),
|
||||
([7, 5, 2], [-1.0, -1.0, -1.0]),
|
||||
([6, 4, 3], [1.0, -1.0, -1.0]),
|
||||
];
|
||||
|
||||
/// Weight of the volume term.
|
||||
const LAMBDA: f64 = 0.1;
|
||||
|
||||
fn chi(d: f64, eps: f64) -> f64 {
|
||||
0.5 * (d + (eps * eps + d * d).sqrt())
|
||||
}
|
||||
|
||||
fn chi_d(d: f64, eps: f64) -> f64 {
|
||||
0.5 * (1.0 + d / (eps * eps + d * d).sqrt())
|
||||
}
|
||||
|
||||
/// The corner Jacobian (columns) of corner `c` of `x` against an ideal cube of edge `s`.
|
||||
fn corner_jacobian(x: &[P; 8], c: usize, s: f64) -> [P; 3] {
|
||||
let (n, sg) = CORNER_FRAMES[c];
|
||||
[0, 1, 2].map(|a| sub(x[n[a]], x[c]).map(|v| v * sg[a] / s))
|
||||
}
|
||||
|
||||
/// The sample points of the energy: the eight corners and the centre. Column `a` of sample `s`'s
|
||||
/// Jacobian is `Σ_k W[s][a][k] x_k` (before dividing by the target edge). The centre keeps
|
||||
/// twisted elements (positive corners, folded inside) out of reach, as the nine-point measure of
|
||||
/// ParaView and of Zhang et al.'s meshers does.
|
||||
const SAMPLES: usize = 9;
|
||||
|
||||
fn sample_weights(s: usize) -> [[f64; 8]; 3] {
|
||||
let mut w = [[0.0; 8]; 3];
|
||||
if s < 8 {
|
||||
let (n, sg) = CORNER_FRAMES[s];
|
||||
for a in 0..3 {
|
||||
w[a][n[a]] += sg[a];
|
||||
w[a][s] -= sg[a];
|
||||
}
|
||||
} else {
|
||||
w = [
|
||||
[-0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25, -0.25],
|
||||
[-0.25, -0.25, 0.25, 0.25, -0.25, -0.25, 0.25, 0.25],
|
||||
[-0.25, -0.25, -0.25, -0.25, 0.25, 0.25, 0.25, 0.25],
|
||||
];
|
||||
}
|
||||
w
|
||||
}
|
||||
|
||||
/// The Jacobian (columns) of sample `smp` of `x` against an ideal cube of edge `s`.
|
||||
fn sample_jacobian(x: &[P; 8], smp: usize, s: f64) -> [P; 3] {
|
||||
let w = sample_weights(smp);
|
||||
[0, 1, 2].map(|a| {
|
||||
let mut c = [0.0; 3];
|
||||
for (k, p) in x.iter().enumerate() {
|
||||
for r in 0..3 {
|
||||
c[r] += w[a][k] * p[r];
|
||||
}
|
||||
}
|
||||
c.map(|v| v / s)
|
||||
})
|
||||
}
|
||||
|
||||
fn det3(j: &[P; 3]) -> f64 {
|
||||
dot(j[0], cross(j[1], j[2]))
|
||||
}
|
||||
|
||||
/// The cofactor columns: `∂ det / ∂ column a`.
|
||||
fn cofactors(j: &[P; 3]) -> [P; 3] {
|
||||
[cross(j[1], j[2]), cross(j[2], j[0]), cross(j[0], j[1])]
|
||||
}
|
||||
|
||||
/// The energy of one corner.
|
||||
fn corner_energy(j: &[P; 3], eps: f64) -> f64 {
|
||||
let d = det3(j);
|
||||
let c = chi(d, eps);
|
||||
let fro: f64 = j.iter().map(|v| dot(*v, *v)).sum();
|
||||
fro * fro.sqrt() / c + LAMBDA * (d * d + 1.0) / c
|
||||
}
|
||||
|
||||
/// The scaled Jacobian of each corner, smallest (`min over corners of det / product of edge
|
||||
/// lengths`), the usual hex quality measure.
|
||||
pub(crate) fn min_scaled_jacobian(x: &[P; 8]) -> f64 {
|
||||
(0..8)
|
||||
.map(|c| {
|
||||
let j = corner_jacobian(x, c, 1.0);
|
||||
let l = norm(j[0]) * norm(j[1]) * norm(j[2]);
|
||||
if l > 0.0 { det3(&j) / l } else { -1.0 }
|
||||
})
|
||||
.fold(f64::INFINITY, f64::min)
|
||||
}
|
||||
|
||||
/// The optimiser's state.
|
||||
pub(crate) struct Optimizer<'a> {
|
||||
pub geom: &'a Geometry<'a>,
|
||||
pub points: Vec<P>,
|
||||
pub hexes: &'a [[u32; 8]],
|
||||
/// Target edge per element.
|
||||
pub size: &'a [f64],
|
||||
pub constraint: Vec<Constraint>,
|
||||
/// Elements around each vertex with the vertex's local position.
|
||||
around: Vec<Vec<(u32, u8)>>,
|
||||
pub eps: f64,
|
||||
/// Weight of the fitting term; 0 holds boundary vertices on their features exactly.
|
||||
pub fit: f64,
|
||||
}
|
||||
|
||||
impl<'a> Optimizer<'a> {
|
||||
pub fn new(geom: &'a Geometry<'a>, points: Vec<P>, hexes: &'a [[u32; 8]], size: &'a [f64], constraint: Vec<Constraint>) -> Optimizer<'a> {
|
||||
let mut around: Vec<Vec<(u32, u8)>> = vec![Vec::new(); points.len()];
|
||||
for (e, h) in hexes.iter().enumerate() {
|
||||
for (k, &v) in h.iter().enumerate() {
|
||||
if let Some(a) = around.get_mut(v as usize) {
|
||||
a.push((e as u32, k as u8));
|
||||
}
|
||||
}
|
||||
}
|
||||
Optimizer { geom, points, hexes, size, constraint, around, eps: 1.0, fit: 0.0 }
|
||||
}
|
||||
|
||||
fn element(&self, e: u32) -> ([P; 8], f64) {
|
||||
let h = self.hexes.get(e as usize).copied().unwrap_or([0; 8]);
|
||||
(h.map(|v| self.points.get(v as usize).copied().unwrap_or([0.0; 3])), self.size.get(e as usize).copied().unwrap_or(1.0))
|
||||
}
|
||||
|
||||
/// The smallest normalised corner determinant over the mesh, and the total energy.
|
||||
pub fn measure(&self) -> (f64, f64) {
|
||||
let (mut m, mut total) = (f64::INFINITY, 0.0);
|
||||
for e in 0..self.hexes.len() as u32 {
|
||||
let (x, s) = self.element(e);
|
||||
for c in 0..SAMPLES {
|
||||
let j = sample_jacobian(&x, c, s);
|
||||
m = m.min(det3(&j));
|
||||
total += s * s * s * corner_energy(&j, self.eps);
|
||||
}
|
||||
}
|
||||
(m, total)
|
||||
}
|
||||
|
||||
/// Where vertex `v`'s feature is closest to `p` (`None` for a free vertex).
|
||||
pub fn target(&self, v: u32, p: P) -> Option<P> {
|
||||
match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) {
|
||||
Constraint::Free => None,
|
||||
Constraint::Corner(c) => Some(self.geom.xyz(c)),
|
||||
Constraint::Curve(c) => self.geom.nearest_on_curve(c, p).map(|x| x.0),
|
||||
Constraint::Patch(pt) => self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)).map(|x| x.1),
|
||||
}
|
||||
}
|
||||
|
||||
/// The fitting weight of vertex `v`: the fit factor times its elements' size (so a distance
|
||||
/// of one element weighs like the shape of its corners), ten times more for a corner.
|
||||
fn fit_weight(&self, v: u32) -> f64 {
|
||||
let a = self.around.get(v as usize).map_or(&[][..], |a| a.as_slice());
|
||||
if a.is_empty() || self.fit == 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
let s = a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(1.0)).sum::<f64>() / a.len() as f64;
|
||||
let corner = if matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_))) { 10.0 } else { 1.0 };
|
||||
8.0 * self.fit * corner * s
|
||||
}
|
||||
|
||||
/// The energy of the corners that involve vertex `v`, with its fitting term.
|
||||
fn local_energy(&self, v: u32) -> f64 {
|
||||
let mut total = 0.0;
|
||||
if self.fit > 0.0 {
|
||||
let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
if let Some(t) = self.target(v, p) {
|
||||
total += self.fit_weight(v) * dist2(p, t);
|
||||
}
|
||||
}
|
||||
for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) {
|
||||
let (x, s) = self.element(e);
|
||||
let (n, _) = CORNER_FRAMES[k as usize];
|
||||
for c in [k as usize, n[0], n[1], n[2], 8] {
|
||||
total += s * s * s * corner_energy(&sample_jacobian(&x, c, s), self.eps);
|
||||
}
|
||||
}
|
||||
total
|
||||
}
|
||||
|
||||
/// The gradient and approximate Hessian of the local energy at vertex `v`.
|
||||
fn derivatives(&self, v: u32) -> (P, [[f64; 3]; 3]) {
|
||||
let mut g = [0.0; 3];
|
||||
let mut h = [[0.0; 3]; 3];
|
||||
for &(e, k) in self.around.get(v as usize).map_or(&[][..], |a| a.as_slice()) {
|
||||
let (x, s) = self.element(e);
|
||||
let k = k as usize;
|
||||
let w = s * s * s;
|
||||
let (nk, _) = CORNER_FRAMES[k];
|
||||
for c in [k, nk[0], nk[1], nk[2], 8] {
|
||||
// ∂ column a / ∂ x_v = coef[a] I.
|
||||
let wt = sample_weights(c);
|
||||
let coef: [f64; 3] = [0, 1, 2].map(|a| wt[a][k] / s);
|
||||
let j = sample_jacobian(&x, c, s);
|
||||
let d = det3(&j);
|
||||
let ch = chi(d, self.eps);
|
||||
let chd = chi_d(d, self.eps);
|
||||
let fro: f64 = j.iter().map(|q| dot(*q, *q)).sum();
|
||||
let cof = cofactors(&j);
|
||||
let root = fro.sqrt();
|
||||
// dE/dJ = A J + B cof(J).
|
||||
let a_coef = 3.0 * root / ch;
|
||||
let f_d = -fro * root * chd / (ch * ch);
|
||||
let g_d = 2.0 * d / ch - (d * d + 1.0) * chd / (ch * ch);
|
||||
let b_coef = f_d + LAMBDA * g_d;
|
||||
let mut dd = [0.0; 3];
|
||||
let mut sum_c2 = 0.0;
|
||||
for a in 0..3 {
|
||||
if coef[a] == 0.0 {
|
||||
continue;
|
||||
}
|
||||
sum_c2 += coef[a] * coef[a];
|
||||
for r in 0..3 {
|
||||
g[r] += w * coef[a] * (a_coef * j[a][r] + b_coef * cof[a][r]);
|
||||
dd[r] += coef[a] * cof[a][r];
|
||||
}
|
||||
}
|
||||
// Positive part of the second derivative along D, χ'' and det'' dropped.
|
||||
let f_dd = 2.0 * fro * root * chd * chd / (ch * ch * ch);
|
||||
let g_dd = 2.0 / ch - 4.0 * d * chd / (ch * ch) + 2.0 * (d * d + 1.0) * chd * chd / (ch * ch * ch);
|
||||
let k_dd = (f_dd + LAMBDA * g_dd).max(0.0);
|
||||
for r in 0..3 {
|
||||
h[r][r] += w * a_coef * sum_c2;
|
||||
for q in 0..3 {
|
||||
h[r][q] += w * k_dd * dd[r] * dd[q];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if self.fit > 0.0 {
|
||||
let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
if let Some(t) = self.target(v, p) {
|
||||
let w = self.fit_weight(v);
|
||||
for r in 0..3 {
|
||||
g[r] += 2.0 * w * (p[r] - t[r]);
|
||||
h[r][r] += 2.0 * w;
|
||||
}
|
||||
}
|
||||
}
|
||||
(g, h)
|
||||
}
|
||||
|
||||
/// Project a point onto vertex `v`'s constraint.
|
||||
fn project(&self, v: u32, p: P) -> P {
|
||||
self.target(v, p).unwrap_or(p)
|
||||
}
|
||||
|
||||
/// The direction a constrained vertex may move along: the step restricted to its tangent.
|
||||
fn restrict(&self, v: u32, d: P) -> P {
|
||||
let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
match self.constraint.get(v as usize).copied().unwrap_or(Constraint::Free) {
|
||||
Constraint::Free => d,
|
||||
Constraint::Corner(_) => [0.0; 3],
|
||||
Constraint::Curve(c) => {
|
||||
let Some((_, i, _)) = self.geom.nearest_on_curve(c, p) else { return [0.0; 3] };
|
||||
let Some(cv) = self.geom.curves.get(c as usize) else { return [0.0; 3] };
|
||||
let (a, b) = (cv.vertices.get(i).copied().unwrap_or(0), cv.vertices.get(i + 1).copied().unwrap_or(0));
|
||||
let t = sub(self.geom.xyz(b), self.geom.xyz(a));
|
||||
let l2 = dot(t, t);
|
||||
if l2 > 0.0 { t.map(|x| x * dot(d, t) / l2) } else { [0.0; 3] }
|
||||
}
|
||||
Constraint::Patch(pt) => {
|
||||
let Some((t, _, _)) = self.geom.nearest(p, |t| self.geom.tri_patch.get(t as usize) == Some(&pt)) else { return [0.0; 3] };
|
||||
let n = self.geom.normals.get(t as usize).copied().unwrap_or([0.0; 3]);
|
||||
let dn = dot(d, n);
|
||||
[0, 1, 2].map(|k| d[k] - dn * n[k])
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// One Newton step of vertex `v`; whether it moved.
|
||||
fn relax(&mut self, v: u32) -> bool {
|
||||
let soft = self.fit > 0.0;
|
||||
if (!soft && matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_))))
|
||||
|| self.around.get(v as usize).is_none_or(|a| a.is_empty())
|
||||
{
|
||||
return false;
|
||||
}
|
||||
let (g, h) = self.derivatives(v);
|
||||
let Some(step) = solve3(h, g.map(|x| -x)) else { return false };
|
||||
let step = if soft { step } else { self.restrict(v, step) };
|
||||
if !step.iter().all(|x| x.is_finite()) || norm(step) == 0.0 {
|
||||
return false;
|
||||
}
|
||||
let old = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
let e0 = self.local_energy(v);
|
||||
let mut tau = 1.0;
|
||||
for _ in 0..12 {
|
||||
let moved = [0, 1, 2].map(|k| old[k] + tau * step[k]);
|
||||
let trial = if soft { moved } else { self.project(v, moved) };
|
||||
if let Some(p) = self.points.get_mut(v as usize) {
|
||||
*p = trial;
|
||||
}
|
||||
let e1 = self.local_energy(v);
|
||||
if e1 < e0 && e1.is_finite() {
|
||||
return true;
|
||||
}
|
||||
tau *= 0.5;
|
||||
}
|
||||
if let Some(p) = self.points.get_mut(v as usize) {
|
||||
*p = old;
|
||||
}
|
||||
false
|
||||
}
|
||||
|
||||
/// Set the regularisation from the smallest corner determinant (Garanzha et al. 2021: ε
|
||||
/// shrinks as the mesh untangles); that determinant and the energy at the new ε.
|
||||
fn update_eps(&mut self) -> (f64, f64) {
|
||||
let (min_d, energy) = self.measure();
|
||||
let old = self.eps;
|
||||
self.eps = (1e-12 + 0.04 * min_d.min(0.0).powi(2)).sqrt();
|
||||
(min_d, if self.eps == old { energy } else { self.measure().1 })
|
||||
}
|
||||
|
||||
/// Sweep all vertices until every corner is positive and the energy settles, or `max_sweeps`;
|
||||
/// the sweeps run.
|
||||
pub fn run(&mut self, max_sweeps: usize, tolerance: f64) -> usize {
|
||||
let n = self.points.len() as u32;
|
||||
let (_, mut energy) = self.update_eps();
|
||||
let mut sweeps = 0;
|
||||
for sweep in 0..max_sweeps {
|
||||
let mut moved = 0usize;
|
||||
for v in 0..n {
|
||||
if self.relax(v) {
|
||||
moved += 1;
|
||||
}
|
||||
}
|
||||
let old_eps = self.eps;
|
||||
let (min_d, e) = self.update_eps();
|
||||
sweeps = sweep + 1;
|
||||
let settled = old_eps == self.eps && (energy - e).abs() <= tolerance * energy.abs();
|
||||
energy = e;
|
||||
if moved == 0 || (min_d > 0.0 && settled) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
sweeps
|
||||
}
|
||||
}
|
||||
|
||||
impl Optimizer<'_> {
|
||||
/// The smallest corner scaled Jacobian of the elements around vertex `v`.
|
||||
fn local_quality(&self, v: u32) -> f64 {
|
||||
self.around
|
||||
.get(v as usize)
|
||||
.map_or(f64::INFINITY, |a| a.iter().map(|&(e, _)| min_scaled_jacobian(&self.element(e).0)).fold(f64::INFINITY, f64::min))
|
||||
}
|
||||
|
||||
/// Move every constrained vertex exactly onto its feature where its elements keep a corner
|
||||
/// scaled Jacobian of at least `floor`; the others stay where the fit left them. Returns
|
||||
/// how many stayed off and the largest distance of a constrained vertex from its feature.
|
||||
pub fn snap(&mut self, floor: f64) -> (usize, f64) {
|
||||
let (mut off, mut worst) = (0usize, 0.0f64);
|
||||
for v in 0..self.points.len() as u32 {
|
||||
let p = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
let Some(t) = self.target(v, p) else { continue };
|
||||
if p == t {
|
||||
continue;
|
||||
}
|
||||
if let Some(x) = self.points.get_mut(v as usize) {
|
||||
*x = t;
|
||||
}
|
||||
if self.local_quality(v) < floor {
|
||||
if let Some(x) = self.points.get_mut(v as usize) {
|
||||
*x = p;
|
||||
}
|
||||
off += 1;
|
||||
worst = worst.max(dist2(p, t).sqrt());
|
||||
}
|
||||
}
|
||||
(off, worst)
|
||||
}
|
||||
|
||||
/// Raise the worst elements: every vertex of an element whose corner scaled Jacobian is
|
||||
/// below `threshold` moves, within its constraint, by a compass search that maximises the
|
||||
/// smallest scaled Jacobian around it (steps shrink by halves; at most `passes` rounds over
|
||||
/// the worst elements, 12 step sizes per vertex).
|
||||
pub fn improve_worst(&mut self, threshold: f64, passes: usize) {
|
||||
for _ in 0..passes {
|
||||
let mut verts: Vec<u32> = Vec::new();
|
||||
for (e, h) in self.hexes.iter().enumerate() {
|
||||
if min_scaled_jacobian(&self.element(e as u32).0) < threshold {
|
||||
verts.extend(h.iter().copied());
|
||||
}
|
||||
}
|
||||
verts.sort_unstable();
|
||||
verts.dedup();
|
||||
if verts.is_empty() {
|
||||
break;
|
||||
}
|
||||
let mut moved = false;
|
||||
for &v in &verts {
|
||||
if matches!(self.constraint.get(v as usize), Some(Constraint::Free)) {
|
||||
moved |= self.compass(v, &|o: &Self, v| o.local_quality(v));
|
||||
}
|
||||
}
|
||||
if !moved {
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Compass search on vertex `v` (within its constraint) for the largest `quality`; steps
|
||||
/// start at 0.15 of its largest element size and halve when no direction improves, 12 times
|
||||
/// at most. Whether it moved.
|
||||
fn compass(&mut self, v: u32, quality: &dyn Fn(&Self, u32) -> f64) -> bool {
|
||||
let p0 = self.points.get(v as usize).copied().unwrap_or([0.0; 3]);
|
||||
let reach = self
|
||||
.around
|
||||
.get(v as usize)
|
||||
.map_or(0.0, |a| a.iter().map(|&(e, _)| self.size.get(e as usize).copied().unwrap_or(0.0)).fold(0.0, f64::max));
|
||||
let mut step = 0.15 * reach;
|
||||
let mut best = quality(self, v);
|
||||
let mut cur = p0;
|
||||
let mut moved = false;
|
||||
for _ in 0..12 {
|
||||
if step <= 1e-9 * reach.max(1e-300) {
|
||||
break;
|
||||
}
|
||||
let mut improved = false;
|
||||
for d in [[1.0, 0.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, 1.0], [0.0, 0.0, -1.0]] {
|
||||
let dir = self.restrict(v, d.map(|x| x * step));
|
||||
if norm(dir) == 0.0 {
|
||||
continue;
|
||||
}
|
||||
let trial = self.project(v, [0, 1, 2].map(|k| cur[k] + dir[k]));
|
||||
if let Some(p) = self.points.get_mut(v as usize) {
|
||||
*p = trial;
|
||||
}
|
||||
let q = quality(self, v);
|
||||
if q > best + 1e-12 {
|
||||
best = q;
|
||||
cur = trial;
|
||||
improved = true;
|
||||
moved = true;
|
||||
} else if let Some(p) = self.points.get_mut(v as usize) {
|
||||
*p = cur;
|
||||
}
|
||||
}
|
||||
if !improved {
|
||||
step *= 0.5;
|
||||
}
|
||||
}
|
||||
if let Some(p) = self.points.get_mut(v as usize) {
|
||||
*p = cur;
|
||||
}
|
||||
moved
|
||||
}
|
||||
|
||||
/// The smallest certified Jacobian ratio of the elements around `v` (negative where the
|
||||
/// trilinear map folds somewhere, even with positive corners).
|
||||
fn certified_quality(&self, v: u32) -> f64 {
|
||||
self.around.get(v as usize).map_or(f64::INFINITY, |a| {
|
||||
a.iter()
|
||||
.map(|&(e, _)| {
|
||||
let x = self.element(e).0;
|
||||
if min_scaled_jacobian(&x) <= 0.0 {
|
||||
return -1e9;
|
||||
}
|
||||
element_validity(ElementKind::Hex8, &x).map_or(-1.0, |c| match c.validity {
|
||||
Validity::Valid => c.ratio.max(1e-300),
|
||||
_ => c.ratio.min(0.0) - 1e-12,
|
||||
})
|
||||
})
|
||||
.fold(f64::INFINITY, f64::min)
|
||||
})
|
||||
}
|
||||
|
||||
/// Repair elements whose trilinear map is not proven positive although their corners are:
|
||||
/// the vertices of each such element (free ones, and boundary ones along their features)
|
||||
/// move by compass search to raise the smallest certified bound around them; `rounds` rounds
|
||||
/// at most. Returns the elements still not proven valid.
|
||||
pub fn repair_certified(&mut self, rounds: usize) -> Vec<u32> {
|
||||
let check = |o: &Self, e: u32| element_validity(ElementKind::Hex8, &o.element(e).0).is_ok_and(|c| c.validity == Validity::Valid);
|
||||
let mut bad: Vec<u32> = (0..self.hexes.len() as u32).filter(|&e| !check(self, e)).collect();
|
||||
for _ in 0..rounds {
|
||||
if bad.is_empty() {
|
||||
break;
|
||||
}
|
||||
let mut verts: Vec<u32> = bad.iter().flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8])).collect();
|
||||
verts.sort_unstable();
|
||||
verts.dedup();
|
||||
for v in verts {
|
||||
if !matches!(self.constraint.get(v as usize), Some(Constraint::Corner(_)) | None) {
|
||||
self.compass(v, &|o: &Self, v| o.certified_quality(v));
|
||||
}
|
||||
}
|
||||
let mut around: Vec<u32> = bad
|
||||
.iter()
|
||||
.flat_map(|&e| self.hexes.get(e as usize).copied().unwrap_or([0; 8]))
|
||||
.flat_map(|v| self.around.get(v as usize).map_or(Vec::new(), |a| a.iter().map(|x| x.0).collect()))
|
||||
.collect();
|
||||
around.sort_unstable();
|
||||
around.dedup();
|
||||
bad = around.into_iter().filter(|&e| !check(self, e)).collect();
|
||||
}
|
||||
bad
|
||||
}
|
||||
}
|
||||
|
||||
/// Solve the symmetric 3 × 3 system `h x = b` (Cramer), `None` when singular.
|
||||
fn solve3(h: [[f64; 3]; 3], b: P) -> Option<P> {
|
||||
let d = dot(h[0], cross(h[1], h[2]));
|
||||
if !(d.abs() > 1e-300) || !d.is_finite() {
|
||||
return None;
|
||||
}
|
||||
let col = |k: usize| [h[0][k], h[1][k], h[2][k]];
|
||||
let (c0, c1, c2) = (col(0), col(1), col(2));
|
||||
Some([dot(b, cross(c1, c2)) / d, dot(c0, cross(b, c2)) / d, dot(c0, cross(c1, b)) / d])
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::octree::geom::tests::box_surface;
|
||||
|
||||
#[test]
|
||||
fn corner_frames_are_right_handed_on_the_unit_cube() {
|
||||
let cube: [P; 8] = [[0., 0., 0.], [1., 0., 0.], [1., 1., 0.], [0., 1., 0.], [0., 0., 1.], [1., 0., 1.], [1., 1., 1.], [0., 1., 1.]];
|
||||
for c in 0..8 {
|
||||
let j = corner_jacobian(&cube, c, 1.0);
|
||||
assert!((det3(&j) - 1.0).abs() < 1e-15, "corner {c}");
|
||||
}
|
||||
assert!((min_scaled_jacobian(&cube) - 1.0).abs() < 1e-15);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_tangled_interior_vertex_is_untangled() {
|
||||
// A 2 × 2 × 2 block of cubes with the centre vertex pushed through a face; the boundary
|
||||
// fixed.
|
||||
let s = box_surface([0.0; 3], [2.0; 3]);
|
||||
let g = Geometry::new(&s, 30.0, 160.0).unwrap();
|
||||
let id = |x: u32, y: u32, z: u32| x + 3 * (y + 3 * z);
|
||||
let mut p = Vec::new();
|
||||
for z in 0..3 {
|
||||
for y in 0..3 {
|
||||
for x in 0..3 {
|
||||
p.push([x as f64, y as f64, z as f64]);
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut h = Vec::new();
|
||||
for z in 0..2 {
|
||||
for y in 0..2 {
|
||||
for x in 0..2 {
|
||||
h.push([
|
||||
id(x, y, z),
|
||||
id(x + 1, y, z),
|
||||
id(x + 1, y + 1, z),
|
||||
id(x, y + 1, z),
|
||||
id(x, y, z + 1),
|
||||
id(x + 1, y, z + 1),
|
||||
id(x + 1, y + 1, z + 1),
|
||||
id(x, y + 1, z + 1),
|
||||
]);
|
||||
}
|
||||
}
|
||||
}
|
||||
let centre = id(1, 1, 1) as usize;
|
||||
p[centre] = [1.7, 1.9, 0.2];
|
||||
let size = vec![1.0; h.len()];
|
||||
// Boundary vertices slide on the box's faces (hard constraints, patch of the nearest
|
||||
// face), the centre is free.
|
||||
let cons: Vec<Constraint> = (0..p.len())
|
||||
.map(|i| {
|
||||
if i == centre {
|
||||
return Constraint::Free;
|
||||
}
|
||||
let t = g.nearest(p[i], |_| true).unwrap().0;
|
||||
Constraint::Patch(g.tri_patch[t as usize])
|
||||
})
|
||||
.collect();
|
||||
let mut o = Optimizer::new(&g, p, &h, &size, cons);
|
||||
assert!(o.measure().0 < 0.0);
|
||||
let st = o.run(200, 1e-9);
|
||||
assert!(o.measure().0 > 0.9, "{st:?}");
|
||||
assert!(sub(o.points[centre], [1.0; 3]).iter().all(|c| c.abs() < 1e-3), "{:?}", o.points[centre]);
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,374 @@
|
||||
//! The 3-refinement transition templates: how a leaf of the balanced tree is cut into
|
||||
//! hexahedra so that it conforms to finer neighbours.
|
||||
//!
|
||||
//! Every vertex of a leaf is *marked* when a finer leaf touches it. Under strong balance the
|
||||
//! marked set of a leaf is one of five patterns (after the closure in [`super::grid`] has split
|
||||
//! every leaf that had another one), each with one template:
|
||||
//!
|
||||
//! | marks | template | hexahedra |
|
||||
//! |---|---|---|
|
||||
//! | none | the leaf itself | 1 |
|
||||
//! | one corner | corner cube and three blocks fanning to the far corner | 4 |
|
||||
//! | the two ends of an edge | three cubes along the edge and two four-hex slabs | 11 |
|
||||
//! | the four corners of a face | nine cubes on the face, nine above them, four on top | 22 |
|
||||
//! | all eight | the 3 × 3 × 3 split | 27 |
|
||||
//!
|
||||
//! A face of a leaf is cut according to its own marks only (none, one corner, one edge, all
|
||||
//! four), and an edge into thirds at its marked ends, so two leaves sharing a face cut it the same
|
||||
//! way: the mesh is conforming by construction. This is the node-marking form of Schneiders'
|
||||
//! 3-refinement (R. Schneiders, *Refining quadrilateral and hexahedral element meshes*, NUMGRID
|
||||
//! 1996; Y. Ito, A. M. Shih, B. K. Soni, *Octree-based reasonable-quality hexahedral mesh
|
||||
//! generation using a new set of refinement templates*, IJNME 77(13):1809–1833, 2009). The
|
||||
//! templates are our own construction from those papers' description: the face template's top
|
||||
//! block is the 3 × 3 → 1 transition traced from its dual sheets, the edge template is two
|
||||
//! copies of the planar 3 → 1 transition swept across the leaf, and the tests prove every
|
||||
//! template's faces, volume and validity.
|
||||
//!
|
||||
//! Coordinates are integers in sixths of the leaf's edge, so every vertex shared with a neighbour
|
||||
//! lies on the global integer lattice of [`super::grid`] and coincident vertices merge exactly.
|
||||
|
||||
/// A point of a leaf in sixths of its edge, `0..=6` per axis.
|
||||
pub(crate) type Local = [i64; 3];
|
||||
/// A hexahedron in the model's node order (corners 1-2-3-4 then 5-6-7-8 above them).
|
||||
pub(crate) type Hex = [Local; 8];
|
||||
|
||||
/// The marked-corner pattern of a leaf, up to symmetry.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub(crate) enum Pattern {
|
||||
None,
|
||||
Corner,
|
||||
Edge,
|
||||
Face,
|
||||
Full,
|
||||
}
|
||||
|
||||
/// A symmetry of the cube: coordinate `k` of the image is coordinate `perm[k]` of the point,
|
||||
/// reflected (`6 − v`) when `flip[k]`.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub(crate) struct Sym {
|
||||
perm: [usize; 3],
|
||||
flip: [bool; 3],
|
||||
}
|
||||
|
||||
impl Sym {
|
||||
pub(crate) fn apply(&self, p: Local) -> Local {
|
||||
let mut q = [0; 3];
|
||||
for (k, out) in q.iter_mut().enumerate() {
|
||||
let v = p[self.perm[k]];
|
||||
*out = if self.flip[k] { 6 - v } else { v };
|
||||
}
|
||||
q
|
||||
}
|
||||
}
|
||||
|
||||
/// The 48 symmetries in a fixed order.
|
||||
pub(crate) fn symmetries() -> Vec<Sym> {
|
||||
let perms = [[0, 1, 2], [0, 2, 1], [1, 0, 2], [1, 2, 0], [2, 0, 1], [2, 1, 0]];
|
||||
let mut out = Vec::with_capacity(48);
|
||||
for perm in perms {
|
||||
for bits in 0..8u8 {
|
||||
out.push(Sym { perm, flip: [bits & 1 != 0, bits & 2 != 0, bits & 4 != 0] });
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// The leaf corner `c` (bit 0: x, bit 1: y, bit 2: z) in sixths.
|
||||
pub(crate) fn corner(c: usize) -> Local {
|
||||
[(c & 1) as i64 * 6, ((c >> 1) & 1) as i64 * 6, ((c >> 2) & 1) as i64 * 6]
|
||||
}
|
||||
|
||||
/// The corner index of a corner point.
|
||||
fn corner_index(p: Local) -> usize {
|
||||
usize::from(p[0] == 6) | (usize::from(p[1] == 6) << 1) | (usize::from(p[2] == 6) << 2)
|
||||
}
|
||||
|
||||
/// The image of a corner mask under a symmetry.
|
||||
fn map_mask(sym: &Sym, mask: u8) -> u8 {
|
||||
(0..8).filter(|c| mask & (1 << c) != 0).fold(0u8, |m, c| m | (1 << corner_index(sym.apply(corner(c)))))
|
||||
}
|
||||
|
||||
/// The pattern of a corner mask, or `None` when no template fits it.
|
||||
pub(crate) fn classify(mask: u8) -> Option<Pattern> {
|
||||
let corners: Vec<Local> = (0..8).filter(|c| mask & (1 << c) != 0).map(corner).collect();
|
||||
let shared = |n: usize| (0..3).filter(|&k| corners.iter().all(|p| p[k] == corners[0][k])).count() == n;
|
||||
match corners.len() {
|
||||
0 => Some(Pattern::None),
|
||||
1 => Some(Pattern::Corner),
|
||||
2 if shared(2) => Some(Pattern::Edge),
|
||||
4 if shared(1) => Some(Pattern::Face),
|
||||
8 => Some(Pattern::Full),
|
||||
_ => None,
|
||||
}
|
||||
}
|
||||
|
||||
/// The canonical mask of a pattern: corner 0; edge 0–1 (along x); face z = 0.
|
||||
fn canonical_mask(p: Pattern) -> u8 {
|
||||
match p {
|
||||
Pattern::None => 0,
|
||||
Pattern::Corner => 0b1,
|
||||
Pattern::Edge => 0b11,
|
||||
Pattern::Face => 0b1111,
|
||||
Pattern::Full => 0xff,
|
||||
}
|
||||
}
|
||||
|
||||
/// The first symmetry (in [`symmetries`] order) taking the pattern's canonical mask to `mask`.
|
||||
pub(crate) fn orientation(syms: &[Sym], pattern: Pattern, mask: u8) -> Option<Sym> {
|
||||
let c = canonical_mask(pattern);
|
||||
syms.iter().find(|s| map_mask(s, c) == mask).copied()
|
||||
}
|
||||
|
||||
/// Thirds to sixths.
|
||||
fn t(x: i64, y: i64, z: i64) -> Local {
|
||||
[2 * x, 2 * y, 2 * z]
|
||||
}
|
||||
|
||||
/// A thirds-aligned cube with lower corner `(x, y, z)` (in thirds) and edge `e` thirds.
|
||||
fn cube(x: i64, y: i64, z: i64, e: i64) -> Hex {
|
||||
let (a, b) = (t(x, y, z), t(x + e, y + e, z + e));
|
||||
[
|
||||
[a[0], a[1], a[2]],
|
||||
[b[0], a[1], a[2]],
|
||||
[b[0], b[1], a[2]],
|
||||
[a[0], b[1], a[2]],
|
||||
[a[0], a[1], b[2]],
|
||||
[b[0], a[1], b[2]],
|
||||
[b[0], b[1], b[2]],
|
||||
[a[0], b[1], b[2]],
|
||||
]
|
||||
}
|
||||
|
||||
/// The planar 3 → 1 transition in `(x, s)`: x in thirds, `s` in halves of the transition height,
|
||||
/// three segments at `s = 0` and one at `s = 2`.
|
||||
const THREE_TO_ONE: [[(i64, i64); 4]; 4] =
|
||||
[[(0, 0), (1, 0), (1, 1), (0, 2)], [(1, 0), (2, 0), (2, 1), (1, 1)], [(2, 0), (3, 0), (3, 2), (2, 1)], [(1, 1), (2, 1), (3, 2), (0, 2)]];
|
||||
|
||||
/// The canonical template of a pattern, in sixths of the leaf.
|
||||
pub(crate) fn template(p: Pattern) -> Vec<Hex> {
|
||||
match p {
|
||||
Pattern::None => vec![cube(0, 0, 0, 3)],
|
||||
Pattern::Full => (0..27).map(|k| cube(k % 3, (k / 3) % 3, k / 9, 1)).collect(),
|
||||
Pattern::Corner => vec![
|
||||
cube(0, 0, 0, 1),
|
||||
[t(1, 0, 0), t(3, 0, 0), t(3, 3, 0), t(1, 1, 0), t(1, 0, 1), t(3, 0, 3), t(3, 3, 3), t(1, 1, 1)],
|
||||
[t(0, 1, 0), t(1, 1, 0), t(3, 3, 0), t(0, 3, 0), t(0, 1, 1), t(1, 1, 1), t(3, 3, 3), t(0, 3, 3)],
|
||||
[t(0, 0, 1), t(1, 0, 1), t(1, 1, 1), t(0, 1, 1), t(0, 0, 3), t(3, 0, 3), t(3, 3, 3), t(0, 3, 3)],
|
||||
],
|
||||
Pattern::Edge => {
|
||||
let mut out: Vec<Hex> = (0..3).map(|i| cube(i, 0, 0, 1)).collect();
|
||||
// Two slabs over the edge column's faces y = 1 and z = 1, each the 3 → 1 transition
|
||||
// swept from the leaf face to the diagonal plane y = z.
|
||||
for q in THREE_TO_ONE {
|
||||
let b = |(x, s): (i64, i64), w: i64| t(x, 1 + s, w * (1 + s));
|
||||
let c = |(x, s): (i64, i64), w: i64| t(x, w * (1 + s), 1 + s);
|
||||
out.push([b(q[0], 0), b(q[1], 0), b(q[2], 0), b(q[3], 0), b(q[0], 1), b(q[1], 1), b(q[2], 1), b(q[3], 1)]);
|
||||
out.push([c(q[0], 0), c(q[1], 0), c(q[2], 0), c(q[3], 0), c(q[0], 1), c(q[1], 1), c(q[2], 1), c(q[3], 1)]);
|
||||
}
|
||||
out
|
||||
}
|
||||
Pattern::Face => {
|
||||
let mut out: Vec<Hex> = (0..9).map(|k| cube(k % 3, k / 3, 0, 1)).collect();
|
||||
// The layer above: its top is the corners raised to the leaf's top, the side nodes of
|
||||
// the 3 → 1 transitions at two thirds, the inside nodes at half height.
|
||||
let top = |a: i64, b: i64| -> Local {
|
||||
let edge = |v: i64| v == 0 || v == 3;
|
||||
match (edge(a), edge(b)) {
|
||||
(true, true) => [2 * a, 2 * b, 6],
|
||||
(false, false) => [2 * a, 2 * b, 3],
|
||||
_ => [2 * a, 2 * b, 4],
|
||||
}
|
||||
};
|
||||
for k in 0..9 {
|
||||
let (i, j) = (k % 3, k / 3);
|
||||
out.push([
|
||||
t(i, j, 1),
|
||||
t(i + 1, j, 1),
|
||||
t(i + 1, j + 1, 1),
|
||||
t(i, j + 1, 1),
|
||||
top(i, j),
|
||||
top(i + 1, j),
|
||||
top(i + 1, j + 1),
|
||||
top(i, j + 1),
|
||||
]);
|
||||
}
|
||||
// The four on top: the y-arch (low), the x-arch (high) and the two side blocks.
|
||||
out.push([[2, 2, 3], [4, 2, 3], [4, 4, 3], [2, 4, 3], [2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4]]);
|
||||
out.push([[2, 0, 4], [4, 0, 4], [4, 6, 4], [2, 6, 4], [0, 0, 6], [6, 0, 6], [6, 6, 6], [0, 6, 6]]);
|
||||
out.push([[2, 2, 3], [2, 4, 3], [2, 6, 4], [2, 0, 4], [0, 2, 4], [0, 4, 4], [0, 6, 6], [0, 0, 6]]);
|
||||
out.push([[4, 2, 3], [4, 4, 3], [4, 6, 4], [4, 0, 4], [6, 2, 4], [6, 4, 4], [6, 6, 6], [6, 0, 6]]);
|
||||
out
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Six times the signed volume of the corner tetrahedron at node 0 (exact on the lattice).
|
||||
pub(crate) fn corner0_orientation(h: &[[i64; 3]; 8]) -> i128 {
|
||||
let e = |j: usize| [0, 1, 2].map(|k| i128::from(h[j][k] - h[0][k]));
|
||||
let (a, b, c) = (e(1), e(3), e(4));
|
||||
a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0]) + a[2] * (b[0] * c[1] - b[1] * c[0])
|
||||
}
|
||||
|
||||
/// The hexahedron with its orientation reversed (the same cell, mirrored node order).
|
||||
pub(crate) fn reversed<T: Copy>(h: [T; 8]) -> [T; 8] {
|
||||
[h[0], h[3], h[2], h[1], h[4], h[7], h[6], h[5]]
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use std::collections::HashMap;
|
||||
|
||||
use super::*;
|
||||
use crate::quality::{Validity, element_validity};
|
||||
use fieldcraft_model::ElementKind;
|
||||
|
||||
/// The face of the leaf a quad lies on (axis, value), if any.
|
||||
fn on_leaf_face(q: &[Local; 4]) -> Option<(usize, i64)> {
|
||||
(0..3).find_map(|k| [0, 6].into_iter().find(|&v| q.iter().all(|p| p[k] == v)).map(|v| (k, v)))
|
||||
}
|
||||
|
||||
fn sorted(mut q: [Local; 4]) -> [Local; 4] {
|
||||
q.sort();
|
||||
q
|
||||
}
|
||||
|
||||
/// The quads a leaf face must carry for its four corner marks, as sorted vertex sets, on the
|
||||
/// face `(axis, value)`.
|
||||
fn face_template(axis: usize, value: i64, marks: [bool; 4]) -> Vec<[Local; 4]> {
|
||||
// Face coordinates (u, v) are the other two axes in increasing order; corners in order
|
||||
// (0,0), (6,0), (0,6), (6,6).
|
||||
let (ua, va) = match axis {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
let lift = |u: i64, v: i64| {
|
||||
let mut p = [0; 3];
|
||||
p[axis] = value;
|
||||
p[ua] = u;
|
||||
p[va] = v;
|
||||
p
|
||||
};
|
||||
// Build in a canonical frame, then map the 2D symmetry that takes the canonical marks to
|
||||
// these marks.
|
||||
let count = marks.iter().filter(|m| **m).count();
|
||||
let canon: Vec<[(i64, i64); 4]> = match count {
|
||||
0 => vec![[(0, 0), (6, 0), (6, 6), (0, 6)]],
|
||||
1 => vec![[(0, 0), (2, 0), (2, 2), (0, 2)], [(2, 0), (6, 0), (6, 6), (2, 2)], [(0, 2), (2, 2), (6, 6), (0, 6)]],
|
||||
2 => {
|
||||
let mut v: Vec<[(i64, i64); 4]> = (0..3).map(|i| [(2 * i, 0), (2 * i + 2, 0), (2 * i + 2, 2), (2 * i, 2)]).collect();
|
||||
v.extend([
|
||||
[(0, 2), (2, 2), (2, 4), (0, 6)],
|
||||
[(2, 2), (4, 2), (4, 4), (2, 4)],
|
||||
[(4, 2), (6, 2), (6, 6), (4, 4)],
|
||||
[(2, 4), (4, 4), (6, 6), (0, 6)],
|
||||
]);
|
||||
v
|
||||
}
|
||||
_ => (0..9)
|
||||
.map(|k| {
|
||||
[(2 * (k % 3), 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3)), (2 * (k % 3) + 2, 2 * (k / 3) + 2), (2 * (k % 3), 2 * (k / 3) + 2)]
|
||||
})
|
||||
.collect(),
|
||||
};
|
||||
let canon_marks: [bool; 4] = match count {
|
||||
0 => [false; 4],
|
||||
1 => [true, false, false, false],
|
||||
2 => [true, true, false, false],
|
||||
_ => [true; 4],
|
||||
};
|
||||
let maps: [fn(i64, i64) -> (i64, i64); 8] = [
|
||||
|u, v| (u, v),
|
||||
|u, v| (6 - u, v),
|
||||
|u, v| (u, 6 - v),
|
||||
|u, v| (6 - u, 6 - v),
|
||||
|u, v| (v, u),
|
||||
|u, v| (6 - v, u),
|
||||
|u, v| (v, 6 - u),
|
||||
|u, v| (6 - v, 6 - u),
|
||||
];
|
||||
let corners = [(0, 0), (6, 0), (0, 6), (6, 6)];
|
||||
let idx = |c: (i64, i64)| corners.iter().position(|&x| x == c).unwrap();
|
||||
let m =
|
||||
maps.iter().find(|m| (0..4).all(|c| canon_marks[c] == marks[idx(m(corners[c].0, corners[c].1))])).expect("a face symmetry for the marks");
|
||||
canon.iter().map(|q| sorted(q.map(|(u, v)| m(u, v)).map(|(u, v)| lift(u, v)))).collect()
|
||||
}
|
||||
|
||||
fn face_corner_marks(axis: usize, value: i64, mask: u8) -> [bool; 4] {
|
||||
let (ua, va) = match axis {
|
||||
0 => (1, 2),
|
||||
1 => (0, 2),
|
||||
_ => (0, 1),
|
||||
};
|
||||
[(0, 0), (6, 0), (0, 6), (6, 6)].map(|(u, v)| {
|
||||
let mut p = [0; 3];
|
||||
p[axis] = value;
|
||||
p[ua] = u;
|
||||
p[va] = v;
|
||||
mask & (1 << corner_index(p)) != 0
|
||||
})
|
||||
}
|
||||
|
||||
/// Every template under every symmetry: valid elements, the leaf's volume, closed (each inner
|
||||
/// face shared by two elements) and, on each leaf face, exactly the face template of that
|
||||
/// face's marks.
|
||||
#[test]
|
||||
fn templates_are_valid_closed_and_conform_to_their_face_marks() {
|
||||
let syms = symmetries();
|
||||
for pattern in [Pattern::None, Pattern::Corner, Pattern::Edge, Pattern::Face, Pattern::Full] {
|
||||
for sym in &syms {
|
||||
let mask = map_mask(sym, canonical_mask(pattern));
|
||||
assert_eq!(classify(mask), Some(pattern));
|
||||
let hexes: Vec<Hex> = template(pattern)
|
||||
.into_iter()
|
||||
.map(|h| {
|
||||
let h = h.map(|p| sym.apply(p));
|
||||
if corner0_orientation(&h) < 0 { reversed(h) } else { h }
|
||||
})
|
||||
.collect();
|
||||
let mut volume = 0.0;
|
||||
let mut faces: HashMap<[Local; 4], usize> = HashMap::new();
|
||||
for h in &hexes {
|
||||
let xyz: Vec<[f64; 3]> = h.iter().map(|p| p.map(|v| v as f64)).collect();
|
||||
let check = element_validity(ElementKind::Hex8, &xyz).unwrap();
|
||||
assert_eq!(check.validity, Validity::Valid, "{pattern:?} {sym:?} {h:?}");
|
||||
volume += check.volume;
|
||||
for f in ElementKind::Hex8.faces() {
|
||||
*faces.entry(sorted([h[f[0]], h[f[1]], h[f[2]], h[f[3]]])).or_default() += 1;
|
||||
}
|
||||
}
|
||||
assert!((volume - 216.0).abs() < 1e-9, "{pattern:?}: volume {volume}");
|
||||
let mut outer: HashMap<(usize, i64), Vec<[Local; 4]>> = HashMap::new();
|
||||
for (q, n) in &faces {
|
||||
match on_leaf_face(q) {
|
||||
Some(f) => {
|
||||
assert_eq!(*n, 1);
|
||||
outer.entry(f).or_default().push(*q);
|
||||
}
|
||||
None => assert_eq!(*n, 2, "{pattern:?} {sym:?}: inner face {q:?} used {n} times"),
|
||||
}
|
||||
}
|
||||
for axis in 0..3 {
|
||||
for value in [0, 6] {
|
||||
let mut got = outer.remove(&(axis, value)).unwrap_or_default();
|
||||
let mut want = face_template(axis, value, face_corner_marks(axis, value, mask));
|
||||
got.sort();
|
||||
want.sort();
|
||||
assert_eq!(got, want, "{pattern:?} {sym:?} face {axis}={value}");
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn patterns_without_a_template_are_recognised() {
|
||||
assert_eq!(classify(0b1001), None, "diagonal of a face");
|
||||
assert_eq!(classify(0b1000_0001), None, "opposite corners");
|
||||
assert_eq!(classify(0b0111), None, "three of a face");
|
||||
assert_eq!(classify(0b0110_1001), None, "four not on a face");
|
||||
assert_eq!(classify(0b0111_1111), None);
|
||||
assert_eq!(classify(0b0011_0011), Some(Pattern::Face));
|
||||
assert_eq!(classify(0b0101_0000), Some(Pattern::Edge));
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,98 @@
|
||||
use super::*;
|
||||
use crate::octree::geom::tests::box_surface;
|
||||
use crate::size::UniformSize;
|
||||
|
||||
#[test]
|
||||
fn a_box_meshes_with_exact_volume_and_tagged_faces() {
|
||||
let s = box_surface([0.0; 3], [10.0, 8.0, 6.0]);
|
||||
let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
|
||||
let r = &m.report;
|
||||
assert!((r.volume - 480.0).abs() < 1e-9 * 480.0, "{r:?}");
|
||||
assert!(r.min_scaled_jacobian > 0.3, "{r:?}");
|
||||
let faces: BTreeSet<u32> = m.boundary.iter().map(|b| b.brep_face).collect();
|
||||
assert_eq!(faces.len(), 6);
|
||||
// Each boundary quad lies in the plane of its face.
|
||||
for b in &m.boundary {
|
||||
let hx = m.hexes[b.element as usize];
|
||||
let q = topo::FACES[(b.face - 1) as usize].map(|i| m.nodes[hx[i] as usize]);
|
||||
let (axis, value) = match b.brep_face {
|
||||
0 => (2, 0.0),
|
||||
1 => (2, 6.0),
|
||||
2 => (1, 0.0),
|
||||
3 => (1, 8.0),
|
||||
4 => (0, 0.0),
|
||||
_ => (0, 10.0),
|
||||
};
|
||||
assert!(q.iter().all(|p| (p[axis] - value).abs() < 1e-9), "{b:?} {q:?}");
|
||||
}
|
||||
}
|
||||
|
||||
/// A closed cylinder of radius `r` and height `h` on the z axis, `n` segments: faces 0 (bottom
|
||||
/// cap), 1 (top cap), 2 (side).
|
||||
pub(crate) fn cylinder_surface(r: f64, h: f64, n: usize) -> BoundarySurface {
|
||||
let mut p = vec![[0.0, 0.0, 0.0], [0.0, 0.0, h]];
|
||||
for z in [0.0, h] {
|
||||
for i in 0..n {
|
||||
let a = std::f64::consts::TAU * i as f64 / n as f64;
|
||||
p.push([r * a.cos(), r * a.sin(), z]);
|
||||
}
|
||||
}
|
||||
let (b, t) = (|i: usize| (2 + i % n) as u32, |i: usize| (2 + n + i % n) as u32);
|
||||
let mut tris = Vec::new();
|
||||
let mut face = Vec::new();
|
||||
for i in 0..n {
|
||||
tris.push([0, b(i + 1), b(i)]);
|
||||
face.push(0);
|
||||
tris.push([1, t(i), t(i + 1)]);
|
||||
face.push(1);
|
||||
tris.push([b(i), b(i + 1), t(i + 1)]);
|
||||
tris.push([b(i), t(i + 1), t(i)]);
|
||||
face.extend([2, 2]);
|
||||
}
|
||||
BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap()
|
||||
}
|
||||
|
||||
/// A UV sphere of radius `r`, one face.
|
||||
pub(crate) fn sphere_surface(r: f64, n: usize) -> BoundarySurface {
|
||||
let rings = n / 2;
|
||||
let mut p = vec![[0.0, 0.0, -r], [0.0, 0.0, r]];
|
||||
for j in 1..rings {
|
||||
let th = std::f64::consts::PI * j as f64 / rings as f64 - std::f64::consts::FRAC_PI_2;
|
||||
for i in 0..n {
|
||||
let a = std::f64::consts::TAU * i as f64 / n as f64;
|
||||
p.push([r * th.cos() * a.cos(), r * th.cos() * a.sin(), r * th.sin()]);
|
||||
}
|
||||
}
|
||||
let id = |j: usize, i: usize| (2 + (j - 1) * n + i % n) as u32;
|
||||
let mut tris = Vec::new();
|
||||
for i in 0..n {
|
||||
tris.push([0, id(1, i + 1), id(1, i)]);
|
||||
tris.push([1, id(rings - 1, i), id(rings - 1, i + 1)]);
|
||||
}
|
||||
for j in 1..rings - 1 {
|
||||
for i in 0..n {
|
||||
tris.push([id(j, i), id(j, i + 1), id(j + 1, i + 1)]);
|
||||
tris.push([id(j, i), id(j + 1, i + 1), id(j + 1, i)]);
|
||||
}
|
||||
}
|
||||
let face = vec![0; tris.len()];
|
||||
BoundarySurface::new(&p, &tris, &face, &crate::surface::SurfaceOptions::default()).unwrap()
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_cylinder_meshes_with_its_rims_captured() {
|
||||
let s = cylinder_surface(5.0, 8.0, 64);
|
||||
let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
|
||||
let r = &m.report;
|
||||
assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}");
|
||||
let faces: BTreeSet<u32> = m.boundary.iter().map(|b| b.brep_face).collect();
|
||||
assert_eq!(faces.len(), 3);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_sphere_meshes() {
|
||||
let s = sphere_surface(5.0, 48);
|
||||
let m = octree_hex_mesh(&s, &UniformSize(2.0), &OctreeOptions::default()).unwrap();
|
||||
let r = &m.report;
|
||||
assert!(((r.volume - r.surface_volume) / r.surface_volume).abs() < 0.01, "{r:?}");
|
||||
}
|
||||
@@ -0,0 +1,404 @@
|
||||
//! Hexahedral mesh topology for the boundary stages: boundary faces, the manifold clean-up of
|
||||
//! the core, padding a boundary with a layer of hexahedra, and pillowing a shrink set.
|
||||
//!
|
||||
//! *Padding* puts one hexahedron on every boundary quad, between the quad and a copy of it, so
|
||||
//! no hexahedron keeps more than one face on the boundary. *Pillowing* a shrink set `S` inserts
|
||||
//! a layer of hexahedra along the faces between `S` and the rest of the mesh; the layer reaches
|
||||
//! the boundary along the loop where `S`'s boundary faces meet the others' (S. A. Mitchell,
|
||||
//! T. J. Tautges, *Pillowing doublets: refining a mesh to ensure that faces share at most one
|
||||
//! edge*, IMR 1995; Gao, Shen, Panozzo 2019, §3.4 use the same insertion around features).
|
||||
|
||||
use std::collections::{BTreeMap, HashMap};
|
||||
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// The six faces of a hex8 as local node positions, right-hand normals outward (the model's
|
||||
/// `ElementKind::Hex8.faces()` order).
|
||||
pub(crate) const FACES: [[usize; 4]; 6] = [[0, 3, 2, 1], [4, 5, 6, 7], [0, 1, 5, 4], [1, 2, 6, 5], [2, 3, 7, 6], [3, 0, 4, 7]];
|
||||
|
||||
/// A face key: the vertices sorted.
|
||||
pub(crate) fn key(q: [u32; 4]) -> [u32; 4] {
|
||||
let mut k = q;
|
||||
k.sort_unstable();
|
||||
k
|
||||
}
|
||||
|
||||
/// The root of `i` in the union-find forest `c`, compressing the path from `i`.
|
||||
pub(crate) fn find(c: &mut [usize], i: usize) -> usize {
|
||||
let mut r = i;
|
||||
while let Some(&p) = c.get(r).filter(|&&p| p != r) {
|
||||
r = p;
|
||||
}
|
||||
if let Some(x) = c.get_mut(i) {
|
||||
*x = r;
|
||||
}
|
||||
r
|
||||
}
|
||||
|
||||
/// A face of an element: (hexahedron, local face 0..6).
|
||||
pub(crate) type FaceRef = (u32, u8);
|
||||
|
||||
/// The vertices of a local face, outward.
|
||||
pub(crate) fn face(h: &[u32; 8], f: u8) -> [u32; 4] {
|
||||
FACES.get(f as usize).map_or([0; 4], |l| l.map(|i| h[i]))
|
||||
}
|
||||
|
||||
/// Every face's elements, keyed by sorted vertices (deterministic map).
|
||||
pub(crate) fn face_map(hexes: &[[u32; 8]]) -> BTreeMap<[u32; 4], Vec<FaceRef>> {
|
||||
let mut m: BTreeMap<[u32; 4], Vec<FaceRef>> = BTreeMap::new();
|
||||
for (i, h) in hexes.iter().enumerate() {
|
||||
for f in 0..6u8 {
|
||||
m.entry(key(face(h, f))).or_default().push((i as u32, f));
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
/// The boundary faces (on one element only), outward, in face-map order.
|
||||
pub(crate) fn boundary(hexes: &[[u32; 8]]) -> Vec<(FaceRef, [u32; 4])> {
|
||||
face_map(hexes)
|
||||
.into_values()
|
||||
.filter(|v| v.len() == 1)
|
||||
.filter_map(|v| v.first().copied())
|
||||
.map(|(h, f)| ((h, f), hexes.get(h as usize).map_or([0; 4], |x| face(x, f))))
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Edges of a quad list with the quads using each (ascending vertex pair → quad indices).
|
||||
pub(crate) fn quad_edges(quads: &[[u32; 4]]) -> BTreeMap<(u32, u32), Vec<u32>> {
|
||||
let mut m: BTreeMap<(u32, u32), Vec<u32>> = BTreeMap::new();
|
||||
for (i, q) in quads.iter().enumerate() {
|
||||
for k in 0..4 {
|
||||
let (a, b) = (q[k], q[(k + 1) % 4]);
|
||||
m.entry((a.min(b), a.max(b))).or_default().push(i as u32);
|
||||
}
|
||||
}
|
||||
m
|
||||
}
|
||||
|
||||
/// The hexahedra to drop so that the boundary of `keep` is a 2-manifold: elements at edges with
|
||||
/// four boundary faces (elements touching along an edge only) and at vertices whose boundary
|
||||
/// faces form more than one fan. Returns the indices to drop (empty when manifold).
|
||||
pub(crate) fn nonmanifold(hexes: &[[u32; 8]], keep: &[bool]) -> Vec<u32> {
|
||||
let kept: Vec<u32> = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect();
|
||||
let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect();
|
||||
let b = boundary(&sub);
|
||||
let quads: Vec<[u32; 4]> = b.iter().map(|(_, q)| *q).collect();
|
||||
let mut bad_vertices: Vec<u32> = Vec::new();
|
||||
for ((a, c), qs) in quad_edges(&quads) {
|
||||
if qs.len() != 2 {
|
||||
bad_vertices.extend([a, c]);
|
||||
}
|
||||
}
|
||||
// Fans: around each vertex, boundary quads connected through shared edges at the vertex.
|
||||
let mut at: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
|
||||
for (i, q) in quads.iter().enumerate() {
|
||||
for &v in q {
|
||||
at.entry(v).or_default().push(i as u32);
|
||||
}
|
||||
}
|
||||
for (v, qs) in &at {
|
||||
let mut comp: Vec<usize> = (0..qs.len()).collect();
|
||||
for i in 0..qs.len() {
|
||||
for j in i + 1..qs.len() {
|
||||
let (qa, qb) = (quads[qs[i] as usize], quads[qs[j] as usize]);
|
||||
// Share an edge through v: another common vertex adjacent to v in both.
|
||||
let nb = |q: [u32; 4]| {
|
||||
let p = q.iter().position(|&x| x == *v).unwrap_or(0);
|
||||
[q[(p + 1) % 4], q[(p + 3) % 4]]
|
||||
};
|
||||
if nb(qa).iter().any(|x| nb(qb).contains(x)) {
|
||||
let (ri, rj) = (find(&mut comp, i), find(&mut comp, j));
|
||||
comp[ri] = rj;
|
||||
}
|
||||
}
|
||||
}
|
||||
let roots = (0..qs.len()).filter(|&i| find(&mut comp, i) == i).count();
|
||||
if roots > 1 {
|
||||
bad_vertices.push(*v);
|
||||
}
|
||||
}
|
||||
bad_vertices.sort_unstable();
|
||||
bad_vertices.dedup();
|
||||
if bad_vertices.is_empty() {
|
||||
return Vec::new();
|
||||
}
|
||||
kept.into_iter().filter(|&i| hexes.get(i as usize).is_some_and(|h| h.iter().any(|v| bad_vertices.binary_search(v).is_ok()))).collect()
|
||||
}
|
||||
|
||||
/// Pad: a hexahedron on each listed boundary face (outward), between it and a copy; the copies
|
||||
/// of the faces' vertices are new points (appended, initially at the original positions).
|
||||
/// Returns the copy of each original vertex.
|
||||
pub(crate) fn pad(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, faces: &[[u32; 4]]) -> Result<BTreeMap<u32, u32>> {
|
||||
let mut copy: BTreeMap<u32, u32> = BTreeMap::new();
|
||||
for q in faces {
|
||||
for &v in q {
|
||||
if let std::collections::btree_map::Entry::Vacant(slot) = copy.entry(v) {
|
||||
let n = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?;
|
||||
points.push(points.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
slot.insert(n);
|
||||
}
|
||||
}
|
||||
}
|
||||
for q in faces {
|
||||
let c = q.map(|v| copy.get(&v).copied().unwrap_or(v));
|
||||
hexes.push([q[0], q[1], q[2], q[3], c[0], c[1], c[2], c[3]]);
|
||||
}
|
||||
Ok(copy)
|
||||
}
|
||||
|
||||
/// One vertex copy made by [`pillow`]: the original, the copy, and the shrink-set elements
|
||||
/// around the original that now use the copy.
|
||||
#[derive(Debug, Clone)]
|
||||
pub(crate) struct PillowCopy {
|
||||
pub original: u32,
|
||||
pub copy: u32,
|
||||
pub elements: Vec<u32>,
|
||||
}
|
||||
|
||||
/// Pillow the shrink set `inside` (one flag per element): a new layer of hexahedra on every face
|
||||
/// between an element of the set and one outside it. Vertices of those faces get one copy per
|
||||
/// group of set elements around them (connected through faces inside the set), which the set's
|
||||
/// elements now use; the new layer joins the originals to the copies. Faces on the boundary are
|
||||
/// not pillowed: where the layer meets the boundary, the set's boundary faces move to the copies,
|
||||
/// so the boundary gains a ring of quads along the loop. New points are appended at the
|
||||
/// originals' positions. Returns the copies made.
|
||||
pub(crate) fn pillow(points: &mut Vec<[f64; 3]>, hexes: &mut Vec<[u32; 8]>, inside: &[bool]) -> Result<Vec<PillowCopy>> {
|
||||
let fm = face_map(hexes);
|
||||
let in_set = |h: u32| inside.get(h as usize).copied().unwrap_or(false);
|
||||
// Faces to pillow, oriented from the outside element into the set.
|
||||
let mut sheet: Vec<(u32, [u32; 4])> = Vec::new();
|
||||
for refs in fm.values() {
|
||||
if let [(a, fa), (b, fb)] = refs.as_slice() {
|
||||
match (in_set(*a), in_set(*b)) {
|
||||
(true, false) => sheet.push((*a, hexes.get(*a as usize).map_or([0; 4], |h| face(h, *fa)))),
|
||||
(false, true) => sheet.push((*b, hexes.get(*b as usize).map_or([0; 4], |h| face(h, *fb)))),
|
||||
_ => {}
|
||||
}
|
||||
}
|
||||
}
|
||||
if sheet.is_empty() {
|
||||
return Ok(Vec::new());
|
||||
}
|
||||
let mut on_sheet: Vec<u32> = sheet.iter().flat_map(|(_, q)| *q).collect();
|
||||
on_sheet.sort_unstable();
|
||||
on_sheet.dedup();
|
||||
// Set elements around each sheet vertex.
|
||||
let mut around: HashMap<u32, Vec<u32>> = HashMap::new();
|
||||
for (i, h) in hexes.iter().enumerate() {
|
||||
if !in_set(i as u32) {
|
||||
continue;
|
||||
}
|
||||
for &v in h {
|
||||
if on_sheet.binary_search(&v).is_ok() {
|
||||
let e = around.entry(v).or_default();
|
||||
if e.last() != Some(&(i as u32)) {
|
||||
e.push(i as u32);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut copies: Vec<PillowCopy> = Vec::new();
|
||||
// (vertex, element) → copy.
|
||||
let mut copy_of: HashMap<(u32, u32), u32> = HashMap::new();
|
||||
for &v in &on_sheet {
|
||||
let mut els = around.get(&v).cloned().unwrap_or_default();
|
||||
els.sort_unstable();
|
||||
els.dedup();
|
||||
// Group by faces inside the set that contain v.
|
||||
let mut group: Vec<usize> = (0..els.len()).collect();
|
||||
for i in 0..els.len() {
|
||||
for j in i + 1..els.len() {
|
||||
let (Some(hi), Some(hj)) = (hexes.get(els[i] as usize), hexes.get(els[j] as usize)) else { continue };
|
||||
let shares = (0..6u8).any(|f| {
|
||||
let q = face(hi, f);
|
||||
q.contains(&v) && (0..6u8).any(|g| key(face(hj, g)) == key(q))
|
||||
});
|
||||
if shares {
|
||||
let (a, b) = (find(&mut group, i), find(&mut group, j));
|
||||
group[a] = b;
|
||||
}
|
||||
}
|
||||
}
|
||||
// Per group: its copy and that copy's index in `copies`.
|
||||
let mut made: BTreeMap<usize, (u32, usize)> = BTreeMap::new();
|
||||
for (i, &e) in els.iter().enumerate() {
|
||||
let r = find(&mut group, i);
|
||||
let (c, at) = match made.get(&r) {
|
||||
Some(&m) => m,
|
||||
None => {
|
||||
let c = u32::try_from(points.len()).map_err(|_| Error::Limit("too many points".into()))?;
|
||||
points.push(points.get(v as usize).copied().unwrap_or([0.0; 3]));
|
||||
copies.push(PillowCopy { original: v, copy: c, elements: Vec::new() });
|
||||
made.insert(r, (c, copies.len() - 1));
|
||||
(c, copies.len() - 1)
|
||||
}
|
||||
};
|
||||
copy_of.insert((v, e), c);
|
||||
if let Some(pc) = copies.get_mut(at) {
|
||||
pc.elements.push(e);
|
||||
}
|
||||
}
|
||||
}
|
||||
// The new layer: original face (seen from outside, normal into the set) below, copies above.
|
||||
let mut layer: Vec<[u32; 8]> = Vec::with_capacity(sheet.len());
|
||||
for (h, q) in &sheet {
|
||||
let r = [q[0], q[3], q[2], q[1]];
|
||||
let c = r.map(|v| copy_of.get(&(v, *h)).copied().unwrap_or(v));
|
||||
layer.push([r[0], r[1], r[2], r[3], c[0], c[1], c[2], c[3]]);
|
||||
}
|
||||
for (i, h) in hexes.iter_mut().enumerate() {
|
||||
if in_set(i as u32) {
|
||||
for v in h.iter_mut() {
|
||||
if let Some(&c) = copy_of.get(&(*v, i as u32)) {
|
||||
*v = c;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
hexes.extend(layer);
|
||||
Ok(copies)
|
||||
}
|
||||
|
||||
/// Hexahedra to drop so that every vertex of the boundary of `keep` can be pushed outward: the
|
||||
/// normals of its boundary quads must lie strictly inside one half-space (their mean has a
|
||||
/// positive product with each, at least `margin`). Where this fails, padding from the vertex
|
||||
/// to the part's surface inverts elements, so for each such vertex the kept element at it with
|
||||
/// the most boundary faces is dropped (the buffer-zone rule of H. Tong, E. Halilaj, Y. J. Zhang,
|
||||
/// *HybridOctree_Hex*, arXiv 2401.05984, 2024, §2.3, in a mean-normal form).
|
||||
pub(crate) fn hidden_vertices(hexes: &[[u32; 8]], keep: &[bool], xyz: &dyn Fn(u32) -> [f64; 3], margin: f64) -> Vec<u32> {
|
||||
let kept: Vec<u32> = (0..hexes.len() as u32).filter(|&i| keep.get(i as usize).copied().unwrap_or(false)).collect();
|
||||
let sub: Vec<[u32; 8]> = kept.iter().filter_map(|&i| hexes.get(i as usize).copied()).collect();
|
||||
let b = boundary(&sub);
|
||||
let mut normals: BTreeMap<u32, Vec<[f64; 3]>> = BTreeMap::new();
|
||||
let mut faces_of: BTreeMap<u32, usize> = BTreeMap::new();
|
||||
for ((h, _), q) in &b {
|
||||
let p = q.map(xyz);
|
||||
let d1 = [0, 1, 2].map(|k| p[2][k] - p[0][k]);
|
||||
let d2 = [0, 1, 2].map(|k| p[3][k] - p[1][k]);
|
||||
let n = [d1[1] * d2[2] - d1[2] * d2[1], d1[2] * d2[0] - d1[0] * d2[2], d1[0] * d2[1] - d1[1] * d2[0]];
|
||||
let l = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
|
||||
if l > 0.0 {
|
||||
for &v in q {
|
||||
normals.entry(v).or_default().push(n.map(|c| c / l));
|
||||
}
|
||||
}
|
||||
*faces_of.entry(*h).or_default() += 1;
|
||||
}
|
||||
let mut drop: Vec<u32> = Vec::new();
|
||||
for (v, ns) in &normals {
|
||||
let mut m = [0.0; 3];
|
||||
for n in ns {
|
||||
for k in 0..3 {
|
||||
m[k] += n[k];
|
||||
}
|
||||
}
|
||||
let l = (m[0] * m[0] + m[1] * m[1] + m[2] * m[2]).sqrt();
|
||||
let ok = l > 0.0 && ns.iter().all(|n| (n[0] * m[0] + n[1] * m[1] + n[2] * m[2]) / l > margin);
|
||||
if ok {
|
||||
continue;
|
||||
}
|
||||
// The kept element at v with the most boundary faces (lowest index on ties).
|
||||
let best = sub
|
||||
.iter()
|
||||
.enumerate()
|
||||
.filter(|(_, h)| h.contains(v))
|
||||
.max_by(|a, b| faces_of.get(&(a.0 as u32)).unwrap_or(&0).cmp(faces_of.get(&(b.0 as u32)).unwrap_or(&0)).then(b.0.cmp(&a.0)))
|
||||
.and_then(|(i, _)| kept.get(i).copied());
|
||||
if let Some(e) = best {
|
||||
drop.push(e);
|
||||
}
|
||||
}
|
||||
drop.sort_unstable();
|
||||
drop.dedup();
|
||||
drop
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn faces_follow_the_model_hex8_face_order() {
|
||||
let model = fieldcraft_model::ElementKind::Hex8.faces();
|
||||
assert_eq!(model.len(), FACES.len());
|
||||
for (m, f) in model.iter().zip(FACES) {
|
||||
assert_eq!(*m, &f[..]);
|
||||
}
|
||||
}
|
||||
|
||||
/// An `n × n × n` block of unit cubes.
|
||||
pub(crate) fn block(n: u32) -> (Vec<[f64; 3]>, Vec<[u32; 8]>) {
|
||||
let id = |x: u32, y: u32, z: u32| x + (n + 1) * (y + (n + 1) * z);
|
||||
let mut p = Vec::new();
|
||||
for z in 0..=n {
|
||||
for y in 0..=n {
|
||||
for x in 0..=n {
|
||||
p.push([x as f64, y as f64, z as f64]);
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut h = Vec::new();
|
||||
for z in 0..n {
|
||||
for y in 0..n {
|
||||
for x in 0..n {
|
||||
h.push([
|
||||
id(x, y, z),
|
||||
id(x + 1, y, z),
|
||||
id(x + 1, y + 1, z),
|
||||
id(x, y + 1, z),
|
||||
id(x, y, z + 1),
|
||||
id(x + 1, y, z + 1),
|
||||
id(x + 1, y + 1, z + 1),
|
||||
id(x, y + 1, z + 1),
|
||||
]);
|
||||
}
|
||||
}
|
||||
}
|
||||
(p, h)
|
||||
}
|
||||
|
||||
fn closed(hexes: &[[u32; 8]]) {
|
||||
for (k, v) in face_map(hexes) {
|
||||
assert!(v.len() <= 2, "face {k:?} on {} elements", v.len());
|
||||
}
|
||||
let b = boundary(hexes);
|
||||
for (_, qs) in quad_edges(&b.iter().map(|x| x.1).collect::<Vec<_>>()) {
|
||||
assert_eq!(qs.len(), 2);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn padding_and_pillowing_keep_a_closed_conforming_mesh() {
|
||||
let (mut p, mut h) = block(3);
|
||||
let b: Vec<[u32; 4]> = boundary(&h).into_iter().map(|x| x.1).collect();
|
||||
assert_eq!(b.len(), 54);
|
||||
let copy = pad(&mut p, &mut h, &b).unwrap();
|
||||
assert_eq!(copy.len(), 56);
|
||||
assert_eq!(h.len(), 27 + 54);
|
||||
closed(&h);
|
||||
// Pillow the padding hexes on the face z = 0 (outer faces down).
|
||||
let pad_start = 27;
|
||||
let inside: Vec<bool> = (0..h.len()).map(|i| i >= pad_start && (0..4).all(|k| p[h[i][k] as usize][2] == 0.0)).collect();
|
||||
assert_eq!(inside.iter().filter(|x| **x).count(), 9);
|
||||
let before = h.len();
|
||||
let copies = pillow(&mut p, &mut h, &inside).unwrap();
|
||||
closed(&h);
|
||||
// Sheet: the 9 inner faces plus the 12 side faces to the neighbouring padding.
|
||||
assert_eq!(h.len(), before + 9 + 12);
|
||||
// The boundary gains the ring: 54 + 12 quads.
|
||||
assert_eq!(boundary(&h).len(), 54 + 12);
|
||||
assert!(!copies.is_empty());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nonmanifold_edges_and_vertices_are_found() {
|
||||
let (_, h) = block(2);
|
||||
// Two cubes touching along an edge.
|
||||
let keep = |ids: &[usize]| (0..8).map(|i| ids.contains(&i)).collect::<Vec<bool>>();
|
||||
assert!(!nonmanifold(&h, &keep(&[0, 3])).is_empty(), "edge contact");
|
||||
assert!(!nonmanifold(&h, &keep(&[0, 7])).is_empty(), "vertex contact");
|
||||
assert!(nonmanifold(&h, &keep(&[0, 1, 3])).is_empty());
|
||||
assert!(nonmanifold(&h, &keep(&[0, 1, 2, 3, 4, 5, 6, 7])).is_empty());
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user