From 96fb3020540553c0897a056c2a5b072b0c87d629 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 14:19:29 +0000 Subject: [PATCH 1/9] Initial plan From bd6b327ea163b2d8a5ca863b8faaf3dc82d0501c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 14:26:24 +0000 Subject: [PATCH 2/9] Add mesh patch splitting functionality and demo Co-authored-by: Ruochun <24469442+Ruochun@users.noreply.github.com> --- src/DEM/BdrsAndObjs.h | 27 +++++++ src/DEM/MeshUtils.cpp | 144 +++++++++++++++++++++++++++++++++ src/demo/CMakeLists.txt | 1 + src/demo/DEMdemo_MeshPatch.cpp | 116 ++++++++++++++++++++++++++ 4 files changed, 288 insertions(+) create mode 100644 src/demo/DEMdemo_MeshPatch.cpp diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index c69e83dd..7e356b87 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -523,6 +523,33 @@ class DEMMeshConnected : public DEMInitializer { Scale(make_float3(s[0], s[1], s[2])); } + //////////////////////////////////////////////////////// + // Mesh patch information for convex patch splitting + //////////////////////////////////////////////////////// + // Patch ID for each triangle facet (-1 means unassigned) + std::vector m_patch_ids; + // Number of patches in this mesh + size_t num_patches = 0; + // Whether patch information has been computed + bool patches_computed = false; + + /// @brief Split the mesh into convex patches based on angle threshold. + /// @param angle_threshold_deg Maximum angle (in degrees) between adjacent face normals to be in same patch. + /// @return Number of patches created. + size_t SplitIntoConvexPatches(float angle_threshold_deg = 30.0f); + + /// @brief Get the patch ID for each triangle. + /// @return Vector of patch IDs (one per triangle). Returns empty vector if patches not computed. + const std::vector& GetPatchIDs() const { return m_patch_ids; } + + /// @brief Get the number of patches in the mesh. + /// @return Number of patches, or 0 if patches not computed. + size_t GetNumPatches() const { return num_patches; } + + /// @brief Check if patch information has been computed. + /// @return True if patches have been computed, false otherwise. + bool ArePatchesComputed() const { return patches_computed; } + //////////////////////////////////////////////////////// // Some geo wildcard-related stuff //////////////////////////////////////////////////////// diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index f0efa3cd..0c586893 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -183,4 +183,148 @@ void DEMMeshConnected::WriteWavefront(const std::string& filename, std::vector 1e-10f) { + normal.x /= length; + normal.y /= length; + normal.z /= length; + } + + return normal; +} + +// Helper function to compute angle between two normals (in degrees) +static float computeAngleBetweenNormals(const float3& n1, const float3& n2) { + // Compute dot product + float dot_product = n1.x * n2.x + n1.y * n2.y + n1.z * n2.z; + + // Clamp to [-1, 1] to avoid numerical issues + dot_product = std::max(-1.0f, std::min(1.0f, dot_product)); + + // Compute angle in radians and convert to degrees + float angle_rad = std::acos(dot_product); + return angle_rad * 180.0f / 3.14159265358979323846f; +} + +// Helper to build adjacency map for triangles (shared edges) +static std::vector> buildAdjacencyMap(const std::vector& face_v_indices) { + size_t num_faces = face_v_indices.size(); + std::vector> adjacency(num_faces); + + // Map from edge (as pair of vertex indices) to faces that share it + std::map, std::vector> edge_to_faces; + + for (size_t i = 0; i < num_faces; ++i) { + const int3& face = face_v_indices[i]; + + // Three edges of the triangle (store with smaller index first for consistency) + std::pair edges[3] = { + {std::min(face.x, face.y), std::max(face.x, face.y)}, + {std::min(face.y, face.z), std::max(face.y, face.z)}, + {std::min(face.z, face.x), std::max(face.z, face.x)} + }; + + for (int e = 0; e < 3; ++e) { + edge_to_faces[edges[e]].push_back(i); + } + } + + // Build adjacency list + for (const auto& entry : edge_to_faces) { + const std::vector& faces = entry.second; + // If two faces share an edge, they are adjacent + if (faces.size() == 2) { + adjacency[faces[0]].push_back(faces[1]); + adjacency[faces[1]].push_back(faces[0]); + } + } + + return adjacency; +} + +size_t DEMMeshConnected::SplitIntoConvexPatches(float angle_threshold_deg) { + if (nTri == 0) { + patches_computed = false; + num_patches = 0; + return 0; + } + + // Initialize patch IDs (all -1 means unassigned) + m_patch_ids.clear(); + m_patch_ids.resize(nTri, -1); + + // Compute face normals for all triangles + std::vector face_normals(nTri); + for (size_t i = 0; i < nTri; ++i) { + const int3& face = m_face_v_indices[i]; + const float3& v0 = m_vertices[face.x]; + const float3& v1 = m_vertices[face.y]; + const float3& v2 = m_vertices[face.z]; + face_normals[i] = computeFaceNormal(v0, v1, v2); + } + + // Build adjacency map (which triangles share edges) + std::vector> adjacency = buildAdjacencyMap(m_face_v_indices); + + // Region growing algorithm to assign patches + int current_patch_id = 0; + std::vector queue; + + for (size_t seed = 0; seed < nTri; ++seed) { + // Skip if already assigned to a patch + if (m_patch_ids[seed] != -1) { + continue; + } + + // Start a new patch from this seed triangle + queue.clear(); + queue.push_back(seed); + m_patch_ids[seed] = current_patch_id; + + // Grow the region + size_t queue_idx = 0; + while (queue_idx < queue.size()) { + size_t current = queue[queue_idx++]; + + // Check all adjacent triangles + for (size_t neighbor : adjacency[current]) { + // Skip if already assigned + if (m_patch_ids[neighbor] != -1) { + continue; + } + + // Check angle between normals + float angle = computeAngleBetweenNormals(face_normals[current], face_normals[neighbor]); + + // If angle is below threshold, add to same patch + if (angle <= angle_threshold_deg) { + m_patch_ids[neighbor] = current_patch_id; + queue.push_back(neighbor); + } + } + } + + // Move to next patch + current_patch_id++; + } + + num_patches = current_patch_id; + patches_computed = true; + + return num_patches; +} + } // end namespace deme diff --git a/src/demo/CMakeLists.txt b/src/demo/CMakeLists.txt index 5463acff..ed6588a6 100644 --- a/src/demo/CMakeLists.txt +++ b/src/demo/CMakeLists.txt @@ -20,6 +20,7 @@ SET(DEMOS DEMdemo_Sieve DEMdemo_SingleSphereCollide DEMdemo_MeshCollide + DEMdemo_MeshPatch DEMdemo_BoxesFalling DEMdemo_TestPack DEMdemo_RotatingDrum diff --git a/src/demo/DEMdemo_MeshPatch.cpp b/src/demo/DEMdemo_MeshPatch.cpp new file mode 100644 index 00000000..58dfe6d4 --- /dev/null +++ b/src/demo/DEMdemo_MeshPatch.cpp @@ -0,0 +1,116 @@ +// Copyright (c) 2021, SBEL GPU Development Team +// Copyright (c) 2021, University of Wisconsin - Madison +// +// SPDX-License-Identifier: BSD-3-Clause + +// ============================================================================= +// A demo that tests mesh patch splitting functionality. +// This demo loads a mesh and splits it into convex patches based on angle +// thresholds, demonstrating the mesh patch splitting utility. +// ============================================================================= + +#include +#include +#include +#include + +#include +#include +#include +#include + +using namespace deme; +using namespace std::filesystem; + +int main() { + std::cout << "========================================" << std::endl; + std::cout << "DEM Mesh Patch Splitting Demo" << std::endl; + std::cout << "========================================" << std::endl; + + // Test with a simple cube mesh + std::cout << "\n--- Test 1: Cube Mesh ---" << std::endl; + auto cube_mesh = std::make_shared(); + bool loaded = cube_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/cube.obj").string()); + + if (loaded) { + std::cout << "Loaded cube mesh successfully" << std::endl; + std::cout << "Number of triangles: " << cube_mesh->GetNumTriangles() << std::endl; + std::cout << "Number of vertices: " << cube_mesh->GetNumNodes() << std::endl; + + // Test different angle thresholds + float thresholds[] = {10.0f, 30.0f, 45.0f, 90.0f}; + + for (float threshold : thresholds) { + size_t num_patches = cube_mesh->SplitIntoConvexPatches(threshold); + std::cout << "\nAngle threshold: " << std::fixed << std::setprecision(1) << threshold << " degrees" << std::endl; + std::cout << "Number of patches: " << num_patches << std::endl; + + // Show patch distribution + if (cube_mesh->ArePatchesComputed()) { + const auto& patch_ids = cube_mesh->GetPatchIDs(); + + // Count triangles per patch + std::map patch_counts; + for (int patch_id : patch_ids) { + patch_counts[patch_id]++; + } + + std::cout << "Patch distribution:" << std::endl; + for (const auto& entry : patch_counts) { + std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; + } + } + } + } else { + std::cout << "Failed to load cube mesh" << std::endl; + } + + // Test with sphere mesh if available + std::cout << "\n--- Test 2: Sphere Mesh ---" << std::endl; + auto sphere_mesh = std::make_shared(); + loaded = sphere_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/sphere.obj").string()); + + if (loaded) { + std::cout << "Loaded sphere mesh successfully" << std::endl; + std::cout << "Number of triangles: " << sphere_mesh->GetNumTriangles() << std::endl; + std::cout << "Number of vertices: " << sphere_mesh->GetNumNodes() << std::endl; + + // Test with 30 degree threshold + size_t num_patches = sphere_mesh->SplitIntoConvexPatches(30.0f); + std::cout << "Split into " << num_patches << " patches (threshold: 30 degrees)" << std::endl; + + if (sphere_mesh->ArePatchesComputed()) { + const auto& patch_ids = sphere_mesh->GetPatchIDs(); + + // Count triangles per patch + std::map patch_counts; + for (int patch_id : patch_ids) { + patch_counts[patch_id]++; + } + + std::cout << "Number of patches with different sizes:" << std::endl; + std::map size_distribution; + for (const auto& entry : patch_counts) { + size_distribution[entry.second]++; + } + for (const auto& entry : size_distribution) { + std::cout << " " << entry.second << " patches with " << entry.first << " triangles each" << std::endl; + } + } + } else { + std::cout << "Sphere mesh not available, skipping" << std::endl; + } + + // Test edge case: empty mesh + std::cout << "\n--- Test 3: Empty Mesh ---" << std::endl; + auto empty_mesh = std::make_shared(); + size_t num_patches = empty_mesh->SplitIntoConvexPatches(30.0f); + std::cout << "Empty mesh patches: " << num_patches << " (expected: 0)" << std::endl; + std::cout << "Patches computed: " << (empty_mesh->ArePatchesComputed() ? "yes" : "no") << " (expected: no)" << std::endl; + + std::cout << "\n========================================" << std::endl; + std::cout << "Demo completed successfully!" << std::endl; + std::cout << "========================================" << std::endl; + + return 0; +} From cd75d6904fa4f6a19bf1c150bd29ad5bd38ab42c Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 14:29:16 +0000 Subject: [PATCH 3/9] Add missing headers and update Clear() method for patches Co-authored-by: Ruochun <24469442+Ruochun@users.noreply.github.com> --- src/DEM/BdrsAndObjs.h | 3 +++ src/DEM/MeshUtils.cpp | 3 +++ 2 files changed, 6 insertions(+) diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index 7e356b87..e56d5224 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -363,6 +363,9 @@ class DEMMeshConnected : public DEMInitializer { this->m_face_n_indices.clear(); this->m_face_uv_indices.clear(); this->m_face_col_indices.clear(); + this->m_patch_ids.clear(); + this->num_patches = 0; + this->patches_computed = false; this->owner = NULL_BODYID; } diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index 0c586893..65ec9e14 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -32,6 +32,9 @@ #include #include #include +#include +#include +#include #include "../kernel/DEMHelperKernels.cuh" #include "BdrsAndObjs.h" From 85405172e77589dbf43b47398c986ac6d4f690d1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 14:32:18 +0000 Subject: [PATCH 4/9] Enhance documentation for mesh patch splitting Co-authored-by: Ruochun <24469442+Ruochun@users.noreply.github.com> --- src/DEM/BdrsAndObjs.h | 5 +++++ src/DEM/MeshUtils.cpp | 3 +++ 2 files changed, 8 insertions(+) diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index e56d5224..9a704ce9 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -537,7 +537,12 @@ class DEMMeshConnected : public DEMInitializer { bool patches_computed = false; /// @brief Split the mesh into convex patches based on angle threshold. + /// @details Uses a region-growing algorithm to group adjacent triangles whose face normals differ by less than + /// the specified angle threshold. Each patch represents a locally convex region of the mesh. Patches are + /// non-overlapping and cover the entire mesh. This is useful for contact force calculations. /// @param angle_threshold_deg Maximum angle (in degrees) between adjacent face normals to be in same patch. + /// Default is 30.0 degrees. Lower values create more patches (stricter convexity), higher values create fewer + /// patches (relaxed convexity). /// @return Number of patches created. size_t SplitIntoConvexPatches(float angle_threshold_deg = 30.0f); diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index 65ec9e14..d5639bef 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -258,6 +258,9 @@ static std::vector> buildAdjacencyMap(const std::vector Date: Sat, 8 Nov 2025 14:34:11 +0000 Subject: [PATCH 5/9] Final implementation complete Co-authored-by: Ruochun <24469442+Ruochun@users.noreply.github.com> --- _codeql_detected_source_root | 1 + 1 file changed, 1 insertion(+) create mode 120000 _codeql_detected_source_root diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root new file mode 120000 index 00000000..945c9b46 --- /dev/null +++ b/_codeql_detected_source_root @@ -0,0 +1 @@ +. \ No newline at end of file From 6383383a1b3e8478e1c3e0be33ff4fce72b836c1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 14:34:45 +0000 Subject: [PATCH 6/9] Remove CodeQL artifact --- _codeql_detected_source_root | 1 - 1 file changed, 1 deletion(-) delete mode 120000 _codeql_detected_source_root diff --git a/_codeql_detected_source_root b/_codeql_detected_source_root deleted file mode 120000 index 945c9b46..00000000 --- a/_codeql_detected_source_root +++ /dev/null @@ -1 +0,0 @@ -. \ No newline at end of file From f1798f90a6ef0b94cacea7c5962cd3b1ee0f3d20 Mon Sep 17 00:00:00 2001 From: Ruochun Date: Sun, 9 Nov 2025 01:10:31 +0800 Subject: [PATCH 7/9] Minor improvement --- src/DEM/MeshUtils.cpp | 67 ++++++++++++++++------------------ src/demo/DEMdemo_MeshPatch.cpp | 42 +++++++++++---------- 2 files changed, 53 insertions(+), 56 deletions(-) diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index d5639bef..5d53799f 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -190,22 +190,19 @@ void DEMMeshConnected::WriteWavefront(const std::string& filename, std::vector 1e-10f) { + if (length > DEME_TINY_FLOAT) { normal.x /= length; normal.y /= length; normal.z /= length; } - + return normal; } @@ -213,38 +210,36 @@ static float3 computeFaceNormal(const float3& v0, const float3& v1, const float3 static float computeAngleBetweenNormals(const float3& n1, const float3& n2) { // Compute dot product float dot_product = n1.x * n2.x + n1.y * n2.y + n1.z * n2.z; - + // Clamp to [-1, 1] to avoid numerical issues dot_product = std::max(-1.0f, std::min(1.0f, dot_product)); - + // Compute angle in radians and convert to degrees float angle_rad = std::acos(dot_product); - return angle_rad * 180.0f / 3.14159265358979323846f; + return angle_rad * 180.0f / deme::PI; } // Helper to build adjacency map for triangles (shared edges) static std::vector> buildAdjacencyMap(const std::vector& face_v_indices) { size_t num_faces = face_v_indices.size(); std::vector> adjacency(num_faces); - + // Map from edge (as pair of vertex indices) to faces that share it std::map, std::vector> edge_to_faces; - + for (size_t i = 0; i < num_faces; ++i) { const int3& face = face_v_indices[i]; - + // Three edges of the triangle (store with smaller index first for consistency) - std::pair edges[3] = { - {std::min(face.x, face.y), std::max(face.x, face.y)}, - {std::min(face.y, face.z), std::max(face.y, face.z)}, - {std::min(face.z, face.x), std::max(face.z, face.x)} - }; - + std::pair edges[3] = {{std::min(face.x, face.y), std::max(face.x, face.y)}, + {std::min(face.y, face.z), std::max(face.y, face.z)}, + {std::min(face.z, face.x), std::max(face.z, face.x)}}; + for (int e = 0; e < 3; ++e) { edge_to_faces[edges[e]].push_back(i); } } - + // Build adjacency list for (const auto& entry : edge_to_faces) { const std::vector& faces = entry.second; @@ -254,7 +249,7 @@ static std::vector> buildAdjacencyMap(const std::vector face_normals(nTri); for (size_t i = 0; i < nTri; ++i) { @@ -281,40 +276,40 @@ size_t DEMMeshConnected::SplitIntoConvexPatches(float angle_threshold_deg) { const float3& v2 = m_vertices[face.z]; face_normals[i] = computeFaceNormal(v0, v1, v2); } - + // Build adjacency map (which triangles share edges) std::vector> adjacency = buildAdjacencyMap(m_face_v_indices); - + // Region growing algorithm to assign patches int current_patch_id = 0; std::vector queue; - + for (size_t seed = 0; seed < nTri; ++seed) { // Skip if already assigned to a patch if (m_patch_ids[seed] != -1) { continue; } - + // Start a new patch from this seed triangle queue.clear(); queue.push_back(seed); m_patch_ids[seed] = current_patch_id; - + // Grow the region size_t queue_idx = 0; while (queue_idx < queue.size()) { size_t current = queue[queue_idx++]; - + // Check all adjacent triangles for (size_t neighbor : adjacency[current]) { // Skip if already assigned if (m_patch_ids[neighbor] != -1) { continue; } - + // Check angle between normals float angle = computeAngleBetweenNormals(face_normals[current], face_normals[neighbor]); - + // If angle is below threshold, add to same patch if (angle <= angle_threshold_deg) { m_patch_ids[neighbor] = current_patch_id; @@ -322,14 +317,14 @@ size_t DEMMeshConnected::SplitIntoConvexPatches(float angle_threshold_deg) { } } } - + // Move to next patch current_patch_id++; } - + num_patches = current_patch_id; patches_computed = true; - + return num_patches; } diff --git a/src/demo/DEMdemo_MeshPatch.cpp b/src/demo/DEMdemo_MeshPatch.cpp index 58dfe6d4..c5bdc5f1 100644 --- a/src/demo/DEMdemo_MeshPatch.cpp +++ b/src/demo/DEMdemo_MeshPatch.cpp @@ -5,7 +5,7 @@ // ============================================================================= // A demo that tests mesh patch splitting functionality. -// This demo loads a mesh and splits it into convex patches based on angle +// This demo loads a mesh and splits it into convex patches based on angle // thresholds, demonstrating the mesh patch splitting utility. // ============================================================================= @@ -26,35 +26,36 @@ int main() { std::cout << "========================================" << std::endl; std::cout << "DEM Mesh Patch Splitting Demo" << std::endl; std::cout << "========================================" << std::endl; - + // Test with a simple cube mesh std::cout << "\n--- Test 1: Cube Mesh ---" << std::endl; auto cube_mesh = std::make_shared(); bool loaded = cube_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/cube.obj").string()); - + if (loaded) { std::cout << "Loaded cube mesh successfully" << std::endl; std::cout << "Number of triangles: " << cube_mesh->GetNumTriangles() << std::endl; std::cout << "Number of vertices: " << cube_mesh->GetNumNodes() << std::endl; - + // Test different angle thresholds - float thresholds[] = {10.0f, 30.0f, 45.0f, 90.0f}; - + float thresholds[] = {10.0f, 45.0f, 90.0f, 300.0f}; + for (float threshold : thresholds) { size_t num_patches = cube_mesh->SplitIntoConvexPatches(threshold); - std::cout << "\nAngle threshold: " << std::fixed << std::setprecision(1) << threshold << " degrees" << std::endl; + std::cout << "\nAngle threshold: " << std::fixed << std::setprecision(1) << threshold << " degrees" + << std::endl; std::cout << "Number of patches: " << num_patches << std::endl; - + // Show patch distribution if (cube_mesh->ArePatchesComputed()) { const auto& patch_ids = cube_mesh->GetPatchIDs(); - + // Count triangles per patch std::map patch_counts; for (int patch_id : patch_ids) { patch_counts[patch_id]++; } - + std::cout << "Patch distribution:" << std::endl; for (const auto& entry : patch_counts) { std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; @@ -64,30 +65,30 @@ int main() { } else { std::cout << "Failed to load cube mesh" << std::endl; } - + // Test with sphere mesh if available std::cout << "\n--- Test 2: Sphere Mesh ---" << std::endl; auto sphere_mesh = std::make_shared(); loaded = sphere_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/sphere.obj").string()); - + if (loaded) { std::cout << "Loaded sphere mesh successfully" << std::endl; std::cout << "Number of triangles: " << sphere_mesh->GetNumTriangles() << std::endl; std::cout << "Number of vertices: " << sphere_mesh->GetNumNodes() << std::endl; - + // Test with 30 degree threshold size_t num_patches = sphere_mesh->SplitIntoConvexPatches(30.0f); std::cout << "Split into " << num_patches << " patches (threshold: 30 degrees)" << std::endl; - + if (sphere_mesh->ArePatchesComputed()) { const auto& patch_ids = sphere_mesh->GetPatchIDs(); - + // Count triangles per patch std::map patch_counts; for (int patch_id : patch_ids) { patch_counts[patch_id]++; } - + std::cout << "Number of patches with different sizes:" << std::endl; std::map size_distribution; for (const auto& entry : patch_counts) { @@ -100,17 +101,18 @@ int main() { } else { std::cout << "Sphere mesh not available, skipping" << std::endl; } - + // Test edge case: empty mesh std::cout << "\n--- Test 3: Empty Mesh ---" << std::endl; auto empty_mesh = std::make_shared(); size_t num_patches = empty_mesh->SplitIntoConvexPatches(30.0f); std::cout << "Empty mesh patches: " << num_patches << " (expected: 0)" << std::endl; - std::cout << "Patches computed: " << (empty_mesh->ArePatchesComputed() ? "yes" : "no") << " (expected: no)" << std::endl; - + std::cout << "Patches computed: " << (empty_mesh->ArePatchesComputed() ? "yes" : "no") << " (expected: no)" + << std::endl; + std::cout << "\n========================================" << std::endl; std::cout << "Demo completed successfully!" << std::endl; std::cout << "========================================" << std::endl; - + return 0; } From 7b7b74092e67b62567caa3035e1d434cb55d8081 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 8 Nov 2025 17:27:52 +0000 Subject: [PATCH 8/9] Add default patch info and manual patch setting capability Co-authored-by: Ruochun <24469442+Ruochun@users.noreply.github.com> --- src/DEM/BdrsAndObjs.h | 29 ++++++++------ src/DEM/MeshUtils.cpp | 46 ++++++++++++++++++++-- src/demo/DEMdemo_MeshPatch.cpp | 72 ++++++++++++++++++++++++---------- 3 files changed, 113 insertions(+), 34 deletions(-) diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index 9a704ce9..07735ad5 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -364,8 +364,8 @@ class DEMMeshConnected : public DEMInitializer { this->m_face_uv_indices.clear(); this->m_face_col_indices.clear(); this->m_patch_ids.clear(); - this->num_patches = 0; - this->patches_computed = false; + this->num_patches = 1; + this->patches_explicitly_set = false; this->owner = NULL_BODYID; } @@ -529,12 +529,12 @@ class DEMMeshConnected : public DEMInitializer { //////////////////////////////////////////////////////// // Mesh patch information for convex patch splitting //////////////////////////////////////////////////////// - // Patch ID for each triangle facet (-1 means unassigned) + // Patch ID for each triangle facet (defaults to 0 for all triangles, assuming convex mesh) std::vector m_patch_ids; // Number of patches in this mesh - size_t num_patches = 0; - // Whether patch information has been computed - bool patches_computed = false; + size_t num_patches = 1; + // Whether patch information has been explicitly set (either computed or manually supplied) + bool patches_explicitly_set = false; /// @brief Split the mesh into convex patches based on angle threshold. /// @details Uses a region-growing algorithm to group adjacent triangles whose face normals differ by less than @@ -546,17 +546,24 @@ class DEMMeshConnected : public DEMInitializer { /// @return Number of patches created. size_t SplitIntoConvexPatches(float angle_threshold_deg = 30.0f); + /// @brief Manually set the patch IDs for each triangle. + /// @details Allows user to manually specify which patch each triangle belongs to. This is useful when + /// the user has pre-computed patch information or wants to define patches based on custom criteria. + /// @param patch_ids Vector of patch IDs, one for each triangle. Must have the same length as the number + /// of triangles in the mesh. Patch IDs should be non-negative integers starting from 0. + void SetPatchIDs(const std::vector& patch_ids); + /// @brief Get the patch ID for each triangle. - /// @return Vector of patch IDs (one per triangle). Returns empty vector if patches not computed. + /// @return Vector of patch IDs (one per triangle). By default, all triangles are in patch 0 (assuming convex mesh). const std::vector& GetPatchIDs() const { return m_patch_ids; } /// @brief Get the number of patches in the mesh. - /// @return Number of patches, or 0 if patches not computed. + /// @return Number of patches. Default is 1 (assuming convex mesh). size_t GetNumPatches() const { return num_patches; } - /// @brief Check if patch information has been computed. - /// @return True if patches have been computed, false otherwise. - bool ArePatchesComputed() const { return patches_computed; } + /// @brief Check if patch information has been explicitly set. + /// @return True if patches have been computed via SplitIntoConvexPatches() or set via SetPatchIDs(), false if using default (single patch). + bool ArePatchesExplicitlySet() const { return patches_explicitly_set; } //////////////////////////////////////////////////////// // Some geo wildcard-related stuff diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index 5d53799f..5b5fb300 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -35,6 +35,7 @@ #include #include #include +#include #include "../kernel/DEMHelperKernels.cuh" #include "BdrsAndObjs.h" @@ -133,6 +134,12 @@ bool DEMMeshConnected::LoadWavefrontMesh(std::string input_file, bool load_norma this->nTri = m_face_v_indices.size(); + // Initialize default patch info: all triangles in patch 0 (assuming convex mesh) + this->m_patch_ids.clear(); + this->m_patch_ids.resize(this->nTri, 0); + this->num_patches = 1; + this->patches_explicitly_set = false; + return true; } @@ -258,8 +265,8 @@ static std::vector> buildAdjacencyMap(const std::vector& patch_ids) { + if (patch_ids.size() != nTri) { + std::stringstream ss; + ss << "SetPatchIDs: Input vector size (" << patch_ids.size() + << ") must match the number of triangles (" << nTri << ") in the mesh." << std::endl; + throw std::runtime_error(ss.str()); + } + + // Validate that all patch IDs are non-negative + for (size_t i = 0; i < patch_ids.size(); ++i) { + if (patch_ids[i] < 0) { + std::stringstream ss; + ss << "SetPatchIDs: Patch ID at index " << i << " is negative (" << patch_ids[i] + << "). All patch IDs must be non-negative integers." << std::endl; + throw std::runtime_error(ss.str()); + } + } + + // Copy the patch IDs + m_patch_ids = patch_ids; + + // Calculate the number of patches (maximum patch ID + 1) + if (!patch_ids.empty()) { + int max_patch_id = *std::max_element(patch_ids.begin(), patch_ids.end()); + num_patches = max_patch_id + 1; + } else { + num_patches = 1; + } + + patches_explicitly_set = true; +} + } // end namespace deme diff --git a/src/demo/DEMdemo_MeshPatch.cpp b/src/demo/DEMdemo_MeshPatch.cpp index c5bdc5f1..13936f0b 100644 --- a/src/demo/DEMdemo_MeshPatch.cpp +++ b/src/demo/DEMdemo_MeshPatch.cpp @@ -28,7 +28,7 @@ int main() { std::cout << "========================================" << std::endl; // Test with a simple cube mesh - std::cout << "\n--- Test 1: Cube Mesh ---" << std::endl; + std::cout << "\n--- Test 1: Cube Mesh with Default Patch Info ---" << std::endl; auto cube_mesh = std::make_shared(); bool loaded = cube_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/cube.obj").string()); @@ -36,8 +36,17 @@ int main() { std::cout << "Loaded cube mesh successfully" << std::endl; std::cout << "Number of triangles: " << cube_mesh->GetNumTriangles() << std::endl; std::cout << "Number of vertices: " << cube_mesh->GetNumNodes() << std::endl; + + // Test default patch info (should be all in patch 0) + std::cout << "\nDefault patch info (assuming convex mesh):" << std::endl; + std::cout << "Number of patches: " << cube_mesh->GetNumPatches() << std::endl; + std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; + const auto& default_patch_ids = cube_mesh->GetPatchIDs(); + std::cout << "All triangles in patch 0: " << (std::all_of(default_patch_ids.begin(), default_patch_ids.end(), + [](int id) { return id == 0; }) ? "yes" : "no") << std::endl; // Test different angle thresholds + std::cout << "\n--- Test 2: Automatic Patch Splitting ---" << std::endl; float thresholds[] = {10.0f, 45.0f, 90.0f, 300.0f}; for (float threshold : thresholds) { @@ -45,29 +54,53 @@ int main() { std::cout << "\nAngle threshold: " << std::fixed << std::setprecision(1) << threshold << " degrees" << std::endl; std::cout << "Number of patches: " << num_patches << std::endl; + std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; // Show patch distribution - if (cube_mesh->ArePatchesComputed()) { - const auto& patch_ids = cube_mesh->GetPatchIDs(); - - // Count triangles per patch - std::map patch_counts; - for (int patch_id : patch_ids) { - patch_counts[patch_id]++; - } - - std::cout << "Patch distribution:" << std::endl; - for (const auto& entry : patch_counts) { - std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; - } + const auto& patch_ids = cube_mesh->GetPatchIDs(); + + // Count triangles per patch + std::map patch_counts; + for (int patch_id : patch_ids) { + patch_counts[patch_id]++; } + + std::cout << "Patch distribution:" << std::endl; + for (const auto& entry : patch_counts) { + std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; + } + } + + // Test manual patch ID setting + std::cout << "\n--- Test 3: Manual Patch ID Setting ---" << std::endl; + size_t num_tris = cube_mesh->GetNumTriangles(); + std::vector manual_patches(num_tris); + // Split triangles into 3 patches based on index + for (size_t i = 0; i < num_tris; ++i) { + manual_patches[i] = i % 3; // Assign patches 0, 1, 2 cyclically + } + + cube_mesh->SetPatchIDs(manual_patches); + std::cout << "Manually set patch IDs (cycling 0, 1, 2)" << std::endl; + std::cout << "Number of patches: " << cube_mesh->GetNumPatches() << std::endl; + std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; + + // Count triangles per patch + const auto& manual_patch_ids = cube_mesh->GetPatchIDs(); + std::map manual_patch_counts; + for (int patch_id : manual_patch_ids) { + manual_patch_counts[patch_id]++; + } + std::cout << "Manual patch distribution:" << std::endl; + for (const auto& entry : manual_patch_counts) { + std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; } } else { std::cout << "Failed to load cube mesh" << std::endl; } // Test with sphere mesh if available - std::cout << "\n--- Test 2: Sphere Mesh ---" << std::endl; + std::cout << "\n--- Test 4: Sphere Mesh ---" << std::endl; auto sphere_mesh = std::make_shared(); loaded = sphere_mesh->LoadWavefrontMesh((GET_DATA_PATH() / "mesh/sphere.obj").string()); @@ -80,7 +113,7 @@ int main() { size_t num_patches = sphere_mesh->SplitIntoConvexPatches(30.0f); std::cout << "Split into " << num_patches << " patches (threshold: 30 degrees)" << std::endl; - if (sphere_mesh->ArePatchesComputed()) { + if (sphere_mesh->ArePatchesExplicitlySet()) { const auto& patch_ids = sphere_mesh->GetPatchIDs(); // Count triangles per patch @@ -103,11 +136,10 @@ int main() { } // Test edge case: empty mesh - std::cout << "\n--- Test 3: Empty Mesh ---" << std::endl; + std::cout << "\n--- Test 5: Empty Mesh ---" << std::endl; auto empty_mesh = std::make_shared(); - size_t num_patches = empty_mesh->SplitIntoConvexPatches(30.0f); - std::cout << "Empty mesh patches: " << num_patches << " (expected: 0)" << std::endl; - std::cout << "Patches computed: " << (empty_mesh->ArePatchesComputed() ? "yes" : "no") << " (expected: no)" + std::cout << "Empty mesh default patches: " << empty_mesh->GetNumPatches() << " (expected: 1)" << std::endl; + std::cout << "Patches explicitly set: " << (empty_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << " (expected: no)" << std::endl; std::cout << "\n========================================" << std::endl; From 7775d88c3737dfc23b1d52db6c1cf9d3b7970d6f Mon Sep 17 00:00:00 2001 From: Ruochun Zhang Date: Mon, 10 Nov 2025 18:23:56 +0800 Subject: [PATCH 9/9] Format and add m_mesh_patch_owner --- src/DEM/API.h | 13 +++++++++---- src/DEM/APIPrivate.cpp | 2 +- src/DEM/BdrsAndObjs.h | 3 ++- src/DEM/MeshUtils.cpp | 6 +++--- src/demo/DEMdemo_MeshPatch.cpp | 22 +++++++++++++--------- 5 files changed, 28 insertions(+), 18 deletions(-) diff --git a/src/DEM/API.h b/src/DEM/API.h index 885354ae..17637b6d 100644 --- a/src/DEM/API.h +++ b/src/DEM/API.h @@ -1749,7 +1749,7 @@ class DEMSolver { // Flattened (analytical) object component definition arrays, potentially jitifiable // These extra analytical entities' owners' ID will be appended to those added thru normal AddClump - std::vector m_anal_owner; + std::vector m_anal_owner; // Material types of these analytical geometries std::vector m_anal_materials; // Initial locations of this obj's components relative to obj's CoM @@ -1769,13 +1769,18 @@ class DEMSolver { // topologically a plane then this param is meaningless, since its normal is determined by its rotation. std::vector m_anal_normals; - // These extra mesh facets' owners' ID will be appended to analytical entities' - std::vector m_mesh_facet_owner; + // These mesh facets' owners' ID, flattened + std::vector m_mesh_facet_owner; // Material types of these mesh facets std::vector m_mesh_facet_materials; - // Material types of these mesh facets + // Three nodes of each triangle, flattened std::vector m_mesh_facets; + // These mesh patches' owners' ID, flattened + std::vector m_mesh_patch_owner; + // Material types of these mesh patches + std::vector m_mesh_patch_materials; + // Clump templates will be flatten and transferred into kernels upon Initialize() std::vector m_template_clump_mass; std::vector m_template_clump_moi; diff --git a/src/DEM/APIPrivate.cpp b/src/DEM/APIPrivate.cpp index 0f93147f..bf5d6a3e 100644 --- a/src/DEM/APIPrivate.cpp +++ b/src/DEM/APIPrivate.cpp @@ -757,7 +757,7 @@ void DEMSolver::preprocessClumps() { void DEMSolver::preprocessTriangleObjs() { nTriMeshes += cached_mesh_objs.size(); - unsigned int thisMeshObj = 0; + bodyID_t thisMeshObj = 0; for (const auto& mesh_obj : cached_mesh_objs) { if (!(mesh_obj->isMaterialSet)) { DEME_ERROR( diff --git a/src/DEM/BdrsAndObjs.h b/src/DEM/BdrsAndObjs.h index 07735ad5..beb81b0a 100644 --- a/src/DEM/BdrsAndObjs.h +++ b/src/DEM/BdrsAndObjs.h @@ -562,7 +562,8 @@ class DEMMeshConnected : public DEMInitializer { size_t GetNumPatches() const { return num_patches; } /// @brief Check if patch information has been explicitly set. - /// @return True if patches have been computed via SplitIntoConvexPatches() or set via SetPatchIDs(), false if using default (single patch). + /// @return True if patches have been computed via SplitIntoConvexPatches() or set via SetPatchIDs(), false if using + /// default (single patch). bool ArePatchesExplicitlySet() const { return patches_explicitly_set; } //////////////////////////////////////////////////////// diff --git a/src/DEM/MeshUtils.cpp b/src/DEM/MeshUtils.cpp index 5b5fb300..ad7bf725 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -339,8 +339,8 @@ size_t DEMMeshConnected::SplitIntoConvexPatches(float angle_threshold_deg) { void DEMMeshConnected::SetPatchIDs(const std::vector& patch_ids) { if (patch_ids.size() != nTri) { std::stringstream ss; - ss << "SetPatchIDs: Input vector size (" << patch_ids.size() - << ") must match the number of triangles (" << nTri << ") in the mesh." << std::endl; + ss << "SetPatchIDs: Input vector size (" << patch_ids.size() << ") must match the number of triangles (" << nTri + << ") in the mesh." << std::endl; throw std::runtime_error(ss.str()); } @@ -348,7 +348,7 @@ void DEMMeshConnected::SetPatchIDs(const std::vector& patch_ids) { for (size_t i = 0; i < patch_ids.size(); ++i) { if (patch_ids[i] < 0) { std::stringstream ss; - ss << "SetPatchIDs: Patch ID at index " << i << " is negative (" << patch_ids[i] + ss << "SetPatchIDs: Patch ID at index " << i << " is negative (" << patch_ids[i] << "). All patch IDs must be non-negative integers." << std::endl; throw std::runtime_error(ss.str()); } diff --git a/src/demo/DEMdemo_MeshPatch.cpp b/src/demo/DEMdemo_MeshPatch.cpp index 13936f0b..60b6076a 100644 --- a/src/demo/DEMdemo_MeshPatch.cpp +++ b/src/demo/DEMdemo_MeshPatch.cpp @@ -36,14 +36,17 @@ int main() { std::cout << "Loaded cube mesh successfully" << std::endl; std::cout << "Number of triangles: " << cube_mesh->GetNumTriangles() << std::endl; std::cout << "Number of vertices: " << cube_mesh->GetNumNodes() << std::endl; - + // Test default patch info (should be all in patch 0) std::cout << "\nDefault patch info (assuming convex mesh):" << std::endl; std::cout << "Number of patches: " << cube_mesh->GetNumPatches() << std::endl; std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; const auto& default_patch_ids = cube_mesh->GetPatchIDs(); - std::cout << "All triangles in patch 0: " << (std::all_of(default_patch_ids.begin(), default_patch_ids.end(), - [](int id) { return id == 0; }) ? "yes" : "no") << std::endl; + std::cout << "All triangles in patch 0: " + << (std::all_of(default_patch_ids.begin(), default_patch_ids.end(), [](int id) { return id == 0; }) + ? "yes" + : "no") + << std::endl; // Test different angle thresholds std::cout << "\n--- Test 2: Automatic Patch Splitting ---" << std::endl; @@ -54,7 +57,8 @@ int main() { std::cout << "\nAngle threshold: " << std::fixed << std::setprecision(1) << threshold << " degrees" << std::endl; std::cout << "Number of patches: " << num_patches << std::endl; - std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; + std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") + << std::endl; // Show patch distribution const auto& patch_ids = cube_mesh->GetPatchIDs(); @@ -70,7 +74,7 @@ int main() { std::cout << " Patch " << entry.first << ": " << entry.second << " triangles" << std::endl; } } - + // Test manual patch ID setting std::cout << "\n--- Test 3: Manual Patch ID Setting ---" << std::endl; size_t num_tris = cube_mesh->GetNumTriangles(); @@ -79,12 +83,12 @@ int main() { for (size_t i = 0; i < num_tris; ++i) { manual_patches[i] = i % 3; // Assign patches 0, 1, 2 cyclically } - + cube_mesh->SetPatchIDs(manual_patches); std::cout << "Manually set patch IDs (cycling 0, 1, 2)" << std::endl; std::cout << "Number of patches: " << cube_mesh->GetNumPatches() << std::endl; std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << std::endl; - + // Count triangles per patch const auto& manual_patch_ids = cube_mesh->GetPatchIDs(); std::map manual_patch_counts; @@ -139,8 +143,8 @@ int main() { std::cout << "\n--- Test 5: Empty Mesh ---" << std::endl; auto empty_mesh = std::make_shared(); std::cout << "Empty mesh default patches: " << empty_mesh->GetNumPatches() << " (expected: 1)" << std::endl; - std::cout << "Patches explicitly set: " << (empty_mesh->ArePatchesExplicitlySet() ? "yes" : "no") << " (expected: no)" - << std::endl; + std::cout << "Patches explicitly set: " << (empty_mesh->ArePatchesExplicitlySet() ? "yes" : "no") + << " (expected: no)" << std::endl; std::cout << "\n========================================" << std::endl; std::cout << "Demo completed successfully!" << std::endl;