Replace SAT with bidirectional projection and axis hygiene check in triangle-triangle collision detection

src/demo/CMakeLists.txt

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

src/demo/DEMdemo_BoxesFalling.cpp

@@ -0,0 +1,136 @@
+//  Copyright (c) 2021, SBEL GPU Development Team
+//  Copyright (c) 2021, University of Wisconsin - Madison
+//
+//	SPDX-License-Identifier: BSD-3-Clause
+
+// =============================================================================
+// A demo featuring multiple boxes falling onto an analytical plane, testing
+// triangle-triangle contacts between boxes and box-plane contacts.
+// =============================================================================
+
+#include <core/ApiVersion.h>
+#include <core/utils/ThreadManager.h>
+#include <DEM/API.h>
+#include <DEM/utils/Samplers.hpp>
+
+#include <filesystem>
+#include <cstdio>
+#include <time.h>
+#include <random>
+
+using namespace deme;
+using namespace std::filesystem;
+
+int main() {
+    DEMSolver DEMSim;
+    DEMSim.SetVerbosity("INFO");
+    DEMSim.SetOutputFormat(OUTPUT_FORMAT::CSV);
+    DEMSim.SetContactOutputContent({"OWNER", "FORCE", "POINT", "NORMAL", "TORQUE", "CNT_WILDCARD"});
+    DEMSim.InstructBoxDomainDimension(20, 20, 15);
+
+    // Enable mesh-mesh contacts so boxes can collide with each other
+    DEMSim.SetMeshUniversalContact(true);
+
+    // Define material properties
+    auto mat_box = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.4}, {"Crr", 0.01}});
+    auto mat_plane = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.5}, {"mu", 0.3}, {"Crr", 0.01}});
+
+    // Add a bottom plane at z = 0
+    DEMSim.AddBCPlane(make_float3(0, 0, 0), make_float3(0, 0, 1), mat_plane);
+
+    // Create multiple boxes in a grid pattern above the plane
+    const int num_boxes_x = 4;
+    const int num_boxes_y = 4;
+    const float box_spacing = 2.0;
+    const float initial_height = 8.0;
+    const float box_size = 0.5;  // Assume cube.obj is roughly 1x1x1, scale to 0.5
+
+    std::vector<std::shared_ptr<DEMTracker>> trackers;
+
+    // Random number generator for small position variations
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<> pos_dist(-0.1, 0.1);
+    std::uniform_real_distribution<> rot_dist(-0.2, 0.2);
+
+    for (int i = 0; i < num_boxes_x; i++) {
+        for (int j = 0; j < num_boxes_y; j++) {
+            float x = (i - num_boxes_x / 2.0 + 0.5) * box_spacing + pos_dist(gen);
+            float y = (j - num_boxes_y / 2.0 + 0.5) * box_spacing + pos_dist(gen);
+            float z = initial_height + (i + j) * 0.5 + pos_dist(gen) * 2;
+
+            auto box = DEMSim.AddWavefrontMeshObject((GET_DATA_PATH() / "mesh/cube.obj").string(), mat_box);
+            box->SetFamily(0);
+            box->SetInitPos(make_float3(x, y, z));
+            box->Scale(box_size);  // Scale the box
+
+            // Set mass and MOI for the box (scaled cubic)
+            float mass = 1000.0 * box_size * box_size * box_size;
+            float moi = mass * box_size * box_size / 6.0;  // MOI for cube
+            box->SetMass(mass);
+            box->SetMOI(make_float3(moi, moi, moi));
+
+            // Add small initial rotation for more interesting dynamics
+            box->SetInitQuat(make_float4(rot_dist(gen), rot_dist(gen), rot_dist(gen), 1.0));
+
+            auto tracker = DEMSim.Track(box);
+            trackers.push_back(tracker);
+        }
+    }
+
+    float step_time = 1e-5;
+    DEMSim.SetInitTimeStep(step_time);
+    DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));
+    DEMSim.SetExpandSafetyType("auto");
+    DEMSim.Initialize();
+
+    // Setup output directory
+    path out_dir = current_path();
+    out_dir /= "DemoOutput_BoxesFalling";
+    create_directory(out_dir);
+
+    float frame_time = 1e-2;
+    int frame = 0;
+    int frame_step = (int)(frame_time / step_time);
+    double final_time = 3.0;  // 3 seconds of simulation
+
+    std::cout << "Starting simulation with " << num_boxes_x * num_boxes_y << " boxes" << std::endl;
+    std::cout << "Frame time: " << frame_time << " s, Time step: " << step_time << " s" << std::endl;
+    std::cout << "Total frames: " << (int)(final_time / frame_time) << std::endl;
+
+    for (int i = 0; i < (int)(final_time / step_time); i++) {
+        if (i % frame_step == 0) {
+            frame++;
+            std::cout << "Frame: " << frame << " (t = " << i * step_time << " s)" << std::endl;
+
+            char meshfilename[100];
+            sprintf(meshfilename, "DEMdemo_boxes_%04d.vtk", frame);
+            DEMSim.WriteMeshFile(out_dir / meshfilename);
+
+            // Print some statistics
+            int boxes_settled = 0;
+            float avg_height = 0.0;
+            for (size_t k = 0; k < trackers.size(); k++) {
+                float3 pos = trackers[k]->Pos();
+                avg_height += pos.z;
+                if (pos.z < 2.0) {  // Consider boxes below 2m as settled
+                    boxes_settled++;
+                }
+            }
+            avg_height /= trackers.size();
+
+            std::cout << "  Average height: " << avg_height << " m" << std::endl;
+            std::cout << "  Boxes settled: " << boxes_settled << " / " << trackers.size() << std::endl;
+
+            DEMSim.ShowMemStats();
+            std::cout << "----------------------------------------" << std::endl;
+        }
+
+        DEMSim.DoDynamics(step_time);
+    }
+
+    DEMSim.ShowThreadCollaborationStats();
+    DEMSim.ShowTimingStats();
+    std::cout << "DEMdemo_BoxesFalling exiting..." << std::endl;
+    return 0;
+}

src/kernel/DEMCalcForceKernels.cu

@@ -50,7 +50,7 @@ __device__ __forceinline__ void calculateContactForcesImpl(deme::DEMSimParams* s
     deme::contact_t ContactType = CONTACT_TYPE;
     // The following quantities are always calculated, regardless of force model
     double3 contactPnt;
-    float3 B2A;  // Unit vector pointing from body B to body A (contact normal)
+    float3 B2A;  // contact normal felt by A, pointing from B to A, unit vector
     // Penetration depth
     // Positive number for overlapping cases, negative for non-overlapping cases
     double overlapDepth = 0.0;
@@ -267,14 +267,10 @@ __device__ __forceinline__ void calculateContactForcesImpl(deme::DEMSimParams* s
             bool shouldDropContact = false;
             bool in_contact = checkTriangleTriangleOverlap<double3, double>(
                 triANode1, triANode2, triANode3, triBNode1, triBNode2, triBNode3, contact_normal, overlapDepth,
-                overlapArea, contactPnt, shouldDropContact, needsNonContactPenetrationCalc);
+                overlapArea, contactPnt, needsNonContactPenetrationCalc);
             B2A = to_float3(contact_normal);
             // Fix ContactType if needed
-            // If the solver thinks this contact is redundant, we drop it directly
-            if (shouldDropContact) {
-                ContactType = deme::NOT_A_CONTACT;
-            }
-            // Then, if the solver says in contact, we do not question it
+            // If the solver says in contact, we do not question it
             if (!in_contact) {
                 // Then, if we have extra margin, we check that if the distance is within the extra margin.
                 // Or, if the solver says no contact and we got overlapDepth > 0, that is by design a flag telling the
@@ -283,9 +279,9 @@ __device__ __forceinline__ void calculateContactForcesImpl(deme::DEMSimParams* s
                     ContactType = deme::NOT_A_CONTACT;
                 }
             }
-            // if (ContactType != deme::NOT_A_CONTACT) {
-            //     printf("Triangle-Triangle contact retained: overlapDepth = %f\n", overlapDepth);
-            // }
+            if (overlapDepth > 0.5) {
+                ContactType = deme::NOT_A_CONTACT;
+            }
         }
 
     } else if constexpr (BType == deme::GEO_T_ANALYTICAL) {

src/kernel/DEMCollisionKernels_SphSph.cuh

@@ -16,22 +16,22 @@
  *
  */
 template <typename T1, typename T2>
-inline __host__ __device__ deme::contact_t checkSpheresOverlap(const T1& XA,
-                                                               const T1& YA,
-                                                               const T1& ZA,
-                                                               const T1& radA,
-                                                               const T1& XB,
-                                                               const T1& YB,
-                                                               const T1& ZB,
-                                                               const T1& radB,
-                                                               T1& CPX,
-                                                               T1& CPY,
-                                                               T1& CPZ,
-                                                               T2& normalX,
-                                                               T2& normalY,
-                                                               T2& normalZ,
-                                                               T1& overlapDepth,
-                                                               T1& overlapArea) {
+__host__ __device__ deme::contact_t checkSpheresOverlap(const T1& XA,
+                                                        const T1& YA,
+                                                        const T1& ZA,
+                                                        const T1& radA,
+                                                        const T1& XB,
+                                                        const T1& YB,
+                                                        const T1& ZB,
+                                                        const T1& radB,
+                                                        T1& CPX,
+                                                        T1& CPY,
+                                                        T1& CPZ,
+                                                        T2& normalX,
+                                                        T2& normalY,
+                                                        T2& normalZ,
+                                                        T1& overlapDepth,
+                                                        T1& overlapArea) {
     T1 centerDist2 = distSquared<T1>(XA, YA, ZA, XB, YB, ZB);
     deme::contact_t contactType;
     if (centerDist2 > (radA + radB) * (radA + radB)) {
@@ -60,20 +60,20 @@ inline __host__ __device__ deme::contact_t checkSpheresOverlap(const T1& XA,
 // normal. Returned contact type is only useful for kT to sort contact types, as for dT's force calculation, the flavor
 // used is determined by type B's actual objType.
 template <typename T1, typename T2, typename T3>
-inline __host__ __device__ deme::contact_t checkSphereEntityOverlap(const T1& A,
-                                                                    const T2& radA,
-                                                                    const deme::objType_t& typeB,
-                                                                    const T1& B,
-                                                                    const float3& dirB,
-                                                                    const float& size1B,
-                                                                    const float& size2B,
-                                                                    const float& size3B,
-                                                                    const float& normal_sign,
-                                                                    const float& beta4Entity,
-                                                                    T1& CP,
-                                                                    float3& cntNormal,
-                                                                    T3& overlapDepth,
-                                                                    T3& overlapArea) {
+__host__ __device__ deme::contact_t checkSphereEntityOverlap(const T1& A,
+                                                             const T2& radA,
+                                                             const deme::objType_t& typeB,
+                                                             const T1& B,
+                                                             const float3& dirB,
+                                                             const float& size1B,
+                                                             const float& size2B,
+                                                             const float& size3B,
+                                                             const float& normal_sign,
+                                                             const float& beta4Entity,
+                                                             T1& CP,
+                                                             float3& cntNormal,
+                                                             T3& overlapDepth,
+                                                             T3& overlapArea) {
     deme::contact_t contactType;
     switch (typeB) {
         case (deme::ANAL_OBJ_TYPE_PLANE): {