ORNL · streeve · Aug 27, 2024
diff --git a/src/CabanaPD.hpp b/src/CabanaPD.hpp
@@ -29,8 +29,10 @@
 #include <CabanaPD_Types.hpp>
 #include <CabanaPD_config.hpp>
 
+#include <force/CabanaPD_ForceModels_Correspondence.hpp>
 #include <force/CabanaPD_ForceModels_LPS.hpp>
 #include <force/CabanaPD_ForceModels_PMB.hpp>
+#include <force/CabanaPD_Force_Correspondence.hpp>
 #include <force/CabanaPD_Force_LPS.hpp>
 #include <force/CabanaPD_Force_PMB.hpp>
 

diff --git a/src/CabanaPD_ForceModels.hpp b/src/CabanaPD_ForceModels.hpp
@@ -31,6 +31,29 @@ struct BaseForceModel
     void thermalStretch( double&, const int, const int ) const {}
 };
 
+struct BaseInfluenceForceModel : public BaseForceModel
+{
+    using base_type = BaseForceModel;
+
+    using base_type::delta;
+    int influence_type;
+
+    BaseInfluenceForceModel( const double _delta, const int _influence = 0 )
+        : base_type( _delta )
+        , influence_type( _influence )
+    {
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    double influence_function( const double xi ) const
+    {
+        if ( influence_type == 1 )
+            return 1.0 / xi;
+        else
+            return 1.0;
+    }
+};
+
 template <typename TemperatureType>
 struct BaseTemperatureModel
 {

diff --git a/src/CabanaPD_Particles.hpp b/src/CabanaPD_Particles.hpp
@@ -652,6 +652,155 @@ class Particles<MemorySpace, LPS, TemperatureIndependent, Dimension>
     using base_type::_timer;
 };
 
+template <class MemorySpace, int Dimension>
+class Particles<MemorySpace, Correspondence, TemperatureIndependent, Dimension>
+    : public Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>
+{
+  public:
+    using self_type = Particles<MemorySpace, Correspondence,
+                                TemperatureIndependent, Dimension>;
+    using base_type =
+        Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>;
+    using memory_space = typename base_type::memory_space;
+    using base_type::dim;
+
+    // Per particle.
+    using base_type::n_ghost;
+    using base_type::n_global;
+    using base_type::n_local;
+    using base_type::size;
+
+    // Split?
+    using tensor_type = Cabana::MemberTypes<double[dim][dim]>;
+    using aosoa_tensor_type = Cabana::AoSoA<tensor_type, memory_space, 1>;
+
+    // Per type.
+    using base_type::n_types;
+
+    // Simulation total domain.
+    using base_type::global_mesh_ext;
+
+    // Simulation sub domain (single MPI rank).
+    using base_type::ghost_mesh_hi;
+    using base_type::ghost_mesh_lo;
+    using base_type::local_mesh_ext;
+    using base_type::local_mesh_hi;
+    using base_type::local_mesh_lo;
+
+    using base_type::dx;
+    using base_type::local_grid;
+
+    using base_type::halo_width;
+
+    // Default constructor.
+    Particles()
+        : base_type()
+    {
+        _init_timer.start();
+        _aosoa_F = aosoa_tensor_type( "Particle Deformation Gradient", 0 );
+        _aosoa_Kinv = aosoa_tensor_type( "Particle Shape Tensor", 0 );
+        _init_timer.stop();
+    }
+
+    // Constructor which initializes particles on regular grid.
+    template <typename... Args>
+    Particles( Args&&... args )
+        : base_type( std::forward<Args>( args )... )
+    {
+        _init_timer.start();
+        _aosoa_F =
+            aosoa_tensor_type( "Particle Deformation Gradient", n_local );
+        _aosoa_Kinv = aosoa_tensor_type( "Particle Shape Tensor", n_local );
+        init_corr();
+        _init_timer.stop();
+    }
+
+    template <typename... Args>
+    void createParticles( Args&&... args )
+    {
+        // Forward arguments to standard or custom particle creation.
+        base_type::createParticles( std::forward<Args>( args )... );
+        _init_timer.start();
+        _aosoa_F.resize( n_local );
+        _aosoa_Kinv.resize( n_local );
+        _init_timer.stop();
+    }
+
+    auto sliceDefGradient()
+    {
+        return Cabana::slice<0>( _aosoa_F, "deformation_gradient" );
+    }
+    auto sliceDefGradient() const { return sliceDefGradient(); }
+    auto sliceShapeTensor()
+    {
+        return Cabana::slice<0>( _aosoa_Kinv, "shape_tensor" );
+    }
+    auto sliceShapeTensor() const { return sliceShapeTensor(); }
+
+    void resize( int new_local, int new_ghost )
+    {
+        base_type::resize( new_local, new_ghost );
+        _timer.start();
+        _aosoa_F.resize( new_local + new_ghost );
+        _aosoa_Kinv.resize( new_local + new_ghost );
+        _timer.stop();
+    }
+
+    void output( [[maybe_unused]] const int output_step,
+                 [[maybe_unused]] const double output_time,
+                 [[maybe_unused]] const bool use_reference = true )
+    {
+        _output_timer.start();
+
+#ifdef Cabana_ENABLE_HDF5
+        Cabana::Experimental::HDF5ParticleOutput::writeTimeStep(
+            h5_config, "particles", MPI_COMM_WORLD, output_step, output_time,
+            n_local, base_type::getPosition( use_reference ),
+            base_type::sliceStrainEnergy(), base_type::sliceForce(),
+            base_type::sliceDisplacement(), base_type::sliceVelocity(),
+            base_type::sliceDamage(), sliceDefGradient() );
+#else
+#ifdef Cabana_ENABLE_SILO
+        Cabana::Grid::Experimental::SiloParticleOutput::
+            writePartialRangeTimeStep(
+                "particles", local_grid->globalGrid(), output_step, output_time,
+                0, n_local, base_type::getPosition( use_reference ),
+                base_type::sliceStrainEnergy(), base_type::sliceForce(),
+                base_type::sliceDisplacement(), base_type::sliceVelocity(),
+                base_type::sliceDamage(), sliceDefGradient() );
+#else
+        log( std::cout, "No particle output enabled." );
+#endif
+#endif
+
+        _output_timer.stop();
+    }
+
+    friend class Comm<self_type, PMB, TemperatureIndependent>;
+    friend class Comm<self_type, LPS, TemperatureIndependent>;
+    friend class Comm<self_type, Correspondence, TemperatureIndependent>;
+
+  protected:
+    void init_corr()
+    {
+        auto F = sliceDefGradient();
+        Cabana::deep_copy( F, 0.0 );
+        auto Kinv = sliceShapeTensor();
+        Cabana::deep_copy( Kinv, 0.0 );
+    }
+
+    aosoa_tensor_type _aosoa_F;
+    aosoa_tensor_type _aosoa_Kinv;
+
+#ifdef Cabana_ENABLE_HDF5
+    using base_type::h5_config;
+#endif
+
+    using base_type::_init_timer;
+    using base_type::_output_timer;
+    using base_type::_timer;
+};
+
 template <class MemorySpace, int Dimension>
 class Particles<MemorySpace, PMB, TemperatureDependent, Dimension>
     : public Particles<MemorySpace, PMB, TemperatureIndependent, Dimension>

diff --git a/src/CabanaPD_Types.hpp b/src/CabanaPD_Types.hpp
@@ -39,6 +39,9 @@ struct LPS
 struct LinearLPS
 {
 };
+struct Correspondence
+{
+};
 
 } // namespace CabanaPD
 #endif
diff --git a/src/force/CabanaPD_ForceModels_Correspondence.hpp b/src/force/CabanaPD_ForceModels_Correspondence.hpp
@@ -0,0 +1,91 @@
+/****************************************************************************
+ * Copyright (c) 2022 by Oak Ridge National Laboratory                      *
+ * All rights reserved.                                                     *
+ *                                                                          *
+ * This file is part of CabanaPD. CabanaPD is distributed under a           *
+ * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
+ * the top-level directory.                                                 *
+ *                                                                          *
+ * SPDX-License-Identifier: BSD-3-Clause                                    *
+ ****************************************************************************/
+
+#ifndef FORCE_MODELS_CORR_H
+#define FORCE_MODELS_CORR_H
+
+#include <CabanaPD_ForceModels.hpp>
+#include <CabanaPD_Types.hpp>
+
+namespace CabanaPD
+{
+template <>
+struct ForceModel<Correspondence, Elastic> : public BaseInfluenceForceModel
+{
+    using base_type = BaseInfluenceForceModel;
+    using base_model = Correspondence;
+    using fracture_type = Elastic;
+    using thermal_type = TemperatureIndependent;
+
+    using base_type::delta;
+    using base_type::influence_type;
+
+    double K;
+    double G;
+    double theta_coeff;
+    double s_coeff;
+
+    ForceModel( const double _delta, const double _K, const double _G,
+                const int _influence = 0 )
+        : base_type( _delta, _influence )
+    {
+        set_param( _delta, _K, _G );
+    }
+
+    void set_param( const double _delta, const double _K, const double _G )
+    {
+        delta = _delta;
+        K = _K;
+        G = _G;
+
+        theta_coeff = 3.0 * K - 5.0 * G;
+        s_coeff = 15.0 * G;
+    }
+
+    template <typename DefGradType>
+    KOKKOS_INLINE_FUNCTION auto getStress( DefGradType F, const int i )
+    {
+        double epsilon[3][3];
+
+        for ( std::size_t d = 0; d < 3; d++ )
+        {
+            epsilon[d][d] = F( i, d, d ) - 1.0;
+            int d2 = getComponent( d );
+            epsilon[d][d2] = 0.5 * ( F( i, d, d2 ) + F( i, d2, d ) );
+            epsilon[d2][d] = epsilon[d][d2];
+        }
+
+        double sigma[3][3];
+        double diag = ( K - 2.0 / 3.0 * G ) *
+                      ( epsilon[0][0] + epsilon[1][1] + epsilon[2][2] );
+        for ( std::size_t d = 0; d < 3; d++ )
+        {
+            sigma[d][d] = diag + 2.0 * G * epsilon[d][d];
+            int d2 = getComponent( d );
+            sigma[d][d2] = 2.0 * G * epsilon[d][d2];
+            sigma[d2][d] = sigma[d][d2];
+        }
+        return sigma;
+    }
+
+    auto getComponent( const int d )
+    {
+        int d2;
+        if ( d < 2 )
+            return d + 1;
+        else
+            return 0;
+    }
+};
+
+} // namespace CabanaPD
+
+#endif
diff --git a/src/force/CabanaPD_ForceModels_LPS.hpp b/src/force/CabanaPD_ForceModels_LPS.hpp
@@ -18,27 +18,24 @@
 namespace CabanaPD
 {
 template <>
-struct ForceModel<LPS, Elastic> : public BaseForceModel
+struct ForceModel<LPS, Elastic> : public BaseInfluenceForceModel
 {
-    using base_type = BaseForceModel;
+    using base_type = BaseInfluenceForceModel;
     using base_model = LPS;
     using fracture_type = Elastic;
     using thermal_type = TemperatureIndependent;
 
     using base_type::delta;
-
-    int influence_type;
+    using base_type::influence_type;
 
     double K;
     double G;
     double theta_coeff;
     double s_coeff;
 
-    ForceModel(){};
     ForceModel( const double _delta, const double _K, const double _G,
                 const int _influence = 0 )
-        : base_type( _delta )
-        , influence_type( _influence )
+        : base_type( _delta, _influence )
     {
         set_param( _delta, _K, _G );
     }
@@ -52,14 +49,6 @@ struct ForceModel<LPS, Elastic> : public BaseForceModel
         theta_coeff = 3.0 * K - 5.0 * G;
         s_coeff = 15.0 * G;
     }
-
-    KOKKOS_INLINE_FUNCTION double influence_function( double xi ) const
-    {
-        if ( influence_type == 1 )
-            return 1.0 / xi;
-        else
-            return 1.0;
-    }
 };
 
 template <>
@@ -80,7 +69,6 @@ struct ForceModel<LPS, Fracture> : public ForceModel<LPS, Elastic>
     double s0;
     double bond_break_coeff;
 
-    ForceModel() {}
     ForceModel( const double _delta, const double _K, const double _G,
                 const double _G0, const int _influence = 0 )
         : base_type( _delta, _K, _G, _influence )