Merge branch 'Xiangyu-Hu:master' into master

.github/workflows/ci.yml

@@ -123,7 +123,7 @@ jobs:
       VCPKG_DEFAULT_TRIPLET: x64-windows
 
     steps:
-      - uses: actions/checkout@v3 # Checks-out the repository under ${{github.workspace}}
+      - uses: actions/checkout@v3.5.2 # Checks-out the repository under ${{github.workspace}}
 
       - name: Update ccache and ninja # For correct caching with ccache on Windows
         shell: bash

CMakeLists.txt

@@ -10,6 +10,7 @@ option(SPHINXSYS_2D "Build sphinxsys_2d library" ON)
 option(SPHINXSYS_3D "Build sphinxsys_3d library" ON)
 option(SPHINXSYS_BUILD_TESTS "Build tests" ON)
 option(SPHINXSYS_DEVELOPER_MODE "Developer mode has more flags active for code quality" ON)
+option(SPHINXSYS_USE_FLOAT "Build using float or not" OFF)
 
 # ------ Global properties (Some cannot be set on INTERFACE targets)
 set(CMAKE_VERBOSE_MAKEFILE OFF CACHE BOOL "Enable verbose compilation commands for Makefile and Ninja" FORCE) # Extra fluff needed for Ninja: https://github.com/ninja-build/ninja/issues/900
@@ -78,6 +79,11 @@ set(Boost_NO_WARN_NEW_VERSIONS TRUE) # In case your CMake version is older than
 find_package(Boost REQUIRED)
 find_path(BOOST_INCLUDE_DIR boost/geometry) # Header-only Boost libraries are not components
 
+if(SPHINXSYS_USE_FLOAT)
+    target_compile_definitions(sphinxsys_core INTERFACE USE_FLOAT)
+    message("-- Use float as default")
+endif()
+
 if(NOT BOOST_INCLUDE_DIR)
     message(FATAL_ERROR "Please install Boost.Geometry library")
 endif()

src/for_2D_build/bodies/solid_body_supplementary.cpp

@@ -40,7 +40,7 @@ namespace SPH
 		Iz /= body_part_volume;
 
 		body_part_mass_properties_ = mass_properties_ptr_keeper_.createPtr<SimTK::MassProperties>(
-			body_part_volume * solid_body_density_, SimTK::Vec3(0), SimTK::UnitInertia(Ix, Iy, Iz));
+			body_part_volume * solid_body_density_, SimTKVec3(0), SimTK::UnitInertia(Ix, Iy, Iz));
 	}
 	//=================================================================================================//
 }

src/for_2D_build/common/data_type.h

@@ -40,7 +40,7 @@ namespace SPH {
 	using AlignedBox = AlignedBox2d;
 	using AngularVecd = Real;
 	using Transformd = Transform2d;
-	using SimTKVecd = SimTK::Vec2;
+	using SimTKVecd = SimTKVec2;
 
 	using  BoundingBox = BaseBoundingBox<Vec2d>;
 

src/for_2D_build/geometries/multi_polygon_shape.cpp

@@ -94,7 +94,7 @@ namespace SPH
 		strategy::buffer::end_round end_strategy;
 		strategy::buffer::side_straight side_strategy;
 		strategy::buffer::point_circle circle_strategy(buffer_res);
-		strategy::buffer::distance_symmetric<double> circle_dist_strategy(buffer_radius);
+		strategy::buffer::distance_symmetric<Real> circle_dist_strategy(buffer_radius);
 
 		// Create the buffer of a multi point
 		model::d2::point_xy<Real> circle_center_pnt;
@@ -154,7 +154,7 @@ namespace SPH
 		std::fstream dataFile(file_path_name);
 		Vecd temp_point;
 		std::vector<Vecd> coordinates;
-		double temp1 = 0.0, temp2 = 0.0;
+		Real temp1 = 0.0, temp2 = 0.0;
 		if (dataFile.fail())
 		{
 			std::cout << "File can not open.\n"

src/for_3D_build/bodies/solid_body_supplementary.cpp

@@ -42,9 +42,9 @@ namespace SPH
 		inertia_products /= body_part_volume;
 
 		body_part_mass_properties_ = mass_properties_ptr_keeper_.createPtr<SimTK::MassProperties>(
-			body_part_volume * solid_body_density_, SimTK::Vec3(0), 
-			SimTK::UnitInertia(SimTK::Vec3(inertia_moments[0],inertia_moments[1],inertia_moments[2]), 
-							   SimTK::Vec3(inertia_products[0],inertia_products[1],inertia_products[2])));
+			body_part_volume * solid_body_density_, SimTKVec3(0), 
+			SimTK::UnitInertia(SimTKVec3(inertia_moments[0],inertia_moments[1],inertia_moments[2]), 
+							   SimTKVec3(inertia_products[0],inertia_products[1],inertia_products[2])));
 	}
 	//=================================================================================================//
 }

src/for_3D_build/common/data_type.h

@@ -39,7 +39,7 @@ namespace SPH {
 	using AlignedBox = AlignedBox3d;
 	using AngularVecd = Vec3d;
 	using Transformd = Transform3d;
-	using SimTKVecd = SimTK::Vec3;
+	using SimTKVecd = SimTKVec3;
 
 	using  BoundingBox = BaseBoundingBox<Vec3d>;
 

src/for_3D_build/geometries/geometric_shape.cpp

@@ -7,7 +7,7 @@ namespace SPH
     {
         SimTK::UnitVec3 normal;
         bool inside = false;
-        contact_geometry_->findNearestPoint(SimTK::Vec3(probe_point[0], probe_point[1], probe_point[2]), inside, normal);
+        contact_geometry_->findNearestPoint(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]), inside, normal);
 
         return inside;
     }
@@ -16,7 +16,7 @@ namespace SPH
     {
         SimTK::UnitVec3 normal;
         bool inside = false;
-        SimTK::Vec3 out_pnt = contact_geometry_->findNearestPoint(SimTK::Vec3(probe_point[0], probe_point[1], probe_point[2]), inside, normal);
+        SimTKVec3 out_pnt = contact_geometry_->findNearestPoint(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]), inside, normal);
 
         return Vecd(out_pnt[0], out_pnt[1], out_pnt[2]);
     }
@@ -30,13 +30,13 @@ namespace SPH
     //=================================================================================================//
     bool GeometricShapeBox::checkContain(const Vec3d &probe_point, bool BOUNDARY_INCLUDED)
     {
-        return brick_.getGeoBox().containsPoint(SimTK::Vec3(probe_point[0], probe_point[1], probe_point[2]));
+        return brick_.getGeoBox().containsPoint(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]));
     }
     //=================================================================================================//
     Vec3d GeometricShapeBox::findClosestPoint(const Vec3d &probe_point)
     {
         bool inside = false;
-        SimTK::Vec3 out_pnt = brick_.getGeoBox().findClosestPointOnSurface(SimTK::Vec3(probe_point[0], probe_point[1], probe_point[2]), inside);
+        SimTKVec3 out_pnt = brick_.getGeoBox().findClosestPointOnSurface(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]), inside);
 
         return Vecd(out_pnt[0], out_pnt[1], out_pnt[2]);
     }

src/for_3D_build/geometries/image_shape.cpp

@@ -45,7 +45,7 @@ namespace SPH
 		ImageShapeSphere(Real radius, Vecd spacings, Vecd center, const std::string &shape_name)
 		: ImageShape(shape_name)
 	{
-		double extend = 1.5;
+		Real extend = 1.5;
 		int length = int(std::ceil(2.0 * extend * radius));
 		Arrayi NxNyNz(length, length, length);
 		image_.reset(new ImageMHD<float, 3>(radius, NxNyNz, spacings));

src/for_3D_build/geometries/triangle_mesh_shape.cpp

@@ -20,11 +20,11 @@ namespace SPH
 	//=================================================================================================//
 	bool TriangleMeshShape::checkContain(const Vec3d &probe_point, bool BOUNDARY_INCLUDED)
 	{
-		SimTK::Vec2 uv_coordinate;
+		SimTKVec2 uv_coordinate;
 		bool inside = false; // note that direct prediction is not reliable sometime.
 		int face_id;
 
-		SimTK::Vec3 closest_pnt = triangle_mesh_->findNearestPoint(EigenToSimTK(probe_point), inside, face_id, uv_coordinate);
+		SimTKVec3 closest_pnt = triangle_mesh_->findNearestPoint(EigenToSimTK(probe_point), inside, face_id, uv_coordinate);
 		Vec3d from_face_to_pnt = probe_point - SimTKToEigen(closest_pnt);
 		Real distance_to_pnt = from_face_to_pnt.norm();
 		Vec3d direction_to_pnt = from_face_to_pnt / (distance_to_pnt + TinyReal);
@@ -57,8 +57,8 @@ namespace SPH
 	{
 		bool inside = false;
 		int face_id;
-		SimTK::Vec2 norm;
-		SimTK::Vec3 closest_pnt = triangle_mesh_->findNearestPoint(SimTK::Vec3(probe_point[0], probe_point[1], probe_point[2]), inside, face_id, norm);
+		SimTKVec2 norm;
+		SimTKVec3 closest_pnt = triangle_mesh_->findNearestPoint(SimTKVec3(probe_point[0], probe_point[1], probe_point[2]), inside, face_id, norm);
 		if (face_id < 0 && face_id > triangle_mesh_->getNumFaces())
 		{
 			std::cout << "\n Error the nearest point is not valid" << std::endl;
@@ -72,8 +72,8 @@ namespace SPH
 	{
 		int number_of_vertices = triangle_mesh_->getNumVertices();
 		// initial reference values
-		Vec3d lower_bound = SimTKToEigen(SimTK::Vec3(Infinity));
-		Vec3d upper_bound = SimTKToEigen(SimTK::Vec3(-Infinity));
+		Vec3d lower_bound = SimTKToEigen(SimTKVec3(Infinity));
+		Vec3d upper_bound = SimTKToEigen(SimTKVec3(-Infinity));
 		for (int i = 0; i != number_of_vertices; ++i)
 		{
 			Vec3d vertex_position = SimTKToEigen(triangle_mesh_->getVertexPosition(i));
@@ -99,7 +99,7 @@ namespace SPH
 		SimTK::PolygonalMesh polymesh;
 		polymesh.loadStlFile(filepathname);
 		polymesh.scaleMesh(scale_factor);
-		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTK::Vec3(translation[0], translation[1], translation[2])));
+		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 	}
 	//=================================================================================================//
 	TriangleMeshShapeSTL::TriangleMeshShapeSTL(const std::string &filepathname, Mat3d rotation,
@@ -116,10 +116,10 @@ namespace SPH
 		polymesh.loadStlFile(filepathname);
 
 		polymesh.scaleMesh(scale_factor);
-		SimTK::Transform_<Real> transform(SimTK::Rotation_<Real>(SimTK::Mat33(rotation(0, 0), rotation(0, 1), rotation(0, 2),
-																			  rotation(1, 0), rotation(1, 1), rotation(1, 2),
-																			  rotation(2, 0), rotation(2, 1), rotation(2, 2))),
-										  SimTK::Vec3(translation[0], translation[1], translation[2]));
+		SimTK::Transform_<double> transform(SimTK::Rotation_<double>(SimTKMat33((double)rotation(0, 0), (double)rotation(0, 1), (double)rotation(0, 2),
+																			  (double)rotation(1, 0), (double)rotation(1, 1), (double)rotation(1, 2),
+																			  (double)rotation(2, 0), (double)rotation(2, 1), (double)rotation(2, 2))),
+										  SimTKVec3((double)translation[0], (double)translation[1], (double)translation[2]));
 		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(transform));
 	}
 //=================================================================================================//
@@ -131,16 +131,16 @@ namespace SPH
 		SimTK::PolygonalMesh polymesh;
 		polymesh.loadStlBuffer(buffer);
 		polymesh.scaleMesh(scale_factor);
-		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTK::Vec3(translation[0], translation[1], translation[2])));
+		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 	}
 #endif
 	//=================================================================================================//
 	TriangleMeshShapeBrick::TriangleMeshShapeBrick(Vecd halfsize, int resolution, Vecd translation,
 												   const std::string &shape_name)
 		: TriangleMeshShape(shape_name)
 	{
-		SimTK::PolygonalMesh polymesh = SimTK::PolygonalMesh::createBrickMesh(SimTK::Vec3(halfsize[0], halfsize[1], halfsize[2]), resolution);
-		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTK::Vec3(translation[0], translation[1], translation[2])));
+		SimTK::PolygonalMesh polymesh = SimTK::PolygonalMesh::createBrickMesh(SimTKVec3(halfsize[0], halfsize[1], halfsize[2]), resolution);
+		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 	}
 	//=================================================================================================//
 	TriangleMeshShapeBrick::TriangleMeshShapeBrick(const TriangleMeshShapeBrick::ShapeParameters &shape_parameters,
@@ -153,7 +153,7 @@ namespace SPH
 		: TriangleMeshShape(shape_name)
 	{
 		SimTK::PolygonalMesh polymesh = SimTK::PolygonalMesh::createSphereMesh(radius, resolution);
-		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTK::Vec3(translation[0], translation[1], translation[2])));
+		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 	}
 	//=================================================================================================//
 	TriangleMeshShapeCylinder::TriangleMeshShapeCylinder(SimTK::UnitVec3 axis, Real radius, Real halflength, int resolution,
@@ -162,7 +162,7 @@ namespace SPH
 	{
 		SimTK::PolygonalMesh polymesh =
 			SimTK::PolygonalMesh::createCylinderMesh(axis, radius, halflength, resolution);
-		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTK::Vec3(translation[0], translation[1], translation[2])));
+		triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(SimTKVec3(translation[0], translation[1], translation[2])));
 	}
 	//=================================================================================================//
 }

src/shared/adaptations/adaptation.cpp

@@ -130,7 +130,7 @@ namespace SPH
 	{
 		SPHAdaptation::registerAdaptationVariables(base_particles);
 
-		base_particles.registerVariable(h_ratio_, "SmoothingLengthRatio", 1.0);
+		base_particles.registerVariable(h_ratio_, "SmoothingLengthRatio", Real(1.0));
 		base_particles.registerSortableVariable<Real>("SmoothingLengthRatio");
 		base_particles.addVariableToReload<Real>("SmoothingLengthRatio");
 	}
@@ -210,12 +210,12 @@ namespace SPH
 	//=================================================================================================//
 	Vec2d ParticleSplitAndMerge::splittingPattern(Vec2d pos, Real particle_spacing, Real delta)
 	{
-		return {pos[0] + 0.5 * particle_spacing * cos(delta), pos[1] + 0.5 * particle_spacing * sin(delta)};
+		return {pos[0] + Real(0.5) * particle_spacing * cos(delta), pos[1] + Real(0.5) * particle_spacing * sin(delta)};
 	}
 	//=================================================================================================//
 	Vec3d ParticleSplitAndMerge::splittingPattern(Vec3d pos, Real particle_spacing, Real delta)
 	{
-		return {pos[0] + 0.5 * particle_spacing * cos(delta), pos[1] + 0.5 * particle_spacing * sin(delta), pos[2]};
+		return {pos[0] + Real(0.5) * particle_spacing * cos(delta), pos[1] + Real(0.5) * particle_spacing * sin(delta), pos[2]};
 	}
 	//=================================================================================================//
 }

src/shared/common/base_data_type.h

@@ -69,7 +69,14 @@ namespace SPH
 	 * 			Using -mtune=native produces code optimized for the local machine under the constraints of the selected instruction set.
 	 */
 
-	using Real = double;
+#if USE_FLOAT
+		using Real = float;
+		using EigMat = Eigen::MatrixXf;
+#else
+		using Real = double;
+		using EigMat = Eigen::MatrixXd;
+#endif
+
 	/** Vector with integers. */
 	using Array2i = Eigen::Array<int, 2, 1>;
 	using Array3i = Eigen::Array<int, 3, 1>;
@@ -82,6 +89,14 @@ namespace SPH
 	/** AlignedBox */
 	using AlignedBox2d = Eigen::AlignedBox<Real, 2>;
 	using AlignedBox3d = Eigen::AlignedBox<Real, 3>;
+	/** SimTK vector and matrix */
+	/**
+	 * Note that SimTK vector and matrix using double precision, whne ENABLE_SINGLE, will transfer to double precision. 
+	*/
+	using SimTKVec2 = SimTK::Vec2;
+	using SimTKVec3 = SimTK::Vec3;
+	using SimTKMat22 = SimTK::Mat22;
+	using SimTKMat33 = SimTK::Mat33;
 	/** Unified initialize to zero for all data type. */
 	/**
 	 * NOTE: Eigen::Matrix<> constexpr constructor?
@@ -142,20 +157,20 @@ namespace SPH
 	};
 	/** Useful float point constants. */
 	constexpr size_t MaxSize_t = std::numeric_limits<size_t>::max();
-	constexpr double MinRealNumber = std::numeric_limits<double>::min();
-	constexpr double MaxRealNumber = std::numeric_limits<double>::max();
+	constexpr Real MinRealNumber = std::numeric_limits<Real>::min();
+	constexpr Real MaxRealNumber = std::numeric_limits<Real>::max();
 	/** Verbal boolean for positive and negative axis directions. */
 	const int xAxis = 0;
 	const int yAxis = 1;
 	const int zAxis = 2;
 	const bool positiveDirection = true;
 	const bool negativeDirection = false;
 	/** Constant parameters. */
-	const Real Pi = Real(M_PI);
-	const Real Eps = 2.22045e-16;
-	const Real TinyReal = 2.71051e-20;
-	const Real Infinity = std::numeric_limits<double>::max();
-	const Real SqrtEps = 1.0e-8;
+	constexpr Real Pi = Real(M_PI);
+	constexpr Real Eps = std::numeric_limits<Real>::epsilon();
+	constexpr Real TinyReal = Real(2.71051e-20);
+	constexpr Real Infinity = std::numeric_limits<Real>::max();
+	constexpr Real SqrtEps = Real(1.0e-8);
 }
 
 #endif // BASE_DATA_TYPE_H

src/shared/common/image_mhd.hpp

@@ -189,7 +189,7 @@ namespace SPH {
 				for (int x = 0; x < width_; x++)
 				{
 					int index = z * width_*height_ + y * width_ + x;
-					double distance = (Vecd(x, y, z) - center).norm() - radius;
+					Real distance = (Vecd(x, y, z) - center).norm() - radius;
 					if (distance < min_value_) min_value_ = distance;
 					if (distance > max_value_) max_value_ = distance;
 					data_[index] = float(distance);
@@ -256,7 +256,7 @@ namespace SPH {
 		int y = (i%sliceSize) / width;
 		int x = (i%sliceSize) % width;
 
-		double gradx = 0.0; double grady = 0.0; double gradz = 0.0;
+		Real gradx = 0.0; Real grady = 0.0; Real gradz = 0.0;
 		//- cds (if inner cell)
 		//- otherwise back/forward scheme
 		if (x == 0)
@@ -367,20 +367,20 @@ namespace SPH {
 		Array3i this_cell = Array3i::Zero();
 		std::vector<int> neighbors = findNeighbors(probe_point, this_cell);
 		Vec3d n_sum = Vecd::Zero();
-		double weight_sum = 0.0;
-		double d_sum = 0.0;
+		Real weight_sum = 0.0;
+		Real d_sum = 0.0;
 		for (const int& i : neighbors)
 		{
 			// checkIndexBound(i);
 			Vec3d nCj = computeNormalAtCell(i);
-			double dCj = float(getValueAtCell(i));
-			double weight_Cj = 1.0 / (fabs(dCj) + Eps);
+			Real dCj = float(getValueAtCell(i));
+			Real weight_Cj = 1.0 / (fabs(dCj) + Eps);
 			n_sum = n_sum + weight_Cj * nCj;
 			weight_sum = weight_sum + weight_Cj;
 			d_sum = d_sum + dCj;
 		}
 		Vec3d n = n_sum / (weight_sum + Eps);
-		double d = d_sum / (weight_sum + Eps);
+		Real d = d_sum / (weight_sum + Eps);
 
 		Vec3d p_image = Vec3d(this_cell[0], this_cell[1], this_cell[2]) + n.normalized() * d;
 		Vec3d p = convertToPhysicalSpace(p_image);
@@ -418,15 +418,15 @@ namespace SPH {
 	{
 		Array3i this_cell;
 		std::vector<int> neighbors = findNeighbors(probe_point, this_cell);
-		double weight_sum = 0.0;
-		double d_sum = 0.0;
+		Real weight_sum = 0.0;
+		Real d_sum = 0.0;
 		if (neighbors.size() > 0)
 		{
 			for (const int& i : neighbors)
 			{
 				// checkIndexBound(i);
-				double dCj = float(getValueAtCell(i));
-				double weight_Cj = 1.0 / (fabs(dCj) + Eps);
+				Real dCj = float(getValueAtCell(i));
+				Real weight_Cj = 1.0 / (fabs(dCj) + Eps);
 				weight_sum = weight_sum + weight_Cj;
 				d_sum = d_sum + dCj;
 			}
@@ -445,16 +445,16 @@ namespace SPH {
 		Array3i this_cell = Array3i::Zero();
 		std::vector<int> neighbors = findNeighbors(probe_point, this_cell);
 		Vec3d n_sum = Vecd::Zero();
-		double weight_sum = 0.0;
-		double d_sum = 0.0;
+		Real weight_sum = 0.0;
+		Real d_sum = 0.0;
 		if (neighbors.size() > 0)
 		{
 			for (const int& i : neighbors)
 			{
 				// checkIndexBound(i);
 				Vec3d nCj = computeNormalAtCell(i);
-				double dCj = float(getValueAtCell(i));
-				double weight_Cj = 1.0 / (fabs(dCj) + Eps);
+				Real dCj = float(getValueAtCell(i));
+				Real weight_Cj = 1.0 / (fabs(dCj) + Eps);
 				n_sum = n_sum + weight_Cj * nCj;
 				weight_sum = weight_sum + weight_Cj;
 				d_sum = d_sum + dCj;