From 68e9bcde9266f5d3c8d63b945ecee2428f892e59 Mon Sep 17 00:00:00 2001 From: jyviko <2501648+jyviko@users.noreply.github.com> Date: Fri, 16 Oct 2026 12:48:57 -0400 Subject: [PATCH] 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 Claude-Session: https://claude.ai/code/session_01UhhZMBTviDQxZvwnyP89fo --- crates/mesh/src/lib.rs | 1 + crates/mesh/src/tet/complex.rs | 24 ++- crates/mesh/src/tet/mod.rs | 14 ++ crates/mesh/src/transition.rs | 357 +++++++++++++++++++++++++++++++++ 4 files changed, 392 insertions(+), 4 deletions(-) create mode 100644 crates/mesh/src/transition.rs diff --git a/crates/mesh/src/lib.rs b/crates/mesh/src/lib.rs index 4cd6626..f74781b 100644 --- a/crates/mesh/src/lib.rs +++ b/crates/mesh/src/lib.rs @@ -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)] diff --git a/crates/mesh/src/tet/complex.rs b/crates/mesh/src/tet/complex.rs index 070ed96..726c41a 100644 --- a/crates/mesh/src/tet/complex.rs +++ b/crates/mesh/src/tet/complex.rs @@ -195,6 +195,8 @@ pub struct TetComplex { /// only at it, so stars never miss a piece). vtets: Vec>, curves: IdMap<(u32, u32), u32>, + /// Boundary face tags no change may remove or add a face of ([`Self::freeze_faces`]). + frozen: Vec, mark: Vec, 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) { 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). diff --git a/crates/mesh/src/tet/mod.rs b/crates/mesh/src/tet/mod.rs index 255e879..72e1a8f 100644 --- a/crates/mesh/src/tet/mod.rs +++ b/crates/mesh/src/tet/mod.rs @@ -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, + /// 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, } 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> = 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), diff --git a/crates/mesh/src/transition.rs b/crates/mesh/src/transition.rs new file mode 100644 index 0000000..b1f776d --- /dev/null +++ b/crates/mesh/src/transition.rs @@ -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, + /// 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 { + if quads.is_empty() { + return Err(Error::Invalid("no interface quadrilaterals".into())); + } + let cut: BTreeSet = 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> = 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 = removed.iter().flatten().copied().collect(); + let rim_nodes: Vec = rim.keys().copied().collect(); + let mut nearest: BTreeMap = 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 = Vec::with_capacity(triangles.capacity()); + for &(t, f) in &kept { + let rim_edges: Vec = (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::() / 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>, from: u32, to: u32) -> Result> { + 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 = 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> { + 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()); + } +}