Enable multi-material

examples/mechanics/kalthoff_winkler.cpp

@@ -76,23 +76,26 @@ void kalthoffWinklerExample( const std::string filename )
     CabanaPD::Prenotch<2> prenotch( v1, v2, notch_positions );
 
     // ====================================================
-    //                    Force model
+    //                    Force model options
     // ====================================================
-    using model_type = CabanaPD::ForceModel<CabanaPD::PMB, CabanaPD::Fracture>;
-    model_type force_model( delta, K, G0 );
-    // using model_type =
-    //     CabanaPD::ForceModel<CabanaPD::LPS, CabanaPD::Fracture>;
-    // model_type force_model( delta, K, G, G0 );
+    using model_type = CabanaPD::PMB;
+    using thermal_type = CabanaPD::TemperatureIndependent;
 
     // ====================================================
     //                 Particle generation
     // ====================================================
     // Does not set displacements, velocities, etc.
     auto particles = std::make_shared<
-        CabanaPD::Particles<memory_space, typename model_type::base_model,
-                            typename model_type::thermal_type>>(
+        CabanaPD::Particles<memory_space, model_type, thermal_type>>(
         exec_space(), low_corner, high_corner, num_cells, halo_width );
 
+    std::array<double, 2> same_d = { delta, delta };
+    std::array<double, 2> fake_K = { K, K * 1.2 };
+    std::array<double, 2> fake_G0 = { G0, G0 * 1.8 };
+    auto force_model =
+        createForceModel( model_type{}, CabanaPD::Fracture{}, thermal_type{},
+                          *particles, same_d, fake_K, fake_G0 );
+
     // ====================================================
     //                Boundary conditions
     // ====================================================
@@ -111,6 +114,7 @@ void kalthoffWinklerExample( const std::string filename )
     auto x = particles->sliceReferencePosition();
     auto v = particles->sliceVelocity();
     auto f = particles->sliceForce();
+    auto type = particles->sliceType();
 
     double v0 = inputs["impactor_velocity"];
     auto init_functor = KOKKOS_LAMBDA( const int pid )
@@ -121,6 +125,8 @@ void kalthoffWinklerExample( const std::string filename )
         if ( x( pid, 1 ) > y_prenotch1 && x( pid, 1 ) < y_prenotch2 &&
              x( pid, 0 ) < -0.5 * height + dx )
             v( pid, 0 ) = v0;
+        if ( x( pid, 1 ) > 0.0 )
+            type( pid ) = 1;
     };
     particles->updateParticles( exec_space{}, init_functor );
 

examples/thermomechanics/inputs/thermal_deformation.json

@@ -1,8 +1,9 @@
 {
     "num_cells"             : {"value": [100, 30, 3]},
     "system_size"           : {"value": [1.0, 0.3, 0.03], "unit": "m"},
+    "num_types"             : {"value": 1},
     "density"               : {"value": 3980,  "unit": "kg/m^3"},
-    "elastic_modulus"       : {"value": 370e+9, "unit": "Pa"},
+    "elastic_modulus"       : {"value": [370e+9, 100e+9, 10e+9], "unit": "Pa", "note": "11, 22, 12"},
     "thermal_coefficient"   : {"value": 7.5E-6, "unit": "oC^{-1}"},
     "reference_temperature" : {"value": 20.0, "unit": "oC"},
     "horizon"               : {"value": 0.03, "unit": "m"},  

examples/thermomechanics/thermal_deformation.cpp

@@ -82,7 +82,8 @@ void thermalDeformationExample( const std::string filename )
     //                    Force model
     // ====================================================
     auto force_model = CabanaPD::createForceModel(
-        model_type{}, CabanaPD::Elastic{}, *particles, delta, K, alpha, temp0 );
+        model_type{}, CabanaPD::Elastic{}, CabanaPD::TemperatureDependent{},
+        *particles, delta, K, alpha, temp0 );
 
     // ====================================================
     //                   Create solver

src/CabanaPD_Comm.hpp

@@ -230,12 +230,13 @@ struct HaloIds
     }
 };
 
-template <class ParticleType, class ModelType, class ThermalType>
+template <class ParticleType, class ModelType, class SpeciesType,
+          class ThermalType>
 class Comm;
 
 // FIXME: extract model from ParticleType instead.
 template <class ParticleType>
-class Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
     int mpi_size = -1;
@@ -320,11 +321,12 @@ class Comm<ParticleType, PMB, TemperatureIndependent>
 };
 
 template <class ParticleType>
-class Comm<ParticleType, LPS, TemperatureIndependent>
-    : public Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, LPS, SingleSpecies, TemperatureIndependent>
+    : public Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
-    using base_type = Comm<ParticleType, PMB, TemperatureIndependent>;
+    using base_type =
+        Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>;
     using memory_space = typename base_type::memory_space;
     using halo_type = typename base_type::halo_type;
     using base_type::gather_u;
@@ -368,12 +370,52 @@ class Comm<ParticleType, LPS, TemperatureIndependent>
     }
 };
 
+template <class ParticleType, class ModelType>
+class Comm<ParticleType, ModelType, MultiSpecies, TemperatureIndependent>
+    : public Comm<ParticleType, ModelType, SingleSpecies,
+                  TemperatureIndependent>
+{
+  public:
+    using base_type =
+        Comm<ParticleType, ModelType, SingleSpecies, TemperatureIndependent>;
+    using memory_space = typename base_type::memory_space;
+    using halo_type = typename base_type::halo_type;
+    using base_type::halo;
+
+    using base_type::_init_timer;
+    using base_type::_timer;
+
+    using gather_species_type =
+        Cabana::Gather<halo_type, typename ParticleType::aosoa_species_type>;
+    std::shared_ptr<gather_species_type> gather_species;
+
+    Comm( ParticleType& particles, int max_export_guess = 100 )
+        : base_type( particles, max_export_guess )
+    {
+        _init_timer.start();
+
+        gather_species = std::make_shared<gather_species_type>(
+            *halo, particles._aosoa_species );
+
+        particles.resize( halo->numLocal(), halo->numGhost() );
+        _init_timer.stop();
+    }
+
+    void gatherSpecies()
+    {
+        _timer.start();
+        gather_species->apply();
+        _timer.stop();
+    }
+};
+
 template <class ParticleType>
-class Comm<ParticleType, PMB, TemperatureDependent>
-    : public Comm<ParticleType, PMB, TemperatureIndependent>
+class Comm<ParticleType, PMB, SingleSpecies, TemperatureDependent>
+    : public Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>
 {
   public:
-    using base_type = Comm<ParticleType, PMB, TemperatureIndependent>;
+    using base_type =
+        Comm<ParticleType, PMB, SingleSpecies, TemperatureIndependent>;
     using memory_space = typename base_type::memory_space;
     using halo_type = typename base_type::halo_type;
     using base_type::halo;

src/CabanaPD_Force.hpp

@@ -129,8 +129,7 @@ getLinearizedDistance( const PosType& x, const PosType& u, const int i,
 template <class ExecutionSpace, class ForceType>
 class Force;
 
-template <class ExecutionSpace>
-class Force<ExecutionSpace, BaseForceModel>
+class ForceBase
 {
   protected:
     bool _half_neigh;
@@ -139,7 +138,7 @@ class Force<ExecutionSpace, BaseForceModel>
     Timer _energy_timer;
 
   public:
-    Force( const bool half_neigh )
+    ForceBase( const bool half_neigh )
         : _half_neigh( half_neigh )
     {
     }

src/CabanaPD_ForceModels.hpp

@@ -16,58 +16,164 @@
 
 namespace CabanaPD
 {
-struct BaseForceModel
+
+template <typename MemorySpace>
+struct BaseForceModel;
+
+template <>
+struct BaseForceModel<SingleSpecies>
 {
+    using species_type = SingleSpecies;
     double delta;
 
     BaseForceModel(){};
     BaseForceModel( const double _delta )
         : delta( _delta ){};
 
-    BaseForceModel( const BaseForceModel& model ) { delta = model.delta; }
+    auto horizon( const int ) { return delta; }
+    auto maxHorizon() { return delta; }
+
+    // No-op for temperature.
+    KOKKOS_INLINE_FUNCTION
+    void thermalStretch( double&, const int, const int ) const {}
+};
+
+template <typename MemorySpace>
+struct BaseForceModel
+{
+    using species_type = MultiSpecies;
+
+    using memory_space = MemorySpace;
+    using view_type_1d = Kokkos::View<double*, memory_space>;
+    view_type_1d delta;
+    double max_delta;
+    std::size_t num_types;
+
+    using view_type_2d = Kokkos::View<double**, memory_space>;
+
+    BaseForceModel( const double _delta )
+        : delta( view_type_1d( "delta", 1 ) )
+        , num_types( 1 )
+    {
+        Kokkos::deep_copy( delta, _delta );
+    }
+
+    template <typename ArrayType>
+    BaseForceModel( const ArrayType& _delta )
+        : delta( view_type_1d( "delta", _delta.size() ) )
+        , num_types( _delta.size() )
+    {
+        setParameters( _delta );
+    }
+
+    template <typename ArrayType>
+    void setParameters( const ArrayType& _delta )
+    {
+        max_delta = 0;
+        auto delta_copy = delta;
+        auto init_func = KOKKOS_LAMBDA( const int i, double& max )
+        {
+            delta_copy( i ) = _delta[i];
+            if ( delta_copy( i ) > max )
+                max = delta_copy( i );
+        };
+        using exec_space = typename memory_space::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, num_types );
+        Kokkos::parallel_reduce( "CabanaPD::Model::Init", policy, init_func,
+                                 max_delta );
+    }
+
+    auto horizon( const int i ) { return delta( i ); }
+    auto maxHorizon() { return max_delta; }
 
     // No-op for temperature.
     KOKKOS_INLINE_FUNCTION
     void thermalStretch( double&, const int, const int ) const {}
 };
 
+template <typename TemperatureType, typename SpeciesType>
+struct BaseTemperatureModel;
+
 template <typename TemperatureType>
-struct BaseTemperatureModel
+struct BaseTemperatureModel<TemperatureType, SingleSpecies>
 {
+    using species_type = SingleSpecies;
+    using memory_space = typename TemperatureType::memory_space;
     double alpha;
     double temp0;
 
     // Temperature field
     TemperatureType temperature;
 
-    BaseTemperatureModel(){};
     BaseTemperatureModel( const TemperatureType _temp, const double _alpha,
                           const double _temp0 )
         : alpha( _alpha )
         , temp0( _temp0 )
         , temperature( _temp ){};
 
-    BaseTemperatureModel( const BaseTemperatureModel& model )
+    void update( const TemperatureType _temp ) { temperature = _temp; }
+
+    // Update stretch with temperature effects.
+    KOKKOS_INLINE_FUNCTION
+    void thermalStretch( double& s, const int i, const int j ) const
     {
-        alpha = model.alpha;
-        temp0 = model.temp0;
-        temperature = model.temperature;
+        double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0;
+        s -= alpha * temp_avg;
     }
+};
 
-    void set_param( const double _alpha, const double _temp0 )
+template <typename TemperatureType, typename ParticleType>
+struct BaseTemperatureModel
+{
+    using species_type = MultiSpecies;
+    using memory_space = typename TemperatureType::memory_space;
+    using view_type_1d = Kokkos::View<double*, memory_space>;
+    view_type_1d alpha;
+    view_type_1d temp0;
+
+    // Particle fields
+    TemperatureType temperature;
+    ParticleType type;
+
+    template <typename ArrayType>
+    BaseTemperatureModel( const TemperatureType& _temp, const ArrayType& _alpha,
+                          const ArrayType& _temp0, const ParticleType& _type )
+        : alpha( view_type_1d( "delta", _alpha.size() ) )
+        , temp0( view_type_1d( "delta", _temp0.size() ) )
+        , temperature( _temp )
+        , type( _type )
     {
-        alpha = _alpha;
-        temp0 = _temp0;
+        setParameters( _alpha, _temp0 );
     }
 
-    void update( const TemperatureType _temp ) { temperature = _temp; }
+    template <typename ArrayType>
+    void setParameters( const ArrayType& _alpha, const ArrayType& _temp0 )
+    {
+        auto alpha_copy = alpha;
+        auto temp0_copy = temp0;
+        auto init_func = KOKKOS_LAMBDA( const int i )
+        {
+            alpha_copy( i ) = _alpha[i];
+            temp0_copy( i ) = _temp0[i];
+        };
+        using exec_space = typename memory_space::execution_space;
+        Kokkos::RangePolicy<exec_space> policy( 0, alpha.size() );
+        Kokkos::parallel_for( "CabanaPD::Model::Init", policy, init_func );
+    }
+
+    void update( const TemperatureType& _temp, const ParticleType& _type )
+    {
+        temperature = _temp;
+        type = _type;
+    }
 
     // Update stretch with temperature effects.
     KOKKOS_INLINE_FUNCTION
     void thermalStretch( double& s, const int i, const int j ) const
     {
-        double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0;
-        s -= alpha * temp_avg;
+        double temp_avg =
+            0.5 * ( temperature( i ) + temperature( j ) ) - temp0( type( i ) );
+        s -= alpha( type( i ) ) * temp_avg;
     }
 };