From 97f6409f40331cbac2edc99c319a6c1fff62740f Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Sun, 28 May 2023 18:57:28 +0200 Subject: [PATCH 1/8] Switch from single to double precision. First Step: test with Real = double by doing the following steps. 1, double -> real " using Real = double; using SimTKVec2 = SimTK::Vec2; using SimTKVec3 = SimTK::Vec3; using SimTKMat22 = SimTK::Mat22; using SimTKMat33 = SimTK::Mat33; using EigMat = Eigen::MatrixXd; " (moved from ZhenxiZhao's implementation) Then, 2, SimTK::Vec"d" -> SimTKVec"d" 3, DoubleVec -> BiVector 4, TripleVec -> TriVector 5, CalculateDoubleDotProduct -> CalculateBiDotProduct (clean the variable name with "double", and change all double to Real in code) 6, double -> Real 7, discards changes in polar_docompostion.h, TriangleMeshDistance.h --- CMakeLists.txt | 5 ++ .../bodies/solid_body_supplementary.cpp | 2 +- src/for_2D_build/common/data_type.h | 2 +- .../geometries/multi_polygon_shape.cpp | 4 +- .../bodies/solid_body_supplementary.cpp | 6 +- src/for_3D_build/common/data_type.h | 2 +- .../geometries/geometric_shape.cpp | 8 +- src/for_3D_build/geometries/image_shape.cpp | 2 +- .../geometries/triangle_mesh_shape.cpp | 28 +++---- src/shared/common/base_data_type.h | 31 +++++-- src/shared/common/image_mhd.hpp | 30 +++---- src/shared/common/large_data_containers.h | 4 +- src/shared/common/scalar_functions.cpp | 84 +++++++++---------- src/shared/common/scalar_functions.h | 8 +- src/shared/common/vector_functions.cpp | 36 ++++---- src/shared/common/vector_functions.h | 22 ++--- src/shared/geometries/level_set.h | 2 +- src/shared/materials/elastic_solid.cpp | 4 +- src/shared/meshes/base_mesh.h | 4 +- .../particle_dynamics_dissipation.hpp | 2 +- .../general_dynamics/general_dynamics.cpp | 2 +- .../solid_dynamics/constraint_dynamics.h | 20 ++--- .../particle_generator_lattice.cpp | 2 +- .../dynamic_time_warping_method.h | 2 +- .../dynamic_time_warping_method.hpp | 4 +- .../ensemble_averaged_method.h | 22 ++--- .../ensemble_averaged_method.hpp | 34 ++++---- .../regression_test/regression_test_base.h | 8 +- .../regression_test/regression_test_base.hpp | 12 +-- .../regression_test/time_averaged_method.h | 24 +++--- .../regression_test/time_averaged_method.hpp | 24 +++--- src/shared/simbody_sphinxsys/state_engine.cpp | 14 ++-- src/shared/simbody_sphinxsys/xml_engine.h | 4 +- tests/2d_examples/test_2d_owsc/owsc.cpp | 12 +-- tests/2d_examples/test_2d_owsc/owsc.h | 4 +- .../src/shell_beam_collision.cpp | 4 +- .../src/tethered_dead_fish_in_flow.cpp | 10 +-- .../heart_volume_change.cpp | 6 +- .../heart_volume_change.h | 8 +- .../src/muscle_solid_contact.cpp | 4 +- .../test_vector_functions.cpp | 10 +-- tests/user_examples/test_2d_stfb/stfb.cpp | 10 +-- tests/user_examples/test_2d_stfb/stfb.h | 4 +- tests/user_examples/test_3d_stfb/stfb.cpp | 10 +-- tests/user_examples/test_3d_stfb/stfb.h | 4 +- .../bernoulli_beam/bernoulli_beam.cpp | 4 +- .../bernoulli_beam/bernoulli_beam.h | 16 ++-- 47 files changed, 292 insertions(+), 272 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 01d0c4a727..98851bfd9c 100644 --- a/CMakeLists.txt +++ b/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_SINGLE "Build using single precision 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,10 @@ 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_SINGLE) + add_definitions(-D ENABLE_SINGLE) +endif() + if(NOT BOOST_INCLUDE_DIR) message(FATAL_ERROR "Please install Boost.Geometry library") endif() diff --git a/src/for_2D_build/bodies/solid_body_supplementary.cpp b/src/for_2D_build/bodies/solid_body_supplementary.cpp index f907980153..fd3303aa6b 100644 --- a/src/for_2D_build/bodies/solid_body_supplementary.cpp +++ b/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( - 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)); } //=================================================================================================// } diff --git a/src/for_2D_build/common/data_type.h b/src/for_2D_build/common/data_type.h index ed226a2c1d..a4d4933d4a 100644 --- a/src/for_2D_build/common/data_type.h +++ b/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; diff --git a/src/for_2D_build/geometries/multi_polygon_shape.cpp b/src/for_2D_build/geometries/multi_polygon_shape.cpp index bf29cdf082..dc30e9994f 100644 --- a/src/for_2D_build/geometries/multi_polygon_shape.cpp +++ b/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 circle_dist_strategy(buffer_radius); + strategy::buffer::distance_symmetric circle_dist_strategy(buffer_radius); // Create the buffer of a multi point model::d2::point_xy circle_center_pnt; @@ -154,7 +154,7 @@ namespace SPH std::fstream dataFile(file_path_name); Vecd temp_point; std::vector 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" diff --git a/src/for_3D_build/bodies/solid_body_supplementary.cpp b/src/for_3D_build/bodies/solid_body_supplementary.cpp index 7adb987351..f30671af64 100644 --- a/src/for_3D_build/bodies/solid_body_supplementary.cpp +++ b/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( - 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]))); } //=================================================================================================// } diff --git a/src/for_3D_build/common/data_type.h b/src/for_3D_build/common/data_type.h index 4ec827ee06..f0d4b7d62f 100644 --- a/src/for_3D_build/common/data_type.h +++ b/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; diff --git a/src/for_3D_build/geometries/geometric_shape.cpp b/src/for_3D_build/geometries/geometric_shape.cpp index 92491cea41..59ad22b26a 100644 --- a/src/for_3D_build/geometries/geometric_shape.cpp +++ b/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]); } diff --git a/src/for_3D_build/geometries/image_shape.cpp b/src/for_3D_build/geometries/image_shape.cpp index 0c2bf1899b..9ad8364161 100644 --- a/src/for_3D_build/geometries/image_shape.cpp +++ b/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(radius, NxNyNz, spacings)); diff --git a/src/for_3D_build/geometries/triangle_mesh_shape.cpp b/src/for_3D_build/geometries/triangle_mesh_shape.cpp index a1b2e291c1..6be47521b8 100644 --- a/src/for_3D_build/geometries/triangle_mesh_shape.cpp +++ b/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_ transform(SimTK::Rotation_(SimTK::Mat33(rotation(0, 0), rotation(0, 1), rotation(0, 2), + SimTK::Transform_ transform(SimTK::Rotation_(SimTKMat33(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])); + SimTKVec3(translation[0], translation[1], translation[2])); triangle_mesh_ = generateTriangleMesh(polymesh.transformMesh(transform)); } //=================================================================================================// @@ -131,7 +131,7 @@ 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 //=================================================================================================// @@ -139,8 +139,8 @@ namespace SPH 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]))); } //=================================================================================================// } diff --git a/src/shared/common/base_data_type.h b/src/shared/common/base_data_type.h index 62b4a48154..6581041639 100644 --- a/src/shared/common/base_data_type.h +++ b/src/shared/common/base_data_type.h @@ -69,7 +69,22 @@ namespace SPH * Using -mtune=native produces code optimized for the local machine under the constraints of the selected instruction set. */ - using Real = double; +#if ENABLE_SINGLE + using Real = float; + using SimTKVec2 = SimTK::fVec2; + using SimTKVec3 = SimTK::fVec3; + using SimTKMat22 = SimTK::fMat22; + using SimTKMat33 = SimTK::fMat33; + using EigMat = Eigen::MatrixXf; +#else + using Real = double; + using SimTKVec2 = SimTK::Vec2; + using SimTKVec3 = SimTK::Vec3; + using SimTKMat22 = SimTK::Mat22; + using SimTKMat33 = SimTK::Mat33; + using EigMat = Eigen::MatrixXd; +#endif + /** Vector with integers. */ using Array2i = Eigen::Array; using Array3i = Eigen::Array; @@ -142,8 +157,8 @@ namespace SPH }; /** Useful float point constants. */ constexpr size_t MaxSize_t = std::numeric_limits::max(); - constexpr double MinRealNumber = std::numeric_limits::min(); - constexpr double MaxRealNumber = std::numeric_limits::max(); + constexpr Real MinRealNumber = std::numeric_limits::min(); + constexpr Real MaxRealNumber = std::numeric_limits::max(); /** Verbal boolean for positive and negative axis directions. */ const int xAxis = 0; const int yAxis = 1; @@ -151,11 +166,11 @@ namespace SPH 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::max(); - const Real SqrtEps = 1.0e-8; + constexpr Real Pi = Real(M_PI); + constexpr Real Eps = std::numeric_limits::epsilon(); + constexpr Real TinyReal = pow(Eps, 1.25); + constexpr Real Infinity = std::numeric_limits::max(); + constexpr Real SqrtEps = sqrt(Eps); } #endif // BASE_DATA_TYPE_H diff --git a/src/shared/common/image_mhd.hpp b/src/shared/common/image_mhd.hpp index 9e204878de..e693ddf058 100644 --- a/src/shared/common/image_mhd.hpp +++ b/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 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 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 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; diff --git a/src/shared/common/large_data_containers.h b/src/shared/common/large_data_containers.h index 89e081f196..5bc5a272c4 100644 --- a/src/shared/common/large_data_containers.h +++ b/src/shared/common/large_data_containers.h @@ -62,10 +62,10 @@ namespace SPH { using StdVec = std::vector; template - using DoubleVec = std::vector>; + using BiVector = std::vector>; template - using TripleVec = std::vector>>; + using TriVector = std::vector>>; } #endif //LARGE_DATA_CONTAINER_H diff --git a/src/shared/common/scalar_functions.cpp b/src/shared/common/scalar_functions.cpp index 9a03b8708f..baef4f0147 100644 --- a/src/shared/common/scalar_functions.cpp +++ b/src/shared/common/scalar_functions.cpp @@ -8,60 +8,60 @@ namespace SPH return SecondAxis(SecondAxis(first_axis)); } //=================================================================================================// - double getLeftStateInWeno(double v_1, double v_2, double v_3, double v_4) + Real getLeftStateInWeno(Real v_1, Real v_2, Real v_3, Real v_4) { - double v1 = v_1; - double v2 = v_2; - double v3 = v_3; - double v4 = v_4; + Real v1 = v_1; + Real v2 = v_2; + Real v3 = v_3; + Real v4 = v_4; - double f1 = 0.5 * v2 + 0.5 * v3; - double f2 = (-0.5) * v1 + 1.5 * v2; - double f3 = v2 / 3.0 + 5.0 * v3 / 6.0 - v4 / 6.0; + Real f1 = 0.5 * v2 + 0.5 * v3; + Real f2 = (-0.5) * v1 + 1.5 * v2; + Real f3 = v2 / 3.0 + 5.0 * v3 / 6.0 - v4 / 6.0; - double epsilon = 1.0e-6; - double s1 = pow(v2 - v3, 2) + epsilon; - double s2 = pow(v2 - v1, 2) + epsilon; - double s3 = pow(3.0 * v2 - 4.0 * v3 + v4, 2) / 4.0 + 13.0 * pow(v2 - 2.0 * v3 + v4, 2) / 12.0 + epsilon; - double s12 = 13.0 * pow(v1 - 2.0 * v2 + v3, 2) / 12.0 + pow(v1 - v3, 2) / 4.0 + epsilon; - double tau_4 = (v1 * (547.0 * v1 - 2522.0 * v2 + 1922.0 * v3 - 494.0 * v4) + v2 * (3443.0 * v2 - 5966.0 * v3 + 1602.0 * v4) + v3 * (2843.0 * v3 - 1642.0 * v4) + 267.0 * v4 * v4) / 240.0; + Real epsilon = 1.0e-6; + Real s1 = pow(v2 - v3, 2) + epsilon; + Real s2 = pow(v2 - v1, 2) + epsilon; + Real s3 = pow(3.0 * v2 - 4.0 * v3 + v4, 2) / 4.0 + 13.0 * pow(v2 - 2.0 * v3 + v4, 2) / 12.0 + epsilon; + Real s12 = 13.0 * pow(v1 - 2.0 * v2 + v3, 2) / 12.0 + pow(v1 - v3, 2) / 4.0 + epsilon; + Real tau_4 = (v1 * (547.0 * v1 - 2522.0 * v2 + 1922.0 * v3 - 494.0 * v4) + v2 * (3443.0 * v2 - 5966.0 * v3 + 1602.0 * v4) + v3 * (2843.0 * v3 - 1642.0 * v4) + 267.0 * v4 * v4) / 240.0; - double alpha_1 = (1.0 + (tau_4 / s1) * (tau_4 / s12)) / 3.0; - double alpha_2 = (1.0 + (tau_4 / s2) * (tau_4 / s12)) / 6.0; - double alpha_3 = (1.0 + tau_4 / s3) / 2.0; - double w_1 = alpha_1 / (alpha_1 + alpha_2 + alpha_3); - double w_2 = alpha_2 / (alpha_1 + alpha_2 + alpha_3); - double w_3 = alpha_3 / (alpha_1 + alpha_2 + alpha_3); - double left_state = w_1 * f1 + w_2 * f2 + w_3 * f3; + Real alpha_1 = (1.0 + (tau_4 / s1) * (tau_4 / s12)) / 3.0; + Real alpha_2 = (1.0 + (tau_4 / s2) * (tau_4 / s12)) / 6.0; + Real alpha_3 = (1.0 + tau_4 / s3) / 2.0; + Real w_1 = alpha_1 / (alpha_1 + alpha_2 + alpha_3); + Real w_2 = alpha_2 / (alpha_1 + alpha_2 + alpha_3); + Real w_3 = alpha_3 / (alpha_1 + alpha_2 + alpha_3); + Real left_state = w_1 * f1 + w_2 * f2 + w_3 * f3; return left_state; } //=================================================================================================// - double getRightStateInWeno(double v_1, double v_2, double v_3, double v_4) + Real getRightStateInWeno(Real v_1, Real v_2, Real v_3, Real v_4) { - double v1 = v_4; - double v2 = v_3; - double v3 = v_2; - double v4 = v_1; + Real v1 = v_4; + Real v2 = v_3; + Real v3 = v_2; + Real v4 = v_1; - double f1 = 0.5 * v2 + 0.5 * v3; - double f2 = (-0.5) * v1 + 1.5 * v2; - double f3 = v2 / 3.0 + 5.0 * v3 / 6.0 - v4 / 6.0; + Real f1 = 0.5 * v2 + 0.5 * v3; + Real f2 = (-0.5) * v1 + 1.5 * v2; + Real f3 = v2 / 3.0 + 5.0 * v3 / 6.0 - v4 / 6.0; - double epsilon = 1.0e-6; - double s1 = pow(v2 - v3, 2) + epsilon; - double s2 = pow(v2 - v1, 2) + epsilon; - double s3 = pow(3.0 * v2 - 4.0 * v3 + v4, 2) / 4.0 + 13.0 * pow(v2 - 2.0 * v3 + v4, 2) / 12.0 + epsilon; - double s12 = 13.0 * pow(v1 - 2.0 * v2 + v3, 2) / 12.0 + pow(v1 - v3, 2) / 4.0 + epsilon; - double tau_4 = (v1 * (547.0 * v1 - 2522.0 * v2 + 1922.0 * v3 - 494.0 * v4) + v2 * (3443.0 * v2 - 5966.0 * v3 + 1602.0 * v4) + v3 * (2843.0 * v3 - 1642.0 * v4) + 267.0 * v4 * v4) / 240.0; + Real epsilon = 1.0e-6; + Real s1 = pow(v2 - v3, 2) + epsilon; + Real s2 = pow(v2 - v1, 2) + epsilon; + Real s3 = pow(3.0 * v2 - 4.0 * v3 + v4, 2) / 4.0 + 13.0 * pow(v2 - 2.0 * v3 + v4, 2) / 12.0 + epsilon; + Real s12 = 13.0 * pow(v1 - 2.0 * v2 + v3, 2) / 12.0 + pow(v1 - v3, 2) / 4.0 + epsilon; + Real tau_4 = (v1 * (547.0 * v1 - 2522.0 * v2 + 1922.0 * v3 - 494.0 * v4) + v2 * (3443.0 * v2 - 5966.0 * v3 + 1602.0 * v4) + v3 * (2843.0 * v3 - 1642.0 * v4) + 267.0 * v4 * v4) / 240.0; - double alpha_1 = (1.0 + (tau_4 / s1) * (tau_4 / s12)) / 3.0; - double alpha_2 = (1.0 + (tau_4 / s2) * (tau_4 / s12)) / 6.0; - double alpha_3 = (1.0 + tau_4 / s3) / 2.0; - double w_1 = alpha_1 / (alpha_1 + alpha_2 + alpha_3); - double w_2 = alpha_2 / (alpha_1 + alpha_2 + alpha_3); - double w_3 = alpha_3 / (alpha_1 + alpha_2 + alpha_3); - double right_state = w_1 * f1 + w_2 * f2 + w_3 * f3; + Real alpha_1 = (1.0 + (tau_4 / s1) * (tau_4 / s12)) / 3.0; + Real alpha_2 = (1.0 + (tau_4 / s2) * (tau_4 / s12)) / 6.0; + Real alpha_3 = (1.0 + tau_4 / s3) / 2.0; + Real w_1 = alpha_1 / (alpha_1 + alpha_2 + alpha_3); + Real w_2 = alpha_2 / (alpha_1 + alpha_2 + alpha_3); + Real w_3 = alpha_3 / (alpha_1 + alpha_2 + alpha_3); + Real right_state = w_1 * f1 + w_2 * f2 + w_3 * f3; return right_state; } diff --git a/src/shared/common/scalar_functions.h b/src/shared/common/scalar_functions.h index 7a4a9c37e1..21dae3e7db 100644 --- a/src/shared/common/scalar_functions.h +++ b/src/shared/common/scalar_functions.h @@ -146,11 +146,11 @@ namespace SPH return (std::isnan(a) || !(std::isfinite(a))) ? true : false; } - inline double rand_norm(double u, double std) + inline Real rand_norm(Real u, Real std) { unsigned seed = (unsigned)std::chrono::system_clock::now().time_since_epoch().count(); std::default_random_engine generator(seed); - std::normal_distribution distribution(u, std); + std::normal_distribution distribution(u, std); return distribution(generator); } /** rotating axis once according to right hand rule. @@ -162,8 +162,8 @@ namespace SPH */ int ThirdAxis(int first_axis); - double getLeftStateInWeno(double v1, double v2, double v3, double v4); - double getRightStateInWeno(double v1, double v2, double v3, double v4); + Real getLeftStateInWeno(Real v1, Real v2, Real v3, Real v4); + Real getRightStateInWeno(Real v1, Real v2, Real v3, Real v4); /** linear heaviside function.*/ Real Heaviside(Real phi, Real half_width); diff --git a/src/shared/common/vector_functions.cpp b/src/shared/common/vector_functions.cpp index 246a0c7e9f..b3694cf549 100644 --- a/src/shared/common/vector_functions.cpp +++ b/src/shared/common/vector_functions.cpp @@ -3,45 +3,45 @@ namespace SPH { //=================================================================================================// - SimTK::Vec2 EigenToSimTK(const Vec2d &eigen_vector) + SimTKVec2 EigenToSimTK(const Vec2d &eigen_vector) { - return SimTK::Vec2(eigen_vector[0], eigen_vector[1]); + return SimTKVec2(eigen_vector[0], eigen_vector[1]); } //=================================================================================================// - SimTK::Vec3 EigenToSimTK(const Vec3d &eigen_vector) + SimTKVec3 EigenToSimTK(const Vec3d &eigen_vector) { - return SimTK::Vec3(eigen_vector[0], eigen_vector[1], eigen_vector[2]); + return SimTKVec3(eigen_vector[0], eigen_vector[1], eigen_vector[2]); } //=================================================================================================// - Vec2d SimTKToEigen(const SimTK::Vec2 &simTK_vector) + Vec2d SimTKToEigen(const SimTKVec2 &simTK_vector) { return Vec2d(simTK_vector[0], simTK_vector[1]); } //=================================================================================================// - Vec3d SimTKToEigen(const SimTK::Vec3 &simTK_vector) + Vec3d SimTKToEigen(const SimTKVec3 &simTK_vector) { return Vec3d(simTK_vector[0], simTK_vector[1], simTK_vector[2]); } - SimTK::Mat22 EigenToSimTK(const Mat2d &eigen_matrix) + SimTKMat22 EigenToSimTK(const Mat2d &eigen_matrix) { - return SimTK::Mat22(eigen_matrix(0, 0), eigen_matrix(0, 1), + return SimTKMat22(eigen_matrix(0, 0), eigen_matrix(0, 1), eigen_matrix(1, 0), eigen_matrix(1, 1)); } //=================================================================================================// - SimTK::Mat33 EigenToSimTK(const Mat3d &eigen_matrix) + SimTKMat33 EigenToSimTK(const Mat3d &eigen_matrix) { - return SimTK::Mat33(eigen_matrix(0, 0), eigen_matrix(0, 1), eigen_matrix(0, 2), + return SimTKMat33(eigen_matrix(0, 0), eigen_matrix(0, 1), eigen_matrix(0, 2), eigen_matrix(1, 0), eigen_matrix(1, 1), eigen_matrix(1, 2), eigen_matrix(2, 0), eigen_matrix(2, 1), eigen_matrix(2, 2)); } - Mat2d SimTKToEigen(const SimTK::Mat22 &simTK_matrix) + Mat2d SimTKToEigen(const SimTKMat22 &simTK_matrix) { return Mat2d{ {simTK_matrix(0, 0), simTK_matrix(0, 1)}, {simTK_matrix(1, 0), simTK_matrix(1, 1)}}; } //=================================================================================================// - Mat3d SimTKToEigen(const SimTK::Mat33 &simTK_matrix) + Mat3d SimTKToEigen(const SimTKMat33 &simTK_matrix) { return Mat3d{ {simTK_matrix(0, 0), simTK_matrix(0, 1), simTK_matrix(0, 2)}, @@ -167,7 +167,7 @@ namespace SPH { for (int j = 0; j < (i + 1); j++) { - double sum = 0; + Real sum = 0; for (int k = 0; k < j; k++) { sum += lower(i, k) * lower(j, k); @@ -195,7 +195,7 @@ namespace SPH { for (int j = 0; j < (i + 1); j++) { - double sum = 0; + Real sum = 0; for (int k = 0; k < j; k++) { sum += lower(i, k) * lower(j, k); @@ -261,7 +261,7 @@ namespace SPH return diag; } //=================================================================================================// - Real CalculateDoubleDotProduct(Mat2d Matrix1, Mat2d Matrix2) + Real CalculateBiDotProduct(Mat2d Matrix1, Mat2d Matrix2) { Real product = 0; for (int i = 0; i < 2; i++) @@ -273,7 +273,7 @@ namespace SPH } return product; } - Real CalculateDoubleDotProduct(Mat3d Matrix1, Mat3d Matrix2) + Real CalculateBiDotProduct(Mat3d Matrix1, Mat3d Matrix2) { Real product = 0; for (int i = 0; i < 3; i++) @@ -369,7 +369,7 @@ namespace SPH //=================================================================================================// Vec2d getPrincipalValuesFromMatrix(const Mat2d &A) { - Eigen::EigenSolver ces(A, /* computeEigenvectors = */ false); + Eigen::EigenSolver ces(A, /* computeEigenvectors = */ false); auto eigen_values = ces.eigenvalues(); std::vector sorted_values = { @@ -384,7 +384,7 @@ namespace SPH //=================================================================================================// Vec3d getPrincipalValuesFromMatrix(const Mat3d &A) { - Eigen::EigenSolver ces(A, /* computeEigenvectors = */ false); + Eigen::EigenSolver ces(A, /* computeEigenvectors = */ false); auto eigen_values = ces.eigenvalues(); std::vector sorted_values = { diff --git a/src/shared/common/vector_functions.h b/src/shared/common/vector_functions.h index 6d32c89468..65fd054d01 100644 --- a/src/shared/common/vector_functions.h +++ b/src/shared/common/vector_functions.h @@ -33,15 +33,15 @@ namespace SPH { - SimTK::Vec2 EigenToSimTK(const Vec2d &eigen_vector); - SimTK::Vec3 EigenToSimTK(const Vec3d &eigen_vector); - Vec2d SimTKToEigen(const SimTK::Vec2 &simTK_vector); - Vec3d SimTKToEigen(const SimTK::Vec3 &simTK_vector); + SimTKVec2 EigenToSimTK(const Vec2d &eigen_vector); + SimTKVec3 EigenToSimTK(const Vec3d &eigen_vector); + Vec2d SimTKToEigen(const SimTKVec2 &simTK_vector); + Vec3d SimTKToEigen(const SimTKVec3 &simTK_vector); - SimTK::Mat22 EigenToSimTK(const Mat2d &eigen_matrix); - SimTK::Mat33 EigenToSimTK(const Mat3d &eigen_matrix); - Mat2d SimTKToEigen(const SimTK::Mat22 &simTK_matrix); - Mat3d SimTKToEigen(const SimTK::Mat33 &simTK_matrix); + SimTKMat22 EigenToSimTK(const Mat2d &eigen_matrix); + SimTKMat33 EigenToSimTK(const Mat3d &eigen_matrix); + Mat2d SimTKToEigen(const SimTKMat22 &simTK_matrix); + Mat3d SimTKToEigen(const SimTKMat33 &simTK_matrix); Vec2d FirstAxisVector(const Vec2d &zero_vector); Vec3d FirstAxisVector(const Vec3d &zero_vector); @@ -66,9 +66,9 @@ namespace SPH Mat2d getDiagonal(const Mat2d &A); Mat3d getDiagonal(const Mat3d &A); - /** double dot product between two matrices, resulting in a scalar value (sum of products of element-wise) */ - Real CalculateDoubleDotProduct(Mat2d Matrix1, Mat2d Matrix2); // calculate double dot - Real CalculateDoubleDotProduct(Mat3d Matrix1, Mat3d Matrix2); // calculate double dot + /** Real dot product between two matrices, resulting in a scalar value (sum of products of element-wise) */ + Real CalculateBiDotProduct(Mat2d Matrix1, Mat2d Matrix2); // calculate Real dot + Real CalculateBiDotProduct(Mat3d Matrix1, Mat3d Matrix2); // calculate Real dot /** get transformation matrix. */ Mat2d getTransformationMatrix(const Vec2d &direction_of_y); diff --git a/src/shared/geometries/level_set.h b/src/shared/geometries/level_set.h index dac7d939ae..3c6b9e6ead 100644 --- a/src/shared/geometries/level_set.h +++ b/src/shared/geometries/level_set.h @@ -138,7 +138,7 @@ namespace SPH /** * @class RefinedLevelSet - * @brief level set which has double resolution of a coarse level set. + * @brief level set which has Real resolution of a coarse level set. */ class RefinedLevelSet : public RefinedMesh { diff --git a/src/shared/materials/elastic_solid.cpp b/src/shared/materials/elastic_solid.cpp index 0cc64ba337..abcc1b849b 100644 --- a/src/shared/materials/elastic_solid.cpp +++ b/src/shared/materials/elastic_solid.cpp @@ -163,8 +163,8 @@ namespace SPH for (int j = 0; j < Dimensions; j++) { // inner sum (b{1-3}) - Summa2 += Lambda_(i, j) * (CalculateDoubleDotProduct(A_[i], strain) * A_[j] + - CalculateDoubleDotProduct(A_[j], strain) * A_[i]); + Summa2 += Lambda_(i, j) * (CalculateBiDotProduct(A_[i], strain) * A_[j] + + CalculateBiDotProduct(A_[j], strain) * A_[i]); } stress_PK2 += Mu_[i] * (((A_[i] * strain) + (strain * A_[i])) + 1 / 2 * (Summa2)); } diff --git a/src/shared/meshes/base_mesh.h b/src/shared/meshes/base_mesh.h index a09624d818..a6b8be028f 100644 --- a/src/shared/meshes/base_mesh.h +++ b/src/shared/meshes/base_mesh.h @@ -147,7 +147,7 @@ namespace SPH /** * @class RefinedMesh - * @brief Abstract base class derived from the coarse mesh but has double resolution. + * @brief Abstract base class derived from the coarse mesh but has Real resolution. * Currently, the design is simple but can be extending for more inter-mesh operations. */ template @@ -166,7 +166,7 @@ namespace SPH /** * @class MultilevelMesh - * @brief Multi-level Meshes with successively double the resolution + * @brief Multi-level Meshes with successively Real the resolution */ template class MultilevelMesh : public MeshFieldType diff --git a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp index 57fd4c45bb..c102cf3ade 100644 --- a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp +++ b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp @@ -436,7 +436,7 @@ namespace SPH template bool DampingWithRandomChoice::RandomChoice() { - return ((double)rand() / (RAND_MAX)) < random_ratio_ ? true : false; + return ((Real)rand() / (RAND_MAX)) < random_ratio_ ? true : false; } //=================================================================================================// template diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp index 47dd8c0f62..94ec91997a 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp @@ -22,7 +22,7 @@ namespace SPH Vecd &pos_n_i = pos_[index_i]; for (int k = 0; k < pos_n_i.size(); ++k) { - pos_n_i[k] += dt * (((double)rand() / (RAND_MAX)) - 0.5) * 2.0 * randomize_scale_; + pos_n_i[k] += dt * (((Real)rand() / (RAND_MAX)) - 0.5) * 2.0 * randomize_scale_; } } //=================================================================================================// diff --git a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h index 60a5bce87b..dfb3d477f4 100644 --- a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h @@ -259,19 +259,19 @@ namespace SPH void update(size_t index_i, Real dt = 0.0) { /** Change to SimTK::Vector. */ - SimTK::Vec3 rr, pos, vel, acc; + SimTKVec3 rr, pos, vel, acc; rr = EigenToSimTK(upgradeToVec3d(this->pos0_[index_i])) - initial_mobod_origin_location_; mobod_.findStationLocationVelocityAndAccelerationInGround(*simbody_state_, rr, pos, vel, acc); /** this is how we calculate the particle position in after transform of MBbody. * const SimTK::Rotation& R_GB = mobod_.getBodyRotation(simbody_state); - * const SimTK::Vec3& p_GB = mobod_.getBodyOriginLocation(simbody_state); - * const SimTK::Vec3 r = R_GB * rr; // re-express station vector p_BS in G (15 flops) + * const SimTKVec3& p_GB = mobod_.getBodyOriginLocation(simbody_state); + * const SimTKVec3 r = R_GB * rr; // re-express station vector p_BS in G (15 flops) * base_particle_data_i.pos_ = (p_GB + r); */ degradeToVecd(SimTKToEigen(pos), this->pos_[index_i]); degradeToVecd(SimTKToEigen(vel), this->vel_[index_i]); - SimTK::Vec3 n = (mobod_.getBodyRotation(*simbody_state_) * EigenToSimTK(upgradeToVec3d(this->n0_[index_i]))); + SimTKVec3 n = (mobod_.getBodyRotation(*simbody_state_) * EigenToSimTK(upgradeToVec3d(this->n0_[index_i]))); degradeToVecd(SimTKToEigen(n), this->n_[index_i]); }; @@ -281,7 +281,7 @@ namespace SPH SimTK::Force::DiscreteForces &force_on_bodies_; SimTK::RungeKuttaMersonIntegrator &integ_; const SimTK::State *simbody_state_; - SimTK::Vec3 initial_mobod_origin_location_; + SimTKVec3 initial_mobod_origin_location_; }; using ConstraintBodyBySimBody = ConstraintBySimBody; using ConstraintBodyPartBySimBody = ConstraintBySimBody; @@ -304,7 +304,7 @@ namespace SPH SimTK::Force::DiscreteForces &force_on_bodies_; SimTK::RungeKuttaMersonIntegrator &integ_; const SimTK::State *simbody_state_; - SimTK::Vec3 current_mobod_origin_location_; + SimTKVec3 current_mobod_origin_location_; public: TotalForceForSimBody(DynamicsIdentifier &identifier, @@ -313,7 +313,7 @@ namespace SPH SimTK::Force::DiscreteForces &force_on_bodies, SimTK::RungeKuttaMersonIntegrator &integ) : BaseLocalDynamicsReduce, DynamicsIdentifier>( - identifier, SimTK::SpatialVec(SimTK::Vec3(0), SimTK::Vec3(0))), + identifier, SimTK::SpatialVec(SimTKVec3(0), SimTKVec3(0))), SolidDataSimple(identifier.getSPHBody()), mass_(particles_->mass_), acc_(particles_->acc_), acc_prior_(particles_->acc_prior_), pos_(particles_->pos_), @@ -334,9 +334,9 @@ namespace SPH SimTK::SpatialVec reduce(size_t index_i, Real dt = 0.0) { Vecd force = (acc_[index_i] + acc_prior_[index_i]) * mass_[index_i]; - SimTK::Vec3 force_from_particle = EigenToSimTK(upgradeToVec3d(force)); - SimTK::Vec3 displacement = EigenToSimTK(upgradeToVec3d(pos_[index_i])) - current_mobod_origin_location_; - SimTK::Vec3 torque_from_particle = SimTK::cross(displacement, force_from_particle); + SimTKVec3 force_from_particle = EigenToSimTK(upgradeToVec3d(force)); + SimTKVec3 displacement = EigenToSimTK(upgradeToVec3d(pos_[index_i])) - current_mobod_origin_location_; + SimTKVec3 torque_from_particle = SimTK::cross(displacement, force_from_particle); return SimTK::SpatialVec(torque_from_particle, force_from_particle); }; diff --git a/src/shared/particle_generator/particle_generator_lattice.cpp b/src/shared/particle_generator/particle_generator_lattice.cpp index e2d2ec4126..0683436458 100644 --- a/src/shared/particle_generator/particle_generator_lattice.cpp +++ b/src/shared/particle_generator/particle_generator_lattice.cpp @@ -40,7 +40,7 @@ namespace SPH { Real local_particle_spacing = particle_adaptation_->getLocalSpacing(target_shape_, position); Real local_particle_volume_ratio = pow(lattice_spacing_ / local_particle_spacing, Dimensions); - if ((double)rand() / (RAND_MAX) < local_particle_volume_ratio) + if ((Real)rand() / (RAND_MAX) < local_particle_volume_ratio) { ParticleGeneratorLattice::initializePositionAndVolumetricMeasure(position, volume / local_particle_volume_ratio); initializeSmoothingLengthRatio(local_particle_spacing); diff --git a/src/shared/regression_test/dynamic_time_warping_method.h b/src/shared/regression_test/dynamic_time_warping_method.h index dad60bb1b8..717bee18ca 100644 --- a/src/shared/regression_test/dynamic_time_warping_method.h +++ b/src/shared/regression_test/dynamic_time_warping_method.h @@ -62,7 +62,7 @@ namespace SPH }; /** the local constrained method used for calculating the dtw distance between two lines. */ - StdVec calculateDTWDistance(DoubleVec dataset_a_, DoubleVec dataset_b_); + StdVec calculateDTWDistance(BiVector dataset_a_, BiVector dataset_b_); public: template diff --git a/src/shared/regression_test/dynamic_time_warping_method.hpp b/src/shared/regression_test/dynamic_time_warping_method.hpp index a5de92f4fd..bb7cbdce12 100644 --- a/src/shared/regression_test/dynamic_time_warping_method.hpp +++ b/src/shared/regression_test/dynamic_time_warping_method.hpp @@ -37,7 +37,7 @@ namespace SPH //=================================================================================================// template StdVec RegressionTestDynamicTimeWarping::calculateDTWDistance - (DoubleVec dataset_a_, DoubleVec dataset_b_) + (BiVector dataset_a_, BiVector dataset_b_) { /* define the container to hold the dtw distance.*/ StdVec dtw_distance; @@ -47,7 +47,7 @@ namespace SPH int a_length = dataset_a_[observation_index].size(); int b_length = dataset_b_[observation_index].size(); /** create a 2D vector with [a_length, b_length] to contain the local DTW distance value. */ - DoubleVec local_dtw_distance(a_length, StdVec(b_length, 0)); + BiVector local_dtw_distance(a_length, StdVec(b_length, 0)); local_dtw_distance[0][0] = calculatePNorm(dataset_a_[observation_index][0], dataset_b_[observation_index][0]); for (int index_i = 1; index_i < a_length; ++index_i) local_dtw_distance[index_i][0] = local_dtw_distance[index_i - 1][0] + calculatePNorm(dataset_a_[observation_index][index_i], dataset_b_[observation_index][0]); diff --git a/src/shared/regression_test/ensemble_averaged_method.h b/src/shared/regression_test/ensemble_averaged_method.h index ea4b0161cb..d388ec7802 100644 --- a/src/shared/regression_test/ensemble_averaged_method.h +++ b/src/shared/regression_test/ensemble_averaged_method.h @@ -43,23 +43,23 @@ namespace SPH using VariableType = decltype(ObserveMethodType::type_indicator_); protected: - DoubleVec meanvalue_, meanvalue_new_; /* the container of (new) mean value. [different from time-averaged]*/ - DoubleVec variance_, variance_new_; /* the container of (new) variance. [different from time-averaged]*/ + BiVector meanvalue_, meanvalue_new_; /* the container of (new) mean value. [different from time-averaged]*/ + BiVector variance_, variance_new_; /* the container of (new) variance. [different from time-averaged]*/ /** the method used for calculating the new variance. */ - void calculateNewVariance(TripleVec &result, DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new); - void calculateNewVariance(TripleVec &result, DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new); - void calculateNewVariance(TripleVec &result, DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new); + void calculateNewVariance(TriVector &result, BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new); + void calculateNewVariance(TriVector &result, BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new); + void calculateNewVariance(TriVector &result, BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new); /** the method used for comparing the meanvalue and variance. */ - int compareParameter(std::string par_name, DoubleVec ¶meter, DoubleVec ¶meter_new, Real &threshold); - int compareParameter(std::string par_name, DoubleVec ¶meter, DoubleVec ¶meter_new, Vecd &threshold); - int compareParameter(std::string par_name, DoubleVec ¶meter, DoubleVec ¶meter_new, Matd &threshold); + int compareParameter(std::string par_name, BiVector ¶meter, BiVector ¶meter_new, Real &threshold); + int compareParameter(std::string par_name, BiVector ¶meter, BiVector ¶meter_new, Vecd &threshold); + int compareParameter(std::string par_name, BiVector ¶meter, BiVector ¶meter_new, Matd &threshold); /** the method used for testing the new result with meanvalue and variance. */ - int testNewResult(int diff, DoubleVec ¤t_result, DoubleVec &meanvalue, DoubleVec &variance); - int testNewResult(int diff, DoubleVec ¤t_result, DoubleVec &meanvalue, DoubleVec &variance); - int testNewResult(int diff, DoubleVec ¤t_result, DoubleVec &meanvalue, DoubleVec &variance); + int testNewResult(int diff, BiVector ¤t_result, BiVector &meanvalue, BiVector &variance); + int testNewResult(int diff, BiVector ¤t_result, BiVector &meanvalue, BiVector &variance); + int testNewResult(int diff, BiVector ¤t_result, BiVector &meanvalue, BiVector &variance); public: template diff --git a/src/shared/regression_test/ensemble_averaged_method.hpp b/src/shared/regression_test/ensemble_averaged_method.hpp index ede9593553..d373066d1c 100644 --- a/src/shared/regression_test/ensemble_averaged_method.hpp +++ b/src/shared/regression_test/ensemble_averaged_method.hpp @@ -35,8 +35,8 @@ namespace SPH { //=================================================================================================// template - void RegressionTestEnsembleAveraged::calculateNewVariance(TripleVec &result, - DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new) + void RegressionTestEnsembleAveraged::calculateNewVariance(TriVector &result, + BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new) { for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -50,8 +50,8 @@ namespace SPH }; //=================================================================================================// template - void RegressionTestEnsembleAveraged::calculateNewVariance(TripleVec &result, - DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new) + void RegressionTestEnsembleAveraged::calculateNewVariance(TriVector &result, + BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new) { for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -66,8 +66,8 @@ namespace SPH }; //=================================================================================================// template - void RegressionTestEnsembleAveraged::calculateNewVariance(TripleVec &result, - DoubleVec &meanvalue_new, DoubleVec &variance, DoubleVec &variance_new) + void RegressionTestEnsembleAveraged::calculateNewVariance(TriVector &result, + BiVector &meanvalue_new, BiVector &variance, BiVector &variance_new) { for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -84,7 +84,7 @@ namespace SPH //=================================================================================================// template int RegressionTestEnsembleAveraged::compareParameter(std::string par_name, - DoubleVec ¶meter, DoubleVec ¶meter_new, Real &threshold) + BiVector ¶meter, BiVector ¶meter_new, Real &threshold) { int count = 0; for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) @@ -104,7 +104,7 @@ namespace SPH ////=================================================================================================// template int RegressionTestEnsembleAveraged::compareParameter(std::string par_name, - DoubleVec ¶meter, DoubleVec ¶meter_new, Vecd &threshold) + BiVector ¶meter, BiVector ¶meter_new, Vecd &threshold) { int count = 0; for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) @@ -125,7 +125,7 @@ namespace SPH ////=================================================================================================// template int RegressionTestEnsembleAveraged::compareParameter(std::string par_name, - DoubleVec ¶meter, DoubleVec ¶meter_new, Matd &threshold) + BiVector ¶meter, BiVector ¶meter_new, Matd &threshold) { int count = 0; for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) @@ -146,8 +146,8 @@ namespace SPH } //=================================================================================================// template - int RegressionTestEnsembleAveraged::testNewResult(int diff, DoubleVec ¤t_result, - DoubleVec &meanvalue, DoubleVec &variance) + int RegressionTestEnsembleAveraged::testNewResult(int diff, BiVector ¤t_result, + BiVector &meanvalue, BiVector &variance) { int count = 0; for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) @@ -168,8 +168,8 @@ namespace SPH } //=================================================================================================// template - int RegressionTestEnsembleAveraged::testNewResult(int diff, DoubleVec ¤t_result, - DoubleVec &meanvalue, DoubleVec &variance) + int RegressionTestEnsembleAveraged::testNewResult(int diff, BiVector ¤t_result, + BiVector &meanvalue, BiVector &variance) { int count = 0; for (int snapshot_index = 0; snapshot_index != SMIN(this->snapshot_, this->number_of_snapshot_old_); ++snapshot_index) @@ -197,8 +197,8 @@ namespace SPH } //=================================================================================================// template - int RegressionTestEnsembleAveraged::testNewResult(int diff, DoubleVec ¤t_result, - DoubleVec &meanvalue, DoubleVec &variance) + int RegressionTestEnsembleAveraged::testNewResult(int diff, BiVector ¤t_result, + BiVector &meanvalue, BiVector &variance) { int count = 0; std::cout << "The current length difference is " << diff << "." << std::endl; @@ -258,7 +258,7 @@ namespace SPH SimTK::Xml::Element mean_element_ = this->mean_variance_xml_engine_in_.getChildElement("Mean_Element"); this->number_of_snapshot_old_ = std::distance(mean_element_.element_begin(), mean_element_.element_end()); - DoubleVec temp(SMAX(this->snapshot_, this->number_of_snapshot_old_), StdVec(this->observation_)); + BiVector temp(SMAX(this->snapshot_, this->number_of_snapshot_old_), StdVec(this->observation_)); meanvalue_ = temp; variance_ = temp; @@ -278,7 +278,7 @@ namespace SPH else if (this->number_of_run_ == 1) { this->number_of_snapshot_old_ = this->snapshot_; - DoubleVec temp(this->snapshot_, StdVec(this->observation_)); + BiVector temp(this->snapshot_, StdVec(this->observation_)); this->result_.push_back(this->current_result_); meanvalue_ = temp; variance_ = temp; diff --git a/src/shared/regression_test/regression_test_base.h b/src/shared/regression_test/regression_test_base.h index bed488fbeb..3c5bfa52f3 100644 --- a/src/shared/regression_test/regression_test_base.h +++ b/src/shared/regression_test/regression_test_base.h @@ -64,13 +64,13 @@ namespace SPH /*< the XmlEngine can operate the node name and elements in xml memory. */ StdVec element_tag_; /*< the container of the tag of current result. */ - DoubleVec current_result_; /*< the container of current run result stored as snapshot * observation. */ - DoubleVec current_result_trans_; /*< the container of current run result with snapshot & observations transposed, + BiVector current_result_; /*< the container of current run result stored as snapshot * observation. */ + BiVector current_result_trans_; /*< the container of current run result with snapshot & observations transposed, because this data structure is required in TA and DTW method. */ - DoubleVec result_in_; /*< the temporary container of each result when reading from .xml memory, and + BiVector result_in_; /*< the temporary container of each result when reading from .xml memory, and it may be snapshot * observations (reading all result for averaged methods), or observations * snapshot (reading specified result for TA and DTW method.) */ - TripleVec result_; /*< the container of results in all runs (run * snapshot * observation) */ + TriVector result_; /*< the container of results in all runs (run * snapshot * observation) */ int snapshot_, observation_; /*< the size of each layer of current result vector. */ int difference_; /*< the length difference of snapshot between old and new result, diff --git a/src/shared/regression_test/regression_test_base.hpp b/src/shared/regression_test/regression_test_base.hpp index ca68d805d2..698537434b 100644 --- a/src/shared/regression_test/regression_test_base.hpp +++ b/src/shared/regression_test/regression_test_base.hpp @@ -66,7 +66,7 @@ namespace SPH size_t number_of_particle_ = this->base_particles_.total_real_particles_; size_t number_of_snapshot_ = std::distance(observe_xml_engine_.root_element_.element_begin(), observe_xml_engine_.root_element_.element_end()); - DoubleVec current_result_temp_(number_of_snapshot_, StdVec(number_of_particle_)); + BiVector current_result_temp_(number_of_snapshot_, StdVec(number_of_particle_)); StdVec element_tag_temp_(number_of_snapshot_); current_result_ = current_result_temp_; element_tag_ = element_tag_temp_; @@ -86,7 +86,7 @@ namespace SPH size_t number_of_particle_ = 1; size_t number_of_snapshot_ = std::distance(observe_xml_engine_.root_element_.element_begin(), observe_xml_engine_.root_element_.element_end()); - DoubleVec current_result_temp_(number_of_snapshot_, StdVec(number_of_particle_)); + BiVector current_result_temp_(number_of_snapshot_, StdVec(number_of_particle_)); StdVec element_tag_temp_(number_of_snapshot_); current_result_ = current_result_temp_; element_tag_ = element_tag_temp_; @@ -103,7 +103,7 @@ namespace SPH { int number_of_snapshot = this->current_result_.size(); int number_of_observation = this->current_result_[0].size(); - DoubleVec temp(number_of_observation, StdVec(number_of_snapshot)); + BiVector temp(number_of_observation, StdVec(number_of_snapshot)); current_result_trans_ = temp; for (int snapshot_index = 0; snapshot_index != number_of_snapshot; ++snapshot_index) for (int observation_index = 0; observation_index != number_of_observation; ++observation_index) @@ -116,14 +116,14 @@ namespace SPH if (number_of_run_ > 1) /*only read the result from the 2nd run, because the 1st run doesn't have previous results. */ { /*Here result_in is a temporary container that reloads each previous result.*/ - DoubleVec result_in_(SMAX(snapshot_, number_of_snapshot_old_), StdVec(observation_)); + BiVector result_in_(SMAX(snapshot_, number_of_snapshot_old_), StdVec(observation_)); for (int run_index_ = 0; run_index_ != number_of_run_ - 1; ++run_index_) { std::string node_name_ = "Round_" + std::to_string(run_index_); SimTK::Xml::Element father_element_ = result_xml_engine_in_.getChildElement(node_name_); for (int observation_index_ = 0; observation_index_ != observation_; ++observation_index_) xmlmemory_io_.readDataFromXmlMemory(result_xml_engine_in_, father_element_, observation_index_, result_in_, this->quantity_name_); - DoubleVec result_temp_ = result_in_; + BiVector result_temp_ = result_in_; for (int delete_ = 0; delete_ != difference_; ++delete_) result_temp_.pop_back(); /* trim the new reading result to unify the length of all results. (number of snapshots) */ result_.push_back(result_temp_); @@ -172,7 +172,7 @@ namespace SPH SimTK::Xml::element_iterator ele_ite = snapshot_element_.element_begin(); result_xml_engine_in_.getRequiredAttributeValue(ele_ite, "number_of_snapshot_for_local_result_", snapshot_); - DoubleVec result_temp_(observation_, StdVec(snapshot_)); + BiVector result_temp_(observation_, StdVec(snapshot_)); result_in_ = result_temp_; SimTK::Xml::Element result_element_ = result_xml_engine_in_.getChildElement("Result_Element"); for (int snapshot_index = 0; snapshot_index != snapshot_; ++snapshot_index) diff --git a/src/shared/regression_test/time_averaged_method.h b/src/shared/regression_test/time_averaged_method.h index 675c8ee16c..a573be091d 100644 --- a/src/shared/regression_test/time_averaged_method.h +++ b/src/shared/regression_test/time_averaged_method.h @@ -55,19 +55,19 @@ namespace SPH StdVec local_meanvalue_; /* the container of meanvalue for current result. */ /** the method used for filtering local extreme values. */ - void filterLocalResult(DoubleVec ¤t_result); - void filterLocalResult(DoubleVec ¤t_result); - void filterLocalResult(DoubleVec ¤t_result); + void filterLocalResult(BiVector ¤t_result); + void filterLocalResult(BiVector ¤t_result); + void filterLocalResult(BiVector ¤t_result); /** the method used for searching steady starting point. */ - void searchSteadyStart(DoubleVec ¤t_result); - void searchSteadyStart(DoubleVec ¤t_result); - void searchSteadyStart(DoubleVec ¤t_result); + void searchSteadyStart(BiVector ¤t_result); + void searchSteadyStart(BiVector ¤t_result); + void searchSteadyStart(BiVector ¤t_result); /** the method used for calculating the new variance. */ - void calculateNewVariance(DoubleVec ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); - void calculateNewVariance(DoubleVec ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); - void calculateNewVariance(DoubleVec ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); + void calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); + void calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); + void calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new); /** the method used for comparing the meanvalue and variance. */ int compareParameter(std::string par_name, StdVec ¶meter, StdVec ¶meter_new, Real &threshold); @@ -75,9 +75,9 @@ namespace SPH int compareParameter(std::string par_name, StdVec ¶meter, StdVec ¶meter_new, Matd &threshold); /** the method used for testing the new result with meanvalue and variance. */ - int testNewResult(DoubleVec ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); - int testNewResult(DoubleVec ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); - int testNewResult(DoubleVec ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); + int testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); + int testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); + int testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance); public: template diff --git a/src/shared/regression_test/time_averaged_method.hpp b/src/shared/regression_test/time_averaged_method.hpp index ff39fb6555..6a408317dd 100644 --- a/src/shared/regression_test/time_averaged_method.hpp +++ b/src/shared/regression_test/time_averaged_method.hpp @@ -35,7 +35,7 @@ namespace SPH { //=================================================================================================// template - void RegressionTestTimeAveraged::filterLocalResult(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::filterLocalResult(BiVector ¤t_result) { int scale = round(this->snapshot_ / 200); std::cout << "The filter scale is " << scale * 2 << "." << std::endl; @@ -69,7 +69,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::filterLocalResult(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::filterLocalResult(BiVector ¤t_result) { int scale = round(this->snapshot_ / 200); std::cout << "The filter scale is " << scale * 2 << "." << std::endl; @@ -106,7 +106,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::filterLocalResult(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::filterLocalResult(BiVector ¤t_result) { int scale = round(this->snapshot_ / 200); std::cout << "The filter scale is " << scale * 2 << "." << std::endl; @@ -146,7 +146,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::searchSteadyStart(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::searchSteadyStart(BiVector ¤t_result) { /* the search is only for one value. */ int scale = round(this->snapshot_ / 20); @@ -170,7 +170,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::searchSteadyStart(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::searchSteadyStart(BiVector ¤t_result) { /* the search is for each value within parameters. */ int scale = round(this->snapshot_ / 20); @@ -194,7 +194,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::searchSteadyStart(DoubleVec ¤t_result) + void RegressionTestTimeAveraged::searchSteadyStart(BiVector ¤t_result) { int scale = round(this->snapshot_ / 20); for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -217,7 +217,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::calculateNewVariance(DoubleVec ¤t_result, + void RegressionTestTimeAveraged::calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new) { for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -230,7 +230,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::calculateNewVariance(DoubleVec ¤t_result, + void RegressionTestTimeAveraged::calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new) { for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -244,7 +244,7 @@ namespace SPH } //=================================================================================================// template - void RegressionTestTimeAveraged::calculateNewVariance(DoubleVec ¤t_result, + void RegressionTestTimeAveraged::calculateNewVariance(BiVector ¤t_result, StdVec &local_meanvalue, StdVec &variance, StdVec &variance_new) { for (int observation_index = 0; observation_index != this->observation_; ++observation_index) @@ -335,7 +335,7 @@ namespace SPH } //=================================================================================================// template - int RegressionTestTimeAveraged::testNewResult(DoubleVec ¤t_result, + int RegressionTestTimeAveraged::testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance) { int count = 0; @@ -365,7 +365,7 @@ namespace SPH } //=================================================================================================// template - int RegressionTestTimeAveraged::testNewResult(DoubleVec ¤t_result, + int RegressionTestTimeAveraged::testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance) { int count = 0; @@ -399,7 +399,7 @@ namespace SPH } //=================================================================================================// template - int RegressionTestTimeAveraged::testNewResult(DoubleVec ¤t_result, + int RegressionTestTimeAveraged::testNewResult(BiVector ¤t_result, StdVec &meanvalue, StdVec &local_meanvalue, StdVec &variance) { int count = 0; diff --git a/src/shared/simbody_sphinxsys/state_engine.cpp b/src/shared/simbody_sphinxsys/state_engine.cpp index 4cfd3723f9..2c401b80c9 100644 --- a/src/shared/simbody_sphinxsys/state_engine.cpp +++ b/src/shared/simbody_sphinxsys/state_engine.cpp @@ -186,17 +186,17 @@ namespace SPH { throw (msg.str()); } //=============================================================================================// - double StateEngine::AddedStateVariable:: + Real StateEngine::AddedStateVariable:: getDerivative() { - //return getCacheVariableValue(state, getName()+"_deriv"); + //return getCacheVariableValue(state, getName()+"_deriv"); return 0.0; } //=============================================================================================// void StateEngine::AddedStateVariable:: setDerivative(Real deriv) { - //return setCacheVariableValue(state, getName()+"_deriv", deriv); + //return setCacheVariableValue(state, getName()+"_deriv", deriv); } //=============================================================================================// void StateEngine::reporter(SimTK::State& state_) @@ -260,7 +260,7 @@ namespace SPH { std::string ele_name = "UIndx_" + std::to_string(i); simbody_xml_engine_.setAttributeToElement(ele_ite, ele_name, mobod_u); } - SimTK::Vec3 transform_ = mobod.getBodyTransform(state_).p(); + SimTKVec3 transform_ = mobod.getBodyTransform(state_).p(); Vecd transform = Vecd::Zero(); if constexpr (Dimensions == 2) @@ -321,13 +321,13 @@ namespace SPH { Vecd transform = Vecd::Zero(); simbody_xml_engine_.getRequiredAttributeValue(ele_ite_, "Transform", transform); - SimTK::Vec3 transform_; + SimTKVec3 transform_; if(Dimensions == 2) { - transform_ = SimTK::Vec3(transform[0],transform[1], 0); + transform_ = SimTKVec3(transform[0],transform[1], 0); }else { - transform_ = SimTK::Vec3(transform[0],transform[1], transform[2]); + transform_ = SimTKVec3(transform[0],transform[1], transform[2]); } mobod.setQToFitTransform(state_, SimTK::Transform(transform_)); diff --git a/src/shared/simbody_sphinxsys/xml_engine.h b/src/shared/simbody_sphinxsys/xml_engine.h index 4bf4c908e9..7e4ee19834 100644 --- a/src/shared/simbody_sphinxsys/xml_engine.h +++ b/src/shared/simbody_sphinxsys/xml_engine.h @@ -144,7 +144,7 @@ namespace SPH virtual ~XmlMemoryIO(){}; template - void writeDataToXmlMemory(XmlEngine &xml_engine, SimTK::Xml::Element &element, const DoubleVec &quantity, + void writeDataToXmlMemory(XmlEngine &xml_engine, SimTK::Xml::Element &element, const BiVector &quantity, int snapshot_, int observation_, const std::string &quantity_name, StdVec &element_tag) { for (int snapshot_index = 0; snapshot_index != snapshot_; ++snapshot_index) @@ -171,7 +171,7 @@ namespace SPH template void readDataFromXmlMemory(XmlEngine &xml_engine, SimTK::Xml::Element &element, - int observation_index, DoubleVec &result_container, const std::string &quantity_name) + int observation_index, BiVector &result_container, const std::string &quantity_name) { int snapshot_index = 0; SimTK::Xml::element_iterator ele_ite = element.element_begin(); diff --git a/tests/2d_examples/test_2d_owsc/owsc.cpp b/tests/2d_examples/test_2d_owsc/owsc.cpp index 88df051d51..4ef3be0990 100644 --- a/tests/2d_examples/test_2d_owsc/owsc.cpp +++ b/tests/2d_examples/test_2d_owsc/owsc.cpp @@ -95,13 +95,13 @@ int main(int ac, char *av[]) * @details Create a Pin mobilizer between an existing parent (inboard) body P and * a new child (outboard) body B created by copying the given bodyInfo into * a privately-owned Body within the constructed MobilizedBody object. - * @param[in] inboard(SimTK::Vec3) Defines the location of the joint point relative to the parent body. - * @param[in] outboard(SimTK::Vec3) Defines the body's origin location to the joint point. - * @note The body's origin location can be the mass center, the the center of mass should be SimTK::Vec3(0) + * @param[in] inboard(SimTKVec3) Defines the location of the joint point relative to the parent body. + * @param[in] outboard(SimTKVec3) Defines the body's origin location to the joint point. + * @note The body's origin location can be the mass center, the the center of mass should be SimTKVec3(0) * in SimTK::MassProperties(mass, com, inertia) */ - SimTK::MobilizedBody::Pin pin_spot(matter.Ground(), SimTK::Transform(SimTK::Vec3(7.92, 0.315, 0.0)), - pin_spot_info, SimTK::Transform(SimTK::Vec3(0.0, 0.0, 0.0))); + SimTK::MobilizedBody::Pin pin_spot(matter.Ground(), SimTK::Transform(SimTKVec3(7.92, 0.315, 0.0)), + pin_spot_info, SimTK::Transform(SimTKVec3(0.0, 0.0, 0.0))); /** set the default angle of the pin. */ pin_spot.setDefaultAngle(0); /** @@ -128,7 +128,7 @@ int main(int ac, char *av[]) * hb=pb*(-d) - hz. Note that this is a signed quantity so the potential energy is * also signed. 0.475 */ - SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(0.0, -gravity_g, 0.0), 0.0); + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTKVec3(0.0, -gravity_g, 0.0), 0.0); /** discrete forces acting on the bodies. */ SimTK::Force::DiscreteForces force_on_bodies(forces, matter); /** diff --git a/tests/2d_examples/test_2d_owsc/owsc.h b/tests/2d_examples/test_2d_owsc/owsc.h index da20a3c693..6aee2d9934 100644 --- a/tests/2d_examples/test_2d_owsc/owsc.h +++ b/tests/2d_examples/test_2d_owsc/owsc.h @@ -247,14 +247,14 @@ class FlapSystemForSimbody : public SolidBodyPartForSimbody : SolidBodyPartForSimbody(sph_body, shape_ptr) { // Vecd mass_center = Vecd(7.92, 0.355); // 0.3355 - // initial_mass_center_ = SimTK::Vec3(mass_center[0], mass_center[1], 0.0); + // initial_mass_center_ = SimTKVec3(mass_center[0], mass_center[1], 0.0); /** UnitInertia_ (const RealP &xx, const RealP &yy, const RealP &zz) * Create a principal unit inertia matrix (only non-zero on diagonal). */ Real Iz = 1.84 / 33.04; body_part_mass_properties_ = mass_properties_ptr_keeper_ - .createPtr(33.04, SimTK::Vec3(0.0), SimTK::UnitInertia(0.0, 0.0, Iz)); + .createPtr(33.04, SimTKVec3(0.0), SimTK::UnitInertia(0.0, 0.0, Iz)); } }; diff --git a/tests/2d_examples/test_2d_shell_beam_collision/src/shell_beam_collision.cpp b/tests/2d_examples/test_2d_shell_beam_collision/src/shell_beam_collision.cpp index 5560eab898..6fe006e5ab 100644 --- a/tests/2d_examples/test_2d_shell_beam_collision/src/shell_beam_collision.cpp +++ b/tests/2d_examples/test_2d_shell_beam_collision/src/shell_beam_collision.cpp @@ -213,9 +213,9 @@ int main(int ac, char *av[]) SolidBodyPartForSimbody shell_multibody(shell, makeShared("Shell")); SimTK::Body::Rigid rigid_info(*shell_multibody.body_part_mass_properties_); SimTK::MobilizedBody::Slider - shellMBody(matter.Ground(), SimTK::Transform(SimTK::Vec3(0)), rigid_info, SimTK::Transform(SimTK::Vec3(0))); + shellMBody(matter.Ground(), SimTK::Transform(SimTKVec3(0)), rigid_info, SimTK::Transform(SimTKVec3(0))); /** Gravity. */ - SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(Real(-150.), 0.0, 0.0)); + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTKVec3(Real(-150.), 0.0, 0.0)); /** discrete forces acting on the bodies. */ SimTK::Force::DiscreteForces force_on_bodies(forces, matter); /** Time stepping method for multibody system.*/ diff --git a/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp b/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp index 1838e7a34e..5679607c35 100644 --- a/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp +++ b/tests/2d_examples/test_2d_tethered_dead_fish_in_flow/src/tethered_dead_fish_in_flow.cpp @@ -353,18 +353,18 @@ int main(int ac, char *av[]) SimTK::GeneralForceSubsystem forces(MBsystem); SimTK::CableTrackerSubsystem cables(MBsystem); /** Mass properties of the fixed spot. */ - SimTK::Body::Rigid fixed_spot_info(SimTK::MassProperties(1.0, SimTK::Vec3(0), SimTK::UnitInertia(1))); + SimTK::Body::Rigid fixed_spot_info(SimTK::MassProperties(1.0, SimTKVec3(0), SimTK::UnitInertia(1))); SolidBodyPartForSimbody fish_head(fish_body, makeShared(createFishHeadShape(fish_body), "FishHead")); /** Mass properties of the constrained spot. */ SimTK::Body::Rigid tethered_spot_info(*fish_head.body_part_mass_properties_); /** Mobility of the fixed spot. */ - SimTK::MobilizedBody::Weld fixed_spot( matter.Ground(), SimTK::Transform( SimTK::Vec3(tethering_point[0], tethering_point[1], 0.0) ), - fixed_spot_info, SimTK::Transform(SimTK::Vec3(0)) ); + SimTK::MobilizedBody::Weld fixed_spot( matter.Ground(), SimTK::Transform( SimTKVec3(tethering_point[0], tethering_point[1], 0.0) ), + fixed_spot_info, SimTK::Transform(SimTKVec3(0)) ); /** Mobility of the tethered spot. * Set the mass center as the origin location of the planar mobilizer */ Vecd disp0 = fish_head.initial_mass_center_ - tethering_point; - SimTK::MobilizedBody::Planar tethered_spot(fixed_spot, SimTK::Transform(SimTK::Vec3(disp0[0], disp0[1], 0.0)), tethered_spot_info, SimTK::Transform( SimTK::Vec3(0) )); + SimTK::MobilizedBody::Planar tethered_spot(fixed_spot, SimTK::Transform(SimTKVec3(disp0[0], disp0[1], 0.0)), tethered_spot_info, SimTK::Transform( SimTKVec3(0) )); /** The tethering line give cable force. * the start point of the cable path is at the origin location of the first mobilizer body, * the end point is the tip of the fish head which has a distance to the origin @@ -372,7 +372,7 @@ int main(int ac, char *av[]) * of the fish head. */ Vecd disp_cable_end = Vecd(cx, cy) - fish_head.initial_mass_center_; - SimTK::CablePath tethering_line(cables, fixed_spot, SimTK::Vec3(0), tethered_spot, SimTK::Vec3(disp_cable_end[0], disp_cable_end[1], 0.0) ); + SimTK::CablePath tethering_line(cables, fixed_spot, SimTKVec3(0), tethered_spot, SimTKVec3(disp_cable_end[0], disp_cable_end[1], 0.0) ); SimTK::CableSpring tethering_spring(forces, tethering_line, 100.0, 3.0, 10.0); // discrete forces acting on the bodies diff --git a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp index 0dac6aee45..26f5c978f5 100644 --- a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp +++ b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp @@ -6,9 +6,9 @@ void write_csv_files( const std::string& parameter_1_name, const std::string& parameter_2_name, const std::string& parameter_3_name, - const std::vector& parameter_1, - const std::vector& parameter_2, - const std::vector& parameter_3 + const std::vector& parameter_1, + const std::vector& parameter_2, + const std::vector& parameter_3 ) { std::ofstream myfile; diff --git a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.h b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.h index a944ec13b5..a76f44cc9a 100644 --- a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.h +++ b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.h @@ -17,9 +17,9 @@ void write_csv_files( const std::string& parameter_1_name, const std::string& parameter_2_name, const std::string& parameter_3_name, - const std::vector& parameter_1, - const std::vector& parameter_2, - const std::vector& parameter_3 + const std::vector& parameter_1, + const std::vector& parameter_2, + const std::vector& parameter_3 ); @@ -43,7 +43,7 @@ class MeshData }; struct Mesh { - std::vector> vertices; + std::vector> vertices; std::vector> faces; }; diff --git a/tests/3d_examples/test_3d_muscle_solid_contact/src/muscle_solid_contact.cpp b/tests/3d_examples/test_3d_muscle_solid_contact/src/muscle_solid_contact.cpp index e8bb11a035..cc474e1335 100644 --- a/tests/3d_examples/test_3d_muscle_solid_contact/src/muscle_solid_contact.cpp +++ b/tests/3d_examples/test_3d_muscle_solid_contact/src/muscle_solid_contact.cpp @@ -115,9 +115,9 @@ int main() */ SimTK::Body::Rigid rigid_info(*plate_multibody.body_part_mass_properties_); SimTK::MobilizedBody::Slider - plateMBody(matter.Ground(), SimTK::Transform(SimTK::Vec3(0)), rigid_info, SimTK::Transform(SimTK::Vec3(0))); + plateMBody(matter.Ground(), SimTK::Transform(SimTKVec3(0)), rigid_info, SimTK::Transform(SimTKVec3(0))); /** Gravity. */ - SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(Real(-100.0), 0.0, 0.0)); + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTKVec3(Real(-100.0), 0.0, 0.0)); /** discrete forces acting on the bodies. */ SimTK::Force::DiscreteForces force_on_bodies(forces, matter); /** Damper. */ diff --git a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp index 74b31156cc..55f745dcb2 100644 --- a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp +++ b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp @@ -3,7 +3,7 @@ using namespace SPH; -TEST(small_vectors, CalculateDoubleDotProduct) +TEST(small_vectors, CalculateBiDotProduct) { Mat3d matrix_1{ {6.5, 7.8, 0.0}, {1.0, 2.0, 3.0}, @@ -13,7 +13,7 @@ TEST(small_vectors, CalculateDoubleDotProduct) {1.0, 5.0, 7.0}, {0.0, 0.0, 0.2}, }; - EXPECT_DOUBLE_EQ(CalculateDoubleDotProduct(matrix_1, matrix_2), 103.58); + EXPECT_EQ(CalculateBiDotProduct(matrix_1, matrix_2), 103.58); } TEST(small_vectors, getCosineOfAngleBetweenTwoVectors) @@ -87,8 +87,8 @@ TEST(TransformationMatrix, getTransformationMatrix) { // Generate 3D unit vectors with a spherically symmetric probability distribution // according to the link: https://mathworld.wolfram.com/SpherePointPicking.html - Real u = (double)rand() / (RAND_MAX); - Real v = (double)rand() / (RAND_MAX); + Real u = (Real)rand() / (RAND_MAX); + Real v = (Real)rand() / (RAND_MAX); Real theta = 2.0 * Pi * u; Real cos_phi = 2.0 * v - 1.0; Real sin_phi = sqrt(1.0 - cos_phi * cos_phi); @@ -118,7 +118,7 @@ TEST(small_vectors, getPrincipalValuesFromMatrix) }; Vec3d principal_stress_2 = getPrincipalValuesFromMatrix(stress_tensor_2); Vec3d principal_stress_2_ref = {15.3, 12.5, -50.4}; - for (size_t i = 0; i < 3; i++) EXPECT_DOUBLE_EQ(principal_stress_2[i], principal_stress_2_ref[i]); + for (size_t i = 0; i < 3; i++) EXPECT_EQ(principal_stress_2[i], principal_stress_2_ref[i]); } //=================================================================================================// diff --git a/tests/user_examples/test_2d_stfb/stfb.cpp b/tests/user_examples/test_2d_stfb/stfb.cpp index e02ef0965a..46df29064e 100644 --- a/tests/user_examples/test_2d_stfb/stfb.cpp +++ b/tests/user_examples/test_2d_stfb/stfb.cpp @@ -91,12 +91,12 @@ int main(int ac, char *av[]) * and a new child (outboard) body B created by copying the given \a bodyInfo * into a privately-owned Body within the constructed %MobilizedBody object. * Specify the mobilizer frames F fixed to parent P and M fixed to child B. - * @param[in] inboard(SimTK::Vec3) Defines the location of the joint point relative to the parent body. - * @param[in] outboard(SimTK::Vec3) Defines the body's origin location to the joint point. - * @note The body's origin location can be the mass center, the the center of mass should be SimTK::Vec3(0) + * @param[in] inboard(SimTKVec3) Defines the location of the joint point relative to the parent body. + * @param[in] outboard(SimTKVec3) Defines the body's origin location to the joint point. + * @note The body's origin location can be the mass center, the the center of mass should be SimTKVec3(0) * in SimTK::MassProperties(mass, com, inertia) */ - SimTK::MobilizedBody::Planar structure_mob(matter.Ground(), SimTK::Transform(SimTK::Vec3(G[0],G[1], 0.0)), structure_info, SimTK::Transform(SimTK::Vec3(0.0, 0.0, 0.0))); + SimTK::MobilizedBody::Planar structure_mob(matter.Ground(), SimTK::Transform(SimTKVec3(G[0],G[1], 0.0)), structure_info, SimTK::Transform(SimTKVec3(0.0, 0.0, 0.0))); /** * @details Add gravity to mb body. * @param[in,out] forces, The subsystem to which this force should be added. @@ -121,7 +121,7 @@ int main(int ac, char *av[]) * hb=pb*(-d) - hz. Note that this is a signed quantity so the potential energy is * also signed. 0.475 */ - SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(0.0, -gravity_g, 0.0), 0.0); + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTKVec3(0.0, -gravity_g, 0.0), 0.0); /** discrete forces acting on the bodies. */ SimTK::Force::DiscreteForces force_on_bodies(forces, matter); /** Time stepping method for multibody system.*/ diff --git a/tests/user_examples/test_2d_stfb/stfb.h b/tests/user_examples/test_2d_stfb/stfb.h index f647a499d2..ea406273eb 100644 --- a/tests/user_examples/test_2d_stfb/stfb.h +++ b/tests/user_examples/test_2d_stfb/stfb.h @@ -96,10 +96,10 @@ class StructureSystemForSimbody : public SolidBodyPartForSimbody : SolidBodyPartForSimbody(sph_body, shape_ptr) { //Vec2d mass_center(G[0], G[1]); - //initial_mass_center_ = SimTK::Vec3(mass_center[0], mass_center[1], 0.0); + //initial_mass_center_ = SimTKVec3(mass_center[0], mass_center[1], 0.0); body_part_mass_properties_ = mass_properties_ptr_keeper_ - .createPtr(StructureMass, SimTK::Vec3(0.0), SimTK::UnitInertia(Ix,Iy,Iz)); + .createPtr(StructureMass, SimTKVec3(0.0), SimTK::UnitInertia(Ix,Iy,Iz)); } }; //---------------------------------------------------------------------- diff --git a/tests/user_examples/test_3d_stfb/stfb.cpp b/tests/user_examples/test_3d_stfb/stfb.cpp index fc85f8b718..0de45f1847 100644 --- a/tests/user_examples/test_3d_stfb/stfb.cpp +++ b/tests/user_examples/test_3d_stfb/stfb.cpp @@ -99,12 +99,12 @@ int main(int ac, char *av[]) * and a new child (outboard) body B created by copying the given \a bodyInfo * into a privately-owned Body within the constructed %MobilizedBody object. * Specify the mobilizer frames F fixed to parent P and M fixed to child B. - * @param[in] inboard(SimTK::Vec3) Defines the location of the joint point relative to the parent body. - * @param[in] outboard(SimTK::Vec3) Defines the body's origin location to the joint point. - * @note The body's origin location can be the mass center, the the center of mass should be SimTK::Vec3(0) + * @param[in] inboard(SimTKVec3) Defines the location of the joint point relative to the parent body. + * @param[in] outboard(SimTKVec3) Defines the body's origin location to the joint point. + * @note The body's origin location can be the mass center, the the center of mass should be SimTKVec3(0) * in SimTK::MassProperties(mass, com, inertia) */ - SimTK::MobilizedBody::Planar structure_mob(matter.Ground(), SimTK::Transform(SimTK::Vec3(G[0],G[1], G[2])), structure_info, SimTK::Transform(SimTK::Vec3(0.0, 0.0, 0.0))); + SimTK::MobilizedBody::Planar structure_mob(matter.Ground(), SimTK::Transform(SimTKVec3(G[0],G[1], G[2])), structure_info, SimTK::Transform(SimTKVec3(0.0, 0.0, 0.0))); /** * @details Add gravity to mb body. * @param[in,out] forces, The subsystem to which this force should be added. @@ -129,7 +129,7 @@ int main(int ac, char *av[]) * hb=pb*(-d) - hz. Note that this is a signed quantity so the potential energy is * also signed. 0.475 */ - SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTK::Vec3(0.0, 0.0, -gravity_g), 0.0); + SimTK::Force::UniformGravity sim_gravity(forces, matter, SimTKVec3(0.0, 0.0, -gravity_g), 0.0); /** discrete forces acting on the bodies. */ SimTK::Force::DiscreteForces force_on_bodies(forces, matter); /** Time stepping method for multibody system.*/ diff --git a/tests/user_examples/test_3d_stfb/stfb.h b/tests/user_examples/test_3d_stfb/stfb.h index 798f53e875..a1e01e60ff 100644 --- a/tests/user_examples/test_3d_stfb/stfb.h +++ b/tests/user_examples/test_3d_stfb/stfb.h @@ -75,10 +75,10 @@ class StructureSystemForSimbody : public SolidBodyPartForSimbody : SolidBodyPartForSimbody(sph_body, shape_ptr) { //Vec2d mass_center(G[0], G[1]); - //initial_mass_center_ = SimTK::Vec3(mass_center[0], mass_center[1], 0.0); + //initial_mass_center_ = SimTKVec3(mass_center[0], mass_center[1], 0.0); body_part_mass_properties_ = mass_properties_ptr_keeper_ - .createPtr(StructureMass, SimTK::Vec3(0.0), SimTK::UnitInertia(Ix,Iy,Iz)); + .createPtr(StructureMass, SimTKVec3(0.0), SimTK::UnitInertia(Ix,Iy,Iz)); } }; //---------------------------------------------------------------------- diff --git a/tests/webassembly_models/bernoulli_beam/bernoulli_beam.cpp b/tests/webassembly_models/bernoulli_beam/bernoulli_beam.cpp index 35ca0f3e4a..742ce967a8 100755 --- a/tests/webassembly_models/bernoulli_beam/bernoulli_beam.cpp +++ b/tests/webassembly_models/bernoulli_beam/bernoulli_beam.cpp @@ -12,7 +12,7 @@ EMSCRIPTEN_BINDINGS(SPHINXSYS) { - emscripten::value_array>("ArrayDouble3") + emscripten::value_array>("ArrayReal3") .element(emscripten::index<0>()) .element(emscripten::index<1>()) .element(emscripten::index<2>()); @@ -22,7 +22,7 @@ EMSCRIPTEN_BINDINGS(SPHINXSYS) .field("ptr", &StlData::ptr); emscripten::register_vector("StringVector"); - emscripten::register_vector("DoubleVector"); + emscripten::register_vector("BiVectortor"); emscripten::register_vector("UIntVector"); emscripten::register_vector("StlsList"); emscripten::register_vector("BodiesList"); diff --git a/tests/webassembly_models/bernoulli_beam/bernoulli_beam.h b/tests/webassembly_models/bernoulli_beam/bernoulli_beam.h index 440ab24af6..a118574c3a 100755 --- a/tests/webassembly_models/bernoulli_beam/bernoulli_beam.h +++ b/tests/webassembly_models/bernoulli_beam/bernoulli_beam.h @@ -15,13 +15,13 @@ using namespace SPH; struct BernoulliBeamInput { - double scale_stl; - std::vector resolution; - double rho_0; - double poisson; - double Youngs_modulus; - double physical_viscosity; - std::array translation; + Real scale_stl; + std::vector resolution; + Real rho_0; + Real poisson; + Real Youngs_modulus; + Real physical_viscosity; + std::array translation; StlList stls; std::string relative_input_path; }; @@ -68,7 +68,7 @@ class BernoulliBeam } public: //C++ Backend functions - void runCompleteSimulation(double endTime) { sim->runSimulation(SPH::Real(endTime)); }; + void runCompleteSimulation(Real endTime) { sim->runSimulation(SPH::Real(endTime)); }; private: std::SharedPtr material_; From 29312bc95e8271d682092005643026c337e9b04e Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Sat, 3 Jun 2023 12:58:21 +0200 Subject: [PATCH 2/8] single-double continue --- src/shared/adaptations/adaptation.cpp | 6 ++-- src/shared/common/base_data_type.h | 16 +++++------ src/shared/common/vector_functions.cpp | 28 +++++++++---------- src/shared/materials/riemann_solver.cpp | 10 +++---- .../fluid_dynamics/fluid_dynamics_complex.hpp | 2 +- .../solid_dynamics/contact_dynamics.h | 10 +++---- .../fluid_structure_interaction.h | 2 +- .../src/T_shaped_pipe.cpp | 2 +- ..._bending_circular_plate_parametric_cvt.cpp | 2 +- .../test_3d_roof_parametric_cvt.cpp | 2 +- .../fluid_dynamics_complex_wkgc.hpp | 2 +- 11 files changed, 41 insertions(+), 41 deletions(-) diff --git a/src/shared/adaptations/adaptation.cpp b/src/shared/adaptations/adaptation.cpp index 13865877b5..493a7d43a9 100644 --- a/src/shared/adaptations/adaptation.cpp +++ b/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("SmoothingLengthRatio"); base_particles.addVariableToReload("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]}; } //=================================================================================================// } diff --git a/src/shared/common/base_data_type.h b/src/shared/common/base_data_type.h index 6581041639..e1bbbe3878 100644 --- a/src/shared/common/base_data_type.h +++ b/src/shared/common/base_data_type.h @@ -71,17 +71,9 @@ namespace SPH #if ENABLE_SINGLE using Real = float; - using SimTKVec2 = SimTK::fVec2; - using SimTKVec3 = SimTK::fVec3; - using SimTKMat22 = SimTK::fMat22; - using SimTKMat33 = SimTK::fMat33; using EigMat = Eigen::MatrixXf; #else using Real = double; - using SimTKVec2 = SimTK::Vec2; - using SimTKVec3 = SimTK::Vec3; - using SimTKMat22 = SimTK::Mat22; - using SimTKMat33 = SimTK::Mat33; using EigMat = Eigen::MatrixXd; #endif @@ -97,6 +89,14 @@ namespace SPH /** AlignedBox */ using AlignedBox2d = Eigen::AlignedBox; using AlignedBox3d = Eigen::AlignedBox; + /** 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? diff --git a/src/shared/common/vector_functions.cpp b/src/shared/common/vector_functions.cpp index b3694cf549..dbcd58e7b0 100644 --- a/src/shared/common/vector_functions.cpp +++ b/src/shared/common/vector_functions.cpp @@ -5,48 +5,48 @@ namespace SPH //=================================================================================================// SimTKVec2 EigenToSimTK(const Vec2d &eigen_vector) { - return SimTKVec2(eigen_vector[0], eigen_vector[1]); + return SimTKVec2((double)eigen_vector[0], (double)eigen_vector[1]); } //=================================================================================================// SimTKVec3 EigenToSimTK(const Vec3d &eigen_vector) { - return SimTKVec3(eigen_vector[0], eigen_vector[1], eigen_vector[2]); + return SimTKVec3((double)eigen_vector[0], (double)eigen_vector[1], (double)eigen_vector[2]); } //=================================================================================================// Vec2d SimTKToEigen(const SimTKVec2 &simTK_vector) { - return Vec2d(simTK_vector[0], simTK_vector[1]); + return Vec2d((Real)simTK_vector[0], (Real)simTK_vector[1]); } //=================================================================================================// Vec3d SimTKToEigen(const SimTKVec3 &simTK_vector) { - return Vec3d(simTK_vector[0], simTK_vector[1], simTK_vector[2]); + return Vec3d((Real)simTK_vector[0], (Real)simTK_vector[1], (Real)simTK_vector[2]); } SimTKMat22 EigenToSimTK(const Mat2d &eigen_matrix) { - return SimTKMat22(eigen_matrix(0, 0), eigen_matrix(0, 1), - eigen_matrix(1, 0), eigen_matrix(1, 1)); + return SimTKMat22((double)eigen_matrix(0, 0), (double)eigen_matrix(0, 1), + (double)eigen_matrix(1, 0), (double)eigen_matrix(1, 1)); } //=================================================================================================// SimTKMat33 EigenToSimTK(const Mat3d &eigen_matrix) { - return SimTKMat33(eigen_matrix(0, 0), eigen_matrix(0, 1), eigen_matrix(0, 2), - eigen_matrix(1, 0), eigen_matrix(1, 1), eigen_matrix(1, 2), - eigen_matrix(2, 0), eigen_matrix(2, 1), eigen_matrix(2, 2)); + return SimTKMat33((double)eigen_matrix(0, 0), (double)eigen_matrix(0, 1), (double)eigen_matrix(0, 2), + (double)eigen_matrix(1, 0), (double)eigen_matrix(1, 1), (double)eigen_matrix(1, 2), + (double)eigen_matrix(2, 0), (double)eigen_matrix(2, 1), (double)eigen_matrix(2, 2)); } Mat2d SimTKToEigen(const SimTKMat22 &simTK_matrix) { return Mat2d{ - {simTK_matrix(0, 0), simTK_matrix(0, 1)}, - {simTK_matrix(1, 0), simTK_matrix(1, 1)}}; + {(Real)simTK_matrix(0, 0), (Real)simTK_matrix(0, 1)}, + {(Real)simTK_matrix(1, 0), (Real)simTK_matrix(1, 1)}}; } //=================================================================================================// Mat3d SimTKToEigen(const SimTKMat33 &simTK_matrix) { return Mat3d{ - {simTK_matrix(0, 0), simTK_matrix(0, 1), simTK_matrix(0, 2)}, - {simTK_matrix(1, 0), simTK_matrix(1, 1), simTK_matrix(1, 2)}, - {simTK_matrix(2, 0), simTK_matrix(2, 1), simTK_matrix(2, 2)}}; + {(Real)simTK_matrix(0, 0), (Real)simTK_matrix(0, 1), (Real)simTK_matrix(0, 2)}, + {(Real)simTK_matrix(1, 0), (Real)simTK_matrix(1, 1), (Real)simTK_matrix(1, 2)}, + {(Real)simTK_matrix(2, 0), (Real)simTK_matrix(2, 1), (Real)simTK_matrix(2, 2)}}; } //=================================================================================================// Vec2d FirstAxisVector(const Vec2d &zero_vector) diff --git a/src/shared/materials/riemann_solver.cpp b/src/shared/materials/riemann_solver.cpp index 8407774826..c419aeb4e8 100644 --- a/src/shared/materials/riemann_solver.cpp +++ b/src/shared/materials/riemann_solver.cpp @@ -27,7 +27,7 @@ namespace SPH //=================================================================================================// Real AcousticRiemannSolver::DissipativePJump(const Real &u_jump) { - return rho0c0_geo_ave_ * u_jump * SMIN(3.0 * SMAX(u_jump * inv_c_ave_, 0.0), 1.0); + return rho0c0_geo_ave_ * u_jump * SMIN( Real(3) * SMAX(u_jump * inv_c_ave_, Real(0)), Real(1)); } //=================================================================================================// Real AcousticRiemannSolver::DissipativeUJump(const Real &p_jump) @@ -100,7 +100,7 @@ namespace SPH Real rho_cl = state_i.rho_ * fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); Real rho_cr = state_j.rho_ * fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); Real rho_clr = (rho_cl * state_i.rho_ + rho_cr * state_j.rho_) / (state_i.rho_ + state_j.rho_); - p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(3.0 * SMAX(rho_ave * (ul - ur), 0.0), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); + p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(Real(3) * SMAX(rho_ave * (ul - ur), Real(0)), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); v_star = state_i.vel_ - e_ij * (s_star - ul); } if (s_star <= 0.0 && 0.0 <= s_r) @@ -109,7 +109,7 @@ namespace SPH Real rho_cl = state_i.rho_ * fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); Real rho_cr = state_j.rho_ * fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); Real rho_clr = (rho_cl * state_i.rho_ + rho_cr * state_j.rho_) / (state_i.rho_ + state_j.rho_); - p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(3.0 * SMAX(rho_ave * (ul - ur), 0.0), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); + p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(Real(3) * SMAX(rho_ave * (ul - ur), Real(0)), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); v_star = state_j.vel_ - e_ij * (s_star - ur); } if (s_r < 0.0) @@ -192,7 +192,7 @@ namespace SPH Real rho_cl = state_i.rho_ * compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); Real rho_cr = state_j.rho_ * compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); Real rho_clr = (rho_cl * state_i.rho_ + rho_cr * state_j.rho_) / (state_i.rho_ + state_j.rho_); - p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(3.0 * SMAX(rho_ave * (ul - ur), 0.0), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); + p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(Real(3) * SMAX(rho_ave * (ul - ur), Real(0)), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); v_star = state_i.vel_ - e_ij * (s_star - ul); rho_star = state_i.rho_ * (s_l - ul) / (s_l - s_star); energy_star = state_i.rho_ * (s_l - ul) / (s_l - s_star) * (state_i.E_ / state_i.rho_ + (s_star - ul) * (s_star + state_i.p_ / state_i.rho_ / (s_l - ul))); @@ -203,7 +203,7 @@ namespace SPH Real rho_cl = state_i.rho_ * compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_); Real rho_cr = state_j.rho_ * compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_); Real rho_clr = (rho_cl * state_i.rho_ + rho_cr * state_j.rho_) / (state_i.rho_ + state_j.rho_); - p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(3.0 * SMAX(rho_ave * (ul - ur), 0.0), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); + p_star = 0.5 * (state_i.p_ + state_j.p_) + 0.5 * (SMIN(Real(3) * SMAX(rho_ave * (ul - ur), Real(0)), rho_clr) * (ul - ur) + s_star * (rho_cr - rho_cl)); v_star = state_j.vel_ - e_ij * (s_star - ur); rho_star = state_j.rho_ * (s_r - ur) / (s_r - s_star); energy_star = state_j.rho_ * (s_r - ur) / (s_r - s_star) * (state_j.E_ / state_j.rho_ + (s_star - ur) * (s_star + state_j.p_ / state_j.rho_ / (s_r - ur))); diff --git a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp index dbc50656a2..c421d97e87 100644 --- a/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp +++ b/src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_complex.hpp @@ -259,7 +259,7 @@ namespace SPH Real r_ij = wall_neighborhood.r_ij_[n]; Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); - Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_wall_external_acceleration); + Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(Real(0), face_wall_external_acceleration); acceleration -= (this->p_[index_i] + p_in_wall) * dW_ijV_j * e_ij; rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; } diff --git a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h index be6f5a2e47..92aba254b4 100644 --- a/src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h @@ -66,7 +66,7 @@ namespace SPH const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { - Real corrected_W_ij = std::max(inner_neighborhood.W_ij_[n] - offset_W_ij_, 0.0); + Real corrected_W_ij = std::max(inner_neighborhood.W_ij_[n] - offset_W_ij_, Real(0)); sigma += corrected_W_ij * mass_[inner_neighborhood.j_[n]]; } contact_density_[index_i] = sigma; @@ -98,7 +98,7 @@ namespace SPH for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { - Real corrected_W_ij = std::max(contact_neighborhood.W_ij_[n] - offset_W_ij_[k], 0.0); + Real corrected_W_ij = std::max(contact_neighborhood.W_ij_[n] - offset_W_ij_[k], Real(0)); sigma += corrected_W_ij * contact_mass_k[contact_neighborhood.j_[n]]; } } @@ -134,7 +134,7 @@ namespace SPH Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { - Real corrected_W_ij = std::max(contact_neighborhood.W_ij_[n] - offset_W_ij_[k], 0.0); + Real corrected_W_ij = std::max(contact_neighborhood.W_ij_[n] - offset_W_ij_[k], Real(0)); sigma += corrected_W_ij * contact_Vol_k[contact_neighborhood.j_[n]]; } constexpr Real heuristic_limiter = 0.1; @@ -359,7 +359,7 @@ namespace SPH // only update particle i Vecd vel_derivative = (vel_i - vel_k[index_j]); Vecd n_j = e_ij.dot(n_k[index_j]) > 0.0 ? n_k[index_j] : -1.0 * n_k[index_j]; - vel_derivative -= SMAX(0.0, vel_derivative.dot(n_j)) * n_j; + vel_derivative -= SMAX(Real(0), vel_derivative.dot(n_j)) * n_j; vel_i += parameter_b[n] * vel_derivative / (mass_i - 2.0 * parameter_b[n]); } // backward sweep @@ -371,7 +371,7 @@ namespace SPH // only update particle i Vecd vel_derivative = (vel_i - vel_k[index_j]); Vecd n_j = e_ij.dot(n_k[index_j]) > 0.0 ? n_k[index_j] : -1.0 * n_k[index_j]; - vel_derivative -= SMAX(0.0, vel_derivative.dot(n_j)) * n_j; + vel_derivative -= SMAX(Real(0), vel_derivative.dot(n_j)) * n_j; vel_i += parameter_b[n - 1] * vel_derivative / (mass_i - 2.0 * parameter_b[n - 1]); } } diff --git a/src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h b/src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h index 7b78bfb074..461ff74e4a 100644 --- a/src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h +++ b/src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h @@ -138,7 +138,7 @@ namespace SPH Vecd e_ij = contact_neighborhood.e_ij_[n]; Real r_ij = contact_neighborhood.r_ij_[n]; Real face_wall_external_acceleration = (acc_prior_k[index_j] - acc_ave_[index_i]).dot(e_ij); - Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(0.0, face_wall_external_acceleration); + Real p_in_wall = p_k[index_j] + rho_n_k[index_j] * r_ij * SMAX(Real(0), face_wall_external_acceleration); Real u_jump = 2.0 * (vel_k[index_j] - vel_ave_[index_i]).dot(n_[index_i]); force += (riemann_solvers_k.DissipativePJump(u_jump) * n_[index_i] - (p_in_wall + p_k[index_j]) * e_ij ) * Vol_[index_i] * contact_neighborhood.dW_ijV_j_[n]; diff --git a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp index 9c213d315c..d0703b1d7f 100644 --- a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp +++ b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp @@ -23,7 +23,7 @@ Real DL_sponge = resolution_ref * 20; /**< Reference size of the emitter buffer Real rho0_f = 1.0; /**< Reference density of fluid. */ Real U_f = 1.0; /**< Characteristic velocity. */ /** Reference sound speed needs to consider the flow speed in the narrow channels. */ -Real c_f = 10.0 * U_f * SMAX(1.0, DH / (2.0 * (DL - DL1))); +Real c_f = 10.0 * U_f * SMAX(Real(1), DH / (2.0 * (DL - DL1))); Real Re = 100.0; /**< Reynolds number. */ Real mu_f = rho0_f * U_f * DH / Re; /**< Dynamics viscosity. */ //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_bending_circular_plate_parametric_cvt/test_3d_bending_circular_plate_parametric_cvt.cpp b/tests/3d_examples/test_3d_bending_circular_plate_parametric_cvt/test_3d_bending_circular_plate_parametric_cvt.cpp index b60d234f9f..5e2bc882b1 100644 --- a/tests/3d_examples/test_3d_bending_circular_plate_parametric_cvt/test_3d_bending_circular_plate_parametric_cvt.cpp +++ b/tests/3d_examples/test_3d_bending_circular_plate_parametric_cvt/test_3d_bending_circular_plate_parametric_cvt.cpp @@ -320,7 +320,7 @@ return_data bending_circular_plate(Real dp_ratio) initialize_external_force.exec(dt); - dt = std::min(thickness/dp, 0.5) * computing_time_step_size.exec(); + dt = std::min(thickness/dp, Real(0.5)) * computing_time_step_size.exec(); {// checking for excessive time step reduction if (dt > max_dt) max_dt = dt; if (dt < max_dt/1e3) throw std::runtime_error("time step decreased too much"); diff --git a/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp b/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp index d0591c8a85..b089355619 100644 --- a/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp +++ b/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp @@ -421,7 +421,7 @@ return_data roof_under_self_weight(Real dp, bool cvt = true, int particle_number initialize_external_force.exec(dt); - dt = std::min(thickness/dp, 0.5) * computing_time_step_size.exec(); + dt = std::min(thickness/dp, Real(0.5)) * computing_time_step_size.exec(); {// checking for excessive time step reduction if (dt > max_dt) max_dt = dt; if (dt < max_dt/1e3) throw std::runtime_error("time step decreased too much"); diff --git a/tests/user_examples/ExtraSources/extra_sources_share/fluid_dynamics_complex_wkgc.hpp b/tests/user_examples/ExtraSources/extra_sources_share/fluid_dynamics_complex_wkgc.hpp index 6c8d1f9d4f..d572ae7d34 100644 --- a/tests/user_examples/ExtraSources/extra_sources_share/fluid_dynamics_complex_wkgc.hpp +++ b/tests/user_examples/ExtraSources/extra_sources_share/fluid_dynamics_complex_wkgc.hpp @@ -51,7 +51,7 @@ namespace SPH Real r_ij = wall_neighborhood.r_ij_[n]; Real face_wall_external_acceleration = (acc_prior_i - acc_ave_k[index_j]).dot(-e_ij); - Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(0.0, face_wall_external_acceleration); + Real p_in_wall = this->p_[index_i] + this->rho_[index_i] * r_ij * SMAX(Real(0), face_wall_external_acceleration); acceleration -= (this->p_[index_i] + p_in_wall) * this->B_[index_i] * e_ij * dW_ijV_j; rho_dissipation += this->riemann_solver_.DissipativeUJump(this->p_[index_i] - p_in_wall) * dW_ijV_j; } From 0bdaf82760643e488b1590e8f0fd438bdd1fe3a8 Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Sun, 4 Jun 2023 10:30:41 +0200 Subject: [PATCH 3/8] 2D compilation pass --- CMakeLists.txt | 6 +++--- src/shared/common/base_data_type.h | 2 +- .../general_dynamics/general_dynamics_refinement.cpp | 2 +- .../general_dynamics/general_dynamics_refinement.h | 2 +- .../particle_dynamics/solid_dynamics/loading_dynamics.cpp | 2 +- 5 files changed, 7 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 98851bfd9c..2896b3dad6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +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_SINGLE "Build using single precision or not" OFF) +option(SPHINXSYS_USE_SINGLE "Build using single precision 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 @@ -79,8 +79,8 @@ 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_SINGLE) - add_definitions(-D ENABLE_SINGLE) +if(SPHINXSYS_USE_SINGLE) + add_definitions(-D USE_SINGLE) endif() if(NOT BOOST_INCLUDE_DIR) diff --git a/src/shared/common/base_data_type.h b/src/shared/common/base_data_type.h index e1bbbe3878..16ba6259dc 100644 --- a/src/shared/common/base_data_type.h +++ b/src/shared/common/base_data_type.h @@ -69,7 +69,7 @@ namespace SPH * Using -mtune=native produces code optimized for the local machine under the constraints of the selected instruction set. */ -#if ENABLE_SINGLE +#if USE_SINGLE using Real = float; using EigMat = Eigen::MatrixXf; #else diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp index 02a945fd84..2cc8fe6f0a 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp @@ -440,7 +440,7 @@ ParticleMergeWithPrescribedArea:: DataDelegateInner(inner_relation), all_particle_data_(particles_->all_particle_data_), vel_n_(particles_->vel_) { - particles_->registerVariable(total_merge_error_, "MergeDensityError", 0.0); + particles_->registerVariable(total_merge_error_, "MergeDensityError", Real(0)); } //=================================================================================================// void ParticleMergeWithPrescribedArea::setupDynamics(Real dt) diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.h b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.h index 3e25862303..068370644c 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.h +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.h @@ -180,7 +180,7 @@ class SplitWithMinimumDensityErrorInner : public ParticleSplitWithPrescribedArea : ParticleSplitWithPrescribedArea(inner_relation.getSPHBody(), refinement_region, body_buffer_width), compute_density_error(inner_relation) { - particles_->registerVariable(total_split_error_, "SplitDensityError", 0.0); + particles_->registerVariable(total_split_error_, "SplitDensityError", Real(0)); }; virtual ~SplitWithMinimumDensityErrorInner(){}; diff --git a/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp index e96e13f637..475d6b17dc 100644 --- a/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/loading_dynamics.cpp @@ -175,7 +175,7 @@ namespace SPH //=================================================================================================// void ForceInBodyRegion::update(size_t index_i, Real dt) { - Real time_factor = SMIN(GlobalStaticVariables::physical_time_ / end_time_, 1.0); + Real time_factor = SMIN(GlobalStaticVariables::physical_time_ / end_time_, Real(1.0)); acc_prior_[index_i] = acceleration_ * time_factor; } //=================================================================================================// From 7fd9f03c3811afc87235b7a0d38be2acfc43437d Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Sun, 4 Jun 2023 11:26:09 +0200 Subject: [PATCH 4/8] All compilation for single precision success --- .../geometries/triangle_mesh_shape.cpp | 8 ++--- .../ensemble_averaged_method.hpp | 30 +++++++++++-------- .../regression_test/time_averaged_method.hpp | 21 +++++++++---- .../src/T_shaped_pipe.cpp | 2 +- .../beam_pulling_pressure_load.cpp | 8 ++--- .../excitation-contraction.cpp | 4 +-- .../excitation_contraction.h | 4 +-- .../heart_volume_change.cpp | 2 +- .../pkj_lv_electrocontraction.h | 4 +-- .../bernoulli_beam_struct_sim.cpp | 6 ++-- 10 files changed, 52 insertions(+), 37 deletions(-) diff --git a/src/for_3D_build/geometries/triangle_mesh_shape.cpp b/src/for_3D_build/geometries/triangle_mesh_shape.cpp index 6be47521b8..7b8dd18e76 100644 --- a/src/for_3D_build/geometries/triangle_mesh_shape.cpp +++ b/src/for_3D_build/geometries/triangle_mesh_shape.cpp @@ -116,10 +116,10 @@ namespace SPH polymesh.loadStlFile(filepathname); polymesh.scaleMesh(scale_factor); - SimTK::Transform_ transform(SimTK::Rotation_(SimTKMat33(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))), - SimTKVec3(translation[0], translation[1], translation[2])); + SimTK::Transform_ transform(SimTK::Rotation_(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)); } //=================================================================================================// diff --git a/src/shared/regression_test/ensemble_averaged_method.hpp b/src/shared/regression_test/ensemble_averaged_method.hpp index d373066d1c..7ee68ec6cc 100644 --- a/src/shared/regression_test/ensemble_averaged_method.hpp +++ b/src/shared/regression_test/ensemble_averaged_method.hpp @@ -42,10 +42,12 @@ namespace SPH for (int observation_index = 0; observation_index != this->observation_; ++observation_index) for (int run_index = 0; run_index != this->number_of_run_; ++run_index) { - variance_new[snapshot_index][observation_index] = SMAX(variance[snapshot_index][observation_index], - variance_new[snapshot_index][observation_index], - pow((result[run_index][snapshot_index][observation_index] - meanvalue_new[snapshot_index][observation_index]), 2), - pow(meanvalue_new[snapshot_index][observation_index] * 1.0e-2, 2) ); + variance_new[snapshot_index][observation_index] = SMAX( + (Real)variance[snapshot_index][observation_index], + (Real)variance_new[snapshot_index][observation_index], + (Real)pow((result[run_index][snapshot_index][observation_index] - meanvalue_new[snapshot_index][observation_index]), 2), + (Real)pow(meanvalue_new[snapshot_index][observation_index] * 1.0e-2, 2) + ); } }; //=================================================================================================// @@ -58,10 +60,12 @@ namespace SPH for (int run_index = 0; run_index != this->number_of_run_; ++run_index) for (int i = 0; i != variance[0][0].size(); ++i) { - variance_new[snapshot_index][observation_index][i] = SMAX(variance[snapshot_index][observation_index][i], - variance_new[snapshot_index][observation_index][i], - pow((result[run_index][snapshot_index][observation_index][i] - meanvalue_new[snapshot_index][observation_index][i]), 2), - pow(meanvalue_new[snapshot_index][observation_index][i] * 1.0e-2, 2) ); + variance_new[snapshot_index][observation_index][i] = SMAX( + (Real)variance[snapshot_index][observation_index][i], + (Real)variance_new[snapshot_index][observation_index][i], + (Real)pow((result[run_index][snapshot_index][observation_index][i] - meanvalue_new[snapshot_index][observation_index][i]), 2), + (Real)pow(meanvalue_new[snapshot_index][observation_index][i] * 1.0e-2, 2) + ); } }; //=================================================================================================// @@ -75,10 +79,12 @@ namespace SPH for (size_t i = 0; i != variance[0][0].size(); ++i) for (size_t j = 0; j != variance[0][0].size(); ++j) { - variance_new[snapshot_index][observation_index](i,j) = SMAX(variance[snapshot_index][observation_index](i,j), - variance_new[snapshot_index][observation_index](i,j), - pow((result[run_index][snapshot_index][observation_index](i,j) - meanvalue_new[snapshot_index][observation_index](i,j)), 2), - pow(meanvalue_new[snapshot_index][observation_index](i,j) * 1.0e-2, 2)); + variance_new[snapshot_index][observation_index](i,j) = SMAX( + (Real)variance[snapshot_index][observation_index](i,j), + (Real)variance_new[snapshot_index][observation_index](i,j), + (Real)pow((result[run_index][snapshot_index][observation_index](i,j) - meanvalue_new[snapshot_index][observation_index](i,j)), 2), + (Real)pow(meanvalue_new[snapshot_index][observation_index](i,j) * Real(0.01), 2) + ); } }; //=================================================================================================// diff --git a/src/shared/regression_test/time_averaged_method.hpp b/src/shared/regression_test/time_averaged_method.hpp index 6a408317dd..6354956e60 100644 --- a/src/shared/regression_test/time_averaged_method.hpp +++ b/src/shared/regression_test/time_averaged_method.hpp @@ -224,8 +224,11 @@ namespace SPH { for (int snapshot_index = snapshot_for_converged_; snapshot_index != this->snapshot_; ++snapshot_index) variance_new[observation_index] += pow((current_result[observation_index][snapshot_index] - local_meanvalue[observation_index]), 2); - variance_new[observation_index] = SMAX( (variance_new[observation_index] / (this->snapshot_ - snapshot_for_converged_)), - variance[observation_index], pow(local_meanvalue[observation_index] * 1.0e-2, 2) ); + variance_new[observation_index] = SMAX( + (Real)(variance_new[observation_index] / (this->snapshot_ - snapshot_for_converged_)), + (Real)variance[observation_index], + (Real)pow(local_meanvalue[observation_index] * Real(0.01), 2) + ); } } //=================================================================================================// @@ -238,8 +241,11 @@ namespace SPH { for (int snapshot_index = snapshot_for_converged_; snapshot_index != this->snapshot_; ++snapshot_index) variance_new[observation_index][i] += pow((current_result[observation_index][snapshot_index][i] - local_meanvalue[observation_index][i]), 2); - variance_new[observation_index][i] = SMAX( (variance_new[observation_index][i] / (this->snapshot_ - snapshot_for_converged_)), - variance[observation_index][i], pow(local_meanvalue[observation_index][i] * 1.0e-2, 2) ); + variance_new[observation_index][i] = SMAX( + (Real)(variance_new[observation_index][i] / (this->snapshot_ - snapshot_for_converged_)), + (Real)variance[observation_index][i], + (Real)pow(local_meanvalue[observation_index][i] * Real(0.01), 2) + ); } } //=================================================================================================// @@ -253,8 +259,11 @@ namespace SPH { for (int snapshot_index = snapshot_for_converged_; snapshot_index != this->snapshot_; ++snapshot_index) variance_new[observation_index](i,j) += pow((current_result[observation_index][snapshot_index](i,j) - local_meanvalue[observation_index](i,j)), 2); - variance_new[observation_index](i,j) = SMAX( (variance_new[observation_index](i,j) / (this->snapshot_ - snapshot_for_converged_)), - variance[observation_index](i,j), pow(local_meanvalue[observation_index](i,j) * 1.0e-2, 2) ); + variance_new[observation_index](i,j) = SMAX( + (Real)(variance_new[observation_index](i,j) / (this->snapshot_ - snapshot_for_converged_)), + (Real)variance[observation_index](i,j), + (Real)pow(local_meanvalue[observation_index](i,j) * Real(0.01), 2) + ); } } //=================================================================================================// diff --git a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp index d0703b1d7f..b7543647b3 100644 --- a/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp +++ b/tests/2d_examples/test_2d_T_shaped_pipe/src/T_shaped_pipe.cpp @@ -23,7 +23,7 @@ Real DL_sponge = resolution_ref * 20; /**< Reference size of the emitter buffer Real rho0_f = 1.0; /**< Reference density of fluid. */ Real U_f = 1.0; /**< Characteristic velocity. */ /** Reference sound speed needs to consider the flow speed in the narrow channels. */ -Real c_f = 10.0 * U_f * SMAX(Real(1), DH / (2.0 * (DL - DL1))); +Real c_f = 10.0 * U_f * SMAX(Real(1), DH / (Real(2.0) * (DL - DL1))); Real Re = 100.0; /**< Reynolds number. */ Real mu_f = rho0_f * U_f * DH / Re; /**< Dynamics viscosity. */ //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp b/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp index c0068b6f83..64982c6e3b 100644 --- a/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp +++ b/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp @@ -151,10 +151,10 @@ int main(int ac, char *av[]) Vecd half_size_0(0.03, 0.03, resolution_ref); BodyRegionByParticle load_surface(beam_body, makeShared(half_size_0, 1, Vecd(0.00, 0.00, 0.1))); StdVec> force_over_time = { - {0.0, 0.0}, - {0.1 * end_time, 0.1 * load_total_force}, - {0.4 * end_time, load_total_force}, - {end_time, load_total_force}}; + {Real(0), Real(0)}, + {Real(0.1) * end_time, Real(0.1) * load_total_force}, + {Real(0.4) * end_time, load_total_force}, + {Real(end_time), Real(load_total_force)}}; SimpleDynamics pull_force(load_surface, force_over_time); std::cout << "load surface particle number: " << load_surface.body_part_particles_.size() << std::endl; diff --git a/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp b/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp index ec8ccdf626..15097dab5d 100644 --- a/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp +++ b/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp @@ -29,8 +29,8 @@ BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); Real rho0_s = 1.06e-3; /** Active stress factor */ Real k_a = 100 * stress_scale; -Real a0[4] = {496.0 * stress_scale, 15196.0 * stress_scale, 3283.0 * stress_scale, 662.0 * stress_scale}; -Real b0[4] = {7.209, 20.417, 11.176, 9.466}; +Real a0[4] = {Real(496.0) * stress_scale, Real(15196.0) * stress_scale, Real(3283.0) * stress_scale, Real(662.0) * stress_scale}; +Real b0[4] = {Real(7.209), Real(20.417), Real(11.176), Real(9.466)}; /** reference stress to achieve weakly compressible condition */ Real poisson = 0.4995; Real bulk_modulus = 2.0 * a0[0] * (1.0 + poisson) / (3.0 * (1.0 - 2.0 * poisson)); diff --git a/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.h b/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.h index 54dac4cc35..c41dce4119 100644 --- a/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.h +++ b/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.h @@ -32,8 +32,8 @@ BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); Real rho0_s = 1.06e-3; /** Active stress factor */ Real k_a = 100 * stress_scale; -Real a0[4] = {496.0 * stress_scale, 15196.0 * stress_scale, 3283.0 * stress_scale, 662.0 * stress_scale}; -Real b0[4] = {7.209, 20.417, 11.176, 9.466}; +Real a0[4] = {Real(496.0) * stress_scale, Real(15196.0) * stress_scale, Real(3283.0) * stress_scale, Real(662.0) * stress_scale}; +Real b0[4] = {Real(7.209), Real(20.417), Real(11.176), Real(9.466)}; /** reference stress to achieve weakly compressible condition */ Real poisson = 0.4995; Real bulk_modulus = 2.0 * a0[0] * (1.0 + poisson) / (3.0 * (1.0 - 2.0 * poisson)); diff --git a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp index 26f5c978f5..ff64108112 100644 --- a/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp +++ b/tests/3d_examples/test_3d_heart_volume_change/heart_volume_change.cpp @@ -30,7 +30,7 @@ void MeshData::load(std::string path_to_mesh, Real scale) for(int i=0;i{ ConstrainedRegionPair(0, fixation) }; StdVec> pressure_over_time = { - {0.0, 0.0}, - {end_time * 0.1, pressure}, - {end_time, pressure } + {Real(0), Real(0)}, + {Real(end_time * 0.1), Real(pressure)}, + {Real(end_time), Real(pressure)} }; input.surface_pressure_tuple_ = StdVec{ PressureTuple(0, specimen, Vec3d(0.1, 0.0, 0.1), pressure_over_time) }; input.particle_relaxation_list_ = { true }; From cf0d25643b7b89dfb87e542ff626bdb3c285895c Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Sun, 4 Jun 2023 17:33:33 +0200 Subject: [PATCH 5/8] Small modification --- CMakeLists.txt | 6 +++--- src/shared/common/base_data_type.h | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2896b3dad6..89a1db8a87 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,7 +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_SINGLE "Build using single precision or not" OFF) +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 @@ -79,8 +79,8 @@ 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_SINGLE) - add_definitions(-D USE_SINGLE) +if(SPHINXSYS_USE_FLOAT) + target_compile_definitions(sphinxsys_core INTERFACE USE_FLOAT) endif() if(NOT BOOST_INCLUDE_DIR) diff --git a/src/shared/common/base_data_type.h b/src/shared/common/base_data_type.h index 16ba6259dc..3691f54778 100644 --- a/src/shared/common/base_data_type.h +++ b/src/shared/common/base_data_type.h @@ -69,7 +69,7 @@ namespace SPH * Using -mtune=native produces code optimized for the local machine under the constraints of the selected instruction set. */ -#if USE_SINGLE +#if USE_FLOAT using Real = float; using EigMat = Eigen::MatrixXf; #else From 61894b613c31d810750a814cb0c0516d003e0f7b Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Tue, 6 Jun 2023 23:18:26 +0200 Subject: [PATCH 6/8] Fix compilation error of clang++. --- .../particle_generator_lattice_supplementary.cpp | 2 +- src/for_3D_build/geometries/triangle_mesh_shape.cpp | 2 +- .../particle_generator_lattice_supplementary.cpp | 2 +- .../particle_generator/particle_generator_network.cpp | 2 +- src/shared/common/base_data_type.h | 4 ++-- src/shared/geometries/base_geometry.cpp | 2 +- src/shared/geometries/level_set.cpp | 2 +- .../dissipation_dynamics/particle_dynamics_dissipation.hpp | 2 +- .../particle_dynamics/general_dynamics/general_dynamics.cpp | 2 +- .../general_dynamics/general_dynamics_refinement.cpp | 2 +- src/shared/particle_generator/particle_generator_lattice.cpp | 2 +- .../common/test_vector_functions/test_vector_functions.cpp | 4 ++-- .../test_thin_structure_dynamics.cpp | 2 +- 13 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp b/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp index b78f797382..15fe13afe5 100644 --- a/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp +++ b/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp @@ -60,7 +60,7 @@ namespace SPH { { if (body_shape_.checkContain(particle_position)) { - Real random_real = (Real)rand() / (RAND_MAX); + Real random_real = Real(rand() / RAND_MAX); // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_) { diff --git a/src/for_3D_build/geometries/triangle_mesh_shape.cpp b/src/for_3D_build/geometries/triangle_mesh_shape.cpp index 7b8dd18e76..d3f8b5cdf1 100644 --- a/src/for_3D_build/geometries/triangle_mesh_shape.cpp +++ b/src/for_3D_build/geometries/triangle_mesh_shape.cpp @@ -36,7 +36,7 @@ namespace SPH { Vec3d jittered = probe_point; // jittering for (int l = 0; l != probe_point.size(); ++l) - jittered[l] = probe_point[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * (SqrtEps + distance_to_pnt * 0.1); + jittered[l] = probe_point[l] + (Real(rand() / RAND_MAX) - 0.5) * (SqrtEps + distance_to_pnt * 0.1); Vec3d from_face_to_jittered = jittered - SimTKToEigen(closest_pnt); Vec3d direction_to_jittered = from_face_to_jittered / (from_face_to_jittered.norm() + TinyReal); cosine_angle = face_normal.dot(direction_to_jittered); diff --git a/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp b/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp index 757ec5ea68..b15dd639b9 100644 --- a/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp +++ b/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp @@ -66,7 +66,7 @@ namespace SPH { if (body_shape_.checkContain(particle_position)) { - Real random_real = (Real)rand() / (RAND_MAX); + Real random_real = Real(rand() / RAND_MAX); // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles. if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_) { diff --git a/src/for_3D_build/particle_generator/particle_generator_network.cpp b/src/for_3D_build/particle_generator/particle_generator_network.cpp index e9e1a2fc45..614773cdb1 100644 --- a/src/for_3D_build/particle_generator/particle_generator_network.cpp +++ b/src/for_3D_build/particle_generator/particle_generator_network.cpp @@ -205,7 +205,7 @@ namespace SPH for (size_t j = 0; j != branches_to_grow.size(); j++) { size_t grow_id = branches_to_grow[j]; - Real rand_num = ((Real)rand() / (RAND_MAX)) - 0.5; + Real rand_num = Real(rand() / RAND_MAX) - 0.5; Real angle_to_use = angle_ + rand_num * 0.05; for (size_t k = 0; k != 2; k++) { diff --git a/src/shared/common/base_data_type.h b/src/shared/common/base_data_type.h index 3691f54778..880c476d3c 100644 --- a/src/shared/common/base_data_type.h +++ b/src/shared/common/base_data_type.h @@ -168,9 +168,9 @@ namespace SPH /** Constant parameters. */ constexpr Real Pi = Real(M_PI); constexpr Real Eps = std::numeric_limits::epsilon(); - constexpr Real TinyReal = pow(Eps, 1.25); + constexpr Real TinyReal = Real(2.71051e-20); constexpr Real Infinity = std::numeric_limits::max(); - constexpr Real SqrtEps = sqrt(Eps); + constexpr Real SqrtEps = Real(1.0e-8); } #endif // BASE_DATA_TYPE_H diff --git a/src/shared/geometries/base_geometry.cpp b/src/shared/geometries/base_geometry.cpp index 16ffb4a991..5e9ee9a60f 100644 --- a/src/shared/geometries/base_geometry.cpp +++ b/src/shared/geometries/base_geometry.cpp @@ -37,7 +37,7 @@ namespace SPH { Vecd jittered = probe_point; for (int l = 0; l != probe_point.size(); ++l) - jittered[l] = probe_point[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * 100.0 * Eps; + jittered[l] = probe_point[l] + (Real(rand() / RAND_MAX) - 0.5) * 100.0 * Eps; if (checkContain(jittered) == is_contain) displacement_to_surface = findClosestPoint(jittered) - jittered; } diff --git a/src/shared/geometries/level_set.cpp b/src/shared/geometries/level_set.cpp index 0a1f01a8fb..4d57f69d23 100644 --- a/src/shared/geometries/level_set.cpp +++ b/src/shared/geometries/level_set.cpp @@ -82,7 +82,7 @@ namespace SPH { Vecd jittered = position; // jittering for (int l = 0; l != position.size(); ++l) - jittered[l] += (((Real)rand() / (RAND_MAX)) - 0.5) * 0.5 * data_spacing_; + jittered[l] += (Real(rand() / RAND_MAX) - 0.5) * 0.5 * data_spacing_; probed_value = probeLevelSetGradient(jittered); } return probed_value.normalized(); diff --git a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp index c102cf3ade..33338830aa 100644 --- a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp +++ b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp @@ -436,7 +436,7 @@ namespace SPH template bool DampingWithRandomChoice::RandomChoice() { - return ((Real)rand() / (RAND_MAX)) < random_ratio_ ? true : false; + return Real(rand() / RAND_MAX) < random_ratio_ ? true : false; } //=================================================================================================// template diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp index 94ec91997a..7b5a8f4106 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp @@ -22,7 +22,7 @@ namespace SPH Vecd &pos_n_i = pos_[index_i]; for (int k = 0; k < pos_n_i.size(); ++k) { - pos_n_i[k] += dt * (((Real)rand() / (RAND_MAX)) - 0.5) * 2.0 * randomize_scale_; + pos_n_i[k] += dt * (Real(rand() / RAND_MAX) - 0.5) * 2.0 * randomize_scale_; } } //=================================================================================================// diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp index 2cc8fe6f0a..f26ed68408 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp @@ -376,7 +376,7 @@ void ParticleSplitWithPrescribedArea::update(size_t index_i, Real dt) Vecd ParticleSplitWithPrescribedArea::getSplittingPosition(const StdVec &new_indices) { srand(int(new_indices[0])); - Real delta_random = 0 + 2.0 * Pi * rand() / RAND_MAX * (2.0 * Pi - 0); + Real delta_random = 0 + 2.0 * Pi * Real(rand() / RAND_MAX) * (2.0 * Pi - 0); Vecd pos_ = particles_->pos_[new_indices[0]]; Real delta = delta_random + Pi; Real Vol_split = Vol_[new_indices[0]] / 2.0; diff --git a/src/shared/particle_generator/particle_generator_lattice.cpp b/src/shared/particle_generator/particle_generator_lattice.cpp index 0683436458..2f179b0fec 100644 --- a/src/shared/particle_generator/particle_generator_lattice.cpp +++ b/src/shared/particle_generator/particle_generator_lattice.cpp @@ -40,7 +40,7 @@ namespace SPH { Real local_particle_spacing = particle_adaptation_->getLocalSpacing(target_shape_, position); Real local_particle_volume_ratio = pow(lattice_spacing_ / local_particle_spacing, Dimensions); - if ((Real)rand() / (RAND_MAX) < local_particle_volume_ratio) + if (Real(rand() / RAND_MAX) < local_particle_volume_ratio) { ParticleGeneratorLattice::initializePositionAndVolumetricMeasure(position, volume / local_particle_volume_ratio); initializeSmoothingLengthRatio(local_particle_spacing); diff --git a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp index 55f745dcb2..fdfd86b0f1 100644 --- a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp +++ b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp @@ -87,8 +87,8 @@ TEST(TransformationMatrix, getTransformationMatrix) { // Generate 3D unit vectors with a spherically symmetric probability distribution // according to the link: https://mathworld.wolfram.com/SpherePointPicking.html - Real u = (Real)rand() / (RAND_MAX); - Real v = (Real)rand() / (RAND_MAX); + Real u = Real(rand() / RAND_MAX); + Real v = Real(rand() / RAND_MAX); Real theta = 2.0 * Pi * u; Real cos_phi = 2.0 * v - 1.0; Real sin_phi = sqrt(1.0 - cos_phi * cos_phi); diff --git a/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp b/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp index 48f97937e6..5f3b428db0 100644 --- a/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp +++ b/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp @@ -215,7 +215,7 @@ int main(int ac, char* av[]) for (int i = 0; i < 10; i++) { - rondom_index.push_back((double)rand() / (RAND_MAX) * shell_particles->total_real_particles_); + rondom_index.push_back(double(rand() / RAND_MAX) * shell_particles->total_real_particles_); von_mises_strain.push_back(shell_particles->getVonMisesStrain(rondom_index[i])); } pseudo_normal = shell_particles->pseudo_n_; From 586d944370e5326aa504ccb611f50759a0213e27 Mon Sep 17 00:00:00 2001 From: Xiangyu Hu Date: Wed, 7 Jun 2023 19:05:43 +0200 Subject: [PATCH 7/8] done --- .../diffusion_dynamics.hpp | 82 +++++++++---------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.hpp b/src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.hpp index 1dba44fb40..e6dcedd093 100644 --- a/src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.hpp +++ b/src/shared/particle_dynamics/diffusion_reaction_dynamics/diffusion_dynamics.hpp @@ -194,7 +194,7 @@ void DiffusionRelaxationDirichlet:: { Real diff_coff_ij = this->all_diffusions_[m]->getInterParticleDiffusionCoff(particle_i, particle_j, e_ij); - Real phi_ij = (*this->diffusion_species_[m])[particle_i] - (*gradient_species_k[m])[particle_j]; + Real phi_ij = (*this->gradient_species_[m])[particle_i] - (*gradient_species_k[m])[particle_j]; (*this->diffusion_dt_[m])[particle_i] += diff_coff_ij * phi_ij * surface_area_ij; } } @@ -280,62 +280,62 @@ void DiffusionRelaxationNeumann:: //=================================================================================================// template DiffusionRelaxationRobin:: - DiffusionRelaxationRobin(BaseContactRelation& contact_relation) - : BaseDiffusionRelaxationContact(contact_relation), - n_(this->particles_->n_) + DiffusionRelaxationRobin(BaseContactRelation &contact_relation) + : BaseDiffusionRelaxationContact(contact_relation), + n_(this->particles_->n_) { - contact_convection_.resize(this->contact_particles_.size()); - contact_T_infinity_.resize(this->all_diffusions_.size()); + contact_convection_.resize(this->contact_particles_.size()); + contact_T_infinity_.resize(this->all_diffusions_.size()); - for (size_t m = 0; m < this->all_diffusions_.size(); ++m) - { - for (size_t k = 0; k != this->contact_particles_.size(); ++k) - { - contact_n_.push_back(&(this->contact_particles_[k]->n_)); - contact_convection_[k] = this->contact_particles_[k]->template registerSharedVariable("Convection"); - contact_T_infinity_[m] = this->contact_particles_[k]->template registerGlobalVariable("T_infinity"); - } - } + for (size_t m = 0; m < this->all_diffusions_.size(); ++m) + { + for (size_t k = 0; k != this->contact_particles_.size(); ++k) + { + contact_n_.push_back(&(this->contact_particles_[k]->n_)); + contact_convection_[k] = this->contact_particles_[k]->template registerSharedVariable("Convection"); + contact_T_infinity_[m] = this->contact_particles_[k]->template registerGlobalVariable("T_infinity"); + } + } } //=================================================================================================// template void DiffusionRelaxationRobin:: - getDiffusionChangeRateRobinContact(size_t particle_i, size_t particle_j, - Real surface_area_ij_Robin, StdLargeVec& convection_k, Real& T_infinity_k) + getDiffusionChangeRateRobinContact(size_t particle_i, size_t particle_j, + Real surface_area_ij_Robin, StdLargeVec &convection_k, Real &T_infinity_k) { - for (size_t m = 0; m < this->all_diffusions_.size(); ++m) - { - Real phi_ij = T_infinity_k - (*this->diffusion_species_[m])[particle_i]; - (*this->diffusion_dt_[m])[particle_i] += convection_k[particle_j] * phi_ij * surface_area_ij_Robin; - } + for (size_t m = 0; m < this->all_diffusions_.size(); ++m) + { + Real phi_ij = T_infinity_k - (*this->diffusion_species_[m])[particle_i]; + (*this->diffusion_dt_[m])[particle_i] += convection_k[particle_j] * phi_ij * surface_area_ij_Robin; + } } //=================================================================================================// template void DiffusionRelaxationRobin:: - interaction(size_t index_i, Real dt) + interaction(size_t index_i, Real dt) { - ParticlesType* particles = this->particles_; + ParticlesType *particles = this->particles_; - for (size_t k = 0; k < this->contact_configuration_.size(); ++k) - { - StdLargeVec& n_k = *(contact_n_[k]); + for (size_t k = 0; k < this->contact_configuration_.size(); ++k) + { + StdLargeVec &n_k = *(contact_n_[k]); - StdLargeVec& convection_k = *(contact_convection_[k]); - Real& T_infinity_k = *(contact_T_infinity_[k]); + StdLargeVec &convection_k = *(contact_convection_[k]); + Real &T_infinity_k = *(contact_T_infinity_[k]); - Neighborhood& contact_neighborhood = (*this->contact_configuration_[k])[index_i]; - for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) - { - size_t index_j = contact_neighborhood.j_[n]; - Real dW_ijV_j_ = contact_neighborhood.dW_ijV_j_[n]; - Vecd& e_ij = contact_neighborhood.e_ij_[n]; + Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Real dW_ijV_j_ = contact_neighborhood.dW_ijV_j_[n]; + Vecd &e_ij = contact_neighborhood.e_ij_[n]; - const Vecd& grad_ijV_j = particles->getKernelGradient(index_i, index_j, dW_ijV_j_, e_ij); - Vecd n_ij = n_[index_i] - n_k[index_j]; - Real area_ij_Robin = grad_ijV_j.dot(n_ij); - getDiffusionChangeRateRobinContact(index_i, index_j, area_ij_Robin, convection_k, T_infinity_k); - } - } + const Vecd &grad_ijV_j = particles->getKernelGradient(index_i, index_j, dW_ijV_j_, e_ij); + Vecd n_ij = n_[index_i] - n_k[index_j]; + Real area_ij_Robin = grad_ijV_j.dot(n_ij); + getDiffusionChangeRateRobinContact(index_i, index_j, area_ij_Robin, convection_k, T_infinity_k); + } + } } //=================================================================================================// template From 8ead465d5a50e5d610fdc29468ea1ef7917282fb Mon Sep 17 00:00:00 2001 From: ChiZhangatTUM Date: Thu, 8 Jun 2023 14:39:00 +0200 Subject: [PATCH 8/8] Fix adaption resoluiton error and update windows action checkout version to 3.5.2 --- .github/workflows/ci.yml | 2 +- .../particle_generator_lattice_supplementary.cpp | 2 +- src/for_3D_build/geometries/triangle_mesh_shape.cpp | 2 +- .../particle_generator_lattice_supplementary.cpp | 2 +- .../particle_generator/particle_generator_network.cpp | 2 +- src/shared/geometries/base_geometry.cpp | 2 +- src/shared/geometries/level_set.cpp | 2 +- .../dissipation_dynamics/particle_dynamics_dissipation.hpp | 2 +- .../particle_dynamics/general_dynamics/general_dynamics.cpp | 2 +- .../general_dynamics/general_dynamics_refinement.cpp | 2 +- src/shared/particle_generator/particle_generator_lattice.cpp | 2 +- .../common/test_vector_functions/test_vector_functions.cpp | 4 ++-- .../test_thin_structure_dynamics.cpp | 2 +- 13 files changed, 14 insertions(+), 14 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ba36b22a65..5a0308d99c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -122,7 +122,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 diff --git a/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp b/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp index 15fe13afe5..b78f797382 100644 --- a/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp +++ b/src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp @@ -60,7 +60,7 @@ namespace SPH { { if (body_shape_.checkContain(particle_position)) { - Real random_real = Real(rand() / RAND_MAX); + Real random_real = (Real)rand() / (RAND_MAX); // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_) { diff --git a/src/for_3D_build/geometries/triangle_mesh_shape.cpp b/src/for_3D_build/geometries/triangle_mesh_shape.cpp index d3f8b5cdf1..7b8dd18e76 100644 --- a/src/for_3D_build/geometries/triangle_mesh_shape.cpp +++ b/src/for_3D_build/geometries/triangle_mesh_shape.cpp @@ -36,7 +36,7 @@ namespace SPH { Vec3d jittered = probe_point; // jittering for (int l = 0; l != probe_point.size(); ++l) - jittered[l] = probe_point[l] + (Real(rand() / RAND_MAX) - 0.5) * (SqrtEps + distance_to_pnt * 0.1); + jittered[l] = probe_point[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * (SqrtEps + distance_to_pnt * 0.1); Vec3d from_face_to_jittered = jittered - SimTKToEigen(closest_pnt); Vec3d direction_to_jittered = from_face_to_jittered / (from_face_to_jittered.norm() + TinyReal); cosine_angle = face_normal.dot(direction_to_jittered); diff --git a/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp b/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp index b15dd639b9..757ec5ea68 100644 --- a/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp +++ b/src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp @@ -66,7 +66,7 @@ namespace SPH { if (body_shape_.checkContain(particle_position)) { - Real random_real = Real(rand() / RAND_MAX); + Real random_real = (Real)rand() / (RAND_MAX); // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles. if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_) { diff --git a/src/for_3D_build/particle_generator/particle_generator_network.cpp b/src/for_3D_build/particle_generator/particle_generator_network.cpp index 614773cdb1..e9e1a2fc45 100644 --- a/src/for_3D_build/particle_generator/particle_generator_network.cpp +++ b/src/for_3D_build/particle_generator/particle_generator_network.cpp @@ -205,7 +205,7 @@ namespace SPH for (size_t j = 0; j != branches_to_grow.size(); j++) { size_t grow_id = branches_to_grow[j]; - Real rand_num = Real(rand() / RAND_MAX) - 0.5; + Real rand_num = ((Real)rand() / (RAND_MAX)) - 0.5; Real angle_to_use = angle_ + rand_num * 0.05; for (size_t k = 0; k != 2; k++) { diff --git a/src/shared/geometries/base_geometry.cpp b/src/shared/geometries/base_geometry.cpp index 5e9ee9a60f..16ffb4a991 100644 --- a/src/shared/geometries/base_geometry.cpp +++ b/src/shared/geometries/base_geometry.cpp @@ -37,7 +37,7 @@ namespace SPH { Vecd jittered = probe_point; for (int l = 0; l != probe_point.size(); ++l) - jittered[l] = probe_point[l] + (Real(rand() / RAND_MAX) - 0.5) * 100.0 * Eps; + jittered[l] = probe_point[l] + (((Real)rand() / (RAND_MAX)) - 0.5) * 100.0 * Eps; if (checkContain(jittered) == is_contain) displacement_to_surface = findClosestPoint(jittered) - jittered; } diff --git a/src/shared/geometries/level_set.cpp b/src/shared/geometries/level_set.cpp index 4d57f69d23..0a1f01a8fb 100644 --- a/src/shared/geometries/level_set.cpp +++ b/src/shared/geometries/level_set.cpp @@ -82,7 +82,7 @@ namespace SPH { Vecd jittered = position; // jittering for (int l = 0; l != position.size(); ++l) - jittered[l] += (Real(rand() / RAND_MAX) - 0.5) * 0.5 * data_spacing_; + jittered[l] += (((Real)rand() / (RAND_MAX)) - 0.5) * 0.5 * data_spacing_; probed_value = probeLevelSetGradient(jittered); } return probed_value.normalized(); diff --git a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp index 33338830aa..c102cf3ade 100644 --- a/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp +++ b/src/shared/particle_dynamics/dissipation_dynamics/particle_dynamics_dissipation.hpp @@ -436,7 +436,7 @@ namespace SPH template bool DampingWithRandomChoice::RandomChoice() { - return Real(rand() / RAND_MAX) < random_ratio_ ? true : false; + return ((Real)rand() / (RAND_MAX)) < random_ratio_ ? true : false; } //=================================================================================================// template diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp index 7b5a8f4106..94ec91997a 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics.cpp @@ -22,7 +22,7 @@ namespace SPH Vecd &pos_n_i = pos_[index_i]; for (int k = 0; k < pos_n_i.size(); ++k) { - pos_n_i[k] += dt * (Real(rand() / RAND_MAX) - 0.5) * 2.0 * randomize_scale_; + pos_n_i[k] += dt * (((Real)rand() / (RAND_MAX)) - 0.5) * 2.0 * randomize_scale_; } } //=================================================================================================// diff --git a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp index f26ed68408..2cc8fe6f0a 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp +++ b/src/shared/particle_dynamics/general_dynamics/general_dynamics_refinement.cpp @@ -376,7 +376,7 @@ void ParticleSplitWithPrescribedArea::update(size_t index_i, Real dt) Vecd ParticleSplitWithPrescribedArea::getSplittingPosition(const StdVec &new_indices) { srand(int(new_indices[0])); - Real delta_random = 0 + 2.0 * Pi * Real(rand() / RAND_MAX) * (2.0 * Pi - 0); + Real delta_random = 0 + 2.0 * Pi * rand() / RAND_MAX * (2.0 * Pi - 0); Vecd pos_ = particles_->pos_[new_indices[0]]; Real delta = delta_random + Pi; Real Vol_split = Vol_[new_indices[0]] / 2.0; diff --git a/src/shared/particle_generator/particle_generator_lattice.cpp b/src/shared/particle_generator/particle_generator_lattice.cpp index 2f179b0fec..0683436458 100644 --- a/src/shared/particle_generator/particle_generator_lattice.cpp +++ b/src/shared/particle_generator/particle_generator_lattice.cpp @@ -40,7 +40,7 @@ namespace SPH { Real local_particle_spacing = particle_adaptation_->getLocalSpacing(target_shape_, position); Real local_particle_volume_ratio = pow(lattice_spacing_ / local_particle_spacing, Dimensions); - if (Real(rand() / RAND_MAX) < local_particle_volume_ratio) + if ((Real)rand() / (RAND_MAX) < local_particle_volume_ratio) { ParticleGeneratorLattice::initializePositionAndVolumetricMeasure(position, volume / local_particle_volume_ratio); initializeSmoothingLengthRatio(local_particle_spacing); diff --git a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp index fdfd86b0f1..55f745dcb2 100644 --- a/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp +++ b/tests/unit_tests_src/shared/common/test_vector_functions/test_vector_functions.cpp @@ -87,8 +87,8 @@ TEST(TransformationMatrix, getTransformationMatrix) { // Generate 3D unit vectors with a spherically symmetric probability distribution // according to the link: https://mathworld.wolfram.com/SpherePointPicking.html - Real u = Real(rand() / RAND_MAX); - Real v = Real(rand() / RAND_MAX); + Real u = (Real)rand() / (RAND_MAX); + Real v = (Real)rand() / (RAND_MAX); Real theta = 2.0 * Pi * u; Real cos_phi = 2.0 * v - 1.0; Real sin_phi = sqrt(1.0 - cos_phi * cos_phi); diff --git a/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp b/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp index 5f3b428db0..48f97937e6 100644 --- a/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp +++ b/tests/unit_tests_src/shared/particle_dynamics/solid_dynamics/test_thin_structure_dynamics/test_thin_structure_dynamics.cpp @@ -215,7 +215,7 @@ int main(int ac, char* av[]) for (int i = 0; i < 10; i++) { - rondom_index.push_back(double(rand() / RAND_MAX) * shell_particles->total_real_particles_); + rondom_index.push_back((double)rand() / (RAND_MAX) * shell_particles->total_real_particles_); von_mises_strain.push_back(shell_particles->getVonMisesStrain(rondom_index[i])); } pseudo_normal = shell_particles->pseudo_n_;