Xiangyu/sycl loop ranges

src/for_2D_build/common/data_type.h

@@ -61,8 +61,11 @@ const Matd reduced_unit_matrix{
 
 /** initial local normal, only works for thin structure dynamics. */
 const Vecd local_pseudo_n_0 = Vecd(0.0, 1.0);
-
 const Vecd ZeroVecd = Vec2d::Zero();
+
+inline Vecd degradeToVecd(const Vec3d &input) { return Vecd(input[0], input[1]); };
+inline Matd degradeToMatd(const Mat3d &input) { return input.block<2, 2>(0, 0); };
+
 } // namespace SPH
 
 #endif // DATA_TYPE_2D_H

src/for_3D_build/common/data_type.h

@@ -61,7 +61,9 @@ const Matd reduced_unit_matrix{
 /** initial local normal, only works for thin structure dynamics. */
 const Vecd local_pseudo_n_0 = Vecd(0.0, 0.0, 1.0);
 const Vecd local_pseudo_b_n_0 = Vecd(0.0, 1.0, 0.0);
-
 const Vecd ZeroVecd = Vec3d::Zero();
+
+inline Vecd degradeToVecd(const Vec3d &input) { return input; };
+inline Matd degradeToMatd(const Mat3d &input) { return input; };
 } // namespace SPH
 #endif // DATA_TYPE_3D_H

src/for_3D_build/common/scalar_functions_3d.cpp

@@ -7,4 +7,5 @@ int NextAxis(int first_axis)
 {
     return first_axis == 2 ? 0 : first_axis + 1;
 }
+//=================================================================================================//
 } // namespace SPH

src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h

@@ -260,7 +260,7 @@ class BarStressRelaxationFirstHalf : public BaseBarRelaxation
 
 /**
  * @class BarStressRelaxationSecondHalf
- * @brief computing stress relaxation process by verlet time stepping
+ * @brief computing stress relaxation process by Verlet time stepping
  * This is the second step
  */
 class BarStressRelaxationSecondHalf : public BaseBarRelaxation

src/shared/bodies/base_body.cpp

@@ -80,11 +80,6 @@ BoundingBox SPHBody::getSPHBodyBounds()
     return is_bound_set_ ? bound_ : initial_shape_->getBounds();
 }
 //=================================================================================================//
-void SPHBody::registerComputingKernel(execution::Implementation<Base> *implementation)
-{
-    all_simple_reduce_computing_kernels_.push_back(implementation);
-}
-//=================================================================================================//
 void SPHBody::defineAdaptationRatios(Real h_spacing_ratio, Real new_system_refinement_ratio)
 {
     sph_adaptation_->resetAdaptationRatios(h_spacing_ratio, new_system_refinement_ratio);

src/shared/bodies/base_body.h

@@ -38,11 +38,11 @@
 #include "adaptation.h"
 #include "all_geometries.h"
 #include "base_data_package.h"
+#include "base_implementation.h"
 #include "base_material.h"
 #include "base_particle_generator.h"
 #include "base_particles.h"
 #include "cell_linked_list.h"
-#include "execution.h"
 #include "sph_system.h"
 #include "sphinxsys_containers.h"
 
@@ -76,9 +76,7 @@ class SPHBody
     bool is_bound_set_;             /**< whether the bounding box is set */
     BoundingBox bound_;             /**< bounding box of the body */
     Shape *initial_shape_;          /**< initial volumetric geometry enclosing the body */
-    int total_body_parts_;
-    StdVec<execution::Implementation<Base> *> all_simple_reduce_computing_kernels_;
-    /**< total number of body parts */
+    int total_body_parts_;          /**< total number of body parts */
 
   public:
     SPHAdaptation *sph_adaptation_;        /**< numerical adaptation policy */
@@ -109,7 +107,6 @@ class SPHBody
     void setSPHBodyBounds(const BoundingBox &bound);
     BoundingBox getSPHBodyBounds();
     BoundingBox getSPHSystemBounds();
-    void registerComputingKernel(execution::Implementation<Base> *implementation);
     int getNewBodyPartID();
     int getTotalBodyParts() { return total_body_parts_; };
     //----------------------------------------------------------------------

src/shared/bodies/base_body_part.cpp

@@ -4,12 +4,28 @@
 namespace SPH
 {
 //=================================================================================================//
+BodyPart::BodyPart(SPHBody &sph_body, const std::string &body_part_name)
+    : sph_body_(sph_body), part_id_(sph_body.getNewBodyPartID()),
+      body_part_name_(body_part_name),
+      base_particles_(sph_body.getBaseParticles()),
+      dv_index_list_(nullptr), sv_range_size_(nullptr),
+      pos_(base_particles_.getVariableDataByName<Vecd>("Position")) {}
+//=================================================================================================//
+BodyPartByParticle::BodyPartByParticle(SPHBody &sph_body, const std::string &body_part_name)
+    : BodyPart(sph_body, body_part_name),
+      body_part_bounds_(Vecd::Zero(), Vecd::Zero()), body_part_bounds_set_(false) {}
+//=================================================================================================//
 void BodyPartByParticle::tagParticles(TaggingParticleMethod &tagging_particle_method)
 {
     for (size_t i = 0; i < base_particles_.TotalRealParticles(); ++i)
     {
         tagging_particle_method(i);
     }
+    dv_index_list_ = base_particles_.addUniqueDiscreteVariableOnly<UnsignedInt>(
+        body_part_name_, body_part_particles_.size(), [&](size_t i) -> Real
+        { return body_part_particles_[i]; });
+    sv_range_size_ = base_particles_.addUniqueSingularVariableOnly<UnsignedInt>(
+        body_part_name_ + "_Size", body_part_particles_.size());
 };
 //=============================================================================================//
 size_t BodyPartByCell::SizeOfLoopRange()
@@ -41,7 +57,6 @@ BodyRegionByParticle::BodyRegionByParticle(SPHBody &sph_body, SharedPtr<Shape> s
 {
     shape_ptr_keeper_.assignRef(shape_ptr);
 }
-//==
 //=================================================================================================//
 void BodyRegionByParticle::tagByContain(size_t particle_index)
 {

src/shared/bodies/base_body_part.h

@@ -44,23 +44,22 @@ using namespace std::placeholders;
 class BodyPart
 {
   public:
-    BodyPart(SPHBody &sph_body, const std::string &body_part_name)
-        : sph_body_(sph_body), part_id_(sph_body.getNewBodyPartID()),
-          body_part_name_(body_part_name),
-          base_particles_(sph_body.getBaseParticles()),
-          pos_(base_particles_.getVariableDataByName<Vecd>("Position")){};
+    BodyPart(SPHBody &sph_body, const std::string &body_part_name);
     virtual ~BodyPart(){};
-
     SPHBody &getSPHBody() { return sph_body_; };
     SPHSystem &getSPHSystem() { return sph_body_.getSPHSystem(); };
     std::string getName() { return body_part_name_; };
     int getPartID() { return part_id_; };
+    DiscreteVariable<UnsignedInt> *dvIndexList() { return dv_index_list_; };
+    SingularVariable<UnsignedInt> *svRangeSize() { return sv_range_size_; };
 
   protected:
     SPHBody &sph_body_;
     int part_id_;
     std::string body_part_name_;
     BaseParticles &base_particles_;
+    DiscreteVariable<UnsignedInt> *dv_index_list_;
+    SingularVariable<UnsignedInt> *sv_range_size_;
     Vecd *pos_;
 };
 
@@ -75,10 +74,7 @@ class BodyPartByParticle : public BodyPart
     BaseParticles &getBaseParticles() { return base_particles_; };
     IndexVector &LoopRange() { return body_part_particles_; };
     size_t SizeOfLoopRange() { return body_part_particles_.size(); };
-
-    BodyPartByParticle(SPHBody &sph_body, const std::string &body_part_name)
-        : BodyPart(sph_body, body_part_name),
-          body_part_bounds_(Vecd::Zero(), Vecd::Zero()), body_part_bounds_set_(false){};
+    BodyPartByParticle(SPHBody &sph_body, const std::string &body_part_name);
     virtual ~BodyPartByParticle(){};
 
     void setBodyPartBounds(BoundingBox bbox)

src/shared/body_relations/base_body_relation.h

@@ -34,7 +34,7 @@
 #include "base_geometry.h"
 #include "base_particles.h"
 #include "cell_linked_list.h"
-#include "execution.h"
+#include "base_implementation.h"
 #include "neighborhood.h"
 
 namespace SPH

src/shared/body_relations/contact_body_relation.cpp

@@ -3,6 +3,8 @@
 #include "base_particle_dynamics.h"
 #include "cell_linked_list.hpp"
 
+#include <numeric>
+
 namespace SPH
 {
 //=================================================================================================//

src/shared/common/base_data_package.h

@@ -72,38 +72,6 @@ using DataContainerAddressAssemble = DataAssemble<DataContainerAddressKeeper, Co
 template <template <typename> typename ContainerType>
 using DataContainerUniquePtrAssemble = DataAssemble<DataContainerUniquePtrKeeper, ContainerType>;
 
-/** a type irrelevant operation on the data assembles  */
-template <template <typename> typename OperationType>
-class DataAssembleOperation
-{
-    OperationType<int> integer_operation;
-    OperationType<Real> scalar_operation;
-    OperationType<Vec2d> vector2d_operation;
-    OperationType<Mat2d> matrix2d_operation;
-    OperationType<Vec3d> vector3d_operation;
-    OperationType<Mat3d> matrix3d_operation;
-
-  public:
-    template <typename... Args>
-    DataAssembleOperation(Args &&...args)
-        : integer_operation(std::forward<Args>(args)...),
-          scalar_operation(std::forward<Args>(args)...),
-          vector2d_operation(std::forward<Args>(args)...),
-          matrix2d_operation(std::forward<Args>(args)...),
-          vector3d_operation(std::forward<Args>(args)...),
-          matrix3d_operation(std::forward<Args>(args)...){};
-    template <typename... OperationArgs>
-    void operator()(OperationArgs &&...operation_args)
-    {
-        integer_operation(std::forward<OperationArgs>(operation_args)...);
-        scalar_operation(std::forward<OperationArgs>(operation_args)...);
-        vector2d_operation(std::forward<OperationArgs>(operation_args)...);
-        matrix2d_operation(std::forward<OperationArgs>(operation_args)...);
-        vector3d_operation(std::forward<OperationArgs>(operation_args)...);
-        matrix3d_operation(std::forward<OperationArgs>(operation_args)...);
-    }
-};
-
 // Please refer: https://www.cppstories.com/2022/tuple-iteration-basics/ for the following code
 template <typename DataAssembleType, typename OperationType>
 class OperationOnDataAssemble

src/shared/common/scalar_functions.h

@@ -147,6 +147,11 @@ inline bool Not_a_number(T a)
     return (std::isnan(a) || !(std::isfinite(a))) ? true : false;
 }
 
+inline Real harmonic_average(const Real &a, const Real &b)
+{
+    return 2.0 * a * b / (a + b);
+}
+
 inline Real rand_normal(Real u, Real std)
 {
     unsigned seed = (unsigned)std::chrono::system_clock::now().time_since_epoch().count();

src/shared/common/sphinxsys_containers.h

@@ -61,6 +61,8 @@ class Wall;            /**< Interaction with wall boundary */
 class Extended;        /**< An extened method of an interaction type */
 class SpatialTemporal; /**< A interaction considering spatial temporal correlations */
 class Dynamic;         /**< A dynamic interaction */
+class Constant;        /**< A constant parameter */
+class Variable;        /**< A variable parameter */
 
 using MaterialVector = StdVec<BaseMaterial *>;
 using SPHBodyVector = StdVec<SPHBody *>;

src/shared/common/sphinxsys_variable.h

@@ -77,10 +77,8 @@ class SingularVariable : public Entity
     ~SingularVariable() { delete value_; };
 
     DataType *ValueAddress() { return delegated_; };
-
     void setValue(const DataType &value) { *delegated_ = value; };
     DataType getValue() { return *delegated_; };
-
     void incrementValue(const DataType &value) { *delegated_ += value; };
 
     template <class ExecutionPolicy>

src/shared/common/vector_functions.cpp

@@ -13,66 +13,6 @@ Vec3d FirstAxisVector(const Vec3d &zero_vector)
     return Vec3d(1.0, 0.0, 0.0);
 };
 //=================================================================================================//
-Vec3d upgradeToVec3d(const Real &input)
-{
-    return Vec3d(input, 0.0, 0.0);
-}
-//=================================================================================================//
-Vec3d upgradeToVec3d(const Vec2d &input)
-{
-    return Vec3d(input[0], input[1], 0.0);
-}
-//=================================================================================================//
-Vec3d upgradeToVec3d(const Vec3d &input)
-{
-    return input;
-}
-//=================================================================================================//
-Mat3d upgradeToMat3d(const Mat2d &input)
-{
-    Mat3d output = Mat3d::Zero();
-    output.block<2, 2>(0, 0) = input;
-    return output;
-}
-//=================================================================================================//
-Mat3d upgradeToMat3d(const Mat3d &input)
-{
-    return input;
-}
-//=================================================================================================//
-Mat2d getInverse(const Mat2d &A)
-{
-    Mat2d minv = Mat2d::Zero();
-    Real det = A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
-    Real invdet = 1.0 / det;
-    minv(0, 0) = A(1, 1) * invdet;
-    minv(0, 1) = -A(0, 1) * invdet;
-    minv(1, 0) = -A(1, 0) * invdet;
-    minv(1, 1) = A(0, 0) * invdet;
-    return minv;
-}
-//=================================================================================================//
-Mat3d getInverse(const Mat3d &A)
-{
-    Real det = A(0, 0) * (A(1, 1) * A(2, 2) - A(2, 1) * A(1, 2)) -
-               A(0, 1) * (A(1, 0) * A(2, 2) - A(1, 2) * A(2, 0)) +
-               A(0, 2) * (A(1, 0) * A(2, 1) - A(1, 1) * A(2, 0));
-
-    Real invdet = 1 / det;
-    Mat3d minv = Mat3d::Zero();
-    minv(0, 0) = (A(1, 1) * A(2, 2) - A(2, 1) * A(1, 2)) * invdet;
-    minv(0, 1) = (A(0, 2) * A(2, 1) - A(0, 1) * A(2, 2)) * invdet;
-    minv(0, 2) = (A(0, 1) * A(1, 2) - A(0, 2) * A(1, 1)) * invdet;
-    minv(1, 0) = (A(1, 2) * A(2, 0) - A(1, 0) * A(2, 2)) * invdet;
-    minv(1, 1) = (A(0, 0) * A(2, 2) - A(0, 2) * A(2, 0)) * invdet;
-    minv(1, 2) = (A(1, 0) * A(0, 2) - A(0, 0) * A(1, 2)) * invdet;
-    minv(2, 0) = (A(1, 0) * A(2, 1) - A(2, 0) * A(1, 1)) * invdet;
-    minv(2, 1) = (A(2, 0) * A(0, 1) - A(0, 0) * A(2, 1)) * invdet;
-    minv(2, 2) = (A(0, 0) * A(1, 1) - A(1, 0) * A(0, 1)) * invdet;
-
-    return minv;
-}
-//=================================================================================================//
 Mat2d getAverageValue(const Mat2d &A, const Mat2d &B)
 {
     Mat2d C = Mat2d::Identity();

src/shared/common/vector_functions.h

@@ -35,17 +35,24 @@ namespace SPH
 Vec2d FirstAxisVector(const Vec2d &zero_vector);
 Vec3d FirstAxisVector(const Vec3d &zero_vector);
 
-Vec3d upgradeToVec3d(const Real &input);
-Vec3d upgradeToVec3d(const Vec2d &input);
-Vec3d upgradeToVec3d(const Vec3d &input);
-Mat3d upgradeToMat3d(const Mat2d &input);
-Mat3d upgradeToMat3d(const Mat3d &input);
+inline Vec3d upgradeToVec3d(const Real &input)
+{
+    return Vec3d(input, 0.0, 0.0);
+};
+inline Vec3d upgradeToVec3d(const Vec2d &input)
+{
+    return Vec3d(input[0], input[1], 0.0);
+};
+inline Vec3d upgradeToVec3d(const Vec3d &input) { return input; };
 
-Vecd degradeToVecd(const Vec3d &input);
-Matd degradeToMatd(const Mat3d &input);
+inline Mat3d upgradeToMat3d(const Mat2d &input)
+{
+    Mat3d output = Mat3d::Zero();
+    output.block<2, 2>(0, 0) = input;
+    return output;
+};
+inline Mat3d upgradeToMat3d(const Mat3d &input) { return input; };
 
-Mat2d getInverse(const Mat2d &A);
-Mat3d getInverse(const Mat3d &A);
 Mat2d getAverageValue(const Mat2d &A, const Mat2d &B);
 Mat3d getAverageValue(const Mat3d &A, const Mat3d &B);
 Mat2d inverseCholeskyDecomposition(const Mat2d &A);