forked from 1b-agent/fieldcraft-upstream
mesh: pyramid transition surface and frozen faces in the tet mesher
transition::pyramid_surface puts a pyramid on every hex interface quad and builds the tet region's surface closed by the pyramid sides, re-stitching the neighbouring faces to the hex rim nodes. QualityOptions::frozen_faces keeps those sides exactly (TetComplex refuses changes that remove or add a frozen face; their vertices are corners), and apex_points checks the result. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo
This commit is contained in:
jyvikoandClaude Opus 5.5 authored and Irakli Kourtis committed
1 parent
3da332be77
commit
68e9bcde92
4 files changed
+392
-4
No files matched your search
@@ -24,6 +24,7 @@ pub mod seeds;
|
||||
pub mod size;
|
||||
pub mod surface;
|
||||
pub mod tet;
|
||||
pub mod transition;
|
||||
mod vec3;
|
||||
|
||||
#[derive(Debug, Clone, PartialEq, thiserror::Error)]
|
||||
|
||||
@@ -195,6 +195,8 @@ pub struct TetComplex {
|
||||
/// only at it, so stars never miss a piece).
|
||||
vtets: Vec<Vec<u32>>,
|
||||
curves: IdMap<(u32, u32), u32>,
|
||||
/// Boundary face tags no change may remove or add a face of ([`Self::freeze_faces`]).
|
||||
frozen: Vec<u32>,
|
||||
mark: Vec<u32>,
|
||||
epoch: u32,
|
||||
live: usize,
|
||||
@@ -276,6 +278,7 @@ impl TetComplex {
|
||||
tags: vec![[NONE; 4]; tets.len()],
|
||||
free: Vec::new(),
|
||||
curves: IdMap::default(),
|
||||
frozen: Vec::new(),
|
||||
mark: vec![0; tets.len()],
|
||||
epoch: 0,
|
||||
live: tets.len(),
|
||||
@@ -406,6 +409,19 @@ impl TetComplex {
|
||||
e
|
||||
}
|
||||
|
||||
/// Keep the boundary faces tagged with `tags` exactly as they are: [`Self::replace`] refuses
|
||||
/// every change that would remove one of them (a split, flip or collapse) or add one. Their
|
||||
/// vertices must also be [`VertexKind::Corner`]s, so nothing moves them.
|
||||
pub fn freeze_faces(&mut self, tags: &[u32]) {
|
||||
self.frozen = tags.to_vec();
|
||||
self.frozen.sort_unstable();
|
||||
self.frozen.dedup();
|
||||
}
|
||||
|
||||
fn is_frozen(&self, tag: u32) -> bool {
|
||||
self.frozen.binary_search(&tag).is_ok()
|
||||
}
|
||||
|
||||
pub fn set_curve(&mut self, a: u32, b: u32, curve: Option<u32>) {
|
||||
let key = sorted2(a, b);
|
||||
let before = match curve {
|
||||
@@ -881,13 +897,13 @@ impl TetComplex {
|
||||
used_declared += usize::from(declared_face.is_some());
|
||||
if let Some((n, back, tag, seen)) = outside.remove(&key) {
|
||||
// The new tetrahedron lies on the side of the face the old one did.
|
||||
if !same_cycle(face, seen) {
|
||||
if !same_cycle(face, seen) || (n == NONE && self.is_frozen(tag) && declared_face.is_some_and(|d| d.1 != tag)) {
|
||||
return Ok(None);
|
||||
}
|
||||
links.push((i, k, if n == NONE { Link::Boundary(declared_face.map_or(tag, |d| d.1)) } else { Link::Old(n, back) }));
|
||||
} else if let Some((oriented, tag)) = declared_face {
|
||||
// A new boundary face, oriented as declared (into the mesh).
|
||||
if !same_cycle(face, oriented) {
|
||||
if !same_cycle(face, oriented) || self.is_frozen(tag) {
|
||||
return Ok(None);
|
||||
}
|
||||
links.push((i, k, Link::Boundary(tag)));
|
||||
@@ -898,8 +914,8 @@ impl TetComplex {
|
||||
if used_declared != declared.len() {
|
||||
return Ok(None);
|
||||
}
|
||||
// What is left outside must be boundary faces the change removes.
|
||||
if outside.values().any(|(n, _, _, _)| *n != NONE) || (!outside.is_empty() && change.new_boundary.is_empty()) {
|
||||
// What is left outside must be boundary faces the change removes, and none frozen.
|
||||
if outside.values().any(|(n, _, tag, _)| *n != NONE || self.is_frozen(*tag)) || (!outside.is_empty() && change.new_boundary.is_empty()) {
|
||||
return Ok(None);
|
||||
}
|
||||
// Feature edges strictly inside the old region must survive (or be the split one).
|
||||
|
||||
@@ -97,6 +97,11 @@ pub struct QualityOptions {
|
||||
/// grows from it at the grading rate. A seed finer than the global size refines along its
|
||||
/// edge; one coarser leaves it to the global size.
|
||||
pub seeds: Vec<crate::seeds::EdgeSeed>,
|
||||
/// B-rep face tags whose boundary triangles are kept exactly: their vertices stay where they
|
||||
/// are and no triangle of theirs is split, flipped or removed. The faces a pyramid
|
||||
/// transition shares with the tetrahedra ([`crate::transition`]) are frozen, so the two
|
||||
/// meshes stay conforming.
|
||||
pub frozen_faces: Vec<u32>,
|
||||
}
|
||||
|
||||
impl Default for QualityOptions {
|
||||
@@ -117,6 +122,7 @@ impl Default for QualityOptions {
|
||||
max_tets: 5_000_000,
|
||||
max_points: 1_000_000,
|
||||
seeds: Vec::new(),
|
||||
frozen_faces: Vec::new(),
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -359,6 +365,8 @@ pub fn complex_from_cdt(surface: &BoundarySurface, cdt: &TetMesh, o: &QualityOpt
|
||||
// Vertex kinds from the boundary faces and curve edges at each vertex.
|
||||
let mut face_tag: Vec<Option<u32>> = vec![None; cdt.points.len()];
|
||||
let mut mixed = vec![false; cdt.points.len()];
|
||||
let mut frozen = vec![false; cdt.points.len()];
|
||||
cx.freeze_faces(&o.frozen_faces);
|
||||
for (t, vs) in cx.live_tets() {
|
||||
let (ns, tg) = (cx.adj(t)?, cx.tags(t)?);
|
||||
for k in 0..4 {
|
||||
@@ -367,6 +375,11 @@ pub fn complex_from_cdt(surface: &BoundarySurface, cdt: &TetMesh, o: &QualityOpt
|
||||
}
|
||||
for v in face_toward(vs, k) {
|
||||
let i = v as usize;
|
||||
if o.frozen_faces.contains(&tg[k])
|
||||
&& let Some(f) = frozen.get_mut(i)
|
||||
{
|
||||
*f = true;
|
||||
}
|
||||
match face_tag.get(i).copied().flatten() {
|
||||
None => {
|
||||
if let Some(s) = face_tag.get_mut(i) {
|
||||
@@ -390,6 +403,7 @@ pub fn complex_from_cdt(surface: &BoundarySurface, cdt: &TetMesh, o: &QualityOpt
|
||||
}
|
||||
for v in 0..cdt.points.len() as u32 {
|
||||
let kind = match (face_tag.get(v as usize).copied().flatten(), curve_count.get(&v).map(Vec::as_slice)) {
|
||||
_ if frozen.get(v as usize).copied().unwrap_or(false) => VertexKind::Corner,
|
||||
(None, _) => VertexKind::Interior,
|
||||
(Some(tag), None) if !mixed.get(v as usize).copied().unwrap_or(true) => VertexKind::Face(tag),
|
||||
(Some(_), Some(&[c1, c2])) if c1 == c2 => VertexKind::Curve(c1),
|
||||
|
||||
@@ -0,0 +1,357 @@
|
||||
//! Conforming hex–tet transitions: a layer of pyramids on a hexahedral mesh's interface faces,
|
||||
//! and the boundary surface of the tetrahedral region that meets them.
|
||||
//!
|
||||
//! A hexahedron's face is a quadrilateral and a tetrahedron's a triangle, so a tet mesh cannot
|
||||
//! share a hex mesh's faces. A pyramid can: its base is the quadrilateral and its four sides are
|
||||
//! triangles. [`pyramid_surface`] puts one pyramid on every interface quadrilateral (its apex
|
||||
//! inside the tet region, about 0.6 of the quadrilateral's mean edge above its centre) and builds
|
||||
//! the surface the tetrahedra fill: the tet region's own surface without the interface faces,
|
||||
//! closed by the pyramids' sides. Where the interface meets the tet region's other faces (the
|
||||
//! rim), those faces' triangles are re-stitched to the hex mesh's rim nodes: each rim vertex of
|
||||
//! the tessellation moves to the nearest hex rim node and each rim edge becomes a fan over the hex
|
||||
//! rim nodes it spans, so the rim of the surface is exactly the hex mesh's. The tet mesher then
|
||||
//! keeps the pyramids' sides unchanged ([`crate::tet::QualityOptions::frozen_faces`]) and
|
||||
//! [`apex_points`] checks that it did.
|
||||
//!
|
||||
//! The result is one conforming mesh: hexahedra, pyramids on their faces, tetrahedra on the
|
||||
//! pyramids' sides and on the rest of the region. The surface is checked by
|
||||
//! [`BoundarySurface::new`] like any other (closed, manifold, no self-intersection): a pyramid
|
||||
//! that reaches through a thin region, or a rim that cannot be stitched, is an error, and the
|
||||
//! caller falls back to another join.
|
||||
//!
|
||||
//! # Sources
|
||||
//! - S. J. Owen, S. A. Canann, S. Saigal, *Pyramid elements for maintaining tetrahedra to
|
||||
//! hexahedra conformability*, AMD Vol. 220, Trends in Unstructured Mesh Generation, ASME 1997,
|
||||
//! pp. 123–129: a pyramid on each quadrilateral of the interface, its apex placed off the face
|
||||
//! inside the tetrahedral region.
|
||||
//! - L. Bedrosian, *Shape functions and integration formulas for three-dimensional finite
|
||||
//! element analysis*, IJNME 35(1):95–108, 1992: the pyramid5 element that makes the
|
||||
//! transition conforming (`fieldcraft_fem::basis`).
|
||||
|
||||
use std::collections::{BTreeMap, BTreeSet, VecDeque};
|
||||
|
||||
use crate::surface::BoundarySurface;
|
||||
use crate::tet::QualityMesh;
|
||||
use crate::tet::complex::face_toward;
|
||||
use crate::{Error, Result};
|
||||
|
||||
/// The apex's height above the quadrilateral's centre, over the quadrilateral's mean edge.
|
||||
pub const APEX_HEIGHT: f64 = 0.6;
|
||||
|
||||
/// The tet region's surface for a pyramid transition ([`pyramid_surface`]).
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct Transition {
|
||||
pub positions: Vec<[f64; 3]>,
|
||||
pub triangles: Vec<[u32; 3]>,
|
||||
pub tri_face: Vec<u32>,
|
||||
/// One apex per interface quadrilateral, in their order.
|
||||
pub apexes: Vec<[f64; 3]>,
|
||||
/// Each pyramid's four side triangles (positions), outward from the tet region, quad by quad.
|
||||
pub sides: Vec<[[f64; 3]; 3]>,
|
||||
}
|
||||
|
||||
fn key(p: [f64; 3]) -> [u64; 3] {
|
||||
p.map(f64::to_bits)
|
||||
}
|
||||
|
||||
fn sub(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
|
||||
[a[0] - b[0], a[1] - b[1], a[2] - b[2]]
|
||||
}
|
||||
|
||||
fn dist(a: [f64; 3], b: [f64; 3]) -> f64 {
|
||||
let d = sub(a, b);
|
||||
(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]).sqrt()
|
||||
}
|
||||
|
||||
/// The tet region's surface with a pyramid on every quadrilateral of `quads` (corners in order,
|
||||
/// their right-hand normal pointing into the tet region, i.e. out of the hexahedra), in place of
|
||||
/// the triangles of `surface` on the B-rep faces `cut`. The pyramids' sides are tagged `tag`.
|
||||
pub fn pyramid_surface(surface: &BoundarySurface, cut: &[u32], quads: &[[[f64; 3]; 4]], tag: u32) -> Result<Transition> {
|
||||
if quads.is_empty() {
|
||||
return Err(Error::Invalid("no interface quadrilaterals".into()));
|
||||
}
|
||||
let cut: BTreeSet<u32> = cut.iter().copied().collect();
|
||||
if cut.contains(&tag) || surface.tri_face().contains(&tag) {
|
||||
return Err(Error::Invalid(format!("the pyramid tag {tag} is a face of the region")));
|
||||
}
|
||||
let mut positions: Vec<[f64; 3]> = surface.points().to_vec();
|
||||
let (mut kept, mut removed): (Vec<([u32; 3], u32)>, Vec<[u32; 3]>) = (Vec::new(), Vec::new());
|
||||
for (t, &f) in surface.triangles().iter().zip(surface.tri_face()) {
|
||||
if cut.contains(&f) {
|
||||
removed.push(*t);
|
||||
} else {
|
||||
kept.push((*t, f));
|
||||
}
|
||||
}
|
||||
if removed.is_empty() {
|
||||
return Err(Error::Invalid("the interface faces have no triangles".into()));
|
||||
}
|
||||
// The hex nodes of the interface, after the surface's points.
|
||||
let mut hex_index: BTreeMap<[u64; 3], u32> = BTreeMap::new();
|
||||
let mut quad_nodes: Vec<[u32; 4]> = Vec::with_capacity(quads.len());
|
||||
for q in quads {
|
||||
if q.iter().flatten().any(|c| !c.is_finite()) {
|
||||
return Err(Error::Invalid("an interface quadrilateral has a non-finite corner".into()));
|
||||
}
|
||||
let mut ids = [0u32; 4];
|
||||
for (slot, p) in ids.iter_mut().zip(q) {
|
||||
let next = u32::try_from(positions.len()).map_err(|_| Error::Limit("too many points".into()))?;
|
||||
*slot = *hex_index.entry(key(*p)).or_insert_with(|| {
|
||||
positions.push(*p);
|
||||
next
|
||||
});
|
||||
}
|
||||
quad_nodes.push(ids);
|
||||
}
|
||||
// The rim of the quadrilaterals: edges of one quadrilateral, as a graph.
|
||||
let mut edge_use: BTreeMap<(u32, u32), usize> = BTreeMap::new();
|
||||
for q in &quad_nodes {
|
||||
for i in 0..4 {
|
||||
let (a, b) = (q[i], q[(i + 1) % 4]);
|
||||
*edge_use.entry((a.min(b), a.max(b))).or_default() += 1;
|
||||
}
|
||||
}
|
||||
if edge_use.values().any(|&n| n > 2) {
|
||||
return Err(Error::Invalid("an edge of the interface is shared by more than two quadrilaterals".into()));
|
||||
}
|
||||
let mut rim: BTreeMap<u32, Vec<u32>> = BTreeMap::new();
|
||||
for (&(a, b), _) in edge_use.iter().filter(|(_, n)| **n == 1) {
|
||||
rim.entry(a).or_default().push(b);
|
||||
rim.entry(b).or_default().push(a);
|
||||
}
|
||||
if rim.values().any(|n| n.len() != 2) {
|
||||
return Err(Error::Unsupported("the interface's rim is not a set of simple loops (it touches itself at a node)".into()));
|
||||
}
|
||||
// The surface's rim: the kept triangles' edges whose reverse is a removed triangle's edge.
|
||||
let removed_edges: BTreeSet<(u32, u32)> = removed.iter().flat_map(|t| [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]).collect();
|
||||
let removed_vertices: BTreeSet<u32> = removed.iter().flatten().copied().collect();
|
||||
let rim_nodes: Vec<u32> = rim.keys().copied().collect();
|
||||
let mut nearest: BTreeMap<u32, u32> = BTreeMap::new();
|
||||
for &(t, _) in &kept {
|
||||
for v in t {
|
||||
if !removed_vertices.contains(&v) || nearest.contains_key(&v) {
|
||||
continue;
|
||||
}
|
||||
let p = positions.get(v as usize).copied().ok_or_else(|| Error::Failed("a surface vertex is out of range".into()))?;
|
||||
let mut best: Option<(f64, u32)> = None;
|
||||
for &h in &rim_nodes {
|
||||
let d = positions.get(h as usize).map_or(f64::INFINITY, |q| dist(p, *q));
|
||||
if best.is_none_or(|b| d < b.0) {
|
||||
best = Some((d, h));
|
||||
}
|
||||
}
|
||||
if let Some((_, h)) = best {
|
||||
nearest.insert(v, h);
|
||||
}
|
||||
}
|
||||
}
|
||||
let at = |v: u32| nearest.get(&v).copied().unwrap_or(v);
|
||||
let mut triangles: Vec<[u32; 3]> = Vec::with_capacity(kept.len() + 4 * quads.len());
|
||||
let mut tri_face: Vec<u32> = Vec::with_capacity(triangles.capacity());
|
||||
for &(t, f) in &kept {
|
||||
let rim_edges: Vec<usize> = (0..3).filter(|&i| removed_edges.contains(&(t[(i + 1) % 3], t[i]))).collect();
|
||||
match rim_edges.as_slice() {
|
||||
[] => {
|
||||
let m = t.map(at);
|
||||
if m[0] == m[1] || m[1] == m[2] || m[2] == m[0] {
|
||||
return Err(Error::Unsupported(
|
||||
"a triangle beside the interface collapses onto its rim (the tessellation is coarser than the hex mesh there)".into(),
|
||||
));
|
||||
}
|
||||
triangles.push(m);
|
||||
tri_face.push(f);
|
||||
}
|
||||
[i] => {
|
||||
let (u, v, w) = (t[*i], t[(i + 1) % 3], t[(i + 2) % 3]);
|
||||
let path = rim_path(&rim, at(u), at(v))?;
|
||||
let apex = at(w);
|
||||
for pair in path.windows(2) {
|
||||
if let [a, b] = *pair {
|
||||
if apex == a || apex == b {
|
||||
return Err(Error::Unsupported("a triangle beside the interface lies along the rim on two sides".into()));
|
||||
}
|
||||
triangles.push([a, b, apex]);
|
||||
tri_face.push(f);
|
||||
}
|
||||
}
|
||||
}
|
||||
_ => return Err(Error::Unsupported("a triangle beside the interface lies along two of its edges".into())),
|
||||
}
|
||||
}
|
||||
// The pyramids: an apex above each quadrilateral, the sides toward the tet region.
|
||||
let mut apexes = Vec::with_capacity(quads.len());
|
||||
let mut sides = Vec::with_capacity(4 * quads.len());
|
||||
for (q, ids) in quads.iter().zip(&quad_nodes) {
|
||||
// Newell normal (into the tet region) and the mean edge.
|
||||
let mut n = [0.0; 3];
|
||||
let mut edge = 0.0;
|
||||
for i in 0..4 {
|
||||
let (p, r) = (q[i], q[(i + 1) % 4]);
|
||||
n = [n[0] + (p[1] - r[1]) * (p[2] + r[2]), n[1] + (p[2] - r[2]) * (p[0] + r[0]), n[2] + (p[0] - r[0]) * (p[1] + r[1])];
|
||||
edge += dist(p, r) / 4.0;
|
||||
}
|
||||
let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
|
||||
if !(len > 0.0 && edge > 0.0 && len.is_finite()) {
|
||||
return Err(Error::Invalid("an interface quadrilateral is degenerate".into()));
|
||||
}
|
||||
let c = [0, 1, 2].map(|k| q.iter().map(|p| p[k]).sum::<f64>() / 4.0);
|
||||
let apex = [0, 1, 2].map(|k| c[k] + n[k] / len * APEX_HEIGHT * edge);
|
||||
let a = u32::try_from(positions.len()).map_err(|_| Error::Limit("too many points".into()))?;
|
||||
positions.push(apex);
|
||||
apexes.push(apex);
|
||||
for i in 0..4 {
|
||||
// The pyramid's side on edge i → i+1 faces out of the pyramid; the region's boundary
|
||||
// there faces into it.
|
||||
let (p, r) = (ids[i], ids[(i + 1) % 4]);
|
||||
triangles.push([r, p, a]);
|
||||
tri_face.push(tag);
|
||||
sides.push([q[(i + 1) % 4], q[i], apex]);
|
||||
}
|
||||
}
|
||||
// Drop the points no triangle uses (the interface's own and the moved rim vertices).
|
||||
let mut used = vec![u32::MAX; positions.len()];
|
||||
let mut compact = Vec::new();
|
||||
for t in &mut triangles {
|
||||
for v in t.iter_mut() {
|
||||
let slot = used.get_mut(*v as usize).ok_or_else(|| Error::Failed("a vertex is out of range".into()))?;
|
||||
if *slot == u32::MAX {
|
||||
*slot = compact.len() as u32;
|
||||
compact.push(positions.get(*v as usize).copied().unwrap_or([0.0; 3]));
|
||||
}
|
||||
*v = *slot;
|
||||
}
|
||||
}
|
||||
Ok(Transition { positions: compact, triangles, tri_face, apexes, sides })
|
||||
}
|
||||
|
||||
/// The rim nodes from `from` to `to` along the shorter way round their loop.
|
||||
fn rim_path(rim: &BTreeMap<u32, Vec<u32>>, from: u32, to: u32) -> Result<Vec<u32>> {
|
||||
if from == to {
|
||||
return Ok(vec![from]);
|
||||
}
|
||||
// Breadth-first: the loop's two ways meet at the far side, so this is the shorter one.
|
||||
let mut previous: BTreeMap<u32, u32> = BTreeMap::new();
|
||||
let mut queue = VecDeque::from([from]);
|
||||
previous.insert(from, from);
|
||||
while let Some(v) = queue.pop_front() {
|
||||
if v == to {
|
||||
break;
|
||||
}
|
||||
for &w in rim.get(&v).map(Vec::as_slice).unwrap_or_default() {
|
||||
if let std::collections::btree_map::Entry::Vacant(e) = previous.entry(w) {
|
||||
e.insert(v);
|
||||
queue.push_back(w);
|
||||
}
|
||||
}
|
||||
}
|
||||
let mut path = vec![to];
|
||||
let mut v = to;
|
||||
while v != from {
|
||||
v = *previous.get(&v).ok_or_else(|| Error::Unsupported("the tessellation's rim and the hex mesh's rim are on different loops".into()))?;
|
||||
path.push(v);
|
||||
if path.len() > rim.len() + 1 {
|
||||
return Err(Error::Failed("the rim walk did not end".into()));
|
||||
}
|
||||
}
|
||||
path.reverse();
|
||||
Ok(path)
|
||||
}
|
||||
|
||||
/// The tet mesh's point index of each pyramid's apex, after checking that the mesh kept every
|
||||
/// pyramid side as one of its boundary faces (so the pyramids and tetrahedra share them).
|
||||
pub fn apex_points(mesh: &QualityMesh, t: &Transition) -> Result<Vec<u32>> {
|
||||
let index: BTreeMap<[u64; 3], u32> = mesh.points.iter().enumerate().map(|(i, p)| (key(*p), i as u32)).collect();
|
||||
let find = |p: [f64; 3]| index.get(&key(p)).copied();
|
||||
let mut faces: BTreeSet<[u32; 3]> = BTreeSet::new();
|
||||
for b in &mesh.boundary {
|
||||
let Some(e) = mesh.elements.get(b.tet as usize) else { continue };
|
||||
let Some(vs) = e.get(..4).and_then(|s| <[u32; 4]>::try_from(s).ok()) else { continue };
|
||||
let mut f = face_toward(vs, usize::from(b.opposite));
|
||||
f.sort_unstable();
|
||||
faces.insert(f);
|
||||
}
|
||||
let mut missing = 0usize;
|
||||
for s in &t.sides {
|
||||
match (find(s[0]), find(s[1]), find(s[2])) {
|
||||
(Some(a), Some(b), Some(c)) => {
|
||||
let mut f = [a, b, c];
|
||||
f.sort_unstable();
|
||||
missing += usize::from(!faces.contains(&f));
|
||||
}
|
||||
_ => missing += 1,
|
||||
}
|
||||
}
|
||||
if missing > 0 {
|
||||
return Err(Error::Failed(format!("the tet mesh changed {missing} of the {} pyramid sides it had to keep", t.sides.len())));
|
||||
}
|
||||
t.apexes.iter().map(|p| find(*p).ok_or_else(|| Error::Failed("a pyramid's apex is not a vertex of the tet mesh".into()))).collect()
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
use crate::cdt::{CdtOptions, tetrahedralize};
|
||||
use crate::surface::SurfaceOptions;
|
||||
use crate::tet::{QualityOptions, quality_mesh};
|
||||
|
||||
/// A box [0, 2] × [0, 2] × [0, 1] as twelve triangles (two per side, outward), sides tagged
|
||||
/// 0..6 (−x, +x, −y, +y, −z, +z).
|
||||
fn slab() -> BoundarySurface {
|
||||
let p: Vec<[f64; 3]> = (0..8).map(|i| [2.0 * (i & 1) as f64, 2.0 * ((i >> 1) & 1) as f64, ((i >> 2) & 1) as f64]).collect();
|
||||
let quads = [([0, 4, 6, 2], 0), ([1, 3, 7, 5], 1), ([0, 1, 5, 4], 2), ([2, 6, 7, 3], 3), ([0, 2, 3, 1], 4), ([4, 5, 7, 6], 5)];
|
||||
let (mut t, mut f) = (Vec::new(), Vec::new());
|
||||
for (q, tag) in quads {
|
||||
t.push([q[0], q[1], q[2]]);
|
||||
t.push([q[0], q[2], q[3]]);
|
||||
f.extend([tag, tag]);
|
||||
}
|
||||
BoundarySurface::new(&p, &t, &f, &SurfaceOptions::default()).unwrap()
|
||||
}
|
||||
|
||||
/// The bottom (z = 0) as a 4 × 4 grid of hex faces, their normals up into the slab.
|
||||
fn bottom_quads() -> Vec<[[f64; 3]; 4]> {
|
||||
let mut q = Vec::new();
|
||||
for j in 0..4 {
|
||||
for i in 0..4 {
|
||||
let (x, y, h) = (0.5 * i as f64, 0.5 * j as f64, 0.5);
|
||||
q.push([[x, y, 0.0], [x + h, y, 0.0], [x + h, y + h, 0.0], [x, y + h, 0.0]]);
|
||||
}
|
||||
}
|
||||
q
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn pyramids_on_a_face_are_kept_by_the_tet_mesh() {
|
||||
let s = slab();
|
||||
let quads = bottom_quads();
|
||||
let t = pyramid_surface(&s, &[4], &quads, 6).unwrap();
|
||||
assert_eq!(t.apexes.len(), 16);
|
||||
let surface = BoundarySurface::new(&t.positions, &t.triangles, &t.tri_face, &SurfaceOptions::default()).unwrap();
|
||||
// The slab less the pyramids (each 0.5² × 0.3 / 3).
|
||||
let pyramids = 16.0 * 0.25 * (APEX_HEIGHT * 0.5) / 3.0;
|
||||
assert!((surface.volume() - (4.0 - pyramids)).abs() < 1e-12, "{}", surface.volume());
|
||||
let cdt = tetrahedralize(&surface, &CdtOptions::default()).unwrap();
|
||||
let options = QualityOptions { size: Some(0.3), frozen_faces: vec![6], ..QualityOptions::default() };
|
||||
let mesh = quality_mesh(&surface, &cdt, Vec::new(), &options).unwrap();
|
||||
let apexes = apex_points(&mesh, &t).unwrap();
|
||||
assert_eq!(apexes.len(), 16);
|
||||
// Without freezing, refinement splits the pyramid sides and the check says so.
|
||||
let free = quality_mesh(&surface, &cdt, Vec::new(), &QualityOptions { size: Some(0.1), ..QualityOptions::default() }).unwrap();
|
||||
assert!(apex_points(&free, &t).is_err());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn a_pyramid_through_a_thin_region_is_refused() {
|
||||
let s = slab();
|
||||
// 1 × 1 quads: apexes 0.6 up, inside the 1 mm slab; one 2 × 2 quad: its apex 1.2 up, outside.
|
||||
let quads: Vec<[[f64; 3]; 4]> = (0..4)
|
||||
.map(|k| ((k % 2) as f64, (k / 2) as f64))
|
||||
.map(|(x, y)| [[x, y, 0.0], [x + 1.0, y, 0.0], [x + 1.0, y + 1.0, 0.0], [x, y + 1.0, 0.0]])
|
||||
.collect();
|
||||
let t = pyramid_surface(&s, &[4], &quads, 6).unwrap();
|
||||
let big: Vec<[[f64; 3]; 4]> = vec![[[0.0, 0.0, 0.0], [2.0, 0.0, 0.0], [2.0, 2.0, 0.0], [0.0, 2.0, 0.0]]];
|
||||
let tall = pyramid_surface(&s, &[4], &big, 6).unwrap();
|
||||
assert!(BoundarySurface::new(&t.positions, &t.triangles, &t.tri_face, &SurfaceOptions::default()).is_ok());
|
||||
assert!(BoundarySurface::new(&tall.positions, &tall.triangles, &tall.tri_face, &SurfaceOptions::default()).is_err());
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user