Add anisotropic material mechanics

examples/mechanics/dogbone_tensile_test.cpp

@@ -36,14 +36,19 @@ void dogboneTensileTestExample( const std::string filename )
     //                Material parameters
     // ====================================================
     double rho0 = inputs["density"];
-    double E = inputs["elastic_modulus"];
-    double nu = 0.25; // Use bond-based model
-    double K = E / ( 3 * ( 1 - 2 * nu ) );
-    double G0 = inputs["fracture_energy"];
-    double sigma_y = inputs["yield_stress"];
     double delta = inputs["horizon"];
     delta += 1e-10;
 
+    double C11 = inputs["C11"];
+    double C13 = inputs["C13"];
+    double C33 = inputs["C33"];
+
+    // Require two values for both yield and failure.
+    std::array<double, 2> sigma_y_input = inputs["yield_stress"];
+    auto sigma_y = CabanaPD::copy( sigma_y_input );
+    std::array<double, 2> sc_input = inputs["critical_stretch"];
+    auto sc = CabanaPD::copy( sc_input );
+
     // ====================================================
     //                  Discretization
     // ====================================================
@@ -154,8 +159,9 @@ void dogboneTensileTestExample( const std::string filename )
     // ====================================================
     //                    Force model
     // ====================================================
-    CabanaPD::ForceModel force_model( model_type{}, mechanics_type{},
-                                      memory_space{}, delta, K, G0, sigma_y );
+    CabanaPD::ForceModel force_model(
+        model_type{}, mechanics_type{}, CabanaPD::TransverselyIsotropic{},
+        memory_space{}, delta, C11, C13, C33, sigma_y, sc );
 
     // ====================================================
     //                   Create solver

examples/mechanics/elastic_wave.cpp

@@ -36,8 +36,10 @@ void elasticWaveExample( const std::string filename )
     //                Material parameters
     // ====================================================
     double rho0 = inputs["density"];
-    auto K = inputs["bulk_modulus"];
-    double G = inputs["shear_modulus"];
+    double C11 = inputs["C11"];
+    // double C12 = inputs["C12"];
+    double C13 = inputs["C13"];
+    double C33 = inputs["C33"];
     double delta = inputs["horizon"];
     delta += 1e-10;
 
@@ -54,9 +56,12 @@ void elasticWaveExample( const std::string filename )
     // ====================================================
     //                    Force model
     // ====================================================
-    using model_type = CabanaPD::LinearLPS;
+    using model_type = CabanaPD::LinearPMB;
+    // CabanaPD::ForceModel force_model( model_type{}, CabanaPD::NoFracture{},
+    //                                  CabanaPD::Cubic{}, delta, C11, C12 );
     CabanaPD::ForceModel force_model( model_type{}, CabanaPD::NoFracture{},
-                                      delta, K, G );
+                                      CabanaPD::TransverselyIsotropic{}, delta,
+                                      C11, C13, C33 );
 
     // ====================================================
     //                 Particle generation

examples/mechanics/inputs/dogbone_tensile_test.json

@@ -7,9 +7,12 @@
     "distance_between_grips" : {"value": 0.115, "unit": "m"},
     "density"                : {"value": 1200, "unit": "kg/m^3"},
     "elastic_modulus"        : {"value": 3e+9, "unit": "Pa"},
-    "fracture_energy"        : {"value": 3000, "unit": "J/m^2"},
+    "C11"        : {"value": 3e+9, "unit": "Pa"},
+    "C13"        : {"value": 3e+7, "unit": "Pa"},
+    "C33"        : {"value": 3e+8, "unit": "Pa"},
+    "critical_stretch"        : {"value": [0.06, 0.01], "unit": "J/m^2"},
+    "yield_stress"           : {"value": [30e6, 1e6], "unit": "Pa"},
     "horizon"                : {"value": 0.002, "unit": "m"},   
-    "yield_stress"           : {"value": 30e6, "unit": "Pa"},
     "grip_velocity"          : {"value": 1e-1, "unit": "m/s"},
     "final_time"             : {"value": 1.5e-2, "unit": "s"},
     "timestep"               : {"value": 3.0e-7, "unit": "s"},

examples/mechanics/inputs/elastic_wave.json

@@ -1,10 +1,14 @@
 {
-    "num_cells"              : {"value": [41, 41, 41]},
+    "num_cells"              : {"value": [100, 100, 100]},
     "system_size"            : {"value": [1.0, 1.0, 1.0], "unit": ""},
-    "density"                : {"value": 100.0,  "unit": ""},
-    "bulk_modulus"           : {"value": 1.0, "unit": ""},
-    "shear_modulus"          : {"value": 0.5, "unit": ""},
-    "horizon"                : {"value": 0.075, "unit": ""},  
+    "density"                : {"value": 3600.0,  "unit": ""},
+    "C11"                    : {"value": 320.0, "unit": "GPa"},
+    "C12"                    : {"value": 150.0, "unit": "GPa"},
+    "C13"                    : {"value": 150.0, "unit": "GPa"},
+    "C33"                    : {"value": 150.0, "unit": "GPa"},
+    "elastic_modulus"        : {"value": 191e+9, "unit": "Pa"},
+    "fracture_energy"        : {"value": 42408, "unit": "J/m^2"},
+    "horizon"                : {"value": 0.03, "unit": ""},  
     "final_time"             : {"value": 0.6, "unit": ""},
     "timestep"               : {"value": 0.01,  "unit": ""},
     "timestep_safety_factor" : {"value": 0.85},

src/CabanaPD_ForceModels.hpp

@@ -72,13 +72,17 @@ struct BaseNoFractureModel
     // model does not support failure.
     KOKKOS_INLINE_FUNCTION
     bool operator()( CriticalStretchTag, const int, const int, const double,
-                     const double ) const
+                     const double, const double ) const
     {
         return false;
     }
 };
 
-struct BaseFractureModel
+template <typename AnisotropyType>
+struct BaseFractureModel;
+
+template <>
+struct BaseFractureModel<Isotropic>
 {
     using fracture_type = Fracture;
 
@@ -118,8 +122,32 @@ struct BaseFractureModel
 
     KOKKOS_INLINE_FUNCTION
     bool operator()( CriticalStretchTag, const int, const int, const double r,
-                     const double xi ) const
+                     const double xi, const double ) const
+    {
+        return r * r >= bond_break_coeff * xi * xi;
+    }
+};
+
+template <>
+struct BaseFractureModel<TransverselyIsotropic>
+{
+    using fracture_type = Fracture;
+
+    Kokkos::Array<double, 2> s0;
+
+    // Constructor to work with plasticity.
+    BaseFractureModel( const Kokkos::Array<double, 2> _s0 )
+        : s0( _s0 )
+    {
+    }
+
+    KOKKOS_INLINE_FUNCTION
+    bool operator()( CriticalStretchTag, const int, const int, const double r,
+                     const double xi, const double xi_z ) const
     {
+        auto s0_o = s0[0] + ( s0[1] - s0[0] ) * xi_z * xi_z / ( xi * xi );
+        auto bond_break_coeff = ( 1.0 + s0_o ) * ( 1.0 + s0_o );
+
         return r * r >= bond_break_coeff * xi * xi;
     }
 };
@@ -187,12 +215,15 @@ struct BaseTemperatureModel<TemperatureDependent, TemperatureType>
     }
 };
 
+template <typename AnisotropyType, typename TemperatureType>
+struct ThermalFractureModel;
+
 template <typename TemperatureType>
-struct ThermalFractureModel
-    : public BaseFractureModel,
+struct ThermalFractureModel<Isotropic, TemperatureType>
+    : public BaseFractureModel<Isotropic>,
       BaseTemperatureModel<TemperatureDependent, TemperatureType>
 {
-    using base_fracture_type = BaseFractureModel;
+    using base_fracture_type = BaseFractureModel<Isotropic>;
     using base_temperature_type =
         BaseTemperatureModel<TemperatureDependent, TemperatureType>;
     using typename base_fracture_type::fracture_type;
@@ -216,7 +247,7 @@ struct ThermalFractureModel
 
     KOKKOS_INLINE_FUNCTION
     bool operator()( CriticalStretchTag, const int i, const int j,
-                     const double r, const double xi ) const
+                     const double r, const double xi, const double ) const
     {
         double temp_avg = 0.5 * ( temperature( i ) + temperature( j ) ) - temp0;
         double bond_break_coeff =
@@ -261,7 +292,7 @@ struct BaseDynamicTemperatureModel
 };
 
 template <typename PeridynamicsModelType, typename MechanicsModelType = Elastic,
-          typename DamageType = Fracture,
+          typename AnisotropyType = Isotropic, typename DamageType = Fracture,
           typename ThermalType = TemperatureIndependent, typename... DataTypes>
 struct ForceModel;
 

src/CabanaPD_Input.hpp

@@ -387,6 +387,15 @@ class Inputs
     nlohmann::json inputs;
 };
 
+template <typename ValueType, std::size_t Dim>
+auto copy( std::array<ValueType, Dim> std )
+{
+    Kokkos::Array<ValueType, Dim> kokkos;
+    for ( std::size_t i = 0; i < Dim; i++ )
+        kokkos[i] = std[i];
+    return kokkos;
+}
+
 } // namespace CabanaPD
 
 #endif

src/CabanaPD_Types.hpp

@@ -53,6 +53,18 @@ struct ElasticPerfectlyPlastic
 {
     using base_type = Plastic;
 };
+//! Fully isotropic.
+struct Isotropic
+{
+};
+//! Cubic anisotropy.
+struct Cubic
+{
+};
+//! Transversely Isotropic anisotropy.
+struct TransverselyIsotropic
+{
+};
 
 // Model category tags.
 struct Pair

src/force/CabanaPD_LPS.hpp

@@ -434,10 +434,12 @@ class Force<MemorySpace, ModelType, LPS, Fracture>
                 // Get the reference positions and displacements.
                 double xi, r, s;
                 double rx, ry, rz;
-                getDistance( x, u, i, j, xi, r, s, rx, ry, rz );
+                double xi_x, xi_y, xi_z;
+                getDistance( x, u, i, j, xi, r, s, rx, ry, rz, xi_x, xi_y,
+                             xi_z );
 
                 // Break if beyond critical stretch unless in no-fail zone.
-                if ( model( CriticalStretchTag{}, i, j, r, xi ) &&
+                if ( model( CriticalStretchTag{}, i, j, r, xi, xi_z ) &&
                      !nofail( i ) && !nofail( i ) )
                 {
                     mu( i, n ) = 0;

src/force/CabanaPD_PMB.hpp

@@ -113,11 +113,13 @@ class Force<MemorySpace, ModelType, PMB, NoFracture>
 
             double xi, r, s;
             double rx, ry, rz;
-            getDistance( x, u, i, j, xi, r, s, rx, ry, rz );
+            double xi_x, xi_y, xi_z;
+            getDistance( x, u, i, j, xi, r, s, rx, ry, rz, xi_x, xi_y, xi_z );
 
             s = model( ThermalStretchTag{}, i, j, s );
 
-            const double coeff = model( ForceCoeffTag{}, i, j, s, vol( j ) );
+            const double coeff = model( ForceCoeffTag{}, i, j, s, vol( j ), xi,
+                                        xi_x, xi_y, xi_z );
             fx_i = coeff * rx / r;
             fy_i = coeff * ry / r;
             fz_i = coeff * rz / r;
@@ -187,8 +189,8 @@ class Force<MemorySpace, ModelType, PMB, NoFracture>
 
             s = model( ThermalStretchTag{}, i, j, s );
 
-            const double coeff =
-                0.5 * model( ForceCoeffTag{}, i, j, s, vol( j ) );
+            const double coeff = 0.5 * model( ForceCoeffTag{}, i, j, s,
+                                              vol( j ), xi, xi_x, xi_y, xi_z );
             const double fx_i = coeff * rx / r;
             const double fy_i = coeff * ry / r;
             const double fz_i = coeff * rz / r;
@@ -273,12 +275,14 @@ class Force<MemorySpace, ModelType, PMB, Fracture>
                 // Get the reference positions and displacements.
                 double xi, r, s;
                 double rx, ry, rz;
-                getDistance( x, u, i, j, xi, r, s, rx, ry, rz );
+                double xi_x, xi_y, xi_z;
+                getDistance( x, u, i, j, xi, r, s, rx, ry, rz, xi_x, xi_y,
+                             xi_z );
 
                 s = model( ThermalStretchTag{}, i, j, s );
 
                 // Break if beyond critical stretch unless in no-fail zone.
-                if ( model( CriticalStretchTag{}, i, j, r, xi ) &&
+                if ( model( CriticalStretchTag{}, i, j, r, xi, xi_z ) &&
                      !nofail( i ) && !nofail( j ) )
                 {
                     mu( i, n ) = 0;
@@ -287,7 +291,8 @@ class Force<MemorySpace, ModelType, PMB, Fracture>
                 else if ( mu( i, n ) > 0 )
                 {
                     const double coeff =
-                        model( ForceCoeffTag{}, i, j, s, vol( j ), n );
+                        model( ForceCoeffTag{}, i, j, s, vol( j ), xi, xi_x,
+                               xi_y, xi_z, n );
 
                     double muij = mu( i, n );
                     fx_i = muij * coeff * rx / r;
@@ -394,7 +399,8 @@ class Force<MemorySpace, ModelType, PMB, Fracture>
                 s = model( ThermalStretchTag{}, i, j, s );
 
                 const double coeff =
-                    0.5 * model( ForceCoeffTag{}, i, j, s, vol( j ), n );
+                    0.5 * model( ForceCoeffTag{}, i, j, s, vol( j ), xi, xi_x,
+                                 xi_y, xi_z, n );
                 const double muij = mu( i, n );
                 const double fx_i = muij * coeff * rx / r;
                 const double fy_i = muij * coeff * ry / r;
@@ -471,8 +477,8 @@ class Force<MemorySpace, ModelType, LinearPMB, NoFracture>
 
             linear_s = model( ThermalStretchTag{}, i, j, linear_s );
 
-            const double coeff =
-                model( ForceCoeffTag{}, i, j, linear_s, vol( j ) );
+            const double coeff = model( ForceCoeffTag{}, i, j, linear_s,
+                                        vol( j ), xi, xi_x, xi_y, xi_z );
             fx_i = coeff * xi_x / xi;
             fy_i = coeff * xi_y / xi;
             fz_i = coeff * xi_z / xi;
@@ -539,8 +545,8 @@ class Force<MemorySpace, ModelType, LinearPMB, NoFracture>
 
             linear_s = model( ThermalStretchTag{}, i, j, linear_s );
 
-            const double coeff =
-                0.5 * model( ForceCoeffTag{}, i, j, linear_s, vol( j ) );
+            const double coeff = 0.5 * model( ForceCoeffTag{}, i, j, linear_s,
+                                              vol( j ), xi, xi_x, xi_y, xi_z );
             const double fx_i = coeff * xi_x / xi;
             const double fy_i = coeff * xi_y / xi;
             const double fz_i = coeff * xi_z / xi;