From 0508b9a04b5a38416830930f90d08e39a39f44e5 Mon Sep 17 00:00:00 2001
From: Sam Reeve <6740307+streeve@users.noreply.github.com>
Date: Fri, 20 Oct 2023 12:09:59 -0400
Subject: [PATCH] Move BC to elastic Will be used for temperature, etc.

---
 examples/elastic_wave.cpp |  5 ++++-
 src/CabanaPD_Boundary.hpp | 32 ++++++++++++++++++++++++++++----
 src/CabanaPD_Solver.hpp   | 35 ++++++++++++++++++++++-------------
 3 files changed, 54 insertions(+), 18 deletions(-)

diff --git a/examples/elastic_wave.cpp b/examples/elastic_wave.cpp
index 61c46ed2..aa2cac54 100644
--- a/examples/elastic_wave.cpp
+++ b/examples/elastic_wave.cpp
@@ -89,8 +89,11 @@ int main( int argc, char* argv[] )
         };
         particles->updateParticles( exec_space{}, init_functor );
 
+        auto bc = CabanaPD::createBoundaryCondition<memory_space>(
+            CabanaPD::ZeroBCTag{} );
+
         auto cabana_pd = CabanaPD::createSolverElastic<memory_space>(
-            inputs, particles, force_model );
+            inputs, particles, force_model, bc );
         cabana_pd->init_force();
         cabana_pd->run();
 
diff --git a/src/CabanaPD_Boundary.hpp b/src/CabanaPD_Boundary.hpp
index 2e45522f..06c0326b 100644
--- a/src/CabanaPD_Boundary.hpp
+++ b/src/CabanaPD_Boundary.hpp
@@ -19,10 +19,6 @@
 
 namespace CabanaPD
 {
-// Empty boundary.
-struct ZeroBoundary
-{
-};
 
 // Define a plane or other rectilinear subset of the system as the boundary.
 struct RegionBoundary
@@ -56,6 +52,9 @@ struct BoundaryIndexSpace<MemorySpace, RegionBoundary>
     index_view_type _view;
     index_view_type _count;
 
+    // Default for empty case.
+    BoundaryIndexSpace() {}
+
     template <class ExecSpace, class Particles>
     BoundaryIndexSpace( ExecSpace exec_space, Particles particles,
                         std::vector<RegionBoundary> planes,
@@ -137,9 +136,27 @@ struct ForceCrackBranchBCTag
 {
 };
 
+struct ZeroBCTag
+{
+};
+
 template <class BCIndexSpace, class BCTag>
 struct BoundaryCondition;
 
+template <class BCIndexSpace>
+struct BoundaryCondition<BCIndexSpace, ZeroBCTag>
+{
+    template <class ExecSpace, class Particles>
+    void update( ExecSpace, Particles, RegionBoundary )
+    {
+    }
+
+    template <class ExecSpace, class ParticleType>
+    void apply( ExecSpace, ParticleType )
+    {
+    }
+};
+
 template <class BCIndexSpace>
 struct BoundaryCondition<BCIndexSpace, ForceValueBCTag>
 {
@@ -262,6 +279,13 @@ auto createBoundaryCondition( BCTag, ExecSpace exec_space, Particles particles,
     return BoundaryCondition<bc_index_type, BCTag>( value, bc_indices );
 }
 
+template <class MemorySpace>
+auto createBoundaryCondition( ZeroBCTag )
+{
+    using bc_index_type = BoundaryIndexSpace<MemorySpace, RegionBoundary>;
+    return BoundaryCondition<bc_index_type, ZeroBCTag>();
+}
+
 } // namespace CabanaPD
 
 #endif
diff --git a/src/CabanaPD_Solver.hpp b/src/CabanaPD_Solver.hpp
index ac718b31..de359b8c 100644
--- a/src/CabanaPD_Solver.hpp
+++ b/src/CabanaPD_Solver.hpp
@@ -90,7 +90,7 @@ class SolverBase
 };
 
 template <class MemorySpace, class InputType, class ParticleType,
-          class ForceModel>
+          class ForceModel, class BoundaryCondition>
 class SolverElastic
 {
   public:
@@ -108,12 +108,14 @@ class SolverElastic
                            Cabana::VerletLayout2D, Cabana::TeamOpTag>;
     using neigh_iter_tag = Cabana::SerialOpTag;
     using input_type = InputType;
+    using bc_type = BoundaryCondition;
 
     SolverElastic( input_type _inputs,
                    std::shared_ptr<particle_type> _particles,
-                   force_model_type force_model )
+                   force_model_type force_model, bc_type bc )
         : inputs( _inputs )
         , particles( _particles )
+        , boundary_condition( bc )
     {
         force_time = 0;
         integrate_time = 0;
@@ -193,6 +195,9 @@ class SolverElastic
         computeForce( *force, *particles, *neighbors, neigh_iter_tag{} );
         computeEnergy( *force, *particles, *neighbors, neigh_iter_tag() );
 
+        // Add boundary condition.
+        boundary_condition.apply( exec_space(), *particles );
+
         particles->output( 0, 0.0, output_reference );
         init_time += init_timer.seconds();
     }
@@ -229,6 +234,9 @@ class SolverElastic
             computeForce( *force, *particles, *neighbors, neigh_iter_tag{} );
             force_time += force_timer.seconds();
 
+            // Add boundary condition.
+            boundary_condition.apply( exec_space(), *particles );
+
             // Integrate - velocity Verlet second half.
             integrate_timer.reset();
             integrator->finalHalfStep( *particles );
@@ -320,6 +328,7 @@ class SolverElastic
     std::shared_ptr<integrator_type> integrator;
     std::shared_ptr<force_type> force;
     std::shared_ptr<neighbor_type> neighbors;
+    bc_type boundary_condition;
 
     std::string output_file;
     std::string error_file;
@@ -343,11 +352,12 @@ class SolverElastic
 template <class MemorySpace, class InputType, class ParticleType,
           class ForceModel, class BoundaryCondition, class PrenotchType>
 class SolverFracture
-    : public SolverElastic<MemorySpace, InputType, ParticleType, ForceModel>
+    : public SolverElastic<MemorySpace, InputType, ParticleType, ForceModel,
+                           BoundaryCondition>
 {
   public:
-    using base_type =
-        SolverElastic<MemorySpace, InputType, ParticleType, ForceModel>;
+    using base_type = SolverElastic<MemorySpace, InputType, ParticleType,
+                                    ForceModel, BoundaryCondition>;
     using exec_space = typename base_type::exec_space;
     using memory_space = typename base_type::memory_space;
 
@@ -366,8 +376,7 @@ class SolverFracture
                     std::shared_ptr<particle_type> _particles,
                     force_model_type force_model, bc_type bc,
                     prenotch_type prenotch )
-        : base_type( _inputs, _particles, force_model )
-        , boundary_condition( bc )
+        : base_type( _inputs, _particles, force_model, bc )
     {
         init_timer.reset();
 
@@ -475,13 +484,13 @@ class SolverFracture
     using base_type::output_reference;
 
   protected:
+    using base_type::boundary_condition;
     using base_type::comm;
     using base_type::force;
     using base_type::inputs;
     using base_type::integrator;
     using base_type::neighbors;
     using base_type::particles;
-    bc_type boundary_condition;
 
     using NeighborView = typename Kokkos::View<int**, memory_space>;
     NeighborView mu;
@@ -505,14 +514,14 @@ class SolverFracture
 };
 
 template <class MemorySpace, class InputsType, class ParticleType,
-          class ForceModel>
+          class ForceModel, class BCType>
 auto createSolverElastic( InputsType inputs,
                           std::shared_ptr<ParticleType> particles,
-                          ForceModel model )
+                          ForceModel model, BCType bc )
 {
-    return std::make_shared<
-        SolverElastic<MemorySpace, InputsType, ParticleType, ForceModel>>(
-        inputs, particles, model );
+    return std::make_shared<SolverElastic<MemorySpace, InputsType, ParticleType,
+                                          ForceModel, BCType>>(
+        inputs, particles, model, bc );
 }
 
 template <class MemorySpace, class InputsType, class ParticleType,