Add mesh patch splitting with default convex assumption and manual patch assignment

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<unsigned int> m_anal_owner;
+    std::vector<bodyID_t> m_anal_owner;
     // Material types of these analytical geometries
     std::vector<materialsOffset_t> 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<float> m_anal_normals;
 
-    // These extra mesh facets' owners' ID will be appended to analytical entities'
-    std::vector<unsigned int> m_mesh_facet_owner;
+    // These mesh facets' owners' ID, flattened
+    std::vector<bodyID_t> m_mesh_facet_owner;
     // Material types of these mesh facets
     std::vector<materialsOffset_t> m_mesh_facet_materials;
-    // Material types of these mesh facets
+    // Three nodes of each triangle, flattened
     std::vector<DEMTriangle> m_mesh_facets;
 
+    // These mesh patches' owners' ID, flattened
+    std::vector<bodyID_t> m_mesh_patch_owner;
+    // Material types of these mesh patches
+    std::vector<materialsOffset_t> m_mesh_patch_materials;
+
     // Clump templates will be flatten and transferred into kernels upon Initialize()
     std::vector<float> m_template_clump_mass;
     std::vector<float3> m_template_clump_moi;

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(

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<int> 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<int>& 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<int>& 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
     ////////////////////////////////////////////////////////

src/DEM/MeshUtils.cpp

@@ -32,6 +32,10 @@
 #include <fstream>
 #include <map>
 #include <unordered_map>
+#include <cmath>
+#include <queue>
+#include <vector>
+#include <sstream>
 
 #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<D
     mf.close();
 }
 
+// Helper function to compute face normal for a triangle
+static float3 computeFaceNormal(const float3& v0, const float3& v1, const float3& v2) {
+    float3 edge1 = make_float3(v1.x - v0.x, v1.y - v0.y, v1.z - v0.z);
+    float3 edge2 = make_float3(v2.x - v0.x, v2.y - v0.y, v2.z - v0.z);
+
+    // Cross product
+    float3 normal = make_float3(edge1.y * edge2.z - edge1.z * edge2.y, edge1.z * edge2.x - edge1.x * edge2.z,
+                                edge1.x * edge2.y - edge1.y * edge2.x);
+
+    // Normalize
+    float length = std::sqrt(normal.x * normal.x + normal.y * normal.y + normal.z * normal.z);
+    if (length > 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<std::vector<size_t>> buildAdjacencyMap(const std::vector<int3>& face_v_indices) {
+    size_t num_faces = face_v_indices.size();
+    std::vector<std::vector<size_t>> adjacency(num_faces);
+
+    // Map from edge (as pair of vertex indices) to faces that share it
+    std::map<std::pair<int, int>, std::vector<size_t>> 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<int, int> 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<size_t>& 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<float3> 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<std::vector<size_t>> adjacency = buildAdjacencyMap(m_face_v_indices);
+
+    // Region growing algorithm to assign patches
+    int current_patch_id = 0;
+    std::vector<size_t> 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<int>& 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

src/demo/CMakeLists.txt

@@ -20,6 +20,7 @@ SET(DEMOS
 		DEMdemo_Sieve
 		DEMdemo_SingleSphereCollide
 		DEMdemo_MeshCollide
+		DEMdemo_MeshPatch
 		DEMdemo_BoxesFalling
 		DEMdemo_TestPack
 		DEMdemo_RotatingDrum