Mesh particles unification #47

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged

Ruochun merged 19 commits into projectchrono:Mesh_Particles_Exp from uwsbel:Mesh_Particles

Feb 4, 2026

cmake/CudaSupportedArchitectures.cmake

-Original file line number
+Diff line change
@@ Expand Up / @@ -20,7 +20,6 @@ @@
     # version of the CUDA Toolkit
     #
     # Minimum CUDA version: 7.0
-    # Maximum CUDA version: 11.6
     function(cuda_supported_architectures)
@@ Expand All / @@ -31,7 +30,8 @@ function(cuda_supported_architectures) @@
     	set(cu10 30 35 50 52 60 61 70 72 75)
     	set(cu11 35 50 52 60 61 70 72 75 80)
     	set(cu11_x 35 50 52 60 61 70 72 75 80 86)
-    	set(cu12_x 50 52 60 61 70 72 75 80 86 89 120)
+    	set(cu12_x 50 52 60 61 70 72 75 80 86 89)
+    	set(cu12_8 50 52 60 61 70 72 75 80 86 89 120)
     	set(cu13_x 75 80 86 89 90 100 120 121)
     	if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 7)
@@ Expand Down Expand Up / @@ -62,6 +62,14 @@ function(cuda_supported_architectures) @@
     		set(CUDASUP_ARCHITECTURES ${cu12_x} CACHE INTERNAL "")
     	endif()
+        if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 13)
+            set(CUDASUP_ARCHITECTURES ${cu13_x} CACHE INTERNAL "")
+        endif()
+    	if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 12.8)
+    		set(CUDASUP_ARCHITECTURES ${cu12_8} CACHE INTERNAL "")
+    	endif()
         if (CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 13)
             set(CUDASUP_ARCHITECTURES ${cu13_x} CACHE INTERNAL "")
         endif()
@@ Expand Down @@

data/mesh/cross.stl

Binary file not shown.

data/mesh/cross_fine.stl

Binary file not shown.

src/DEM/API.h

-Original file line number
+Diff line change
@@ Expand Up / @@ -10,6 +10,7 @@ @@
     #include <set>
     #include <cfloat>
     #include <functional>
+    #include <thread>
     #include "kT.h"
     #include "dT.h"
@@ Expand Down Expand Up / @@ -1203,6 +1204,7 @@ class DEMSolver { @@
         /// Remove host-side cached vectors (so you can re-define them, and then re-initialize system)
         void ClearCache();
+        /// Output methods enqueue asynchronous writes; call WaitForPendingOutput() to block for completion.
         /// Write the current status of clumps to a file
         void WriteClumpFile(const std::string& outfilename, unsigned int accuracy = 10) const;
         void WriteClumpFile(const std::filesystem::path& outfilename, unsigned int accuracy = 10) const {
@@ Expand Down Expand Up / @@ -1231,6 +1233,8 @@ class DEMSolver { @@
         /// Write the current status of all meshes to a file.
         void WriteMeshFile(const std::string& outfilename) const;
         void WriteMeshFile(const std::filesystem::path& outfilename) const { WriteMeshFile(outfilename.string()); }
+        /// Wait for any in-flight async output to finish.
+        void WaitForPendingOutput() const;
         /// @brief Read 3 columns of your choice from a CSV filem and group them by clump_header.
         /// @param infilename CSV filename.
@@ Expand Down Expand Up / @@ -1471,6 +1475,8 @@ class DEMSolver { @@
         /// @brief Specify the output file format of meshes.
         /// @param format A choice between "VTK", "OBJ", "STL", "PLY".
         void SetMeshOutputFormat(const std::string& format);
+        /// @brief Enable/disable per-patch face colors in PLY mesh output (for testing auto patch splitting only).
+        void EnableMeshPatchColorOutput(bool enable = true);
         /// @brief Clear stored solver logs (errors, warnings, messages).
         void ClearLog() { Logger::GetInstance().Clear(); }
         /// @brief Show error and warnings.
@@ Expand Down Expand Up / @@ -1566,10 +1572,13 @@ class DEMSolver { @@
                                          CNT_OUTPUT_CONTENT::CNT_WILDCARD;
         // The output file format for meshes
         MESH_FORMAT m_mesh_out_format = MESH_FORMAT::VTK;
+        // If PLY mesh output should include per-patch face colors
+        bool m_mesh_out_ply_patch_colors = false;
         // If the solver should output wildcards to file
         bool m_is_out_owner_wildcards = false;
         bool m_is_out_cnt_wildcards = false;
         bool m_is_out_geo_wildcards = false;
+        mutable std::thread m_output_thread;
         // User-instructed simulation `world' size. Note it is an approximate of the true size and we will generate a world
         // not smaller than this. This is useful if the user want to automatically add BCs enclosing this user-defined
@@ Expand Down Expand Up / @@ -1741,6 +1750,8 @@ class DEMSolver { @@
         size_t nSpheresGM = 0;
         // Total number of triangle facets
         size_t nTriGM = 0;
+        // Total number of triangles that need neighbor info (compact neighbor array size)
+        size_t nTriNeighbors = 0;
         // Total number of mesh patches
         size_t nMeshPatches = 0;
         // Number of analytical entites (as components of some external objects)
@@ Expand Down Expand Up / @@ -1896,6 +1907,8 @@ class DEMSolver { @@
         std::vector<float3> m_input_mesh_obj_xyz;
         std::vector<float4> m_input_mesh_obj_rot;
         std::vector<unsigned int> m_input_mesh_obj_family;
+        std::vector<notStupidBool_t> m_input_mesh_obj_convex;
+        std::vector<notStupidBool_t> m_input_mesh_obj_never_winner;
         // Processed unique family prescription info
         std::vector<familyPrescription_t> m_unique_family_prescription;
@@ Expand Down Expand Up / @@ -1931,6 +1944,10 @@ class DEMSolver { @@
         std::vector<bodyID_t> m_mesh_facet_owner;
         // Patch ID for each mesh facet, flattened
         std::vector<bodyID_t> m_mesh_facet_patch;
+        // Per-facet edge neighbors (global triangle indices, NULL_BODYID if boundary)
+        std::vector<bodyID_t> m_mesh_facet_neighbor1;
+        std::vector<bodyID_t> m_mesh_facet_neighbor2;
+        std::vector<bodyID_t> m_mesh_facet_neighbor3;
         // Three nodes of each triangle, flattened
         std::vector<DEMTriangle> m_mesh_facets;
@@ Expand Down Expand Up / @@ -2065,6 +2082,7 @@ class DEMSolver { @@
                                    size_t nSpheres,
                                    size_t nTriMesh,
                                    size_t nFacets,
+                                   size_t nTriNeighbors,
                                    size_t nMeshPatches,
                                    unsigned int nExtObj_old,
                                    unsigned int nAnalGM_old);
@@ Expand Down @@

src/DEM/APIPrivate.cpp

            
                      Original file line number
                      Diff line number
                      Diff line change
                  
    @@ -13,13 +13,93 @@
  
    #include <thread>

    #include <chrono>

    #include <cstring>

    #include <cstdint>

    #include <cmath>

    #include <limits>

    #include <algorithm>

    #include <map>

    #include <tuple>

    #include <unordered_map>

    #include <array>

    namespace deme {

    namespace {

    struct EdgeInfo {

        size_t tri = 0;

        int edge = 0;

    };

    inline uint64_t makeEdgeKey(int a, int b) {

        const uint32_t lo = static_cast<uint32_t>(std::min(a, b));

        const uint32_t hi = static_cast<uint32_t>(std::max(a, b));

        return (static_cast<uint64_t>(lo) << 32) | static_cast<uint64_t>(hi);

    }

    std::vector<std::array<bodyID_t, 3>> buildTriangleEdgeNeighbors(const std::vector<int3>& face_v_indices,

                                                                     const std::vector<float3>& vertices) {

        const size_t n_faces = face_v_indices.size();

        std::vector<std::array<bodyID_t, 3>> neighbors(n_faces, {NULL_BODYID, NULL_BODYID, NULL_BODYID});

        if (n_faces == 0) {

            return neighbors;

        }

        std::vector<size_t> canon;

        if (!vertices.empty()) {

            const double eps = computeVertexQuantEps(vertices);

            canon = buildCanonicalVertexMap(vertices, eps);

        }

        std::unordered_map<uint64_t, std::vector<EdgeInfo>> edge_map;

        edge_map.reserve(n_faces * 3);

        for (size_t i = 0; i < n_faces; ++i) {

            const int3& face = face_v_indices[i];

            const int v0_raw = face.x;

            const int v1_raw = face.y;

            const int v2_raw = face.z;

            if (v0_raw < 0 || v1_raw < 0 || v2_raw < 0) {

                continue;

            }

            int v0 = v0_raw;

            int v1 = v1_raw;

            int v2 = v2_raw;

            if (!canon.empty()) {

                if (static_cast<size_t>(v0_raw) >= canon.size() || static_cast<size_t>(v1_raw) >= canon.size() ||

                    static_cast<size_t>(v2_raw) >= canon.size()) {

                    continue;

                }

                v0 = static_cast<int>(canon[static_cast<size_t>(v0_raw)]);

                v1 = static_cast<int>(canon[static_cast<size_t>(v1_raw)]);

                v2 = static_cast<int>(canon[static_cast<size_t>(v2_raw)]);

            }

            if (v0 == v1 || v1 == v2 || v2 == v0) {

                continue;

            }

            const uint64_t e0 = makeEdgeKey(v0, v1);

            const uint64_t e1 = makeEdgeKey(v1, v2);

            const uint64_t e2 = makeEdgeKey(v2, v0);

            edge_map[e0].push_back(EdgeInfo{i, 0});

            edge_map[e1].push_back(EdgeInfo{i, 1});

            edge_map[e2].push_back(EdgeInfo{i, 2});

        }

        for (const auto& entry : edge_map) {

            const auto& info = entry.second;

            if (info.size() == 2) {

                const EdgeInfo& a = info[0];

                const EdgeInfo& b = info[1];

                neighbors[a.tri][a.edge] = static_cast<bodyID_t>(b.tri);

                neighbors[b.tri][b.edge] = static_cast<bodyID_t>(a.tri);

            }

        }

        return neighbors;

    }

    }  // namespace

    void DEMSolver::assertSysInit(const std::string& method_name) {

        if (!sys_initialized) {

            DEME_ERROR("DEMSolver's method %s can only be called after calling Initialize()", method_name.c_str());

    @@ -717,6 +797,11 @@ void DEMSolver::preprocessAnalyticalObjs() {
  
                        addAnalCompTemplate(ANAL_OBJ_TYPE_CYL_INF, comp_mat.at(i), thisLoadExtObj, param.cyl.center,

                                            param.cyl.dir, param.cyl.radius, 0, 0, param.cyl.normal);

                        break;

                    case OBJ_COMPONENT::PLANAR_CYL:

                        addAnalCompTemplate(ANAL_OBJ_TYPE_PLANAR_CYL, comp_mat.at(i), thisLoadExtObj, param.cyl.center,

                                            param.cyl.dir, param.cyl.radius, 0, 0,

                                            param.cyl.normal);

                        break;

                    default:

                        DEME_ERROR(std::string("There is at least one analytical boundary that has a type not supported."));

                }

    @@ -835,17 +920,43 @@ void DEMSolver::preprocessTriangleObjs() {
  
            m_input_mesh_obj_xyz.push_back(mesh_obj->init_pos);

            m_input_mesh_obj_rot.push_back(mesh_obj->init_oriQ);

            m_input_mesh_obj_family.push_back(mesh_obj->family_code);

            m_input_mesh_obj_convex.push_back(mesh_obj->is_convex ? 1 : 0);

            m_input_mesh_obj_never_winner.push_back(mesh_obj->never_winner ? 1 : 0);

            m_mesh_facet_owner.insert(m_mesh_facet_owner.end(), mesh_obj->GetNumTriangles(), thisMeshObj);

            // Initialize patch IDs if not already set (default: all facets in patch 0)

            if (!mesh_obj->patches_explicitly_set && mesh_obj->m_patch_ids.empty()) {

                mesh_obj->SetPatchIDs({0});

            const bodyID_t tri_offset = static_cast<bodyID_t>(m_mesh_facets.size());

            std::vector<std::array<bodyID_t, 3>> local_neighbors;

            if (mesh_obj->IsConvex() && mesh_obj->IsNeverWinner()) {

                local_neighbors.assign(mesh_obj->GetNumTriangles(), {NULL_BODYID, NULL_BODYID, NULL_BODYID});

            } else {

                local_neighbors = buildTriangleEdgeNeighbors(mesh_obj->m_face_v_indices, mesh_obj->m_vertices);

            }

            // Force single-patch semantics: one patch per mesh (all facets in patch 0)

            if (mesh_obj->patches_explicitly_set || mesh_obj->GetNumPatches() > 1) {

                DEME_WARNING(

                    "Mesh patch IDs were provided or computed, but single-patch mode is enabled; all facets will be "

                    "assigned to one patch.");

            }

            if (mesh_obj->GetNumTriangles() > 0) {

                mesh_obj->m_patch_ids.assign(mesh_obj->GetNumTriangles(), 0);

            } else {

                mesh_obj->m_patch_ids.clear();

            }

            mesh_obj->nPatches = 1;

            mesh_obj->patches_explicitly_set = true;

            mesh_obj->m_patch_locations.clear();

            mesh_obj->patch_locations_explicitly_set = false;

            if (mesh_obj->materials.size() != 1 && !mesh_obj->materials.empty()) {

                auto mat = mesh_obj->materials[0];

                mesh_obj->materials.assign(1, mat);

                DEME_WARNING("Mesh provided multiple patch materials; single-patch mode keeps only the first material.");

            }

            // Populate patch owner and material arrays (one entry per patch in this mesh)

            // Note patch_id in a mesh is always 0-based, and contiguous

            std::vector<materialsOffset_t> patch_materials(mesh_obj->GetNumPatches());

            for (size_t facet_idx = 0; facet_idx < mesh_obj->GetNumPatches(); facet_idx++) {

            for (size_t facet_idx = 0; facet_idx < mesh_obj->GetNumTriangles(); facet_idx++) {

                // patch_id is per-triangle

                bodyID_t patch_id = mesh_obj->m_patch_ids.at(facet_idx);

                // Assign this facet's material to its patch (will overwrite for each facet, but they should be consistent

    @@ -881,6 +992,11 @@ void DEMSolver::preprocessTriangleObjs() {
  
                    }

                }

                m_mesh_facets.push_back(tri);

                const auto& nb = local_neighbors[i];

                m_mesh_facet_neighbor1.push_back(nb[0] == NULL_BODYID ? NULL_BODYID : nb[0] + tri_offset);

                m_mesh_facet_neighbor2.push_back(nb[1] == NULL_BODYID ? NULL_BODYID : nb[1] + tri_offset);

                m_mesh_facet_neighbor3.push_back(nb[2] == NULL_BODYID ? NULL_BODYID : nb[2] + tri_offset);

            }

            thisLoadPatchCount += mesh_obj->GetNumPatches();

    @@ -1310,17 +1426,30 @@ void DEMSolver::setSimParams() {
  
    }

    void DEMSolver::allocateGPUArrays() {

        size_t tri_neighbors = 0;

        for (const auto& mesh_obj : cached_mesh_objs) {

            if (!mesh_obj) {

                continue;

            }

            if (!(mesh_obj->IsConvex() && mesh_obj->IsNeverWinner())) {

                tri_neighbors += mesh_obj->GetNumTriangles();

            }

        }

        nTriNeighbors = tri_neighbors;

        // Resize arrays based on the statistical data we have

        std::thread dThread = std::move(std::thread([this]() {

            this->dT->allocateGPUArrays(this->nOwnerBodies, this->nOwnerClumps, this->nExtObj, this->nTriMeshes,

                                        this->nSpheresGM, this->nTriGM, this->nMeshPatches, this->nAnalGM,

                                        this->nSpheresGM, this->nTriGM, this->nTriNeighbors, this->nMeshPatches,

                                        this->nAnalGM,

                                        this->nExtraContacts, this->nDistinctMassProperties,

                                        this->nDistinctClumpBodyTopologies, this->nDistinctClumpComponents,

                                        this->nJitifiableClumpComponents, this->nMatTuples);

        }));

        std::thread kThread = std::move(std::thread([this]() {

            this->kT->allocateGPUArrays(this->nOwnerBodies, this->nOwnerClumps, this->nExtObj, this->nTriMeshes,

                                        this->nSpheresGM, this->nTriGM, this->nAnalGM, this->nExtraContacts,

                                        this->nSpheresGM, this->nTriGM, this->nTriNeighbors, this->nAnalGM,

                                        this->nExtraContacts,

                                        this->nDistinctMassProperties, this->nDistinctClumpBodyTopologies,

                                        this->nDistinctClumpComponents, this->nJitifiableClumpComponents, this->nMatTuples);

        }));

    @@ -1340,8 +1469,10 @@ void DEMSolver::initializeGPUArrays() {
  
            // Analytical objects' initial stats

            m_input_ext_obj_xyz, m_input_ext_obj_rot, m_input_ext_obj_family,

            // Meshed objects' initial stats

            cached_mesh_objs, m_input_mesh_obj_xyz, m_input_mesh_obj_rot, m_input_mesh_obj_family, m_mesh_facet_owner,

            m_mesh_facet_patch, m_mesh_facets, m_mesh_patch_owner, m_mesh_patch_materials,

            cached_mesh_objs, m_input_mesh_obj_xyz, m_input_mesh_obj_rot, m_input_mesh_obj_family,

            m_input_mesh_obj_convex, m_input_mesh_obj_never_winner, m_mesh_facet_owner, m_mesh_facet_patch,

            m_mesh_facet_neighbor1, m_mesh_facet_neighbor2, m_mesh_facet_neighbor3, m_mesh_facets, m_mesh_patch_owner,

            m_mesh_patch_materials,

            // Clump template name mapping

            m_template_number_name_map,

            // Clump template info (mass, sphere components, materials etc.)

    @@ -1363,7 +1494,8 @@ void DEMSolver::initializeGPUArrays() {
  
            // Analytical objects' initial stats

            m_input_ext_obj_family,

            // Meshed objects' initial stats

            m_input_mesh_obj_family, m_mesh_facet_owner, m_mesh_facet_patch, m_mesh_facets,

            m_input_mesh_obj_family, m_input_mesh_obj_convex, m_input_mesh_obj_never_winner, m_mesh_facet_owner,

            m_mesh_facet_patch, m_mesh_facet_neighbor1, m_mesh_facet_neighbor2, m_mesh_facet_neighbor3, m_mesh_facets,

            // Analytical obj physics properties

            m_ext_obj_comp_num,

            // Family mask

    @@ -1380,6 +1512,7 @@ void DEMSolver::updateClumpMeshArrays(size_t nOwners,
  
                                          size_t nSpheres,

                                          size_t nTriMesh,

                                          size_t nFacets,

                                          size_t nTriNeighbors,

                                          size_t nMeshPatches,

                                          unsigned int nExtObj,

                                          unsigned int nAnalGM) {

    @@ -1393,8 +1526,10 @@ void DEMSolver::updateClumpMeshArrays(size_t nOwners,
  
            // Analytical objects' initial stats

            m_input_ext_obj_xyz, m_input_ext_obj_rot, m_input_ext_obj_family,

            // Meshed objects' initial stats

            cached_mesh_objs, m_input_mesh_obj_xyz, m_input_mesh_obj_rot, m_input_mesh_obj_family, m_mesh_facet_owner,

            m_mesh_facet_patch, m_mesh_facets, m_mesh_patch_owner, m_mesh_patch_materials,

            cached_mesh_objs, m_input_mesh_obj_xyz, m_input_mesh_obj_rot, m_input_mesh_obj_family,

            m_input_mesh_obj_convex, m_input_mesh_obj_never_winner, m_mesh_facet_owner, m_mesh_facet_patch,

            m_mesh_facet_neighbor1, m_mesh_facet_neighbor2, m_mesh_facet_neighbor3, m_mesh_facets, m_mesh_patch_owner,

            m_mesh_patch_materials,

            // Clump template info (mass, sphere components, materials etc.)

            flattened_clump_templates,

            // Analytical obj physics properties

    @@ -1408,22 +1543,23 @@ void DEMSolver::updateClumpMeshArrays(size_t nOwners,
  
            // I/O and misc.

            m_no_output_families, m_tracked_objs,

            // Number of entities, old

            nOwners, nClumps, nSpheres, nTriMesh, nFacets, nMeshPatches, nExtObj, nAnalGM);

            nOwners, nClumps, nSpheres, nTriMesh, nFacets, nTriNeighbors, nMeshPatches, nExtObj, nAnalGM);

        kT->updateClumpMeshArrays(

            // Clump batchs' initial stats

            cached_input_clump_batches,

            // Analytical objects' initial stats

            m_input_ext_obj_family,

            // Meshed objects' initial stats

            m_input_mesh_obj_family, m_mesh_facet_owner, m_mesh_facet_patch, m_mesh_facets,

            m_input_mesh_obj_family, m_input_mesh_obj_convex, m_input_mesh_obj_never_winner, m_mesh_facet_owner,

            m_mesh_facet_patch, m_mesh_facet_neighbor1, m_mesh_facet_neighbor2, m_mesh_facet_neighbor3, m_mesh_facets,

            // Analytical obj physics properties

            m_ext_obj_comp_num,

            // Family mask

            m_family_mask_matrix,

            // Templates and misc.

            flattened_clump_templates,

            // Number of entities, old

            nOwners, nClumps, nSpheres, nTriMesh, nFacets, nMeshPatches, nExtObj, nAnalGM);

            nOwners, nClumps, nSpheres, nTriMesh, nFacets, nTriNeighbors, nMeshPatches, nExtObj, nAnalGM);

    }

    void DEMSolver::packDataPointers() {

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Mesh particles unification #47

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!