From 9cbb85f703f4cb432dd7b82ec02f80281eba6a55 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 09:50:11 +0200 Subject: [PATCH 01/14] feat(volume_mesh): scaffold prism/pyramid cell types and shape data --- .../src/volume_mesh/cell_data.rs | 109 ++++++++++++++++++ .../src/volume_mesh/mod.rs | 31 ++++- 2 files changed, 138 insertions(+), 2 deletions(-) create mode 100644 crates/polyscope-structures/src/volume_mesh/cell_data.rs diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs new file mode 100644 index 0000000..c0c4f16 --- /dev/null +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -0,0 +1,109 @@ +// Items are wired in by subsequent tasks (tet decomposition + face dispatch). +// The allow is removed at that point. +#![allow(dead_code)] + +//! Static shape data for each `VolumeCellType`. +//! +//! Defines, for each cell type: +//! - the unique vertex indices on each face (used for canonical face hashing) +//! - the per-face triangulation stencil (used for rendering and edge detection) +//! - the tet-decomposition pattern (used for slicing and isosurface extraction) +//! +//! Index conventions mirror upstream C++ Polyscope (`src/volume_mesh.cpp`, +//! `stencilTet`, `stencilHex`, `stencilPrism`, `stencilPyramid`): +//! - Tet: 4 verts in slots 0..3, sentinels in 4..7. 4 triangular faces. +//! - Hex: 8 verts in 0..7. 6 quadrilateral faces, each split into 2 triangles. +//! - Prism: 6 verts in 0..5, sentinels in 6..7. Bottom tri (0,1,2), top tri (3,4,5). +//! 5 faces: 1 tri (bottom) + 3 quads (sides) + 1 tri (top) = 8 triangles. +//! - Pyramid: 5 verts in 0..4, sentinels in 5..7. Base quad (0,1,2,3), apex (4). +//! 5 faces: 1 quad (base) + 4 tris (sides) = 6 triangles. + +use super::VolumeCellType; + +/// A face is described by (unique vertex slot list, triangulation). +/// +/// `polygon` lists the unique cell-local vertex slots in CCW order around the +/// face (3 for triangles, 4 for quads). `triangulation` is the list of triangle +/// stencils used to render the face. +pub struct FaceData { + pub polygon: &'static [usize], + pub triangulation: &'static [[usize; 3]], +} + +// ===== Tet ===== +const TET_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, + FaceData { polygon: &[0, 1, 3], triangulation: &[[0, 1, 3]] }, + FaceData { polygon: &[0, 3, 2], triangulation: &[[0, 3, 2]] }, + FaceData { polygon: &[1, 2, 3], triangulation: &[[1, 2, 3]] }, +]; + +// ===== Hex ===== +// Numbered like in the VTK file-formats diagram, with slots 6 and 7 swapped +// to match upstream (see polyscope/src/volume_mesh.cpp:43). +const HEX_FACES: &[FaceData] = &[ + FaceData { polygon: &[2, 1, 0, 3], triangulation: &[[2, 1, 0], [2, 0, 3]] }, // Bottom + FaceData { polygon: &[4, 0, 1, 5], triangulation: &[[4, 0, 1], [4, 1, 5]] }, // Front + FaceData { polygon: &[5, 1, 2, 6], triangulation: &[[5, 1, 2], [5, 2, 6]] }, // Right + FaceData { polygon: &[7, 3, 0, 4], triangulation: &[[7, 3, 0], [7, 0, 4]] }, // Left + FaceData { polygon: &[6, 2, 3, 7], triangulation: &[[6, 2, 3], [6, 3, 7]] }, // Back + FaceData { polygon: &[7, 4, 5, 6], triangulation: &[[7, 4, 5], [7, 5, 6]] }, // Top +]; + +// ===== Prism (wedge) ===== +// Slots 0,1,2 = bottom triangle; 3,4,5 = top triangle (slots 3,4,5 align with 0,1,2). +const PRISM_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, // Bottom tri + FaceData { polygon: &[0, 3, 5, 2], triangulation: &[[0, 5, 2], [0, 3, 5]] }, // Side quad 1 + FaceData { polygon: &[2, 5, 4, 1], triangulation: &[[2, 4, 5], [2, 1, 4]] }, // Side quad 2 + FaceData { polygon: &[0, 1, 4, 3], triangulation: &[[3, 0, 4], [0, 1, 4]] }, // Side quad 3 + FaceData { polygon: &[3, 4, 5], triangulation: &[[3, 4, 5]] }, // Top tri +]; + +// ===== Pyramid ===== +// Slots 0..3 = base quad (CCW from outside, looking from -apex toward base); +// Slot 4 = apex. +const PYRAMID_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 1, 2, 3], triangulation: &[[0, 3, 2], [0, 2, 1]] }, // Base quad + FaceData { polygon: &[0, 1, 4], triangulation: &[[0, 1, 4]] }, // Side 1 + FaceData { polygon: &[1, 2, 4], triangulation: &[[1, 2, 4]] }, // Side 2 + FaceData { polygon: &[2, 3, 4], triangulation: &[[2, 3, 4]] }, // Side 3 + FaceData { polygon: &[3, 0, 4], triangulation: &[[3, 0, 4]] }, // Side 4 +]; + +/// Returns the face data table for a given cell type. +#[must_use] +pub fn face_data_for(cell_type: VolumeCellType) -> &'static [FaceData] { + match cell_type { + VolumeCellType::Tet => TET_FACES, + VolumeCellType::Hex => HEX_FACES, + VolumeCellType::Prism => PRISM_FACES, + VolumeCellType::Pyramid => PYRAMID_FACES, + } +} + +/// Number of vertices used by a cell type (non-sentinel slots in the `[u32; 8]`). +#[must_use] +pub fn num_verts_in_cell(cell_type: VolumeCellType) -> usize { + match cell_type { + VolumeCellType::Tet => 4, + VolumeCellType::Hex => 8, + VolumeCellType::Prism => 6, + VolumeCellType::Pyramid => 5, + } +} + +/// Builds a canonical (sorted) face key for hashing. +/// +/// Face polygons can have 3 or 4 unique vertices. Triangular faces leave slot 3 +/// as `u32::MAX` so that triangle and quad keys never collide. +#[must_use] +pub fn canonical_face_key(cell: &[u32; 8], polygon: &[usize]) -> [u32; 4] { + let mut key = [u32::MAX; 4]; + debug_assert!(polygon.len() == 3 || polygon.len() == 4); + for (i, &slot) in polygon.iter().enumerate() { + key[i] = cell[slot]; + } + key.sort_unstable(); + key +} diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index c5cabe3..2bc43bb 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -41,6 +41,7 @@ //! mesh.add_vertex_scalar_quantity("temperature", vec![0.0, 0.5, 1.0, 0.25]); //! ``` +mod cell_data; mod color_quantity; mod scalar_quantity; pub mod slice_geometry; @@ -64,10 +65,14 @@ use polyscope_render::{ /// Cell type for volume meshes. #[derive(Debug, Clone, Copy, PartialEq, Eq)] pub enum VolumeCellType { - /// Tetrahedron (4 vertices) + /// Tetrahedron (4 vertices, 4 triangular faces) Tet, - /// Hexahedron (8 vertices) + /// Hexahedron (8 vertices, 6 quadrilateral faces) Hex, + /// Triangular prism / wedge (6 vertices, 2 tri + 3 quad faces) + Prism, + /// Square pyramid (5 vertices, 1 quad + 4 tri faces) + Pyramid, } /// A volume mesh structure (tetrahedral or hexahedral). @@ -1199,6 +1204,8 @@ impl VolumeMesh { std::array::from_fn(|i| self.vertices[cell[i] as usize]); slice_hex(hex_verts, plane_origin, plane_normal) } + // Prism/Pyramid slice geometry wired in Task 5 of upstream-port plan. + VolumeCellType::Prism | VolumeCellType::Pyramid => CellSliceResult::empty(), }; if slice.has_intersection() { @@ -1452,6 +1459,26 @@ const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ mod tests { use super::*; + #[test] + fn test_cell_type_enum_has_prism_and_pyramid() { + // Compile-time check: pattern match must be exhaustive over all 4 variants. + let types = [ + VolumeCellType::Tet, + VolumeCellType::Hex, + VolumeCellType::Prism, + VolumeCellType::Pyramid, + ]; + for t in types { + let label = match t { + VolumeCellType::Tet => "tet", + VolumeCellType::Hex => "hex", + VolumeCellType::Prism => "prism", + VolumeCellType::Pyramid => "pyramid", + }; + assert!(!label.is_empty()); + } + } + #[test] fn test_interior_face_detection() { // Two tets sharing a face From a6697d15bc99d3df54a0098d1181de3583e6d83d Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 09:50:59 +0200 Subject: [PATCH 02/14] feat(volume_mesh): use sentinel-count to classify cell type --- .../src/volume_mesh/mod.rs | 64 +++++++++++++++++-- 1 file changed, 60 insertions(+), 4 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 2bc43bb..7e1c2e5 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -192,12 +192,29 @@ impl VolumeMesh { } /// Returns the cell type of the given cell. + /// + /// Cell type is determined by the number of sentinel (`u32::MAX`) indices + /// in the 8-slot cell array (matches upstream C++ Polyscope): + /// - 0 sentinels → `Hex` (8 verts) + /// - 2 sentinels → `Prism` (6 verts) + /// - 3 sentinels → `Pyramid` (5 verts) + /// - 4 sentinels → `Tet` (4 verts) + /// + /// # Panics + /// Panics if `cell_idx` is out of range or the sentinel count is invalid + /// (1, 5, 6, 7, or 8 sentinels). #[must_use] pub fn cell_type(&self, cell_idx: usize) -> VolumeCellType { - if self.cells[cell_idx][4] == u32::MAX { - VolumeCellType::Tet - } else { - VolumeCellType::Hex + let sentinels = self.cells[cell_idx] + .iter() + .filter(|&&v| v == u32::MAX) + .count(); + match sentinels { + 0 => VolumeCellType::Hex, + 2 => VolumeCellType::Prism, + 3 => VolumeCellType::Pyramid, + 4 => VolumeCellType::Tet, + n => panic!("VolumeMesh cell {cell_idx}: invalid sentinel count {n} (expected 0/2/3/4)"), } } @@ -1459,6 +1476,45 @@ const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ mod tests { use super::*; + #[test] + fn test_cell_type_detection_tet() { + let mesh = VolumeMesh::new( + "t", + vec![Vec3::ZERO, Vec3::X, Vec3::Y, Vec3::Z], + vec![[0, 1, 2, 3, u32::MAX, u32::MAX, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Tet); + } + + #[test] + fn test_cell_type_detection_hex() { + let verts = (0..8).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new("h", verts, vec![[0, 1, 2, 3, 4, 5, 6, 7]]); + assert_eq!(mesh.cell_type(0), VolumeCellType::Hex); + } + + #[test] + fn test_cell_type_detection_prism() { + let verts = (0..6).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new( + "p", + verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Prism); + } + + #[test] + fn test_cell_type_detection_pyramid() { + let verts = (0..5).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new( + "py", + verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Pyramid); + } + #[test] fn test_cell_type_enum_has_prism_and_pyramid() { // Compile-time check: pattern match must be exhaustive over all 4 variants. From 57fb27de73964bb4d25ef92238c8fcecda710792 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 09:52:51 +0200 Subject: [PATCH 03/14] feat(volume_mesh): tet-decomposition + centroid for prism/pyramid --- .../src/volume_mesh/cell_data.rs | 104 ++++++++++++++ .../src/volume_mesh/mod.rs | 129 ++++++++++++------ 2 files changed, 190 insertions(+), 43 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index c0c4f16..e7fe483 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -107,3 +107,107 @@ pub fn canonical_face_key(cell: &[u32; 8], polygon: &[usize]) -> [u32; 4] { key.sort_unstable(); key } + +/// Decomposes a hex cell into 5 tetrahedra using a fixed diagonal pattern. +/// +/// Central-diagonal pattern from Dompierre et al.; the 5-tet split works for +/// any convex hex. Matches the previous polyscope-rs decomposition. +const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ + [0, 1, 2, 5], + [0, 2, 7, 5], + [0, 2, 3, 7], + [0, 5, 7, 4], + [2, 7, 5, 6], +]; + +/// Decomposes a triangular prism into 3 tetrahedra. +/// +/// Picks a consistent diagonal split on the quad face opposite the +/// lowest-numbered vertex, matching the upstream algorithm (`decomposePrism` in +/// `polyscope/src/volume_mesh.cpp`). Consistency across adjacent cells matters +/// so that shared faces are tessellated identically and no gaps appear in +/// slice caps for mixed meshes. +#[must_use] +pub fn decompose_prism(cell: &[u32; 8]) -> [[u32; 4]; 3] { + let mut p: [u32; 6] = [cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]]; + + let min_idx = (0..6).min_by_key(|&i| p[i]).unwrap(); + + if min_idx < 3 { + let rot = match min_idx { + 0 => 0, + 1 => 2, + _ => 1, + }; + rotate_prism_in_place(&mut p, rot); + } else { + let top_pos = min_idx - 3; + let rot = match top_pos { + 0 => 0, + 1 => 2, + _ => 1, + }; + p.swap(0, 3); + p.swap(1, 4); + p.swap(2, 5); + rotate_prism_in_place(&mut p, rot); + } + + if p[2].min(p[4]) < p[1].min(p[5]) { + [ + [p[0], p[5], p[4], p[3]], + [p[0], p[4], p[5], p[2]], + [p[0], p[4], p[2], p[1]], + ] + } else { + [ + [p[0], p[5], p[4], p[3]], + [p[0], p[1], p[5], p[2]], + [p[0], p[5], p[1], p[4]], + ] + } +} + +fn rotate_prism_in_place(p: &mut [u32; 6], rot: usize) { + const BOTTOM_ROT: [[usize; 3]; 3] = [[0, 1, 2], [1, 2, 0], [2, 0, 1]]; + const TOP_ROT: [[usize; 3]; 3] = [[3, 4, 5], [4, 5, 3], [5, 3, 4]]; + let src = *p; + for i in 0..3 { + p[i] = src[BOTTOM_ROT[rot][i]]; + p[i + 3] = src[TOP_ROT[rot][i]]; + } +} + +/// Decomposes a square pyramid into 2 tetrahedra by splitting the base quad +/// along the diagonal containing the smaller of {p[0], p[2]} vs {p[1], p[3]}. +/// Consistent split ensures adjacent cells tessellate the shared face the same way. +#[must_use] +pub fn decompose_pyramid(cell: &[u32; 8]) -> [[u32; 4]; 2] { + let p: [u32; 5] = [cell[0], cell[1], cell[2], cell[3], cell[4]]; + + if p[0].min(p[2]) < p[1].min(p[3]) { + [ + [p[0], p[2], p[4], p[1]], + [p[0], p[4], p[2], p[3]], + ] + } else { + [ + [p[1], p[3], p[4], p[2]], + [p[1], p[4], p[3], p[0]], + ] + } +} + +/// Returns the tet-decomposition for any cell, dispatching on cell type. +#[must_use] +pub fn decompose_cell_to_tets(cell: &[u32; 8], cell_type: VolumeCellType) -> Vec<[u32; 4]> { + match cell_type { + VolumeCellType::Tet => vec![[cell[0], cell[1], cell[2], cell[3]]], + VolumeCellType::Hex => HEX_TO_TET_PATTERN + .iter() + .map(|t| [cell[t[0]], cell[t[1]], cell[t[2]], cell[t[3]]]) + .collect(), + VolumeCellType::Prism => decompose_prism(cell).to_vec(), + VolumeCellType::Pyramid => decompose_pyramid(cell).to_vec(), + } +} diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 7e1c2e5..6f54808 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -279,29 +279,19 @@ impl VolumeMesh { } /// Decomposes all cells into tetrahedra. - /// Tets pass through unchanged, hexes are decomposed into 5 tets. + /// + /// Decomposition counts per cell type: + /// - Tet → 1 tet (passthrough) + /// - Hex → 5 tets (fixed diagonal pattern) + /// - Prism → 3 tets (consistent diagonal split) + /// - Pyramid → 2 tets (consistent diagonal split) #[must_use] pub fn decompose_to_tets(&self) -> Vec<[u32; 4]> { let mut tets = Vec::new(); - - for cell in &self.cells { - if cell[4] == u32::MAX { - // Already a tet - tets.push([cell[0], cell[1], cell[2], cell[3]]); - } else { - // Hex - decompose using diagonal pattern (5 tets) - for tet_local in &HEX_TO_TET_PATTERN { - let tet = [ - cell[tet_local[0]], - cell[tet_local[1]], - cell[tet_local[2]], - cell[tet_local[3]], - ]; - tets.push(tet); - } - } + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + tets.extend(cell_data::decompose_cell_to_tets(cell, ct)); } - tets } @@ -339,22 +329,17 @@ impl VolumeMesh { face_counts } - /// Computes the centroid of a cell. - fn cell_centroid(&self, cell: &[u32; 8]) -> Vec3 { - if cell[4] == u32::MAX { - // Tetrahedron: average of 4 vertices - let sum = self.vertices[cell[0] as usize] - + self.vertices[cell[1] as usize] - + self.vertices[cell[2] as usize] - + self.vertices[cell[3] as usize]; - sum / 4.0 - } else { - // Hexahedron: average of 8 vertices - let sum = (0..8) - .map(|i| self.vertices[cell[i] as usize]) - .fold(Vec3::ZERO, |a, b| a + b); - sum / 8.0 + /// Computes the centroid of a cell as the mean of its real (non-sentinel) vertices. + pub(crate) fn cell_centroid(&self, cell: &[u32; 8]) -> Vec3 { + let mut sum = Vec3::ZERO; + let mut count = 0u32; + for &v in cell { + if v != u32::MAX { + sum += self.vertices[v as usize]; + count += 1; + } } + sum / count as f32 } /// Tests if a cell should be visible based on slice planes. @@ -1463,19 +1448,77 @@ pub struct VolumeMeshRenderGeometry { pub vertex_colors: Option>, } -/// Diagonal decomposition patterns (5 tets). -const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ - [0, 1, 2, 5], - [0, 2, 7, 5], - [0, 2, 3, 7], - [0, 5, 7, 4], - [2, 7, 5, 6], -]; - #[cfg(test)] mod tests { use super::*; + #[test] + fn test_decompose_prism_to_tets() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let mesh = VolumeMesh::new( + "prism_only", + verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + let tets = mesh.decompose_to_tets(); + assert_eq!(tets.len(), 3, "prism should decompose to 3 tets"); + for tet in &tets { + for &v in tet { + assert!(v < 6, "tet vertex index {v} out of range for prism"); + } + } + } + + #[test] + fn test_decompose_pyramid_to_tets() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let mesh = VolumeMesh::new( + "pyr_only", + verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + let tets = mesh.decompose_to_tets(); + assert_eq!(tets.len(), 2, "pyramid should decompose to 2 tets"); + for tet in &tets { + for &v in tet { + assert!(v < 5, "tet vertex index {v} out of range for pyramid"); + } + } + } + + #[test] + fn test_cell_centroid_prism() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(2.0, 0.0, 0.0), + Vec3::new(1.0, 2.0, 0.0), + Vec3::new(0.0, 0.0, 4.0), + Vec3::new(2.0, 0.0, 4.0), + Vec3::new(1.0, 2.0, 4.0), + ]; + let mesh = VolumeMesh::new( + "p", + verts.clone(), + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + let expected: Vec3 = verts.iter().copied().sum::() / 6.0; + let centroid = mesh.cell_centroid(&mesh.cells()[0]); + assert!((centroid - expected).length() < 1e-5); + } + #[test] fn test_cell_type_detection_tet() { let mesh = VolumeMesh::new( From 9010c0bb47c18c6955b78e8ce9947e57ce5c9b66 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 09:55:35 +0200 Subject: [PATCH 04/14] refactor(volume_mesh): dispatch face iteration via cell_data table --- .../src/volume_mesh/cell_data.rs | 15 - .../src/volume_mesh/mod.rs | 345 +++++++----------- 2 files changed, 135 insertions(+), 225 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index e7fe483..d44e38b 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -1,7 +1,3 @@ -// Items are wired in by subsequent tasks (tet decomposition + face dispatch). -// The allow is removed at that point. -#![allow(dead_code)] - //! Static shape data for each `VolumeCellType`. //! //! Defines, for each cell type: @@ -82,17 +78,6 @@ pub fn face_data_for(cell_type: VolumeCellType) -> &'static [FaceData] { } } -/// Number of vertices used by a cell type (non-sentinel slots in the `[u32; 8]`). -#[must_use] -pub fn num_verts_in_cell(cell_type: VolumeCellType) -> usize { - match cell_type { - VolumeCellType::Tet => 4, - VolumeCellType::Hex => 8, - VolumeCellType::Prism => 6, - VolumeCellType::Pyramid => 5, - } -} - /// Builds a canonical (sorted) face key for hashing. /// /// Face polygons can have 3 or 4 unique vertices. Triangular faces leave slot 3 diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 6f54808..22a42ae 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -304,28 +304,13 @@ impl VolumeMesh { /// Computes face counts for interior/exterior detection. fn compute_face_counts(&self) -> HashMap<[u32; 4], usize> { let mut face_counts: HashMap<[u32; 4], usize> = HashMap::new(); - - for cell in &self.cells { - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - *face_counts.entry(key).or_insert(0) += 1; - } - } else { - // Hexahedron - each quad face uses same 4 vertices - for quad in HEX_FACE_STENCIL { - // Get the 4 unique vertices of this quad face - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; // The fourth vertex - let key = canonical_face_key(v0, v1, v2, Some(v3)); - *face_counts.entry(key).or_insert(0) += 1; - } + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + *face_counts.entry(key).or_insert(0) += 1; } } - face_counts } @@ -360,38 +345,22 @@ impl VolumeMesh { true } - /// Computes face counts for interior/exterior detection, only for visible cells. + /// Same as `compute_face_counts` but skips cells culled by slice planes. fn compute_face_counts_with_culling( &self, planes: &[(Vec3, Vec3)], ) -> HashMap<[u32; 4], usize> { let mut face_counts: HashMap<[u32; 4], usize> = HashMap::new(); - - for cell in &self.cells { - // Skip cells culled by slice planes + for (cell_idx, cell) in self.cells.iter().enumerate() { if !self.is_cell_visible(cell, planes) { continue; } - - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - *face_counts.entry(key).or_insert(0) += 1; - } - } else { - // Hexahedron - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - *face_counts.entry(key).or_insert(0) += 1; - } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + *face_counts.entry(key).or_insert(0) += 1; } } - face_counts } @@ -401,38 +370,19 @@ impl VolumeMesh { let mut positions = Vec::new(); let mut faces = Vec::new(); - for cell in &self.cells { - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - if face_counts[&key] == 1 { - // Exterior face - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts[&key] != 1 { + continue; } - } else { - // Hexahedron - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - if face_counts[&key] == 1 { - // Exterior face - emit both triangles - for [a, b, c] in quad { - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } - } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + faces.push([base_idx, base_idx + 1, base_idx + 2]); } } } @@ -441,52 +391,30 @@ impl VolumeMesh { } /// Generates triangulated exterior faces with cell culling based on slice planes. - /// Only cells whose centroid is on the positive side of all planes are rendered. fn generate_render_geometry_with_culling( &self, planes: &[(Vec3, Vec3)], ) -> (Vec, Vec<[u32; 3]>) { - // Compute face counts only for visible cells let face_counts = self.compute_face_counts_with_culling(planes); let mut positions = Vec::new(); let mut faces = Vec::new(); - for cell in &self.cells { - // Skip cells culled by slice planes + for (cell_idx, cell) in self.cells.iter().enumerate() { if !self.is_cell_visible(cell, planes) { continue; } - - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - // Render face if it's exterior among visible cells - if face_counts.get(&key) == Some(&1) { - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; } - } else { - // Hexahedron - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - if face_counts.get(&key) == Some(&1) { - for [a, b, c] in quad { - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } - } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + faces.push([base_idx, base_idx + 1, base_idx + 2]); } } } @@ -505,47 +433,24 @@ impl VolumeMesh { // First pass: generate geometry and track indices for (cell_idx, cell) in self.cells.iter().enumerate() { - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - if face_counts[&key] == 1 { - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - vertex_indices.push(cell[a] as usize); - vertex_indices.push(cell[b] as usize); - vertex_indices.push(cell[c] as usize); - cell_indices.push(cell_idx); - cell_indices.push(cell_idx); - cell_indices.push(cell_idx); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts[&key] != 1 { + continue; } - } else { - // Hexahedron - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - if face_counts[&key] == 1 { - for [a, b, c] in quad { - let base_idx = positions.len() as u32; - positions.push(self.vertices[cell[a] as usize]); - positions.push(self.vertices[cell[b] as usize]); - positions.push(self.vertices[cell[c] as usize]); - vertex_indices.push(cell[a] as usize); - vertex_indices.push(cell[b] as usize); - vertex_indices.push(cell[c] as usize); - cell_indices.push(cell_idx); - cell_indices.push(cell_idx); - cell_indices.push(cell_idx); - faces.push([base_idx, base_idx + 1, base_idx + 2]); - } - } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + vertex_indices.push(cell[a] as usize); + vertex_indices.push(cell[b] as usize); + vertex_indices.push(cell[c] as usize); + cell_indices.push(cell_idx); + cell_indices.push(cell_idx); + cell_indices.push(cell_idx); + faces.push([base_idx, base_idx + 1, base_idx + 2]); } } } @@ -770,33 +675,18 @@ impl VolumeMesh { fn generate_cell_index_per_triangle(&self) -> Vec { let face_counts = self.compute_face_counts(); let mut cell_indices = Vec::new(); - for (cell_idx, cell) in self.cells.iter().enumerate() { - if cell[4] == u32::MAX { - // Tetrahedron - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - if face_counts[&key] == 1 { - cell_indices.push(cell_idx as u32); - } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; } - } else { - // Hexahedron - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - if face_counts[&key] == 1 { - // 2 triangles per quad face - cell_indices.push(cell_idx as u32); - cell_indices.push(cell_idx as u32); - } + for _tri in face.triangulation { + cell_indices.push(cell_idx as u32); } } } - cell_indices } @@ -804,34 +694,21 @@ impl VolumeMesh { fn generate_cell_index_per_triangle_with_culling(&self, planes: &[(Vec3, Vec3)]) -> Vec { let face_counts = self.compute_face_counts_with_culling(planes); let mut cell_indices = Vec::new(); - for (cell_idx, cell) in self.cells.iter().enumerate() { if !self.is_cell_visible(cell, planes) { continue; } - - if cell[4] == u32::MAX { - for [a, b, c] in TET_FACE_STENCIL { - let key = canonical_face_key(cell[a], cell[b], cell[c], None); - if face_counts.get(&key) == Some(&1) { - cell_indices.push(cell_idx as u32); - } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; } - } else { - for quad in HEX_FACE_STENCIL { - let v0 = cell[quad[0][0]]; - let v1 = cell[quad[0][1]]; - let v2 = cell[quad[0][2]]; - let v3 = cell[quad[1][2]]; - let key = canonical_face_key(v0, v1, v2, Some(v3)); - if face_counts.get(&key) == Some(&1) { - cell_indices.push(cell_idx as u32); - cell_indices.push(cell_idx as u32); - } + for _tri in face.triangulation { + cell_indices.push(cell_idx as u32); } } } - cell_indices } @@ -1416,27 +1293,6 @@ impl HasQuantities for VolumeMesh { use std::collections::HashMap; -/// Generates a canonical (sorted) face key for hashing. -/// For triangular faces, the fourth element is `u32::MAX`. -fn canonical_face_key(v0: u32, v1: u32, v2: u32, v3: Option) -> [u32; 4] { - let mut key = [v0, v1, v2, v3.unwrap_or(u32::MAX)]; - key.sort_unstable(); - key -} - -/// Face stencil for tetrahedra: 4 triangular faces -const TET_FACE_STENCIL: [[usize; 3]; 4] = [[0, 2, 1], [0, 1, 3], [0, 3, 2], [1, 2, 3]]; - -/// Face stencil for hexahedra: 6 quad faces (each as 2 triangles sharing diagonal) -const HEX_FACE_STENCIL: [[[usize; 3]; 2]; 6] = [ - [[2, 1, 0], [2, 0, 3]], // Bottom - [[4, 0, 1], [4, 1, 5]], // Front - [[5, 1, 2], [5, 2, 6]], // Right - [[7, 3, 0], [7, 0, 4]], // Left - [[6, 2, 3], [6, 3, 7]], // Back - [[7, 4, 5], [7, 5, 6]], // Top -]; - /// Render geometry data with optional quantity values. pub struct VolumeMeshRenderGeometry { pub positions: Vec, @@ -1452,6 +1308,75 @@ pub struct VolumeMeshRenderGeometry { mod tests { use super::*; + #[test] + fn test_single_prism_all_exterior() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let mesh = VolumeMesh::new( + "p", + verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + let (_, faces) = mesh.generate_render_geometry(); + assert_eq!(faces.len(), 8, "single prism should have 8 triangles (2 tri + 3*2 quad)"); + } + + #[test] + fn test_single_pyramid_all_exterior() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let mesh = VolumeMesh::new( + "py", + verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + let (_, faces) = mesh.generate_render_geometry(); + assert_eq!(faces.len(), 6, "single pyramid should have 6 triangles (2 base + 4 sides)"); + } + + #[test] + fn test_mixed_cell_mesh_renders_all() { + let verts = vec![ + // Tet + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + // Prism (translated +5 in x) + Vec3::new(5.0, 0.0, 0.0), + Vec3::new(6.0, 0.0, 0.0), + Vec3::new(5.5, 1.0, 0.0), + Vec3::new(5.0, 0.0, 1.0), + Vec3::new(6.0, 0.0, 1.0), + Vec3::new(5.5, 1.0, 1.0), + // Pyramid (translated +10 in x) + Vec3::new(10.0, 0.0, 0.0), + Vec3::new(11.0, 0.0, 0.0), + Vec3::new(11.0, 1.0, 0.0), + Vec3::new(10.0, 1.0, 0.0), + Vec3::new(10.5, 0.5, 1.0), + ]; + let cells = vec![ + [0, 1, 2, 3, u32::MAX, u32::MAX, u32::MAX, u32::MAX], + [4, 5, 6, 7, 8, 9, u32::MAX, u32::MAX], + [10, 11, 12, 13, 14, u32::MAX, u32::MAX, u32::MAX], + ]; + let mesh = VolumeMesh::new("m", verts, cells); + let (_, faces) = mesh.generate_render_geometry(); + assert_eq!(faces.len(), 4 + 8 + 6, "mixed mesh should sum per-cell triangle counts"); + } + #[test] fn test_decompose_prism_to_tets() { let verts = vec![ From fd2256f3f5415c7e73c44e171b3a5c28478cbf65 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 09:57:01 +0200 Subject: [PATCH 05/14] feat(volume_mesh): slice_prism / slice_pyramid via tet decomposition --- .../src/volume_mesh/mod.rs | 14 +- .../src/volume_mesh/slice_geometry.rs | 165 ++++++++++++++++++ 2 files changed, 176 insertions(+), 3 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 22a42ae..a6dfa98 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -49,7 +49,7 @@ mod vector_quantity; pub use color_quantity::*; pub use scalar_quantity::*; -pub use slice_geometry::{CellSliceResult, slice_hex, slice_tet}; +pub use slice_geometry::{CellSliceResult, slice_hex, slice_prism, slice_pyramid, slice_tet}; pub use vector_quantity::*; // Re-export SliceMeshData from this module @@ -1083,8 +1083,16 @@ impl VolumeMesh { std::array::from_fn(|i| self.vertices[cell[i] as usize]); slice_hex(hex_verts, plane_origin, plane_normal) } - // Prism/Pyramid slice geometry wired in Task 5 of upstream-port plan. - VolumeCellType::Prism | VolumeCellType::Pyramid => CellSliceResult::empty(), + VolumeCellType::Prism => { + let prism_verts: [Vec3; 6] = + std::array::from_fn(|i| self.vertices[cell[i] as usize]); + slice_prism(prism_verts, plane_origin, plane_normal) + } + VolumeCellType::Pyramid => { + let pyr_verts: [Vec3; 5] = + std::array::from_fn(|i| self.vertices[cell[i] as usize]); + slice_pyramid(pyr_verts, plane_origin, plane_normal) + } }; if slice.has_intersection() { diff --git a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs index 71fa174..e4453ad 100644 --- a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs +++ b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs @@ -148,6 +148,105 @@ pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> } } +/// Slice a triangular prism by decomposing into 3 tetrahedra. +/// +/// # Arguments +/// * `vertices` - The 6 vertices of the prism (slots 0..2 = bottom tri, 3..5 = top tri) +/// * `plane_origin` - A point on the plane +/// * `plane_normal` - The plane normal (points toward kept geometry) +#[must_use] +pub fn slice_prism( + vertices: [Vec3; 6], + plane_origin: Vec3, + plane_normal: Vec3, +) -> CellSliceResult { + // Symmetric 3-tet decomposition. Cross-cell consistency is not required for + // isolated slicing — the choice doesn't affect correctness. + let tet_indices = [ + [0usize, 5, 4, 3], + [0, 4, 5, 2], + [0, 4, 2, 1], + ]; + + let mut all_vertices = Vec::new(); + let mut all_interp = Vec::new(); + + for tet in &tet_indices { + let r = slice_tet( + vertices[tet[0]], + vertices[tet[1]], + vertices[tet[2]], + vertices[tet[3]], + plane_origin, + plane_normal, + ); + for (local_a, local_b, t) in r.interpolation { + let pa = tet[local_a as usize] as u32; + let pb = tet[local_b as usize] as u32; + all_interp.push((pa, pb, t)); + } + all_vertices.extend(r.vertices); + } + + merge_slice_vertices(&mut all_vertices, &mut all_interp); + if all_vertices.len() >= 3 { + order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); + } + + CellSliceResult { + vertices: all_vertices, + interpolation: all_interp, + } +} + +/// Slice a square pyramid by decomposing into 2 tetrahedra. +/// +/// # Arguments +/// * `vertices` - The 5 vertices of the pyramid (slots 0..3 = base quad, 4 = apex) +/// * `plane_origin` - A point on the plane +/// * `plane_normal` - The plane normal (points toward kept geometry) +#[must_use] +pub fn slice_pyramid( + vertices: [Vec3; 5], + plane_origin: Vec3, + plane_normal: Vec3, +) -> CellSliceResult { + let tet_indices = [ + [0usize, 2, 4, 1], + [0, 4, 2, 3], + ]; + + let mut all_vertices = Vec::new(); + let mut all_interp = Vec::new(); + + for tet in &tet_indices { + let r = slice_tet( + vertices[tet[0]], + vertices[tet[1]], + vertices[tet[2]], + vertices[tet[3]], + plane_origin, + plane_normal, + ); + for (local_a, local_b, t) in r.interpolation { + let pa = tet[local_a as usize] as u32; + let pb = tet[local_b as usize] as u32; + all_interp.push((pa, pb, t)); + } + all_vertices.extend(r.vertices); + } + + merge_slice_vertices(&mut all_vertices, &mut all_interp); + if all_vertices.len() >= 3 { + order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); + } + + CellSliceResult { + vertices: all_vertices, + interpolation: all_interp, + } +} + /// Orders polygon vertices in counter-clockwise order around the centroid. /// /// This ensures the resulting polygon is suitable for rendering with correct face winding. @@ -376,6 +475,72 @@ mod tests { } } + #[test] + fn test_slice_prism_through_middle() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let result = slice_prism(verts, Vec3::new(0.0, 0.0, 0.5), Vec3::Z); + assert!(result.has_intersection()); + assert!( + result.vertices.len() >= 3, + "expected at least 3 verts, got {}", + result.vertices.len() + ); + for v in &result.vertices { + assert!((v.z - 0.5).abs() < 1e-4, "vertex z={} should be 0.5", v.z); + } + } + + #[test] + fn test_slice_prism_no_intersection() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let result = slice_prism(verts, Vec3::new(0.0, 0.0, 2.0), Vec3::Z); + assert!(!result.has_intersection()); + } + + #[test] + fn test_slice_pyramid_through_middle() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let result = slice_pyramid(verts, Vec3::new(0.0, 0.0, 0.5), Vec3::Z); + assert!(result.has_intersection()); + assert!(result.vertices.len() >= 3); + for v in &result.vertices { + assert!((v.z - 0.5).abs() < 1e-4); + } + } + + #[test] + fn test_slice_pyramid_no_intersection() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let result = slice_pyramid(verts, Vec3::new(0.0, 0.0, -1.0), Vec3::Z); + assert!(!result.has_intersection()); + } + #[test] fn test_polygon_ordering() { // Verify vertices are ordered correctly (counter-clockwise) From 8937f2b73a312169c75d6af5f31b07bb08305739 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 10:00:31 +0200 Subject: [PATCH 06/14] feat(volume_mesh): public register_prism_mesh / register_pyramid_mesh --- .../src/volume_mesh/mod.rs | 66 +++++++++++++++++++ crates/polyscope/src/volume_mesh.rs | 44 +++++++++++++ 2 files changed, 110 insertions(+) diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index a6dfa98..2b1094a 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -179,6 +179,41 @@ impl VolumeMesh { Self::new(name, vertices, hexes) } + /// Creates a triangular-prism (wedge) mesh. + /// + /// Each prism has 6 vertices: slots 0..2 form the bottom triangle and + /// slots 3..5 form the top triangle (slot `i+3` should be the vertex + /// above slot `i`). Cells are stored as 8-index arrays with the last two + /// slots set to `u32::MAX` for sentinel detection. + pub fn new_prism_mesh( + name: impl Into, + vertices: Vec, + prisms: Vec<[u32; 6]>, + ) -> Self { + let cells: Vec<[u32; 8]> = prisms + .into_iter() + .map(|p| [p[0], p[1], p[2], p[3], p[4], p[5], u32::MAX, u32::MAX]) + .collect(); + Self::new(name, vertices, cells) + } + + /// Creates a square-pyramid mesh. + /// + /// Each pyramid has 5 vertices: slots 0..3 form the base quad (in CCW + /// order viewed from outside the cell) and slot 4 is the apex. Cells are + /// stored as 8-index arrays with the last three slots set to `u32::MAX`. + pub fn new_pyramid_mesh( + name: impl Into, + vertices: Vec, + pyramids: Vec<[u32; 5]>, + ) -> Self { + let cells: Vec<[u32; 8]> = pyramids + .into_iter() + .map(|p| [p[0], p[1], p[2], p[3], p[4], u32::MAX, u32::MAX, u32::MAX]) + .collect(); + Self::new(name, vertices, cells) + } + /// Returns the number of vertices. #[must_use] pub fn num_vertices(&self) -> usize { @@ -1316,6 +1351,37 @@ pub struct VolumeMeshRenderGeometry { mod tests { use super::*; + #[test] + fn test_new_prism_mesh_constructor() { + let verts = vec![ + Vec3::ZERO, + Vec3::X, + Vec3::Y, + Vec3::Z, + Vec3::X + Vec3::Z, + Vec3::Y + Vec3::Z, + ]; + let prisms = vec![[0u32, 1, 2, 3, 4, 5]]; + let mesh = VolumeMesh::new_prism_mesh("p", verts, prisms); + assert_eq!(mesh.num_cells(), 1); + assert_eq!(mesh.cell_type(0), VolumeCellType::Prism); + } + + #[test] + fn test_new_pyramid_mesh_constructor() { + let verts = vec![ + Vec3::ZERO, + Vec3::X, + Vec3::X + Vec3::Y, + Vec3::Y, + Vec3::splat(0.5) + Vec3::Z, + ]; + let pyramids = vec![[0u32, 1, 2, 3, 4]]; + let mesh = VolumeMesh::new_pyramid_mesh("py", verts, pyramids); + assert_eq!(mesh.num_cells(), 1); + assert_eq!(mesh.cell_type(0), VolumeCellType::Pyramid); + } + #[test] fn test_single_prism_all_exterior() { let verts = vec![ diff --git a/crates/polyscope/src/volume_mesh.rs b/crates/polyscope/src/volume_mesh.rs index 72780d6..1a2a15f 100644 --- a/crates/polyscope/src/volume_mesh.rs +++ b/crates/polyscope/src/volume_mesh.rs @@ -68,6 +68,50 @@ pub fn register_hex_mesh( VolumeMeshHandle { name } } +/// Registers a triangular-prism (wedge) mesh with polyscope. +/// +/// Each entry in `prisms` is 6 vertex indices: slots 0..2 = bottom triangle, +/// slots 3..5 = top triangle (with slot `i+3` directly above slot `i`). +pub fn register_prism_mesh( + name: impl Into, + vertices: Vec, + prisms: Vec<[u32; 6]>, +) -> VolumeMeshHandle { + let name = name.into(); + let mesh = VolumeMesh::new_prism_mesh(name.clone(), vertices, prisms); + + with_context_mut(|ctx| { + ctx.registry + .register(Box::new(mesh)) + .expect("failed to register prism mesh"); + ctx.update_extents(); + }); + + VolumeMeshHandle { name } +} + +/// Registers a square-pyramid mesh with polyscope. +/// +/// Each entry in `pyramids` is 5 vertex indices: slots 0..3 = base quad (CCW +/// from outside), slot 4 = apex. +pub fn register_pyramid_mesh( + name: impl Into, + vertices: Vec, + pyramids: Vec<[u32; 5]>, +) -> VolumeMeshHandle { + let name = name.into(); + let mesh = VolumeMesh::new_pyramid_mesh(name.clone(), vertices, pyramids); + + with_context_mut(|ctx| { + ctx.registry + .register(Box::new(mesh)) + .expect("failed to register pyramid mesh"); + ctx.update_extents(); + }); + + VolumeMeshHandle { name } +} + /// Registers a generic volume mesh with polyscope. /// /// Cells are stored as 8-index arrays. For tetrahedra, indices 4-7 should be `u32::MAX`. From 7cee690f41fe11ece2e15e1633ad6a1f00f9a2a1 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 10:01:59 +0200 Subject: [PATCH 07/14] docs(examples): add mixed-cell-type volume mesh demo --- crates/polyscope/Cargo.toml | 4 ++ examples/volume_mesh_mixed_cells_demo.rs | 69 ++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 examples/volume_mesh_mixed_cells_demo.rs diff --git a/crates/polyscope/Cargo.toml b/crates/polyscope/Cargo.toml index 8b76708..cd260e3 100644 --- a/crates/polyscope/Cargo.toml +++ b/crates/polyscope/Cargo.toml @@ -52,6 +52,10 @@ path = "../../examples/volume_grid_demo.rs" name = "volume_mesh_demo" path = "../../examples/volume_mesh_demo.rs" +[[example]] +name = "volume_mesh_mixed_cells_demo" +path = "../../examples/volume_mesh_mixed_cells_demo.rs" + [[example]] name = "slice_plane_demo" path = "../../examples/slice_plane_demo.rs" diff --git a/examples/volume_mesh_mixed_cells_demo.rs b/examples/volume_mesh_mixed_cells_demo.rs new file mode 100644 index 0000000..09024c3 --- /dev/null +++ b/examples/volume_mesh_mixed_cells_demo.rs @@ -0,0 +1,69 @@ +//! Demonstrates all four volume-mesh cell types side by side: tet, hex, +//! prism (wedge), and pyramid. Each cell is rendered as its own mesh so +//! the structure list shows them separately. +//! +//! Run with: cargo run --example `volume_mesh_mixed_cells_demo` + +use glam::Vec3; +use polyscope_rs::{ + Result, init, register_hex_mesh, register_prism_mesh, register_pyramid_mesh, + register_tet_mesh, show, +}; + +fn main() -> Result<()> { + env_logger::init(); + init()?; + + // ----- Tet ----- + let tet_verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let tet = register_tet_mesh("tet", tet_verts, vec![[0, 1, 2, 3]]); + tet.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 1.0]); + + // ----- Hex (translated +3 in x) ----- + let dx = Vec3::new(3.0, 0.0, 0.0); + let hex_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 1.0, 1.0) + dx, + Vec3::new(0.0, 1.0, 1.0) + dx, + ]; + let hex = register_hex_mesh("hex", hex_verts, vec![[0, 1, 2, 3, 4, 5, 6, 7]]); + hex.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]); + + // ----- Prism (translated +6 in x) ----- + let dx = Vec3::new(6.0, 0.0, 0.0); + let prism_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(0.5, 1.0, 0.0) + dx, + Vec3::new(0.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 0.0, 1.0) + dx, + Vec3::new(0.5, 1.0, 1.0) + dx, + ]; + let prism = register_prism_mesh("prism", prism_verts, vec![[0, 1, 2, 3, 4, 5]]); + prism.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0]); + + // ----- Pyramid (translated +9 in x) ----- + let dx = Vec3::new(9.0, 0.0, 0.0); + let pyr_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 1.0, 0.0) + dx, + Vec3::new(0.5, 0.5, 1.0) + dx, + ]; + let pyr = register_pyramid_mesh("pyramid", pyr_verts, vec![[0, 1, 2, 3, 4]]); + pyr.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 0.0, 1.0]); + + show(); + Ok(()) +} From be474a0d4afe075190c77bb5a93dc787855bffb3 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 10:02:36 +0200 Subject: [PATCH 08/14] docs: prism/pyramid cell support (upstream PR #353) --- CHANGELOG.md | 19 +++++++++++++++++++ docs/feature-status.md | 3 ++- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 628a6ce..9114ceb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,25 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added +- Volume Mesh now supports prism (wedge) and pyramid cell types in addition to + tetrahedra and hexahedra. New constructors: `register_prism_mesh`, + `register_pyramid_mesh`, `VolumeMesh::new_prism_mesh`, + `VolumeMesh::new_pyramid_mesh`. Mixed-cell meshes are supported by using the + 8-slot cell array with sentinel `u32::MAX` indices in unused slots + (matches upstream Polyscope PR #353). +- `slice_prism` and `slice_pyramid` helpers in + `polyscope_structures::volume_mesh::slice_geometry`. +- Example `volume_mesh_mixed_cells_demo` showing all four cell types. + +### Changed +- `VolumeMesh::cell_type` now classifies cells by sentinel count (0/2/3/4 → + Hex/Prism/Pyramid/Tet) instead of by `cell[4] == u32::MAX`. Mixed meshes + produced by external pipelines must place sentinels in the correct trailing + slots; tet meshes built with `new_tet_mesh` continue to work unchanged. + ## [0.5.9] - 2026-03-02 ### Changed diff --git a/docs/feature-status.md b/docs/feature-status.md index e22eb7f..91588f0 100644 --- a/docs/feature-status.md +++ b/docs/feature-status.md @@ -9,7 +9,7 @@ Feature parity tracking between polyscope-rs and C++ Polyscope 2.x. | Point Cloud | Full | Full | Sphere impostors via instanced rendering | | Surface Mesh | Full | Full | Triangles + arbitrary polygons, full quantity support | | Curve Network | Full | Full | Lines + tubes via compute shaders | -| Volume Mesh | Full | Full | Tet/hex, interior face detection, slice capping | +| Volume Mesh | Full | Full | Tet/hex/prism/pyramid, interior face detection, slice capping | | Volume Grid | Full | Full | Node/cell scalars, gridcube + isosurface (marching cubes) | | Camera View | Full | Full | Frustum visualization | | Floating Quantities | Full | Full | Scalar/color images, depth/color/raw render images | @@ -77,6 +77,7 @@ Feature parity tracking between polyscope-rs and C++ Polyscope 2.x. - [x] Degenerate bounding box tolerance (upstream commit 3198ab5) - [x] `remove_everything()` / `remove_all_groups()` scene reset (upstream commit f34f403) - [x] Improved camera flight interpolation via inverse view matrix (upstream commit 067f760) +- [x] Volume Mesh prism + pyramid cell support (upstream PR #353, commit dcbaedb) --- From 57a3321033ee5db557b9b599258c82b510358d57 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 10:03:52 +0200 Subject: [PATCH 09/14] style(volume_mesh): apply rustfmt to new prism/pyramid code --- .../src/volume_mesh/cell_data.rs | 110 +++++++++++++----- .../src/volume_mesh/mod.rs | 34 +++--- .../src/volume_mesh/slice_geometry.rs | 17 +-- examples/volume_mesh_mixed_cells_demo.rs | 4 +- 4 files changed, 107 insertions(+), 58 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index d44e38b..907a097 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -28,43 +28,103 @@ pub struct FaceData { // ===== Tet ===== const TET_FACES: &[FaceData] = &[ - FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, - FaceData { polygon: &[0, 1, 3], triangulation: &[[0, 1, 3]] }, - FaceData { polygon: &[0, 3, 2], triangulation: &[[0, 3, 2]] }, - FaceData { polygon: &[1, 2, 3], triangulation: &[[1, 2, 3]] }, + FaceData { + polygon: &[0, 2, 1], + triangulation: &[[0, 2, 1]], + }, + FaceData { + polygon: &[0, 1, 3], + triangulation: &[[0, 1, 3]], + }, + FaceData { + polygon: &[0, 3, 2], + triangulation: &[[0, 3, 2]], + }, + FaceData { + polygon: &[1, 2, 3], + triangulation: &[[1, 2, 3]], + }, ]; // ===== Hex ===== // Numbered like in the VTK file-formats diagram, with slots 6 and 7 swapped // to match upstream (see polyscope/src/volume_mesh.cpp:43). const HEX_FACES: &[FaceData] = &[ - FaceData { polygon: &[2, 1, 0, 3], triangulation: &[[2, 1, 0], [2, 0, 3]] }, // Bottom - FaceData { polygon: &[4, 0, 1, 5], triangulation: &[[4, 0, 1], [4, 1, 5]] }, // Front - FaceData { polygon: &[5, 1, 2, 6], triangulation: &[[5, 1, 2], [5, 2, 6]] }, // Right - FaceData { polygon: &[7, 3, 0, 4], triangulation: &[[7, 3, 0], [7, 0, 4]] }, // Left - FaceData { polygon: &[6, 2, 3, 7], triangulation: &[[6, 2, 3], [6, 3, 7]] }, // Back - FaceData { polygon: &[7, 4, 5, 6], triangulation: &[[7, 4, 5], [7, 5, 6]] }, // Top + FaceData { + polygon: &[2, 1, 0, 3], + triangulation: &[[2, 1, 0], [2, 0, 3]], + }, // Bottom + FaceData { + polygon: &[4, 0, 1, 5], + triangulation: &[[4, 0, 1], [4, 1, 5]], + }, // Front + FaceData { + polygon: &[5, 1, 2, 6], + triangulation: &[[5, 1, 2], [5, 2, 6]], + }, // Right + FaceData { + polygon: &[7, 3, 0, 4], + triangulation: &[[7, 3, 0], [7, 0, 4]], + }, // Left + FaceData { + polygon: &[6, 2, 3, 7], + triangulation: &[[6, 2, 3], [6, 3, 7]], + }, // Back + FaceData { + polygon: &[7, 4, 5, 6], + triangulation: &[[7, 4, 5], [7, 5, 6]], + }, // Top ]; // ===== Prism (wedge) ===== // Slots 0,1,2 = bottom triangle; 3,4,5 = top triangle (slots 3,4,5 align with 0,1,2). const PRISM_FACES: &[FaceData] = &[ - FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, // Bottom tri - FaceData { polygon: &[0, 3, 5, 2], triangulation: &[[0, 5, 2], [0, 3, 5]] }, // Side quad 1 - FaceData { polygon: &[2, 5, 4, 1], triangulation: &[[2, 4, 5], [2, 1, 4]] }, // Side quad 2 - FaceData { polygon: &[0, 1, 4, 3], triangulation: &[[3, 0, 4], [0, 1, 4]] }, // Side quad 3 - FaceData { polygon: &[3, 4, 5], triangulation: &[[3, 4, 5]] }, // Top tri + FaceData { + polygon: &[0, 2, 1], + triangulation: &[[0, 2, 1]], + }, // Bottom tri + FaceData { + polygon: &[0, 3, 5, 2], + triangulation: &[[0, 5, 2], [0, 3, 5]], + }, // Side quad 1 + FaceData { + polygon: &[2, 5, 4, 1], + triangulation: &[[2, 4, 5], [2, 1, 4]], + }, // Side quad 2 + FaceData { + polygon: &[0, 1, 4, 3], + triangulation: &[[3, 0, 4], [0, 1, 4]], + }, // Side quad 3 + FaceData { + polygon: &[3, 4, 5], + triangulation: &[[3, 4, 5]], + }, // Top tri ]; // ===== Pyramid ===== // Slots 0..3 = base quad (CCW from outside, looking from -apex toward base); // Slot 4 = apex. const PYRAMID_FACES: &[FaceData] = &[ - FaceData { polygon: &[0, 1, 2, 3], triangulation: &[[0, 3, 2], [0, 2, 1]] }, // Base quad - FaceData { polygon: &[0, 1, 4], triangulation: &[[0, 1, 4]] }, // Side 1 - FaceData { polygon: &[1, 2, 4], triangulation: &[[1, 2, 4]] }, // Side 2 - FaceData { polygon: &[2, 3, 4], triangulation: &[[2, 3, 4]] }, // Side 3 - FaceData { polygon: &[3, 0, 4], triangulation: &[[3, 0, 4]] }, // Side 4 + FaceData { + polygon: &[0, 1, 2, 3], + triangulation: &[[0, 3, 2], [0, 2, 1]], + }, // Base quad + FaceData { + polygon: &[0, 1, 4], + triangulation: &[[0, 1, 4]], + }, // Side 1 + FaceData { + polygon: &[1, 2, 4], + triangulation: &[[1, 2, 4]], + }, // Side 2 + FaceData { + polygon: &[2, 3, 4], + triangulation: &[[2, 3, 4]], + }, // Side 3 + FaceData { + polygon: &[3, 0, 4], + triangulation: &[[3, 0, 4]], + }, // Side 4 ]; /// Returns the face data table for a given cell type. @@ -171,15 +231,9 @@ pub fn decompose_pyramid(cell: &[u32; 8]) -> [[u32; 4]; 2] { let p: [u32; 5] = [cell[0], cell[1], cell[2], cell[3], cell[4]]; if p[0].min(p[2]) < p[1].min(p[3]) { - [ - [p[0], p[2], p[4], p[1]], - [p[0], p[4], p[2], p[3]], - ] + [[p[0], p[2], p[4], p[1]], [p[0], p[4], p[2], p[3]]] } else { - [ - [p[1], p[3], p[4], p[2]], - [p[1], p[4], p[3], p[0]], - ] + [[p[1], p[3], p[4], p[2]], [p[1], p[4], p[3], p[0]]] } } diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 2b1094a..79a5b6a 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -249,7 +249,9 @@ impl VolumeMesh { 2 => VolumeCellType::Prism, 3 => VolumeCellType::Pyramid, 4 => VolumeCellType::Tet, - n => panic!("VolumeMesh cell {cell_idx}: invalid sentinel count {n} (expected 0/2/3/4)"), + n => { + panic!("VolumeMesh cell {cell_idx}: invalid sentinel count {n} (expected 0/2/3/4)") + } } } @@ -1392,13 +1394,13 @@ mod tests { Vec3::new(1.0, 0.0, 1.0), Vec3::new(0.5, 1.0, 1.0), ]; - let mesh = VolumeMesh::new( - "p", - verts, - vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], - ); + let mesh = VolumeMesh::new("p", verts, vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]]); let (_, faces) = mesh.generate_render_geometry(); - assert_eq!(faces.len(), 8, "single prism should have 8 triangles (2 tri + 3*2 quad)"); + assert_eq!( + faces.len(), + 8, + "single prism should have 8 triangles (2 tri + 3*2 quad)" + ); } #[test] @@ -1416,7 +1418,11 @@ mod tests { vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], ); let (_, faces) = mesh.generate_render_geometry(); - assert_eq!(faces.len(), 6, "single pyramid should have 6 triangles (2 base + 4 sides)"); + assert_eq!( + faces.len(), + 6, + "single pyramid should have 6 triangles (2 base + 4 sides)" + ); } #[test] @@ -1448,7 +1454,11 @@ mod tests { ]; let mesh = VolumeMesh::new("m", verts, cells); let (_, faces) = mesh.generate_render_geometry(); - assert_eq!(faces.len(), 4 + 8 + 6, "mixed mesh should sum per-cell triangle counts"); + assert_eq!( + faces.len(), + 4 + 8 + 6, + "mixed mesh should sum per-cell triangle counts" + ); } #[test] @@ -1538,11 +1548,7 @@ mod tests { #[test] fn test_cell_type_detection_prism() { let verts = (0..6).map(|i| Vec3::splat(i as f32)).collect(); - let mesh = VolumeMesh::new( - "p", - verts, - vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], - ); + let mesh = VolumeMesh::new("p", verts, vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]]); assert_eq!(mesh.cell_type(0), VolumeCellType::Prism); } diff --git a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs index e4453ad..7d84d3d 100644 --- a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs +++ b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs @@ -155,18 +155,10 @@ pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> /// * `plane_origin` - A point on the plane /// * `plane_normal` - The plane normal (points toward kept geometry) #[must_use] -pub fn slice_prism( - vertices: [Vec3; 6], - plane_origin: Vec3, - plane_normal: Vec3, -) -> CellSliceResult { +pub fn slice_prism(vertices: [Vec3; 6], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { // Symmetric 3-tet decomposition. Cross-cell consistency is not required for // isolated slicing — the choice doesn't affect correctness. - let tet_indices = [ - [0usize, 5, 4, 3], - [0, 4, 5, 2], - [0, 4, 2, 1], - ]; + let tet_indices = [[0usize, 5, 4, 3], [0, 4, 5, 2], [0, 4, 2, 1]]; let mut all_vertices = Vec::new(); let mut all_interp = Vec::new(); @@ -211,10 +203,7 @@ pub fn slice_pyramid( plane_origin: Vec3, plane_normal: Vec3, ) -> CellSliceResult { - let tet_indices = [ - [0usize, 2, 4, 1], - [0, 4, 2, 3], - ]; + let tet_indices = [[0usize, 2, 4, 1], [0, 4, 2, 3]]; let mut all_vertices = Vec::new(); let mut all_interp = Vec::new(); diff --git a/examples/volume_mesh_mixed_cells_demo.rs b/examples/volume_mesh_mixed_cells_demo.rs index 09024c3..c54c7e3 100644 --- a/examples/volume_mesh_mixed_cells_demo.rs +++ b/examples/volume_mesh_mixed_cells_demo.rs @@ -6,8 +6,8 @@ use glam::Vec3; use polyscope_rs::{ - Result, init, register_hex_mesh, register_prism_mesh, register_pyramid_mesh, - register_tet_mesh, show, + Result, init, register_hex_mesh, register_prism_mesh, register_pyramid_mesh, register_tet_mesh, + show, }; fn main() -> Result<()> { From c454089040a2a1c832fd7f3836b5ba295d9024ca Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 10:14:02 +0200 Subject: [PATCH 10/14] refactor(volume_mesh): simplify cell-type dispatch and slice helpers - cell_type_of(): 1-3 slot checks instead of iter-filter-count, removes per-cell sentinel scan from 9 hot loops. - decompose_to_tets / num_tets: drop per-cell Vec allocation via for_each_tet callback; num_tets is now O(n) instead of O(decompose+alloc). - cell_centroid: iterate only real vert count from cell type. - slice_hex/prism/pyramid: extract shared slice_via_tet_decomposition. - Fix UI label that miscounted prisms/pyramids as hexes. - Update stale "tetrahedral or hexahedral" doc on VolumeMesh. --- .../src/volume_mesh/cell_data.rs | 71 +++++++-- .../src/volume_mesh/mod.rs | 75 ++++----- .../src/volume_mesh/slice_geometry.rs | 149 +++++------------- 3 files changed, 132 insertions(+), 163 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index 907a097..e680d0f 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -138,6 +138,46 @@ pub fn face_data_for(cell_type: VolumeCellType) -> &'static [FaceData] { } } +/// Classifies a cell by examining its trailing sentinel slots. +/// +/// Sentinels are always placed at the end of the 8-slot array, so the cell +/// type is fully determined by which of slots 4/5/6 are sentinel: this +/// requires at most 3 comparisons (Tet hits the first). +#[must_use] +pub fn cell_type_of(cell: &[u32; 8]) -> VolumeCellType { + if cell[4] == u32::MAX { + VolumeCellType::Tet // 4 sentinels in slots 4..7 + } else if cell[5] == u32::MAX { + VolumeCellType::Pyramid // 3 sentinels in slots 5..7 + } else if cell[6] == u32::MAX { + VolumeCellType::Prism // 2 sentinels in slots 6..7 + } else { + VolumeCellType::Hex // 0 sentinels + } +} + +/// Number of real (non-sentinel) vertices in each cell type. +#[must_use] +pub fn num_real_verts(cell_type: VolumeCellType) -> usize { + match cell_type { + VolumeCellType::Tet => 4, + VolumeCellType::Pyramid => 5, + VolumeCellType::Prism => 6, + VolumeCellType::Hex => 8, + } +} + +/// Number of tetrahedra produced by `decompose_cell_to_tets` for each cell type. +#[must_use] +pub fn num_tets_in_cell(cell_type: VolumeCellType) -> usize { + match cell_type { + VolumeCellType::Tet => 1, + VolumeCellType::Pyramid => 2, + VolumeCellType::Prism => 3, + VolumeCellType::Hex => 5, + } +} + /// Builds a canonical (sorted) face key for hashing. /// /// Face polygons can have 3 or 4 unique vertices. Triangular faces leave slot 3 @@ -157,7 +197,7 @@ pub fn canonical_face_key(cell: &[u32; 8], polygon: &[usize]) -> [u32; 4] { /// /// Central-diagonal pattern from Dompierre et al.; the 5-tet split works for /// any convex hex. Matches the previous polyscope-rs decomposition. -const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ +pub(super) const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ [0, 1, 2, 5], [0, 2, 7, 5], [0, 2, 3, 7], @@ -237,16 +277,25 @@ pub fn decompose_pyramid(cell: &[u32; 8]) -> [[u32; 4]; 2] { } } -/// Returns the tet-decomposition for any cell, dispatching on cell type. -#[must_use] -pub fn decompose_cell_to_tets(cell: &[u32; 8], cell_type: VolumeCellType) -> Vec<[u32; 4]> { +/// Invokes `f` once per tet in the decomposition of `cell`, dispatching on +/// cell type. Avoids the per-cell `Vec` allocation of returning a collection. +pub fn for_each_tet(cell: &[u32; 8], cell_type: VolumeCellType, mut f: F) { match cell_type { - VolumeCellType::Tet => vec![[cell[0], cell[1], cell[2], cell[3]]], - VolumeCellType::Hex => HEX_TO_TET_PATTERN - .iter() - .map(|t| [cell[t[0]], cell[t[1]], cell[t[2]], cell[t[3]]]) - .collect(), - VolumeCellType::Prism => decompose_prism(cell).to_vec(), - VolumeCellType::Pyramid => decompose_pyramid(cell).to_vec(), + VolumeCellType::Tet => f([cell[0], cell[1], cell[2], cell[3]]), + VolumeCellType::Hex => { + for t in &HEX_TO_TET_PATTERN { + f([cell[t[0]], cell[t[1]], cell[t[2]], cell[t[3]]]); + } + } + VolumeCellType::Prism => { + for tet in decompose_prism(cell) { + f(tet); + } + } + VolumeCellType::Pyramid => { + for tet in decompose_pyramid(cell) { + f(tet); + } + } } } diff --git a/crates/polyscope-structures/src/volume_mesh/mod.rs b/crates/polyscope-structures/src/volume_mesh/mod.rs index 79a5b6a..87ccf7d 100644 --- a/crates/polyscope-structures/src/volume_mesh/mod.rs +++ b/crates/polyscope-structures/src/volume_mesh/mod.rs @@ -75,10 +75,11 @@ pub enum VolumeCellType { Pyramid, } -/// A volume mesh structure (tetrahedral or hexahedral). +/// A volume mesh structure with mixed cell types (tet / hex / prism / pyramid). /// -/// Cells are stored as arrays of 8 vertex indices. For tetrahedra, -/// only the first 4 indices are used (indices 4-7 are set to `u32::MAX`). +/// Cells are stored as arrays of 8 vertex indices. Unused trailing slots are +/// set to `u32::MAX`: tets use 4 slots, pyramids 5, prisms 6, hexes all 8. +/// The cell type is recovered from the sentinel pattern via [`Self::cell_type`]. pub struct VolumeMesh { name: String, @@ -228,31 +229,12 @@ impl VolumeMesh { /// Returns the cell type of the given cell. /// - /// Cell type is determined by the number of sentinel (`u32::MAX`) indices - /// in the 8-slot cell array (matches upstream C++ Polyscope): - /// - 0 sentinels → `Hex` (8 verts) - /// - 2 sentinels → `Prism` (6 verts) - /// - 3 sentinels → `Pyramid` (5 verts) - /// - 4 sentinels → `Tet` (4 verts) - /// - /// # Panics - /// Panics if `cell_idx` is out of range or the sentinel count is invalid - /// (1, 5, 6, 7, or 8 sentinels). + /// Sentinels are placed at the end of the 8-slot array, so the type is + /// determined by which of slots 4/5/6 holds `u32::MAX` (matches upstream + /// C++ Polyscope's sentinel-count classification). #[must_use] pub fn cell_type(&self, cell_idx: usize) -> VolumeCellType { - let sentinels = self.cells[cell_idx] - .iter() - .filter(|&&v| v == u32::MAX) - .count(); - match sentinels { - 0 => VolumeCellType::Hex, - 2 => VolumeCellType::Prism, - 3 => VolumeCellType::Pyramid, - 4 => VolumeCellType::Tet, - n => { - panic!("VolumeMesh cell {cell_idx}: invalid sentinel count {n} (expected 0/2/3/4)") - } - } + cell_data::cell_type_of(&self.cells[cell_idx]) } /// Returns the vertices. @@ -324,18 +306,20 @@ impl VolumeMesh { /// - Pyramid → 2 tets (consistent diagonal split) #[must_use] pub fn decompose_to_tets(&self) -> Vec<[u32; 4]> { - let mut tets = Vec::new(); - for (cell_idx, cell) in self.cells.iter().enumerate() { - let ct = self.cell_type(cell_idx); - tets.extend(cell_data::decompose_cell_to_tets(cell, ct)); + let mut tets = Vec::with_capacity(self.cells.len() * 2); + for cell in &self.cells { + cell_data::for_each_tet(cell, cell_data::cell_type_of(cell), |t| tets.push(t)); } tets } - /// Returns the number of tetrahedra (including decomposed hexes). + /// Returns the number of tetrahedra (including decomposed cells). #[must_use] pub fn num_tets(&self) -> usize { - self.decompose_to_tets().len() + self.cells + .iter() + .map(|c| cell_data::num_tets_in_cell(cell_data::cell_type_of(c))) + .sum() } /// Computes face counts for interior/exterior detection. @@ -353,15 +337,12 @@ impl VolumeMesh { /// Computes the centroid of a cell as the mean of its real (non-sentinel) vertices. pub(crate) fn cell_centroid(&self, cell: &[u32; 8]) -> Vec3 { + let n = cell_data::num_real_verts(cell_data::cell_type_of(cell)); let mut sum = Vec3::ZERO; - let mut count = 0u32; - for &v in cell { - if v != u32::MAX { - sum += self.vertices[v as usize]; - count += 1; - } + for &v in &cell[..n] { + sum += self.vertices[v as usize]; } - sum / count as f32 + sum / n as f32 } /// Tests if a cell should be visible based on slice planes. @@ -949,15 +930,19 @@ impl VolumeMesh { /// Builds the egui UI for this volume mesh. pub fn build_egui_ui(&mut self, ui: &mut egui::Ui) { - // Info - let num_tets = self.cells.iter().filter(|c| c[4] == u32::MAX).count(); - let num_hexes = self.num_cells() - num_tets; + // Info — count each cell type + let mut counts = [0usize; 4]; // Tet, Hex, Prism, Pyramid + for cell in &self.cells { + counts[cell_data::cell_type_of(cell) as usize] += 1; + } ui.label(format!( - "{} verts, {} cells ({} tets, {} hexes)", + "{} verts, {} cells ({} tet, {} hex, {} prism, {} pyramid)", self.num_vertices(), self.num_cells(), - num_tets, - num_hexes + counts[VolumeCellType::Tet as usize], + counts[VolumeCellType::Hex as usize], + counts[VolumeCellType::Prism as usize], + counts[VolumeCellType::Pyramid as usize], )); // Color diff --git a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs index 7d84d3d..94d9cdd 100644 --- a/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs +++ b/crates/polyscope-structures/src/volume_mesh/slice_geometry.rs @@ -88,35 +88,22 @@ pub fn slice_tet( } } -/// Slice a hexahedron by decomposing into 5 tetrahedra. -/// -/// Hexahedra are sliced by treating them as 5 tetrahedra (using the standard -/// symmetric decomposition), then merging the resulting polygons. +/// Slice a polyhedral cell by decomposing it into tetrahedra, slicing each, +/// and merging the resulting polygons. /// -/// # Arguments -/// * `vertices` - The 8 vertices of the hexahedron in standard ordering -/// * `plane_origin` - A point on the plane -/// * `plane_normal` - The plane normal (points toward kept geometry) -/// -/// # Returns -/// A `CellSliceResult` containing 0, 3-6 vertices depending on the intersection. -#[must_use] -pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { - // Standard decomposition of a hex into 5 tets - // This decomposition is symmetric and works for any hex orientation - let tet_indices = [ - [0, 1, 3, 4], - [1, 2, 3, 6], - [1, 4, 5, 6], - [3, 4, 6, 7], - [1, 3, 4, 6], // Central tet connecting all others - ]; - +/// Used by `slice_hex`, `slice_prism`, and `slice_pyramid`. The tet table's +/// local indices (0..3) are remapped back to cell-local indices via the table. +fn slice_via_tet_decomposition( + vertices: &[Vec3], + tet_indices: &[[usize; 4]], + plane_origin: Vec3, + plane_normal: Vec3, +) -> CellSliceResult { let mut all_vertices = Vec::new(); let mut all_interp = Vec::new(); - for tet in &tet_indices { - let result = slice_tet( + for tet in tet_indices { + let r = slice_tet( vertices[tet[0]], vertices[tet[1]], vertices[tet[2]], @@ -124,20 +111,17 @@ pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> plane_origin, plane_normal, ); - - // Remap interpolation indices from local tet indices to hex indices - for (local_a, local_b, t) in result.interpolation { - let hex_a = tet[local_a as usize] as u32; - let hex_b = tet[local_b as usize] as u32; - all_interp.push((hex_a, hex_b, t)); + for (local_a, local_b, t) in r.interpolation { + all_interp.push(( + tet[local_a as usize] as u32, + tet[local_b as usize] as u32, + t, + )); } - all_vertices.extend(result.vertices); + all_vertices.extend(r.vertices); } - // Merge and deduplicate vertices that are close together merge_slice_vertices(&mut all_vertices, &mut all_interp); - - // Order vertices to form valid polygon if all_vertices.len() >= 3 { order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); } @@ -148,92 +132,43 @@ pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> } } +/// Slice a hexahedron by decomposing into 5 tetrahedra. +/// +/// # Returns +/// A `CellSliceResult` containing 0, 3-6 vertices depending on the intersection. +#[must_use] +pub fn slice_hex(vertices: [Vec3; 8], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { + // Symmetric 5-tet decomposition that works for any hex orientation. + const TETS: [[usize; 4]; 5] = [ + [0, 1, 3, 4], + [1, 2, 3, 6], + [1, 4, 5, 6], + [3, 4, 6, 7], + [1, 3, 4, 6], // Central tet connecting all others + ]; + slice_via_tet_decomposition(&vertices, &TETS, plane_origin, plane_normal) +} + /// Slice a triangular prism by decomposing into 3 tetrahedra. /// -/// # Arguments -/// * `vertices` - The 6 vertices of the prism (slots 0..2 = bottom tri, 3..5 = top tri) -/// * `plane_origin` - A point on the plane -/// * `plane_normal` - The plane normal (points toward kept geometry) +/// Slots 0..2 form the bottom triangle, 3..5 the top. #[must_use] pub fn slice_prism(vertices: [Vec3; 6], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { - // Symmetric 3-tet decomposition. Cross-cell consistency is not required for - // isolated slicing — the choice doesn't affect correctness. - let tet_indices = [[0usize, 5, 4, 3], [0, 4, 5, 2], [0, 4, 2, 1]]; - - let mut all_vertices = Vec::new(); - let mut all_interp = Vec::new(); - - for tet in &tet_indices { - let r = slice_tet( - vertices[tet[0]], - vertices[tet[1]], - vertices[tet[2]], - vertices[tet[3]], - plane_origin, - plane_normal, - ); - for (local_a, local_b, t) in r.interpolation { - let pa = tet[local_a as usize] as u32; - let pb = tet[local_b as usize] as u32; - all_interp.push((pa, pb, t)); - } - all_vertices.extend(r.vertices); - } - - merge_slice_vertices(&mut all_vertices, &mut all_interp); - if all_vertices.len() >= 3 { - order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); - } - - CellSliceResult { - vertices: all_vertices, - interpolation: all_interp, - } + const TETS: [[usize; 4]; 3] = [[0, 5, 4, 3], [0, 4, 5, 2], [0, 4, 2, 1]]; + slice_via_tet_decomposition(&vertices, &TETS, plane_origin, plane_normal) } /// Slice a square pyramid by decomposing into 2 tetrahedra. /// -/// # Arguments -/// * `vertices` - The 5 vertices of the pyramid (slots 0..3 = base quad, 4 = apex) -/// * `plane_origin` - A point on the plane -/// * `plane_normal` - The plane normal (points toward kept geometry) +/// Slots 0..3 form the base quad, slot 4 is the apex. #[must_use] pub fn slice_pyramid( vertices: [Vec3; 5], plane_origin: Vec3, plane_normal: Vec3, ) -> CellSliceResult { - let tet_indices = [[0usize, 2, 4, 1], [0, 4, 2, 3]]; - - let mut all_vertices = Vec::new(); - let mut all_interp = Vec::new(); - - for tet in &tet_indices { - let r = slice_tet( - vertices[tet[0]], - vertices[tet[1]], - vertices[tet[2]], - vertices[tet[3]], - plane_origin, - plane_normal, - ); - for (local_a, local_b, t) in r.interpolation { - let pa = tet[local_a as usize] as u32; - let pb = tet[local_b as usize] as u32; - all_interp.push((pa, pb, t)); - } - all_vertices.extend(r.vertices); - } - - merge_slice_vertices(&mut all_vertices, &mut all_interp); - if all_vertices.len() >= 3 { - order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); - } - - CellSliceResult { - vertices: all_vertices, - interpolation: all_interp, - } + const TETS: [[usize; 4]; 2] = [[0, 2, 4, 1], [0, 4, 2, 3]]; + slice_via_tet_decomposition(&vertices, &TETS, plane_origin, plane_normal) } /// Orders polygon vertices in counter-clockwise order around the centroid. From 2f210742b92a945bc18711f1242d66498f57f84c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 19 May 2026 08:30:46 +0000 Subject: [PATCH 11/14] Merge origin/main and resolve changelog/docs conflicts Agent-Logs-Url: https://github.com/xarthurx/polyscope-rs/sessions/ac8c73f6-5b86-422b-a2ef-a39b2b4091ec Co-authored-by: xarthurx <1921878+xarthurx@users.noreply.github.com> --- .github/workflows/ci.yml | 12 +- .github/workflows/publish.yml | 6 +- CHANGELOG.md | 10 + CLAUDE.md | 2 +- Cargo.toml | 10 +- crates/polyscope-render/src/camera.rs | 82 ++++++++ .../src/curve_network_render.rs | 125 ++++++++++++ crates/polyscope-render/src/engine/mod.rs | 50 +++-- .../src/shaders/curve_network_tube.wgsl | 48 ++++- .../src/shaders/pick_curve_tube.wgsl | 43 +++- .../shaders/reflected_curve_network_tube.wgsl | 41 +++- .../src/volume_grid/scalar_quantity.rs | 190 ++++++++++++++++-- crates/polyscope/src/slice_plane.rs | 38 ++++ crates/polyscope/tests/api_coverage_test.rs | 38 ++++ docs/feature-status.md | 3 + 15 files changed, 642 insertions(+), 56 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1ae8145..ce137d0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -9,13 +9,15 @@ on: env: CARGO_TERM_COLOR: always + # Run JavaScript actions on Node.js 24 (Node 20 deprecates Sep 2026). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true jobs: check: name: Check runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: dtolnay/rust-toolchain@stable - uses: Swatinem/rust-cache@v2 - name: cargo check @@ -25,7 +27,7 @@ jobs: name: Clippy runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: dtolnay/rust-toolchain@stable with: components: clippy @@ -37,7 +39,7 @@ jobs: name: Format runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: dtolnay/rust-toolchain@stable with: components: rustfmt @@ -48,7 +50,7 @@ jobs: name: Test runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 with: lfs: true - uses: dtolnay/rust-toolchain@stable @@ -66,7 +68,7 @@ jobs: env: RUSTDOCFLAGS: -D warnings steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - uses: dtolnay/rust-toolchain@stable - uses: Swatinem/rust-cache@v2 - name: cargo doc diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml index cb1c27a..c79b2dd 100644 --- a/.github/workflows/publish.yml +++ b/.github/workflows/publish.yml @@ -7,19 +7,21 @@ on: env: CARGO_TERM_COLOR: always + # Run JavaScript actions on Node.js 24 (Node 20 deprecates Sep 2026). + FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true jobs: publish: name: Publish runs-on: ubuntu-latest steps: - - uses: actions/checkout@v4 + - uses: actions/checkout@v5 - name: Install Rust toolchain uses: dtolnay/rust-toolchain@stable - name: Cache cargo registry - uses: actions/cache@v4 + uses: actions/cache@v5 with: path: | ~/.cargo/registry diff --git a/CHANGELOG.md b/CHANGELOG.md index 9114ceb..3377521 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 produced by external pipelines must place sentinels in the correct trailing slots; tet meshes built with `new_tet_mesh` continue to work unchanged. +## [0.5.10] - 2026-05-04 + +### Fixed +- VolumeGrid isosurface indexing bug — VolumeGrid stores values as `i + j*nx + k*nx*ny` (z-slowest) but `marching_cubes` indexes as `(i*ny+j)*nz + k` (x-slowest), so isosurfaces were silently X/Z-transposed on uniform grids and visibly broken on non-uniform grids. Fixed by passing dims as `(nz, ny, nx)`, swizzling output positions/normals, and reversing triangle winding (the X↔Z swap is a reflection that flips handedness — without re-winding, `front_facing` shader checks and registered-mesh face normals were inverted). Added regression tests for non-uniform dimensions, anisotropic spacing, and winding consistency (upstream commit e91a709) +- Curve network tube rendering in orthographic projection — `curve_network_tube.wgsl`, `reflected_curve_network_tube.wgsl`, and `pick_curve_tube.wgsl` constructed perspective rays unconditionally (`normalize(world_pos - camera_pos)`), distorting tubes in ortho mode. The shaders now branch on a new `is_orthographic` flag in `CameraUniforms` and emit parallel rays along the world-space view forward direction in ortho mode. The ray-cylinder intersection routine also now handles parallel rays via end-cap intersection — without this, ortho viewing straight down a tube would produce NaN and the tube would disappear or fail to pick (upstream commit 51953c2) + +### Added +- `add_slice_plane_auto()` — convenience constructor that picks the smallest unused "Scene Slice Plane N" index, mirroring C++ Polyscope's no-args `addSlicePlane()` (upstream commit 24ec7e3) +- `SlicePlaneHandle::remove(self)` — consume-self method for removing a slice plane, mirroring C++ `SlicePlane::remove()` + ## [0.5.9] - 2026-03-02 ### Changed diff --git a/CLAUDE.md b/CLAUDE.md index fc77097..5c9d767 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -4,7 +4,7 @@ This file provides guidance to Claude Code when working with this repository. ## Project Overview -polyscope-rs is a Rust-native 3D visualization library for geometric data, ported from C++ [Polyscope](https://polyscope.run). **Core paradigm**: Structures (geometric objects) + Quantities (data on structures). Version 0.5.9, ~100% feature parity with C++ Polyscope 2.x. +polyscope-rs is a Rust-native 3D visualization library for geometric data, ported from C++ [Polyscope](https://polyscope.run). **Core paradigm**: Structures (geometric objects) + Quantities (data on structures). Version 0.5.10, ~100% feature parity with C++ Polyscope 2.x. ## Build Commands diff --git a/Cargo.toml b/Cargo.toml index b679961..dd66dda 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,7 +9,7 @@ members = [ ] [workspace.package] -version = "0.5.9" +version = "0.5.10" edition = "2024" rust-version = "1.85" license = "MIT" @@ -18,10 +18,10 @@ authors = ["polyscope-rs contributors"] [workspace.dependencies] # Internal crates -polyscope-core = { version = "0.5.9", path = "crates/polyscope-core" } -polyscope-render = { version = "0.5.9", path = "crates/polyscope-render" } -polyscope-ui = { version = "0.5.9", path = "crates/polyscope-ui" } -polyscope-structures = { version = "0.5.9", path = "crates/polyscope-structures" } +polyscope-core = { version = "0.5.10", path = "crates/polyscope-core" } +polyscope-render = { version = "0.5.10", path = "crates/polyscope-render" } +polyscope-ui = { version = "0.5.10", path = "crates/polyscope-ui" } +polyscope-structures = { version = "0.5.10", path = "crates/polyscope-structures" } # Math glam = { version = "0.29", features = ["serde"] } diff --git a/crates/polyscope-render/src/camera.rs b/crates/polyscope-render/src/camera.rs index 2ec739b..f7340ce 100644 --- a/crates/polyscope-render/src/camera.rs +++ b/crates/polyscope-render/src/camera.rs @@ -890,4 +890,86 @@ mod tests { "Orthographic zoom in should decrease scale" ); } + + /// Mirrors the WGSL formula used by the curve tube ortho fix: + /// `forward = -vec3(view[0].z, view[1].z, view[2].z)`. + /// glam's `view.x_axis` is column 0, `view.x_axis.z` is `view[0].z` in WGSL. + fn shader_forward_from_view(view: Mat4) -> Vec3 { + -Vec3::new(view.x_axis.z, view.y_axis.z, view.z_axis.z) + } + + /// The shader's view-forward extraction must match `Camera::forward()` for + /// every camera pose, otherwise ortho ray casting in `curve_network_tube.wgsl`, + /// `reflected_curve_network_tube.wgsl`, and `pick_curve_tube.wgsl` will use + /// the wrong ray direction. + #[test] + fn shader_view_forward_matches_camera_forward() { + let cases: &[(Vec3, Vec3, Vec3)] = &[ + // (eye, target, up) + (Vec3::new(0.0, 0.0, 5.0), Vec3::ZERO, Vec3::Y), + (Vec3::new(5.0, 0.0, 0.0), Vec3::ZERO, Vec3::Y), + (Vec3::new(0.0, 5.0, 0.0), Vec3::ZERO, Vec3::Z), + (Vec3::new(3.0, 4.0, 5.0), Vec3::ZERO, Vec3::Y), + ( + Vec3::new(-2.0, 1.0, -3.0), + Vec3::new(1.0, 1.0, 1.0), + Vec3::Y, + ), + ( + Vec3::new(7.0, -2.0, 0.5), + Vec3::new(-1.0, 0.0, 2.0), + Vec3::Z, + ), + ]; + for &(eye, target, up) in cases { + let view = Mat4::look_at_rh(eye, target, up); + let extracted = shader_forward_from_view(view); + let expected = (target - eye).normalize(); + assert!( + extracted.distance(expected) < 1e-5, + "case eye={eye:?} target={target:?} up={up:?}: \ + extracted={extracted:?} expected={expected:?}" + ); + } + } + + /// `CameraUniforms` must be exactly 272 bytes: four mat4x4 (256) + vec3 + + /// f32 flag. If anyone changes the layout, every shader that aliases + /// `camera_pos` as `vec4` and reads `.w` as the ortho flag breaks + /// silently. + #[test] + fn camera_uniforms_layout_is_stable() { + use crate::engine::CameraUniforms; + assert_eq!(std::mem::size_of::(), 272); + assert_eq!(std::mem::align_of::(), 4); + } + + /// Setting `Camera::projection_mode` to `Orthographic` must propagate to + /// the GPU uniform's `is_orthographic` field as 1.0; perspective gives 0.0. + /// The three tube shaders branch on `camera.camera_pos.w > 0.5`, so the + /// exact float values matter. + /// + /// Calls the same `CameraUniforms::from_camera` that `update_camera_uniforms` + /// uses on the render path, so this test catches drift if either side changes. + #[test] + fn ortho_flag_propagates_to_uniform() { + use crate::engine::CameraUniforms; + + let mut camera = Camera::new(1.0); + camera.projection_mode = ProjectionMode::Perspective; + let u_persp = CameraUniforms::from_camera(&camera); + assert!( + u_persp.is_orthographic < 0.5, + "perspective should produce flag below the 0.5 shader threshold, got {}", + u_persp.is_orthographic + ); + + camera.projection_mode = ProjectionMode::Orthographic; + let u_ortho = CameraUniforms::from_camera(&camera); + assert!( + u_ortho.is_orthographic > 0.5, + "orthographic must clear the `> 0.5` shader threshold, got {}", + u_ortho.is_orthographic + ); + } } diff --git a/crates/polyscope-render/src/curve_network_render.rs b/crates/polyscope-render/src/curve_network_render.rs index f175cb5..a4a3339 100644 --- a/crates/polyscope-render/src/curve_network_render.rs +++ b/crates/polyscope-render/src/curve_network_render.rs @@ -421,4 +421,129 @@ mod tests { // Must be 16-byte aligned for GPU uniform buffers assert_eq!(size % 16, 0, "CurveNetworkUniforms must be 16-byte aligned"); } + + // ======================================================================== + // Shader-math validation: parallel-ray cylinder intersection. + // + // The three tube shaders (`curve_network_tube.wgsl`, + // `reflected_curve_network_tube.wgsl`, `pick_curve_tube.wgsl`) include a + // parallel-ray branch needed for ortho viewing straight down a tube. This + // module mirrors that branch in Rust so its math can be unit-tested + // (no GPU required). If the WGSL diverges from this Rust port, the + // ortho head-on tube case will silently regress. + // ======================================================================== + + use glam::Vec3; + + /// Mirrors the parallel-ray branch in the three tube shaders. + /// Returns Some((t, hit_point)) on hit, None on miss. + /// Ignores normal direction since pick variant doesn't compute one. + fn ray_cylinder_parallel_intersect( + ray_origin: Vec3, + ray_dir: Vec3, + cyl_start: Vec3, + cyl_end: Vec3, + cyl_radius: f32, + ) -> Option<(f32, Vec3)> { + let cyl_axis = cyl_end - cyl_start; + let cyl_dir = cyl_axis.normalize(); + let delta = ray_origin - cyl_start; + let delta_perp = delta - cyl_dir.dot(delta) * cyl_dir; + + if delta_perp.length_squared() > cyl_radius * cyl_radius { + return None; + } + let ray_dot_cyl = ray_dir.dot(cyl_dir); + if ray_dot_cyl.abs() < 1e-8 { + return None; + } + let t_start = (cyl_start - ray_origin).dot(cyl_dir) / ray_dot_cyl; + let t_end = (cyl_end - ray_origin).dot(cyl_dir) / ray_dot_cyl; + let mut t_cap = t_start.min(t_end); + if t_cap < 0.001 { + t_cap = t_start.max(t_end); + if t_cap < 0.001 { + return None; + } + } + Some((t_cap, ray_origin + t_cap * ray_dir)) + } + + /// Camera looks straight down a Z-aligned tube in ortho mode. Ray origin + /// is pushed back so t > 0. Should hit the front end cap (closer to camera). + #[test] + fn parallel_ray_through_axis_hits_front_cap() { + let cyl_start = Vec3::new(0.0, 0.0, 0.0); + let cyl_end = Vec3::new(0.0, 0.0, 5.0); + let radius = 0.1_f32; + // Camera at +Z looking down -Z (toward origin) + let ray_dir = Vec3::new(0.0, 0.0, -1.0); + let world_position = Vec3::new(0.0, 0.0, 5.5); // on bbox front (toward camera) + let extent = (cyl_end - cyl_start).length() + 2.0 * radius; + let ray_origin = world_position - extent * ray_dir; + + let hit = ray_cylinder_parallel_intersect(ray_origin, ray_dir, cyl_start, cyl_end, radius); + let (t, p) = hit.expect("parallel ray through axis should hit cylinder cap"); + assert!(t > 0.001, "t must be positive, got {t}"); + // Front cap is cyl_end (z=5.0); hit point z must equal cyl_end.z + assert!( + (p.z - cyl_end.z).abs() < 1e-4, + "expected hit at z={}, got {p:?}", + cyl_end.z + ); + } + + /// Same setup, ray offset within radius — must still hit the cap (the + /// disk is filled, not just the rim). + #[test] + fn parallel_ray_offset_within_radius_hits() { + let cyl_start = Vec3::ZERO; + let cyl_end = Vec3::new(0.0, 0.0, 5.0); + let radius = 0.1_f32; + let ray_dir = Vec3::new(0.0, 0.0, -1.0); + // Offset 0.05 from axis (within radius=0.1) + let world_position = Vec3::new(0.05, 0.0, 5.5); + let extent = (cyl_end - cyl_start).length() + 2.0 * radius; + let ray_origin = world_position - extent * ray_dir; + + let hit = ray_cylinder_parallel_intersect(ray_origin, ray_dir, cyl_start, cyl_end, radius); + assert!(hit.is_some(), "ray within radius should hit cap"); + } + + /// Ray offset beyond radius — must miss. + #[test] + fn parallel_ray_offset_beyond_radius_misses() { + let cyl_start = Vec3::ZERO; + let cyl_end = Vec3::new(0.0, 0.0, 5.0); + let radius = 0.1_f32; + let ray_dir = Vec3::new(0.0, 0.0, -1.0); + let world_position = Vec3::new(0.5, 0.0, 5.5); // 5x radius from axis + let ray_origin = world_position - 10.0 * ray_dir; + + let hit = ray_cylinder_parallel_intersect(ray_origin, ray_dir, cyl_start, cyl_end, radius); + assert!(hit.is_none(), "ray outside radius must miss"); + } + + /// Reverse-direction ray (looking from -Z toward +Z) must hit the OTHER + /// cap (cyl_start side, since that's now nearer to the camera). + #[test] + fn parallel_ray_reverse_direction_hits_other_cap() { + let cyl_start = Vec3::new(0.0, 0.0, 0.0); + let cyl_end = Vec3::new(0.0, 0.0, 5.0); + let radius = 0.1_f32; + let ray_dir = Vec3::new(0.0, 0.0, 1.0); // looking down +Z now + let world_position = Vec3::new(0.0, 0.0, -0.5); + let extent = (cyl_end - cyl_start).length() + 2.0 * radius; + let ray_origin = world_position - extent * ray_dir; + + let (t, p) = + ray_cylinder_parallel_intersect(ray_origin, ray_dir, cyl_start, cyl_end, radius) + .expect("reverse-direction parallel ray should hit"); + assert!(t > 0.001); + assert!( + (p.z - cyl_start.z).abs() < 1e-4, + "expected hit at z={}, got {p:?}", + cyl_start.z + ); + } } diff --git a/crates/polyscope-render/src/engine/mod.rs b/crates/polyscope-render/src/engine/mod.rs index 864f1a0..24b28aa 100644 --- a/crates/polyscope-render/src/engine/mod.rs +++ b/crates/polyscope-render/src/engine/mod.rs @@ -24,16 +24,20 @@ use crate::slice_plane_render::SlicePlaneRenderData; use crate::tone_mapping::ToneMapPass; /// Camera uniforms for GPU. +/// +/// `is_orthographic` is encoded as a float (0.0 = perspective, 1.0 = ortho) so that +/// shaders aliasing `camera_pos` as `vec4` can read it as the `.w` component +/// without type punning. Ray-cast primitive shaders (curve tubes) branch on this to +/// emit parallel rays in ortho mode instead of perspective rays. #[repr(C)] #[derive(Debug, Clone, Copy, bytemuck::Pod, bytemuck::Zeroable)] -#[allow(clippy::pub_underscore_fields)] pub struct CameraUniforms { pub view: [[f32; 4]; 4], pub proj: [[f32; 4]; 4], pub view_proj: [[f32; 4]; 4], pub inv_proj: [[f32; 4]; 4], pub camera_pos: [f32; 3], - pub _padding: f32, + pub is_orthographic: f32, } impl Default for CameraUniforms { @@ -44,7 +48,32 @@ impl Default for CameraUniforms { view_proj: glam::Mat4::IDENTITY.to_cols_array_2d(), inv_proj: glam::Mat4::IDENTITY.to_cols_array_2d(), camera_pos: [0.0, 0.0, 5.0], - _padding: 0.0, + is_orthographic: 0.0, + } + } +} + +impl CameraUniforms { + /// Packs a `Camera` into the GPU uniform layout. Single source of truth for + /// matrix derivation and the `is_orthographic` flag — `update_camera_uniforms` + /// uses this, and tests assert against this directly to catch field drift. + #[must_use] + pub fn from_camera(camera: &crate::camera::Camera) -> Self { + let view = camera.view_matrix(); + let proj = camera.projection_matrix(); + let view_proj = proj * view; + let inv_proj = proj.inverse(); + let is_orthographic = match camera.projection_mode { + crate::camera::ProjectionMode::Orthographic => 1.0, + crate::camera::ProjectionMode::Perspective => 0.0, + }; + Self { + view: view.to_cols_array_2d(), + proj: proj.to_cols_array_2d(), + view_proj: view_proj.to_cols_array_2d(), + inv_proj: inv_proj.to_cols_array_2d(), + camera_pos: camera.position.to_array(), + is_orthographic, } } } @@ -980,20 +1009,7 @@ impl RenderEngine { /// Updates camera uniforms. pub fn update_camera_uniforms(&self) { - let view = self.camera.view_matrix(); - let proj = self.camera.projection_matrix(); - let view_proj = proj * view; - let inv_proj = proj.inverse(); - - let uniforms = CameraUniforms { - view: view.to_cols_array_2d(), - proj: proj.to_cols_array_2d(), - view_proj: view_proj.to_cols_array_2d(), - inv_proj: inv_proj.to_cols_array_2d(), - camera_pos: self.camera.position.to_array(), - _padding: 0.0, - }; - + let uniforms = CameraUniforms::from_camera(&self.camera); self.queue .write_buffer(&self.camera_buffer, 0, bytemuck::cast_slice(&[uniforms])); } diff --git a/crates/polyscope-render/src/shaders/curve_network_tube.wgsl b/crates/polyscope-render/src/shaders/curve_network_tube.wgsl index af32230..c225715 100644 --- a/crates/polyscope-render/src/shaders/curve_network_tube.wgsl +++ b/crates/polyscope-render/src/shaders/curve_network_tube.wgsl @@ -104,6 +104,37 @@ fn ray_cylinder_intersect( // Quadratic coefficients for intersection with infinite cylinder let a = dot(ray_dir_perp, ray_dir_perp); + + // Parallel-ray case: ray parallel (or nearly) to cylinder axis. Quadratic + // degenerates (a≈0); intersect against the two end-cap disks instead. + // Triggered by ortho viewing straight down a tube — without this, all + // fragments produce NaN and the tube disappears or shows garbage pixels. + if (a < 1e-8) { + // Ray must be inside the infinite cylinder to hit either cap. + if (dot(delta_perp, delta_perp) > cyl_radius * cyl_radius) { + return false; + } + let ray_dot_cyl = dot(ray_dir, cyl_dir); + if (abs(ray_dot_cyl) < 1e-8) { + return false; + } + let t_start = dot(cyl_start - ray_origin, cyl_dir) / ray_dot_cyl; + let t_end = dot(cyl_end - ray_origin, cyl_dir) / ray_dot_cyl; + var t_cap = min(t_start, t_end); + if (t_cap < 0.001) { + t_cap = max(t_start, t_end); + if (t_cap < 0.001) { + return false; + } + } + // Outward cap normal: cyl_start cap faces -cyl_dir, cyl_end cap faces +cyl_dir. + let cap_normal = select(cyl_dir, -cyl_dir, t_start < t_end); + *t_hit = t_cap; + *hit_point = ray_origin + t_cap * ray_dir; + *hit_normal = cap_normal; + return true; + } + let b = 2.0 * dot(ray_dir_perp, delta_perp); let c = dot(delta_perp, delta_perp) - cyl_radius * cyl_radius; @@ -159,9 +190,20 @@ fn fs_main(in: VertexOutput) -> FragmentOutput { let tip = edge_vertices[in.edge_id * 2u + 1u].xyz; let radius = uniforms.radius; - // Setup ray from camera through this fragment - let ray_origin = camera.camera_pos.xyz; - let ray_dir = normalize(in.world_position - ray_origin); + // Setup ray. Perspective: from camera through fragment. Orthographic: parallel + // along world-space view forward, pushed back behind the cylinder so the + // intersection routine's t > 0 check always sees the front face of the cylinder + // regardless of which face of the impostor bounding box this fragment came from. + var ray_origin: vec3; + var ray_dir: vec3; + if (camera.camera_pos.w > 0.5) { + ray_dir = -vec3(camera.view[0].z, camera.view[1].z, camera.view[2].z); + let cyl_extent = length(tip - tail) + 2.0 * radius; + ray_origin = in.world_position - cyl_extent * ray_dir; + } else { + ray_origin = camera.camera_pos.xyz; + ray_dir = normalize(in.world_position - ray_origin); + } // Ray-cylinder intersection var t_hit: f32; diff --git a/crates/polyscope-render/src/shaders/pick_curve_tube.wgsl b/crates/polyscope-render/src/shaders/pick_curve_tube.wgsl index 76c7192..27cddc8 100644 --- a/crates/polyscope-render/src/shaders/pick_curve_tube.wgsl +++ b/crates/polyscope-render/src/shaders/pick_curve_tube.wgsl @@ -78,10 +78,32 @@ fn ray_cylinder_intersect( // Quadratic coefficients for intersection with infinite cylinder let a = dot(ray_dir_perp, ray_dir_perp); - if (a < 0.0001) { - // Ray parallel to cylinder axis - return false; + + // Parallel-ray case: ray parallel to cylinder axis (e.g. ortho viewing + // straight down a tube). Quadratic degenerates; intersect with end caps. + // Without this, ortho-mode picks miss tubes viewed end-on. + if (a < 1e-8) { + if (dot(delta_perp, delta_perp) > cyl_radius * cyl_radius) { + return false; + } + let ray_dot_cyl = dot(ray_dir, cyl_dir); + if (abs(ray_dot_cyl) < 1e-8) { + return false; + } + let t_start = dot(cyl_start - ray_origin, cyl_dir) / ray_dot_cyl; + let t_end = dot(cyl_end - ray_origin, cyl_dir) / ray_dot_cyl; + var t_cap = min(t_start, t_end); + if (t_cap < 0.001) { + t_cap = max(t_start, t_end); + if (t_cap < 0.001) { + return false; + } + } + *t_hit = t_cap; + *hit_point = ray_origin + t_cap * ray_dir; + return true; } + let b = 2.0 * dot(ray_dir_perp, delta_perp); let c = dot(delta_perp, delta_perp) - cyl_radius * cyl_radius; @@ -134,9 +156,18 @@ fn fs_main(in: VertexOutput) -> FragmentOutput { // Use at least the minimum pick radius for easier selection let radius = max(pick.radius, pick.min_pick_radius); - // Setup ray from camera through this fragment - let ray_origin = camera.camera_pos.xyz; - let ray_dir = normalize(in.world_position - ray_origin); + // Setup ray. Perspective: from camera through fragment. Orthographic: parallel + // along world-space view forward, pushed back behind the cylinder so t > 0. + var ray_origin: vec3; + var ray_dir: vec3; + if (camera.camera_pos.w > 0.5) { + ray_dir = -vec3(camera.view[0].z, camera.view[1].z, camera.view[2].z); + let cyl_extent = length(tip - tail) + 2.0 * radius; + ray_origin = in.world_position - cyl_extent * ray_dir; + } else { + ray_origin = camera.camera_pos.xyz; + ray_dir = normalize(in.world_position - ray_origin); + } // Ray-cylinder intersection var t_hit: f32; diff --git a/crates/polyscope-render/src/shaders/reflected_curve_network_tube.wgsl b/crates/polyscope-render/src/shaders/reflected_curve_network_tube.wgsl index 5e27e04..13c3b12 100644 --- a/crates/polyscope-render/src/shaders/reflected_curve_network_tube.wgsl +++ b/crates/polyscope-render/src/shaders/reflected_curve_network_tube.wgsl @@ -111,6 +111,32 @@ fn ray_cylinder_intersect( let delta_perp = delta - dot(delta, cyl_dir) * cyl_dir; let a = dot(ray_dir_perp, ray_dir_perp); + + // Parallel-ray case: see curve_network_tube.wgsl for rationale. + if (a < 1e-8) { + if (dot(delta_perp, delta_perp) > cyl_radius * cyl_radius) { + return false; + } + let ray_dot_cyl = dot(ray_dir, cyl_dir); + if (abs(ray_dot_cyl) < 1e-8) { + return false; + } + let t_start = dot(cyl_start - ray_origin, cyl_dir) / ray_dot_cyl; + let t_end = dot(cyl_end - ray_origin, cyl_dir) / ray_dot_cyl; + var t_cap = min(t_start, t_end); + if (t_cap < 0.001) { + t_cap = max(t_start, t_end); + if (t_cap < 0.001) { + return false; + } + } + let cap_normal = select(cyl_dir, -cyl_dir, t_start < t_end); + *t_hit = t_cap; + *hit_point = ray_origin + t_cap * ray_dir; + *hit_normal = cap_normal; + return true; + } + let b = 2.0 * dot(ray_dir_perp, delta_perp); let c = dot(delta_perp, delta_perp) - cyl_radius * cyl_radius; @@ -168,9 +194,18 @@ fn fs_main(in: VertexOutput) -> FragmentOutput { discard; } - // Setup ray from camera through this fragment - let ray_origin = camera.camera_pos.xyz; - let ray_dir = normalize(in.world_position - ray_origin); + // Setup ray. Perspective: from camera through fragment. Orthographic: parallel + // along world-space view forward, pushed back behind the cylinder so t > 0. + var ray_origin: vec3; + var ray_dir: vec3; + if (camera.camera_pos.w > 0.5) { + ray_dir = -vec3(camera.view[0].z, camera.view[1].z, camera.view[2].z); + let cyl_extent = length(tip - tail) + 2.0 * radius; + ray_origin = in.world_position - cyl_extent * ray_dir; + } else { + ray_origin = camera.camera_pos.xyz; + ray_dir = normalize(in.world_position - ray_origin); + } // Ray-cylinder intersection var t_hit: f32; diff --git a/crates/polyscope-structures/src/volume_grid/scalar_quantity.rs b/crates/polyscope-structures/src/volume_grid/scalar_quantity.rs index f0036bf..f03eda2 100644 --- a/crates/polyscope-structures/src/volume_grid/scalar_quantity.rs +++ b/crates/polyscope-structures/src/volume_grid/scalar_quantity.rs @@ -208,19 +208,25 @@ impl VolumeGridNodeScalarQuantity { /// Extracts the isosurface mesh using marching cubes. /// - /// MC output vertices are in grid index space: vertex (i,j,k) has coords - /// that need swizzle(z,y,x) * `grid_spacing` + `bound_min` to transform to world space. + /// `VolumeGrid` stores values as `idx = i + j*nx + k*nx*ny` (z-slowest, x-fastest), + /// but `marching_cubes` indexes as `(i_mc*ny + j_mc)*nz + k_mc` (x-slowest, z-fastest). + /// We pass dims in `(nz, ny, nx)` order so MC's inner loop steps over our X axis, + /// then map output vertices/normals back to world space: + /// `world_x = k_mc, world_y = j_mc, world_z = i_mc` + /// (i.e. swap the X and Z components of every position and normal). + /// + /// That swap has determinant -1, which flips triangle handedness. We restore + /// CCW winding by swapping two indices in every triangle — required for the + /// `front_facing` test in `simple_mesh.wgsl` and for `register_isosurface_as_mesh` + /// which recomputes per-face normals from winding. pub fn extract_isosurface(&mut self) -> &McmMesh { if self.isosurface_mesh_cache.is_none() || self.isosurface_dirty { let nx = self.node_dim.x; let ny = self.node_dim.y; let nz = self.node_dim.z; - let mut mesh = marching_cubes(&self.values, self.isosurface_level, nx, ny, nz); + let mut mesh = marching_cubes(&self.values, self.isosurface_level, nz, ny, nx); - // Transform from MC index space to world space - // MC uses indexing (i * ny + j) * nz + k, output coords are in (i,j,k) space - // Need to map: x_world = x_mc * spacing_z + bound_min.z (swizzle z,y,x) let cell_dim = Vec3::new( (nx - 1).max(1) as f32, (ny - 1).max(1) as f32, @@ -229,26 +235,25 @@ impl VolumeGridNodeScalarQuantity { let spacing = (self.bound_max - self.bound_min) / cell_dim; for v in &mut mesh.vertices { - // MC output: v.x is in i-dimension, v.y in j-dimension, v.z in k-dimension - // Grid layout: i maps to x, j maps to y, k maps to z (no swizzle needed - // since our MC uses same indexing as the grid) *v = Vec3::new( - v.x * spacing.x + self.bound_min.x, + v.z * spacing.x + self.bound_min.x, v.y * spacing.y + self.bound_min.y, - v.z * spacing.z + self.bound_min.z, + v.x * spacing.z + self.bound_min.z, ); } - // Transform normals (only need to scale, then renormalize) for n in &mut mesh.normals { - // Scale normals by inverse spacing to account for non-uniform grid - *n = Vec3::new(n.x / spacing.x, n.y / spacing.y, n.z / spacing.z); + *n = Vec3::new(n.z / spacing.x, n.y / spacing.y, n.x / spacing.z); let len = n.length(); if len > 0.0 { *n /= len; } } + for tri in mesh.indices.chunks_exact_mut(3) { + tri.swap(1, 2); + } + self.isosurface_mesh_cache = Some(mesh); self.isosurface_dirty = false; } @@ -938,3 +943,160 @@ impl Quantity for VolumeGridCellScalarQuantity { self } } + +#[cfg(test)] +mod tests { + use super::*; + + /// Build a non-uniform 3x4x5 grid with field value = `world_x` and verify the + /// isosurface lies on the plane `world_x` = level. Regression test for an + /// indexing-order mismatch between `VolumeGrid` storage (z-slowest) and + /// `marching_cubes` (x-slowest) — equivalent to upstream C++ commit e91a709. + /// Also verifies that triangle winding stays consistent after the + /// handedness-flipping swizzle. + #[test] + fn isosurface_x_aligned_non_uniform_grid() { + let nx: u32 = 3; + let ny: u32 = 4; + let nz: u32 = 5; + let bound_min = Vec3::new(0.0, 0.0, 0.0); + let bound_max = Vec3::new(2.0, 3.0, 4.0); // spacing = (1, 1, 1) + + let mut values = Vec::with_capacity((nx * ny * nz) as usize); + for _k in 0..nz { + for _j in 0..ny { + for i in 0..nx { + values.push(i as f32); + } + } + } + + let mut q = VolumeGridNodeScalarQuantity::new( + "test", + "grid", + values, + UVec3::new(nx, ny, nz), + bound_min, + bound_max, + ); + q.set_isosurface_level(1.5); + + let mesh = q.extract_isosurface(); + assert!(!mesh.vertices.is_empty(), "isosurface should not be empty"); + for v in &mesh.vertices { + assert!( + (v.x - 1.5).abs() < 1e-4, + "vertex {v:?} should lie on plane world_x = 1.5" + ); + assert!( + v.y >= bound_min.y - 1e-4 && v.y <= bound_max.y + 1e-4, + "vertex {v:?} world_y outside [0, 3]" + ); + assert!( + v.z >= bound_min.z - 1e-4 && v.z <= bound_max.z + 1e-4, + "vertex {v:?} world_z outside [0, 4]" + ); + } + + // Field gradient is +X everywhere, so outward normals point in +X. + // Stored normals must agree, AND face normals computed from triangle + // winding must agree — otherwise `front_facing` and registered-mesh + // normals will be inverted. + for n in &mesh.normals { + assert!( + n.x > 0.5, + "stored normal {n:?} should point in +X direction" + ); + } + for tri in mesh.indices.chunks_exact(3) { + let v0 = mesh.vertices[tri[0] as usize]; + let v1 = mesh.vertices[tri[1] as usize]; + let v2 = mesh.vertices[tri[2] as usize]; + let geom_normal = (v1 - v0).cross(v2 - v0); + assert!( + geom_normal.x > 0.0, + "triangle winding gives normal {geom_normal:?}, expected +X" + ); + } + } + + /// Anisotropic-spacing variant: 3x4x5 grid with `bound_max = (20, 3, 4)` so + /// `spacing = (10, 1, 1)`. Field is `world_x / 10` (i.e. integer i), level + /// is 1.5 — vertices should land on plane `world_x = 15`. Catches a bug + /// where someone swizzles indices but uses the wrong `spacing` axis. + #[test] + fn isosurface_anisotropic_spacing_x_axis() { + let nx: u32 = 3; + let ny: u32 = 4; + let nz: u32 = 5; + let bound_min = Vec3::new(0.0, 0.0, 0.0); + let bound_max = Vec3::new(20.0, 3.0, 4.0); + + let mut values = Vec::with_capacity((nx * ny * nz) as usize); + for _k in 0..nz { + for _j in 0..ny { + for i in 0..nx { + values.push(i as f32); + } + } + } + + let mut q = VolumeGridNodeScalarQuantity::new( + "test", + "grid", + values, + UVec3::new(nx, ny, nz), + bound_min, + bound_max, + ); + q.set_isosurface_level(1.5); + + let mesh = q.extract_isosurface(); + assert!(!mesh.vertices.is_empty()); + for v in &mesh.vertices { + assert!( + (v.x - 15.0).abs() < 1e-3, + "vertex {v:?} should lie on plane world_x = 15.0" + ); + } + } + + /// Same setup but field value = `world_z`, so the isosurface should lie on a + /// plane `world_z` = level. Catches axis-confusion regressions on Z specifically. + #[test] + fn isosurface_z_aligned_non_uniform_grid() { + let nx: u32 = 3; + let ny: u32 = 4; + let nz: u32 = 5; + let bound_min = Vec3::new(0.0, 0.0, 0.0); + let bound_max = Vec3::new(2.0, 3.0, 4.0); + + let mut values = Vec::with_capacity((nx * ny * nz) as usize); + for k in 0..nz { + for _j in 0..ny { + for _i in 0..nx { + values.push(k as f32); + } + } + } + + let mut q = VolumeGridNodeScalarQuantity::new( + "test", + "grid", + values, + UVec3::new(nx, ny, nz), + bound_min, + bound_max, + ); + q.set_isosurface_level(2.5); + + let mesh = q.extract_isosurface(); + assert!(!mesh.vertices.is_empty(), "isosurface should not be empty"); + for v in &mesh.vertices { + assert!( + (v.z - 2.5).abs() < 1e-4, + "vertex {v:?} should lie on plane world_z = 2.5" + ); + } + } +} diff --git a/crates/polyscope/src/slice_plane.rs b/crates/polyscope/src/slice_plane.rs index 0138bc7..ead716b 100644 --- a/crates/polyscope/src/slice_plane.rs +++ b/crates/polyscope/src/slice_plane.rs @@ -64,6 +64,36 @@ pub fn add_slice_plane_with_pose( SlicePlaneHandle { name } } +/// Adds a slice plane with an auto-generated name like "Scene Slice Plane 0". +/// +/// Mirrors C++ Polyscope's `addSlicePlane()` (no-args overload). Returns a +/// handle to the newly created plane. The chosen index is the smallest +/// non-negative integer N for which "Scene Slice Plane N" is not already in +/// use, so removing a middle plane and re-adding will reclaim that index. +/// +/// Search and creation happen under a single `with_context_mut` lock — without +/// that, two concurrent callers could pick the same name and one would receive +/// a handle to the other's plane (since core's `add_slice_plane` is an upsert). +pub fn add_slice_plane_auto() -> SlicePlaneHandle { + let name = with_context_mut(|ctx| { + let mut i = 0usize; + let candidate = loop { + let c = format!("Scene Slice Plane {i}"); + if !ctx.has_slice_plane(&c) { + break c; + } + i += 1; + }; + let length_scale = ctx.length_scale; + let center = (ctx.bounding_box.0 + ctx.bounding_box.1) * 0.5; + let plane = ctx.add_slice_plane(&candidate); + plane.set_plane_size(length_scale * 0.25); + plane.set_origin(center); + candidate + }); + SlicePlaneHandle { name } +} + /// Gets an existing slice plane by name. #[must_use] pub fn get_slice_plane(name: &str) -> Option { @@ -116,6 +146,14 @@ impl SlicePlaneHandle { &self.name } + /// Removes this slice plane from the scene. + /// + /// Consumes the handle since the underlying plane is gone afterwards. + /// Mirrors C++ Polyscope's `SlicePlane::remove()`. + pub fn remove(self) { + remove_slice_plane(&self.name); + } + /// Sets the pose (origin and normal) of the slice plane. pub fn set_pose(&self, origin: Vec3, normal: Vec3) -> &Self { with_context_mut(|ctx| { diff --git a/crates/polyscope/tests/api_coverage_test.rs b/crates/polyscope/tests/api_coverage_test.rs index 47ce4c8..121c8c3 100644 --- a/crates/polyscope/tests/api_coverage_test.rs +++ b/crates/polyscope/tests/api_coverage_test.rs @@ -148,6 +148,44 @@ fn api_coverage_tests() { assert!(get_all_slice_planes().is_empty()); } + // --- Test: add_slice_plane_auto generates "Scene Slice Plane N" names --- + { + remove_all_slice_planes(); + + let p0 = add_slice_plane_auto(); + let p1 = add_slice_plane_auto(); + let p2 = add_slice_plane_auto(); + + assert_eq!(p0.name(), "Scene Slice Plane 0"); + assert_eq!(p1.name(), "Scene Slice Plane 1"); + assert_eq!(p2.name(), "Scene Slice Plane 2"); + assert_eq!(get_all_slice_planes().len(), 3); + } + + // --- Test: add_slice_plane_auto skips already-used indices --- + { + remove_all_slice_planes(); + + // Manually claim indices 0 and 2; auto should pick 1 first, then 3. + add_slice_plane("Scene Slice Plane 0"); + add_slice_plane("Scene Slice Plane 2"); + let auto1 = add_slice_plane_auto(); + let auto2 = add_slice_plane_auto(); + assert_eq!(auto1.name(), "Scene Slice Plane 1"); + assert_eq!(auto2.name(), "Scene Slice Plane 3"); + } + + // --- Test: SlicePlaneHandle::remove() consumes handle and removes plane --- + { + remove_all_slice_planes(); + + let plane = add_slice_plane("handle_remove_test"); + assert!(get_slice_plane("handle_remove_test").is_some()); + + plane.remove(); + assert!(get_slice_plane("handle_remove_test").is_none()); + } + // ======================================================================== // GROUP TESTS // ======================================================================== diff --git a/docs/feature-status.md b/docs/feature-status.md index 91588f0..8031fd2 100644 --- a/docs/feature-status.md +++ b/docs/feature-status.md @@ -77,6 +77,9 @@ Feature parity tracking between polyscope-rs and C++ Polyscope 2.x. - [x] Degenerate bounding box tolerance (upstream commit 3198ab5) - [x] `remove_everything()` / `remove_all_groups()` scene reset (upstream commit f34f403) - [x] Improved camera flight interpolation via inverse view matrix (upstream commit 067f760) +- [x] VolumeGrid isosurface indexing fix for non-uniform dimensions — corrects silent X/Z transpose; adds regression tests (upstream commit e91a709) +- [x] Curve network tube ortho-projection ray casting — parallel rays in ortho mode via new `is_orthographic` camera flag (upstream commit 51953c2) +- [x] Slice plane API improvements — `add_slice_plane_auto()` (auto-named) and `SlicePlaneHandle::remove(self)` (upstream commit 24ec7e3) - [x] Volume Mesh prism + pyramid cell support (upstream PR #353, commit dcbaedb) --- From e375485fd27022769644b2cfb0ae379ff8021b16 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 11:28:45 +0200 Subject: [PATCH 12/14] fix(volume_mesh): correct prism side quad 2 triangle winding MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The triangulation [[2,4,5], [2,1,4]] for prism side quad 2 was wound opposite to its polygon order [2,5,4,1], producing inward face normals in init_render_data (which uses (p1-p0).cross(p2-p0)). Other prism faces and tet/hex/pyramid stencils are consistent. Switch to [[2,5,4], [4,1,2]] (same diagonal split, CCW). Upstream C++ Polyscope ships the same inward stencil — likely latent there too. --- crates/polyscope-structures/src/volume_mesh/cell_data.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index e680d0f..7fe15ac 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -89,8 +89,8 @@ const PRISM_FACES: &[FaceData] = &[ }, // Side quad 1 FaceData { polygon: &[2, 5, 4, 1], - triangulation: &[[2, 4, 5], [2, 1, 4]], - }, // Side quad 2 + triangulation: &[[2, 5, 4], [4, 1, 2]], + }, // Side quad 2 (winding fixed vs. upstream stencilPrism, which is inward here) FaceData { polygon: &[0, 1, 4, 3], triangulation: &[[3, 0, 4], [0, 1, 4]], From 28172c984a58cf237b411089931717ac9eac5280 Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 11:32:38 +0200 Subject: [PATCH 13/14] docs(plan): add prism/pyramid implementation plan --- .../2026-05-19-volume-mesh-prism-pyramid.md | 1606 +++++++++++++++++ 1 file changed, 1606 insertions(+) create mode 100644 docs/plans/2026-05-19-volume-mesh-prism-pyramid.md diff --git a/docs/plans/2026-05-19-volume-mesh-prism-pyramid.md b/docs/plans/2026-05-19-volume-mesh-prism-pyramid.md new file mode 100644 index 0000000..a259a1b --- /dev/null +++ b/docs/plans/2026-05-19-volume-mesh-prism-pyramid.md @@ -0,0 +1,1606 @@ +# Volume Mesh Prism & Pyramid Cell Support Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Add prism (6-vertex wedge) and pyramid (5-vertex) cell support to `VolumeMesh`, achieving feature parity with upstream C++ Polyscope PR #353. + +**Architecture:** Extract cell-shape data (stencils, polygon faces, decomposition patterns) into a new `cell_data.rs` submodule under `polyscope-structures/src/volume_mesh/`. Switch `VolumeMesh::cell_type()` from the current sentinel-slot-4 check to a sentinel-count classifier (0=Hex, 2=Prism, 3=Pyramid, 4=Tet — matches upstream). Refactor each per-cell iteration site (`compute_face_counts`, `generate_render_geometry`, `generate_render_geometry_with_culling`, `generate_render_geometry_with_quantities`, `generate_cell_index_per_triangle`, `decompose_to_tets`, `cell_centroid`) to dispatch through a unified `face_data_for(cell_type)` table instead of branching on tet-vs-hex. Add `slice_prism` / `slice_pyramid` to `slice_geometry.rs` using tet-decomposition (same strategy as existing `slice_hex`). Expose `new_prism_mesh` / `new_pyramid_mesh` on `VolumeMesh` and matching `register_*` functions on the public crate. + +**Tech Stack:** Rust, glam (math), wgpu/egui (rendering — unchanged), cargo test, cargo clippy. + +**File-size note:** `volume_mesh/mod.rs` is currently 1569 lines. Adding four-way matches everywhere without refactoring would push it past the 2000-line project limit. The plan extracts ~150 lines into a new `cell_data.rs` and replaces if/else branches with helper dispatch — net growth in `mod.rs` should be small. + +--- + +## File Structure + +- **New:** `crates/polyscope-structures/src/volume_mesh/cell_data.rs` — static tables for each cell type (face polygons, triangulation stencils, tet-decomposition patterns), plus `canonical_face_key` and dispatch helpers `face_data_for`, `face_polygon_for`, `decompose_cell_to_tets`. Single source of truth for shape-data. +- **Modify:** `crates/polyscope-structures/src/volume_mesh/mod.rs` — extend `VolumeCellType` enum, switch `cell_type()` to sentinel-count, replace per-cell tet-vs-hex branches with calls into `cell_data` helpers. Delete the now-relocated constants (`TET_FACE_STENCIL`, `HEX_FACE_STENCIL`, `HEX_TO_TET_PATTERN`, `canonical_face_key`). +- **Modify:** `crates/polyscope-structures/src/volume_mesh/slice_geometry.rs` — add `slice_prism` and `slice_pyramid` alongside `slice_tet`/`slice_hex`. +- **Modify:** `crates/polyscope/src/volume_mesh.rs` — add `register_prism_mesh` and `register_pyramid_mesh` next to the existing `register_tet_mesh`/`register_hex_mesh`. +- **Modify:** `crates/polyscope/src/lib.rs` — re-export the two new register functions if needed (check existing export pattern in the same file). +- **New:** `examples/volume_mesh_mixed_cells_demo.rs` — visual demo registering one mesh of each cell type plus one mixed mesh. +- **Modify:** `docs/feature-status.md` — mark prism/pyramid support as complete under "Completed Features". +- **Modify:** `CHANGELOG.md` — add entry under a new `## [Unreleased]` (or next version) section. + +--- + +### Task 1: Create `cell_data` submodule with new enum variants + +**Files:** +- Create: `crates/polyscope-structures/src/volume_mesh/cell_data.rs` +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs` (declare submodule, extend enum) + +- [ ] **Step 1: Write the failing test** + +Append to the `#[cfg(test)] mod tests` block at the bottom of `crates/polyscope-structures/src/volume_mesh/mod.rs`: + +```rust +#[test] +fn test_cell_type_enum_has_prism_and_pyramid() { + // Compile-time check: pattern match must be exhaustive over all 4 variants. + let types = [ + VolumeCellType::Tet, + VolumeCellType::Hex, + VolumeCellType::Prism, + VolumeCellType::Pyramid, + ]; + for t in types { + let label = match t { + VolumeCellType::Tet => "tet", + VolumeCellType::Hex => "hex", + VolumeCellType::Prism => "prism", + VolumeCellType::Pyramid => "pyramid", + }; + assert!(!label.is_empty()); + } +} +``` + +- [ ] **Step 2: Run test to verify it fails** + +Run: `cargo test -p polyscope-structures test_cell_type_enum_has_prism_and_pyramid -- --nocapture` +Expected: FAIL with "no variant or associated item named `Prism`". + +- [ ] **Step 3: Extend the enum** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, replace lines 64-71 (the `VolumeCellType` enum) with: + +```rust +/// Cell type for volume meshes. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +pub enum VolumeCellType { + /// Tetrahedron (4 vertices, 4 triangular faces) + Tet, + /// Hexahedron (8 vertices, 6 quadrilateral faces) + Hex, + /// Triangular prism / wedge (6 vertices, 2 tri + 3 quad faces) + Prism, + /// Square pyramid (5 vertices, 1 quad + 4 tri faces) + Pyramid, +} +``` + +- [ ] **Step 4: Create the new submodule** + +Create `crates/polyscope-structures/src/volume_mesh/cell_data.rs` with the following content: + +```rust +//! Static shape data for each `VolumeCellType`. +//! +//! Defines, for each cell type: +//! - the unique vertex indices on each face (used for canonical face hashing) +//! - the per-face triangulation stencil (used for rendering and edge detection) +//! - the tet-decomposition pattern (used for slicing and isosurface extraction) +//! +//! Index conventions mirror upstream C++ Polyscope (`src/volume_mesh.cpp`, +//! `stencilTet`, `stencilHex`, `stencilPrism`, `stencilPyramid`): +//! - Tet: 4 verts in slots 0..3, sentinels in 4..7. 4 triangular faces. +//! - Hex: 8 verts in 0..7. 6 quadrilateral faces, each split into 2 triangles. +//! - Prism: 6 verts in 0..5, sentinels in 6..7. Bottom tri (0,1,2), top tri (3,4,5). +//! 5 faces: 1 tri (bottom) + 3 quads (sides) + 1 tri (top) = 8 triangles. +//! - Pyramid: 5 verts in 0..4, sentinels in 5..7. Base quad (0,1,2,3), apex (4). +//! 5 faces: 1 quad (base) + 4 tris (sides) = 6 triangles. + +use super::VolumeCellType; + +/// A face is described by (unique vertex slot list, triangulation). +/// +/// `polygon` lists the unique cell-local vertex slots in CCW order around the +/// face (3 for triangles, 4 for quads). `triangulation` is the list of triangle +/// stencils used to render the face. +pub struct FaceData { + pub polygon: &'static [usize], + pub triangulation: &'static [[usize; 3]], +} + +// ===== Tet ===== +const TET_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, + FaceData { polygon: &[0, 1, 3], triangulation: &[[0, 1, 3]] }, + FaceData { polygon: &[0, 3, 2], triangulation: &[[0, 3, 2]] }, + FaceData { polygon: &[1, 2, 3], triangulation: &[[1, 2, 3]] }, +]; + +// ===== Hex ===== +// Numbered like in the VTK file-formats diagram, with slots 6 and 7 swapped +// to match upstream (see polyscope/src/volume_mesh.cpp:43). +const HEX_FACES: &[FaceData] = &[ + FaceData { polygon: &[2, 1, 0, 3], triangulation: &[[2, 1, 0], [2, 0, 3]] }, // Bottom + FaceData { polygon: &[4, 0, 1, 5], triangulation: &[[4, 0, 1], [4, 1, 5]] }, // Front + FaceData { polygon: &[5, 1, 2, 6], triangulation: &[[5, 1, 2], [5, 2, 6]] }, // Right + FaceData { polygon: &[7, 3, 0, 4], triangulation: &[[7, 3, 0], [7, 0, 4]] }, // Left + FaceData { polygon: &[6, 2, 3, 7], triangulation: &[[6, 2, 3], [6, 3, 7]] }, // Back + FaceData { polygon: &[7, 4, 5, 6], triangulation: &[[7, 4, 5], [7, 5, 6]] }, // Top +]; + +// ===== Prism (wedge) ===== +// Slots 0,1,2 = bottom triangle; 3,4,5 = top triangle (slots 3,4,5 align with 0,1,2). +const PRISM_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 2, 1], triangulation: &[[0, 2, 1]] }, // Bottom tri + FaceData { polygon: &[0, 3, 5, 2], triangulation: &[[0, 5, 2], [0, 3, 5]] }, // Side quad 1 + FaceData { polygon: &[2, 5, 4, 1], triangulation: &[[2, 4, 5], [2, 1, 4]] }, // Side quad 2 + FaceData { polygon: &[0, 1, 4, 3], triangulation: &[[3, 0, 4], [0, 1, 4]] }, // Side quad 3 + FaceData { polygon: &[3, 4, 5], triangulation: &[[3, 4, 5]] }, // Top tri +]; + +// ===== Pyramid ===== +// Slots 0..3 = base quad (CCW from outside, looking from -apex toward base); +// Slot 4 = apex. +const PYRAMID_FACES: &[FaceData] = &[ + FaceData { polygon: &[0, 1, 2, 3], triangulation: &[[0, 3, 2], [0, 2, 1]] }, // Base quad + FaceData { polygon: &[0, 1, 4], triangulation: &[[0, 1, 4]] }, // Side 1 + FaceData { polygon: &[1, 2, 4], triangulation: &[[1, 2, 4]] }, // Side 2 + FaceData { polygon: &[2, 3, 4], triangulation: &[[2, 3, 4]] }, // Side 3 + FaceData { polygon: &[3, 0, 4], triangulation: &[[3, 0, 4]] }, // Side 4 +]; + +/// Returns the face data table for a given cell type. +#[must_use] +pub fn face_data_for(cell_type: VolumeCellType) -> &'static [FaceData] { + match cell_type { + VolumeCellType::Tet => TET_FACES, + VolumeCellType::Hex => HEX_FACES, + VolumeCellType::Prism => PRISM_FACES, + VolumeCellType::Pyramid => PYRAMID_FACES, + } +} + +/// Number of vertices used by a cell type (non-sentinel slots in the `[u32; 8]`). +#[must_use] +pub fn num_verts_in_cell(cell_type: VolumeCellType) -> usize { + match cell_type { + VolumeCellType::Tet => 4, + VolumeCellType::Hex => 8, + VolumeCellType::Prism => 6, + VolumeCellType::Pyramid => 5, + } +} + +/// Builds a canonical (sorted) face key for hashing. +/// +/// Face polygons can have 3 or 4 unique vertices. Triangular faces leave slot 3 +/// as `u32::MAX` so that triangle and quad keys never collide. +#[must_use] +pub fn canonical_face_key(cell: &[u32; 8], polygon: &[usize]) -> [u32; 4] { + let mut key = [u32::MAX; 4]; + debug_assert!(polygon.len() == 3 || polygon.len() == 4); + for (i, &slot) in polygon.iter().enumerate() { + key[i] = cell[slot]; + } + key.sort_unstable(); + key +} +``` + +Then in `crates/polyscope-structures/src/volume_mesh/mod.rs`, add the submodule declaration. Look for the existing `mod ...;` declarations near line 44, and add: + +```rust +mod cell_data; +``` + +between `mod color_quantity;` and `mod scalar_quantity;`. Keep alphabetical order if other declarations follow it. + +- [ ] **Step 5: Run test to verify it passes** + +Run: `cargo test -p polyscope-structures test_cell_type_enum_has_prism_and_pyramid -- --nocapture` +Expected: PASS. + +- [ ] **Step 6: Run clippy** + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 7: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/cell_data.rs crates/polyscope-structures/src/volume_mesh/mod.rs +git commit -m "feat(volume_mesh): scaffold prism/pyramid cell types and shape data" +``` + +--- + +### Task 2: Switch `cell_type()` to sentinel-count detection + +**Files:** +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs:189-197` (the `cell_type` method) + +- [ ] **Step 1: Write the failing tests** + +Append to the `#[cfg(test)] mod tests` block: + +```rust +#[test] +fn test_cell_type_detection_tet() { + let mesh = VolumeMesh::new( + "t", + vec![Vec3::ZERO, Vec3::X, Vec3::Y, Vec3::Z], + vec![[0, 1, 2, 3, u32::MAX, u32::MAX, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Tet); +} + +#[test] +fn test_cell_type_detection_hex() { + // 8 distinct vertices, no sentinels + let verts = (0..8).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new("h", verts, vec![[0, 1, 2, 3, 4, 5, 6, 7]]); + assert_eq!(mesh.cell_type(0), VolumeCellType::Hex); +} + +#[test] +fn test_cell_type_detection_prism() { + // 6 verts in slots 0..5, two sentinels in slots 6..7 + let verts = (0..6).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new( + "p", + verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Prism); +} + +#[test] +fn test_cell_type_detection_pyramid() { + // 5 verts in slots 0..4, three sentinels in slots 5..7 + let verts = (0..5).map(|i| Vec3::splat(i as f32)).collect(); + let mesh = VolumeMesh::new( + "py", + verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + assert_eq!(mesh.cell_type(0), VolumeCellType::Pyramid); +} +``` + +- [ ] **Step 2: Run tests to verify they fail** + +Run: `cargo test -p polyscope-structures test_cell_type_detection -- --nocapture` +Expected: FAIL — prism and pyramid currently misclassified as `Hex` (any non-tet returns `Hex` under the old check). + +- [ ] **Step 3: Replace the detection implementation** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, replace the `cell_type` method body (currently lines 189-197): + +```rust + /// Returns the cell type of the given cell. + #[must_use] + pub fn cell_type(&self, cell_idx: usize) -> VolumeCellType { + if self.cells[cell_idx][4] == u32::MAX { + VolumeCellType::Tet + } else { + VolumeCellType::Hex + } + } +``` + +with: + +```rust + /// Returns the cell type of the given cell. + /// + /// Cell type is determined by the number of sentinel (`u32::MAX`) indices + /// in the 8-slot cell array (matches upstream C++ Polyscope): + /// - 0 sentinels → `Hex` (8 verts) + /// - 2 sentinels → `Prism` (6 verts) + /// - 3 sentinels → `Pyramid` (5 verts) + /// - 4 sentinels → `Tet` (4 verts) + /// + /// # Panics + /// Panics if `cell_idx` is out of range or the sentinel count is invalid + /// (1, 5, 6, 7, or 8 sentinels). + #[must_use] + pub fn cell_type(&self, cell_idx: usize) -> VolumeCellType { + let sentinels = self.cells[cell_idx] + .iter() + .filter(|&&v| v == u32::MAX) + .count(); + match sentinels { + 0 => VolumeCellType::Hex, + 2 => VolumeCellType::Prism, + 3 => VolumeCellType::Pyramid, + 4 => VolumeCellType::Tet, + n => panic!("VolumeMesh cell {cell_idx}: invalid sentinel count {n} (expected 0/2/3/4)"), + } + } +``` + +- [ ] **Step 4: Run tests to verify they pass** + +Run: `cargo test -p polyscope-structures test_cell_type_detection -- --nocapture` +Expected: All 4 PASS. + +- [ ] **Step 5: Re-run the full structure test suite to check for regressions** + +Run: `cargo test -p polyscope-structures -- --nocapture` +Expected: All existing tests still PASS. + +- [ ] **Step 6: Run clippy** + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 7: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/mod.rs +git commit -m "feat(volume_mesh): use sentinel-count to classify cell type" +``` + +--- + +### Task 3: Wire prism / pyramid into `decompose_to_tets` and `cell_centroid` + +**Files:** +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs:259-336` (decompose_to_tets, cell_centroid) +- Modify: `crates/polyscope-structures/src/volume_mesh/cell_data.rs` (add tet-decomposition helpers) + +- [ ] **Step 1: Write the failing tests** + +Append to the test module: + +```rust +#[test] +fn test_decompose_prism_to_tets() { + // Unit triangular prism: bottom tri at z=0, top tri at z=1 + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), // 0 + Vec3::new(1.0, 0.0, 0.0), // 1 + Vec3::new(0.5, 1.0, 0.0), // 2 + Vec3::new(0.0, 0.0, 1.0), // 3 + Vec3::new(1.0, 0.0, 1.0), // 4 + Vec3::new(0.5, 1.0, 1.0), // 5 + ]; + let mesh = VolumeMesh::new( + "prism_only", + verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + let tets = mesh.decompose_to_tets(); + // A prism decomposes into exactly 3 tetrahedra + assert_eq!(tets.len(), 3, "prism should decompose to 3 tets"); + for tet in &tets { + for &v in tet { + assert!(v < 6, "tet vertex index {v} out of range for prism"); + } + } +} + +#[test] +fn test_decompose_pyramid_to_tets() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), // 0 + Vec3::new(1.0, 0.0, 0.0), // 1 + Vec3::new(1.0, 1.0, 0.0), // 2 + Vec3::new(0.0, 1.0, 0.0), // 3 + Vec3::new(0.5, 0.5, 1.0), // 4 (apex) + ]; + let mesh = VolumeMesh::new( + "pyr_only", + verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + let tets = mesh.decompose_to_tets(); + // A pyramid decomposes into exactly 2 tetrahedra + assert_eq!(tets.len(), 2, "pyramid should decompose to 2 tets"); + for tet in &tets { + for &v in tet { + assert!(v < 5, "tet vertex index {v} out of range for pyramid"); + } + } +} + +#[test] +fn test_cell_centroid_prism() { + // Build a prism with known centroid + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(2.0, 0.0, 0.0), + Vec3::new(1.0, 2.0, 0.0), + Vec3::new(0.0, 0.0, 4.0), + Vec3::new(2.0, 0.0, 4.0), + Vec3::new(1.0, 2.0, 4.0), + ]; + let mesh = VolumeMesh::new( + "p", verts.clone(), + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + // Centroid is mean of the 6 vertices + let expected: Vec3 = verts.iter().copied().sum::() / 6.0; + let centroid = mesh.cell_centroid(&mesh.cells()[0]); + assert!((centroid - expected).length() < 1e-5); +} +``` + +- [ ] **Step 2: Run tests to verify they fail** + +Run: `cargo test -p polyscope-structures test_decompose_prism_to_tets test_decompose_pyramid_to_tets test_cell_centroid_prism -- --nocapture` +Expected: FAIL — `decompose_to_tets` currently only handles tet/hex (returns 0 tets for prism/pyramid) and `cell_centroid` averages 8 vertices for non-tet cells (incorrect for prism/pyramid since slots 6-7 hold `u32::MAX`). + +- [ ] **Step 3: Add tet-decomposition helpers to `cell_data.rs`** + +Append to `crates/polyscope-structures/src/volume_mesh/cell_data.rs`: + +```rust +/// Decomposes a hex cell into 5 tetrahedra using a fixed diagonal pattern. +/// +/// This matches the existing polyscope-rs decomposition (the central-diagonal +/// pattern from Dompierre et al.). The 5-tet split works for any convex hex. +const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ + [0, 1, 2, 5], + [0, 2, 7, 5], + [0, 2, 3, 7], + [0, 5, 7, 4], + [2, 7, 5, 6], +]; + +/// Decomposes a triangular prism into 3 tetrahedra. +/// +/// Picks a consistent diagonal split on the quad face opposite the +/// lowest-numbered vertex, matching the upstream algorithm (decomposePrism in +/// polyscope/src/volume_mesh.cpp). Consistency across adjacent cells matters +/// for mixed meshes so that shared faces are tessellated identically and no +/// gaps appear in slice caps. +pub fn decompose_prism(cell: &[u32; 8]) -> [[u32; 4]; 3] { + let mut p: [u32; 6] = [cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]]; + + // Find index of smallest vertex + let min_idx = (0..6).min_by_key(|&i| p[i]).unwrap(); + + if min_idx < 3 { + // Smallest vertex is in the bottom triangle. Rotate bottom so it lands in slot 0. + let rot = match min_idx { + 0 => 0, + 1 => 2, + _ => 1, // min_idx == 2 + }; + rotate_prism_in_place(&mut p, rot); + } else { + // Smallest is in the top triangle. Reflect top<->bottom, then rotate. + let top_pos = min_idx - 3; + let rot = match top_pos { + 0 => 0, + 1 => 2, + _ => 1, // top_pos == 2 + }; + p.swap(0, 3); + p.swap(1, 4); + p.swap(2, 5); + rotate_prism_in_place(&mut p, rot); + } + + // Split the quad opposite V0 along the diagonal containing the smaller of + // {p[2], p[4]} vs {p[1], p[5]}, so adjacent cells choose the same split. + if p[2].min(p[4]) < p[1].min(p[5]) { + // Diagonal 2-4 + [ + [p[0], p[5], p[4], p[3]], + [p[0], p[4], p[5], p[2]], + [p[0], p[4], p[2], p[1]], + ] + } else { + // Diagonal 1-5 + [ + [p[0], p[5], p[4], p[3]], + [p[0], p[1], p[5], p[2]], + [p[0], p[5], p[1], p[4]], + ] + } +} + +fn rotate_prism_in_place(p: &mut [u32; 6], rot: usize) { + const BOTTOM_ROT: [[usize; 3]; 3] = [[0, 1, 2], [1, 2, 0], [2, 0, 1]]; + const TOP_ROT: [[usize; 3]; 3] = [[3, 4, 5], [4, 5, 3], [5, 3, 4]]; + let src = *p; + for i in 0..3 { + p[i] = src[BOTTOM_ROT[rot][i]]; + p[i + 3] = src[TOP_ROT[rot][i]]; + } +} + +/// Decomposes a square pyramid into 2 tetrahedra by splitting the base quad +/// along the diagonal containing the smaller of {p[0], p[2]} vs {p[1], p[3]}. +/// Consistent split ensures adjacent cells tessellate the shared face the same way. +pub fn decompose_pyramid(cell: &[u32; 8]) -> [[u32; 4]; 2] { + let p: [u32; 5] = [cell[0], cell[1], cell[2], cell[3], cell[4]]; + + if p[0].min(p[2]) < p[1].min(p[3]) { + // Diagonal 0-2 + [ + [p[0], p[2], p[4], p[1]], + [p[0], p[4], p[2], p[3]], + ] + } else { + // Diagonal 1-3 + [ + [p[1], p[3], p[4], p[2]], + [p[1], p[4], p[3], p[0]], + ] + } +} + +/// Returns the tet-decomposition for any cell, dispatching on cell type. +pub fn decompose_cell_to_tets(cell: &[u32; 8], cell_type: VolumeCellType) -> Vec<[u32; 4]> { + match cell_type { + VolumeCellType::Tet => vec![[cell[0], cell[1], cell[2], cell[3]]], + VolumeCellType::Hex => HEX_TO_TET_PATTERN + .iter() + .map(|t| [cell[t[0]], cell[t[1]], cell[t[2]], cell[t[3]]]) + .collect(), + VolumeCellType::Prism => decompose_prism(cell).to_vec(), + VolumeCellType::Pyramid => decompose_pyramid(cell).to_vec(), + } +} +``` + +- [ ] **Step 4: Replace `decompose_to_tets` in `mod.rs`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, replace `decompose_to_tets` (currently lines 259-284): + +```rust + /// Decomposes all cells into tetrahedra. + /// Tets pass through unchanged, hexes are decomposed into 5 tets. + #[must_use] + pub fn decompose_to_tets(&self) -> Vec<[u32; 4]> { + let mut tets = Vec::new(); + + for cell in &self.cells { + if cell[4] == u32::MAX { + // Already a tet + tets.push([cell[0], cell[1], cell[2], cell[3]]); + } else { + // Hex - decompose using diagonal pattern (5 tets) + for tet_local in &HEX_TO_TET_PATTERN { + let tet = [ + cell[tet_local[0]], + cell[tet_local[1]], + cell[tet_local[2]], + cell[tet_local[3]], + ]; + tets.push(tet); + } + } + } + + tets + } +``` + +with: + +```rust + /// Decomposes all cells into tetrahedra. + /// + /// Decomposition counts per cell type: + /// - Tet → 1 tet (passthrough) + /// - Hex → 5 tets (fixed diagonal pattern) + /// - Prism → 3 tets (consistent diagonal split) + /// - Pyramid → 2 tets (consistent diagonal split) + #[must_use] + pub fn decompose_to_tets(&self) -> Vec<[u32; 4]> { + let mut tets = Vec::new(); + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + tets.extend(cell_data::decompose_cell_to_tets(cell, ct)); + } + tets + } +``` + +Add `use cell_data;` (or `use super::cell_data;` depending on context) near the top of `mod.rs` — check the existing `use` statements around lines 56-62. If those are absolute paths, just refer to `cell_data::...` directly since the submodule is in scope. + +- [ ] **Step 5: Replace `cell_centroid` in `mod.rs`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, replace `cell_centroid` (currently lines 320-336): + +```rust + /// Computes the centroid of a cell. + fn cell_centroid(&self, cell: &[u32; 8]) -> Vec3 { + if cell[4] == u32::MAX { + // Tetrahedron: average of 4 vertices + let sum = self.vertices[cell[0] as usize] + + self.vertices[cell[1] as usize] + + self.vertices[cell[2] as usize] + + self.vertices[cell[3] as usize]; + sum / 4.0 + } else { + // Hexahedron: average of 8 vertices + let sum = (0..8) + .map(|i| self.vertices[cell[i] as usize]) + .fold(Vec3::ZERO, |a, b| a + b); + sum / 8.0 + } + } +``` + +with: + +```rust + /// Computes the centroid of a cell as the mean of its real (non-sentinel) vertices. + pub(crate) fn cell_centroid(&self, cell: &[u32; 8]) -> Vec3 { + let mut sum = Vec3::ZERO; + let mut count = 0u32; + for &v in cell { + if v != u32::MAX { + sum += self.vertices[v as usize]; + count += 1; + } + } + sum / count as f32 + } +``` + +(Note the visibility change to `pub(crate)` to allow the test to call it directly. If the existing visibility was private, leave it `fn cell_centroid` and instead expose a test helper or use a non-direct-access test; otherwise `pub(crate)` is fine for crate-internal tests in the same module.) + +- [ ] **Step 6: Delete the now-orphaned `HEX_TO_TET_PATTERN` constant in `mod.rs`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, delete the constant block at line 1442-1449: + +```rust +/// Diagonal decomposition patterns (5 tets). +const HEX_TO_TET_PATTERN: [[usize; 4]; 5] = [ + [0, 1, 2, 5], + [0, 2, 7, 5], + [0, 2, 3, 7], + [0, 5, 7, 4], + [2, 7, 5, 6], +]; +``` + +It now lives in `cell_data.rs`. + +- [ ] **Step 7: Run tests to verify they pass** + +Run: `cargo test -p polyscope-structures test_decompose_prism_to_tets test_decompose_pyramid_to_tets test_cell_centroid_prism -- --nocapture` +Expected: PASS. + +Run the existing hex/tet tests too: + +Run: `cargo test -p polyscope-structures test_hex_to_tet_decomposition test_single_tet_all_exterior -- --nocapture` +Expected: PASS (regression check). + +- [ ] **Step 8: Run clippy** + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 9: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/mod.rs crates/polyscope-structures/src/volume_mesh/cell_data.rs +git commit -m "feat(volume_mesh): tet-decomposition + centroid for prism/pyramid" +``` + +--- + +### Task 4: Refactor face-counting and render geometry to dispatch through `cell_data` + +This task consolidates the six near-duplicate per-cell loops (`compute_face_counts`, `compute_face_counts_with_culling`, `generate_render_geometry`, `generate_render_geometry_with_culling`, `generate_render_geometry_with_quantities`, `generate_cell_index_per_triangle`, `generate_cell_index_per_triangle_with_culling`) so each handles all 4 cell types uniformly. + +**Files:** +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs:292-829` (the seven functions listed above) + +- [ ] **Step 1: Write the failing tests** + +Append to the test module: + +```rust +#[test] +fn test_single_prism_all_exterior() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let mesh = VolumeMesh::new( + "p", verts, + vec![[0, 1, 2, 3, 4, 5, u32::MAX, u32::MAX]], + ); + let (_, faces) = mesh.generate_render_geometry(); + // Prism has 5 faces total: 2 tris + 3 quads = 8 triangles after triangulation + assert_eq!(faces.len(), 8, "single prism should have 8 triangles (2 tri + 3*2 quad)"); +} + +#[test] +fn test_single_pyramid_all_exterior() { + let verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let mesh = VolumeMesh::new( + "py", verts, + vec![[0, 1, 2, 3, 4, u32::MAX, u32::MAX, u32::MAX]], + ); + let (_, faces) = mesh.generate_render_geometry(); + // Pyramid: 1 quad base + 4 tri sides = 2 + 4 = 6 triangles + assert_eq!(faces.len(), 6, "single pyramid should have 6 triangles (2 base + 4 sides)"); +} + +#[test] +fn test_mixed_cell_mesh_renders_all() { + // One tet + one prism + one pyramid sharing zero faces -> all faces exterior + let verts = vec![ + // Tet + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + // Prism (translated +5 in x) + Vec3::new(5.0, 0.0, 0.0), + Vec3::new(6.0, 0.0, 0.0), + Vec3::new(5.5, 1.0, 0.0), + Vec3::new(5.0, 0.0, 1.0), + Vec3::new(6.0, 0.0, 1.0), + Vec3::new(5.5, 1.0, 1.0), + // Pyramid (translated +10 in x) + Vec3::new(10.0, 0.0, 0.0), + Vec3::new(11.0, 0.0, 0.0), + Vec3::new(11.0, 1.0, 0.0), + Vec3::new(10.0, 1.0, 0.0), + Vec3::new(10.5, 0.5, 1.0), + ]; + let cells = vec![ + [0, 1, 2, 3, u32::MAX, u32::MAX, u32::MAX, u32::MAX], // tet (4 tris) + [4, 5, 6, 7, 8, 9, u32::MAX, u32::MAX], // prism (8 tris) + [10, 11, 12, 13, 14, u32::MAX, u32::MAX, u32::MAX], // pyramid (6 tris) + ]; + let mesh = VolumeMesh::new("m", verts, cells); + let (_, faces) = mesh.generate_render_geometry(); + assert_eq!(faces.len(), 4 + 8 + 6, "mixed mesh should sum per-cell triangle counts"); +} +``` + +- [ ] **Step 2: Run tests to verify they fail** + +Run: `cargo test -p polyscope-structures test_single_prism_all_exterior test_single_pyramid_all_exterior test_mixed_cell_mesh_renders_all -- --nocapture` +Expected: FAIL — prism/pyramid currently go down the `else` (hex) branch and emit 12 triangles each. + +- [ ] **Step 3: Replace `compute_face_counts`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, replace `compute_face_counts` (currently lines 293-318): + +```rust + /// Computes face counts for interior/exterior detection. + fn compute_face_counts(&self) -> HashMap<[u32; 4], usize> { + let mut face_counts: HashMap<[u32; 4], usize> = HashMap::new(); + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + *face_counts.entry(key).or_insert(0) += 1; + } + } + face_counts + } +``` + +- [ ] **Step 4: Replace `compute_face_counts_with_culling`** + +Replace the function at lines 356-389 with: + +```rust + /// Same as `compute_face_counts` but skips cells culled by slice planes. + fn compute_face_counts_with_culling( + &self, + planes: &[(Vec3, Vec3)], + ) -> HashMap<[u32; 4], usize> { + let mut face_counts: HashMap<[u32; 4], usize> = HashMap::new(); + for (cell_idx, cell) in self.cells.iter().enumerate() { + if !self.is_cell_visible(cell, planes) { + continue; + } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + *face_counts.entry(key).or_insert(0) += 1; + } + } + face_counts + } +``` + +- [ ] **Step 5: Replace `generate_render_geometry`** + +Replace the function at lines 391-434 with: + +```rust + /// Generates triangulated exterior faces for rendering. + fn generate_render_geometry(&self) -> (Vec, Vec<[u32; 3]>) { + let face_counts = self.compute_face_counts(); + let mut positions = Vec::new(); + let mut faces = Vec::new(); + + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts[&key] != 1 { + continue; + } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + faces.push([base_idx, base_idx + 1, base_idx + 2]); + } + } + } + + (positions, faces) + } +``` + +- [ ] **Step 6: Replace `generate_render_geometry_with_culling`** + +Replace the function at lines 436-488 with: + +```rust + /// Generates triangulated exterior faces with cell culling based on slice planes. + fn generate_render_geometry_with_culling( + &self, + planes: &[(Vec3, Vec3)], + ) -> (Vec, Vec<[u32; 3]>) { + let face_counts = self.compute_face_counts_with_culling(planes); + let mut positions = Vec::new(); + let mut faces = Vec::new(); + + for (cell_idx, cell) in self.cells.iter().enumerate() { + if !self.is_cell_visible(cell, planes) { + continue; + } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; + } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + faces.push([base_idx, base_idx + 1, base_idx + 2]); + } + } + } + + (positions, faces) + } +``` + +- [ ] **Step 7: Replace `generate_render_geometry_with_quantities` per-cell loop** + +The body of this function spans roughly lines 492-... and has two passes (geometry + normals + per-vertex quantity values). The first pass iterates cells; replace the per-cell `if cell[4] == u32::MAX { tet branch } else { hex branch }` block with: + +```rust + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts[&key] != 1 { + continue; + } + for &[a, b, c] in face.triangulation { + let base_idx = positions.len() as u32; + positions.push(self.vertices[cell[a] as usize]); + positions.push(self.vertices[cell[b] as usize]); + positions.push(self.vertices[cell[c] as usize]); + vertex_indices.push(cell[a] as usize); + vertex_indices.push(cell[b] as usize); + vertex_indices.push(cell[c] as usize); + cell_indices.push(cell_idx); + cell_indices.push(cell_idx); + cell_indices.push(cell_idx); + faces.push([base_idx, base_idx + 1, base_idx + 2]); + } + } +``` + +Read the surrounding context with `cargo run -- ` or by opening the file: confirm that `vertex_indices`, `cell_indices`, `positions`, `faces`, and `face_counts` are the bindings in scope and adjust names if any differ. The rest of the function (normals computation, quantity sampling) does not depend on cell type and stays unchanged. + +- [ ] **Step 8: Replace `generate_cell_index_per_triangle` and `generate_cell_index_per_triangle_with_culling`** + +Replace the body of each function with the same dispatch pattern. For `generate_cell_index_per_triangle` (around lines 763-794): + +```rust + fn generate_cell_index_per_triangle(&self) -> Vec { + let face_counts = self.compute_face_counts(); + let mut cell_indices = Vec::new(); + for (cell_idx, cell) in self.cells.iter().enumerate() { + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; + } + for _tri in face.triangulation { + cell_indices.push(cell_idx as u32); + } + } + } + cell_indices + } +``` + +And the `_with_culling` variant (around lines 796-829): + +```rust + fn generate_cell_index_per_triangle_with_culling(&self, planes: &[(Vec3, Vec3)]) -> Vec { + let face_counts = self.compute_face_counts_with_culling(planes); + let mut cell_indices = Vec::new(); + for (cell_idx, cell) in self.cells.iter().enumerate() { + if !self.is_cell_visible(cell, planes) { + continue; + } + let ct = self.cell_type(cell_idx); + for face in cell_data::face_data_for(ct) { + let key = cell_data::canonical_face_key(cell, face.polygon); + if face_counts.get(&key) != Some(&1) { + continue; + } + for _tri in face.triangulation { + cell_indices.push(cell_idx as u32); + } + } + } + cell_indices + } +``` + +- [ ] **Step 9: Delete the now-orphaned constants in `mod.rs`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, delete these blocks (now in `cell_data.rs`): + +- `canonical_face_key` (lines ~1410-1416) +- `TET_FACE_STENCIL` (line ~1419) +- `HEX_FACE_STENCIL` (lines ~1421-1429) + +The `VolumeMeshRenderGeometry` struct (lines ~1431-1440) stays. + +- [ ] **Step 10: Run tests** + +Run: `cargo test -p polyscope-structures -- --nocapture` +Expected: All tests PASS — both the new prism/pyramid/mixed tests and every existing tet/hex test (regression check). + +- [ ] **Step 11: Build the full workspace and run clippy** + +Run: `cargo build --workspace` +Expected: Builds clean. + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 12: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/mod.rs crates/polyscope-structures/src/volume_mesh/cell_data.rs +git commit -m "refactor(volume_mesh): dispatch face iteration via cell_data table" +``` + +--- + +### Task 5: Slice geometry for prism and pyramid + +**Files:** +- Modify: `crates/polyscope-structures/src/volume_mesh/slice_geometry.rs` (add two new functions) +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs` (update the `pub use slice_geometry::...` re-export) + +- [ ] **Step 1: Write the failing tests** + +Append to the `#[cfg(test)] mod tests` block in `slice_geometry.rs`: + +```rust +#[test] +fn test_slice_prism_through_middle() { + // Unit triangular prism z ∈ [0,1]; slice horizontally at y=... no, let's use z=0.5 + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + let result = slice_prism(verts, Vec3::new(0.0, 0.0, 0.5), Vec3::Z); + assert!(result.has_intersection()); + // Slice through the prism mid-height should yield a triangular cross-section + assert!(result.vertices.len() >= 3, "expected at least 3 verts, got {}", result.vertices.len()); + for v in &result.vertices { + assert!((v.z - 0.5).abs() < 1e-4, "vertex z={} should be 0.5", v.z); + } +} + +#[test] +fn test_slice_prism_no_intersection() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.0, 0.0, 1.0), + Vec3::new(1.0, 0.0, 1.0), + Vec3::new(0.5, 1.0, 1.0), + ]; + // Plane above the prism + let result = slice_prism(verts, Vec3::new(0.0, 0.0, 2.0), Vec3::Z); + assert!(!result.has_intersection()); +} + +#[test] +fn test_slice_pyramid_through_middle() { + // Unit pyramid: base z=0, apex at (0.5,0.5,1) + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let result = slice_pyramid(verts, Vec3::new(0.0, 0.0, 0.5), Vec3::Z); + assert!(result.has_intersection()); + // Slice at half height yields a square cross-section + assert!(result.vertices.len() >= 3); + for v in &result.vertices { + assert!((v.z - 0.5).abs() < 1e-4); + } +} + +#[test] +fn test_slice_pyramid_no_intersection() { + let verts = [ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(1.0, 1.0, 0.0), + Vec3::new(0.0, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let result = slice_pyramid(verts, Vec3::new(0.0, 0.0, -1.0), Vec3::Z); + assert!(!result.has_intersection()); +} +``` + +- [ ] **Step 2: Run tests to verify they fail** + +Run: `cargo test -p polyscope-structures slice_prism slice_pyramid -- --nocapture` +Expected: FAIL — `slice_prism` / `slice_pyramid` not yet defined. + +- [ ] **Step 3: Implement `slice_prism`** + +In `crates/polyscope-structures/src/volume_mesh/slice_geometry.rs`, append after `slice_hex` (after line 149): + +```rust +/// Slice a triangular prism by decomposing into 3 tetrahedra. +/// +/// # Arguments +/// * `vertices` - The 6 vertices of the prism (slots 0..2 = bottom tri, 3..5 = top tri) +/// * `plane_origin` - A point on the plane +/// * `plane_normal` - The plane normal (points toward kept geometry) +#[must_use] +pub fn slice_prism(vertices: [Vec3; 6], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { + // Symmetric 3-tet decomposition. Matches the upstream prism split for V0 + // in the bottom-left corner with no rotation (default case): tets + // {0, 5, 4, 3}, {0, 4, 5, 2}, {0, 4, 2, 1} + // (See cell_data::decompose_prism for the consistent variant; here we don't + // care about cross-cell consistency because slicing operates on isolated cells.) + let tet_indices = [ + [0usize, 5, 4, 3], + [0, 4, 5, 2], + [0, 4, 2, 1], + ]; + + let mut all_vertices = Vec::new(); + let mut all_interp = Vec::new(); + + for tet in &tet_indices { + let r = slice_tet( + vertices[tet[0]], + vertices[tet[1]], + vertices[tet[2]], + vertices[tet[3]], + plane_origin, + plane_normal, + ); + for (local_a, local_b, t) in r.interpolation { + let pa = tet[local_a as usize] as u32; + let pb = tet[local_b as usize] as u32; + all_interp.push((pa, pb, t)); + } + all_vertices.extend(r.vertices); + } + + merge_slice_vertices(&mut all_vertices, &mut all_interp); + if all_vertices.len() >= 3 { + order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); + } + + CellSliceResult { vertices: all_vertices, interpolation: all_interp } +} + +/// Slice a square pyramid by decomposing into 2 tetrahedra. +/// +/// # Arguments +/// * `vertices` - The 5 vertices of the pyramid (slots 0..3 = base quad, 4 = apex) +/// * `plane_origin` - A point on the plane +/// * `plane_normal` - The plane normal (points toward kept geometry) +#[must_use] +pub fn slice_pyramid(vertices: [Vec3; 5], plane_origin: Vec3, plane_normal: Vec3) -> CellSliceResult { + // Split base quad along diagonal 0-2 (consistent with cell_data::decompose_pyramid + // when p[0].min(p[2]) is the smaller pair). For isolated slicing the split + // choice doesn't affect correctness. + let tet_indices = [ + [0usize, 2, 4, 1], + [0, 4, 2, 3], + ]; + + let mut all_vertices = Vec::new(); + let mut all_interp = Vec::new(); + + for tet in &tet_indices { + let r = slice_tet( + vertices[tet[0]], + vertices[tet[1]], + vertices[tet[2]], + vertices[tet[3]], + plane_origin, + plane_normal, + ); + for (local_a, local_b, t) in r.interpolation { + let pa = tet[local_a as usize] as u32; + let pb = tet[local_b as usize] as u32; + all_interp.push((pa, pb, t)); + } + all_vertices.extend(r.vertices); + } + + merge_slice_vertices(&mut all_vertices, &mut all_interp); + if all_vertices.len() >= 3 { + order_polygon_vertices(&mut all_vertices, &mut all_interp, plane_normal); + } + + CellSliceResult { vertices: all_vertices, interpolation: all_interp } +} +``` + +- [ ] **Step 4: Re-export the new functions** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, find the existing line (around line 51): + +```rust +pub use slice_geometry::{CellSliceResult, slice_hex, slice_tet}; +``` + +and replace with: + +```rust +pub use slice_geometry::{CellSliceResult, slice_hex, slice_prism, slice_pyramid, slice_tet}; +``` + +- [ ] **Step 5: Run tests to verify they pass** + +Run: `cargo test -p polyscope-structures slice_prism slice_pyramid -- --nocapture` +Expected: PASS. + +Run the full slice geometry suite: + +Run: `cargo test -p polyscope-structures slice_geometry -- --nocapture` +Expected: All PASS. + +- [ ] **Step 6: Run clippy** + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 7: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/slice_geometry.rs crates/polyscope-structures/src/volume_mesh/mod.rs +git commit -m "feat(volume_mesh): slice_prism / slice_pyramid via tet decomposition" +``` + +--- + +### Task 6: Public constructors and registration functions + +**Files:** +- Modify: `crates/polyscope-structures/src/volume_mesh/mod.rs` (add `new_prism_mesh`, `new_pyramid_mesh`) +- Modify: `crates/polyscope/src/volume_mesh.rs` (add `register_prism_mesh`, `register_pyramid_mesh`) +- Modify: `crates/polyscope/src/lib.rs` if `register_*` functions are explicitly re-exported (verify by grep first) + +- [ ] **Step 1: Write the failing test** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, append: + +```rust +#[test] +fn test_new_prism_mesh_constructor() { + let verts = vec![ + Vec3::ZERO, Vec3::X, Vec3::Y, + Vec3::Z, Vec3::X + Vec3::Z, Vec3::Y + Vec3::Z, + ]; + let prisms = vec![[0u32, 1, 2, 3, 4, 5]]; + let mesh = VolumeMesh::new_prism_mesh("p", verts, prisms); + assert_eq!(mesh.num_cells(), 1); + assert_eq!(mesh.cell_type(0), VolumeCellType::Prism); +} + +#[test] +fn test_new_pyramid_mesh_constructor() { + let verts = vec![ + Vec3::ZERO, Vec3::X, Vec3::X + Vec3::Y, Vec3::Y, + Vec3::splat(0.5) + Vec3::Z, + ]; + let pyramids = vec![[0u32, 1, 2, 3, 4]]; + let mesh = VolumeMesh::new_pyramid_mesh("py", verts, pyramids); + assert_eq!(mesh.num_cells(), 1); + assert_eq!(mesh.cell_type(0), VolumeCellType::Pyramid); +} +``` + +- [ ] **Step 2: Run tests to verify they fail** + +Run: `cargo test -p polyscope-structures test_new_prism_mesh_constructor test_new_pyramid_mesh_constructor -- --nocapture` +Expected: FAIL — methods not defined. + +- [ ] **Step 3: Add `new_prism_mesh` and `new_pyramid_mesh`** + +In `crates/polyscope-structures/src/volume_mesh/mod.rs`, after `new_hex_mesh` (after line 175), insert: + +```rust + /// Creates a triangular-prism (wedge) mesh. + /// + /// Each prism has 6 vertices: slots 0..2 form the bottom triangle and + /// slots 3..5 form the top triangle (slot `i+3` should be the vertex + /// above slot `i`). Cells are stored as 8-index arrays with the last two + /// slots set to `u32::MAX` for sentinel detection. + pub fn new_prism_mesh( + name: impl Into, + vertices: Vec, + prisms: Vec<[u32; 6]>, + ) -> Self { + let cells: Vec<[u32; 8]> = prisms + .into_iter() + .map(|p| [p[0], p[1], p[2], p[3], p[4], p[5], u32::MAX, u32::MAX]) + .collect(); + Self::new(name, vertices, cells) + } + + /// Creates a square-pyramid mesh. + /// + /// Each pyramid has 5 vertices: slots 0..3 form the base quad (in CCW + /// order viewed from outside the cell) and slot 4 is the apex. Cells are + /// stored as 8-index arrays with the last three slots set to `u32::MAX`. + pub fn new_pyramid_mesh( + name: impl Into, + vertices: Vec, + pyramids: Vec<[u32; 5]>, + ) -> Self { + let cells: Vec<[u32; 8]> = pyramids + .into_iter() + .map(|p| [p[0], p[1], p[2], p[3], p[4], u32::MAX, u32::MAX, u32::MAX]) + .collect(); + Self::new(name, vertices, cells) + } +``` + +- [ ] **Step 4: Add the registration wrappers** + +In `crates/polyscope/src/volume_mesh.rs`, after `register_hex_mesh` (after line 69), insert: + +```rust +/// Registers a triangular-prism (wedge) mesh with polyscope. +/// +/// Each entry in `prisms` is 6 vertex indices: slots 0..2 = bottom triangle, +/// slots 3..5 = top triangle (with slot `i+3` directly above slot `i`). +pub fn register_prism_mesh( + name: impl Into, + vertices: Vec, + prisms: Vec<[u32; 6]>, +) -> VolumeMeshHandle { + let name = name.into(); + let mesh = VolumeMesh::new_prism_mesh(name.clone(), vertices, prisms); + + with_context_mut(|ctx| { + ctx.registry + .register(Box::new(mesh)) + .expect("failed to register prism mesh"); + ctx.update_extents(); + }); + + VolumeMeshHandle { name } +} + +/// Registers a square-pyramid mesh with polyscope. +/// +/// Each entry in `pyramids` is 5 vertex indices: slots 0..3 = base quad (CCW +/// from outside), slot 4 = apex. +pub fn register_pyramid_mesh( + name: impl Into, + vertices: Vec, + pyramids: Vec<[u32; 5]>, +) -> VolumeMeshHandle { + let name = name.into(); + let mesh = VolumeMesh::new_pyramid_mesh(name.clone(), vertices, pyramids); + + with_context_mut(|ctx| { + ctx.registry + .register(Box::new(mesh)) + .expect("failed to register pyramid mesh"); + ctx.update_extents(); + }); + + VolumeMeshHandle { name } +} +``` + +- [ ] **Step 5: Check if re-exports are needed** + +Run: `grep -n "register_tet_mesh\|register_hex_mesh" crates/polyscope/src/lib.rs` + +If the existing `register_tet_mesh` and `register_hex_mesh` appear in an explicit `pub use` list, add `register_prism_mesh` and `register_pyramid_mesh` to the same list. If they are auto-exported via `pub use crate::volume_mesh::*;` or similar wildcard, no change is needed. + +If a `pub use` exists like: +```rust +pub use crate::volume_mesh::{register_hex_mesh, register_tet_mesh, register_volume_mesh, VolumeMeshHandle, ...}; +``` +extend it to: +```rust +pub use crate::volume_mesh::{register_hex_mesh, register_prism_mesh, register_pyramid_mesh, register_tet_mesh, register_volume_mesh, VolumeMeshHandle, ...}; +``` + +- [ ] **Step 6: Run all tests** + +Run: `cargo test --workspace -- --nocapture` +Expected: All PASS. + +- [ ] **Step 7: Run clippy** + +Run: `cargo clippy --workspace -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 8: Commit** + +```bash +git add crates/polyscope-structures/src/volume_mesh/mod.rs crates/polyscope/src/volume_mesh.rs crates/polyscope/src/lib.rs +git commit -m "feat(volume_mesh): public register_prism_mesh / register_pyramid_mesh" +``` + +--- + +### Task 7: Visual demo example + +**Files:** +- Create: `examples/volume_mesh_mixed_cells_demo.rs` + +- [ ] **Step 1: Write the example** + +Look at `examples/volume_mesh_demo.rs` (13K, has the existing layout) to confirm the import style and structure. Then create `examples/volume_mesh_mixed_cells_demo.rs`: + +```rust +//! Demonstrates all four volume-mesh cell types side by side: tet, hex, +//! prism (wedge), and pyramid. Each cell is rendered as its own mesh so +//! the structure list shows them separately. + +use polyscope_rs::*; + +fn main() -> Result<()> { + init()?; + + // ----- Tet ----- + let tet_verts = vec![ + Vec3::new(0.0, 0.0, 0.0), + Vec3::new(1.0, 0.0, 0.0), + Vec3::new(0.5, 1.0, 0.0), + Vec3::new(0.5, 0.5, 1.0), + ]; + let tet = register_tet_mesh("tet", tet_verts, vec![[0, 1, 2, 3]]); + tet.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 1.0]); + + // ----- Hex (translated +3 in x) ----- + let dx = Vec3::new(3.0, 0.0, 0.0); + let hex_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 1.0, 1.0) + dx, + Vec3::new(0.0, 1.0, 1.0) + dx, + ]; + let hex = register_hex_mesh("hex", hex_verts, vec![[0, 1, 2, 3, 4, 5, 6, 7]]); + hex.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]); + + // ----- Prism (translated +6 in x) ----- + let dx = Vec3::new(6.0, 0.0, 0.0); + let prism_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(0.5, 1.0, 0.0) + dx, + Vec3::new(0.0, 0.0, 1.0) + dx, + Vec3::new(1.0, 0.0, 1.0) + dx, + Vec3::new(0.5, 1.0, 1.0) + dx, + ]; + let prism = register_prism_mesh("prism", prism_verts, vec![[0, 1, 2, 3, 4, 5]]); + prism.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0]); + + // ----- Pyramid (translated +9 in x) ----- + let dx = Vec3::new(9.0, 0.0, 0.0); + let pyr_verts = vec![ + Vec3::new(0.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 0.0, 0.0) + dx, + Vec3::new(1.0, 1.0, 0.0) + dx, + Vec3::new(0.0, 1.0, 0.0) + dx, + Vec3::new(0.5, 0.5, 1.0) + dx, + ]; + let pyr = register_pyramid_mesh("pyramid", pyr_verts, vec![[0, 1, 2, 3, 4]]); + pyr.add_vertex_scalar_quantity("z", vec![0.0, 0.0, 0.0, 0.0, 1.0]); + + show(); + Ok(()) +} +``` + +- [ ] **Step 2: Build the example to verify it compiles** + +Run: `cargo build --example volume_mesh_mixed_cells_demo` +Expected: Builds clean. + +- [ ] **Step 3: Run clippy on the example** + +Run: `cargo clippy --example volume_mesh_mixed_cells_demo -- -D warnings` +Expected: Zero warnings. + +- [ ] **Step 4: Run the demo manually and verify visually** + +Run: `cargo run --release --example volume_mesh_mixed_cells_demo` + +Visual checks: +- Four distinct shapes appear in a row along the x axis. +- The structure list on the right panel shows four entries: `tet`, `hex`, `prism`, `pyramid`. +- Each has a "z" scalar quantity selectable in its UI. +- Toggling visibility hides only that one cell. +- Hovering / clicking each cell highlights it (picking should work — confirms the per-triangle cell-index mapping is correct). + +If the visual check passes, proceed. If anything looks wrong (especially mismatched faces, gaps, or wrong shapes), debug before committing — the most likely culprit is a wrong vertex order in `PRISM_FACES` / `PYRAMID_FACES` in `cell_data.rs`. + +- [ ] **Step 5: Commit** + +```bash +git add examples/volume_mesh_mixed_cells_demo.rs +git commit -m "docs(examples): add mixed-cell-type volume mesh demo" +``` + +--- + +### Task 8: Update documentation + +**Files:** +- Modify: `docs/feature-status.md` +- Modify: `CHANGELOG.md` + +- [ ] **Step 1: Update `feature-status.md`** + +In `docs/feature-status.md`, find the Volume Mesh row in the Structures table (around line 12): + +```markdown +| Volume Mesh | Full | Full | Tet/hex, interior face detection, slice capping | +``` + +Update it to: + +```markdown +| Volume Mesh | Full | Full | Tet/hex/prism/pyramid, interior face detection, slice capping | +``` + +Find the `## Completed Features` section (around line 60) and add a new entry under the existing upstream-port lines (after the entry for inverse-view interpolation, line 79): + +```markdown +- [x] Volume Mesh prism + pyramid cell support (upstream PR #353, commit dcbaedb) +``` + +Find the `### Upstream Ports (Medium-Term)` subsection under Planned Work — no removal needed, the prism/pyramid item was not previously listed there. + +- [ ] **Step 2: Update `CHANGELOG.md`** + +In `CHANGELOG.md`, add a new entry at the top of the file (above `## [0.5.9]`): + +```markdown +## [Unreleased] + +### Added +- Volume Mesh now supports prism (wedge) and pyramid cell types in addition to + tetrahedra and hexahedra. New constructors: `register_prism_mesh`, + `register_pyramid_mesh`, `VolumeMesh::new_prism_mesh`, + `VolumeMesh::new_pyramid_mesh`. Mixed-cell meshes are supported by using the + 8-slot cell array with sentinel `u32::MAX` indices in unused slots + (matches upstream Polyscope PR #353). +- `slice_prism` and `slice_pyramid` helpers in + `polyscope_structures::volume_mesh::slice_geometry`. +- Example `volume_mesh_mixed_cells_demo` showing all four cell types. + +### Changed +- `VolumeMesh::cell_type` now classifies cells by sentinel count (0/2/3/4 → + Hex/Prism/Pyramid/Tet) instead of by `cell[4] == u32::MAX`. Mixed meshes + produced by external pipelines must place sentinels in the correct trailing + slots; tet meshes built with `new_tet_mesh` continue to work unchanged. +``` + +- [ ] **Step 3: Verify the docs build / render correctly** + +Run: `head -30 CHANGELOG.md` +Expected: New `## [Unreleased]` section sits cleanly above `## [0.5.9]`. + +- [ ] **Step 4: Commit** + +```bash +git add docs/feature-status.md CHANGELOG.md +git commit -m "docs: prism/pyramid cell support (upstream PR #353)" +``` + +--- + +## Self-Review + +### Spec coverage + +Upstream PR #353 adds: +1. PRISM and PYRAMID enum variants — Task 1 ✓ +2. Sentinel-count cell-type detection — Task 2 ✓ +3. Per-cell-type face stencils + face polygons + real-edge stencils — Task 1 (stencils + polygons); real-edge stencils are only used for wireframe rendering, which polyscope-rs does not currently expose for volume meshes; not required for parity at the geometry level. +4. `decomposePrism` / `decomposePyramid` — Task 3 ✓ +5. Slice integration through `computeTets` — Tasks 3 + 5 ✓ +6. Constructors taking mixed cells — Task 6 (`new_prism_mesh` / `new_pyramid_mesh` plus pre-existing `register_volume_mesh` for mixed) ✓ +7. UI integration / menu / picking — already cell-type-agnostic in polyscope-rs because dispatch goes through `cell_type()`; verified visually in Task 7. + +Edge cases left intentionally out of scope: +- Wireframe edges along "real" cell edges (would require porting the `realEdgeStencil`). Tracked separately if/when polyscope-rs adds volume-mesh edge rendering. +- Categorical-data isosurface on prism/pyramid (the polyscope-rs port does isosurface on volume **grids**, not on volume **meshes** with level-set quantities, so the upstream `activeLevelSetQuantity` path is not exercised the same way). + +### Placeholder scan + +No "TBD" / "implement later" / "similar to Task N" / "add appropriate handling" — every step has a complete code block or an exact command. + +### Type consistency + +- `VolumeCellType` variants `Tet`, `Hex`, `Prism`, `Pyramid` used consistently across Tasks 1–8. +- `decompose_prism` returns `[[u32; 4]; 3]` and `decompose_pyramid` returns `[[u32; 4]; 2]`; the dispatcher `decompose_cell_to_tets` wraps both in `Vec<[u32; 4]>`. Consistent. +- `slice_prism` takes `[Vec3; 6]`, `slice_pyramid` takes `[Vec3; 5]` (matching upstream and the new constructors). +- `register_prism_mesh(prisms: Vec<[u32; 6]>)` and `register_pyramid_mesh(pyramids: Vec<[u32; 5]>)` mirror the inner `new_prism_mesh` / `new_pyramid_mesh` signatures. +- `canonical_face_key` signature: `fn(&[u32; 8], &[usize]) -> [u32; 4]`. Same call site uses `face.polygon` (`&'static [usize]`) — type-compatible. + +--- + +## Execution Handoff + +**Plan complete and saved to `docs/plans/2026-05-19-volume-mesh-prism-pyramid.md`. Two execution options:** + +**1. Subagent-Driven (recommended)** — I dispatch a fresh subagent per task, review between tasks, fast iteration. + +**2. Inline Execution** — Execute tasks in this session using executing-plans, batch execution with checkpoints. + +**Which approach?** From f2cf90ff51f41890530a07e8ff03d1228d4e929e Mon Sep 17 00:00:00 2001 From: "ZM.TreeWSL" Date: Tue, 19 May 2026 11:36:57 +0200 Subject: [PATCH 14/14] fix(doc): escape array indexing in decompose_pyramid doc comment Rustdoc parsed p[0], p[2], p[1], p[3] as broken intra-doc links, failing the workspace doc build with -D warnings. Wrap the expressions in backticks so rustdoc treats them as code spans. --- crates/polyscope-structures/src/volume_mesh/cell_data.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/polyscope-structures/src/volume_mesh/cell_data.rs b/crates/polyscope-structures/src/volume_mesh/cell_data.rs index 7fe15ac..cc8fcd8 100644 --- a/crates/polyscope-structures/src/volume_mesh/cell_data.rs +++ b/crates/polyscope-structures/src/volume_mesh/cell_data.rs @@ -264,7 +264,7 @@ fn rotate_prism_in_place(p: &mut [u32; 6], rot: usize) { } /// Decomposes a square pyramid into 2 tetrahedra by splitting the base quad -/// along the diagonal containing the smaller of {p[0], p[2]} vs {p[1], p[3]}. +/// along the diagonal containing the smaller of `{p[0], p[2]}` vs `{p[1], p[3]}`. /// Consistent split ensures adjacent cells tessellate the shared face the same way. #[must_use] pub fn decompose_pyramid(cell: &[u32; 8]) -> [[u32; 4]; 2] {