Successfully Running **ck_version** and **sycl_version** of Soil Mechanics

src/shared/common/base_data_type.h

@@ -56,9 +56,11 @@ namespace SPH
 #if SPHINXSYS_USE_FLOAT
 using Real = float;
 using UnsignedInt = u_int32_t;
+namespace math = sycl;
 #else
 using Real = double;
 using UnsignedInt = size_t;
+namespace math = std;
 #endif // SPHINXSYS_USE_FLOAT
 
 /** Vector with integers. */

src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h

@@ -32,6 +32,7 @@
 
 #include "base_general_dynamics.h"
 #include "general_continuum.h"
+#include "general_continuum.hpp"
 
 namespace SPH
 {

src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h

@@ -34,6 +34,7 @@
 #include "constraint_dynamics.h"
 #include "fluid_integration.hpp"
 #include "general_continuum.h"
+#include "general_continuum.hpp"
 namespace SPH
 {
 namespace continuum_dynamics

src/shared/physical_closure/materials/all_materials.h

@@ -36,5 +36,6 @@
 #include "diffusion_reaction.h"
 #include "elastic_solid.h"
 #include "general_continuum.h"
+#include "general_continuum.hpp"
 #include "inelastic_solid.h"
 #include "weakly_compressible_fluid.h"

src/shared/physical_closure/materials/general_continuum.cpp

@@ -1,4 +1,5 @@
 #include "general_continuum.h"
+#include "general_continuum.hpp"
 
 namespace SPH
 {

src/shared/physical_closure/materials/general_continuum.h

@@ -65,6 +65,30 @@ class GeneralContinuum : public WeaklyCompressibleFluid
     Real ContactStiffness() { return contact_stiffness_; };
 
     virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress);
+
+
+    class GeneralContinuumKernel
+    {
+      public:
+        GeneralContinuumKernel(GeneralContinuum &encloser):
+        E_(encloser.E_), G_(encloser.G_),K_(encloser.K_),
+        nu_(encloser.nu_),contact_stiffness_(encloser.contact_stiffness_),
+        rho0_(encloser.rho0_){};
+
+        inline Real getYoungsModulus() { return E_; };
+        inline Real getPoissonRatio() { return nu_; };
+        inline Real getDensity() { return rho0_; };
+        inline Real getBulkModulus(Real youngs_modulus, Real poisson_ratio);
+        inline Real getShearModulus(Real youngs_modulus, Real poisson_ratio);
+        inline Real getLambda(Real youngs_modulus, Real poisson_ratio);
+      protected:
+        Real E_;                 /* Youngs or tensile modules  */
+        Real G_;                 /* shear modules  */
+        Real K_;                 /* bulk modules  */
+        Real nu_;                /* Poisson ratio  */
+        Real contact_stiffness_; /* contact-force stiffness related to bulk modulus*/
+        Real rho0_; /* contact-force stiffness related to bulk modulus*/
+    };
 };
 
 class PlasticContinuum : public GeneralContinuum
@@ -93,6 +117,32 @@ class PlasticContinuum : public GeneralContinuum
 
     virtual Mat3d ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor);
     virtual Mat3d ReturnMapping(Mat3d &stress_tensor);
+
+
+    class PlasticKernel: public GeneralContinuum::GeneralContinuumKernel
+    {
+      public:
+
+        PlasticKernel(PlasticContinuum &encloser) : GeneralContinuum::GeneralContinuumKernel(encloser),
+        c_(encloser.c_),phi_(encloser.phi_),
+        psi_(encloser.psi_),alpha_phi_(encloser.alpha_phi_),k_c_(encloser.k_c_){};
+
+        inline Real getDPConstantsA(Real friction_angle);
+        inline Mat3d ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor);  
+        inline Mat3d ReturnMapping(Mat3d &stress_tensor);
+        inline Real getFrictionAngle() { return phi_; };
+
+
+      protected:
+          Real c_;                            /* cohesion  */
+          Real phi_;                          /* friction angle  */
+          Real psi_;                          /* dilatancy angle  */
+          Real alpha_phi_;                    /* Drucker-Prager's constants */
+          Real k_c_;                          /* Drucker-Prager's constants */
+          Real stress_dimension_ = 3.0; /* plain strain condition */   //Temporarily cancel const --need to check
+
+
+    };
 };
 
 class J2Plasticity : public GeneralContinuum
@@ -120,4 +170,4 @@ class J2Plasticity : public GeneralContinuum
     virtual Real HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor);
 };
 } // namespace SPH
-#endif // GENERAL_CONTINUUM_H
+#endif // GENERAL_CONTINUUM_H

src/shared/physical_closure/materials/general_continuum.hpp

@@ -0,0 +1,67 @@
+#ifndef GENERAL_CONTINUUM_HPP
+#define GENERAL_CONTINUUM_HPP
+
+#include "general_continuum.h"
+
+namespace SPH
+{
+Real GeneralContinuum::GeneralContinuumKernel::getBulkModulus(Real youngs_modulus, Real poisson_ratio)
+{
+    return youngs_modulus / 3.0 / (1.0 - 2.0 * poisson_ratio);
+}
+//=================================================================================================//
+Real GeneralContinuum::GeneralContinuumKernel::getShearModulus(Real youngs_modulus, Real poisson_ratio)
+{
+    return 0.5 * youngs_modulus / (1.0 + poisson_ratio);
+}
+//=================================================================================================//
+Real GeneralContinuum::GeneralContinuumKernel::getLambda(Real youngs_modulus, Real poisson_ratio)
+{
+    return nu_ * youngs_modulus / (1.0 + poisson_ratio) / (1.0 - 2.0 * poisson_ratio);
+}
+//=================================================================================================//
+//=================================================================================================//
+Real PlasticContinuum::PlasticKernel::getDPConstantsA(Real friction_angle)
+{
+    return tan(friction_angle) / sqrt(9.0 + 12.0 * tan(friction_angle) * tan(friction_angle));
+};
+//=================================================================================================//
+Mat3d PlasticContinuum::PlasticKernel::ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor)
+{
+    Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());
+    Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose());
+    Mat3d deviatoric_strain_rate = strain_rate - (1.0 / stress_dimension_) * strain_rate.trace() * Mat3d::Identity();
+    Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor;
+    Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity();
+    Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();
+    Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor.trace() - k_c_;
+    Real lambda_dot_ = 0;
+    Mat3d g = Mat3d::Zero();
+    if (f >= TinyReal)
+    {
+        Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum();
+        // non-associate flow rule
+        lambda_dot_ = (3.0 * alpha_phi_ * K_ * strain_rate.trace() + (G_ / (sqrt(stress_tensor_J2)+TinyReal)) * deviatoric_stress_times_strain_rate) / (9.0 * alpha_phi_ * K_ * getDPConstantsA(psi_) + G_);
+        g = lambda_dot_ * (3.0 * K_ * getDPConstantsA(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2+ TinyReal)));
+    }
+    Mat3d stress_rate_temp = stress_rate_elastic - g;
+    return stress_rate_temp;
+}; 
+//=================================================================================================//
+Mat3d PlasticContinuum::PlasticKernel::ReturnMapping(Mat3d &stress_tensor)
+{
+    Real stress_tensor_I1 = stress_tensor.trace();
+    if (-alpha_phi_ * stress_tensor_I1 + k_c_ < 0)
+    stress_tensor -= (1.0 / stress_dimension_) * (stress_tensor_I1 - k_c_ / alpha_phi_) * Mat3d::Identity();
+    stress_tensor_I1 = stress_tensor.trace();
+    Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity();
+    volatile Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();
+    if (-alpha_phi_ * stress_tensor_I1 + k_c_ < sqrt(stress_tensor_J2))
+    {
+        Real r = (-alpha_phi_ * stress_tensor_I1 + k_c_) / (sqrt(stress_tensor_J2) + TinyReal);
+        stress_tensor = r * deviatoric_stress_tensor + (1.0 / stress_dimension_) * stress_tensor_I1 * Mat3d::Identity();
+    }
+    return stress_tensor;
+}
+}// namespace SPH
+#endif //GENERAL_CONTINUUM_HPP

src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h

@@ -40,5 +40,6 @@
 #include "interaction_algorithms_ck.hpp"
 #include "particle_sort_ck.hpp"
 #include "simple_algorithms_ck.h"
+#include "all_continum_dynamics.h"
 
 #endif // ALL_SHARED_PHYSICAL_DYNAMICS_CK_H

src/shared/shared_ck/particle_dynamics/continum_dynamics/all_continum_dynamics.h

@@ -0,0 +1,39 @@
+/* ------------------------------------------------------------------------- *
+ *                                SPHinXsys                                  *
+ * ------------------------------------------------------------------------- *
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for    *
+ * physical accurate simulation and aims to model coupled industrial dynamic *
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH   *
+ * (smoothed particle hydrodynamics), a meshless computational method using  *
+ * particle discretization.                                                  *
+ *                                                                           *
+ * SPHinXsys is partially funded by German Research Foundation               *
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,            *
+ *  HU1527/12-1 and HU1527/12-4.                                             *
+ *                                                                           *
+ * Portions copyright (c) 2017-2023 Technical University of Munich and       *
+ * the authors' affiliations.                                                *
+ *                                                                           *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may   *
+ * not use this file except in compliance with the License. You may obtain a *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+/**
+ * @file    all_continum_dynamics.h
+ */
+
+
+#pragma once
+
+#include "continuum_integration_1st_ck.h"
+#include "continuum_integration_1st_ck.hpp"
+#include "continuum_integration_2nd_ck.h"
+#include "continuum_integration_2nd_ck.hpp"
+#include "continuum_dynamics_variable_ck.h"
+#include "continuum_dynamics_variable_ck.hpp"
+#include "stress_diffusion_ck.h"
+#include "stress_diffusion_ck.hpp"
+#include "initilization_dynamics_ck.h"
+#include "initilization_dynamics_ck.hpp"

src/shared/shared_ck/particle_dynamics/continum_dynamics/continuum_dynamics_variable_ck.h

@@ -0,0 +1,99 @@
+/* ------------------------------------------------------------------------- *
+ *                                SPHinXsys                                  *
+ * ------------------------------------------------------------------------- *
+ * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
+ * Hydrodynamics for industrial compleX systems. It provides C++ APIs for    *
+ * physical accurate simulation and aims to model coupled industrial dynamic *
+ * systems including fluid, solid, multi-body dynamics and beyond with SPH   *
+ * (smoothed particle hydrodynamics), a meshless computational method using  *
+ * particle discretization.                                                  *
+ *                                                                           *
+ * SPHinXsys is partially funded by German Research Foundation               *
+ * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1,            *
+ *  HU1527/12-1 and HU1527/12-4.                                             *
+ *                                                                           *
+ * Portions copyright (c) 2017-2023 Technical University of Munich and       *
+ * the authors' affiliations.                                                *
+ *                                                                           *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may   *
+ * not use this file except in compliance with the License. You may obtain a *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.        *
+ *                                                                           *
+ * ------------------------------------------------------------------------- */
+/**
+ * @file continuum_dynamics_variable_ck.h
+ * @brief Here, we define the algorithm classes for computing derived solid dynamics variables.
+ * @details These variable can be added into variable list for state output.
+ * @author Shuang Li, Xiangyu Hu and Xiangyu Hu
+ */
+
+#ifndef CONTINUUM_DYNAMICS_VARIABLE_CK_H
+#define CONTINUUM_DYNAMICS_VARIABLE_CK_H
+
+#include "base_general_dynamics.h"
+#include "general_continuum.h"
+#include "general_continuum.hpp"
+
+namespace SPH
+{
+namespace continuum_dynamics
+{
+/**
+ * @class VerticalStress
+ */
+class VerticalStressCK : public BaseDerivedVariable<Real>
+{
+  public:
+    explicit VerticalStressCK(SPHBody &sph_body);
+    virtual ~VerticalStressCK(){};
+    class UpdateKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        UpdateKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser);
+        void update(size_t index_i, Real dt = 0.0);
+
+      protected:
+        Mat3d *stress_tensor_3D_;
+        Real *derived_variable_;
+    };
+  protected:
+    DiscreteVariable<Mat3d> *dv_stress_tensor_3D_;
+    DiscreteVariable<Real> *dv_derived_variable_;
+};
+
+/**
+ * @class AccumulatedDeviatoricPlasticStrain
+ */
+class AccDeviatoricPlasticStrainCK : public BaseDerivedVariable<Real>
+{
+  using PlasticKernel = typename PlasticContinuum::PlasticKernel;
+  public:
+    explicit AccDeviatoricPlasticStrainCK(SPHBody &sph_body);
+    virtual ~AccDeviatoricPlasticStrainCK(){};
+
+    class UpdateKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        UpdateKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser);
+        void update(size_t index_i, Real dt = 0.0);
+
+      protected:
+        PlasticKernel plastic_kernel_;
+        Mat3d *stress_tensor_3D_, *strain_tensor_3D_;
+        Real *derived_variable_;
+        Real E_, nu_;
+    };
+
+  protected:
+    PlasticContinuum &plastic_continuum_;
+    DiscreteVariable<Mat3d> *dv_stress_tensor_3D_, *dv_strain_tensor_3D_;
+    DiscreteVariable<Real> *dv_derived_variable_;
+    Real E_, nu_;
+};
+
+} // namespace continuum_dynamics
+} // namespace SPH
+
+#endif // CONTINUUM_DYNAMICS_VARIABLE_CK_H

src/shared/shared_ck/particle_dynamics/continum_dynamics/continuum_dynamics_variable_ck.hpp

@@ -0,0 +1,57 @@
+#ifndef CONTINUUM_DYNAMICS_VARIABLE_CK_HPP
+#define CONTINUUM_DYNAMICS_VARIABLE_CK_HPP
+
+
+#include "continuum_dynamics_variable_ck.h"
+
+namespace SPH
+{
+namespace continuum_dynamics
+{
+    //=============================================================================================//
+    VerticalStressCK::VerticalStressCK(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "VerticalStress"),
+    dv_stress_tensor_3D_(this->particles_->template registerStateVariableOnly<Mat3d>("StressTensor3D")),
+    dv_derived_variable_(this->particles_->template registerStateVariableOnly<Real>("VerticalStress")){}
+    //=================================================================================================//
+    template <class ExecutionPolicy, class EncloserType>
+    VerticalStressCK::UpdateKernel::UpdateKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser):
+    stress_tensor_3D_(encloser.dv_stress_tensor_3D_->DelegatedData(ex_policy)),
+    derived_variable_(encloser.dv_derived_variable_->DelegatedData(ex_policy)){}
+    //=================================================================================================//   
+    void VerticalStressCK::UpdateKernel::update(size_t index_i, Real dt)
+    {
+        derived_variable_[index_i] = stress_tensor_3D_[index_i](1, 1);
+    }
+    //=============================================================================================//
+    AccDeviatoricPlasticStrainCK::AccDeviatoricPlasticStrainCK(SPHBody &sph_body)
+    : BaseDerivedVariable<Real>(sph_body, "AccDeviatoricPlasticStrain"),
+    plastic_continuum_(DynamicCast<PlasticContinuum>(this, this->sph_body_.getBaseMaterial())),
+    dv_stress_tensor_3D_(this->particles_->template registerStateVariableOnly<Mat3d>("StressTensor3D")),
+    dv_strain_tensor_3D_(this->particles_->template registerStateVariableOnly<Mat3d>("StrainTensor3D")),
+    dv_derived_variable_(this->particles_->template registerStateVariableOnly<Real>("AccDeviatoricPlasticStrain")),
+    E_(plastic_continuum_.getYoungsModulus()),nu_(plastic_continuum_.getPoissonRatio()){}
+    //=================================================================================================//
+    template <class ExecutionPolicy, class EncloserType>
+    AccDeviatoricPlasticStrainCK::UpdateKernel::UpdateKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser):
+    plastic_kernel_(encloser.plastic_continuum_),
+    stress_tensor_3D_(encloser.dv_stress_tensor_3D_->DelegatedData(ex_policy)),
+    strain_tensor_3D_(encloser.dv_strain_tensor_3D_->DelegatedData(ex_policy)),
+    derived_variable_(encloser.dv_derived_variable_->DelegatedData(ex_policy)),
+    E_(encloser.E_),nu_(encloser.nu_){}
+    //=================================================================================================//   
+    void AccDeviatoricPlasticStrainCK::UpdateKernel::update(size_t index_i, Real dt)
+    {
+        Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity();
+        Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace();
+        Mat3d elastic_strain_tensor_3D = deviatoric_stress / (2.0 * plastic_kernel_.getShearModulus(E_, nu_)) +
+                                        hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_kernel_.getBulkModulus(E_, nu_));
+        Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D;
+        Mat3d deviatoric_strain_tensor = plastic_strain_tensor_3D - (1.0 / (Real)Dimensions) * plastic_strain_tensor_3D.trace() * Mat3d::Identity();
+        Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum();
+        derived_variable_[index_i] = sqrt(sum * 2.0 / 3.0);
+    }
+} // namespace continuum_dynamics
+} // namespace SPH
+
+#endif // CONTINUUM_DYNAMICS_VARIABLE_CK_HPP