Merge pull request #5 from projectchrono/Mesh_Particles_Exp

Temporarily fix maxTriTriPenetration

Author: fl0rianr

src/DEM/APIPrivate.cpp

@@ -1117,8 +1117,8 @@ void DEMSolver::setSolverParams() {
     // dT->solverFlags.should_sort_pairs = should_sort_contacts;
 
     // Error out policies
-    kT->solverFlags.errOutAvgSphCnts = threshold_error_out_num_cnts;
-    dT->solverFlags.errOutAvgSphCnts = threshold_error_out_num_cnts;
+    kT->solverFlags.errOutAvgPrimitiveCnts = threshold_error_out_num_cnts;
+    dT->solverFlags.errOutAvgPrimitiveCnts = threshold_error_out_num_cnts;
     // simParams-stored variables need to be sync-ed to device
     kT->simParams->errOutBinSphNum = threshold_too_many_spheres_in_bin;
     dT->simParams->errOutBinSphNum = threshold_too_many_spheres_in_bin;

src/DEM/Structs.h

@@ -556,7 +556,7 @@ struct SolverFlags {
     // The max number of average contacts per sphere has before the solver errors out. The reason why I didn't use the
     // number of contacts for the sphere that has the most is that, well, we can have a huge sphere and it just will
     // have more contacts. But if avg cnt is high, that means probably the contact margin is out of control now.
-    float errOutAvgSphCnts = 300.;
+    float errOutAvgPrimitiveCnts = 300.;
 
     // Whether there are contacts that can never be removed.
     bool hasPersistentContacts = false;

src/algorithms/DEMContactDetection.cu

@@ -1110,7 +1110,7 @@ void contactDetection(std::shared_ptr<JitHelper::CachedProgram>& bin_sphere_kern
                 (*pNumUniqueNewA > 0) ? (float)(*scratchPad.numPrimitiveContacts) / (float)(*pNumUniqueNewA) : 0.0;
 
             DEME_DEBUG_PRINTF("Average number of contacts for each geometry: %.7g", stateParams.avgCntsPerPrimitive);
-            if (stateParams.avgCntsPerPrimitive > solverFlags.errOutAvgSphCnts) {
+            if (stateParams.avgCntsPerPrimitive > solverFlags.errOutAvgPrimitiveCnts) {
                 DEME_ERROR(
                     "On average a primitive (spheres, triangle facets) has %.7g contacts with other primitives, more "
                     "than the max allowance (%.7g).\nIf you believe this number is expected with the physics you are "
@@ -1123,7 +1123,7 @@ void contactDetection(std::shared_ptr<JitHelper::CachedProgram>& bin_sphere_kern
                     "a.k.a. elements initialized inside walls.\nIf none works and you are going to discuss this on "
                     "forum https://groups.google.com/g/projectchrono, please include a visual rendering of the "
                     "simulation before crash in the post.\n",
-                    stateParams.avgCntsPerPrimitive, solverFlags.errOutAvgSphCnts);
+                    stateParams.avgCntsPerPrimitive, solverFlags.errOutAvgPrimitiveCnts);
             }
 
             scratchPad.finishUsingTempVector("new_idA_runlength");

src/demo/DEMdemo_DrumCubes.cpp

@@ -30,8 +30,8 @@ int main() {
     DEMSim.SetNoForceRecord();
     DEMSim.SetMeshUniversalContact(true);
 
-    auto mat_type_cube = DEMSim.LoadMaterial({{"E", 1e9}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.5}, {"Crr", 0.01}});
-    auto mat_type_drum = DEMSim.LoadMaterial({{"E", 2e9}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.5}, {"Crr", 0.01}});
+    auto mat_type_cube = DEMSim.LoadMaterial({{"E", 1e6}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.5}, {"Crr", 0.01}});
+    auto mat_type_drum = DEMSim.LoadMaterial({{"E", 2e6}, {"nu", 0.3}, {"CoR", 0.6}, {"mu", 0.5}, {"Crr", 0.01}});
     DEMSim.SetMaterialPropertyPair("mu", mat_type_cube, mat_type_drum, 0.5);
 
     const float cube_size = 0.01f;
@@ -54,13 +54,17 @@ int main() {
     float IZZ = CylMass * CylRad * CylRad / 2;
     float IYY = (CylMass / 12) * (3 * CylRad * CylRad + CylHeight * CylHeight);
     auto Drum = DEMSim.AddExternalObject();
-    Drum->AddCylinder(CylCenter, CylAxis, CylRad, mat_type_drum, 0);
+    // Drum->AddCylinder(CylCenter, CylAxis, CylRad, mat_type_drum, 0);
+    Drum->AddPlane(make_float3(CylRad, 0, 0), make_float3(-1, 0, 0), mat_type_drum);
+    Drum->AddPlane(make_float3(-CylRad, 0, 0), make_float3(1, 0, 0), mat_type_drum);
+    Drum->AddPlane(make_float3(0, CylRad, 0), make_float3(0, -1, 0), mat_type_drum);
+    Drum->AddPlane(make_float3(0, -CylRad, 0), make_float3(0, 1, 0), mat_type_drum);
     Drum->SetMass(CylMass);
     Drum->SetMOI(make_float3(IYY, IYY, IZZ));
     auto Drum_tracker = DEMSim.Track(Drum);
     unsigned int drum_family = 100;
     Drum->SetFamily(drum_family);
-    const float rpm = 20.0f;
+    const float rpm = 200.0f;
     const float drum_ang_vel = rpm * 2.0f * PI / 60.0f;
     DEMSim.SetFamilyPrescribedAngVel(drum_family, "0", "0", to_string_with_precision(drum_ang_vel));
     auto top_bot_planes = DEMSim.AddExternalObject();
@@ -97,7 +101,7 @@ int main() {
     float max_v;
 
     DEMSim.InstructBoxDomainDimension(0.4, 0.4, 0.4);
-    float step_size = 5e-6f;
+    float step_size = 1e-4f;
     DEMSim.SetInitTimeStep(step_size);
     DEMSim.SetGPUTimersEnabled(true);
     DEMSim.SetGravitationalAcceleration(make_float3(0, 0, -9.81));
@@ -134,7 +138,7 @@ int main() {
                   << drum_torque.z << std::endl;
 
         float3 force_on_BC = planes_tracker->ContactAcc() * planes_tracker->Mass();
-        std::cout << "Contact force on bottom plane is " << force_on_BC.z << std::endl;
+        std::cout << "Contact force on bottom plane is " << std::abs(force_on_BC.z) << std::endl;
 
         DEMSim.DoDynamics(frame_time);
     }

src/kernel/DEMKinematicMisc.cu

@@ -82,7 +82,8 @@ __global__ void computeMarginFromAbsv_implTri(deme::DEMSimParams* simParams,
         // as our meshed particle representation is surface only, so we need to account for existing penetration length
         // in our future-proof contact detection, always.
         double penetrationMargin = *maxTriTriPenetration;
-        penetrationMargin = (meshUniversalContact && penetrationMargin > 0.0) ? penetrationMargin : 0.0;
+        //// TODO: Temporary measure
+        penetrationMargin = 0.;  // (meshUniversalContact && penetrationMargin > 0.0) ? penetrationMargin : 0.0;
         // Clamp penetration margin to the maximum allowed value to prevent super large margins
         if (penetrationMargin > simParams->capTriTriPenetration) {
             penetrationMargin = simParams->capTriTriPenetration;