Contact physics

examples/mechanics/crack_branching.cpp

@@ -114,7 +114,7 @@ void crackBranchingExample( const std::string filename )
     //                   Create solver
     // ====================================================
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
-        inputs, particles, force_model, prenotch );
+        inputs, particles, force_model );
 
     // ====================================================
     //                Boundary conditions
@@ -137,7 +137,7 @@ void crackBranchingExample( const std::string filename )
     // ====================================================
     //                   Simulation run
     // ====================================================
-    cabana_pd->init();
+    cabana_pd->init( bc, prenotch );
     cabana_pd->run( bc );
 }
 

examples/mechanics/fragmenting_cylinder.cpp

@@ -85,7 +85,8 @@ void fragmentingCylinderExample( const std::string filename )
     {
         double rsq = ( x[0] - x_center ) * ( x[0] - x_center ) +
                      ( x[1] - y_center ) * ( x[1] - y_center );
-        if ( rsq < Rin * Rin || rsq > Rout * Rout )
+        if ( rsq < Rin * Rin || rsq > Rout * Rout || x[2] > 0.05 ||
+             x[2] < -0.05 )
             return false;
         return true;
     };
@@ -137,12 +138,29 @@ void fragmentingCylinderExample( const std::string filename )
     particles->updateParticles( exec_space{}, init_functor );
 
     // ====================================================
-    //                   Simulation run
+    //  Simulation run with contact physics
     // ====================================================
-    auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
-        inputs, particles, force_model );
-    cabana_pd->init();
-    cabana_pd->run();
+    if ( inputs["use_contact"] )
+    {
+        double r_c = inputs["contact_horizon_factor"];
+        r_c *= dx[0];
+        CabanaPD::NormalRepulsionModel contact_model( delta, r_c, K );
+
+        auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
+            inputs, particles, force_model, contact_model );
+        cabana_pd->init();
+        cabana_pd->run();
+    }
+    // ====================================================
+    //  Simulation run without contact
+    // ====================================================
+    else
+    {
+        auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
+            inputs, particles, force_model );
+        cabana_pd->init();
+        cabana_pd->run();
+    }
 }
 
 // Initialize MPI+Kokkos.

examples/mechanics/inputs/fragmenting_cylinder.json

@@ -1,12 +1,15 @@
 {
     "num_cells"                : {"value": [51, 51, 101]},
-    "system_size"              : {"value": [0.05, 0.05, 0.1], "unit": "m"},
+    "dx":           {"value": [0.001, 0.001, 0.001]},
+    "system_size"              : {"value": [0.1, 0.1, 0.15], "unit": "m"},
     "grid_perturbation_factor" : {"value": 0.1},
     "density"                  : {"value": 7800, "unit": "kg/m^3"},
     "bulk_modulus"             : {"value": 130e+9, "unit": "Pa"},
     "shear_modulus"            : {"value": 78e+9, "unit": "Pa"},
     "critical_stretch"         : {"value": 0.02},
     "horizon"                  : {"value": 0.00417462, "unit": "m"},  
+    "use_contact"              : {"value": true},
+    "contact_horizon_factor"   : {"value": 0.9},  
     "cylinder_outer_radius"    : {"value": 0.025, "unit": "m"},
     "cylinder_inner_radius"    : {"value": 0.02, "unit": "m"},
     "max_radial_velocity"      : {"value": 200, "unit": "m/s"}, 

examples/mechanics/kalthoff_winkler.cpp

@@ -116,7 +116,7 @@ void kalthoffWinklerExample( const std::string filename )
     //                   Create solver
     // ====================================================
     auto cabana_pd = CabanaPD::createSolverFracture<memory_space>(
-        inputs, particles, force_model, prenotch );
+        inputs, particles, force_model );
 
     // ====================================================
     //                Boundary conditions
@@ -132,7 +132,7 @@ void kalthoffWinklerExample( const std::string filename )
     // ====================================================
     //                   Simulation run
     // ====================================================
-    cabana_pd->init();
+    cabana_pd->init( bc, prenotch );
     cabana_pd->run( bc );
 }
 

src/CabanaPD.hpp

@@ -29,8 +29,10 @@
 #include <CabanaPD_Types.hpp>
 #include <CabanaPD_config.hpp>
 
+#include <force/CabanaPD_ContactModels.hpp>
 #include <force/CabanaPD_ForceModels_LPS.hpp>
 #include <force/CabanaPD_ForceModels_PMB.hpp>
+#include <force/CabanaPD_Force_Contact.hpp>
 #include <force/CabanaPD_Force_LPS.hpp>
 #include <force/CabanaPD_Force_PMB.hpp>
 

src/CabanaPD_Comm.hpp

@@ -286,7 +286,9 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
 
         _init_timer.stop();
     }
-    ~Comm() {}
+
+    auto size() { return mpi_size; }
+    auto rank() { return mpi_rank; }
 
     // Determine which particles should be ghosted, reallocating and recounting
     // if needed.

src/CabanaPD_Force.hpp

@@ -156,7 +156,20 @@ class Force<MemorySpace, BaseForceModel>
     {
     }
 
-    // Primary constructor: use positions and construct neighbors.
+    // General constructor (necessary for contact, but could be used by any
+    // force routine).
+    template <class PositionType>
+    Force( const bool half_neigh, const double delta,
+           const PositionType& positions, const std::size_t num_local,
+           const double mesh_min[3], const double mesh_max[3],
+           const double tol = 1e-14 )
+        : _half_neigh( half_neigh )
+        , _neigh_list( neighbor_list_type( positions, 0, num_local, delta + tol,
+                                           1.0, mesh_min, mesh_max ) )
+    {
+    }
+
+    // Constructor which stores existing neighbors.
     template <class NeighborListType>
     Force( const bool half_neigh, const NeighborListType& neighbors )
         : _half_neigh( half_neigh )
@@ -224,7 +237,7 @@ class Force<MemorySpace, BaseForceModel>
 ******************************************************************************/
 template <class ForceType, class ParticleType, class ParallelType>
 void computeForce( ForceType& force, ParticleType& particles,
-                   const ParallelType& neigh_op_tag )
+                   const ParallelType& neigh_op_tag, const bool reset = true )
 {
     auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
@@ -233,7 +246,8 @@ void computeForce( ForceType& force, ParticleType& particles,
     auto f_a = particles.sliceForceAtomic();
 
     // Reset force.
-    Cabana::deep_copy( f, 0.0 );
+    if ( reset )
+        Cabana::deep_copy( f, 0.0 );
 
     // if ( half_neigh )
     // Forces must be atomic for half list
@@ -278,7 +292,7 @@ double computeEnergy( ForceType& force, ParticleType& particles,
 template <class ForceType, class ParticleType, class NeighborView,
           class ParallelType>
 void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
-                   const ParallelType& neigh_op_tag )
+                   const ParallelType& neigh_op_tag, const bool reset = true )
 {
     auto n_local = particles.n_local;
     auto x = particles.sliceReferencePosition();
@@ -287,7 +301,8 @@ void computeForce( ForceType& force, ParticleType& particles, NeighborView& mu,
     auto f_a = particles.sliceForceAtomic();
 
     // Reset force.
-    Cabana::deep_copy( f, 0.0 );
+    if ( reset )
+        Cabana::deep_copy( f, 0.0 );
 
     // if ( half_neigh )
     // Forces must be atomic for half list

src/CabanaPD_Input.hpp

@@ -32,22 +32,7 @@ class Inputs
         inputs = parse( filename );
 
         // Add additional derived inputs to json. System size.
-        auto size = inputs["system_size"]["value"];
-        std::string size_unit = inputs["system_size"]["unit"];
-        for ( std::size_t d = 0; d < size.size(); d++ )
-        {
-            double s = size[d];
-            double low = -0.5 * s;
-            double high = 0.5 * s;
-            inputs["low_corner"]["value"][d] = low;
-            inputs["high_corner"]["value"][d] = high;
-
-            double nc = inputs["num_cells"]["value"][d];
-            inputs["dx"]["value"][d] = ( high - low ) / nc;
-        }
-        inputs["low_corner"]["unit"] = size_unit;
-        inputs["high_corner"]["unit"] = size_unit;
-        inputs["dx"]["unit"] = size_unit;
+        setupSize();
 
         // Number of steps.
         double tf = inputs["final_time"]["value"];
@@ -104,7 +89,73 @@ class Inputs
         // Not yet a user option.
         inputs["half_neigh"]["value"] = false;
     }
-    ~Inputs() {}
+
+    void setupSize()
+    {
+        if ( inputs.contains( "system_size" ) )
+        {
+            auto size = inputs["system_size"]["value"];
+            std::string size_unit = inputs["system_size"]["unit"];
+            for ( std::size_t d = 0; d < size.size(); d++ )
+            {
+                double s = size[d];
+                double low = -0.5 * s;
+                double high = 0.5 * s;
+                inputs["low_corner"]["value"][d] = low;
+                inputs["high_corner"]["value"][d] = high;
+            }
+            inputs["low_corner"]["unit"] = size_unit;
+            inputs["high_corner"]["unit"] = size_unit;
+        }
+        else if ( inputs.contains( "low_corner" ) &&
+                  ( inputs.contains( "high_corner" ) ) )
+        {
+            auto low_corner = inputs["low_corner"]["value"];
+            auto high_corner = inputs["high_corner"]["value"];
+            std::string size_unit = inputs["low_corner"]["unit"];
+            for ( std::size_t d = 0; d < low_corner.size(); d++ )
+            {
+                double low = low_corner[d];
+                double high = high_corner[d];
+                inputs["system_size"]["value"][d] = high - low;
+            }
+            inputs["system_size"]["unit"] = size_unit;
+        }
+        else
+        {
+            throw std::runtime_error( "Must input either system_size or "
+                                      "both low_corner and high_corner." );
+        }
+
+        if ( inputs.contains( "dx" ) )
+        {
+            auto size = inputs["system_size"]["value"];
+            for ( std::size_t d = 0; d < size.size(); d++ )
+            {
+                double system_size = inputs["system_size"]["value"][d];
+                double dx = inputs["dx"]["value"][d];
+                inputs["num_cells"]["value"][d] =
+                    static_cast<int>( system_size / dx );
+            }
+        }
+        else if ( inputs.contains( "num_cells" ) )
+        {
+            auto size = inputs["system_size"]["value"];
+            for ( std::size_t d = 0; d < size.size(); d++ )
+            {
+                double low = inputs["low_corner"]["value"][d];
+                double high = inputs["low_corner"]["value"][d];
+                double nc = inputs["num_cells"]["value"][d];
+                inputs["dx"]["value"][d] = ( high - low ) / nc;
+            }
+            std::string size_unit = inputs["system_size"]["unit"];
+            inputs["dx"]["unit"] = size_unit;
+        }
+        else
+        {
+            throw std::runtime_error( "Must input either num_cells or dx." );
+        }
+    }
 
     void computeCriticalTimeStep()
     {