Update/dem

RELEASE_NOTES.md

@@ -28,7 +28,11 @@ Notable changes include:
   * Bug Fixes / improvements:
     * Fixed numerous compiler warnings with newer compilers such as G++ 9.4.
     * r-path for additional TPLs can be propogated to Spheral libraries with `SPHERAL_ADDITIONAL_RPATHS`.
-
+    * The DEM package has received a significant updated.
+      * Added simple analytic solid boundaries (planes,cylinder,sphere).
+      * User given more control over DEM fast time stepping.
+      * Simplified initialization and improved robustness.
+
 Version v2023.03.0 -- Release date 2023-03-29
 ==============================================
 

src/DEM/CMakeLists.txt

@@ -2,6 +2,13 @@ include_directories(.)
 
 set(DEM_inst 
     DEMBase
+    SolidBoundary/SolidBoundaryBase
+    SolidBoundary/DEMBoundaryPolicy
+    SolidBoundary/InfinitePlaneSolidBoundary
+    SolidBoundary/CircularPlaneSolidBoundary
+    SolidBoundary/RectangularPlaneSolidBoundary
+    SolidBoundary/CylinderSolidBoundary
+    SolidBoundary/SphereSolidBoundary
     ReplacePairFieldList
     IncrementPairFieldList
     ReplaceAndIncrementPairFieldList
@@ -14,15 +21,26 @@ set(DEM_sources
 
 instantiate(DEM_inst DEM_sources)
 
-set(DEM_headers 
+set(DEM_headers
+    SolidBoundary/SolidBoundaryBase.hh
+    SolidBoundary/DEMBoundaryPolicy.hh
+    SolidBoundary/InfinitePlaneSolidBoundary.hh
+    SolidBoundary/CircularPlaneSolidBoundary.hh
+    SolidBoundary/RectangularPlaneSolidBoundary.hh
+    SolidBoundary/CylinderSolidBoundary.hh
+    SolidBoundary/SphereSolidBoundary.hh
     ReplacePairFieldList.hh
     ReplaceAndIncrementPairFieldList.hh
     IncrementPairFieldList.hh
     LinearSpringDEM.hh
     DEMBase.hh
     DEMDimension.hh
     DEMFieldNames.hh
-    ContactStorageLocation.hh)
+    ContactStorageLocation.hh
+    setUniqueNodeIDs.hh
+    setUniqueNodeIDsInline.hh)
+
+add_subdirectory(SolidBoundary)
 
 spheral_install_python_files(DEM.py)
 

src/DEM/ContactStorageLocation.cc

@@ -1,15 +1,27 @@
-#include <vector>
+//---------------------------------Spheral++----------------------------------//
+// ContactIndex -- Simple structure used to track where pairwise variables 
+//                 are stored in the pairFieldLists. A -1 indicates an invalid
+//                 entry. The contact index has two configurations. The first
+//                 is for particle-particle contacts, for which store* and 
+//                 pair* properties are used. The second is for particle-
+//                 solid boundary interactions and there store* and 
+//                 solidBoundary are used. In both case, unused properties are
+//                 flagged as invalid (-1).
+//
+// J.M. Pearl 2022
+//----------------------------------------------------------------------------//
 
 #include "DEM/ContactStorageLocation.hh"
 
 namespace Spheral {
 
 ContactIndex::ContactIndex():
-  storeNodeList(0),
-  storeNode(0),
-  storeContact(0),
-  pairNodeList(0),
-  pairNode(0){
+  storeNodeList(-1),
+  storeNode(-1),
+  storeContact(-1),
+  pairNodeList(-1),
+  pairNode(-1),
+  solidBoundary(-1){
 }
 
 ContactIndex::ContactIndex(const int storeNodeListIndex, 
@@ -21,8 +33,22 @@ ContactIndex::ContactIndex(const int storeNodeListIndex,
   storeNode(storeNodeIndex),
   storeContact(storeContactIndex),
   pairNodeList(pairNodeListIndex),
-  pairNode(pairNodeIndex){
+  pairNode(pairNodeIndex),
+  solidBoundary(-1){
 }
 
+ContactIndex::ContactIndex(const int storeNodeListIndex, 
+                           const int storeNodeIndex, 
+                           const int storeContactIndex,
+                           const int solidBoundaryIndex):
+  storeNodeList(storeNodeListIndex),
+  storeNode(storeNodeIndex),
+  storeContact(storeContactIndex),
+  pairNodeList(-1),
+  pairNode(-1),
+  solidBoundary(solidBoundaryIndex){
+}
+
+
 
 }

src/DEM/ContactStorageLocation.hh

@@ -1,29 +1,41 @@
-
-
 //---------------------------------Spheral++----------------------------------//
-// DEM type alias for the particle rotation
+// ContactIndex -- Simple structure used to track where pairwise variables 
+//                 are stored in the pairFieldLists. A -1 indicates an invalid
+//                 entry. The contact index has two configurations. The first
+//                 is for particle-particle contacts, for which store* and 
+//                 pair* properties are used. The second is for particle-
+//                 solid boundary interactions and there store* and 
+//                 solidBoundary are used. In both case, unused properties are
+//                 flagged as invalid (-1).
+//
+// J.M. Pearl 2022
 //----------------------------------------------------------------------------//
 #ifndef __Spheral_ContactStorageLocation_hh__
 #define __Spheral_ContactStorageLocation_hh__
 
 namespace Spheral{
-
   struct ContactIndex {
     ContactIndex();
     ContactIndex(const int storeNodeListIndex, 
                  const int storeNodeIndex, 
                  const int storeContactIndex,
                  const int pairNodeListIndex,
                  const int pairNodeIndex);
+    ContactIndex(const int storeNodeListIndex, 
+                 const int storeNodeIndex, 
+                 const int storeContactIndex,
+                 const int solidBoundaryIndex);
 
       int storeNodeList; 
       int storeNode; 
       int storeContact;
 
       int pairNodeList;
       int pairNode;
-  };
 
+      int solidBoundary;
+  };
+}
 
 #else
 
@@ -33,4 +45,3 @@ namespace Spheral {
 }
 
 #endif
-}

src/DEM/DEM.py

@@ -23,6 +23,17 @@ def LinearSpringDEM(dataBase,
 
     assert dataBase.numDEMNodeLists == dataBase.numNodeLists, "all nodelists must be dem nodelists"
     assert stepsPerCollision > 1, "stepsPerCollision too low, reccomended is 25-50"
+    assert cohesiveTensileStrength >= 0, "cohesiveTensileStrength must be positive"
+    assert normalSpringConstant >= 0, "normalSpringConstant must be positive"
+    assert normalRestitutionCoefficient >= 0 and normalRestitutionCoefficient <= 1, "normalSpringConstant must be between 1 and 0"
+    assert tangentialSpringConstant >= 0, "normalSpringConstant must be positive"
+    assert tangentialRestitutionCoefficient >= 0 and tangentialRestitutionCoefficient <= 1, "normalSpringConstant must be between 1 and 0"
+    assert dynamicFrictionCoefficient >= 0, "dynamicFrictionCoefficient must be positive"
+    assert dynamicFrictionCoefficient >= 0, "dynamicFrictionCoefficient must be positive"
+    assert rollingFrictionCoefficient >= 0, "rollingFrictionCoefficient must be positive"
+    assert torsionalFrictionCoefficient >= 0, "torsionalFrictionCoefficient must be positive"
+
+    #if (stepsPerCollision < 10) print("WARNING: stepsPerCollision is very low, reccomended is 25-50")
 
     # we might want to allow the user to set less parameters with reasonable defaults
     #if tangentialSpringConstant is None: