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

src/shared/materials/general_continuum.h

@@ -67,6 +67,17 @@ class GeneralContinuum : public WeaklyCompressibleFluid
     virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress);
 
     virtual GeneralContinuum *ThisObjectPtr() override { return this; };
+
+    class GeneralContinuumKernel
+    {
+      public:
+        GeneralContinuumKernel(GeneralContinuum &encloser):
+        E_(encloser.E_), G_(encloser.G_),K_(encloser.K_){};
+      protected:
+        Real E_;                 /* Youngs or tensile modules  */
+        Real G_;                 /* shear modules  */
+        Real K_;                 /* bulk modules  */
+    };
 };
 
 class PlasticContinuum : public GeneralContinuum
@@ -97,6 +108,67 @@ class PlasticContinuum : public GeneralContinuum
     virtual Mat3d ReturnMapping(Mat3d &stress_tensor);
 
     virtual GeneralContinuum *ThisObjectPtr() override { return this; };
+
+    class PlasticKernel:GeneralContinuum::GeneralContinuumKernel
+    {
+      public:
+
+        PlasticKernel(PlasticContinuum &encloser) : GeneralContinuum::GeneralContinuumKernel(encloser),
+        psi_(encloser.psi_),alpha_phi_(encloser.alpha_phi_),k_c_(encloser.k_c_)
+        {};
+
+        Real getDPConstantsA(Real friction_angle)
+        {
+          return tan(friction_angle) / sqrt(9.0 + 12.0 * tan(friction_angle) * tan(friction_angle));
+        };
+
+        Mat3d 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 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;
+        };
+
+
+      protected:
+          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
@@ -125,4 +197,4 @@ class J2Plasticity : public GeneralContinuum
     virtual J2Plasticity *ThisObjectPtr() override { return this; };
 };
 } // namespace SPH
-#endif // GENERAL_CONTINUUM_H
+#endif // GENERAL_CONTINUUM_H

src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h

@@ -40,4 +40,9 @@
 #include "particle_sort_ck.hpp"
 #include "simple_algorithms_ck.h"
 
+//soil mechanics
+#include "continuum_integration_1st_ck.h"
+#include "continuum_integration_1st_ck.hpp"
+#include "continuum_integration_2nd_ck.h"
+#include "continuum_integration_2nd_ck.hpp"
 #endif // ALL_SHARED_PHYSICAL_DYNAMICS_CK_H

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

@@ -0,0 +1,156 @@
+/* ------------------------------------------------------------------------- *
+ *                                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_integration_1st_ck.h
+ * @brief 	Here, we define the algorithm classes for continuum dynamics within the body.
+ * @details CK and SYCL version.
+ * @author	Shuang Li,Xiangyu Hu and Shuaihao Zhang
+ */
+#ifndef CONTINUUM_INTEGRATION_1ST_CK_H
+#define CONTINUUM_INTEGRATION_1ST_CK_H
+
+#include "constraint_dynamics.h"
+#include "acoustic_step_1st_half.h"
+#include "general_continuum.h"
+namespace SPH
+{
+namespace continuum_dynamics
+{
+
+template <class BaseInteractionType>
+class PlasticAcousticStep : public fluid_dynamics::AcousticStep<BaseInteractionType>
+{
+
+  public:
+    template <class DynamicsIdentifier>
+    explicit PlasticAcousticStep(DynamicsIdentifier &identifier);
+    virtual ~PlasticAcousticStep(){};
+
+  protected:
+    PlasticContinuum &plastic_continuum_;
+    DiscreteVariable<Mat3d> *dv_stress_tensor_3D_, *dv_strain_tensor_3D_, *dv_stress_rate_3D_, *dv_strain_rate_3D_;
+    DiscreteVariable<Matd> *dv_velocity_gradient_;
+
+};
+
+//Step1-Inner
+template <typename...>
+class PlasticAcousticStep1stHalf;
+
+template <class RiemannSolverType, class KernelCorrectionType, typename... Parameters>
+class PlasticAcousticStep1stHalf<Inner<OneLevel, RiemannSolverType, KernelCorrectionType, Parameters...>>
+    : public PlasticAcousticStep<Interaction<Inner<Parameters...>>>
+{
+    using BaseInteraction = PlasticAcousticStep<Interaction<Inner<Parameters...>>>;
+
+  public:
+    explicit PlasticAcousticStep1stHalf(Relation<Inner<Parameters...>> &inner_relation);
+    virtual ~PlasticAcousticStep1stHalf(){};
+
+    class InitializeKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        InitializeKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser);
+        void initialize(size_t index_i, Real dt = 0.0);
+
+      protected:
+        Real *rho_, *p_, *drho_dt_;
+        Vecd *vel_, *dpos_;
+        Mat3d *stress_tensor_3D_;
+    };
+
+    class InteractKernel : public BaseInteraction::InteractKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser);
+        void interact(size_t index_i, Real dt = 0.0);
+
+      protected:
+        KernelCorrectionType correction_;
+        RiemannSolverType riemann_solver_;
+        Real *Vol_, *rho_, *p_, *drho_dt_, *mass_;
+        Vecd *force_;
+        Mat3d *stress_tensor_3D_;
+    };
+
+    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:
+        Real *mass_;
+        Vecd *vel_, *force_, *force_prior_;
+    };
+
+  protected:
+    KernelCorrectionType correction_;
+    RiemannSolverType riemann_solver_;
+};
+
+//Step1-wall
+template <class RiemannSolverType, class KernelCorrectionType, typename... Parameters>
+class PlasticAcousticStep1stHalf<Contact<Wall, RiemannSolverType, KernelCorrectionType, Parameters...>>
+    : public PlasticAcousticStep<Interaction<Contact<Wall, Parameters...>>>
+{
+    using BaseInteraction = PlasticAcousticStep<Interaction<Contact<Wall, Parameters...>>>;
+
+  public:
+    explicit PlasticAcousticStep1stHalf(Relation<Contact<Parameters...>> &wall_contact_relation);
+    virtual ~PlasticAcousticStep1stHalf(){};
+
+    class InteractKernel : public BaseInteraction::InteractKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index);
+        void interact(size_t index_i, Real dt = 0.0);
+
+      protected:
+        KernelCorrectionType correction_;
+        RiemannSolverType riemann_solver_;
+        Real *Vol_, *rho_, *mass_, *p_, *drho_dt_;
+        Vecd *force_, *force_prior_;
+        Mat3d *stress_tensor_3D_;
+
+        Real *wall_Vol_;
+        Vecd *wall_acc_ave_;
+
+    };
+
+  protected:
+    KernelCorrectionType correction_;
+    RiemannSolverType riemann_solver_;
+};
+
+
+using PlasticAcousticStep1stHalfWithWallRiemannCK =
+    PlasticAcousticStep1stHalf<Inner<OneLevel, AcousticRiemannSolver, NoKernelCorrection>,
+                        Contact<Wall, AcousticRiemannSolver, NoKernelCorrection>>;
+} // namespace continuum_dynamics
+} // namespace SPH
+#endif // CONTINUUM_INTEGRATION_1ST_CK_H