Sync/2022.03.08 finalization merge from upstream

SPHINXsys/src/for_2D_build/common/data_type.h

@@ -15,6 +15,8 @@ namespace SPH {
 	using SymMatd = SymMat2d;
 	using AngularVecd = Real;
 
+	using  BoundingBox = BoundingBox2d;
+
 	using Transformd = Transform2d;
 
 	template<class DataType, int ARRAY_SIZE>

SPHINXsys/src/for_2D_build/meshes/base_mesh_supplementary.cpp

@@ -29,21 +29,21 @@ namespace SPH {
 			}, ap);
 	}
 	//=============================================================================================//
-	Vecu BaseMesh::transfer1DtoMeshIndex(const Vecu &number_of_grid_points, size_t i)
+	Vecu BaseMesh::transfer1DtoMeshIndex(const Vecu &number_of_mesh_indexes, size_t i)
 	{
-		size_t row_size = number_of_grid_points[1];
+		size_t row_size = number_of_mesh_indexes[1];
 		size_t column = i / row_size;
 		return Vec2u(column, i - column * row_size);
 	}
 	//=============================================================================================//
-	size_t BaseMesh::transferMeshIndexTo1D(const Vecu &number_of_grid_points, const Vecu &grid_index)
+	size_t BaseMesh::transferMeshIndexTo1D(const Vecu &number_of_mesh_indexes, const Vecu &mesh_index)
 	{
-		return grid_index[0] * number_of_grid_points[1] + grid_index[1];
+		return mesh_index[0] * number_of_mesh_indexes[1] + mesh_index[1];
 	}
     //=============================================================================================//
-    size_t BaseMesh::transferMeshIndexToMortonOrder(const Vecu &grid_index)
+    size_t BaseMesh::transferMeshIndexToMortonOrder(const Vecu &mesh_index)
     {
-        return MortonCode(grid_index[0]) | (MortonCode(grid_index[1]) << 1);
+        return MortonCode(mesh_index[0]) | (MortonCode(mesh_index[1]) << 1);
     }
 }
 //=============================================================================================//

SPHINXsys/src/for_2D_build/meshes/cell_linked_list_supplementary.cpp

@@ -7,6 +7,7 @@
 #include "base_kernel.h"
 #include "base_body.h"
 #include "base_particles.h"
+#include "base_particle_dynamics.h"
 #include "neighbor_relation.h"
 
 namespace SPH

SPHINXsys/src/for_2D_build/meshes/mesh_with_data_packages.hpp

@@ -10,32 +10,32 @@
 #include "mesh_with_data_packages.h"
 
 //=================================================================================================//
-namespace SPH {
+namespace SPH
+{
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<class DataType>
-	DataType BaseDataPackage<PKG_SIZE, ADDRS_SIZE>
-		::probeDataPackage(PackageDataAddress<DataType>& pkg_data_addrs, const Vecd& position)
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <class DataType>
+	DataType BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
+		probeDataPackage(PackageDataAddress<DataType> &pkg_data_addrs, const Vecd &position)
 	{
 		Vecu grid_idx = CellIndexFromPosition(position);
 		Vecd grid_pos = GridPositionFromIndex(grid_idx);
 		Vecd alpha = (position - grid_pos) / grid_spacing_;
 		Vecd beta = Vec2d(1.0) - alpha;
 
-		DataType bilinear
-			= *pkg_data_addrs[grid_idx[0]][grid_idx[1]] * beta[0] * beta[1]
-			+ *pkg_data_addrs[grid_idx[0] + 1][grid_idx[1]] * alpha[0] * beta[1]
-			+ *pkg_data_addrs[grid_idx[0]][grid_idx[1] + 1] * beta[0] * alpha[1]
-			+ *pkg_data_addrs[grid_idx[0] + 1][grid_idx[1] + 1] * alpha[0] * alpha[1];
+		DataType bilinear = *pkg_data_addrs[grid_idx[0]][grid_idx[1]] * beta[0] * beta[1] +
+							*pkg_data_addrs[grid_idx[0] + 1][grid_idx[1]] * alpha[0] * beta[1] +
+							*pkg_data_addrs[grid_idx[0]][grid_idx[1] + 1] * beta[0] * alpha[1] +
+							*pkg_data_addrs[grid_idx[0] + 1][grid_idx[1] + 1] * alpha[0] * alpha[1];
 
-		return  bilinear;
+		return bilinear;
 	}
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<typename InDataType, typename OutDataType>
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <typename InDataType, typename OutDataType>
 	void BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
-		computeGradient(PackageDataAddress<InDataType>& in_pkg_data_addrs,
-			PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt)
+		computeGradient(PackageDataAddress<InDataType> &in_pkg_data_addrs,
+						PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt)
 	{
 		for (int i = 1; i != PKG_SIZE + 1; ++i)
 			for (int j = 1; j != PKG_SIZE + 1; ++j)
@@ -46,11 +46,11 @@ namespace SPH {
 			}
 	}
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<typename InDataType, typename OutDataType>
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <typename InDataType, typename OutDataType>
 	void BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
-		computeNormalizedGradient(PackageDataAddress<InDataType>& in_pkg_data_addrs,
-			PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt)
+		computeNormalizedGradient(PackageDataAddress<InDataType> &in_pkg_data_addrs,
+								  PackageDataAddress<OutDataType> out_pkg_data_addrs, Real dt)
 	{
 		for (int i = 1; i != PKG_SIZE + 1; ++i)
 			for (int j = 1; j != PKG_SIZE + 1; ++j)
@@ -62,11 +62,11 @@ namespace SPH {
 			}
 	}
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<typename DataType>
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <typename DataType>
 	void BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
-		initializePackageDataAddress(PackageData<DataType>& pkg_data,
-			PackageDataAddress<DataType>& pkg_data_addrs)
+		initializePackageDataAddress(PackageData<DataType> &pkg_data,
+									 PackageDataAddress<DataType> &pkg_data_addrs)
 	{
 		for (int i = 0; i != ADDRS_SIZE; ++i)
 			for (int j = 0; j != ADDRS_SIZE; ++j)
@@ -75,10 +75,10 @@ namespace SPH {
 			}
 	}
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<typename DataType>
-	DataType  BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
-		CornerAverage(PackageDataAddress<DataType>& pkg_data_addrs, Veci addrs_index, Veci corner_direction)
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <typename DataType>
+	DataType BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
+		CornerAverage(PackageDataAddress<DataType> &pkg_data_addrs, Veci addrs_index, Veci corner_direction)
 	{
 		DataType average(0);
 		for (int i = 0; i != 2; ++i)
@@ -91,19 +91,19 @@ namespace SPH {
 		return average * 0.25;
 	}
 	//=================================================================================================//
-	template<int PKG_SIZE, int ADDRS_SIZE>
-	template<typename DataType>
+	template <int PKG_SIZE, int ADDRS_SIZE>
+	template <typename DataType>
 	void BaseDataPackage<PKG_SIZE, ADDRS_SIZE>::
-		assignPackageDataAddress(PackageDataAddress<DataType>& pkg_data_addrs, Vecu& addrs_index,
-			PackageData<DataType>& pkg_data, Vecu& data_index)
+		assignPackageDataAddress(PackageDataAddress<DataType> &pkg_data_addrs, Vecu &addrs_index,
+								 PackageData<DataType> &pkg_data, Vecu &data_index)
 	{
 		pkg_data_addrs[addrs_index[0]][addrs_index[1]] = &pkg_data[data_index[0]][data_index[1]];
 	}
 	//=================================================================================================//
-	template<class MeshFieldType, class DataPackageType> 
-	template<typename DataType, typename PackageDataType, PackageDataType DataPackageType:: * MemPtr>
+	template <class MeshFieldType, class DataPackageType>
+	template <typename DataType, typename PackageDataType, PackageDataType DataPackageType::*MemPtr>
 	DataType MeshWithDataPackages<MeshFieldType, DataPackageType>::
-		DataValueFromGlobalIndex(Vecu global_grid_index)
+		DataValueFromGlobalIndex(const Vecu &global_grid_index)
 	{
 		Vecu pkg_index_(0);
 		Vecu local_data_index(0);
@@ -113,54 +113,71 @@ namespace SPH {
 			pkg_index_[n] = cell_index_in_this_direction;
 			local_data_index[n] = global_grid_index[n] - cell_index_in_this_direction * pkg_size_;
 		}
-		PackageDataType& data = data_pkg_addrs_[pkg_index_[0]][pkg_index_[1]]->*MemPtr;
+		PackageDataType &data = data_pkg_addrs_[pkg_index_[0]][pkg_index_[1]]->*MemPtr;
 		return data[local_data_index[0]][local_data_index[1]];
 	}
 	//=================================================================================================//
-	template<class MeshFieldType, class DataPackageType> 
-	void MeshWithDataPackages<MeshFieldType, DataPackageType>::initializePackageAddressesInACell(Vecu cell_index)
+	template <class MeshFieldType, class DataPackageType>
+	void MeshWithDataPackages<MeshFieldType, DataPackageType>::
+		initializePackageAddressesInACell(const Vecu &cell_index)
 	{
 		int i = (int)cell_index[0];
 		int j = (int)cell_index[1];
 
-		DataPackageType* data_pkg = data_pkg_addrs_[i][j];
-		if (data_pkg->is_inner_pkg_) {
+		DataPackageType *data_pkg = data_pkg_addrs_[i][j];
+		if (data_pkg->is_inner_pkg_)
+		{
 			for (int l = 0; l != pkg_addrs_size_; ++l)
-				for (int m = 0; m != pkg_addrs_size_; ++m) {
-					std::pair<int, int>  x_pair = CellShiftAndDataIndex(l);
-					std::pair<int, int>  y_pair = CellShiftAndDataIndex(m);
-					data_pkg->assignAllPackageDataAddress(Vecu(l, m), 
+				for (int m = 0; m != pkg_addrs_size_; ++m)
+				{
+					std::pair<int, int> x_pair = CellShiftAndDataIndex(l);
+					std::pair<int, int> y_pair = CellShiftAndDataIndex(m);
+					data_pkg->assignAllPackageDataAddress(
+						Vecu(l, m),
 						data_pkg_addrs_[i + x_pair.first][j + y_pair.first],
 						Vecu(x_pair.second, y_pair.second));
 				}
 		}
 	}
 	//=================================================================================================//
-	template<class MeshFieldType, class DataPackageType>
+	template <class MeshFieldType, class DataPackageType>
 	void MeshWithDataPackages<MeshFieldType, DataPackageType>::allocateMeshDataMatrix()
 	{
 		Allocate2dArray(data_pkg_addrs_, number_of_cells_);
 	}
 	//=================================================================================================//
-	template<class MeshFieldType, class DataPackageType>
+	template <class MeshFieldType, class DataPackageType>
 	void MeshWithDataPackages<MeshFieldType, DataPackageType>::deleteMeshDataMatrix()
 	{
 		Delete2dArray(data_pkg_addrs_, number_of_cells_);
 	}
 	//=================================================================================================//
-	template<class MeshFieldType, class DataPackageType>
-	template<class DataType, typename PackageDataAddressType, PackageDataAddressType DataPackageType:: * MemPtr>
-	DataType MeshWithDataPackages<MeshFieldType, DataPackageType>::probeMesh(const Vecd& position)
+	template <class MeshFieldType, class DataPackageType>
+	void MeshWithDataPackages<MeshFieldType, DataPackageType>::
+		assignDataPackageAddress(const Vecu &cell_index, DataPackageType *data_pkg)
+	{
+		data_pkg_addrs_[cell_index[0]][cell_index[1]] = data_pkg;
+	}
+	//=================================================================================================//
+	template <class MeshFieldType, class DataPackageType>
+	DataPackageType *MeshWithDataPackages<MeshFieldType, DataPackageType>::
+		DataPackageFromCellIndex(const Vecu &cell_index)
+	{
+		return data_pkg_addrs_[cell_index[0]][cell_index[1]];
+	}
+	//=================================================================================================//
+	template <class MeshFieldType, class DataPackageType>
+	template <class DataType, typename PackageDataAddressType, PackageDataAddressType DataPackageType::*MemPtr>
+	DataType MeshWithDataPackages<MeshFieldType, DataPackageType>::probeMesh(const Vecd &position)
 	{
-		Vecu grid_index =CellIndexFromPosition(position);
+		Vecu grid_index = CellIndexFromPosition(position);
 		size_t i = grid_index[0];
 		size_t j = grid_index[1];
 
-		DataPackageType* data_pkg = data_pkg_addrs_[i][j];
-		PackageDataAddressType& pkg_data_addrs = data_pkg->*MemPtr;
-		return data_pkg->is_inner_pkg_ ?
-			data_pkg->DataPackageType::template probeDataPackage<DataType>(pkg_data_addrs, position)
-			: *pkg_data_addrs[0][0];
+		DataPackageType *data_pkg = data_pkg_addrs_[i][j];
+		PackageDataAddressType &pkg_data_addrs = data_pkg->*MemPtr;
+		return data_pkg->is_inner_pkg_ ? data_pkg->DataPackageType::template probeDataPackage<DataType>(pkg_data_addrs, position)
+									   : *pkg_data_addrs[0][0];
 	}
 	//=================================================================================================//
 }

SPHINXsys/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp

@@ -9,6 +9,7 @@
 #include "base_mesh.h"
 #include "base_body.h"
 #include "base_particles.h"
+#include "base_particle_dynamics.h"
 
 namespace SPH {
 	//=================================================================================================//

SPHINXsys/src/for_2D_build/particles/solid_particles_supplementary.cpp

@@ -4,11 +4,13 @@
  */
 
 #include "solid_particles.h"
+#include "solid_particles_variable.h"
 #include "base_body.h"
 
-namespace SPH {
+namespace SPH
+{
 	//=============================================================================================//
-	void SolidParticles::ParticleTranslationAndRotation(Transformd& transform) 
+	void SolidParticles::ParticleTranslationAndRotation(Transformd &transform)
 	{
 		for (size_t i = 0; i != total_real_particles_; ++i)
 		{
@@ -102,20 +104,50 @@ namespace SPH {
 		return getVonMisesStressFromMatrix(sigma);
 	}
 	//=================================================================================================//
-	Vecd ElasticSolidParticles::displacement(size_t particle_i) //not tested in 2D
+	Real ElasticSolidParticles::von_Mises_strain(size_t particle_i) //not tested in 2D
 	{
-		Vecd disp = pos_n_[particle_i] - pos_0_[particle_i];
-		return disp;
+
+		Mat2d F = F_[particle_i];
+		Mat2d epsilon = 0.5 * (~F * F - Mat2d(1.0)); //calculation of the Green-Lagrange strain tensor
+
+		Real epsilonxx = epsilon(0, 0);
+		Real epsilonyy = epsilon(1, 1);
+		Real epsilonzz = 0; //z-components zero for 2D measures
+		Real epsilonxy = epsilon(0, 1);
+		Real epsilonxz = 0; //z-components zero for 2D measures
+		Real epsilonyz = 0; //z-components zero for 2D measures
+
+		return sqrt((1.0 / 3.0) * (std::pow(epsilonxx - epsilonyy, 2.0) + std::pow(epsilonyy - epsilonzz, 2.0) +
+								   std::pow(epsilonzz - epsilonxx, 2.0)) +
+					2.0 * (std::pow(epsilonxy, 2.0) + std::pow(epsilonyz, 2.0) + std::pow(epsilonxz, 2.0)));
 	}
-	//=================================================================================================//
-	Vecd ElasticSolidParticles::normal(size_t particle_i) //not tested in 2D
+	//=============================================================================================//
+	void VonMisesStress::operator()(size_t index_i, Real dt)
 	{
-		Vecd normal_vec = n_[particle_i];
-		return normal_vec;
+		Real J = rho0_ / rho_n_[index_i];
+		Mat2d F = F_[index_i];
+		Mat2d stress = stress_PK1_[index_i];
+		Mat2d sigma = (stress * ~F) / J;
+
+		Real sigmaxx = sigma(0, 0);
+		Real sigmayy = sigma(1, 1);
+		Real sigmaxy = sigma(0, 1);
+
+		derived_variable_[index_i] =
+			sqrt(sigmaxx * sigmaxx + sigmayy * sigmayy - sigmaxx * sigmayy + 3.0 * sigmaxy * sigmaxy);
 	}
-	//=================================================================================================//
-	Real ElasticSolidParticles::von_Mises_strain(size_t particle_i) //not tested in 2D
+	//=============================================================================================//
+	void VonMisesStrain::operator()(size_t index_i, Real dt)
 	{
-		return von_Mises_strain_static(particle_i);
+		Mat2d F = F_[index_i];
+		Mat2d epsilon = 0.5 * (~F * F - Mat2d(1.0)); //calculation of the Green-Lagrange strain tensor
+
+		Real epsilonxx = epsilon(0, 0);
+		Real epsilonyy = epsilon(1, 1);
+		Real epsilonxy = epsilon(0, 1);
+
+		derived_variable_[index_i] = sqrt(0.5 * ((epsilonxx - epsilonyy) * (epsilonxx - epsilonyy) +
+												 epsilonyy * epsilonyy + epsilonxx * epsilonxx +
+												 2.0 * epsilonxy * epsilonxy));
 	}
 }

SPHINXsys/src/for_3D_build/common/data_type.h

@@ -16,6 +16,8 @@ namespace SPH {
 	using SymMatd = SymMat3d;
 	using AngularVecd = Vec3d;
 
+	using  BoundingBox = BoundingBox3d;
+
 	using Transformd = Transform3d;
 
 	template<class DataType, int ARRAY_SIZE>

SPHINXsys/src/for_3D_build/geometries/all_geometries.h

@@ -6,6 +6,7 @@
 This is the header file that user code should include to pick up all
 geometry classes used in SPHinXsys. **/
 
+#include "geometric_shape.h"
 #include "triangle_mesh_shape.h"
 #include "image_shape.h"
 #include "level_set_shape.h"