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 c69e83dd..beb81b0a 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 = 1; + this->patches_explicitly_set = false; this->owner = NULL_BODYID; } @@ -523,6 +526,46 @@ 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 (defaults to 0 for all triangles, assuming convex mesh) + std::vector m_patch_ids; + // Number of patches in this mesh + 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 + /// 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); + + /// @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). 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. Default is 1 (assuming convex mesh). + 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). + 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 f0efa3cd..ad7bf725 100644 --- a/src/DEM/MeshUtils.cpp +++ b/src/DEM/MeshUtils.cpp @@ -32,6 +32,10 @@ #include #include #include +#include +#include +#include +#include #include "../kernel/DEMHelperKernels.cuh" #include "BdrsAndObjs.h" @@ -130,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; } @@ -183,4 +193,179 @@ void DEMMeshConnected::WriteWavefront(const std::string& filename, std::vector DEME_TINY_FLOAT) { + 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 / 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)}}; + + 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; +} + +// Split mesh into convex patches using region-growing algorithm. +// The algorithm groups adjacent triangles (sharing an edge) if the angle between their +// face normals is below the threshold. Each patch represents a locally convex region. +size_t DEMMeshConnected::SplitIntoConvexPatches(float angle_threshold_deg) { + if (nTri == 0) { + patches_explicitly_set = false; + num_patches = 1; + 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_explicitly_set = true; + + return num_patches; +} + +// Manually set patch IDs for each triangle +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; + 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/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..60b6076a --- /dev/null +++ b/src/demo/DEMdemo_MeshPatch.cpp @@ -0,0 +1,154 @@ +// 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 with Default Patch Info ---" << 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 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) { + 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; + std::cout << "Patches explicitly set: " << (cube_mesh->ArePatchesExplicitlySet() ? "yes" : "no") + << std::endl; + + // Show patch distribution + 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 4: 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->ArePatchesExplicitlySet()) { + 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 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 << "\n========================================" << std::endl; + std::cout << "Demo completed successfully!" << std::endl; + std::cout << "========================================" << std::endl; + + return 0; +}