From 827d3b22aa1b2576ad314c304485879f31f14a82 Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Sat, 9 Aug 2025 11:19:27 -0500 Subject: [PATCH 01/15] Start on supporting prisms and pyramids. Edges still need work --- include/polyscope/types.h | 2 +- include/polyscope/volume_mesh.h | 2 + src/volume_mesh.cpp | 106 ++++++++++++++++++++++++++++++-- test/src/volume_mesh_test.cpp | 55 ++++++++++++++--- 4 files changed, 152 insertions(+), 13 deletions(-) diff --git a/include/polyscope/types.h b/include/polyscope/types.h index 894ca864..841ed0ed 100644 --- a/include/polyscope/types.h +++ b/include/polyscope/types.h @@ -22,7 +22,7 @@ enum class MeshShadeStyle { Smooth = 0, Flat, TriFlat }; enum class MeshSelectionMode { Auto = 0, VerticesOnly, FacesOnly }; enum class CurveNetworkElement { NODE = 0, EDGE }; enum class VolumeMeshElement { VERTEX = 0, EDGE, FACE, CELL }; -enum class VolumeCellType { TET = 0, HEX }; +enum class VolumeCellType { TET = 0, HEX, PRISM, PYRAMID }; enum class VolumeGridElement { NODE = 0, CELL }; enum class IsolineStyle { Stripe = 0, Contour }; diff --git a/include/polyscope/volume_mesh.h b/include/polyscope/volume_mesh.h index 83e4d30a..d9d1bb8a 100644 --- a/include/polyscope/volume_mesh.h +++ b/include/polyscope/volume_mesh.h @@ -262,6 +262,8 @@ class VolumeMesh : public QuantityStructure { // clang-format off static const std::vector>> stencilTet; static const std::vector>> stencilHex; + static const std::vector>> stencilPrism; + static const std::vector>> stencilPyramid; static const std::array, 8> rotationMap; static const std::array, 6>, 4> diagonalMap; diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index 3e4bf589..0f107525 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -93,8 +93,87 @@ const std::vector>> VolumeMesh::stencilHex = {{6,2,3}, {6,3,7}}, {{7,4,5}, {7,5,6}}, }; + + const std::vector>> VolumeMesh::stencilPrism = + { + {{0,2,1}}, + {{0,1,4}, {0,4,3}}, + {{1,5,4}, {1,2,5}}, + {{0,3,5}, {0,5,2}}, + {{3,4,5}}, + }; + + const std::vector>> VolumeMesh::stencilPyramid = + { + {{0,3,2}, {0,2,1}}, + {{0,1,4}}, + {{1,2,4}}, + {{2,3,4}}, + {{3,0,4}}, + }; // clang-format on + +namespace detail { +using Tet = std::array; +using Prism = std::array; +using Pyramid = std::array; +using Cell = std::array; +using Tets = std::vector; + +Prism rotatePrism(const Prism& p, int rot) { + static constexpr int bottomRot[3][3] = { + {0, 1, 2}, // identity + {1, 2, 0}, // rot1 + {2, 0, 1} // rot2 + }; + static constexpr int topRot[3][3] = {{3, 4, 5}, {4, 5, 3}, {5, 3, 4}}; + Prism r; + for (int i = 0; i < 3; ++i) r[i] = p[bottomRot[rot][i]]; + for (int i = 0; i < 3; ++i) r[i + 3] = p[topRot[rot][i]]; + return r; +} +Tets decomposePrism(const Cell& cell) { + Prism p = {cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]}; + int minIdx = std::min_element(p.begin(), p.end()) - p.begin(); + int rot = 0; + // If smallest vertex is in bottom triangle + if (minIdx < 3) { + // rotate bottom triangle so min vertex is V0 + int bPos = (minIdx == 0 ? 0 : (minIdx == 1 ? 2 : 1)); // mapping from input idx to rot + rot = bPos; + } else { + // smallest vertex in top triangle, reflect wedge + int topPos = minIdx - 3; + rot = (topPos == 0 ? 0 : (topPos == 1 ? 2 : 1)); + // swap top/bottom + std::swap(p[0], p[3]); + std::swap(p[1], p[4]); + std::swap(p[2], p[5]); + } + // Rotate prism to have V0 in bottom left corner + p = rotatePrism(p, rot); + // Now, decide how to split the quad opposite V0 + if (std::min(p[2], p[4]) < std::min(p[1], p[5])) { + // Split along diagonal 2-4 + return {{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 { + // Split along diagonal 1-5 + return {{p[0], p[5], p[4], p[3]}, {p[0], p[1], p[5], p[2]}, {p[0], p[5], p[1], p[4]}}; + } +} +Tets decomposePyramid(const Cell& cell) { + Pyramid p = {cell[0], cell[1], cell[2], cell[3], cell[4]}; + if (std::min(p[0], p[2]) < std::min(p[1], p[3])) { + // Split along diagonal 0-2 + return {{p[0], p[2], p[4], p[1]}, {p[0], p[4], p[2], p[3]}}; + } else { + // Split along diagonal 1-3 + return {{p[1], p[3], p[4], p[2]}, {p[1], p[4], p[3], p[0]}}; + } +} +} // namespace detail + VolumeMesh::VolumeMesh(std::string name, const std::vector& vertexPositions_, const std::vector>& cellIndices_) : QuantityStructure(name, typeName()), @@ -264,14 +343,26 @@ void VolumeMesh::computeTets() { case VolumeCellType::TET: tetCount += 1; break; + case VolumeCellType::PRISM: + tetCount += 3; + break; + case VolumeCellType::PYRAMID: + tetCount += 2; + break; } } - // Mark each edge as real or not (in the original mesh) - std::vector> realEdges; + // Each hex can make up to 6 tets tets.resize(tetCount); - realEdges.resize(tetCount); size_t tetIdx = 0; + auto addTets = [&](const detail::Tets& newTets) { + for (const auto& tet : newTets) { + for (size_t i = 0; i < 4; i++) { + tets[tetIdx][i] = tet[i]; + } + tetIdx++; + } + }; for (size_t iC = 0; iC < nCells(); iC++) { switch (cellType(iC)) { case VolumeCellType::HEX: { @@ -335,13 +426,20 @@ void VolumeMesh::computeTets() { } break; } - case VolumeCellType::TET: + case VolumeCellType::TET: { for (size_t i = 0; i < 4; i++) { tets[tetIdx][i] = cells[iC][i]; } tetIdx++; break; } + case VolumeCellType::PRISM: + addTets(detail::decomposePrism(cells[iC])); + break; + case VolumeCellType::PYRAMID: + addTets(detail::decomposePyramid(cells[iC])); + break; + } } } diff --git a/test/src/volume_mesh_test.cpp b/test/src/volume_mesh_test.cpp index 1a81a786..9aa18760 100644 --- a/test/src/volume_mesh_test.cpp +++ b/test/src/volume_mesh_test.cpp @@ -12,12 +12,12 @@ TEST_F(PolyscopeTest, ShowVolumeMesh) { // Tets only std::vector tet_verts = { {0, 0, 0}, - {0, 0, 1}, + {1, 0, 0}, {0, 1, 0}, - {0, 1, 1}, + {0, 0, 1}, }; std::vector> tet_cells = { - {0,1,2,4} + {0,1,2,3} }; polyscope::registerTetMesh("tet", tet_verts, tet_cells); @@ -25,13 +25,13 @@ TEST_F(PolyscopeTest, ShowVolumeMesh) { // Hexes only std::vector hex_verts = { {0, 0, 0}, - {0, 0, 1}, - {0, 1, 0}, - {0, 1, 1}, {1, 0, 0}, - {1, 0, 1}, {1, 1, 0}, + {0, 1, 0}, + {0, 0, 1}, + {1, 0, 1}, {1, 1, 1}, + {0, 1, 1}, }; std::vector> hex_cells = { {0,1,2,3,4,5,6,7}, @@ -55,7 +55,7 @@ TEST_F(PolyscopeTest, ShowVolumeMesh) { // Mixed elements, shared array std::vector> combined_cells = { - {0, 1, 3, 4, -1, -1, -1, -1}, + {0, 1, 2, 3, -1, -1, -1, -1}, {4, 5, 6, 7, 8, 9, 10, 11}, }; polyscope::registerVolumeMesh("tet hex mix combined", combined_verts, combined_cells); @@ -63,6 +63,45 @@ TEST_F(PolyscopeTest, ShowVolumeMesh) { polyscope::removeAllStructures(); } + +TEST_F(PolyscopeTest, ShowVolumeMeshHexWedgePyramidTet) { + // clang-format off + std::vector vertices = {// Base hex vertices + {0, 0, 0}, // V0 + {1, 0, 0}, // V1 + {1, 1, 0}, // V2 + {0, 1, 0}, // V3 + {0, 0, 1}, // V4 + {1, 0, 1}, // V5 + {1, 1, 1}, // V6 + {0, 1, 1}, // V7 + // Top Prism Vertices + {0.0, 0.5, 1.5}, // V8 + {1.0, 0.5, 1.5}, // V9 + // Side Prism Vertices + {1.5, 0.5, 0.0}, // V10 + {1.5, 0.5, 1.0}, // V11 + // Bottom Pyramid Vertex + {0.5, 0.5, -0.5} // V12 + }; + // clang-format on + std::vector> cells = { + // Base Hex cell + {0, 1, 2, 3, 4, 5, 6, 7}, + // Top Prism cell + {4, 7, 8, 5, 6, 9, -1, -1}, + // Side Prism cell + {1, 10, 2, 5, 11, 6, -1, -1}, + // Bottom Pyramid cell + {0, 3, 2, 1, 12, -1, -1, -1}, + // Tet Connecting side and top prisms + {5, 11, 6, 9, -1, -1, -1, -1}, + }; + polyscope::registerVolumeMesh("hex wedge pyramid tet", vertices, cells); + polyscope::show(); + polyscope::removeAllStructures(); +} + TEST_F(PolyscopeTest, VolumeMeshUpdatePositions) { std::vector verts; std::vector> cells; From dcf540013be89f5c7d8a035cec8c0f4e53221915 Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Sat, 9 Aug 2025 15:08:07 -0500 Subject: [PATCH 02/15] Fix --- src/volume_mesh.cpp | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index 0f107525..c3ac6b22 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -96,16 +96,16 @@ const std::vector>> VolumeMesh::stencilHex = const std::vector>> VolumeMesh::stencilPrism = { - {{0,2,1}}, - {{0,1,4}, {0,4,3}}, - {{1,5,4}, {1,2,5}}, - {{0,3,5}, {0,5,2}}, - {{3,4,5}}, + {{0,1,2}}, + {{0,2,3}, {3,2,5}}, + {{5,2,1}, {5,1,4}}, + {{4,1,3}, {3,1,0}}, + {{5,4,3}}, }; const std::vector>> VolumeMesh::stencilPyramid = { - {{0,3,2}, {0,2,1}}, + {{0,3,1}, {1,3,2}}, {{0,1,4}}, {{1,2,4}}, {{2,3,4}}, @@ -833,6 +833,10 @@ const std::vector>>& VolumeMesh::cellStencil(V return stencilTet; case VolumeCellType::HEX: return stencilHex; + case VolumeCellType::PRISM: + return stencilPrism; + case VolumeCellType::PYRAMID: + return stencilPyramid; } // unreachable @@ -1073,11 +1077,18 @@ void VolumeMesh::recomputeGeometryIfPopulated() { } VolumeCellType VolumeMesh::cellType(size_t i) const { - bool isHex = cells[i][4] < INVALID_IND_32; - if (isHex) { + auto sentinelIndices = std::count(cells[i].begin(), cells[i].end(), INVALID_IND_32); + switch (sentinelIndices) { + case 0: // no sentinel indices, must be a hex return VolumeCellType::HEX; - } else { + case 2: // two sentinel indices, must be a prism + return VolumeCellType::PRISM; + case 3: // three sentinel indices, must be a pyramid + return VolumeCellType::PYRAMID; + case 4: // four sentinel indices, must be a tet return VolumeCellType::TET; + default: + throw std::runtime_error("VolumeMesh::cellType: Invalid number of sentinel indices in cell"); } }; From 73ab829c99ebaade357bd80a1ff631976aeaa4c0 Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Sat, 9 Aug 2025 15:09:45 -0500 Subject: [PATCH 03/15] Some test cleanup --- test/src/volume_mesh_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/src/volume_mesh_test.cpp b/test/src/volume_mesh_test.cpp index 9aa18760..781ed181 100644 --- a/test/src/volume_mesh_test.cpp +++ b/test/src/volume_mesh_test.cpp @@ -97,8 +97,8 @@ TEST_F(PolyscopeTest, ShowVolumeMeshHexWedgePyramidTet) { // Tet Connecting side and top prisms {5, 11, 6, 9, -1, -1, -1, -1}, }; - polyscope::registerVolumeMesh("hex wedge pyramid tet", vertices, cells); - polyscope::show(); + polyscope::registerVolumeMesh("hex prism pyramid tet", vertices, cells); + polyscope::show(3); polyscope::removeAllStructures(); } From be3c7d97c53ba77e878ca0e5e47698d5e9290d2e Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Mon, 11 Aug 2025 08:40:52 -0500 Subject: [PATCH 04/15] Fix for the edges --- test/src/volume_mesh_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/src/volume_mesh_test.cpp b/test/src/volume_mesh_test.cpp index 781ed181..3a12f531 100644 --- a/test/src/volume_mesh_test.cpp +++ b/test/src/volume_mesh_test.cpp @@ -98,7 +98,7 @@ TEST_F(PolyscopeTest, ShowVolumeMeshHexWedgePyramidTet) { {5, 11, 6, 9, -1, -1, -1, -1}, }; polyscope::registerVolumeMesh("hex prism pyramid tet", vertices, cells); - polyscope::show(3); + polyscope::show(); polyscope::removeAllStructures(); } From a3f97bfa60a5e2aeb8d5af02c25a0645c7914700 Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Mon, 11 Aug 2025 09:01:53 -0500 Subject: [PATCH 05/15] Some more work on the edges --- include/polyscope/volume_mesh.h | 7 +- src/volume_mesh.cpp | 124 ++++++++++++++++++++++++-------- 2 files changed, 102 insertions(+), 29 deletions(-) diff --git a/include/polyscope/volume_mesh.h b/include/polyscope/volume_mesh.h index d9d1bb8a..70ca84df 100644 --- a/include/polyscope/volume_mesh.h +++ b/include/polyscope/volume_mesh.h @@ -180,7 +180,7 @@ class VolumeMesh : public QuantityStructure { void fillGeometryBuffers(render::ShaderProgram& p); void fillSliceGeometryBuffers(render::ShaderProgram& p); static const std::vector>>& cellStencil(VolumeCellType type); - + static const std::vector>& cellFaces(VolumeCellType type); // Slice plane listeners std::vector volumeSlicePlaneListeners; void addSlicePlaneListener(polyscope::SlicePlane* sp); @@ -264,6 +264,11 @@ class VolumeMesh : public QuantityStructure { static const std::vector>> stencilHex; static const std::vector>> stencilPrism; static const std::vector>> stencilPyramid; + static const std::vector> facesTet; + static const std::vector> facesHex; + static const std::vector> facesPrism; + static const std::vector> facesPyramid; + static const std::array, 8> rotationMap; static const std::array, 6>, 4> diagonalMap; diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index c3ac6b22..2463c250 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -30,7 +30,13 @@ const std::vector>> VolumeMesh::stencilTet = {{0,3,2}}, {{1,2,3}}, }; - +const std::vector> VolumeMesh::facesTet = +{ + {0,2,1}, + {0,1,3}, + {0,3,2}, + {1,2,3} +}; // Indirection to place vertex 0 always in the bottom left corner const std::array, 8> VolumeMesh::rotationMap = {{ @@ -94,23 +100,51 @@ const std::vector>> VolumeMesh::stencilHex = {{7,4,5}, {7,5,6}}, }; + const std::vector> VolumeMesh::facesHex = + { + {2,1,0,3}, + {4,0,1,5}, + {5,1,2,6}, + {7,3,0,4}, + {6,2,3,7}, + {7,4,5,6} + }; + const std::vector>> VolumeMesh::stencilPrism = { {{0,1,2}}, - {{0,2,3}, {3,2,5}}, - {{5,2,1}, {5,1,4}}, - {{4,1,3}, {3,1,0}}, - {{5,4,3}}, + {{0,4,1}, {0,3,4}}, + {{1,5,4}, {1,2,5}}, + {{3,0,5}, {0,2,5}}, + {{3,5,4}}, + }; + + const std::vector> VolumeMesh::facesPrism = + { + {0,1,2}, + {0,3,4,1}, + {1,4,5,2}, + {0,2,5,3}, + {3,5,4} }; const std::vector>> VolumeMesh::stencilPyramid = { - {{0,3,1}, {1,3,2}}, + {{0,3,2}, {0,2,1}}, {{0,1,4}}, {{1,2,4}}, {{2,3,4}}, {{3,0,4}}, }; + + const std::vector> VolumeMesh::facesPyramid = + { + {0,1,2,3}, + {0,1,4}, + {1,2,4}, + {2,3,4}, + {3,0,4} + }; // clang-format on @@ -745,7 +779,6 @@ void VolumeMesh::fillGeometryBuffers(render::ShaderProgram& p) { } void VolumeMesh::computeConnectivityData() { - // NOTE: If we were to fill buffers naively via a loop over cells, we get pretty bad z-fighting artifacts where // interior edges ever-so-slightly show through the exterior boundary (more generally, any place 3 faces meet at an // edge, which happens everywhere in a tet mesh). @@ -753,6 +786,10 @@ void VolumeMesh::computeConnectivityData() { // To mitigate this issue, we fill the buffer such that all exterior faces come first, then all interior faces, so // that exterior faces always win depth ties. This doesn't totally eliminate the problem, but greatly improves the // most egregious cases. + auto makeEdgeKey = [](uint32_t a, uint32_t b) -> std::pair { + if (a > b) std::swap(a, b); + return {a, b}; + }; // == Allocate buffers triangleVertexInds.data.clear(); @@ -761,8 +798,6 @@ void VolumeMesh::computeConnectivityData() { triangleFaceInds.data.resize(3 * nFacesTriangulation()); triangleCellInds.data.clear(); triangleCellInds.data.resize(3 * nFacesTriangulation()); - triangleCellInds.data.clear(); - triangleCellInds.data.resize(3 * nFacesTriangulation()); baryCoord.data.clear(); baryCoord.data.resize(3 * nFacesTriangulation()); edgeIsReal.data.clear(); @@ -770,58 +805,77 @@ void VolumeMesh::computeConnectivityData() { faceType.data.clear(); faceType.data.resize(nFaces()); - size_t iF = 0; + size_t iF = 0; // face counter size_t iFront = 0; size_t iBack = nFacesTriangulation() - 1; + + // Loop over all cells for (size_t iC = 0; iC < nCells(); iC++) { const std::array& cell = cells[iC]; VolumeCellType cellT = cellType(iC); + const auto& faceTris = cellStencil(cellT); // triangulated faces + const auto& faceVerts = cellFaces(cellT); // original polygon faces + // Loop over all faces of the cell - for (const std::vector>& face : cellStencil(cellT)) { + for (size_t f = 0; f < faceTris.size(); f++) { - // Loop over the face's triangulation - for (size_t j = 0; j < face.size(); j++) { - const std::array& tri = face[j]; + const auto& triList = faceTris[f]; + const auto& poly = faceVerts[f]; - // Enumerate exterior faces in the front of the draw buffer, and interior faces in the back. - // (see note above) + // Build a set of "real" edges from the polygon + std::set> realEdges; + for (size_t p = 0; p < poly.size(); p++) { + realEdges.insert(makeEdgeKey(cell[poly[p]], cell[poly[(p + 1) % poly.size()]])); + } + + // Loop over triangles in the triangulation of this face + for (size_t j = 0; j < triList.size(); j++) { + const auto& tri = triList[j]; + + // Decide where to store (front for exterior, back for interior) size_t iData; if (faceIsInterior[iF]) { - iData = iBack; - iBack--; + iData = iBack--; } else { - iData = iFront; - iFront++; + iData = iFront++; } + // Store triangle vertices for (size_t k = 0; k < 3; k++) { triangleVertexInds.data[3 * iData + k] = cell[tri[k]]; } + + // Face & cell indices for (size_t k = 0; k < 3; k++) triangleFaceInds.data[3 * iData + k] = iF; for (size_t k = 0; k < 3; k++) triangleCellInds.data[3 * iData + k] = iC; + // Barycentric coords baryCoord.data[3 * iData + 0] = glm::vec3{1., 0., 0.}; baryCoord.data[3 * iData + 1] = glm::vec3{0., 1., 0.}; baryCoord.data[3 * iData + 2] = glm::vec3{0., 0., 1.}; - glm::vec3 edgeRealV{0., 1., 0.}; - if (j == 0) edgeRealV.x = 1.; - if (j + 1 == face.size()) edgeRealV.z = 1.; - for (int k = 0; k < 3; k++) edgeIsReal.data[3 * iData + k] = edgeRealV; + // Mark edges as real or not + glm::vec3 realFlag(0.0f); + for (size_t k = 0; k < 3; k++) { + bool isReal = realEdges.count(makeEdgeKey(cell[tri[k]], cell[tri[(k + 1) % 3]])) > 0; + realFlag[k] = isReal ? 1.0f : 0.0f; + edgeIsReal.data[3 * iData + k] = realFlag; + } + for (int k = 0; k < 3; k++) edgeIsReal.data[3 * iData + k] = realFlag; } - float faceTypeFloat = faceIsInterior[iF] ? 1. : 0.; - for (int k = 0; k < 3; k++) faceType.data[iF] = faceTypeFloat; + // Face type: 1 for interior, 0 for exterior + faceType.data[iF] = faceIsInterior[iF] ? 1.f : 0.f; iF++; } } + // Mark buffers updated triangleVertexInds.markHostBufferUpdated(); triangleFaceInds.markHostBufferUpdated(); triangleCellInds.markHostBufferUpdated(); - triangleCellInds.markHostBufferUpdated(); baryCoord.markHostBufferUpdated(); edgeIsReal.markHostBufferUpdated(); faceType.markHostBufferUpdated(); @@ -838,11 +892,25 @@ const std::vector>>& VolumeMesh::cellStencil(V case VolumeCellType::PYRAMID: return stencilPyramid; } - // unreachable return stencilTet; } +const std::vector>& VolumeMesh::cellFaces(VolumeCellType type) { + switch (type) { + case VolumeCellType::TET: + return facesTet; + case VolumeCellType::HEX: + return facesHex; + case VolumeCellType::PRISM: + return facesPrism; + case VolumeCellType::PYRAMID: + return facesPyramid; + } + // unreachable + return facesTet; +} + void VolumeMesh::computeFaceNormals() { From 1db6fee2d0cba8343421d620ed59457161a9bd93 Mon Sep 17 00:00:00 2001 From: Matt Sutton Date: Mon, 11 Aug 2025 09:02:26 -0500 Subject: [PATCH 06/15] Set the number of frames on show --- test/src/volume_mesh_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/src/volume_mesh_test.cpp b/test/src/volume_mesh_test.cpp index 3a12f531..781ed181 100644 --- a/test/src/volume_mesh_test.cpp +++ b/test/src/volume_mesh_test.cpp @@ -98,7 +98,7 @@ TEST_F(PolyscopeTest, ShowVolumeMeshHexWedgePyramidTet) { {5, 11, 6, 9, -1, -1, -1, -1}, }; polyscope::registerVolumeMesh("hex prism pyramid tet", vertices, cells); - polyscope::show(); + polyscope::show(3); polyscope::removeAllStructures(); } From 62531de3d600ac33f1549ee6e786ea023db48bbe Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 17:46:09 -0700 Subject: [PATCH 07/15] refactor and clean up computeTets() --- src/volume_mesh.cpp | 296 ++++++++++++++++++-------------------------- 1 file changed, 121 insertions(+), 175 deletions(-) diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index b4fbd77d..8927572a 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -37,55 +37,6 @@ const std::vector> VolumeMesh::facesTet = {0,3,2}, {1,2,3} }; -// Indirection to place vertex 0 always in the bottom left corner -const std::array, 8> VolumeMesh::rotationMap = - {{ - {0, 1, 2, 3, 4, 5, 7, 6}, - {1, 0, 4, 5, 2, 3, 6, 7}, - {2, 1, 5, 6, 3, 0, 7, 4}, - {3, 0, 1, 2, 7, 4, 6, 5}, - {4, 0, 3, 7, 5, 1, 6, 2}, - {5, 1, 0, 4, 7, 2, 6, 3}, - {7, 3, 2, 6, 4, 0, 5, 1}, - {6, 2, 1, 5, 7, 3, 4, 0} - }}; - -// Map indirected cube to tets -const std::array, 6>, 4> VolumeMesh::diagonalMap = - {{ - {{ - {0, 1, 2, 5}, - {0, 2, 6, 5}, - {0, 2, 3, 6}, - {0, 5, 6, 4}, - {2, 6, 5, 7}, - {0, 0, 0, 0} - }}, - {{ - {0, 5, 6, 4}, - {0, 1, 6, 5}, - {1, 7, 6, 5}, - {0, 6, 2, 3}, - {0, 6, 1, 2}, - {1, 6, 7, 2} - }}, - {{ - {0, 4, 5, 7}, - {0, 3, 6, 7}, - {0, 6, 4, 7}, - {0, 1, 2, 5}, - {0, 3, 7, 2}, - {0, 7, 5, 2} - }}, - {{ - {0, 2, 3, 7}, - {0, 3, 6, 7}, - {0, 6, 4, 7}, - {0, 5, 7, 4}, - {1, 5, 7, 0}, - {1, 7, 2, 0} - }} - }}; const std::vector>> VolumeMesh::stencilHex = @@ -206,6 +157,121 @@ Tets decomposePyramid(const Cell& cell) { return {{p[1], p[3], p[4], p[2]}, {p[1], p[4], p[3], p[0]}}; } } +Tets decomposeHex(const Cell& cell) { + // The challenge is to make sure we tet-decompse adjacent cells in a consistent way, so there are not little gaps + // between mismatched tetrahedralizations, for non-planar faces. This approach is adapted from + // https://www.researchgate.net/profile/Julien-Dompierre/publication/221561839_How_to_Subdivide_Pyramids_Prisms_and_Hexahedra_into_Tetrahedra/links/0912f509c0b7294059000000/How-to-Subdivide-Pyramids-Prisms-and-Hexahedra-into-Tetrahedra.pdf?origin=publication_detail + // It's a bit hard to look at but it works. + // The resulting tet counts are either 5 or 6 per hex. Like the rest of this class, this function only really works + // for convex cells, there is no attempt to handle the nonconvex case. + + // clang-format off + + // Indirection to place vertex 0 always in the bottom left corner + const std::array, 8> rotationMap = {{{0, 1, 2, 3, 4, 5, 7, 6}, + {1, 0, 4, 5, 2, 3, 6, 7}, + {2, 1, 5, 6, 3, 0, 7, 4}, + {3, 0, 1, 2, 7, 4, 6, 5}, + {4, 0, 3, 7, 5, 1, 6, 2}, + {5, 1, 0, 4, 7, 2, 6, 3}, + {7, 3, 2, 6, 4, 0, 5, 1}, + {6, 2, 1, 5, 7, 3, 4, 0}}}; + + + // Map indirected cube to tets + const std::array, 6>, 4> diagonalMap = + {{ + {{ + {0, 1, 2, 5}, + {0, 2, 6, 5}, + {0, 2, 3, 6}, + {0, 5, 6, 4}, + {2, 6, 5, 7}, + {0, 0, 0, 0} + }}, + {{ + {0, 5, 6, 4}, + {0, 1, 6, 5}, + {1, 7, 6, 5}, + {0, 6, 2, 3}, + {0, 6, 1, 2}, + {1, 6, 7, 2} + }}, + {{ + {0, 4, 5, 7}, + {0, 3, 6, 7}, + {0, 6, 4, 7}, + {0, 1, 2, 5}, + {0, 3, 7, 2}, + {0, 7, 5, 2} + }}, + {{ + {0, 2, 3, 7}, + {0, 3, 6, 7}, + {0, 6, 4, 7}, + {0, 5, 7, 4}, + {1, 5, 7, 0}, + {1, 7, 2, 0} + }} + }}; + // clang-format on + + + std::array sortedNumbering; + std::iota(sortedNumbering.begin(), sortedNumbering.end(), 0); + std::sort(sortedNumbering.begin(), sortedNumbering.end(), + [&cell](size_t a, size_t b) -> bool { return cell[a] < cell[b]; }); + std::array rotatedNumbering; + std::copy(rotationMap[sortedNumbering[0]].begin(), rotationMap[sortedNumbering[0]].end(), rotatedNumbering.begin()); + size_t n = 0; + size_t diagCount = 0; + auto checkDiagonal = [&cell, &rotatedNumbering](size_t a1, size_t a2, size_t b1, size_t b2) { + return (cell[rotatedNumbering[a1]] < cell[rotatedNumbering[b1]] && + cell[rotatedNumbering[a1]] < cell[rotatedNumbering[b2]]) || + (cell[rotatedNumbering[a2]] < cell[rotatedNumbering[b1]] && + cell[rotatedNumbering[a2]] < cell[rotatedNumbering[b2]]); + }; + if (checkDiagonal(1, 7, 2, 5)) { + n += 4; + diagCount++; + } + if (checkDiagonal(3, 7, 2, 6)) { + n += 2; + diagCount++; + } + if (checkDiagonal(4, 7, 5, 6)) { + n += 1; + diagCount++; + } + if (n == 1 || n == 6) { + size_t temp = rotatedNumbering[1]; + rotatedNumbering[1] = rotatedNumbering[4]; + rotatedNumbering[4] = rotatedNumbering[3]; + rotatedNumbering[3] = temp; + temp = rotatedNumbering[5]; + rotatedNumbering[5] = rotatedNumbering[6]; + rotatedNumbering[6] = rotatedNumbering[2]; + rotatedNumbering[2] = temp; + } else if (n == 2 || n == 5) { + size_t temp = rotatedNumbering[1]; + rotatedNumbering[1] = rotatedNumbering[3]; + rotatedNumbering[3] = rotatedNumbering[4]; + rotatedNumbering[4] = temp; + temp = rotatedNumbering[5]; + rotatedNumbering[5] = rotatedNumbering[2]; + rotatedNumbering[2] = rotatedNumbering[6]; + rotatedNumbering[6] = temp; + } + size_t tetCount = (diagCount == 0 ? 5 : 6); + std::array, 6> tetMap = diagonalMap[diagCount]; + Tets result(tetCount); + for (size_t k = 0; k < tetCount; k++) { + for (size_t i = 0; i < 4; i++) { + result[k][i] = cell[rotatedNumbering[tetMap[k][i]]]; + } + } + return result; +} } // namespace detail VolumeMesh::VolumeMesh(std::string name, const std::vector& vertexPositions_, @@ -335,138 +401,18 @@ void VolumeMesh::computeCounts() { void VolumeMesh::computeTets() { - // Algorithm from - // https://www.researchgate.net/profile/Julien-Dompierre/publication/221561839_How_to_Subdivide_Pyramids_Prisms_and_Hexahedra_into_Tetrahedra/links/0912f509c0b7294059000000/How-to-Subdivide-Pyramids-Prisms-and-Hexahedra-into-Tetrahedra.pdf?origin=publication_detail - // It's a bit hard to look at but it works - // Uses vertex numberings to ensure consistent diagonals between faces, and keeps tet counts to 5 or 6 per hex - size_t tetCount = 0; - // Get number of tets first - for (size_t iC = 0; iC < nCells(); iC++) { - switch (cellType(iC)) { - case VolumeCellType::HEX: { - std::array sortedNumbering; - std::iota(sortedNumbering.begin(), sortedNumbering.end(), 0); - std::sort(sortedNumbering.begin(), sortedNumbering.end(), - [this, iC](size_t a, size_t b) -> bool { return cells[iC][a] < cells[iC][b]; }); - std::array rotatedNumbering; - std::copy(rotationMap[sortedNumbering[0]].begin(), rotationMap[sortedNumbering[0]].end(), - rotatedNumbering.begin()); - size_t diagCount = 0; - auto checkDiagonal = [this, rotatedNumbering, iC](size_t a1, size_t a2, size_t b1, size_t b2) { - return (cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b1]] && - cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b2]]) || - (cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b1]] && - cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b2]]); - }; - if (checkDiagonal(1, 7, 2, 5)) { - diagCount++; - } - if (checkDiagonal(3, 7, 2, 6)) { - diagCount++; - } - if (checkDiagonal(4, 7, 5, 6)) { - diagCount++; - } - if (diagCount == 0) { - tetCount += 5; - } else { - tetCount += 6; - } - break; - } - case VolumeCellType::TET: - tetCount += 1; - break; - case VolumeCellType::PRISM: - tetCount += 3; - break; - case VolumeCellType::PYRAMID: - tetCount += 2; - break; - } - } - - // Each hex can make up to 6 tets - tets.resize(tetCount); - size_t tetIdx = 0; + tets.clear(); auto addTets = [&](const detail::Tets& newTets) { - for (const auto& tet : newTets) { - for (size_t i = 0; i < 4; i++) { - tets[tetIdx][i] = tet[i]; - } - tetIdx++; - } + for (const auto& tet : newTets) tets.push_back(tet); }; for (size_t iC = 0; iC < nCells(); iC++) { switch (cellType(iC)) { - case VolumeCellType::HEX: { - std::array sortedNumbering; - std::iota(sortedNumbering.begin(), sortedNumbering.end(), 0); - std::sort(sortedNumbering.begin(), sortedNumbering.end(), - [this, iC](size_t a, size_t b) -> bool { return cells[iC][a] < cells[iC][b]; }); - std::array rotatedNumbering; - std::copy(rotationMap[sortedNumbering[0]].begin(), rotationMap[sortedNumbering[0]].end(), - rotatedNumbering.begin()); - size_t n = 0; - size_t diagCount = 0; - // Diagonal exists on the pair of vertices which contain the minimum vertex number - auto checkDiagonal = [this, rotatedNumbering, iC](size_t a1, size_t a2, size_t b1, size_t b2) { - return (cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b1]] && - cells[iC][rotatedNumbering[a1]] < cells[iC][rotatedNumbering[b2]]) || - (cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b1]] && - cells[iC][rotatedNumbering[a2]] < cells[iC][rotatedNumbering[b2]]); - }; - // Minimum vertex will always have 3 diagonals, check other three faces - if (checkDiagonal(1, 7, 2, 5)) { - n += 4; - diagCount++; - } - if (checkDiagonal(3, 7, 2, 6)) { - n += 2; - diagCount++; - } - if (checkDiagonal(4, 7, 5, 6)) { - n += 1; - diagCount++; - } - // Rotate by 120 or 240 degrees depending on diagonal positions - if (n == 1 || n == 6) { - size_t temp = rotatedNumbering[1]; - rotatedNumbering[1] = rotatedNumbering[4]; - rotatedNumbering[4] = rotatedNumbering[3]; - rotatedNumbering[3] = temp; - temp = rotatedNumbering[5]; - rotatedNumbering[5] = rotatedNumbering[6]; - rotatedNumbering[6] = rotatedNumbering[2]; - rotatedNumbering[2] = temp; - } else if (n == 2 || n == 5) { - size_t temp = rotatedNumbering[1]; - rotatedNumbering[1] = rotatedNumbering[3]; - rotatedNumbering[3] = rotatedNumbering[4]; - rotatedNumbering[4] = temp; - temp = rotatedNumbering[5]; - rotatedNumbering[5] = rotatedNumbering[2]; - rotatedNumbering[2] = rotatedNumbering[6]; - rotatedNumbering[6] = temp; - } - - // Map final tets according to diagonalMap and the number of diagonals not incident to V_0 - std::array, 6> tetMap = diagonalMap[diagCount]; - for (size_t k = 0; k < (diagCount == 0 ? 5 : 6); k++) { - for (size_t i = 0; i < 4; i++) { - tets[tetIdx][i] = cells[iC][rotatedNumbering[tetMap[k][i]]]; - } - tetIdx++; - } + case VolumeCellType::HEX: + addTets(detail::decomposeHex(cells[iC])); break; - } - case VolumeCellType::TET: { - for (size_t i = 0; i < 4; i++) { - tets[tetIdx][i] = cells[iC][i]; - } - tetIdx++; + case VolumeCellType::TET: + tets.push_back({cells[iC][0], cells[iC][1], cells[iC][2], cells[iC][3]}); break; - } case VolumeCellType::PRISM: addTets(detail::decomposePrism(cells[iC])); break; From 0c06462d6f117c748ce4d7702ca677db169ed428 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 20:34:27 -0700 Subject: [PATCH 08/15] cache isReal labels, treat all cell types consistently Co-authored-by: Copilot --- include/polyscope/volume_mesh.h | 12 +++++ src/volume_mesh.cpp | 93 ++++++++++++++++++++++++++------- 2 files changed, 85 insertions(+), 20 deletions(-) diff --git a/include/polyscope/volume_mesh.h b/include/polyscope/volume_mesh.h index 0bd700bc..a575c52e 100644 --- a/include/polyscope/volume_mesh.h +++ b/include/polyscope/volume_mesh.h @@ -175,6 +175,9 @@ class VolumeMesh : public Structure { void fillSliceGeometryBuffers(render::ShaderProgram& p); static const std::vector>>& cellStencil(VolumeCellType type); static const std::vector>& cellFaces(VolumeCellType type); + // For each (face, triangle, edge-within-triangle), whether that edge is a real polygon boundary edge. + // Indexed as [face][tri][k], paralleling cellStencil(). + static const std::vector>>& cellRealEdgeStencil(VolumeCellType type); // Slice plane listeners std::vector volumeSlicePlaneListeners; void addSlicePlaneListener(polyscope::SlicePlane* sp); @@ -266,6 +269,15 @@ class VolumeMesh : public Structure { static const std::array, 8> rotationMap; static const std::array, 6>, 4> diagonalMap; + // precomputed real-edge flags: [face][tri][k] = is edge k of that triangle a polygon boundary edge? + static std::vector>> realEdgeStencilTet; + static std::vector>> realEdgeStencilHex; + static std::vector>> realEdgeStencilPrism; + static std::vector>> realEdgeStencilPyramid; + + static bool constantDataInitialized; + static void initializeConstantData(); + // clang-format off // === Quantity adders diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index 8927572a..03213354 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -98,6 +98,52 @@ const std::vector>> VolumeMesh::stencilHex = }; // clang-format on +namespace { +// Computes, for each face/tri/edge in a stencil, whether that triangle edge lies on the polygon boundary. +std::vector>> +buildRealEdgeStencil(const std::vector>>& stencil, + const std::vector>& faces) { + std::vector>> result; + result.reserve(stencil.size()); + for (size_t f = 0; f < stencil.size(); f++) { + const auto& poly = faces[f]; + std::set> realEdges; + for (size_t p = 0; p < poly.size(); p++) { + size_t a = poly[p], b = poly[(p + 1) % poly.size()]; + if (a > b) std::swap(a, b); + realEdges.insert({a, b}); + } + std::vector> faceResult; + faceResult.reserve(stencil[f].size()); + for (const auto& tri : stencil[f]) { + std::array edgeReal; + for (size_t k = 0; k < 3; k++) { + size_t a = tri[k], b = tri[(k + 1) % 3]; + if (a > b) std::swap(a, b); + edgeReal[k] = realEdges.count({a, b}) > 0; + } + faceResult.push_back(edgeReal); + } + result.push_back(std::move(faceResult)); + } + return result; +} +} // namespace + +std::vector>> VolumeMesh::realEdgeStencilTet; +std::vector>> VolumeMesh::realEdgeStencilHex; +std::vector>> VolumeMesh::realEdgeStencilPrism; +std::vector>> VolumeMesh::realEdgeStencilPyramid; +bool VolumeMesh::constantDataInitialized = false; + +void VolumeMesh::initializeConstantData() { + realEdgeStencilTet = buildRealEdgeStencil(stencilTet, facesTet); + realEdgeStencilHex = buildRealEdgeStencil(stencilHex, facesHex); + realEdgeStencilPrism = buildRealEdgeStencil(stencilPrism, facesPrism); + realEdgeStencilPyramid = buildRealEdgeStencil(stencilPyramid, facesPyramid); + constantDataInitialized = true; +} + namespace detail { using Tet = std::array; @@ -312,8 +358,14 @@ edgeWidth(uniquePrefix() + "edgeWidth", 0.), // == misc values activeLevelSetQuantity(nullptr) + { // clang-format on + + if (!constantDataInitialized) { + initializeConstantData(); + } + vertexPositions.checkInvalidValues(); cullWholeElements.setPassive(true); @@ -747,11 +799,6 @@ void VolumeMesh::computeConnectivityData() { // To mitigate this issue, we fill the buffer such that all exterior faces come first, then all interior faces, so // that exterior faces always win depth ties. This doesn't totally eliminate the problem, but greatly improves the // most egregious cases. - auto makeEdgeKey = [](uint32_t a, uint32_t b) -> std::pair { - if (a > b) std::swap(a, b); - return {a, b}; - }; - // == Allocate buffers triangleVertexInds.data.clear(); triangleVertexInds.data.resize(3 * nFacesTriangulation()); @@ -775,20 +822,13 @@ void VolumeMesh::computeConnectivityData() { const std::array& cell = cells[iC]; VolumeCellType cellT = cellType(iC); - const auto& faceTris = cellStencil(cellT); // triangulated faces - const auto& faceVerts = cellFaces(cellT); // original polygon faces + const auto& faceTris = cellStencil(cellT); // triangulated faces + const auto& faceRealEdges = cellRealEdgeStencil(cellT); // real-edge flags per face/tri/edge // Loop over all faces of the cell for (size_t f = 0; f < faceTris.size(); f++) { const auto& triList = faceTris[f]; - const auto& poly = faceVerts[f]; - - // Build a set of "real" edges from the polygon - std::set> realEdges; - for (size_t p = 0; p < poly.size(); p++) { - realEdges.insert(makeEdgeKey(cell[poly[p]], cell[poly[(p + 1) % poly.size()]])); - } // Loop over triangles in the triangulation of this face for (size_t j = 0; j < triList.size(); j++) { @@ -817,13 +857,11 @@ void VolumeMesh::computeConnectivityData() { baryCoord.data[3 * iData + 2] = glm::vec3{0., 0., 1.}; // Mark edges as real or not - glm::vec3 realFlag(0.0f); - for (size_t k = 0; k < 3; k++) { - bool isReal = realEdges.count(makeEdgeKey(cell[tri[k]], cell[tri[(k + 1) % 3]])) > 0; - realFlag[k] = isReal ? 1.0f : 0.0f; - edgeIsReal.data[3 * iData + k] = realFlag; + for (int k = 0; k < 3; k++) { + for (int c = 0; c < 3; c++) { + edgeIsReal.data[3 * iData + k][c] = faceRealEdges[f][j][c] ? 1.0f : 0.0f; + } } - for (int k = 0; k < 3; k++) edgeIsReal.data[3 * iData + k] = realFlag; } // Face type: 1 for interior, 0 for exterior @@ -872,6 +910,21 @@ const std::vector>& VolumeMesh::cellFaces(VolumeCellType typ return facesTet; } +const std::vector>>& VolumeMesh::cellRealEdgeStencil(VolumeCellType type) { + switch (type) { + case VolumeCellType::TET: + return realEdgeStencilTet; + case VolumeCellType::HEX: + return realEdgeStencilHex; + case VolumeCellType::PRISM: + return realEdgeStencilPrism; + case VolumeCellType::PYRAMID: + return realEdgeStencilPyramid; + } + // unreachable + return realEdgeStencilTet; +} + void VolumeMesh::computeFaceNormals() { From 6064df3111d88f602eadcdbac50bb812e6b9c657 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 20:54:10 -0700 Subject: [PATCH 09/15] make sure the slicing case is tested --- test/include/polyscope_test.h | 28 +++++++++++++++ test/src/volume_mesh_test.cpp | 66 ++++++++++++++++++----------------- 2 files changed, 62 insertions(+), 32 deletions(-) diff --git a/test/include/polyscope_test.h b/test/include/polyscope_test.h index aa6548f9..a84c2ff9 100644 --- a/test/include/polyscope_test.h +++ b/test/include/polyscope_test.h @@ -166,3 +166,31 @@ inline std::tuple, std::vector>> getVo return std::make_tuple(combined_verts, combined_cells); }; + +inline std::tuple, std::vector>> getHexWedgePyramidTetMeshData() { + // clang-format off + std::vector vertices = { + {0, 0, 0}, // V0 (base hex) + {1, 0, 0}, // V1 + {1, 1, 0}, // V2 + {0, 1, 0}, // V3 + {0, 0, 1}, // V4 + {1, 0, 1}, // V5 + {1, 1, 1}, // V6 + {0, 1, 1}, // V7 + {0.0, 0.5, 1.5}, // V8 (top prism) + {1.0, 0.5, 1.5}, // V9 + {1.5, 0.5, 0.0}, // V10 (side prism) + {1.5, 0.5, 1.0}, // V11 + {0.5, 0.5, -0.5}, // V12 (pyramid apex) + }; + std::vector> cells = { + {0, 1, 2, 3, 4, 5, 6, 7}, // hex + {4, 7, 8, 5, 6, 9, -1, -1}, // top prism + {1, 10, 2, 5, 11, 6, -1, -1},// side prism + {0, 3, 2, 1, 12, -1, -1, -1},// pyramid + {5, 11, 6, 9, -1, -1, -1, -1},// tet + }; + // clang-format on + return std::make_tuple(vertices, cells); +}; diff --git a/test/src/volume_mesh_test.cpp b/test/src/volume_mesh_test.cpp index 3bbe6a1b..a99734ad 100644 --- a/test/src/volume_mesh_test.cpp +++ b/test/src/volume_mesh_test.cpp @@ -65,43 +65,45 @@ TEST_F(PolyscopeTest, ShowVolumeMesh) { TEST_F(PolyscopeTest, ShowVolumeMeshHexWedgePyramidTet) { - // clang-format off - std::vector vertices = {// Base hex vertices - {0, 0, 0}, // V0 - {1, 0, 0}, // V1 - {1, 1, 0}, // V2 - {0, 1, 0}, // V3 - {0, 0, 1}, // V4 - {1, 0, 1}, // V5 - {1, 1, 1}, // V6 - {0, 1, 1}, // V7 - // Top Prism Vertices - {0.0, 0.5, 1.5}, // V8 - {1.0, 0.5, 1.5}, // V9 - // Side Prism Vertices - {1.5, 0.5, 0.0}, // V10 - {1.5, 0.5, 1.0}, // V11 - // Bottom Pyramid Vertex - {0.5, 0.5, -0.5} // V12 - }; - // clang-format on - std::vector> cells = { - // Base Hex cell - {0, 1, 2, 3, 4, 5, 6, 7}, - // Top Prism cell - {4, 7, 8, 5, 6, 9, -1, -1}, - // Side Prism cell - {1, 10, 2, 5, 11, 6, -1, -1}, - // Bottom Pyramid cell - {0, 3, 2, 1, 12, -1, -1, -1}, - // Tet Connecting side and top prisms - {5, 11, 6, 9, -1, -1, -1, -1}, - }; + std::vector vertices; + std::vector> cells; + std::tie(vertices, cells) = getHexWedgePyramidTetMeshData(); polyscope::registerVolumeMesh("hex prism pyramid tet", vertices, cells); polyscope::show(3); polyscope::removeAllStructures(); } +TEST_F(PolyscopeTest, VolumeMeshInspectHexWedgePyramidTet) { + std::vector vertices; + std::vector> cells; + std::tie(vertices, cells) = getHexWedgePyramidTetMeshData(); + polyscope::VolumeMesh* psVol = polyscope::registerVolumeMesh("hex prism pyramid tet", vertices, cells); + + polyscope::SlicePlane* p = polyscope::addSlicePlane(); + p->setVolumeMeshToInspect("hex prism pyramid tet"); + polyscope::show(3); + + // with a scalar quantity + std::vector vals(vertices.size(), 0.44); + auto q1 = psVol->addVertexScalarQuantity("vals", vals); + q1->setEnabled(true); + polyscope::show(3); + + // with a categorical quantity + auto q1Cat = psVol->addVertexScalarQuantity("vals cat", vals, polyscope::DataType::CATEGORICAL); + q1Cat->setEnabled(true); + polyscope::show(3); + + // try orthographic rendering + polyscope::view::setProjectionMode(polyscope::ProjectionMode::Orthographic); + polyscope::show(3); + polyscope::view::setProjectionMode(polyscope::ProjectionMode::Perspective); + polyscope::show(3); + + polyscope::removeAllStructures(); + polyscope::removeLastSceneSlicePlane(); +} + TEST_F(PolyscopeTest, VolumeMeshUpdatePositions) { std::vector verts; std::vector> cells; From 3ceacc66f3490cecffe40541109926215562ca69 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 20:54:27 -0700 Subject: [PATCH 10/15] comment fixes, use anon namespace --- src/volume_mesh.cpp | 57 ++++++++++++++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 19 deletions(-) diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index 03213354..88dbdf1c 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -144,27 +144,32 @@ void VolumeMesh::initializeConstantData() { constantDataInitialized = true; } +// some internal helpers for connectivity management +namespace { -namespace detail { using Tet = std::array; using Prism = std::array; using Pyramid = std::array; using Cell = std::array; using Tets = std::vector; -Prism rotatePrism(const Prism& p, int rot) { - static constexpr int bottomRot[3][3] = { - {0, 1, 2}, // identity - {1, 2, 0}, // rot1 - {2, 0, 1} // rot2 - }; - static constexpr int topRot[3][3] = {{3, 4, 5}, {4, 5, 3}, {5, 3, 4}}; - Prism r; - for (int i = 0; i < 3; ++i) r[i] = p[bottomRot[rot][i]]; - for (int i = 0; i < 3; ++i) r[i + 3] = p[topRot[rot][i]]; - return r; -} Tets decomposePrism(const Cell& cell) { + // A helper to decompose a prism cell into tets. + // + // The challenge is to make sure we tet-decompse adjacent cells in a consistent way, so there are not little gaps + // between mismatched tetrahedralizations for non-planar cell faces. + // + // WARNING: I'm not sure if we handle all mixed-element cases correctly + + auto rotatePrism = [](const Prism& p, int rot) -> Prism { + static constexpr int bottomRot[3][3] = {{0, 1, 2}, {1, 2, 0}, {2, 0, 1}}; + static constexpr int topRot[3][3] = {{3, 4, 5}, {4, 5, 3}, {5, 3, 4}}; + Prism r; + for (int i = 0; i < 3; ++i) r[i] = p[bottomRot[rot][i]]; + for (int i = 0; i < 3; ++i) r[i + 3] = p[topRot[rot][i]]; + return r; + }; + Prism p = {cell[0], cell[1], cell[2], cell[3], cell[4], cell[5]}; int minIdx = std::min_element(p.begin(), p.end()) - p.begin(); int rot = 0; @@ -194,6 +199,13 @@ Tets decomposePrism(const Cell& cell) { } } Tets decomposePyramid(const Cell& cell) { + // A helper to decompose a pyramid cell into tets. + // + // The challenge is to make sure we tet-decompse adjacent cells in a consistent way, so there are not little gaps + // between mismatched tetrahedralizations for non-planar cell faces. + // + // WARNING: I'm not sure if we handle all mixed-element cases correctly + Pyramid p = {cell[0], cell[1], cell[2], cell[3], cell[4]}; if (std::min(p[0], p[2]) < std::min(p[1], p[3])) { // Split along diagonal 0-2 @@ -204,8 +216,15 @@ Tets decomposePyramid(const Cell& cell) { } } Tets decomposeHex(const Cell& cell) { + // A helper to decompose a hex cell into tets. + // // The challenge is to make sure we tet-decompse adjacent cells in a consistent way, so there are not little gaps - // between mismatched tetrahedralizations, for non-planar faces. This approach is adapted from + // between mismatched tetrahedralizations for non-planar cell faces. This implements a strategy that gets it right, at + // least for pure hex meshes. + // + // WARNING: I'm not sure if we handle all mixed-element cases correctly + // + // This approach is adapted from // https://www.researchgate.net/profile/Julien-Dompierre/publication/221561839_How_to_Subdivide_Pyramids_Prisms_and_Hexahedra_into_Tetrahedra/links/0912f509c0b7294059000000/How-to-Subdivide-Pyramids-Prisms-and-Hexahedra-into-Tetrahedra.pdf?origin=publication_detail // It's a bit hard to look at but it works. // The resulting tet counts are either 5 or 6 per hex. Like the rest of this class, this function only really works @@ -318,7 +337,7 @@ Tets decomposeHex(const Cell& cell) { } return result; } -} // namespace detail +} // namespace VolumeMesh::VolumeMesh(std::string name, const std::vector& vertexPositions_, const std::vector>& cellIndices_) @@ -454,22 +473,22 @@ void VolumeMesh::computeCounts() { void VolumeMesh::computeTets() { tets.clear(); - auto addTets = [&](const detail::Tets& newTets) { + auto addTets = [&](const Tets& newTets) { for (const auto& tet : newTets) tets.push_back(tet); }; for (size_t iC = 0; iC < nCells(); iC++) { switch (cellType(iC)) { case VolumeCellType::HEX: - addTets(detail::decomposeHex(cells[iC])); + addTets(decomposeHex(cells[iC])); break; case VolumeCellType::TET: tets.push_back({cells[iC][0], cells[iC][1], cells[iC][2], cells[iC][3]}); break; case VolumeCellType::PRISM: - addTets(detail::decomposePrism(cells[iC])); + addTets(decomposePrism(cells[iC])); break; case VolumeCellType::PYRAMID: - addTets(detail::decomposePyramid(cells[iC])); + addTets(decomposePyramid(cells[iC])); break; } } From 91f14169d6cecca22c6596f84dce773a7f15a2d7 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 20:58:48 -0700 Subject: [PATCH 11/15] add some docs comments --- include/polyscope/volume_mesh.h | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/include/polyscope/volume_mesh.h b/include/polyscope/volume_mesh.h index a575c52e..ffb8eddc 100644 --- a/include/polyscope/volume_mesh.h +++ b/include/polyscope/volume_mesh.h @@ -297,12 +297,22 @@ class VolumeMesh : public Structure { }; // Register functions + +// Register a volume mesh with all tetrahedral cells. tetIndices has shape (nTets, 4). template VolumeMesh* registerTetMesh(std::string name, const V& vertexPositions, const C& tetIndices); + +// Register a volume mesh with all hexahedral cells. hexIndices has shape (nHexes, 8). template VolumeMesh* registerHexMesh(std::string name, const V& vertexPositions, const C& hexIndices); + +// Register a volume mesh with general mixed cell types (tet, hex, prism, pyramid). cellIndices has shape (nCells, 8); +// cells with fewer than 8 vertices are right-padded with negative indices to indicate unused slots. template -VolumeMesh* registerVolumeMesh(std::string name, const V& vertexPositions, const C& hexIndices); +VolumeMesh* registerVolumeMesh(std::string name, const V& vertexPositions, const C& cellIndices); + +// Register a volume mesh from separate tet and hex index arrays. +// Cells are ordered with all tets first, then all hexes. Note that pyrams and prisms are also supported, but not from this function. template VolumeMesh* registerTetHexMesh(std::string name, const V& vertexPositions, const Ct& tetIndices, const Ch& hexIndices); From c81b4b71fcd6536ce48d98d26e200f142abd2f61 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 22:40:31 -0700 Subject: [PATCH 12/15] change prism indexing order --- src/volume_mesh.cpp | 20 ++++++++++---------- test/include/polyscope_test.h | 4 ++-- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index 88dbdf1c..a42cdf4a 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -63,20 +63,20 @@ const std::vector>> VolumeMesh::stencilHex = const std::vector>> VolumeMesh::stencilPrism = { - {{0,1,2}}, - {{0,4,1}, {0,3,4}}, - {{1,5,4}, {1,2,5}}, - {{3,0,5}, {0,2,5}}, - {{3,5,4}}, + {{0,2,1}}, + {{0,5,2}, {0,3,5}}, + {{2,4,5}, {2,1,4}}, + {{3,0,4}, {0,1,4}}, + {{3,4,5}}, }; const std::vector> VolumeMesh::facesPrism = { - {0,1,2}, - {0,3,4,1}, - {1,4,5,2}, - {0,2,5,3}, - {3,5,4} + {0,2,1}, + {0,3,5,2}, + {2,5,4,1}, + {0,1,4,3}, + {3,4,5} }; const std::vector>> VolumeMesh::stencilPyramid = diff --git a/test/include/polyscope_test.h b/test/include/polyscope_test.h index a84c2ff9..e8749d6a 100644 --- a/test/include/polyscope_test.h +++ b/test/include/polyscope_test.h @@ -186,8 +186,8 @@ inline std::tuple, std::vector>> getHe }; std::vector> cells = { {0, 1, 2, 3, 4, 5, 6, 7}, // hex - {4, 7, 8, 5, 6, 9, -1, -1}, // top prism - {1, 10, 2, 5, 11, 6, -1, -1},// side prism + {4, 8, 7, 5, 9, 6, -1, -1}, // top prism + {1, 2, 10, 5, 6, 11, -1, -1}, // side prism {0, 3, 2, 1, 12, -1, -1, -1},// pyramid {5, 11, 6, 9, -1, -1, -1, -1},// tet }; From 94d436e9c7e685504fdf4de584be8d7ced8ea097 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 22:46:59 -0700 Subject: [PATCH 13/15] add slice plane inspect to volume mesh menu --- src/volume_mesh.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index a42cdf4a..fe128d36 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -1151,6 +1151,23 @@ void VolumeMesh::buildCustomOptionsUI() { material.manuallyChanged(); setMaterial(material.get()); // trigger the other updates that happen on set() } + + if (ImGui::BeginMenu("Inspect with slice plane")) { + if (state::slicePlanes.empty()) { + if (ImGui::Button("Add slice plane")) { + SlicePlane* sp = addSlicePlane(); + sp->setVolumeMeshToInspect(getName()); + } + } else { + for (auto& sp : state::slicePlanes) { + bool isInspecting = sp->getVolumeMeshToInspect() == getName(); + if (ImGui::MenuItem(sp->name.c_str(), NULL, isInspecting)) { + sp->setVolumeMeshToInspect(isInspecting ? "" : getName()); + } + } + } + ImGui::EndMenu(); + } } From 06819ad82345d934742bb56b3c40d064dbba8263 Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 22:47:29 -0700 Subject: [PATCH 14/15] minor --- src/slice_plane.cpp | 1 + src/volume_mesh.cpp | 2 +- src/volume_mesh_scalar_quantity.cpp | 2 ++ 3 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/slice_plane.cpp b/src/slice_plane.cpp index 204b3756..f84f4185 100644 --- a/src/slice_plane.cpp +++ b/src/slice_plane.cpp @@ -174,6 +174,7 @@ void SlicePlane::setSliceGeomUniforms(render::ShaderProgram& p) { void SlicePlane::setVolumeMeshToInspect(std::string meshname) { + requestRedraw(); VolumeMesh* oldMeshToInspect = inspectedMeshName == "" ? nullptr : polyscope::getVolumeMesh(inspectedMeshName); if (oldMeshToInspect != nullptr) { oldMeshToInspect->removeSlicePlaneListener(this); diff --git a/src/volume_mesh.cpp b/src/volume_mesh.cpp index fe128d36..dcd6bc1d 100644 --- a/src/volume_mesh.cpp +++ b/src/volume_mesh.cpp @@ -61,7 +61,7 @@ const std::vector>> VolumeMesh::stencilHex = {7,4,5,6} }; - const std::vector>> VolumeMesh::stencilPrism = + const std::vector>> VolumeMesh::stencilPrism = { {{0,2,1}}, {{0,5,2}, {0,3,5}}, diff --git a/src/volume_mesh_scalar_quantity.cpp b/src/volume_mesh_scalar_quantity.cpp index e3bf2a9c..67085a30 100644 --- a/src/volume_mesh_scalar_quantity.cpp +++ b/src/volume_mesh_scalar_quantity.cpp @@ -152,6 +152,7 @@ void VolumeMeshVertexScalarQuantity::setEnabledLevelSet(bool v) { isDrawingLevelSet = false; parent.setLevelSetQuantity(nullptr); } + requestRedraw(); } void VolumeMeshVertexScalarQuantity::drawSlice(polyscope::SlicePlane* sp) { @@ -198,6 +199,7 @@ void VolumeMeshVertexScalarQuantity::setLevelSetVisibleQuantity(std::string name render::engine->setMaterial(*levelSetProgram, parent.getMaterial()); fillLevelSetData(*levelSetProgram); setLevelSetUniforms(*levelSetProgram); + requestRedraw(); showQuantity = q; } From 395b17a0967ce1fa98a8ff8f15920549000f6cbb Mon Sep 17 00:00:00 2001 From: Nicholas Sharp Date: Sat, 2 May 2026 23:25:58 -0700 Subject: [PATCH 15/15] naming --- include/polyscope/volume_mesh.ipp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/polyscope/volume_mesh.ipp b/include/polyscope/volume_mesh.ipp index d06f3640..3bff6e18 100644 --- a/include/polyscope/volume_mesh.ipp +++ b/include/polyscope/volume_mesh.ipp @@ -43,11 +43,11 @@ VolumeMesh* registerHexMesh(std::string name, const V& vertexPositions, const F& } template -VolumeMesh* registerVolumeMesh(std::string name, const V& vertexPositions, const F& faceIndices) { +VolumeMesh* registerVolumeMesh(std::string name, const V& vertexPositions, const F& cellIndices) { checkInitialized(); VolumeMesh* s = new VolumeMesh(name, standardizeVectorArray(vertexPositions), - standardizeVectorArray, 8>(faceIndices)); + standardizeVectorArray, 8>(cellIndices)); bool success = registerStructure(s); if (!success) {