diff --git a/src/for_2D_build/common/data_type.h b/src/for_2D_build/common/data_type.h
index f1a425efa4..4f93dc7720 100644
--- a/src/for_2D_build/common/data_type.h
+++ b/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
\ No newline at end of file
diff --git a/src/for_2D_build/common/vector_functions_2d.cpp b/src/for_2D_build/common/vector_functions_2d.cpp
deleted file mode 100644
index 0bc0ec9ed9..0000000000
--- a/src/for_2D_build/common/vector_functions_2d.cpp
+++ /dev/null
@@ -1,23 +0,0 @@
-#include "vector_functions.h"
-
-namespace SPH
-{
-//=================================================================================================//
-Vec2d degradeToVecd(const Vec3d &input)
-{
- Vec2d output = Vec2d::Zero();
- for (int i = 0; i != Dimensions; i++)
- output[i] = input[i];
- return output;
-}
-//=================================================================================================//
-Mat2d degradeToMatd(const Mat3d &input)
-{
- Mat2d output = Mat2d::Zero();
- for (int i = 0; i != Dimensions; i++)
- for (int j = 0; j != Dimensions; j++)
- output(i, j) = input(i, j);
- return output;
-}
-//=================================================================================================//
-} // namespace SPH
diff --git a/src/for_3D_build/common/data_type.h b/src/for_3D_build/common/data_type.h
index 288fc2d852..540ceb3975 100644
--- a/src/for_3D_build/common/data_type.h
+++ b/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
\ No newline at end of file
diff --git a/src/for_3D_build/common/scalar_functions_3d.cpp b/src/for_3D_build/common/scalar_functions_3d.cpp
index 13687003a4..47c495481c 100644
--- a/src/for_3D_build/common/scalar_functions_3d.cpp
+++ b/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
diff --git a/src/for_3D_build/common/vector_functions_3d.cpp b/src/for_3D_build/common/vector_functions_3d.cpp
deleted file mode 100644
index c77196c40e..0000000000
--- a/src/for_3D_build/common/vector_functions_3d.cpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#include "vector_functions.h"
-
-namespace SPH
-{
-//=================================================================================================//
-Vec3d degradeToVecd(const Vec3d &input)
-{
- return input;
-}
-//=================================================================================================//
-Mat3d degradeToMatd(const Mat3d &input)
-{
- return input;
-}
-//=================================================================================================//
-} // namespace SPH
diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h
index 75cd7f9ef4..cdf03b2e6a 100644
--- a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h
+++ b/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
diff --git a/src/shared/bodies/base_body.cpp b/src/shared/bodies/base_body.cpp
index 71edccef28..40a1a81e44 100644
--- a/src/shared/bodies/base_body.cpp
+++ b/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 *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);
diff --git a/src/shared/bodies/base_body.h b/src/shared/bodies/base_body.h
index 81639414a7..1fe128b969 100644
--- a/src/shared/bodies/base_body.h
+++ b/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 *> 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 *implementation);
int getNewBodyPartID();
int getTotalBodyParts() { return total_body_parts_; };
//----------------------------------------------------------------------
diff --git a/src/shared/bodies/base_body_part.cpp b/src/shared/bodies/base_body_part.cpp
index ec793229fc..b524348b69 100644
--- a/src/shared/bodies/base_body_part.cpp
+++ b/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("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(
+ body_part_name_, body_part_particles_.size(), [&](size_t i) -> Real
+ { return body_part_particles_[i]; });
+ sv_range_size_ = base_particles_.addUniqueSingularVariableOnly(
+ body_part_name_ + "_Size", body_part_particles_.size());
};
//=============================================================================================//
size_t BodyPartByCell::SizeOfLoopRange()
@@ -41,7 +57,6 @@ BodyRegionByParticle::BodyRegionByParticle(SPHBody &sph_body, SharedPtr s
{
shape_ptr_keeper_.assignRef(shape_ptr);
}
-//==
//=================================================================================================//
void BodyRegionByParticle::tagByContain(size_t particle_index)
{
diff --git a/src/shared/bodies/base_body_part.h b/src/shared/bodies/base_body_part.h
index 9e9812ffaf..83acfac919 100644
--- a/src/shared/bodies/base_body_part.h
+++ b/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("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 *dvIndexList() { return dv_index_list_; };
+ SingularVariable *svRangeSize() { return sv_range_size_; };
protected:
SPHBody &sph_body_;
int part_id_;
std::string body_part_name_;
BaseParticles &base_particles_;
+ DiscreteVariable *dv_index_list_;
+ SingularVariable *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)
diff --git a/src/shared/body_relations/base_body_relation.h b/src/shared/body_relations/base_body_relation.h
index a348dddc59..c4f85c6c73 100644
--- a/src/shared/body_relations/base_body_relation.h
+++ b/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
diff --git a/src/shared/body_relations/contact_body_relation.cpp b/src/shared/body_relations/contact_body_relation.cpp
index 62a5b68fed..769bbf9f77 100644
--- a/src/shared/body_relations/contact_body_relation.cpp
+++ b/src/shared/body_relations/contact_body_relation.cpp
@@ -3,6 +3,8 @@
#include "base_particle_dynamics.h"
#include "cell_linked_list.hpp"
+#include
+
namespace SPH
{
//=================================================================================================//
diff --git a/src/shared/common/base_data_package.h b/src/shared/common/base_data_package.h
index d86d55636f..7eae97a163 100644
--- a/src/shared/common/base_data_package.h
+++ b/src/shared/common/base_data_package.h
@@ -72,38 +72,6 @@ using DataContainerAddressAssemble = DataAssemble typename ContainerType>
using DataContainerUniquePtrAssemble = DataAssemble;
-/** a type irrelevant operation on the data assembles */
-template typename OperationType>
-class DataAssembleOperation
-{
- OperationType integer_operation;
- OperationType scalar_operation;
- OperationType vector2d_operation;
- OperationType matrix2d_operation;
- OperationType vector3d_operation;
- OperationType matrix3d_operation;
-
- public:
- template
- DataAssembleOperation(Args &&...args)
- : integer_operation(std::forward(args)...),
- scalar_operation(std::forward(args)...),
- vector2d_operation(std::forward(args)...),
- matrix2d_operation(std::forward(args)...),
- vector3d_operation(std::forward(args)...),
- matrix3d_operation(std::forward(args)...){};
- template
- void operator()(OperationArgs &&...operation_args)
- {
- integer_operation(std::forward(operation_args)...);
- scalar_operation(std::forward(operation_args)...);
- vector2d_operation(std::forward(operation_args)...);
- matrix2d_operation(std::forward(operation_args)...);
- vector3d_operation(std::forward(operation_args)...);
- matrix3d_operation(std::forward(operation_args)...);
- }
-};
-
// Please refer: https://www.cppstories.com/2022/tuple-iteration-basics/ for the following code
template
class OperationOnDataAssemble
diff --git a/src/shared/common/scalar_functions.h b/src/shared/common/scalar_functions.h
index a2a47d394c..44dc525556 100644
--- a/src/shared/common/scalar_functions.h
+++ b/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();
diff --git a/src/shared/common/sphinxsys_containers.h b/src/shared/common/sphinxsys_containers.h
index a4f705baea..7e7fc946ff 100644
--- a/src/shared/common/sphinxsys_containers.h
+++ b/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;
using SPHBodyVector = StdVec;
diff --git a/src/shared/common/sphinxsys_variable.h b/src/shared/common/sphinxsys_variable.h
index ce75a8403e..8ae0184cc9 100644
--- a/src/shared/common/sphinxsys_variable.h
+++ b/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
diff --git a/src/shared/common/vector_functions.cpp b/src/shared/common/vector_functions.cpp
index 57a777c59d..e5046d1cfa 100644
--- a/src/shared/common/vector_functions.cpp
+++ b/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();
diff --git a/src/shared/common/vector_functions.h b/src/shared/common/vector_functions.h
index 9e288ef7c4..3d37853fd3 100644
--- a/src/shared/common/vector_functions.h
+++ b/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);
diff --git a/src/shared/particle_dynamics/base_local_dynamics.h b/src/shared/particle_dynamics/base_local_dynamics.h
index 4d89feaaf9..a6426829f1 100644
--- a/src/shared/particle_dynamics/base_local_dynamics.h
+++ b/src/shared/particle_dynamics/base_local_dynamics.h
@@ -33,8 +33,11 @@
#include "base_data_package.h"
#include "base_particle_dynamics.h"
#include "execution_policy.h"
+#include "reduce_functors.h"
#include "sphinxsys_containers.h"
+#include
+
namespace SPH
{
/**
@@ -59,14 +62,11 @@ class BaseLocalDynamics
sph_body_(identifier.getSPHBody()),
particles_(&sph_body_.getBaseParticles()){};
virtual ~BaseLocalDynamics(){};
+ typedef DynamicsIdentifier Identifier;
DynamicsIdentifier &getDynamicsIdentifier() { return identifier_; };
SPHBody &getSPHBody() { return sph_body_; };
BaseParticles *getParticles() { return particles_; };
virtual void setupDynamics(Real dt = 0.0){}; // setup global parameters
- void registerComputingKernel(Implementation *implementation)
- {
- sph_body_.registerComputingKernel(implementation);
- };
protected:
DynamicsIdentifier &identifier_;
@@ -84,19 +84,22 @@ template
class BaseLocalDynamicsReduce : public BaseLocalDynamics
{
public:
+ typedef Operation OperationType;
+ using ReturnType = typename Operation::ReturnType;
explicit BaseLocalDynamicsReduce(DynamicsIdentifier &identifier)
: BaseLocalDynamics(identifier),
+ reference_(ReduceReference::value),
quantity_name_("ReducedQuantity"){};
virtual ~BaseLocalDynamicsReduce(){};
- using ReturnType = decltype(Operation::reference_);
- ReturnType Reference() { return operation_.reference_; };
+ ReturnType Reference() { return reference_; };
std::string QuantityName() { return quantity_name_; };
Operation &getOperation() { return operation_; };
virtual ReturnType outputResult(ReturnType reduced_value) { return reduced_value; }
protected:
Operation operation_;
+ ReturnType reference_;
std::string quantity_name_;
};
template
diff --git a/src/shared/particle_dynamics/execution/base_implementation.h b/src/shared/particle_dynamics/execution/base_implementation.h
new file mode 100644
index 0000000000..63819c4b4e
--- /dev/null
+++ b/src/shared/particle_dynamics/execution/base_implementation.h
@@ -0,0 +1,59 @@
+/* ------------------------------------------------------------------------- *
+ * 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 base_implementation.h
+ * @brief Here we define the execution policy relevant to parallel computing.
+ * @details This analog of the standard library on the same functions.
+ * @author Alberto Guarnieri and Xiangyu Hu
+ */
+
+#ifndef BASE_IMPLEMENTATION_H
+#define BASE_IMPLEMENTATION_H
+
+#include "sphinxsys_containers.h"
+
+namespace SPH
+{
+namespace execution
+{
+template
+class Implementation;
+
+template <>
+class Implementation
+{
+ public:
+ explicit Implementation() {}
+ ~Implementation() {}
+
+ bool isUpdated() { return is_updated_; };
+ void resetUpdated() { is_updated_ = false; };
+
+ protected:
+ bool is_updated_ = false;
+ void setUpdated() { is_updated_ = true; };
+};
+
+} // namespace execution
+} // namespace SPH
+#endif // BASE_IMPLEMENTATION_H
diff --git a/src/shared/particle_dynamics/execution/execution.h b/src/shared/particle_dynamics/execution/implementation.h
similarity index 75%
rename from src/shared/particle_dynamics/execution/execution.h
rename to src/shared/particle_dynamics/execution/implementation.h
index 80b61413fe..4313b1ef87 100644
--- a/src/shared/particle_dynamics/execution/execution.h
+++ b/src/shared/particle_dynamics/execution/implementation.h
@@ -21,41 +21,25 @@
* *
* ------------------------------------------------------------------------- */
/**
- * @file execution.h
+ * @file implementation.h
* @brief Here we define the execution policy relevant to parallel computing.
* @details This analog of the standard library on the same functions.
* @author Alberto Guarnieri and Xiangyu Hu
*/
-#ifndef EXECUTION_H
-#define EXECUTION_H
+#ifndef IMPLEMENTATION_H
+#define IMPLEMENTATION_H
#include "base_data_type.h"
+#include "base_implementation.h"
#include "execution_policy.h"
+#include "loop_range.h"
#include "ownership.h"
-#include "sphinxsys_containers.h"
namespace SPH
{
namespace execution
{
-template
-class Implementation;
-
-template <>
-class Implementation
-{
- public:
- explicit Implementation() {}
- ~Implementation() {}
-
- bool isUpdated() { return is_updated_; };
- void resetUpdated() { is_updated_ = false; };
-
- protected:
- bool is_updated_ = false;
- void setUpdated() { is_updated_ = true; };
-};
template
class Implementation
@@ -63,25 +47,24 @@ class Implementation
{
public:
explicit Implementation(LocalDynamicsType &local_dynamics)
- : Implementation(),
- local_dynamics_(local_dynamics), computing_kernel_(nullptr) {}
+ : Implementation(), local_dynamics_(local_dynamics),
+ computing_kernel_(nullptr) {}
~Implementation()
{
delete computing_kernel_;
}
template
- ComputingKernelType *getComputingKernel(Args &&... args)
+ ComputingKernelType *getComputingKernel(Args &&...args)
{
if (computing_kernel_ == nullptr)
{
- local_dynamics_.registerComputingKernel(this, std::forward(args)...);
computing_kernel_ = new ComputingKernelType(
- ExecutionPolicy{}, local_dynamics_, std::forward(args)...);
- setUpdated();
+ ExecutionPolicy{}, this->local_dynamics_, std::forward(args)...);
+ this->setUpdated();
}
- if (!isUpdated())
+ if (!this->isUpdated())
{
overwriteComputingKernel(std::forward(args)...);
}
@@ -90,17 +73,17 @@ class Implementation
}
template
- void overwriteComputingKernel(Args &&... args)
+ void overwriteComputingKernel(Args &&...args)
{
*computing_kernel_ = ComputingKernelType(
- ExecutionPolicy{}, local_dynamics_, std::forward(args)...);
- setUpdated();
+ ExecutionPolicy{}, this->local_dynamics_, std::forward(args)...);
+ this->setUpdated();
}
- private:
+ protected:
LocalDynamicsType &local_dynamics_;
ComputingKernelType *computing_kernel_;
};
} // namespace execution
} // namespace SPH
-#endif // EXECUTION_H
+#endif // IMPLEMENTATION_H
diff --git a/src/shared/particle_dynamics/particle_functors.h b/src/shared/particle_dynamics/particle_functors.h
index b9737b9ebe..4c168d33a4 100644
--- a/src/shared/particle_dynamics/particle_functors.h
+++ b/src/shared/particle_dynamics/particle_functors.h
@@ -30,6 +30,7 @@
#define PARTICLE_FUNCTORS_H
#include "base_particles.hpp"
+#include "reduce_functors.h"
namespace SPH
{
@@ -53,7 +54,7 @@ class WithinScope
class AllParticles : public WithinScope
{
public:
- explicit AllParticles(BaseParticles *base_particles) : WithinScope(){};
+ explicit AllParticles(BaseParticles *base_particles) : WithinScope() {};
bool operator()(size_t index_i)
{
return true;
@@ -68,7 +69,7 @@ class IndicatedParticles : public WithinScope
public:
explicit IndicatedParticles(BaseParticles *base_particles)
: WithinScope(),
- indicator_(base_particles->getVariableDataByName("Indicator")){};
+ indicator_(base_particles->getVariableDataByName("Indicator")) {};
bool operator()(size_t index_i)
{
return indicator_[index_i] == INDICATOR;
@@ -85,7 +86,7 @@ class NotIndicatedParticles : public WithinScope
public:
explicit NotIndicatedParticles(BaseParticles *base_particles)
: WithinScope(),
- indicator_(base_particles->getVariableDataByName("Indicator")){};
+ indicator_(base_particles->getVariableDataByName("Indicator")) {};
bool operator()(size_t index_i)
{
return indicator_[index_i] != INDICATOR;
@@ -105,7 +106,7 @@ class ParameterFixed : public ParticleParameter
T parameter_;
public:
- explicit ParameterFixed(const T &c) : ParticleParameter(), parameter_(c){};
+ explicit ParameterFixed(const T &c) : ParticleParameter(), parameter_(c) {};
T operator()(size_t index_i)
{
return parameter_;
@@ -118,7 +119,7 @@ class ParameterVariable : public ParticleParameter
T *v_;
public:
- explicit ParameterVariable(T *v) : ParticleParameter(), v_(v){};
+ explicit ParameterVariable(T *v) : ParticleParameter(), v_(v) {};
T operator()(size_t index_i)
{
return v_[index_i];
@@ -134,13 +135,13 @@ class ParticleAverage // base class to indicate the concept of particle average
template
class PairAverageFixed : public ParticleAverage
{
- const T average_;
+ T average_;
public:
PairAverageFixed(const T &c1, const T &c2)
- : ParticleAverage(), average_(0.5 * (c1 + c2)){};
+ : ParticleAverage(), average_(0.5 * (c1 + c2)) {};
explicit PairAverageFixed(const T &c)
- : PairAverageFixed(c, c){};
+ : PairAverageFixed(c, c) {};
T operator()(size_t index_i, size_t index_j)
{
return average_;
@@ -150,7 +151,7 @@ class PairAverageFixed : public ParticleAverage
class GeomAverage : public ParticleAverage
{
public:
- GeomAverage() : ParticleAverage(){};
+ GeomAverage() : ParticleAverage() {};
protected:
Real inverse(const Real &x) { return 1.0 / x; };
@@ -166,9 +167,9 @@ class PairGeomAverageFixed : public GeomAverage
public:
PairGeomAverageFixed(const T &c1, const T &c2)
- : GeomAverage(), geom_average_(2.0 * c1 * c2 * inverse(c1 + c2)){};
+ : GeomAverage(), geom_average_(2.0 * c1 * c2 * inverse(c1 + c2)) {};
explicit PairGeomAverageFixed(const T &c)
- : PairGeomAverageFixed(c, c){};
+ : PairGeomAverageFixed(c, c) {};
T operator()(size_t index_i, size_t index_j)
{
return geom_average_;
@@ -182,9 +183,9 @@ class PairAverageVariable : public ParticleAverage
public:
PairAverageVariable(T *v1, T *v2)
- : ParticleAverage(), v1_(v1), v2_(v2){};
+ : ParticleAverage(), v1_(v1), v2_(v2) {};
explicit PairAverageVariable(T *v)
- : PairAverageVariable(v, v){};
+ : PairAverageVariable(v, v) {};
T operator()(size_t index_i, size_t index_j)
{
return 0.5 * (v1_[index_i] + v2_[index_j]);
@@ -198,9 +199,9 @@ class PairGeomAverageVariable : public GeomAverage
public:
PairGeomAverageVariable(T *v1, T *v2)
- : GeomAverage(), v1_(v1), v2_(v2){};
+ : GeomAverage(), v1_(v1), v2_(v2) {};
explicit PairGeomAverageVariable(T *v)
- : PairGeomAverageVariable(v, v){};
+ : PairGeomAverageVariable(v, v) {};
T operator()(size_t index_i, size_t index_j)
{
@@ -217,7 +218,7 @@ class KernelCorrection // base class to indicate the concept of kernel correctio
class NoKernelCorrection : public KernelCorrection
{
public:
- NoKernelCorrection(BaseParticles *particles) : KernelCorrection(){};
+ NoKernelCorrection(BaseParticles *particles) : KernelCorrection() {};
Real operator()(size_t index_i)
{
return 1.0;
@@ -229,7 +230,7 @@ class LinearGradientCorrection : public KernelCorrection
public:
LinearGradientCorrection(BaseParticles *particles)
: KernelCorrection(),
- B_(particles->getVariableDataByName("LinearGradientCorrectionMatrix")){};
+ B_(particles->getVariableDataByName("LinearGradientCorrectionMatrix")) {};
Matd operator()(size_t index_i)
{
@@ -243,7 +244,7 @@ class LinearGradientCorrection : public KernelCorrection
class SingleResolution
{
public:
- SingleResolution(BaseParticles *particles){};
+ SingleResolution(BaseParticles *particles) {};
Real operator()(size_t index_i)
{
return 1.0;
@@ -256,7 +257,7 @@ class AdaptiveResolution
{
public:
AdaptiveResolution(BaseParticles *particles)
- : h_ratio_(particles->getVariableDataByName("SmoothingLengthRatio")){};
+ : h_ratio_(particles->getVariableDataByName("SmoothingLengthRatio")) {};
Real operator()(size_t index_i)
{
@@ -266,62 +267,5 @@ class AdaptiveResolution
protected:
Real *h_ratio_;
};
-//----------------------------------------------------------------------
-// Particle reduce functors
-//----------------------------------------------------------------------
-template
-struct ReduceSum
-{
- ReturnType reference_ = ZeroData::value;
- ReturnType operator()(const ReturnType &x, const ReturnType &y) const { return x + y; };
-};
-
-struct ReduceMax
-{
- Real reference_ = MinReal;
- Real operator()(Real x, Real y) const { return SMAX(x, y); };
-};
-
-struct ReduceMin
-{
- Real reference_ = MaxReal;
- Real operator()(Real x, Real y) const { return SMIN(x, y); };
-};
-
-struct ReduceOR
-{
- bool reference_ = false;
- bool operator()(bool x, bool y) const { return x || y; };
-};
-
-struct ReduceAND
-{
- bool reference_ = true;
- bool operator()(bool x, bool y) const { return x && y; };
-};
-
-struct ReduceLowerBound
-{
- Vecd reference_ = MaxReal * Vecd::Ones();
- Vecd operator()(const Vecd &x, const Vecd &y) const
- {
- Vecd lower_bound;
- for (int i = 0; i < lower_bound.size(); ++i)
- lower_bound[i] = SMIN(x[i], y[i]);
- return lower_bound;
- };
-};
-
-struct ReduceUpperBound
-{
- Vecd reference_ = MinReal * Vecd::Ones();
- Vecd operator()(const Vecd &x, const Vecd &y) const
- {
- Vecd upper_bound;
- for (int i = 0; i < upper_bound.size(); ++i)
- upper_bound[i] = SMAX(x[i], y[i]);
- return upper_bound;
- };
-};
} // namespace SPH
#endif // PARTICLE_FUNCTORS_H
diff --git a/src/shared/particle_dynamics/particle_iterators.h b/src/shared/particle_dynamics/particle_iterators.h
index 013e8c790e..ebf57dae57 100644
--- a/src/shared/particle_dynamics/particle_iterators.h
+++ b/src/shared/particle_dynamics/particle_iterators.h
@@ -30,11 +30,9 @@
#define PARTICLE_ITERATORS_H
#include "base_data_package.h"
-#include "execution.h"
+#include "implementation.h"
#include "sphinxsys_containers.h"
-#include
-
namespace SPH
{
using namespace execution;
@@ -245,6 +243,7 @@ inline ReturnType particle_reduce(const ParallelPolicy &par, const IndexVector &
return operation(x, y);
});
};
+
/**
* BodypartByCell-wise reduce iterators (for sequential and parallel computing).
*/
@@ -288,42 +287,5 @@ inline ReturnType particle_reduce(const ParallelPolicy &par, const ConcurrentCel
[&](const ReturnType &x, const ReturnType &y) -> ReturnType
{ return operation(x, y); });
}
-
-template
-T exclusive_scan(const SequencedPolicy &seq_policy, T *first, T *d_first, UnsignedInt d_size, Op op)
-{
- UnsignedInt scan_size = d_size - 1;
- std::exclusive_scan(first, first + d_size, d_first, T{0}, op);
- return d_first[scan_size];
-}
-
-template
-T exclusive_scan(const ParallelPolicy &par_policy, T *first, T *d_first, UnsignedInt d_size, Op op)
-{
- // Exclusive scan is the same as inclusive, but shifted by one
- UnsignedInt scan_size = d_size - 1;
- d_first[0] = T{0};
- using range_type = tbb::blocked_range;
- tbb::parallel_scan(
- range_type(0, scan_size), d_first[0],
- [=](const range_type &r, T sum, bool is_final_scan) -> T
- {
- T tmp = sum;
- for (UnsignedInt i = r.begin(); i < r.end(); ++i)
- {
- tmp = op(tmp, first[i]);
- if (is_final_scan)
- {
- d_first[i + 1] = tmp;
- }
- }
- return tmp;
- },
- [&](const T &a, const T &b)
- {
- return op(a, b);
- });
- return d_first[scan_size];
-}
} // namespace SPH
#endif // PARTICLE_ITERATORS_H
diff --git a/src/shared/particle_dynamics/reduce_functors.h b/src/shared/particle_dynamics/reduce_functors.h
new file mode 100644
index 0000000000..c10b51156c
--- /dev/null
+++ b/src/shared/particle_dynamics/reduce_functors.h
@@ -0,0 +1,134 @@
+/* ------------------------------------------------------------------------- *
+ * 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 reduce_functors.h
+ * @brief TBD.
+ * @author Xiangyu Hu
+ */
+
+#ifndef REDUCE_FUNCTORS_H
+#define REDUCE_FUNCTORS_H
+
+#include "base_data_package.h"
+
+namespace SPH
+{
+template
+struct ReduceReference;
+
+template
+struct ReduceBase
+{
+ typedef DataType ReturnType;
+};
+
+template
+struct ReduceSum : ReduceBase
+{
+ DataType operator()(const DataType &x, const DataType &y) const { return x + y; };
+};
+
+template
+struct ReduceReference>
+{
+ static inline DataType value = ZeroData::value;
+};
+
+struct ReduceMax : ReduceBase
+{
+ Real operator()(Real x, Real y) const { return SMAX(x, y); };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline Real value = MinReal;
+};
+
+struct ReduceMin : ReduceBase
+{
+ Real operator()(Real x, Real y) const { return SMIN(x, y); };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline Real value = MaxReal;
+};
+struct ReduceOR : ReduceBase
+{
+ bool operator()(bool x, bool y) const { return x || y; };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline bool value = false;
+};
+
+struct ReduceAND : ReduceBase
+{
+ bool operator()(bool x, bool y) const { return x && y; };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline bool value = true;
+};
+struct ReduceLowerBound : ReduceBase
+{
+ Vecd operator()(const Vecd &x, const Vecd &y) const
+ {
+ Vecd lower_bound;
+ for (int i = 0; i < lower_bound.size(); ++i)
+ lower_bound[i] = SMIN(x[i], y[i]);
+ return lower_bound;
+ };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline Vecd value = MaxReal * Vecd::Ones();
+};
+
+struct ReduceUpperBound : ReduceBase
+{
+ Vecd operator()(const Vecd &x, const Vecd &y) const
+ {
+ Vecd upper_bound;
+ for (int i = 0; i < upper_bound.size(); ++i)
+ upper_bound[i] = SMAX(x[i], y[i]);
+ return upper_bound;
+ };
+};
+
+template <>
+struct ReduceReference
+{
+ static inline Vecd value = MinReal * Vecd::Ones();
+};
+
+} // namespace SPH
+#endif // REDUCE_FUNCTORS_H
diff --git a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h
index 2b57bf6173..9587d91e73 100644
--- a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h
+++ b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h
@@ -148,41 +148,6 @@ class ConstrainSolidBodyMassCenter : public MotionConstraint
void update(size_t index_i, Real dt = 0.0);
};
-struct SimbodyState
-{
- Vec3d initial_origin_location_;
- Vec3d origin_location_, origin_velocity_, origin_acceleration_;
- Vec3d angular_velocity_, angular_acceleration_;
- Mat3d rotation_;
-
- SimbodyState() : initial_origin_location_(ZeroData::value),
- origin_location_(ZeroData::value),
- origin_velocity_(ZeroData::value),
- origin_acceleration_(ZeroData::value),
- angular_velocity_(ZeroData::value),
- angular_acceleration_(ZeroData::value),
- rotation_(Mat3d::Identity()) {}
-
- // implemented according to the Simbody API function with the same name
- void findStationLocationVelocityAndAccelerationInGround(const Vec3d &initial_location,
- const Vec3d &initial_normal,
- Vec3d &locationOnGround,
- Vec3d &velocityInGround,
- Vec3d &accelerationInGround,
- Vec3d &normalInGround)
- {
- Vec3d temp_location = rotation_ * (initial_location - initial_origin_location_);
- locationOnGround = origin_location_ + temp_location;
-
- Vec3d temp_velocity = angular_velocity_.cross(temp_location);
- velocityInGround = origin_velocity_ + temp_velocity;
- accelerationInGround = origin_acceleration_ +
- angular_acceleration_.cross(temp_location) +
- angular_velocity_.cross(temp_velocity);
- normalInGround = rotation_ * initial_normal;
- };
-};
-
/**
* @class ConstraintBySimBody
* @brief Constrain by the motion computed from Simbody.
diff --git a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.hpp b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.hpp
index 424bcabf7f..b51dea279b 100644
--- a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.hpp
+++ b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.hpp
@@ -79,7 +79,6 @@ void ConstraintBySimBody::updateSimbodyState(const SimTK::St
template
void ConstraintBySimBody::update(size_t index_i, Real dt)
{
- /** Change to SimTK::Vector. */
Vec3d pos, vel, acc, n;
simbody_state_.findStationLocationVelocityAndAccelerationInGround(
upgradeToVec3d(this->pos0_[index_i]), upgradeToVec3d(n0_[index_i]), pos, vel, acc, n);
diff --git a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.h
index da19abf5ad..8d2a9ff31a 100644
--- a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.h
+++ b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.h
@@ -147,7 +147,7 @@ class BaseElasticIntegration : public LocalDynamics, public DataDelegateInner
/**
* @class BaseIntegration1stHalf
- * @brief computing stress relaxation process by verlet time stepping
+ * @brief computing stress relaxation process by Verlet time stepping
* This is the first step
*/
class BaseIntegration1stHalf : public BaseElasticIntegration
@@ -167,7 +167,7 @@ class BaseIntegration1stHalf : public BaseElasticIntegration
/**
* @class Integration1stHalf
- * @brief computing stress relaxation process by verlet time stepping
+ * @brief computing stress relaxation process by Verlet time stepping
* This is the first step
*/
class Integration1stHalf : public BaseIntegration1stHalf
@@ -321,7 +321,7 @@ class Integration1stHalfPK2RightCauchy : public Integration1stHalfPK2
/**
* @class Integration2ndHalf
- * @brief computing stress relaxation process by verlet time stepping
+ * @brief computing stress relaxation process by Verlet time stepping
* This is the second step
*/
class Integration2ndHalf : public BaseElasticIntegration
diff --git a/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp
index 96e03942a3..2cc601a7b2 100644
--- a/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp
+++ b/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp
@@ -1,6 +1,8 @@
#include "loading_dynamics.h"
#include "base_general_dynamics.h"
+#include
+
namespace SPH
{
namespace solid_dynamics
diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h
index bd9affc9da..a8f0f4faa0 100644
--- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h
+++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h
@@ -259,7 +259,7 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
/**
* @class ShellStressRelaxationSecondHalf
- * @brief computing stress relaxation process by verlet time stepping
+ * @brief computing stress relaxation process by Verlet time stepping
* This is the second step
*/
class ShellStressRelaxationSecondHalf : public BaseShellRelaxation
diff --git a/src/shared/particles/base_particles.cpp b/src/shared/particles/base_particles.cpp
index 9135c88b1c..462451a003 100644
--- a/src/shared/particles/base_particles.cpp
+++ b/src/shared/particles/base_particles.cpp
@@ -10,7 +10,7 @@ namespace SPH
{
//=================================================================================================//
BaseParticles::BaseParticles(SPHBody &sph_body, BaseMaterial *base_material)
- : v_total_real_particles_(nullptr), real_particles_bound_(0), particles_bound_(0),
+ : sv_total_real_particles_(nullptr), real_particles_bound_(0), particles_bound_(0),
original_id_(nullptr), sorted_id_(nullptr),
pos_(nullptr), Vol_(nullptr), rho_(nullptr), mass_(nullptr),
sph_body_(sph_body), body_name_(sph_body.getName()),
@@ -23,7 +23,7 @@ BaseParticles::BaseParticles(SPHBody &sph_body, BaseMaterial *base_material)
read_restart_variable_from_xml_(variables_to_restart_, restart_xml_parser_)
{
sph_body.assignBaseParticles(this);
- v_total_real_particles_ = registerSingularVariable("TotalRealParticles");
+ sv_total_real_particles_ = registerSingularVariable("TotalRealParticles");
}
//=================================================================================================//
void BaseParticles::initializeBasicParticleVariables()
@@ -58,7 +58,7 @@ void BaseParticles::registerPositionAndVolumetricMeasureFromReload()
//=================================================================================================//
void BaseParticles::initializeAllParticlesBounds(size_t number_of_particles)
{
- UnsignedInt *total_real_particles = v_total_real_particles_->ValueAddress();
+ UnsignedInt *total_real_particles = sv_total_real_particles_->ValueAddress();
*total_real_particles = number_of_particles;
real_particles_bound_ = number_of_particles;
particles_bound_ = real_particles_bound_;
diff --git a/src/shared/particles/base_particles.h b/src/shared/particles/base_particles.h
index 7fcbb9f3b7..7b6968cdfd 100644
--- a/src/shared/particles/base_particles.h
+++ b/src/shared/particles/base_particles.h
@@ -99,7 +99,7 @@ class BaseParticles
// particles_bound_ gives the total number of particles in all groups.
//----------------------------------------------------------------------
protected:
- SingularVariable *v_total_real_particles_;
+ SingularVariable *sv_total_real_particles_;
UnsignedInt real_particles_bound_;
UnsignedInt particles_bound_;
@@ -109,9 +109,10 @@ class BaseParticles
//----------------------------------------------------------------------
// Generalized particle manipulation
//----------------------------------------------------------------------
- UnsignedInt TotalRealParticles() { return *v_total_real_particles_->ValueAddress(); };
- void incrementTotalRealParticles(UnsignedInt increment = 1) { *v_total_real_particles_->ValueAddress() += increment; };
- void decrementTotalRealParticles(UnsignedInt decrement = 1) { *v_total_real_particles_->ValueAddress() -= decrement; };
+ SingularVariable *svTotalRealParticles() { return sv_total_real_particles_; };
+ UnsignedInt TotalRealParticles() { return *sv_total_real_particles_->ValueAddress(); };
+ void incrementTotalRealParticles(UnsignedInt increment = 1) { *sv_total_real_particles_->ValueAddress() += increment; };
+ void decrementTotalRealParticles(UnsignedInt decrement = 1) { *sv_total_real_particles_->ValueAddress() -= decrement; };
UnsignedInt RealParticlesBound() { return real_particles_bound_; };
UnsignedInt ParticlesBound() { return particles_bound_; };
void initializeAllParticlesBounds(size_t total_real_particles);
@@ -136,6 +137,8 @@ class BaseParticles
public:
template
DataType *addUniqueDiscreteVariable(const std::string &name, size_t data_size, Args &&...args);
+ template
+ DiscreteVariable *addUniqueDiscreteVariableOnly(const std::string &name, size_t data_size, Args &&...args);
template
DataType *registerDiscreteVariable(const std::string &name, size_t data_size, Args &&...args);
@@ -158,7 +161,11 @@ class BaseParticles
DiscreteVariable *registerDiscreteVariableOnly(const std::string &name, size_t data_size, Args &&...args);
template
DiscreteVariable *registerStateVariableOnly(const std::string &name, Args &&...args);
+ template
+ DiscreteVariable *registerStateVariableOnlyFrom(const std::string &new_name, const std::string &old_name);
+ template
+ SingularVariable *addUniqueSingularVariableOnly(const std::string &name, DataType initial_value = ZeroData::value);
template
SingularVariable *registerSingularVariable(const std::string &name, DataType initial_value = ZeroData::value);
template
diff --git a/src/shared/particles/base_particles.hpp b/src/shared/particles/base_particles.hpp
index fc6cefe2c2..14131bf73c 100644
--- a/src/shared/particles/base_particles.hpp
+++ b/src/shared/particles/base_particles.hpp
@@ -64,6 +64,17 @@ DataType *BaseParticles::
return variable->DataField();
}
//=================================================================================================//
+template
+DiscreteVariable *BaseParticles::
+ addUniqueDiscreteVariableOnly(const std::string &name, size_t data_size, Args &&...args)
+{
+
+ DiscreteVariable *variable =
+ unique_variable_ptrs_.createPtr>(name, data_size);
+ initializeVariable(variable, std::forward(args)...);
+ return variable;
+}
+//=================================================================================================//
template
DataType *BaseParticles::registerDiscreteVariable(const std::string &name,
size_t data_size, Args &&...args)
@@ -143,6 +154,24 @@ DataType *BaseParticles::registerStateVariableFrom(
}
//=================================================================================================//
template
+DiscreteVariable *BaseParticles::registerStateVariableOnlyFrom(
+ const std::string &new_name, const std::string &old_name)
+{
+ DiscreteVariable *variable = findVariableByName(all_discrete_variables_, old_name);
+
+ if (variable == nullptr)
+ {
+ std::cout << "\nError: the old variable '" << old_name << "' is not registered!\n";
+ std::cout << __FILE__ << ':' << __LINE__ << std::endl;
+ exit(1);
+ }
+
+ DataType *old_data_field = variable->DataField();
+ return registerStateVariableOnly(new_name, [&](size_t index)
+ { return old_data_field[index]; });
+}
+//=================================================================================================//
+template
DataType *BaseParticles::registerStateVariableFrom(
const std::string &name, const StdLargeVec &geometric_data)
{
@@ -198,6 +227,16 @@ DataType *BaseParticles::getVariableDataByName(const std::string &name)
return variable->DataField();
}
//=================================================================================================//
+template
+SingularVariable *BaseParticles::
+ addUniqueSingularVariableOnly(const std::string &name, DataType initial_value)
+{
+
+ SingularVariable *variable =
+ unique_variable_ptrs_.createPtr>(name, initial_value);
+ return variable;
+}
+//=================================================================================================//
template
SingularVariable *BaseParticles::
registerSingularVariable(const std::string &name, DataType initial_value)
diff --git a/src/shared/shared_ck/body_relation/relation_ck.h b/src/shared/shared_ck/body_relation/relation_ck.h
index 7f804158f5..cf6d957e69 100644
--- a/src/shared/shared_ck/body_relation/relation_ck.h
+++ b/src/shared/shared_ck/body_relation/relation_ck.h
@@ -31,7 +31,7 @@
#include "base_body.h"
#include "base_particles.h"
-#include "execution.h"
+#include "implementation.h"
namespace SPH
{
diff --git a/src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h b/src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h
index 9e0e553eef..da656a6c1e 100644
--- a/src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h
+++ b/src/shared/shared_ck/particle_dynamics/all_shared_physical_dynamics_ck.h
@@ -31,12 +31,11 @@
#ifndef ALL_SHARED_PHYSICAL_DYNAMICS_CK_H
#define ALL_SHARED_PHYSICAL_DYNAMICS_CK_H
-#include "acoustic_step_1st_half.hpp"
-#include "acoustic_step_2nd_half.hpp"
+#include "all_fluid_structure_interactions.h"
#include "all_general_dynamics_ck.h"
+#include "all_shared_fluid_dynamics_ck.h"
+#include "all_solid_dynamics_ck.h"
#include "complex_algorithms_ck.h"
-#include "density_regularization.hpp"
-#include "fluid_time_step_ck.hpp"
#include "interaction_algorithms_ck.hpp"
#include "particle_sort_ck.hpp"
#include "simple_algorithms_ck.h"
diff --git a/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_body_relation.hpp b/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_body_relation.hpp
index a5e4c1b70c..819df491ba 100644
--- a/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_body_relation.hpp
+++ b/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_body_relation.hpp
@@ -4,6 +4,7 @@
#include "update_body_relation.h"
#include "cell_linked_list.hpp"
+#include "particle_iterators_ck.h"
namespace SPH
{
diff --git a/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_cell_linked_list.hpp b/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_cell_linked_list.hpp
index f422e6d026..84db456503 100644
--- a/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_cell_linked_list.hpp
+++ b/src/shared/shared_ck/particle_dynamics/configuration_dynamics/update_cell_linked_list.hpp
@@ -6,6 +6,7 @@
#include "adaptation.hpp"
#include "base_particles.hpp"
#include "mesh_iterators.hpp"
+#include "particle_iterators_ck.h"
namespace SPH
{
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/acoustic_step_2nd_half.h b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/acoustic_step_2nd_half.h
index 25ef020009..d45aa80e9b 100644
--- a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/acoustic_step_2nd_half.h
+++ b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/acoustic_step_2nd_half.h
@@ -88,6 +88,9 @@ class AcousticStep2ndHalf
+ friend class FSI::PressureForceFromFluid;
+
KernelCorrectionType kernel_correction_;
RiemannSolverType riemann_solver_;
};
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/all_shared_fluid_dynamics_ck.h b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/all_shared_fluid_dynamics_ck.h
new file mode 100644
index 0000000000..a99086e10d
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/all_shared_fluid_dynamics_ck.h
@@ -0,0 +1,40 @@
+/* ------------------------------------------------------------------------- *
+ * 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_shared_fluid_dynamics_ck.h
+ * @brief Head file for all shared physics dynamics for both 2- and 3D build.
+ * This is the header file that user code should include to pick up all
+ particle dynamics capabilities.
+ * @author Chi Zhang and Xiangyu Hu
+ */
+
+#ifndef ALL_SHARED_FLUID_DYNAMICS_CK_H
+#define ALL_SHARED_FLUID_DYNAMICS_CK_H
+
+#include "acoustic_step_1st_half.hpp"
+#include "acoustic_step_2nd_half.hpp"
+#include "density_regularization.hpp"
+#include "fluid_time_step_ck.hpp"
+#include "viscous_force.hpp"
+
+#endif // ALL_SHARED_FLUID_DYNAMICS_CK_H
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/density_regularization.h b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/density_regularization.h
index 5e312dd809..4c76643188 100644
--- a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/density_regularization.h
+++ b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/density_regularization.h
@@ -187,6 +187,7 @@ class DensityRegularization>
StdVec *> dv_contact_mass_;
};
+using DensityRegularizationComplex = DensityRegularization, Contact<>>;
using DensityRegularizationComplexFreeSurface = DensityRegularization, Contact<>>;
} // namespace fluid_dynamics
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.h b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.h
new file mode 100644
index 0000000000..6304720735
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.h
@@ -0,0 +1,178 @@
+/* ------------------------------------------------------------------------- *
+ * 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 viscous_force.h
+ * @brief Here, we define the algorithm classes for computing viscous forces in fluids.
+ * @details TBD.
+ * @author Xiangyu Hu
+ */
+
+#ifndef VISCOUS_FORCE_H
+#define VISCOUS_FORCE_H
+
+#include "base_fluid_dynamics.h"
+#include "force_prior_ck.h"
+#include "interaction_ck.hpp"
+#include "kernel_correction_ck.hpp"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+template
+class ViscousForceCK;
+
+template class RelationType, typename... Parameters>
+class ViscousForceCK>
+ : public Interaction>, public ForcePriorCK
+{
+ using ViscosityKernel = typename ViscosityType::ComputingKernel;
+ using CorrectionKernel = typename KernelCorrectionType::ComputingKernel;
+
+ public:
+ template
+ explicit ViscousForceCK(BaseRelationType &base_relation);
+ virtual ~ViscousForceCK() {};
+
+ class InteractKernel
+ : public Interaction>::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, Args &&...args);
+
+ protected:
+ ViscosityKernel viscosity_;
+ CorrectionKernel correction_;
+ Real *Vol_;
+ Vecd *vel_, *viscous_force_;
+ Real smoothing_length_sq_;
+ };
+
+ protected:
+ template
+ friend class FSI::ViscousForceFromFluid;
+
+ ViscosityType viscosity_method_;
+ KernelCorrectionType kernel_correction_;
+ DiscreteVariable *dv_Vol_;
+ DiscreteVariable *dv_vel_, *dv_viscous_force_;
+ Real smoothing_length_sq_;
+};
+
+template
+class ViscousForceCK>
+ : public ViscousForceCK>
+{
+ using BaseViscousForceType = ViscousForceCK>;
+
+ public:
+ explicit ViscousForceCK(Relation> &inner_relation)
+ : BaseViscousForceType(inner_relation) {};
+ virtual ~ViscousForceCK() {};
+
+ class InteractKernel : public BaseViscousForceType::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser)
+ : BaseViscousForceType::InteractKernel(ex_policy, encloser){};
+ void interact(size_t index_i, Real dt = 0.0);
+ };
+};
+
+template
+class ViscousForceCK>
+ : public ViscousForceCK>
+{
+ using BaseViscousForceType = ViscousForceCK>;
+
+ public:
+ explicit ViscousForceCK(Relation> &contact_relation)
+ : BaseViscousForceType(contact_relation) {};
+ virtual ~ViscousForceCK() {};
+
+ class InteractKernel : public BaseViscousForceType::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, size_t contact_index)
+ : BaseViscousForceType::InteractKernel(ex_policy, encloser, contact_index),
+ wall_Vol_(encloser.dv_wall_Vol_[contact_index]->DelegatedDataField(ex_policy)),
+ wall_vel_ave_(encloser.dv_wall_vel_ave_[contact_index]->DelegatedDataField(ex_policy)){};
+ void interact(size_t index_i, Real dt = 0.0);
+
+ protected:
+ Real *wall_Vol_;
+ Vecd *wall_vel_ave_;
+ };
+};
+
+template
+class Viscosity;
+
+template <>
+class Viscosity
+{
+ public:
+ Viscosity(BaseParticles *particles)
+ : mu_(DynamicCast(this, particles->getBaseMaterial()).ReferenceViscosity()) {};
+
+ class ComputingKernel : public ParameterFixed
+ {
+ public:
+ template
+ ComputingKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser)
+ : ParameterFixed(encloser.mu_){};
+ };
+
+ protected:
+ Real mu_;
+};
+
+template <>
+class Viscosity
+{
+ public:
+ Viscosity(BaseParticles *particles)
+ : dv_mu_(particles->getVariableByName("VariableViscosity")) {};
+
+ class ComputingKernel : public ParameterVariable
+ {
+ public:
+ template
+ ComputingKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser)
+ : ParameterVariable(encloser.dv_mu_->DelegatedDataField(ex_policy)){};
+ };
+
+ protected:
+ DiscreteVariable *dv_mu_;
+};
+
+using ViscousForceInnerCK = ViscousForceCK, NoKernelCorrectionCK>>;
+using ViscousForceWithWallCK = ViscousForceCK, NoKernelCorrectionCK>,
+ Contact, NoKernelCorrectionCK>>;
+} // namespace fluid_dynamics
+} // namespace SPH
+#endif // VISCOUS_FORCE_H
\ No newline at end of file
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.hpp b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.hpp
new file mode 100644
index 0000000000..84ae9aae4b
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_dynamics/viscous_force.hpp
@@ -0,0 +1,77 @@
+#pragma once
+
+#include "viscous_force.h"
+
+namespace SPH
+{
+namespace fluid_dynamics
+{
+//=================================================================================================//
+template class RelationType, typename... Parameters>
+template
+ViscousForceCK>::
+ ViscousForceCK(BaseRelationType &base_relation)
+ : Interaction>(base_relation),
+ ForcePriorCK(this->particles_, "ViscousForce"),
+ viscosity_method_(this->particles_), kernel_correction_(this->particles_),
+ dv_Vol_(this->particles_->template getVariableByName("VolumetricMeasure")),
+ dv_vel_(this->particles_->template getVariableByName("Velocity")),
+ dv_viscous_force_(this->dv_current_force_),
+ smoothing_length_sq_(pow(this->sph_body_.sph_adaptation_->ReferenceSmoothingLength(), 2)) {}
+//=================================================================================================//
+template class RelationType, typename... Parameters>
+template
+ViscousForceCK>::InteractKernel::
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, Args &&...args)
+ : Interaction>::InteractKernel(ex_policy, encloser, std::forward(args)...),
+ viscosity_(ex_policy, encloser.viscosity_method_),
+ correction_(ex_policy, encloser.kernel_correction_),
+ Vol_(encloser.dv_Vol_->DelegatedDataField(ex_policy)),
+ vel_(encloser.dv_vel_->DelegatedDataField(ex_policy)),
+ viscous_force_(encloser.dv_viscous_force_->DelegatedDataField(ex_policy)),
+ smoothing_length_sq_(encloser.smoothing_length_sq_) {}
+//=================================================================================================//
+template
+void ViscousForceCK>::
+ InteractKernel::interact(size_t index_i, Real dt)
+{
+ Vecd force = Vecd::Zero();
+ for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+ {
+ UnsignedInt index_j = this->neighbor_index_[n];
+ Vecd e_ij = this->e_ij(index_i, index_j);
+ Vecd vec_r_ij = this->vec_r_ij(index_i, index_j);
+ Vecd vel_derivative = (this->vel_[index_i] - this->vel_[index_j]) /
+ (vec_r_ij.squaredNorm() + 0.01 * this->smoothing_length_sq_);
+
+ force += vec_r_ij.dot((this->correction_(index_i) + this->correction_(index_j)) * e_ij) *
+ harmonic_average(this->viscosity_(index_i), this->viscosity_(index_j)) *
+ vel_derivative * this->dW_ij(index_i, index_j) * this->Vol_[index_j];
+ }
+ this->viscous_force_[index_i] = force * this->Vol_[index_i];
+}
+//=================================================================================================//
+template
+void ViscousForceCK>::
+ InteractKernel::interact(size_t index_i, Real dt)
+{
+ Vecd force = Vecd::Zero();
+ for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+ {
+ UnsignedInt index_j = this->neighbor_index_[n];
+ Vecd e_ij = this->e_ij(index_i, index_j);
+ Vecd vec_r_ij = this->vec_r_ij(index_i, index_j);
+ Vecd vel_derivative = (this->vel_[index_i] - this->wall_vel_ave_[index_j]) /
+ (vec_r_ij.squaredNorm() + 0.01 * this->smoothing_length_sq_);
+
+ force += vec_r_ij.dot(this->correction_(index_i) * e_ij) *
+ this->viscosity_(index_i) * vel_derivative *
+ this->dW_ij(index_i, index_j) * this->wall_Vol_[index_j];
+ }
+ this->viscous_force_[index_i] += force * this->Vol_[index_i];
+}
+//=================================================================================================//
+} // namespace fluid_dynamics
+} // namespace SPH
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/all_fluid_structure_interactions.h b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/all_fluid_structure_interactions.h
new file mode 100644
index 0000000000..75876cb11c
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/all_fluid_structure_interactions.h
@@ -0,0 +1,29 @@
+/* ------------------------------------------------------------------------- *
+ * 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_fluid_structure_interactions.h
+ */
+
+#pragma once
+
+#include "force_on_structure.hpp"
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.h b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.h
new file mode 100644
index 0000000000..684dbb5b20
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.h
@@ -0,0 +1,152 @@
+/* ------------------------------------------------------------------------- *
+ * 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 force_on_structure.h
+ * @brief Here, we define the algorithm classes for fluid structure interaction.
+ * @author Xiangyu Hu
+ */
+
+#ifndef FORCE_ON_STRUCTRUCTURE_H
+#define FORCE_ON_STRUCTRUCTURE_H
+
+#include "base_material.h"
+#include "force_prior_ck.hpp"
+#include "interaction_ck.hpp"
+#include "riemann_solver.h"
+
+namespace SPH
+{
+namespace FSI
+{
+template
+class ForceFromFluid : public Interaction>, public ForcePriorCK
+{
+ using CorrectionKernel = typename KernelCorrectionType::ComputingKernel;
+
+ public:
+ template
+ explicit ForceFromFluid(ContactRelationType &contact_relation, const std::string &force_name);
+ virtual ~ForceFromFluid(){};
+
+ class InteractKernel
+ : public Interaction>::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index);
+
+ protected:
+ Real *Vol_;
+ Vecd *force_from_fluid_, *vel_ave_;
+ CorrectionKernel contact_correction_;
+ Real *contact_Vol_;
+ Vecd *contact_vel_;
+ };
+
+ protected:
+ Solid &solid_;
+ DiscreteVariable *dv_Vol_;
+ DiscreteVariable *dv_force_from_fluid_, *dv_vel_ave_;
+
+ StdVec contact_fluid_;
+ StdVec contact_kernel_correction_;
+ StdVec *> contact_Vol_;
+ StdVec *> contact_vel_;
+};
+
+template
+class ViscousForceFromFluid;
+
+template
+class ViscousForceFromFluid>
+ : public ForceFromFluid
+{
+
+ using ViscosityType = decltype(ViscousForceType::viscosity_method_);
+ using ViscosityKernel = typename ViscosityType::ComputingKernel;
+ using BaseForceFromFluid = ForceFromFluid;
+
+ public:
+ template
+ explicit ViscousForceFromFluid(ContactRelationType &contact_relation);
+ virtual ~ViscousForceFromFluid(){};
+ class InteractKernel : public BaseForceFromFluid::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index);
+ void interact(size_t index_i, Real dt = 0.0);
+
+ protected:
+ ViscosityKernel viscosity_;
+ Real smoothing_length_sq_;
+ };
+
+ protected:
+ StdVec contact_viscosity_method_;
+ StdVec contact_smoothing_length_sq_;
+};
+template
+using ViscousForceOnStructure = ViscousForceFromFluid>;
+
+template
+class PressureForceFromFluid;
+
+template
+class PressureForceFromFluid>
+ : public ForceFromFluid
+{
+ using RiemannSolverType = decltype(AcousticStep2ndHalfType::riemann_solver_);
+ using KernelCorrectionType = decltype(AcousticStep2ndHalfType::kernel_correction_);
+ using BaseForceFromFluid = ForceFromFluid;
+
+ public:
+ template
+ explicit PressureForceFromFluid(ContactRelationType &contact_relation);
+ virtual ~PressureForceFromFluid(){};
+
+ class InteractKernel : public BaseForceFromFluid::InteractKernel
+ {
+ public:
+ template
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index);
+ void interact(size_t index_i, Real dt = 0.0);
+
+ protected:
+ Vecd *acc_ave_, *n_;
+ RiemannSolverType riemann_solver_;
+ Real *contact_rho_, *contact_mass_, *contact_p_;
+ Vecd *contact_force_prior_;
+ };
+
+ protected:
+ DiscreteVariable *dv_acc_ave_, *dv_n_;
+ StdVec contact_riemann_solver_;
+ StdVec *> dv_contact_rho_, dv_contact_mass_, dv_contact_p_;
+ StdVec *> dv_contact_force_prior_;
+};
+template
+using PressureForceOnStructure = PressureForceFromFluid>;
+} // namespace FSI
+} // namespace SPH
+#endif // FORCE_ON_STRUCTRUCTURE_H
\ No newline at end of file
diff --git a/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.hpp b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.hpp
new file mode 100644
index 0000000000..a39bec2504
--- /dev/null
+++ b/src/shared/shared_ck/particle_dynamics/fluid_structure_interaction/force_on_structure.hpp
@@ -0,0 +1,140 @@
+#ifndef FORCE_ON_STRUCTURE_HPP
+#define FORCE_ON_STRUCTURE_HPP
+
+#include "force_on_structure.h"
+
+namespace SPH
+{
+namespace FSI
+{
+//=================================================================================================//
+template
+template
+ForceFromFluid::
+ ForceFromFluid(ContactRelationType &contact_relation, const std::string &force_name)
+ : Interaction>(contact_relation), ForcePriorCK(this->particles_, force_name),
+ solid_(DynamicCast(this, this->sph_body_.getBaseMaterial())),
+ dv_Vol_(this->particles_->template getVariableByName("VolumetricMeasure")),
+ dv_force_from_fluid_(this->dv_current_force_),
+ dv_vel_ave_(solid_.AverageVelocityVariable(this->particles_))
+{
+ for (size_t k = 0; k != this->contact_particles_.size(); ++k)
+ {
+ contact_fluid_.push_back(DynamicCast(this, &this->contact_particles_[k]->getBaseMaterial()));
+ contact_kernel_correction_.push_back(KernelCorrectionType(this->contact_particles_[k]));
+ contact_Vol_.push_back(this->contact_particles_[k]->template getVariableByName("VolumetricMeasure"));
+ contact_vel_.push_back(this->contact_particles_[k]->template getVariableByName("Velocity"));
+ }
+}
+//=================================================================================================//
+template
+template
+ForceFromFluid::InteractKernel::
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index)
+ : Interaction>::InteractKernel(ex_policy, encloser, contact_index),
+ Vol_(encloser.dv_Vol_->DelegatedDataField(ex_policy)),
+ force_from_fluid_(encloser.dv_force_from_fluid_->DelegatedDataField(ex_policy)),
+ vel_ave_(encloser.dv_vel_ave_->DelegatedDataField(ex_policy)),
+ contact_correction_(ex_policy, encloser.contact_kernel_correction_[contact_index]),
+ contact_Vol_(encloser.contact_Vol_[contact_index]->DelegatedDataField(ex_policy)),
+ contact_vel_(encloser.contact_vel_[contact_index]->DelegatedDataField(ex_policy)) {}
+//=================================================================================================//
+template
+template
+ViscousForceFromFluid>::
+ ViscousForceFromFluid(ContactRelationType &contact_relation)
+ : BaseForceFromFluid(contact_relation, "ViscousForceFromFluid")
+{
+ for (size_t k = 0; k != this->contact_particles_.size(); ++k)
+ {
+ contact_viscosity_method_.push_back(ViscosityType(this->contact_particles_[k]));
+ contact_smoothing_length_sq_.push_back(pow(this->contact_bodies_[k]->sph_adaptation_->ReferenceSmoothingLength(), 2));
+ }
+}
+//=================================================================================================//
+template
+template
+ViscousForceFromFluid>::InteractKernel::
+ InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index)
+ : BaseForceFromFluid::InteractKernel(ex_policy, encloser, contact_index),
+ viscosity_(ex_policy, encloser.contact_viscosity_method_[contact_index]),
+ smoothing_length_sq_(encloser.contact_smoothing_length_sq_[contact_index]) {}
+//=================================================================================================//
+template
+void ViscousForceFromFluid>::
+ InteractKernel::interact(size_t index_i, Real dt)
+{
+ Vecd force = Vecd::Zero();
+ for (UnsignedInt n = this->FirstNeighbor(index_i); n != this->LastNeighbor(index_i); ++n)
+ {
+ UnsignedInt index_j = this->neighbor_index_[n];
+ Vecd e_ij = this->e_ij(index_i, index_j);
+ Vecd vec_r_ij = this->vec_r_ij(index_i, index_j);
+ Vecd vel_derivative = (this->vel_ave_[index_i] - this->contact_vel_[index_j]) /
+ (vec_r_ij.squaredNorm() + 0.01 * smoothing_length_sq_);
+
+ force += vec_r_ij.dot(this->contact_correction_(index_j) * e_ij) *
+ viscosity_(index_j) * vel_derivative *
+ this->dW_ij(index_i, index_j) * this->contact_Vol_[index_j];
+ }
+
+ this->force_from_fluid_[index_i] = force * this->Vol_[index_i];
+}
+//=================================================================================================//
+template
+template
+PressureForceFromFluid>::
+ PressureForceFromFluid(ContactRelationType &contact_relation)
+ : BaseForceFromFluid(contact_relation, "PressureForceFromFluid"),
+ dv_acc_ave_(this->solid_.AverageAccelerationVariable(this->particles_)),
+ dv_n_(this->particles_->template getVariableByName("NormalDirection"))
+{
+ for (size_t k = 0; k != this->contact_particles_.size(); ++k)
+ {
+ contact_riemann_solver_.push_back(RiemannSolverType(*this->contact_fluid_[k], *this->contact_fluid_[k]));
+ dv_contact_rho_.push_back(this->contact_particles_[k]->template getVariableByName("Density"));
+ dv_contact_mass_.push_back(this->contact_particles_[k]->template getVariableByName