From daa0878faef3d11bf2918d997d3c76eabe75991c Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Sun, 8 Sep 2024 10:06:56 +0200 Subject: [PATCH 01/21] Generalized non-hourglass formulation --- src/shared/materials/general_continuum.cpp | 58 +++++ src/shared/materials/general_continuum.h | 26 ++ .../continuum_integration.cpp | 123 +++++++++- .../continuum_integration.h | 45 ++++ .../oscillating_beam_UL.cpp | 1 - .../CMakeLists.txt | 21 ++ .../oscillating_beam_ULGNOG.cpp | 229 ++++++++++++++++++ .../BeamBody_TotalKineticEnergy.xml | 103 ++++++++ ...amBody_TotalKineticEnergy_Run_0_result.xml | 9 + ...amBody_TotalKineticEnergy_Run_1_result.xml | 9 + ...amBody_TotalKineticEnergy_Run_2_result.xml | 9 + ...amBody_TotalKineticEnergy_Run_3_result.xml | 9 + ...amBody_TotalKineticEnergy_Run_4_result.xml | 9 + ...amBody_TotalKineticEnergy_Run_5_result.xml | 9 + ...eamBody_TotalKineticEnergy_dtwdistance.xml | 4 + .../BeamBody_TotalKineticEnergy_runtimes.dat | 3 + .../regression_test_tool.py | 37 +++ .../test_3d_taylor_bar_UL/CMakeLists.txt | 26 ++ .../MyObserver_Position.xml | 63 +++++ .../MyObserver_Position_Run_0_result.xml | 9 + .../MyObserver_Position_Run_1_result.xml | 9 + .../MyObserver_Position_Run_2_result.xml | 9 + .../MyObserver_Position_Run_3_result.xml | 9 + .../MyObserver_Position_Run_4_result.xml | 9 + .../MyObserver_Position_Run_5_result.xml | 9 + .../MyObserver_Position_dtwdistance.xml | 4 + .../MyObserver_Position_runtimes.dat | 3 + .../regression_test_tool.py | 38 +++ .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 181 ++++++++++++++ .../test_3d_taylor_bar_UL/taylor_bar_UL.h | 132 ++++++++++ 30 files changed, 1203 insertions(+), 2 deletions(-) create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat create mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp create mode 100644 tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index 22afa89d4e..5a63c4b52f 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -77,4 +77,62 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor) } return stress_tensor; } +//=================================================================================================// +Matd J2Plasticity::ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor) +{ + Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); + Matd deviatoric_strain_rate = strain_rate - (1.0 / (Real)Dimensions) * strain_rate.trace() * Matd::Identity(); + Matd shear_stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate; + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + Real lambda_dot_ = 0; + Matd g = Matd::Zero(); + if (f > TinyReal) + { + Real deviatoric_stress_times_strain_rate = (shear_stress.cwiseProduct(strain_rate)).sum(); + lambda_dot_ = deviatoric_stress_times_strain_rate / (sqrt(2.0 * stress_tensor_J2) * (1.0 + hardening_modulus_ / (3.0 * G_))); + g = lambda_dot_ * (sqrt(2.0) * G_ * shear_stress / (sqrt(stress_tensor_J2))); + } + return shear_stress_rate_elastic - g; +} +//=================================================================================================// +Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + if (f > TinyReal) + { + + Real r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); + shear_stress = r * shear_stress; + } + return shear_stress; +} +//=================================================================================================// +Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real r = 1.0; + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + if (f > TinyReal) + { + + r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); + } + return r; +} +//=================================================================================================// +Real J2Plasticity::HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor) +{ + Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); + Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + if (f > TinyReal) + { + return 0.5 * f / (G_ + hardening_modulus_ / 3.0); + } + else + { + return 0.0; + } +} } // namespace SPH \ No newline at end of file diff --git a/src/shared/materials/general_continuum.h b/src/shared/materials/general_continuum.h index 115b4eeba2..dcf59e260d 100644 --- a/src/shared/materials/general_continuum.h +++ b/src/shared/materials/general_continuum.h @@ -98,5 +98,31 @@ class PlasticContinuum : public GeneralContinuum virtual GeneralContinuum *ThisObjectPtr() override { return this; }; }; + +class J2Plasticity : public GeneralContinuum +{ + protected: + Real yield_stress_; + Real hardening_modulus_; + const Real sqrt_2_over_3_ = sqrt(2.0 / 3.0); + + public: + explicit J2Plasticity(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio, Real yield_stress, Real hardening_modulus = 0.0) + : GeneralContinuum(rho0, c0, youngs_modulus, poisson_ratio), + yield_stress_(yield_stress), hardening_modulus_(hardening_modulus) + { + material_type_name_ = "J2Plasticity"; + }; + virtual ~J2Plasticity(){}; + + Real YieldStress() { return yield_stress_; }; + Real HardeningModulus() { return hardening_modulus_; }; + + virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress, Real &hardening_factor); + virtual Matd ReturnMappingShearStress(Matd &shear_stress, Real &hardening_factor); + virtual Real ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor); + virtual Real HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor); + virtual J2Plasticity *ThisObjectPtr() override { return this; }; +}; } // namespace SPH #endif // GENERAL_CONTINUUM_H diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index c72958a5e7..65c07d5236 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -11,6 +11,25 @@ ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body) vel_(*particles_->registerSharedVariable("Velocity")), stress_tensor_3D_(*particles_->registerSharedVariable("StressTensor3D")) {} //=================================================================================================// +AcousticTimeStepSize::AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL) + : LocalDynamicsReduce(sph_body), + DataDelegateSimple(sph_body), fluid_(DynamicCast(this, particles_->getBaseMaterial())), + rho_(*particles_->getVariableDataByName("Density")), + p_(*particles_->getVariableDataByName("Pressure")), + vel_(*particles_->getVariableDataByName("Velocity")), + smoothing_length_min_(sph_body.sph_adaptation_->MinimumSmoothingLength()), + acousticCFL_(acousticCFL) {} +//=================================================================================================// +Real AcousticTimeStepSize::reduce(size_t index_i, Real dt) +{ + return fluid_.getSoundSpeed(p_[index_i], rho_[index_i]) + vel_[index_i].norm(); +} +//=================================================================================================// +Real AcousticTimeStepSize::outputResult(Real reduced_value) +{ + return acousticCFL_ * smoothing_length_min_ / (reduced_value + TinyReal); +} +//=================================================================================================// ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation &inner_relation) : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), @@ -125,5 +144,107 @@ void StressDiffusion::interaction(size_t index_i, Real dt) stress_rate_3D_[index_i] = diffusion_stress_rate_; } //====================================================================================// +ShearAccelerationRelaxationHourglassControl :: + ShearAccelerationRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi) + : fluid_dynamics::BaseIntegration(inner_relation), + continuum_(DynamicCast(this, particles_->getBaseMaterial())), + shear_stress_(*particles_->registerSharedVariable("ShearStress")), + shear_stress_rate_(*particles_->registerSharedVariable("ShearStressRate")), + velocity_gradient_(*particles_->registerSharedVariable("VelocityGradient")), + strain_tensor_(*particles_->registerSharedVariable("StrainTensor")), + strain_tensor_rate_(*particles_->registerSharedVariable("StrainTensorRate")), + B_(*particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), + von_mises_stress_(*particles_->registerSharedVariable("VonMisesStress")), + von_mises_strain_(*particles_->registerSharedVariable("VonMisesStrain")), + scale_penalty_force_(*particles_->registerSharedVariable("ScalePenaltyForce", Real(1.0))), + acc_shear_(*particles_->registerSharedVariable("AccelerationByShear")), + acc_hourglass_(*particles_->registerSharedVariable("AccelerationHourglass")), + G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), xi_(xi) +{ + particles_->addVariableToSort("ShearStress"); + particles_->addVariableToSort("ShearStressRate"); + particles_->addVariableToSort("VelocityGradient"); + particles_->addVariableToSort("StrainTensor"); + particles_->addVariableToSort("StrainTensorRate"); + particles_->addVariableToSort("VonMisesStress"); + particles_->addVariableToSort("VonMisesStrain"); + particles_->addVariableToSort("ScalePenaltyForce"); + particles_->addVariableToSort("AccelerationByShear"); + particles_->addVariableToSort("AccelerationHourglass"); +} +//====================================================================================// +void ShearAccelerationRelaxationHourglassControl::initialization(size_t index_i, Real dt) +{ + Matd velocity_gradient = Matd::Zero(); + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + Vecd vel_i = vel_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; + Vecd &e_ij = inner_neighborhood.e_ij_[n]; + Matd velocity_gradient_ij; + velocity_gradient_ij = -(vel_i - vel_[index_j]) * (B_[index_i] * e_ij * dW_ijV_j).transpose(); + velocity_gradient += velocity_gradient_ij; + } + velocity_gradient_[index_i] = velocity_gradient; +} +//====================================================================================// +void ShearAccelerationRelaxationHourglassControl::interaction(size_t index_i, Real dt) +{ + shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); + shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; + Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); + von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); + strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); + strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt; + von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); +} +//====================================================================================// +void ShearAccelerationRelaxationHourglassControl::update(size_t index_i, Real dt) +{ + Real rho_i = rho_[index_i]; + Matd shear_stress_i = shear_stress_[index_i]; + Vecd acceleration = Vecd::Zero(); + Vecd acceleration_hourglass = Vecd::Zero(); + Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; + Vecd &e_ij = inner_neighborhood.e_ij_[n]; + acceleration += ((shear_stress_i + shear_stress_[index_j]) / rho_i) * dW_ijV_j * e_ij; + Vecd v_ij = vel_[index_i] - vel_[index_j]; + Real r_ij = inner_neighborhood.r_ij_[n]; + Vecd v_ij_correction = v_ij - 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_j]) * r_ij * e_ij; + acceleration_hourglass += 0.5 * xi_ * (scale_penalty_force_[index_i] + scale_penalty_force_[index_j]) * G_ * v_ij_correction * dW_ijV_j * dt / (rho_i * r_ij); + } + acc_hourglass_[index_i] += acceleration_hourglass; + acc_shear_[index_i] = acceleration + acc_hourglass_[index_i]; +} +//====================================================================================// +ShearStressRelaxationHourglassControlJ2Plasticity :: + ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) + : ShearAccelerationRelaxationHourglassControl(inner_relation, xi), + J2_plasticity_(DynamicCast(this, particles_->getBaseMaterial())), + hardening_factor_(*particles_->registerSharedVariable("HardeningFactor")) +{ + particles_->addVariableToSort("HardeningFactor"); +} +//====================================================================================// +void ShearStressRelaxationHourglassControlJ2Plasticity::interaction(size_t index_i, Real dt) +{ + shear_stress_rate_[index_i] = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]); + Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate_[index_i] * dt; + Real hardening_factor_increment = J2_plasticity_.HardeningFactorRate(shear_stress_try, hardening_factor_[index_i]); + hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment; + scale_penalty_force_[index_i] = J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); + shear_stress_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_try, hardening_factor_[index_i]); + Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); + von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); + Matd strain_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); + strain_tensor_[index_i] += strain_rate * dt; + von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); +} } // namespace continuum_dynamics -} // namespace SPH +} // namespace SPH \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 126ba22393..6ae13dcc23 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -49,6 +49,22 @@ class ContinuumInitialCondition : public LocalDynamics, public DataDelegateSimpl StdLargeVec &stress_tensor_3D_; }; +class AcousticTimeStepSize : public LocalDynamicsReduce, public DataDelegateSimple +{ + public: + explicit AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL = 0.6); + virtual ~AcousticTimeStepSize(){}; + Real reduce(size_t index_i, Real dt = 0.0); + virtual Real outputResult(Real reduced_value) override; + + protected: + Fluid &fluid_; + StdLargeVec &rho_, &p_; + StdLargeVec &vel_; + Real smoothing_length_min_; + Real acousticCFL_; +}; + template class BaseIntegration1stHalf : public FluidDynamicsType { @@ -203,6 +219,35 @@ class StressDiffusion : public BasePlasticIntegration Real zeta_ = 0.1, fai_; /*diffusion coefficient*/ Real smoothing_length_, sound_speed_; }; + +class ShearAccelerationRelaxationHourglassControl : public fluid_dynamics::BaseIntegration +{ + public: + explicit ShearAccelerationRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi_ = 4.0); + virtual ~ShearAccelerationRelaxationHourglassControl(){}; + void initialization(size_t index_i, Real dt = 0.0); + void interaction(size_t index_i, Real dt = 0.0); + void update(size_t index_i, Real dt = 0.0); + + protected: + GeneralContinuum &continuum_; + StdLargeVec &shear_stress_, &shear_stress_rate_, &velocity_gradient_, &strain_tensor_, &strain_tensor_rate_, &B_; + StdLargeVec &von_mises_stress_, &von_mises_strain_, &scale_penalty_force_; + StdLargeVec &acc_shear_, &acc_hourglass_; + Real G_, xi_; +}; + +class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearAccelerationRelaxationHourglassControl +{ + public: + explicit ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi_ = 0.2); + virtual ~ShearStressRelaxationHourglassControlJ2Plasticity(){}; + void interaction(size_t index_i, Real dt = 0.0); + + protected: + J2Plasticity &J2_plasticity_; + StdLargeVec &hardening_factor_; +}; } // namespace continuum_dynamics } // namespace SPH #endif // CONTINUUM_INTEGRATION_H \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index 56f021df34..9708a04ad5 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -26,7 +26,6 @@ BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus Real poisson = 0.3975; // Poisson ratio -// Real poisson = 0.4; //Poisson ratio Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); Real gravity_g = 0.0; //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt new file mode 100644 index 0000000000..36b6b1e42d --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt @@ -0,0 +1,21 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ + DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_2d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME} + COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp new file mode 100644 index 0000000000..dc5dbc03ea --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -0,0 +1,229 @@ +/* ----------------------------------------------------------------------------* + * SPHinXsys: 2D oscillation beam * + * ----------------------------------------------------------------------------* + * This is the one of the basic test cases for understanding SPH method for * + * solid simulation based on updated Lagrangian method. * + * A generalized hourglass control method is used here. * + * In this case, the constraint of the beam is implemented with * + * internal constrained subregion. * + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu * + * ----------------------------------------------------------------------------*/ +#include "sphinxsys.h" +using namespace SPH; +//------------------------------------------------------------------------------ +// global parameters for the case +//------------------------------------------------------------------------------ +Real PL = 0.2; // beam length +Real PH = 0.02; // for thick plate +Real SL = 0.06; // depth of the insert +Real resolution_ref = PH / 10; +Real BW = resolution_ref * 4; // boundary width, at least three particles +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), + Vec2d(PL + 3.0 * BW, PL / 2.0)); +//---------------------------------------------------------------------- +// Material properties of the fluid. +//---------------------------------------------------------------------- +Real rho0_s = 1.0e3; // reference density +Real Youngs_modulus = 2.0e6; // reference Youngs modulus +Real poisson = 0.3975; // Poisson ratio +Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); +Real gravity_g = 0.0; +//---------------------------------------------------------------------- +// Parameters for initial condition on velocity +//---------------------------------------------------------------------- +Real kl = 1.875; +Real M = sin(kl) + sinh(kl); +Real N = cos(kl) + cosh(kl); +Real Q = 2.0 * (cos(kl) * sinh(kl) - sin(kl) * cosh(kl)); +Real vf = 0.05; +Real R = PL / (0.5 * Pi); +Real U_ref = vf * c0 * (M * (cos(kl) - cosh(kl)) - N * (sin(kl) - sinh(kl))) / Q; +//---------------------------------------------------------------------- +// Geometric shapes used in the system. +//---------------------------------------------------------------------- +// a beam base shape +std::vector beam_base_shape{ + Vecd(-SL - BW, -PH / 2 - BW), Vecd(-SL - BW, PH / 2 + BW), Vecd(0.0, PH / 2 + BW), + Vecd(0.0, -PH / 2 - BW), Vecd(-SL - BW, -PH / 2 - BW)}; +// a beam shape +std::vector beam_shape{ + Vecd(-SL, -PH / 2), Vecd(-SL, PH / 2), Vecd(PL, PH / 2), Vecd(PL, -PH / 2), Vecd(-SL, -PH / 2)}; +// Beam observer location +StdVec observation_location = {Vecd(PL, 0.0)}; +//---------------------------------------------------------------------- +// Define the beam body +//---------------------------------------------------------------------- +class Beam : public MultiPolygonShape +{ + public: + explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(beam_base_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); + } +}; +//---------------------------------------------------------------------- +// application dependent initial condition +//---------------------------------------------------------------------- +class BeamInitialCondition + : public fluid_dynamics::FluidInitialCondition +{ + public: + explicit BeamInitialCondition(RealBody &beam_column) + : fluid_dynamics::FluidInitialCondition(beam_column){}; + + protected: + void update(size_t index_i, Real dt) + { + /** initial velocity profile */ + Real x = pos_[index_i][0] / PL; + if (x > 0.0) + { + vel_[index_i][1] = vf * c0 * + (M * (cos(kl * x) - cosh(kl * x)) - N * (sin(kl * x) - sinh(kl * x))) / Q; + } + }; +}; +//---------------------------------------------------------------------- +// define the beam base which will be constrained. +//---------------------------------------------------------------------- +MultiPolygon createBeamConstrainShape() +{ + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(beam_base_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(beam_shape, ShapeBooleanOps::sub); + return multi_polygon; +}; +//------------------------------------------------------------------------------ +// the main program +//------------------------------------------------------------------------------ +int main(int ac, char *av[]) +{ + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem with global controls. + //---------------------------------------------------------------------- + SPHSystem sph_system(system_domain_bounds, resolution_ref); + sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + RealBody beam_body(sph_system, makeShared("BeamBody")); + beam_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); + beam_body.generateParticles(); + + ObserverBody beam_observer(sph_system, "BeamObserver"); + beam_observer.sph_adaptation_->resetAdaptationRatios(1.15, 2.0); + beam_observer.generateParticles(observation_location); + //---------------------------------------------------------------------- + // Define body relation map. + // The contact map gives the topological connections between the bodies. + // Basically the the range of bodies to build neighbor particle lists. + // Generally, we first define all the inner relations, then the contact relations. + // At last, we define the complex relaxations by combining previous defined + // inner and contact relations. + //---------------------------------------------------------------------- + ContactRelation beam_observer_contact(beam_observer, {&beam_body}); + InnerRelation beam_body_inner(beam_body); + //----------------------------------------------------------------------------- + // this section define all numerical methods will be used in this case + //----------------------------------------------------------------------------- + /** initial condition */ + InteractionWithUpdate correction_matrix(beam_body_inner); + SimpleDynamics beam_initial_velocity(beam_body); + Dynamics1Level + beam_shear_acceleration_relaxation(beam_body_inner); + Dynamics1Level beam_pressure_relaxation(beam_body_inner); + Dynamics1Level beam_density_relaxation(beam_body_inner); + SimpleDynamics beam_volume_update(beam_body); + ReduceDynamics fluid_advection_time_step(beam_body, U_ref, 0.2); + ReduceDynamics fluid_acoustic_time_step(beam_body, 0.4); + // clamping a solid body part. + BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); + SimpleDynamics constraint_beam_base(beam_base); + //----------------------------------------------------------------------------- + // outputs + //----------------------------------------------------------------------------- + BodyStatesRecordingToVtp write_beam_states(beam_body); + write_beam_states.addToWrite(beam_body, "VonMisesStress"); + write_beam_states.addToWrite(beam_body, "VonMisesStrain"); + write_beam_states.addToWrite(beam_body, "Density"); + write_beam_states.addToWrite(beam_body, "Pressure"); + ObservedQuantityRecording write_beam_tip_displacement("Position", beam_observer_contact); + RegressionTestDynamicTimeWarping> write_beam_kinetic_energy(beam_body); + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- + sph_system.initializeSystemCellLinkedLists(); + sph_system.initializeSystemConfigurations(); + beam_initial_velocity.exec(); + correction_matrix.exec(); + //---------------------------------------------------------------------- + // Setup computing time-step controls. + //---------------------------------------------------------------------- + int ite = 0; + Real T0 = 1.0; + Real End_Time = T0; + Real D_Time = End_Time / 100; /**< Time period for data observing */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + //----------------------------------------------------------------------------- + // from here the time stepping begins + //----------------------------------------------------------------------------- + write_beam_states.writeToFile(0); + write_beam_tip_displacement.writeToFile(0); + write_beam_kinetic_energy.writeToFile(0); + // computation loop starts + while (GlobalStaticVariables::physical_time_ < End_Time) + { + Real integration_time = 0.0; + while (integration_time < D_Time) + { + Real relaxation_time = 0.0; + Real advection_dt = fluid_advection_time_step.exec(); + beam_volume_update.exec(); + while (relaxation_time < advection_dt) + { + Real acoustic_dt = fluid_acoustic_time_step.exec(); + beam_pressure_relaxation.exec(acoustic_dt); + constraint_beam_base.exec(); + beam_shear_acceleration_relaxation.exec(acoustic_dt); + beam_density_relaxation.exec(acoustic_dt); + ite++; + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + if (ite % 500 == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " advection_dt: " + << advection_dt << " acoustic_dt: " + << acoustic_dt << "\n"; + } + } + beam_body.updateCellLinkedList(); + beam_body_inner.updateConfiguration(); + correction_matrix.exec(); + } + write_beam_tip_displacement.writeToFile(ite); + write_beam_kinetic_energy.writeToFile(ite); + TickCount t2 = TickCount::now(); + write_beam_states.writeToFile(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (sph_system.GenerateRegressionData()) + { + write_beam_kinetic_energy.generateDataBase(1.0e-3); + } + else + { + write_beam_kinetic_energy.testResult(); + } + return 0; +} \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml new file mode 100644 index 0000000000..c6b0a3834f --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml new file mode 100644 index 0000000000..7a8cacca58 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml new file mode 100644 index 0000000000..4d9716f718 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml new file mode 100644 index 0000000000..8de5e83ed5 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml new file mode 100644 index 0000000000..4f75fb4c37 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml new file mode 100644 index 0000000000..9ada01ea11 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml new file mode 100644 index 0000000000..0313a9973c --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml new file mode 100644 index 0000000000..e59f0fa6a5 --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat new file mode 100644 index 0000000000..de07de18db --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +4 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..4ad6ca7edc --- /dev/null +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py @@ -0,0 +1,37 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_2d_oscillating_beam_ULGNOG +""" + +case_name = "test_2d_oscillating_beam_ULGNOG" +body_name = "BeamBody" +parameter_name = "TotalKineticEnergy" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_case() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if converged == "true": + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt new file mode 100644 index 0000000000..a7924438eb --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt @@ -0,0 +1,26 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + +set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") + +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "solid_dynamics; contact") diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml new file mode 100644 index 0000000000..da81bfd22e --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml new file mode 100644 index 0000000000..60bdb4418f --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml new file mode 100644 index 0000000000..d23ad8ee5c --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml new file mode 100644 index 0000000000..4764515cdb --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml new file mode 100644 index 0000000000..0d57d75632 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml new file mode 100644 index 0000000000..3886ec8fdb --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml new file mode 100644 index 0000000000..2e8750dac5 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml new file mode 100644 index 0000000000..d8b65a8b22 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat new file mode 100644 index 0000000000..de07de18db --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +4 \ No newline at end of file diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..2612da0e6a --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/regression_test_tool.py @@ -0,0 +1,38 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_3d_taylor_bar +""" + +case_name = "test_3d_taylor_bar_UL" +body_name = "MyObserver" +parameter_name = "Position" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_particle_relaxation() + sphinxsys.run_case_with_reload() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true"): + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp new file mode 100644 index 0000000000..4be91ff9c5 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -0,0 +1,181 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#include "taylor_bar_UL.h" +#include "sphinxsys.h" +using namespace SPH; + +int main(int ac, char *av[]) +{ + /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ + SPHSystem system(system_domain_bounds, particle_spacing_ref); + system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + system.setRunParticleRelaxation(false); + system.setReloadParticles(true); + + RealBody column(system, makeShared("Column")); + column.defineBodyLevelSetShape(); + column.defineMaterial(rho0_s, c0, Youngs_modulus, poisson, yield_stress); + (!system.RunParticleRelaxation() && system.ReloadParticles()) + ? column.generateParticles(column.getName()) + : column.generateParticles(); + + SolidBody wall_boundary(system, makeShared("Wall")); + wall_boundary.defineMaterial(rho0_s, Youngs_modulus, poisson); + wall_boundary.generateParticles(); + + ObserverBody my_observer(system, "MyObserver"); + StdVec observation_location = {Vecd(0.0, 0.0, PW)}; + my_observer.generateParticles(observation_location); + + /**body relation topology */ + InnerRelation column_inner(column); + ContactRelation my_observer_contact(my_observer, {&column}); + SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (system.RunParticleRelaxation()) + { + using namespace relax_dynamics; + /** Random reset the insert body particle position. */ + SimpleDynamics random_column_particles(column); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_column_to_vtp(column); + BodyStatesRecordingToVtp write_states(system); + /** Write the particle reload files. */ + ReloadParticleIO write_particle_reload_files(column); + /** A Physics relaxation step. */ + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); + /** + * @brief Particle relaxation starts here. + */ + random_column_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_states.writeToFile(0.0); + + /** relax particles of the insert body. */ + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 100 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; + write_column_to_vtp.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of cylinder body finish !" << std::endl; + /** Output results. */ + write_particle_reload_files.writeToFile(0.0); + return 0; + } + //---------------------------------------------------------------------- + // All numerical methods will be used in this case. + //---------------------------------------------------------------------- + SimpleDynamics initial_condition(column); + SimpleDynamics wall_normal_direction(wall_boundary); + InteractionWithUpdate corrected_configuration(column_inner); + Dynamics1Level column_shear_stress_relaxation(column_inner); + Dynamics1Level column_pressure_relaxation(column_inner); + Dynamics1Level column_density_relaxation(column_inner); + SimpleDynamics beam_volume_update(column); + ReduceDynamics advection_time_step(column, U_max, 0.2); + ReduceDynamics acoustic_time_step(column, 0.4); + InteractionDynamics column_wall_contact_force(column_wall_contact); + /**define simple data file input and outputs functions. */ + BodyStatesRecordingToVtp write_states(system); + write_states.addToWrite(column, "VonMisesStress"); + write_states.addToWrite(column, "VonMisesStrain"); + write_states.addToWrite(column, "Pressure"); + write_states.addToWrite(column, "Density"); + //---------------------------------------------------------------------- + // Output + //---------------------------------------------------------------------- + RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); + //---------------------------------------------------------------------- + // From here the time stepping begins. + //---------------------------------------------------------------------- + GlobalStaticVariables::physical_time_ = 0.0; + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + wall_normal_direction.exec(); + corrected_configuration.exec(); + initial_condition.exec(); + //---------------------------------------------------------------------- + // Setup time-stepping related simulation parameters. + //---------------------------------------------------------------------- + int ite = 0; + Real end_time = 6.0e-5; + int screen_output_interval = 100; + Real output_period = end_time / 60; + /** Statistics for computing time. */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + write_states.writeToFile(); + write_displacement.writeToFile(0); + //---------------------------------------------------------------------- + // Main time-stepping loop. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + while (integration_time < output_period) + { + Real relaxation_time = 0.0; + Real advection_dt = advection_time_step.exec(); + beam_volume_update.exec(); + while (relaxation_time < advection_dt) + { + Real acoustic_dt = acoustic_time_step.exec(); + if (ite % screen_output_interval == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " advection_dt: " + << advection_dt << " acoustic_dt: " + << acoustic_dt << "\n"; + } + column_wall_contact_force.exec(acoustic_dt); + + column_pressure_relaxation.exec(acoustic_dt); + column_shear_stress_relaxation.exec(acoustic_dt); + column_density_relaxation.exec(acoustic_dt); + + ite++; + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + } + column.updateCellLinkedList(); + column_inner.updateConfiguration(); + column_wall_contact.updateConfiguration(); + corrected_configuration.exec(); + } + TickCount t2 = TickCount::now(); + write_states.writeToFile(ite); + write_displacement.writeToFile(ite); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall_boundary time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (system.GenerateRegressionData()) + { + write_displacement.generateDataBase(0.1); + } + else + { + write_displacement.testResult(); + } + + return 0; +} diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h new file mode 100644 index 0000000000..64d776d673 --- /dev/null +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h @@ -0,0 +1,132 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#pragma once +#include "sphinxsys.h" +using namespace SPH; +//---------------------------------------------------------------------- +// Global geometry parameters. +//---------------------------------------------------------------------- +Real PL = 0.00391; /**< X-direction domain. */ +Real PW = 0.02346; /**< Z-direction domain. */ +Real particle_spacing_ref = PL / 12.0; +Real SL = particle_spacing_ref * 4.0; +Real inner_circle_radius = PL; +Vec3d domain_lower_bound(-4.0 * PL, -4.0 * PL, -SL); +Vec3d domain_upper_bound(4.0 * PL, 4.0 * PL, 2.0 * PW); +BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); +int resolution(20); +//---------------------------------------------------------------------- +// Material properties and global parameters +//---------------------------------------------------------------------- +Real rho0_s = 2700.0; /**< Reference density. */ +Real poisson = 0.3; /**< Poisson ratio. */ +Real Youngs_modulus = 78.2e9; +Real yield_stress = 0.29e9; +Real vel_0 = 373.0; +Real U_max = vel_0; +Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); + +class WallBoundary : public ComplexShape +{ + public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + Vecd halfsize_holder(3.0 * PL, 3.0 * PL, 0.5 * SL); + Vecd translation_holder(0.0, 0.0, -0.5 * SL); + add(halfsize_holder, resolution, translation_holder); + } +}; +/** Define the body. */ +class Column : public ComplexShape +{ +public: + explicit Column(const std::string& shape_name) : ComplexShape(shape_name) + { + Vecd translation_column(0.0, 0.0, 0.5 * PW + particle_spacing_ref); + add(SimTK::UnitVec3(0, 0, 1.0), inner_circle_radius, + 0.5 * PW, resolution, translation_column); + } +}; +/** + * application dependent initial condition + */ +class InitialCondition + : public fluid_dynamics::FluidInitialCondition +{ + public: + explicit InitialCondition(SPHBody &sph_body) + : fluid_dynamics::FluidInitialCondition(sph_body){}; + + void update(size_t index_i, Real dt) + { + vel_[index_i][2] = -vel_0; + } +}; + +class DynamicContactForceWithWall : public LocalDynamics, + public DataDelegateContact +{ + public: + explicit DynamicContactForceWithWall(SurfaceContactRelation &solid_body_contact_relation, Real penalty_strength = 1.0) + : LocalDynamics(solid_body_contact_relation.getSPHBody()), + DataDelegateContact(solid_body_contact_relation), + continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), + Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), + vel_(*particles_->getVariableDataByName("Velocity")), + force_prior_(*particles_->getVariableDataByName("ForcePrior")), + penalty_strength_(penalty_strength) + { + impedance_ = continuum_.ReferenceDensity() * sqrt(continuum_.ContactStiffness()); + reference_pressure_ = continuum_.ReferenceDensity() * continuum_.ContactStiffness(); + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); + contact_vel_.push_back(contact_particles_[k]->registerSharedVariable("Velocity")); + contact_n_.push_back(contact_particles_[k]->template getVariableDataByName("NormalDirection")); + } + }; + virtual ~DynamicContactForceWithWall(){}; + void interaction(size_t index_i, Real dt = 0.0) + { + Vecd force = Vecd::Zero(); + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Real particle_spacing_j1 = 1.0 / contact_bodies_[k]->sph_adaptation_->ReferenceSpacing(); + Real particle_spacing_ratio2 = + 1.0 / (sph_body_.sph_adaptation_->ReferenceSpacing() * particle_spacing_j1); + particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2; + + StdLargeVec &n_k = *(contact_n_[k]); + StdLargeVec &vel_n_k = *(contact_vel_[k]); + StdLargeVec &Vol_k = *(contact_Vol_[k]); + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Vecd e_ij = contact_neighborhood.e_ij_[n]; + Vecd n_k_j = n_k[index_j]; + Real impedance_p = 0.5 * impedance_ * (vel_[index_i] - vel_n_k[index_j]).dot(-n_k_j); + Real overlap = contact_neighborhood.r_ij_[n] * n_k_j.dot(e_ij); + Real delta = 2.0 * overlap * particle_spacing_j1; + Real beta = delta < 1.0 ? (1.0 - delta) * (1.0 - delta) * particle_spacing_ratio2 : 0.0; + Real penalty_p = penalty_strength_ * beta * fabs(overlap) * reference_pressure_; + // force due to pressure + force -= 2.0 * (impedance_p + penalty_p) * e_ij.dot(n_k_j) * + n_k_j * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; + } + } + force_prior_[index_i] += force * Vol_[index_i]; + }; + + protected: + GeneralContinuum &continuum_; + StdLargeVec &Vol_; + StdLargeVec &vel_, &force_prior_; // note that prior force directly used here + StdVec *> contact_Vol_; + StdVec *> contact_vel_, contact_n_; + Real penalty_strength_; + Real impedance_, reference_pressure_; +}; From 0efef079ee16bfef26e7e11b40ee5691f1f5c222 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Sun, 8 Sep 2024 13:06:36 +0200 Subject: [PATCH 02/21] fix bug in regression test --- .../3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt | 11 ++++------- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 3 +-- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt index a7924438eb..bcbd334aa6 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt +++ b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt @@ -7,8 +7,8 @@ SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ + DESTINATION ${BUILD_INPUT_PATH}) add_executable(${PROJECT_NAME}) aux_source_directory(. DIR_SRCS) @@ -17,10 +17,7 @@ target_link_libraries(${PROJECT_NAME} sphinxsys_3d) set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) - + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") - -set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "solid_dynamics; contact") diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 4be91ff9c5..ce65c29e0c 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -44,7 +44,6 @@ int main(int ac, char *av[]) SimpleDynamics random_column_particles(column); /** Write the body state to Vtp file. */ BodyStatesRecordingToVtp write_column_to_vtp(column); - BodyStatesRecordingToVtp write_states(system); /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ @@ -62,7 +61,7 @@ int main(int ac, char *av[]) { relaxation_step_inner.exec(); ite_p += 1; - if (ite_p % 100 == 0) + if (ite_p % 200 == 0) { std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; write_column_to_vtp.writeToFile(ite_p); From 856aa9dae5f031fdc657241a01562c5e5295d65a Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Sun, 8 Sep 2024 13:19:31 +0200 Subject: [PATCH 03/21] fix bug --- tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index ce65c29e0c..c30479c5c6 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -34,6 +34,8 @@ int main(int ac, char *av[]) InnerRelation column_inner(column); ContactRelation my_observer_contact(my_observer, {&column}); SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); + /**define simple data file input and outputs functions. */ + BodyStatesRecordingToVtp write_states(system); //---------------------------------------------------------------------- // Run particle relaxation for body-fitted distribution if chosen. //---------------------------------------------------------------------- @@ -85,8 +87,6 @@ int main(int ac, char *av[]) ReduceDynamics advection_time_step(column, U_max, 0.2); ReduceDynamics acoustic_time_step(column, 0.4); InteractionDynamics column_wall_contact_force(column_wall_contact); - /**define simple data file input and outputs functions. */ - BodyStatesRecordingToVtp write_states(system); write_states.addToWrite(column, "VonMisesStress"); write_states.addToWrite(column, "VonMisesStrain"); write_states.addToWrite(column, "Pressure"); From fc5c7f77297007fcb7732945ae8fd88d795dd4f0 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Sun, 8 Sep 2024 13:22:27 +0200 Subject: [PATCH 04/21] adjust format --- tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index c30479c5c6..2bf421bd53 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -36,6 +36,10 @@ int main(int ac, char *av[]) SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); /**define simple data file input and outputs functions. */ BodyStatesRecordingToVtp write_states(system); + write_states.addToWrite(column, "VonMisesStress"); + write_states.addToWrite(column, "VonMisesStrain"); + write_states.addToWrite(column, "Pressure"); + write_states.addToWrite(column, "Density"); //---------------------------------------------------------------------- // Run particle relaxation for body-fitted distribution if chosen. //---------------------------------------------------------------------- @@ -87,10 +91,6 @@ int main(int ac, char *av[]) ReduceDynamics advection_time_step(column, U_max, 0.2); ReduceDynamics acoustic_time_step(column, 0.4); InteractionDynamics column_wall_contact_force(column_wall_contact); - write_states.addToWrite(column, "VonMisesStress"); - write_states.addToWrite(column, "VonMisesStrain"); - write_states.addToWrite(column, "Pressure"); - write_states.addToWrite(column, "Density"); //---------------------------------------------------------------------- // Output //---------------------------------------------------------------------- From a32f5883019bf552c4bcfb298ffa62e9ab25323e Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Sun, 8 Sep 2024 16:16:37 +0200 Subject: [PATCH 05/21] fix bug for regression test --- .../3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt | 11 +++++++---- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 10 +++++----- .../3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h | 4 ++-- 3 files changed, 14 insertions(+), 11 deletions(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt index bcbd334aa6..acb72b326d 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt +++ b/tests/3d_examples/test_3d_taylor_bar_UL/CMakeLists.txt @@ -7,8 +7,8 @@ SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ - DESTINATION ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) add_executable(${PROJECT_NAME}) aux_source_directory(. DIR_SRCS) @@ -17,7 +17,10 @@ target_link_libraries(${PROJECT_NAME} sphinxsys_3d) set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") + +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "continuum_dynamics") diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 2bf421bd53..38b1aba9e7 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -16,7 +16,7 @@ int main(int ac, char *av[]) system.setReloadParticles(true); RealBody column(system, makeShared("Column")); - column.defineBodyLevelSetShape(); + column.defineBodyLevelSetShape()->writeLevelSet(system); column.defineMaterial(rho0_s, c0, Youngs_modulus, poisson, yield_stress); (!system.RunParticleRelaxation() && system.ReloadParticles()) ? column.generateParticles(column.getName()) @@ -36,10 +36,6 @@ int main(int ac, char *av[]) SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); /**define simple data file input and outputs functions. */ BodyStatesRecordingToVtp write_states(system); - write_states.addToWrite(column, "VonMisesStress"); - write_states.addToWrite(column, "VonMisesStrain"); - write_states.addToWrite(column, "Pressure"); - write_states.addToWrite(column, "Density"); //---------------------------------------------------------------------- // Run particle relaxation for body-fitted distribution if chosen. //---------------------------------------------------------------------- @@ -94,6 +90,10 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // Output //---------------------------------------------------------------------- + write_states.addToWrite(column, "VonMisesStress"); + write_states.addToWrite(column, "VonMisesStrain"); + write_states.addToWrite(column, "Pressure"); + write_states.addToWrite(column, "Density"); RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); //---------------------------------------------------------------------- // From here the time stepping begins. diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h index 64d776d673..d7a1f200f7 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h @@ -28,7 +28,7 @@ Real yield_stress = 0.29e9; Real vel_0 = 373.0; Real U_max = vel_0; Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); - +/** Define the wall. */ class WallBoundary : public ComplexShape { public: @@ -65,7 +65,7 @@ class InitialCondition vel_[index_i][2] = -vel_0; } }; - +/** Contact force. */ class DynamicContactForceWithWall : public LocalDynamics, public DataDelegateContact { From 49fb0e64e6e87954347e715f3643e7b2b3492ea1 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Mon, 9 Sep 2024 22:25:38 +0200 Subject: [PATCH 06/21] test-1 --- .../test_3d_aaa_taylor_bar/CMakeLists.txt | 26 +++ .../MyObserver_Position.xml | 63 ++++++ .../MyObserver_Position_Run_0_result.xml | 9 + .../MyObserver_Position_Run_1_result.xml | 9 + .../MyObserver_Position_Run_2_result.xml | 9 + .../MyObserver_Position_Run_3_result.xml | 9 + .../MyObserver_Position_Run_4_result.xml | 9 + .../MyObserver_Position_Run_5_result.xml | 9 + .../MyObserver_Position_dtwdistance.xml | 4 + .../MyObserver_Position_runtimes.dat | 3 + .../regression_test_tool.py | 38 ++++ .../test_3d_aaa_taylor_bar/taylor_bar_UL.cpp | 180 ++++++++++++++++++ .../test_3d_aaa_taylor_bar/taylor_bar_UL.h | 132 +++++++++++++ 13 files changed, 500 insertions(+) create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp create mode 100644 tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt b/tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt new file mode 100644 index 0000000000..a7924438eb --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt @@ -0,0 +1,26 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) +execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + +set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") + +set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "solid_dynamics; contact") diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml new file mode 100644 index 0000000000..da81bfd22e --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml new file mode 100644 index 0000000000..60bdb4418f --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml new file mode 100644 index 0000000000..d23ad8ee5c --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml new file mode 100644 index 0000000000..4764515cdb --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml new file mode 100644 index 0000000000..0d57d75632 --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml new file mode 100644 index 0000000000..3886ec8fdb --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml new file mode 100644 index 0000000000..2e8750dac5 --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml new file mode 100644 index 0000000000..d8b65a8b22 --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat new file mode 100644 index 0000000000..de07de18db --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +4 \ No newline at end of file diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..2612da0e6a --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py @@ -0,0 +1,38 @@ +# !/usr/bin/env python3 +import os +import sys + +path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') +sys.path.append(path) +from regression_test_base_tool import SphinxsysRegressionTest + +""" +case name: test_3d_taylor_bar +""" + +case_name = "test_3d_taylor_bar_UL" +body_name = "MyObserver" +parameter_name = "Position" + +number_of_run_times = 0 +converged = 0 +sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) + + +while True: + print("Now start a new run......") + sphinxsys.run_particle_relaxation() + sphinxsys.run_case_with_reload() + number_of_run_times += 1 + converged = sphinxsys.read_dat_file() + print("Please note: This is the", number_of_run_times, "run!") + if number_of_run_times <= 200: + if (converged == "true"): + print("The tested parameters of all variables are converged, and the run will stop here!") + break + elif converged != "true": + print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") + continue + else: + print("It's too many runs but still not converged, please try again!") + break diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp new file mode 100644 index 0000000000..38b1aba9e7 --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp @@ -0,0 +1,180 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#include "taylor_bar_UL.h" +#include "sphinxsys.h" +using namespace SPH; + +int main(int ac, char *av[]) +{ + /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ + SPHSystem system(system_domain_bounds, particle_spacing_ref); + system.handleCommandlineOptions(ac, av)->setIOEnvironment(); + system.setRunParticleRelaxation(false); + system.setReloadParticles(true); + + RealBody column(system, makeShared("Column")); + column.defineBodyLevelSetShape()->writeLevelSet(system); + column.defineMaterial(rho0_s, c0, Youngs_modulus, poisson, yield_stress); + (!system.RunParticleRelaxation() && system.ReloadParticles()) + ? column.generateParticles(column.getName()) + : column.generateParticles(); + + SolidBody wall_boundary(system, makeShared("Wall")); + wall_boundary.defineMaterial(rho0_s, Youngs_modulus, poisson); + wall_boundary.generateParticles(); + + ObserverBody my_observer(system, "MyObserver"); + StdVec observation_location = {Vecd(0.0, 0.0, PW)}; + my_observer.generateParticles(observation_location); + + /**body relation topology */ + InnerRelation column_inner(column); + ContactRelation my_observer_contact(my_observer, {&column}); + SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); + /**define simple data file input and outputs functions. */ + BodyStatesRecordingToVtp write_states(system); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (system.RunParticleRelaxation()) + { + using namespace relax_dynamics; + /** Random reset the insert body particle position. */ + SimpleDynamics random_column_particles(column); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_column_to_vtp(column); + /** Write the particle reload files. */ + ReloadParticleIO write_particle_reload_files(column); + /** A Physics relaxation step. */ + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); + /** + * @brief Particle relaxation starts here. + */ + random_column_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_states.writeToFile(0.0); + + /** relax particles of the insert body. */ + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; + write_column_to_vtp.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of cylinder body finish !" << std::endl; + /** Output results. */ + write_particle_reload_files.writeToFile(0.0); + return 0; + } + //---------------------------------------------------------------------- + // All numerical methods will be used in this case. + //---------------------------------------------------------------------- + SimpleDynamics initial_condition(column); + SimpleDynamics wall_normal_direction(wall_boundary); + InteractionWithUpdate corrected_configuration(column_inner); + Dynamics1Level column_shear_stress_relaxation(column_inner); + Dynamics1Level column_pressure_relaxation(column_inner); + Dynamics1Level column_density_relaxation(column_inner); + SimpleDynamics beam_volume_update(column); + ReduceDynamics advection_time_step(column, U_max, 0.2); + ReduceDynamics acoustic_time_step(column, 0.4); + InteractionDynamics column_wall_contact_force(column_wall_contact); + //---------------------------------------------------------------------- + // Output + //---------------------------------------------------------------------- + write_states.addToWrite(column, "VonMisesStress"); + write_states.addToWrite(column, "VonMisesStrain"); + write_states.addToWrite(column, "Pressure"); + write_states.addToWrite(column, "Density"); + RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); + //---------------------------------------------------------------------- + // From here the time stepping begins. + //---------------------------------------------------------------------- + GlobalStaticVariables::physical_time_ = 0.0; + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + wall_normal_direction.exec(); + corrected_configuration.exec(); + initial_condition.exec(); + //---------------------------------------------------------------------- + // Setup time-stepping related simulation parameters. + //---------------------------------------------------------------------- + int ite = 0; + Real end_time = 6.0e-5; + int screen_output_interval = 100; + Real output_period = end_time / 60; + /** Statistics for computing time. */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + write_states.writeToFile(); + write_displacement.writeToFile(0); + //---------------------------------------------------------------------- + // Main time-stepping loop. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + while (integration_time < output_period) + { + Real relaxation_time = 0.0; + Real advection_dt = advection_time_step.exec(); + beam_volume_update.exec(); + while (relaxation_time < advection_dt) + { + Real acoustic_dt = acoustic_time_step.exec(); + if (ite % screen_output_interval == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " advection_dt: " + << advection_dt << " acoustic_dt: " + << acoustic_dt << "\n"; + } + column_wall_contact_force.exec(acoustic_dt); + + column_pressure_relaxation.exec(acoustic_dt); + column_shear_stress_relaxation.exec(acoustic_dt); + column_density_relaxation.exec(acoustic_dt); + + ite++; + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + } + column.updateCellLinkedList(); + column_inner.updateConfiguration(); + column_wall_contact.updateConfiguration(); + corrected_configuration.exec(); + } + TickCount t2 = TickCount::now(); + write_states.writeToFile(ite); + write_displacement.writeToFile(ite); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall_boundary time for computation: " << tt.seconds() << " seconds." << std::endl; + + if (system.GenerateRegressionData()) + { + write_displacement.generateDataBase(0.1); + } + else + { + write_displacement.testResult(); + } + + return 0; +} diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h b/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h new file mode 100644 index 0000000000..d7a1f200f7 --- /dev/null +++ b/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h @@ -0,0 +1,132 @@ +/** + * @file taylor_bar.cpp + * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. + * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu + */ +#pragma once +#include "sphinxsys.h" +using namespace SPH; +//---------------------------------------------------------------------- +// Global geometry parameters. +//---------------------------------------------------------------------- +Real PL = 0.00391; /**< X-direction domain. */ +Real PW = 0.02346; /**< Z-direction domain. */ +Real particle_spacing_ref = PL / 12.0; +Real SL = particle_spacing_ref * 4.0; +Real inner_circle_radius = PL; +Vec3d domain_lower_bound(-4.0 * PL, -4.0 * PL, -SL); +Vec3d domain_upper_bound(4.0 * PL, 4.0 * PL, 2.0 * PW); +BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); +int resolution(20); +//---------------------------------------------------------------------- +// Material properties and global parameters +//---------------------------------------------------------------------- +Real rho0_s = 2700.0; /**< Reference density. */ +Real poisson = 0.3; /**< Poisson ratio. */ +Real Youngs_modulus = 78.2e9; +Real yield_stress = 0.29e9; +Real vel_0 = 373.0; +Real U_max = vel_0; +Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); +/** Define the wall. */ +class WallBoundary : public ComplexShape +{ + public: + explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) + { + Vecd halfsize_holder(3.0 * PL, 3.0 * PL, 0.5 * SL); + Vecd translation_holder(0.0, 0.0, -0.5 * SL); + add(halfsize_holder, resolution, translation_holder); + } +}; +/** Define the body. */ +class Column : public ComplexShape +{ +public: + explicit Column(const std::string& shape_name) : ComplexShape(shape_name) + { + Vecd translation_column(0.0, 0.0, 0.5 * PW + particle_spacing_ref); + add(SimTK::UnitVec3(0, 0, 1.0), inner_circle_radius, + 0.5 * PW, resolution, translation_column); + } +}; +/** + * application dependent initial condition + */ +class InitialCondition + : public fluid_dynamics::FluidInitialCondition +{ + public: + explicit InitialCondition(SPHBody &sph_body) + : fluid_dynamics::FluidInitialCondition(sph_body){}; + + void update(size_t index_i, Real dt) + { + vel_[index_i][2] = -vel_0; + } +}; +/** Contact force. */ +class DynamicContactForceWithWall : public LocalDynamics, + public DataDelegateContact +{ + public: + explicit DynamicContactForceWithWall(SurfaceContactRelation &solid_body_contact_relation, Real penalty_strength = 1.0) + : LocalDynamics(solid_body_contact_relation.getSPHBody()), + DataDelegateContact(solid_body_contact_relation), + continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), + Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), + vel_(*particles_->getVariableDataByName("Velocity")), + force_prior_(*particles_->getVariableDataByName("ForcePrior")), + penalty_strength_(penalty_strength) + { + impedance_ = continuum_.ReferenceDensity() * sqrt(continuum_.ContactStiffness()); + reference_pressure_ = continuum_.ReferenceDensity() * continuum_.ContactStiffness(); + for (size_t k = 0; k != contact_particles_.size(); ++k) + { + contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); + contact_vel_.push_back(contact_particles_[k]->registerSharedVariable("Velocity")); + contact_n_.push_back(contact_particles_[k]->template getVariableDataByName("NormalDirection")); + } + }; + virtual ~DynamicContactForceWithWall(){}; + void interaction(size_t index_i, Real dt = 0.0) + { + Vecd force = Vecd::Zero(); + for (size_t k = 0; k < contact_configuration_.size(); ++k) + { + Real particle_spacing_j1 = 1.0 / contact_bodies_[k]->sph_adaptation_->ReferenceSpacing(); + Real particle_spacing_ratio2 = + 1.0 / (sph_body_.sph_adaptation_->ReferenceSpacing() * particle_spacing_j1); + particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2; + + StdLargeVec &n_k = *(contact_n_[k]); + StdLargeVec &vel_n_k = *(contact_vel_[k]); + StdLargeVec &Vol_k = *(contact_Vol_[k]); + Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; + for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) + { + size_t index_j = contact_neighborhood.j_[n]; + Vecd e_ij = contact_neighborhood.e_ij_[n]; + Vecd n_k_j = n_k[index_j]; + Real impedance_p = 0.5 * impedance_ * (vel_[index_i] - vel_n_k[index_j]).dot(-n_k_j); + Real overlap = contact_neighborhood.r_ij_[n] * n_k_j.dot(e_ij); + Real delta = 2.0 * overlap * particle_spacing_j1; + Real beta = delta < 1.0 ? (1.0 - delta) * (1.0 - delta) * particle_spacing_ratio2 : 0.0; + Real penalty_p = penalty_strength_ * beta * fabs(overlap) * reference_pressure_; + // force due to pressure + force -= 2.0 * (impedance_p + penalty_p) * e_ij.dot(n_k_j) * + n_k_j * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; + } + } + force_prior_[index_i] += force * Vol_[index_i]; + }; + + protected: + GeneralContinuum &continuum_; + StdLargeVec &Vol_; + StdLargeVec &vel_, &force_prior_; // note that prior force directly used here + StdVec *> contact_Vol_; + StdVec *> contact_vel_, contact_n_; + Real penalty_strength_; + Real impedance_, reference_pressure_; +}; From 9e54b645bfbfbd4da032eb48828bce4817289a6c Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 10 Sep 2024 01:03:05 +0200 Subject: [PATCH 07/21] test-2 --- .../CMakeLists.txt | 0 .../regression_test_tool/MyObserver_Position.xml | 0 .../MyObserver_Position_Run_0_result.xml | 0 .../MyObserver_Position_Run_1_result.xml | 0 .../MyObserver_Position_Run_2_result.xml | 0 .../MyObserver_Position_Run_3_result.xml | 0 .../MyObserver_Position_Run_4_result.xml | 0 .../MyObserver_Position_Run_5_result.xml | 0 .../MyObserver_Position_dtwdistance.xml | 0 .../regression_test_tool/MyObserver_Position_runtimes.dat | 0 .../regression_test_tool/regression_test_tool.py | 0 .../taylor_bar_UL.cpp | 6 ++---- .../taylor_bar_UL.h | 1 - tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 5 +---- 14 files changed, 3 insertions(+), 9 deletions(-) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/CMakeLists.txt (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_0_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_1_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_2_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_3_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_4_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_Run_5_result.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_dtwdistance.xml (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/MyObserver_Position_runtimes.dat (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/regression_test_tool/regression_test_tool.py (100%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/taylor_bar_UL.cpp (98%) rename tests/3d_examples/{test_3d_aaa_taylor_bar => test_0_taylor_bar_UL}/taylor_bar_UL.h (99%) diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt b/tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/CMakeLists.txt rename to tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_0_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_1_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_2_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_3_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_4_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_Run_5_result.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_dtwdistance.xml rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/MyObserver_Position_runtimes.dat rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py similarity index 100% rename from tests/3d_examples/test_3d_aaa_taylor_bar/regression_test_tool/regression_test_tool.py rename to tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp similarity index 98% rename from tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp rename to tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp index 38b1aba9e7..47b1eadc2b 100644 --- a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp @@ -4,16 +4,14 @@ * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu */ #include "taylor_bar_UL.h" -#include "sphinxsys.h" -using namespace SPH; int main(int ac, char *av[]) { /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ SPHSystem system(system_domain_bounds, particle_spacing_ref); - system.handleCommandlineOptions(ac, av)->setIOEnvironment(); system.setRunParticleRelaxation(false); system.setReloadParticles(true); + system.handleCommandlineOptions(ac, av)->setIOEnvironment(); RealBody column(system, makeShared("Column")); column.defineBodyLevelSetShape()->writeLevelSet(system); @@ -49,7 +47,7 @@ int main(int ac, char *av[]) /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ - RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); + RelaxationStepLevelSetInner relaxation_step_inner(column_inner); /** * @brief Particle relaxation starts here. */ diff --git a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h similarity index 99% rename from tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h rename to tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h index d7a1f200f7..19d4130366 100644 --- a/tests/3d_examples/test_3d_aaa_taylor_bar/taylor_bar_UL.h +++ b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h @@ -3,7 +3,6 @@ * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu */ -#pragma once #include "sphinxsys.h" using namespace SPH; //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 38b1aba9e7..bc3d1d6165 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -4,9 +4,6 @@ * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu */ #include "taylor_bar_UL.h" -#include "sphinxsys.h" -using namespace SPH; - int main(int ac, char *av[]) { /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ @@ -49,7 +46,7 @@ int main(int ac, char *av[]) /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ - RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); + RelaxationStepInner relaxation_step_inner(column_inner); /** * @brief Particle relaxation starts here. */ From 7fa96dbf28ca444127a2186e73b869523d097972 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 10 Sep 2024 05:09:43 +0200 Subject: [PATCH 08/21] test-3 --- tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp index 47b1eadc2b..e95a7bc606 100644 --- a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp @@ -47,7 +47,7 @@ int main(int ac, char *av[]) /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ - RelaxationStepLevelSetInner relaxation_step_inner(column_inner); + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); /** * @brief Particle relaxation starts here. */ From f68f50b4e141bdb77d0e0b7fab9b1fbfe9f54f49 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 10 Sep 2024 08:41:48 +0200 Subject: [PATCH 09/21] solve the bug --- .../test_0_taylor_bar_UL/CMakeLists.txt | 26 --- .../MyObserver_Position.xml | 63 ------- .../MyObserver_Position_Run_0_result.xml | 9 - .../MyObserver_Position_Run_1_result.xml | 9 - .../MyObserver_Position_Run_2_result.xml | 9 - .../MyObserver_Position_Run_3_result.xml | 9 - .../MyObserver_Position_Run_4_result.xml | 9 - .../MyObserver_Position_Run_5_result.xml | 9 - .../MyObserver_Position_dtwdistance.xml | 4 - .../MyObserver_Position_runtimes.dat | 3 - .../regression_test_tool.py | 38 ---- .../test_0_taylor_bar_UL/taylor_bar_UL.cpp | 178 ------------------ .../test_0_taylor_bar_UL/taylor_bar_UL.h | 131 ------------- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 4 +- 14 files changed, 2 insertions(+), 499 deletions(-) delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp delete mode 100644 tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h diff --git a/tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt b/tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt deleted file mode 100644 index a7924438eb..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/CMakeLists.txt +++ /dev/null @@ -1,26 +0,0 @@ -STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) -PROJECT("${CURRENT_FOLDER}") - -SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) -SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") -SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") -SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") - -file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -execute_process(COMMAND ${CMAKE_COMMAND} -E make_directory ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ DESTINATION ${BUILD_INPUT_PATH}) - -add_executable(${PROJECT_NAME}) -aux_source_directory(. DIR_SRCS) -target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) -target_link_libraries(${PROJECT_NAME} sphinxsys_3d) -set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") - -add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --relax=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) -add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --reload=true --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) - -set_tests_properties(${PROJECT_NAME} PROPERTIES DEPENDS "${PROJECT_NAME}_particle_relaxation") - -set_tests_properties(${PROJECT_NAME} PROPERTIES LABELS "solid_dynamics; contact") diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml deleted file mode 100644 index da81bfd22e..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position.xml +++ /dev/null @@ -1,63 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml deleted file mode 100644 index 60bdb4418f..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_0_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml deleted file mode 100644 index d23ad8ee5c..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_1_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml deleted file mode 100644 index 4764515cdb..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_2_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml deleted file mode 100644 index 0d57d75632..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_3_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml deleted file mode 100644 index 3886ec8fdb..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_4_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml deleted file mode 100644 index 2e8750dac5..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_Run_5_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml deleted file mode 100644 index d8b65a8b22..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml +++ /dev/null @@ -1,4 +0,0 @@ - - - - diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat deleted file mode 100644 index de07de18db..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/MyObserver_Position_runtimes.dat +++ /dev/null @@ -1,3 +0,0 @@ -true -6 -4 \ No newline at end of file diff --git a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py deleted file mode 100644 index 2612da0e6a..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/regression_test_tool/regression_test_tool.py +++ /dev/null @@ -1,38 +0,0 @@ -# !/usr/bin/env python3 -import os -import sys - -path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') -sys.path.append(path) -from regression_test_base_tool import SphinxsysRegressionTest - -""" -case name: test_3d_taylor_bar -""" - -case_name = "test_3d_taylor_bar_UL" -body_name = "MyObserver" -parameter_name = "Position" - -number_of_run_times = 0 -converged = 0 -sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) - - -while True: - print("Now start a new run......") - sphinxsys.run_particle_relaxation() - sphinxsys.run_case_with_reload() - number_of_run_times += 1 - converged = sphinxsys.read_dat_file() - print("Please note: This is the", number_of_run_times, "run!") - if number_of_run_times <= 200: - if (converged == "true"): - print("The tested parameters of all variables are converged, and the run will stop here!") - break - elif converged != "true": - print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") - continue - else: - print("It's too many runs but still not converged, please try again!") - break diff --git a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp deleted file mode 100644 index e95a7bc606..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.cpp +++ /dev/null @@ -1,178 +0,0 @@ -/** - * @file taylor_bar.cpp - * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. - * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu - */ -#include "taylor_bar_UL.h" - -int main(int ac, char *av[]) -{ - /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ - SPHSystem system(system_domain_bounds, particle_spacing_ref); - system.setRunParticleRelaxation(false); - system.setReloadParticles(true); - system.handleCommandlineOptions(ac, av)->setIOEnvironment(); - - RealBody column(system, makeShared("Column")); - column.defineBodyLevelSetShape()->writeLevelSet(system); - column.defineMaterial(rho0_s, c0, Youngs_modulus, poisson, yield_stress); - (!system.RunParticleRelaxation() && system.ReloadParticles()) - ? column.generateParticles(column.getName()) - : column.generateParticles(); - - SolidBody wall_boundary(system, makeShared("Wall")); - wall_boundary.defineMaterial(rho0_s, Youngs_modulus, poisson); - wall_boundary.generateParticles(); - - ObserverBody my_observer(system, "MyObserver"); - StdVec observation_location = {Vecd(0.0, 0.0, PW)}; - my_observer.generateParticles(observation_location); - - /**body relation topology */ - InnerRelation column_inner(column); - ContactRelation my_observer_contact(my_observer, {&column}); - SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); - /**define simple data file input and outputs functions. */ - BodyStatesRecordingToVtp write_states(system); - //---------------------------------------------------------------------- - // Run particle relaxation for body-fitted distribution if chosen. - //---------------------------------------------------------------------- - if (system.RunParticleRelaxation()) - { - using namespace relax_dynamics; - /** Random reset the insert body particle position. */ - SimpleDynamics random_column_particles(column); - /** Write the body state to Vtp file. */ - BodyStatesRecordingToVtp write_column_to_vtp(column); - /** Write the particle reload files. */ - ReloadParticleIO write_particle_reload_files(column); - /** A Physics relaxation step. */ - RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); - /** - * @brief Particle relaxation starts here. - */ - random_column_particles.exec(0.25); - relaxation_step_inner.SurfaceBounding().exec(); - write_states.writeToFile(0.0); - - /** relax particles of the insert body. */ - int ite_p = 0; - while (ite_p < 1000) - { - relaxation_step_inner.exec(); - ite_p += 1; - if (ite_p % 200 == 0) - { - std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the column body N = " << ite_p << "\n"; - write_column_to_vtp.writeToFile(ite_p); - } - } - std::cout << "The physics relaxation process of cylinder body finish !" << std::endl; - /** Output results. */ - write_particle_reload_files.writeToFile(0.0); - return 0; - } - //---------------------------------------------------------------------- - // All numerical methods will be used in this case. - //---------------------------------------------------------------------- - SimpleDynamics initial_condition(column); - SimpleDynamics wall_normal_direction(wall_boundary); - InteractionWithUpdate corrected_configuration(column_inner); - Dynamics1Level column_shear_stress_relaxation(column_inner); - Dynamics1Level column_pressure_relaxation(column_inner); - Dynamics1Level column_density_relaxation(column_inner); - SimpleDynamics beam_volume_update(column); - ReduceDynamics advection_time_step(column, U_max, 0.2); - ReduceDynamics acoustic_time_step(column, 0.4); - InteractionDynamics column_wall_contact_force(column_wall_contact); - //---------------------------------------------------------------------- - // Output - //---------------------------------------------------------------------- - write_states.addToWrite(column, "VonMisesStress"); - write_states.addToWrite(column, "VonMisesStrain"); - write_states.addToWrite(column, "Pressure"); - write_states.addToWrite(column, "Density"); - RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); - //---------------------------------------------------------------------- - // From here the time stepping begins. - //---------------------------------------------------------------------- - GlobalStaticVariables::physical_time_ = 0.0; - system.initializeSystemCellLinkedLists(); - system.initializeSystemConfigurations(); - wall_normal_direction.exec(); - corrected_configuration.exec(); - initial_condition.exec(); - //---------------------------------------------------------------------- - // Setup time-stepping related simulation parameters. - //---------------------------------------------------------------------- - int ite = 0; - Real end_time = 6.0e-5; - int screen_output_interval = 100; - Real output_period = end_time / 60; - /** Statistics for computing time. */ - TickCount t1 = TickCount::now(); - TimeInterval interval; - //---------------------------------------------------------------------- - // First output before the main loop. - //---------------------------------------------------------------------- - write_states.writeToFile(); - write_displacement.writeToFile(0); - //---------------------------------------------------------------------- - // Main time-stepping loop. - //---------------------------------------------------------------------- - while (GlobalStaticVariables::physical_time_ < end_time) - { - Real integration_time = 0.0; - while (integration_time < output_period) - { - Real relaxation_time = 0.0; - Real advection_dt = advection_time_step.exec(); - beam_volume_update.exec(); - while (relaxation_time < advection_dt) - { - Real acoustic_dt = acoustic_time_step.exec(); - if (ite % screen_output_interval == 0) - { - std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ << " advection_dt: " - << advection_dt << " acoustic_dt: " - << acoustic_dt << "\n"; - } - column_wall_contact_force.exec(acoustic_dt); - - column_pressure_relaxation.exec(acoustic_dt); - column_shear_stress_relaxation.exec(acoustic_dt); - column_density_relaxation.exec(acoustic_dt); - - ite++; - relaxation_time += acoustic_dt; - integration_time += acoustic_dt; - GlobalStaticVariables::physical_time_ += acoustic_dt; - } - column.updateCellLinkedList(); - column_inner.updateConfiguration(); - column_wall_contact.updateConfiguration(); - corrected_configuration.exec(); - } - TickCount t2 = TickCount::now(); - write_states.writeToFile(ite); - write_displacement.writeToFile(ite); - TickCount t3 = TickCount::now(); - interval += t3 - t2; - } - TickCount t4 = TickCount::now(); - TimeInterval tt; - tt = t4 - t1 - interval; - std::cout << "Total wall_boundary time for computation: " << tt.seconds() << " seconds." << std::endl; - - if (system.GenerateRegressionData()) - { - write_displacement.generateDataBase(0.1); - } - else - { - write_displacement.testResult(); - } - - return 0; -} diff --git a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h b/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h deleted file mode 100644 index 19d4130366..0000000000 --- a/tests/3d_examples/test_0_taylor_bar_UL/taylor_bar_UL.h +++ /dev/null @@ -1,131 +0,0 @@ -/** - * @file taylor_bar.cpp - * @brief This is the case setup for plastic taylor bar using updated Lagragian SPH. - * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu - */ -#include "sphinxsys.h" -using namespace SPH; -//---------------------------------------------------------------------- -// Global geometry parameters. -//---------------------------------------------------------------------- -Real PL = 0.00391; /**< X-direction domain. */ -Real PW = 0.02346; /**< Z-direction domain. */ -Real particle_spacing_ref = PL / 12.0; -Real SL = particle_spacing_ref * 4.0; -Real inner_circle_radius = PL; -Vec3d domain_lower_bound(-4.0 * PL, -4.0 * PL, -SL); -Vec3d domain_upper_bound(4.0 * PL, 4.0 * PL, 2.0 * PW); -BoundingBox system_domain_bounds(domain_lower_bound, domain_upper_bound); -int resolution(20); -//---------------------------------------------------------------------- -// Material properties and global parameters -//---------------------------------------------------------------------- -Real rho0_s = 2700.0; /**< Reference density. */ -Real poisson = 0.3; /**< Poisson ratio. */ -Real Youngs_modulus = 78.2e9; -Real yield_stress = 0.29e9; -Real vel_0 = 373.0; -Real U_max = vel_0; -Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -/** Define the wall. */ -class WallBoundary : public ComplexShape -{ - public: - explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name) - { - Vecd halfsize_holder(3.0 * PL, 3.0 * PL, 0.5 * SL); - Vecd translation_holder(0.0, 0.0, -0.5 * SL); - add(halfsize_holder, resolution, translation_holder); - } -}; -/** Define the body. */ -class Column : public ComplexShape -{ -public: - explicit Column(const std::string& shape_name) : ComplexShape(shape_name) - { - Vecd translation_column(0.0, 0.0, 0.5 * PW + particle_spacing_ref); - add(SimTK::UnitVec3(0, 0, 1.0), inner_circle_radius, - 0.5 * PW, resolution, translation_column); - } -}; -/** - * application dependent initial condition - */ -class InitialCondition - : public fluid_dynamics::FluidInitialCondition -{ - public: - explicit InitialCondition(SPHBody &sph_body) - : fluid_dynamics::FluidInitialCondition(sph_body){}; - - void update(size_t index_i, Real dt) - { - vel_[index_i][2] = -vel_0; - } -}; -/** Contact force. */ -class DynamicContactForceWithWall : public LocalDynamics, - public DataDelegateContact -{ - public: - explicit DynamicContactForceWithWall(SurfaceContactRelation &solid_body_contact_relation, Real penalty_strength = 1.0) - : LocalDynamics(solid_body_contact_relation.getSPHBody()), - DataDelegateContact(solid_body_contact_relation), - continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), - Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), - vel_(*particles_->getVariableDataByName("Velocity")), - force_prior_(*particles_->getVariableDataByName("ForcePrior")), - penalty_strength_(penalty_strength) - { - impedance_ = continuum_.ReferenceDensity() * sqrt(continuum_.ContactStiffness()); - reference_pressure_ = continuum_.ReferenceDensity() * continuum_.ContactStiffness(); - for (size_t k = 0; k != contact_particles_.size(); ++k) - { - contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); - contact_vel_.push_back(contact_particles_[k]->registerSharedVariable("Velocity")); - contact_n_.push_back(contact_particles_[k]->template getVariableDataByName("NormalDirection")); - } - }; - virtual ~DynamicContactForceWithWall(){}; - void interaction(size_t index_i, Real dt = 0.0) - { - Vecd force = Vecd::Zero(); - for (size_t k = 0; k < contact_configuration_.size(); ++k) - { - Real particle_spacing_j1 = 1.0 / contact_bodies_[k]->sph_adaptation_->ReferenceSpacing(); - Real particle_spacing_ratio2 = - 1.0 / (sph_body_.sph_adaptation_->ReferenceSpacing() * particle_spacing_j1); - particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2; - - StdLargeVec &n_k = *(contact_n_[k]); - StdLargeVec &vel_n_k = *(contact_vel_[k]); - StdLargeVec &Vol_k = *(contact_Vol_[k]); - Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; - for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) - { - size_t index_j = contact_neighborhood.j_[n]; - Vecd e_ij = contact_neighborhood.e_ij_[n]; - Vecd n_k_j = n_k[index_j]; - Real impedance_p = 0.5 * impedance_ * (vel_[index_i] - vel_n_k[index_j]).dot(-n_k_j); - Real overlap = contact_neighborhood.r_ij_[n] * n_k_j.dot(e_ij); - Real delta = 2.0 * overlap * particle_spacing_j1; - Real beta = delta < 1.0 ? (1.0 - delta) * (1.0 - delta) * particle_spacing_ratio2 : 0.0; - Real penalty_p = penalty_strength_ * beta * fabs(overlap) * reference_pressure_; - // force due to pressure - force -= 2.0 * (impedance_p + penalty_p) * e_ij.dot(n_k_j) * - n_k_j * contact_neighborhood.dW_ij_[n] * Vol_k[index_j]; - } - } - force_prior_[index_i] += force * Vol_[index_i]; - }; - - protected: - GeneralContinuum &continuum_; - StdLargeVec &Vol_; - StdLargeVec &vel_, &force_prior_; // note that prior force directly used here - StdVec *> contact_Vol_; - StdVec *> contact_vel_, contact_n_; - Real penalty_strength_; - Real impedance_, reference_pressure_; -}; diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index bc3d1d6165..7c4f61dc0e 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -8,9 +8,9 @@ int main(int ac, char *av[]) { /** Setup the system. Please the make sure the global domain bounds are correctly defined. */ SPHSystem system(system_domain_bounds, particle_spacing_ref); - system.handleCommandlineOptions(ac, av)->setIOEnvironment(); system.setRunParticleRelaxation(false); system.setReloadParticles(true); + system.handleCommandlineOptions(ac, av)->setIOEnvironment(); RealBody column(system, makeShared("Column")); column.defineBodyLevelSetShape()->writeLevelSet(system); @@ -46,7 +46,7 @@ int main(int ac, char *av[]) /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ - RelaxationStepInner relaxation_step_inner(column_inner); + RelaxationStepLevelSetCorrectionInner relaxation_step_inner(column_inner); /** * @brief Particle relaxation starts here. */ From d18fb68299c6096c15072145ff53780aba41b537 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 10 Sep 2024 08:55:19 +0200 Subject: [PATCH 10/21] small changes --- src/shared/materials/general_continuum.cpp | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index 5a63c4b52f..d728d54b81 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -102,7 +102,6 @@ Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_ Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); if (f > TinyReal) { - Real r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); shear_stress = r * shear_stress; } @@ -116,7 +115,6 @@ Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); if (f > TinyReal) { - r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); } return r; @@ -126,13 +124,6 @@ Real J2Plasticity::HardeningFactorRate(const Matd &shear_stress, Real &hardening { Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); - if (f > TinyReal) - { - return 0.5 * f / (G_ + hardening_modulus_ / 3.0); - } - else - { - return 0.0; - } + return (f > TinyReal) ? 0.5 * f / (G_ + hardening_modulus_ / 3.0) : 0.0; } } // namespace SPH \ No newline at end of file From c2c4dc7eb4689476b2d20d18141db284c41ba8a3 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Mon, 16 Sep 2024 18:58:26 +0200 Subject: [PATCH 11/21] naming changes --- src/shared/materials/general_continuum.cpp | 2 -- src/shared/materials/general_continuum.h | 2 +- .../continuum_integration.cpp | 18 +++++++++--------- .../continuum_dynamics/continuum_integration.h | 10 +++++----- .../oscillating_beam_ULGNOG.cpp | 2 +- 5 files changed, 16 insertions(+), 18 deletions(-) diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index d728d54b81..92c94430d8 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -114,9 +114,7 @@ Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) Real r = 1.0; Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); if (f > TinyReal) - { r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); - } return r; } //=================================================================================================// diff --git a/src/shared/materials/general_continuum.h b/src/shared/materials/general_continuum.h index dcf59e260d..e809ef7f5a 100644 --- a/src/shared/materials/general_continuum.h +++ b/src/shared/materials/general_continuum.h @@ -43,7 +43,7 @@ class GeneralContinuum : public WeaklyCompressibleFluid Real contact_stiffness_; /* contact-force stiffness related to bulk modulus*/ public: explicit GeneralContinuum(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio) - : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(c0 * c0) + : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(rho0_ * c0 * c0) { material_type_name_ = "GeneralContinuum"; E_ = youngs_modulus; diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 65c07d5236..f9b0888e01 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -116,7 +116,7 @@ void ShearStressRelaxation::update(size_t index_i, Real dt) //====================================================================================// StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation) : BasePlasticIntegration(inner_relation), - fai_(DynamicCast(this, plastic_continuum_).getFrictionAngle()), + phi_(DynamicCast(this, plastic_continuum_).getFrictionAngle()), smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), sound_speed_(plastic_continuum_.ReferenceSoundSpeed()) {} //====================================================================================// @@ -135,17 +135,17 @@ void StressDiffusion::interaction(size_t index_i, Real dt) Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; - diffusion_stress_(0, 0) -= (1 - sin(fai_)) * density * gravity * y_ij; + diffusion_stress_(0, 0) -= (1 - sin(phi_)) * density * gravity * y_ij; diffusion_stress_(1, 1) -= density * gravity * y_ij; - diffusion_stress_(2, 2) -= (1 - sin(fai_)) * density * gravity * y_ij; + diffusion_stress_(2, 2) -= (1 - sin(phi_)) * density * gravity * y_ij; diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); } stress_rate_3D_[index_i] = diffusion_stress_rate_; } //====================================================================================// -ShearAccelerationRelaxationHourglassControl :: - ShearAccelerationRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi) +ShearStressRelaxationHourglassControl :: + ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi) : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), shear_stress_(*particles_->registerSharedVariable("ShearStress")), @@ -173,7 +173,7 @@ ShearAccelerationRelaxationHourglassControl :: particles_->addVariableToSort("AccelerationHourglass"); } //====================================================================================// -void ShearAccelerationRelaxationHourglassControl::initialization(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl::initialization(size_t index_i, Real dt) { Matd velocity_gradient = Matd::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; @@ -190,7 +190,7 @@ void ShearAccelerationRelaxationHourglassControl::initialization(size_t index_i, velocity_gradient_[index_i] = velocity_gradient; } //====================================================================================// -void ShearAccelerationRelaxationHourglassControl::interaction(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl::interaction(size_t index_i, Real dt) { shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; @@ -201,7 +201,7 @@ void ShearAccelerationRelaxationHourglassControl::interaction(size_t index_i, Re von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); } //====================================================================================// -void ShearAccelerationRelaxationHourglassControl::update(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl::update(size_t index_i, Real dt) { Real rho_i = rho_[index_i]; Matd shear_stress_i = shear_stress_[index_i]; @@ -225,7 +225,7 @@ void ShearAccelerationRelaxationHourglassControl::update(size_t index_i, Real dt //====================================================================================// ShearStressRelaxationHourglassControlJ2Plasticity :: ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) - : ShearAccelerationRelaxationHourglassControl(inner_relation, xi), + : ShearStressRelaxationHourglassControl(inner_relation, xi), J2_plasticity_(DynamicCast(this, particles_->getBaseMaterial())), hardening_factor_(*particles_->registerSharedVariable("HardeningFactor")) { diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 6ae13dcc23..d33bc04289 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -216,15 +216,15 @@ class StressDiffusion : public BasePlasticIntegration void interaction(size_t index_i, Real dt = 0.0); protected: - Real zeta_ = 0.1, fai_; /*diffusion coefficient*/ + Real zeta_ = 0.1, phi_; /*diffusion coefficient*/ Real smoothing_length_, sound_speed_; }; -class ShearAccelerationRelaxationHourglassControl : public fluid_dynamics::BaseIntegration +class ShearStressRelaxationHourglassControl : public fluid_dynamics::BaseIntegration { public: - explicit ShearAccelerationRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi_ = 4.0); - virtual ~ShearAccelerationRelaxationHourglassControl(){}; + explicit ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi_ = 4.0); + virtual ~ShearStressRelaxationHourglassControl(){}; void initialization(size_t index_i, Real dt = 0.0); void interaction(size_t index_i, Real dt = 0.0); void update(size_t index_i, Real dt = 0.0); @@ -237,7 +237,7 @@ class ShearAccelerationRelaxationHourglassControl : public fluid_dynamics::BaseI Real G_, xi_; }; -class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearAccelerationRelaxationHourglassControl +class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearStressRelaxationHourglassControl { public: explicit ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi_ = 0.2); diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp index dc5dbc03ea..64cd3460bd 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -131,7 +131,7 @@ int main(int ac, char *av[]) /** initial condition */ InteractionWithUpdate correction_matrix(beam_body_inner); SimpleDynamics beam_initial_velocity(beam_body); - Dynamics1Level + Dynamics1Level beam_shear_acceleration_relaxation(beam_body_inner); Dynamics1Level beam_pressure_relaxation(beam_body_inner); Dynamics1Level beam_density_relaxation(beam_body_inner); From 2d72b9892a47aca76f3444dce76c512c5e34d2ec Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Mon, 16 Sep 2024 20:33:19 +0200 Subject: [PATCH 12/21] solve the conflicts and update to the latest version --- src/shared/materials/general_continuum.cpp | 10 ++--- src/shared/materials/general_continuum.h | 2 +- .../continuum_integration.cpp | 38 +++++++++---------- .../continuum_integration.h | 18 ++++----- .../oscillating_beam_ULGNOG.cpp | 15 ++++---- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 12 +++--- .../test_3d_taylor_bar_UL/taylor_bar_UL.h | 22 +++++------ 7 files changed, 57 insertions(+), 60 deletions(-) diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index 92c94430d8..65e291f1f9 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -64,9 +64,7 @@ Mat3d PlasticContinuum::ReturnMapping(Mat3d &stress_tensor) { Real stress_tensor_I1 = stress_tensor.trace(); if (-alpha_phi_ * stress_tensor_I1 + k_c_ < 0) - { stress_tensor -= (1.0 / stress_dimension_) * (stress_tensor_I1 - k_c_ / alpha_phi_) * Mat3d::Identity(); - } stress_tensor_I1 = stress_tensor.trace(); Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity(); Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); @@ -100,12 +98,10 @@ Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_ { Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); + Real r = 1.0; if (f > TinyReal) - { - Real r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); - shear_stress = r * shear_stress; - } - return shear_stress; + r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); + return r * shear_stress; } //=================================================================================================// Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) diff --git a/src/shared/materials/general_continuum.h b/src/shared/materials/general_continuum.h index e809ef7f5a..681d652324 100644 --- a/src/shared/materials/general_continuum.h +++ b/src/shared/materials/general_continuum.h @@ -22,7 +22,7 @@ * ------------------------------------------------------------------------- */ /** * @file general_continuum.h - * @brief Describe the linear elastic and Drucker-Prager's plastic model + * @brief Describe the linear elastic, J2 plasticity, and Drucker-Prager's plastic model * @author Shuaihao Zhang and Xiangyu Hu */ diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index d05bb69372..b524c5ec0c 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -11,21 +11,21 @@ ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body) vel_(particles_->registerStateVariable("Velocity")), stress_tensor_3D_(particles_->registerStateVariable("StressTensor3D")) {} //=================================================================================================// -AcousticTimeStepSize::AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL) +AcousticTimeStep::AcousticTimeStep(SPHBody &sph_body, Real acousticCFL) : LocalDynamicsReduce(sph_body), - DataDelegateSimple(sph_body), fluid_(DynamicCast(this, particles_->getBaseMaterial())), - rho_(*particles_->getVariableDataByName("Density")), - p_(*particles_->getVariableDataByName("Pressure")), - vel_(*particles_->getVariableDataByName("Velocity")), + fluid_(DynamicCast(this, particles_->getBaseMaterial())), + rho_(particles_->getVariableDataByName("Density")), + p_(particles_->getVariableDataByName("Pressure")), + vel_(particles_->getVariableDataByName("Velocity")), smoothing_length_min_(sph_body.sph_adaptation_->MinimumSmoothingLength()), acousticCFL_(acousticCFL) {} //=================================================================================================// -Real AcousticTimeStepSize::reduce(size_t index_i, Real dt) +Real AcousticTimeStep::reduce(size_t index_i, Real dt) { return fluid_.getSoundSpeed(p_[index_i], rho_[index_i]) + vel_[index_i].norm(); } //=================================================================================================// -Real AcousticTimeStepSize::outputResult(Real reduced_value) +Real AcousticTimeStep::outputResult(Real reduced_value) { return acousticCFL_ * smoothing_length_min_ / (reduced_value + TinyReal); } @@ -148,17 +148,17 @@ ShearStressRelaxationHourglassControl :: ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi) : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), - shear_stress_(*particles_->registerSharedVariable("ShearStress")), - shear_stress_rate_(*particles_->registerSharedVariable("ShearStressRate")), - velocity_gradient_(*particles_->registerSharedVariable("VelocityGradient")), - strain_tensor_(*particles_->registerSharedVariable("StrainTensor")), - strain_tensor_rate_(*particles_->registerSharedVariable("StrainTensorRate")), - B_(*particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), - von_mises_stress_(*particles_->registerSharedVariable("VonMisesStress")), - von_mises_strain_(*particles_->registerSharedVariable("VonMisesStrain")), - scale_penalty_force_(*particles_->registerSharedVariable("ScalePenaltyForce", Real(1.0))), - acc_shear_(*particles_->registerSharedVariable("AccelerationByShear")), - acc_hourglass_(*particles_->registerSharedVariable("AccelerationHourglass")), + shear_stress_(particles_->registerStateVariable("ShearStress")), + shear_stress_rate_(particles_->registerStateVariable("ShearStressRate")), + velocity_gradient_(particles_->registerStateVariable("VelocityGradient")), + strain_tensor_(particles_->registerStateVariable("StrainTensor")), + strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), + B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), + von_mises_stress_(particles_->registerStateVariable("VonMisesStress")), + von_mises_strain_(particles_->registerStateVariable("VonMisesStrain")), + scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce", Real(1.0))), + acc_shear_(particles_->registerStateVariable("AccelerationByShear")), + acc_hourglass_(particles_->registerStateVariable("AccelerationHourglass")), G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), xi_(xi) { particles_->addVariableToSort("ShearStress"); @@ -227,7 +227,7 @@ ShearStressRelaxationHourglassControlJ2Plasticity :: ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) : ShearStressRelaxationHourglassControl(inner_relation, xi), J2_plasticity_(DynamicCast(this, particles_->getBaseMaterial())), - hardening_factor_(*particles_->registerSharedVariable("HardeningFactor")) + hardening_factor_(particles_->registerStateVariable("HardeningFactor")) { particles_->addVariableToSort("HardeningFactor"); } diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 432966ad22..87613e62ae 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -49,18 +49,18 @@ class ContinuumInitialCondition : public LocalDynamics Mat3d *stress_tensor_3D_; }; -class AcousticTimeStepSize : public LocalDynamicsReduce, public DataDelegateSimple +class AcousticTimeStep : public LocalDynamicsReduce { public: - explicit AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL = 0.6); - virtual ~AcousticTimeStepSize(){}; + explicit AcousticTimeStep(SPHBody &sph_body, Real acousticCFL = 0.6); + virtual ~AcousticTimeStep(){}; Real reduce(size_t index_i, Real dt = 0.0); virtual Real outputResult(Real reduced_value) override; protected: Fluid &fluid_; - StdLargeVec &rho_, &p_; - StdLargeVec &vel_; + Real *rho_, *p_; + Vecd *vel_; Real smoothing_length_min_; Real acousticCFL_; }; @@ -231,9 +231,9 @@ class ShearStressRelaxationHourglassControl : public fluid_dynamics::BaseIntegra protected: GeneralContinuum &continuum_; - StdLargeVec &shear_stress_, &shear_stress_rate_, &velocity_gradient_, &strain_tensor_, &strain_tensor_rate_, &B_; - StdLargeVec &von_mises_stress_, &von_mises_strain_, &scale_penalty_force_; - StdLargeVec &acc_shear_, &acc_hourglass_; + Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_; + Real *von_mises_stress_, *von_mises_strain_, *scale_penalty_force_; + Vecd *acc_shear_, *acc_hourglass_; Real G_, xi_; }; @@ -246,7 +246,7 @@ class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearStressRela protected: J2Plasticity &J2_plasticity_; - StdLargeVec &hardening_factor_; + Real *hardening_factor_; }; } // namespace continuum_dynamics } // namespace SPH diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp index 64cd3460bd..d939a0a66b 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -136,8 +136,8 @@ int main(int ac, char *av[]) Dynamics1Level beam_pressure_relaxation(beam_body_inner); Dynamics1Level beam_density_relaxation(beam_body_inner); SimpleDynamics beam_volume_update(beam_body); - ReduceDynamics fluid_advection_time_step(beam_body, U_ref, 0.2); - ReduceDynamics fluid_acoustic_time_step(beam_body, 0.4); + ReduceDynamics advection_time_step(beam_body, U_ref, 0.2); + ReduceDynamics acoustic_time_step(beam_body, 0.4); // clamping a solid body part. BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); SimpleDynamics constraint_beam_base(beam_base); @@ -161,6 +161,7 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // Setup computing time-step controls. //---------------------------------------------------------------------- + Real &physical_time = *sph_system.getSystemVariableDataByName("PhysicalTime"); int ite = 0; Real T0 = 1.0; Real End_Time = T0; @@ -174,17 +175,17 @@ int main(int ac, char *av[]) write_beam_tip_displacement.writeToFile(0); write_beam_kinetic_energy.writeToFile(0); // computation loop starts - while (GlobalStaticVariables::physical_time_ < End_Time) + while (physical_time < End_Time) { Real integration_time = 0.0; while (integration_time < D_Time) { Real relaxation_time = 0.0; - Real advection_dt = fluid_advection_time_step.exec(); + Real advection_dt = advection_time_step.exec(); beam_volume_update.exec(); while (relaxation_time < advection_dt) { - Real acoustic_dt = fluid_acoustic_time_step.exec(); + Real acoustic_dt = acoustic_time_step.exec(); beam_pressure_relaxation.exec(acoustic_dt); constraint_beam_base.exec(); beam_shear_acceleration_relaxation.exec(acoustic_dt); @@ -192,11 +193,11 @@ int main(int ac, char *av[]) ite++; relaxation_time += acoustic_dt; integration_time += acoustic_dt; - GlobalStaticVariables::physical_time_ += acoustic_dt; + physical_time += acoustic_dt; if (ite % 500 == 0) { std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ << " advection_dt: " + << physical_time << " advection_dt: " << advection_dt << " acoustic_dt: " << acoustic_dt << "\n"; } diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 7c4f61dc0e..3bd1e8b4f4 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -81,8 +81,8 @@ int main(int ac, char *av[]) Dynamics1Level column_pressure_relaxation(column_inner); Dynamics1Level column_density_relaxation(column_inner); SimpleDynamics beam_volume_update(column); - ReduceDynamics advection_time_step(column, U_max, 0.2); - ReduceDynamics acoustic_time_step(column, 0.4); + ReduceDynamics advection_time_step(column, U_max, 0.2); + ReduceDynamics acoustic_time_step(column, 0.4); InteractionDynamics column_wall_contact_force(column_wall_contact); //---------------------------------------------------------------------- // Output @@ -95,7 +95,6 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // From here the time stepping begins. //---------------------------------------------------------------------- - GlobalStaticVariables::physical_time_ = 0.0; system.initializeSystemCellLinkedLists(); system.initializeSystemConfigurations(); wall_normal_direction.exec(); @@ -104,6 +103,7 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // Setup time-stepping related simulation parameters. //---------------------------------------------------------------------- + Real &physical_time = *system.getSystemVariableDataByName("PhysicalTime"); int ite = 0; Real end_time = 6.0e-5; int screen_output_interval = 100; @@ -119,7 +119,7 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // Main time-stepping loop. //---------------------------------------------------------------------- - while (GlobalStaticVariables::physical_time_ < end_time) + while (physical_time < end_time) { Real integration_time = 0.0; while (integration_time < output_period) @@ -133,7 +133,7 @@ int main(int ac, char *av[]) if (ite % screen_output_interval == 0) { std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ << " advection_dt: " + << physical_time << " advection_dt: " << advection_dt << " acoustic_dt: " << acoustic_dt << "\n"; } @@ -146,7 +146,7 @@ int main(int ac, char *av[]) ite++; relaxation_time += acoustic_dt; integration_time += acoustic_dt; - GlobalStaticVariables::physical_time_ += acoustic_dt; + physical_time += acoustic_dt; } column.updateCellLinkedList(); column_inner.updateConfiguration(); diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h index d7a1f200f7..2e3e3df898 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.h @@ -74,9 +74,9 @@ class DynamicContactForceWithWall : public LocalDynamics, : LocalDynamics(solid_body_contact_relation.getSPHBody()), DataDelegateContact(solid_body_contact_relation), continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), - Vol_(*particles_->getVariableDataByName("VolumetricMeasure")), - vel_(*particles_->getVariableDataByName("Velocity")), - force_prior_(*particles_->getVariableDataByName("ForcePrior")), + Vol_(particles_->getVariableDataByName("VolumetricMeasure")), + vel_(particles_->getVariableDataByName("Velocity")), + force_prior_(particles_->getVariableDataByName("ForcePrior")), penalty_strength_(penalty_strength) { impedance_ = continuum_.ReferenceDensity() * sqrt(continuum_.ContactStiffness()); @@ -84,7 +84,7 @@ class DynamicContactForceWithWall : public LocalDynamics, for (size_t k = 0; k != contact_particles_.size(); ++k) { contact_Vol_.push_back(contact_particles_[k]->getVariableDataByName("VolumetricMeasure")); - contact_vel_.push_back(contact_particles_[k]->registerSharedVariable("Velocity")); + contact_vel_.push_back(contact_particles_[k]->registerStateVariable("Velocity")); contact_n_.push_back(contact_particles_[k]->template getVariableDataByName("NormalDirection")); } }; @@ -99,9 +99,9 @@ class DynamicContactForceWithWall : public LocalDynamics, 1.0 / (sph_body_.sph_adaptation_->ReferenceSpacing() * particle_spacing_j1); particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2; - StdLargeVec &n_k = *(contact_n_[k]); - StdLargeVec &vel_n_k = *(contact_vel_[k]); - StdLargeVec &Vol_k = *(contact_Vol_[k]); + Vecd *n_k = contact_n_[k]; + Vecd *vel_n_k = contact_vel_[k]; + Real *Vol_k = contact_Vol_[k]; Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i]; for (size_t n = 0; n != contact_neighborhood.current_size_; ++n) { @@ -123,10 +123,10 @@ class DynamicContactForceWithWall : public LocalDynamics, protected: GeneralContinuum &continuum_; - StdLargeVec &Vol_; - StdLargeVec &vel_, &force_prior_; // note that prior force directly used here - StdVec *> contact_Vol_; - StdVec *> contact_vel_, contact_n_; + Real *Vol_; + Vecd *vel_, *force_prior_; // note that prior force directly used here + StdVec contact_Vol_; + StdVec contact_vel_, contact_n_; Real penalty_strength_; Real impedance_, reference_pressure_; }; From 73848c258528b9c81006856244edd196eb7e3f39 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 00:03:04 +0200 Subject: [PATCH 13/21] fix regression test --- .../regression_test_tool/MyObserver_Position_dtwdistance.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml index d8b65a8b22..a585fb7aee 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml @@ -1,4 +1,4 @@ - + From 64c701977968554fe273ccd7a7c6f7aed570fab7 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 05:51:53 +0200 Subject: [PATCH 14/21] fix regression test --- .../regression_test_tool/MyObserver_Position_dtwdistance.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml index a585fb7aee..bcf6137195 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml +++ b/tests/3d_examples/test_3d_taylor_bar_UL/regression_test_tool/MyObserver_Position_dtwdistance.xml @@ -1,4 +1,4 @@ - + From d54b327ab9c6f4dae73352ca4ccc43f059728723 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 14:22:15 +0200 Subject: [PATCH 15/21] delete unnecessary variables --- src/shared/materials/general_continuum.cpp | 10 +++------- .../continuum_dynamics/continuum_integration.cpp | 6 ++---- .../oscillating_beam_UL.cpp | 1 - .../oscillating_beam_ULGNOG.cpp | 1 - .../oscillating_plate_UL.cpp | 3 --- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 9 ++++----- 6 files changed, 9 insertions(+), 21 deletions(-) diff --git a/src/shared/materials/general_continuum.cpp b/src/shared/materials/general_continuum.cpp index 65e291f1f9..d799cdfec8 100644 --- a/src/shared/materials/general_continuum.cpp +++ b/src/shared/materials/general_continuum.cpp @@ -43,16 +43,15 @@ Mat3d PlasticContinuum::ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &st Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose()); Mat3d deviatoric_strain_rate = strain_rate - (1.0 / stress_dimension_) * strain_rate.trace() * Mat3d::Identity(); Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor; - Real stress_tensor_I1 = stress_tensor.trace(); Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity(); Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum(); - Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor_I1 - k_c_; + Real f = sqrt(stress_tensor_J2) + alpha_phi_ * stress_tensor.trace() - k_c_; Real lambda_dot_ = 0; Mat3d g = Mat3d::Zero(); if (f >= TinyReal) { Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum(); - //non-associate flow rule + // non-associate flow rule lambda_dot_ = (3.0 * alpha_phi_ * K_ * strain_rate.trace() + (G_ / sqrt(stress_tensor_J2)) * deviatoric_stress_times_strain_rate) / (9.0 * alpha_phi_ * K_ * getDPConstantsA(psi_) + G_); g = lambda_dot_ * (3.0 * K_ * getDPConstantsA(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2))); } @@ -107,11 +106,8 @@ Matd J2Plasticity::ReturnMappingShearStress(Matd &shear_stress, Real &hardening_ Real J2Plasticity::ScalePenaltyForce(Matd &shear_stress, Real &hardening_factor) { Real stress_tensor_J2 = 0.5 * (shear_stress.cwiseProduct(shear_stress.transpose())).sum(); - Real r = 1.0; Real f = sqrt(2.0 * stress_tensor_J2) - sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_); - if (f > TinyReal) - r = (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal); - return r; + return (f > TinyReal) ? (sqrt_2_over_3_ * (hardening_modulus_ * hardening_factor + yield_stress_)) / (sqrt(2.0 * stress_tensor_J2) + TinyReal) : 1.0; } //=================================================================================================// Real J2Plasticity::HardeningFactorRate(const Matd &shear_stress, Real &hardening_factor) diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index b524c5ec0c..f5df025a98 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -99,11 +99,9 @@ void ShearStressRelaxation::interaction(size_t index_i, Real dt) } velocity_gradient_[index_i] = velocity_gradient; /*calculate strain*/ - Matd strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose()); - strain_tensor_rate_[index_i] = strain_rate; + strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; - Matd strain_i = strain_tensor_[index_i]; - von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_i); + von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); } //====================================================================================// void ShearStressRelaxation::update(size_t index_i, Real dt) diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index bfdbec1cee..92aa21433b 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -27,7 +27,6 @@ Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus Real poisson = 0.3975; // Poisson ratio Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -Real gravity_g = 0.0; //---------------------------------------------------------------------- // Parameters for initial condition on velocity //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp index d939a0a66b..76beeaca52 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -28,7 +28,6 @@ Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus Real poisson = 0.3975; // Poisson ratio Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -Real gravity_g = 0.0; //---------------------------------------------------------------------- // Parameters for initial condition on velocity //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp index dceabddb3d..01a4153ded 100644 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp +++ b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp @@ -27,7 +27,6 @@ Real rho0_s = 1000.0; /** Normalized density. */ Real Youngs_modulus = 100.0e6; /** Normalized Youngs Modulus. */ Real poisson = 0.3; /** Poisson ratio. */ Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -Real gravity_g = 0.0; Real governing_vibration_integer_x = 2.0; Real governing_vibration_integer_y = 2.0; @@ -147,12 +146,10 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- SimpleDynamics initial_velocity(plate_body); InteractionWithUpdate corrected_configuration(plate_body_inner); - InteractionDynamics plate_shear_acceleration(plate_body_inner); Dynamics1Level plate_pressure_relaxation(plate_body_inner); Dynamics1Level plate_density_relaxation(plate_body_inner); Dynamics1Level plate_shear_stress_relaxation(plate_body_inner); - ReduceDynamics fluid_advection_time_step(plate_body, U_ref, 0.2); ReduceDynamics fluid_acoustic_time_step(plate_body, 0.4); BoundaryGeometry boundary_geometry(plate_body, "BoundaryGeometry"); diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 3bd1e8b4f4..464f0a3535 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -31,8 +31,6 @@ int main(int ac, char *av[]) InnerRelation column_inner(column); ContactRelation my_observer_contact(my_observer, {&column}); SurfaceContactRelation column_wall_contact(column, {&wall_boundary}); - /**define simple data file input and outputs functions. */ - BodyStatesRecordingToVtp write_states(system); //---------------------------------------------------------------------- // Run particle relaxation for body-fitted distribution if chosen. //---------------------------------------------------------------------- @@ -42,7 +40,7 @@ int main(int ac, char *av[]) /** Random reset the insert body particle position. */ SimpleDynamics random_column_particles(column); /** Write the body state to Vtp file. */ - BodyStatesRecordingToVtp write_column_to_vtp(column); + BodyStatesRecordingToVtp write_column_to_vtp(system); /** Write the particle reload files. */ ReloadParticleIO write_particle_reload_files(column); /** A Physics relaxation step. */ @@ -52,8 +50,7 @@ int main(int ac, char *av[]) */ random_column_particles.exec(0.25); relaxation_step_inner.SurfaceBounding().exec(); - write_states.writeToFile(0.0); - + write_column_to_vtp.writeToFile(0); /** relax particles of the insert body. */ int ite_p = 0; while (ite_p < 1000) @@ -87,6 +84,8 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- // Output //---------------------------------------------------------------------- + /**define simple data file input and outputs functions. */ + BodyStatesRecordingToVtp write_states(system); write_states.addToWrite(column, "VonMisesStress"); write_states.addToWrite(column, "VonMisesStrain"); write_states.addToWrite(column, "Pressure"); From 1c95bd1f6f20391f13f1048ce42bfaae8c7a8d94 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 17:45:08 +0200 Subject: [PATCH 16/21] Put the variables that will not be involved in the calculation into derived variables. --- .../all_continuum_dynamics.h | 1 + .../continuum_dynamics_variable.cpp | 57 +++++++++++ .../continuum_dynamics_variable.h | 99 +++++++++++++++++++ .../continuum_integration.cpp | 22 +---- .../continuum_integration.h | 12 +-- .../continuum_integration.hpp | 32 +----- .../column_collapse.cpp | 7 +- .../oscillating_beam_UL.cpp | 5 +- .../oscillating_beam_ULGNOG.cpp | 5 +- .../oscillating_plate_UL.cpp | 5 +- .../test_3d_repose_angle/repose_angle.cpp | 4 + .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 9 +- 12 files changed, 186 insertions(+), 72 deletions(-) create mode 100644 src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp create mode 100644 src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h diff --git a/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h b/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h index dbcc5ee06e..328c11e76a 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h +++ b/src/shared/particle_dynamics/continuum_dynamics/all_continuum_dynamics.h @@ -30,3 +30,4 @@ #include "base_continuum_dynamics.h" #include "continuum_integration.hpp" +#include "continuum_dynamics_variable.h" diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp new file mode 100644 index 0000000000..c084c36ba4 --- /dev/null +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.cpp @@ -0,0 +1,57 @@ +#include "continuum_dynamics_variable.h" + +namespace SPH +{ +namespace continuum_dynamics +{ +//=============================================================================================// +VonMisesStress::VonMisesStress(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VonMisesStress"), + p_(particles_->getVariableDataByName("Pressure")), + shear_stress_(particles_->getVariableDataByName("ShearStress")) {} +//=============================================================================================// +void VonMisesStress::update(size_t index_i, Real dt) +{ + Matd stress_tensor = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); + derived_variable_[index_i] = getVonMisesStressFromMatrix(stress_tensor); +} +//=============================================================================================// +VonMisesStrain::VonMisesStrain(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VonMisesStrain"), + strain_tensor_(particles_->getVariableDataByName("StrainTensor")) {} +//=============================================================================================// +void VonMisesStrain::update(size_t index_i, Real dt) +{ + derived_variable_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); +} +//=============================================================================================// +VerticalStress::VerticalStress(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "VerticalStress"), + stress_tensor_3D_(particles_->getVariableDataByName("StressTensor3D")) {} +//=============================================================================================// +void VerticalStress::update(size_t index_i, Real dt) +{ + derived_variable_[index_i] = stress_tensor_3D_[index_i](1, 1); +} +//=============================================================================================// +AccDeviatoricPlasticStrain::AccDeviatoricPlasticStrain(SPHBody &sph_body) + : BaseDerivedVariable(sph_body, "AccDeviatoricPlasticStrain"), + plastic_continuum_(DynamicCast(this, sph_body_.getBaseMaterial())), + stress_tensor_3D_(particles_->getVariableDataByName("StressTensor3D")), + strain_tensor_3D_(particles_->getVariableDataByName("StrainTensor3D")), + E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) {} +//=============================================================================================// +void AccDeviatoricPlasticStrain::update(size_t index_i, Real dt) +{ + Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity(); + Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace(); + Mat3d elastic_strain_tensor_3D = deviatoric_stress / (2.0 * plastic_continuum_.getShearModulus(E_, nu_)) + + hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_continuum_.getBulkModulus(E_, nu_)); + Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D; + Mat3d deviatoric_strain_tensor = plastic_strain_tensor_3D - (1.0 / (Real)Dimensions) * plastic_strain_tensor_3D.trace() * Mat3d::Identity(); + Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum(); + derived_variable_[index_i] = sqrt(sum * 2.0 / 3.0); +} +//=================================================================================================// +} // namespace continuum_dynamics +} // namespace SPH diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h new file mode 100644 index 0000000000..54fd194b37 --- /dev/null +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_dynamics_variable.h @@ -0,0 +1,99 @@ +/* ------------------------------------------------------------------------- * + * SPHinXsys * + * ------------------------------------------------------------------------- * + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle * + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic * + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4. * + * * + * Portions copyright (c) 2017-2023 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a * + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------- */ +/** + * @file continuum_dynamics_variable.h + * @brief Here, we define the algorithm classes for computing derived solid dynamics variables. + * @details These variable can be added into variable list for state output. + * @author Shuaihao Zhang and Xiangyu Hu + */ + +#ifndef CONTINUUM_DYNAMICS_VARIABLE_H +#define CONTINUUM_DYNAMICS_VARIABLE_H + +#include "base_general_dynamics.h" +#include "general_continuum.h" + +namespace SPH +{ +namespace continuum_dynamics +{ +/** + * @class VonMisesStress + * @brief computing von_Mises_stress + */ +class VonMisesStress : public BaseDerivedVariable +{ + public: + explicit VonMisesStress(SPHBody &sph_body); + virtual ~VonMisesStress(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Real *p_; + Matd *shear_stress_; +}; +/** + * @class VonMisesStrain + * @brief computing von_Mises_strain + */ +class VonMisesStrain : public BaseDerivedVariable +{ + public: + explicit VonMisesStrain(SPHBody &sph_body); + virtual ~VonMisesStrain(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Matd *strain_tensor_; +}; +/** + * @class VerticalStress + */ +class VerticalStress : public BaseDerivedVariable +{ + public: + explicit VerticalStress(SPHBody &sph_body); + virtual ~VerticalStress(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + Mat3d *stress_tensor_3D_; +}; +/** + * @class AccumulatedDeviatoricPlasticStrain + */ +class AccDeviatoricPlasticStrain : public BaseDerivedVariable +{ + public: + explicit AccDeviatoricPlasticStrain(SPHBody &sph_body); + virtual ~AccDeviatoricPlasticStrain(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + PlasticContinuum &plastic_continuum_; + Mat3d *stress_tensor_3D_, *strain_tensor_3D_; + Real E_, nu_; +}; +} // namespace continuum_dynamics +} // namespace SPH +#endif // CONTINUUM_DYNAMICS_VARIABLE_H \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index f5df025a98..5921d97634 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -65,15 +65,11 @@ ShearStressRelaxation::ShearStressRelaxation(BaseInnerRelation &inner_relation) velocity_gradient_(particles_->registerStateVariable("VelocityGradient")), strain_tensor_(particles_->registerStateVariable("StrainTensor")), strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), - von_mises_stress_(particles_->registerStateVariable("VonMisesStress")), - von_mises_strain_(particles_->registerStateVariable("VonMisesStrain")), - Vol_(particles_->getVariableDataByName("VolumetricMeasure")), - B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")) + B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), + Vol_(particles_->getVariableDataByName("VolumetricMeasure")) { particles_->addVariableToSort("ShearStress"); particles_->addVariableToSort("ShearStressRate"); - particles_->addVariableToSort("VonMisesStress"); - particles_->addVariableToSort("VonMisesStrain"); particles_->addVariableToSort("VelocityGradient"); particles_->addVariableToSort("StrainTensor"); particles_->addVariableToSort("StrainTensorRate"); @@ -98,18 +94,14 @@ void ShearStressRelaxation::interaction(size_t index_i, Real dt) velocity_gradient -= v_ij * (B_[index_i] * e_ij * dW_ijV_j).transpose(); } velocity_gradient_[index_i] = velocity_gradient; - /*calculate strain*/ strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; - von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); } //====================================================================================// void ShearStressRelaxation::update(size_t index_i, Real dt) { shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); shear_stress_[index_i] += shear_stress_rate_[index_i] * dt * 0.5; - Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); - von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); } //====================================================================================// StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation) @@ -152,8 +144,6 @@ ShearStressRelaxationHourglassControl :: strain_tensor_(particles_->registerStateVariable("StrainTensor")), strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), - von_mises_stress_(particles_->registerStateVariable("VonMisesStress")), - von_mises_strain_(particles_->registerStateVariable("VonMisesStrain")), scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce", Real(1.0))), acc_shear_(particles_->registerStateVariable("AccelerationByShear")), acc_hourglass_(particles_->registerStateVariable("AccelerationHourglass")), @@ -164,8 +154,6 @@ ShearStressRelaxationHourglassControl :: particles_->addVariableToSort("VelocityGradient"); particles_->addVariableToSort("StrainTensor"); particles_->addVariableToSort("StrainTensorRate"); - particles_->addVariableToSort("VonMisesStress"); - particles_->addVariableToSort("VonMisesStrain"); particles_->addVariableToSort("ScalePenaltyForce"); particles_->addVariableToSort("AccelerationByShear"); particles_->addVariableToSort("AccelerationHourglass"); @@ -192,11 +180,8 @@ void ShearStressRelaxationHourglassControl::interaction(size_t index_i, Real dt) { shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; - Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); - von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt; - von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); } //====================================================================================// void ShearStressRelaxationHourglassControl::update(size_t index_i, Real dt) @@ -238,11 +223,8 @@ void ShearStressRelaxationHourglassControlJ2Plasticity::interaction(size_t index hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment; scale_penalty_force_[index_i] = J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); shear_stress_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_try, hardening_factor_[index_i]); - Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); - von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); Matd strain_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); strain_tensor_[index_i] += strain_rate * dt; - von_mises_strain_[index_i] = getVonMisesStressFromMatrix(strain_tensor_[index_i]); } } // namespace continuum_dynamics } // namespace SPH \ No newline at end of file diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 87613e62ae..cc1464c839 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -103,9 +103,8 @@ class ShearStressRelaxation : public fluid_dynamics::BaseIntegration @@ -119,7 +118,6 @@ class BasePlasticIntegration : public fluid_dynamics::BaseIntegration, RiemannSolverType> protected: RiemannSolverType riemann_solver_; - Real *acc_deviatoric_plastic_strain_, *vertical_stress_; Real *Vol_, *mass_; - Real E_, nu_; - - Real getDeviatoricPlasticStrain(Mat3d &strain_tensor); }; using PlasticIntegration2ndHalfInnerNoRiemann = PlasticIntegration2ndHalf, NoRiemannSolver>; using PlasticIntegration2ndHalfInnerRiemann = PlasticIntegration2ndHalf, AcousticRiemannSolver>; @@ -232,7 +226,7 @@ class ShearStressRelaxationHourglassControl : public fluid_dynamics::BaseIntegra protected: GeneralContinuum &continuum_; Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_; - Real *von_mises_stress_, *von_mises_strain_, *scale_penalty_force_; + Real *scale_penalty_force_; Vecd *acc_shear_, *acc_hourglass_; Real G_, xi_; }; diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp index fa9ca9243b..87e50b6d69 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp @@ -28,12 +28,8 @@ BasePlasticIntegration::BasePlasticIntegration(BaseRelationT strain_tensor_3D_(this->particles_->template registerStateVariable("StrainTensor3D")), stress_rate_3D_(this->particles_->template registerStateVariable("StressRate3D")), strain_rate_3D_(this->particles_->template registerStateVariable("StrainRate3D")), - elastic_strain_tensor_3D_(this->particles_->template registerStateVariable("ElasticStrainTensor3D")), - elastic_strain_rate_3D_(this->particles_->template registerStateVariable("ElasticStrainRate3D")), velocity_gradient_(this->particles_->template registerStateVariable("VelocityGradient")) { - this->particles_->template addVariableToSort("ElasticStrainTensor3D"); - this->particles_->template addVariableToSort("ElasticStrainRate3D"); this->particles_->template addVariableToSort("StrainTensor3D"); this->particles_->template addVariableToSort("StressTensor3D"); this->particles_->template addVariableToSort("StrainRate3D"); @@ -142,15 +138,8 @@ template PlasticIntegration2ndHalf, RiemannSolverType>::PlasticIntegration2ndHalf(BaseInnerRelation &inner_relation) : BasePlasticIntegration(inner_relation), riemann_solver_(plastic_continuum_, plastic_continuum_, 20.0 * (Real)Dimensions), - acc_deviatoric_plastic_strain_(particles_->registerStateVariable("AccDeviatoricPlasticStrain")), - vertical_stress_(particles_->registerStateVariable("VerticalStress")), Vol_(particles_->getVariableDataByName("VolumetricMeasure")), - mass_(particles_->getVariableDataByName("Mass")), - E_(plastic_continuum_.getYoungsModulus()), nu_(plastic_continuum_.getPoissonRatio()) -{ - particles_->addVariableToSort("AccDeviatoricPlasticStrain"); - particles_->addVariableToSort("VerticalStress"); -} + mass_(particles_->getVariableDataByName("Mass")) {} //=================================================================================================// template void PlasticIntegration2ndHalf, RiemannSolverType>::initialization(size_t index_i, Real dt) @@ -181,14 +170,6 @@ void PlasticIntegration2ndHalf, RiemannSolverType>::interaction(size_t i } //=================================================================================================// template -Real PlasticIntegration2ndHalf, RiemannSolverType>::getDeviatoricPlasticStrain(Mat3d &strain_tensor) -{ - Mat3d deviatoric_strain_tensor = strain_tensor - (1.0 / (Real)Dimensions) * strain_tensor.trace() * Mat3d::Identity(); - Real sum = (deviatoric_strain_tensor.cwiseProduct(deviatoric_strain_tensor)).sum(); - return sqrt(sum * 2.0 / 3.0); -} -//=================================================================================================// -template void PlasticIntegration2ndHalf, RiemannSolverType>::update(size_t index_i, Real dt) { rho_[index_i] += drho_dt_[index_i] * dt * 0.5; @@ -198,18 +179,9 @@ void PlasticIntegration2ndHalf, RiemannSolverType>::update(size_t index_ stress_rate_3D_[index_i] += stress_tensor_rate_3D_; stress_tensor_3D_[index_i] += stress_rate_3D_[index_i] * dt; /*return mapping*/ - Mat3d stress_tensor_ = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); - stress_tensor_3D_[index_i] = stress_tensor_; - vertical_stress_[index_i] = stress_tensor_3D_[index_i](1, 1); + stress_tensor_3D_[index_i] = plastic_continuum_.ReturnMapping(stress_tensor_3D_[index_i]); strain_rate_3D_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); strain_tensor_3D_[index_i] += strain_rate_3D_[index_i] * dt; - /*calculate elastic strain*/ - Mat3d deviatoric_stress = stress_tensor_3D_[index_i] - (1.0 / 3.0) * stress_tensor_3D_[index_i].trace() * Mat3d::Identity(); - Real hydrostatic_pressure = (1.0 / 3.0) * stress_tensor_3D_[index_i].trace(); - elastic_strain_tensor_3D_[index_i] = deviatoric_stress / (2.0 * plastic_continuum_.getShearModulus(E_, nu_)) + - hydrostatic_pressure * Mat3d::Identity() / (9.0 * plastic_continuum_.getBulkModulus(E_, nu_)); - Mat3d plastic_strain_tensor_3D = strain_tensor_3D_[index_i] - elastic_strain_tensor_3D_[index_i]; - acc_deviatoric_plastic_strain_[index_i] = getDeviatoricPlasticStrain(plastic_strain_tensor_3D); } //=================================================================================================// template diff --git a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp index 650ed6d87c..65f74ec306 100644 --- a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp +++ b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp @@ -99,12 +99,10 @@ int main(int ac, char *av[]) Gravity gravity(Vecd(0.0, -gravity_g)); SimpleDynamics> constant_gravity(soil_block, gravity); SimpleDynamics wall_boundary_normal_direction(wall_boundary); - Dynamics1Level granular_stress_relaxation(soil_block_inner, soil_block_contact); Dynamics1Level granular_density_relaxation(soil_block_inner, soil_block_contact); InteractionWithUpdate soil_density_by_summation(soil_block_inner, soil_block_contact); InteractionDynamics stress_diffusion(soil_block_inner); - ReduceDynamics soil_acoustic_time_step(soil_block, 0.4); //---------------------------------------------------------------------- // Define the methods for I/O operations, observations @@ -113,9 +111,10 @@ int main(int ac, char *av[]) BodyStatesRecordingToVtp body_states_recording(sph_system); body_states_recording.addToWrite(soil_block, "Pressure"); body_states_recording.addToWrite(soil_block, "Density"); + SimpleDynamics vertical_stress(soil_block); body_states_recording.addToWrite(soil_block, "VerticalStress"); + SimpleDynamics accumulated_deviatoric_plastic_strain(soil_block); body_states_recording.addToWrite(soil_block, "AccDeviatoricPlasticStrain"); - body_states_recording.addToWrite(wall_boundary, "NormalDirection"); RestartIO restart_io(sph_system); RegressionTestDynamicTimeWarping> write_mechanical_energy(soil_block, gravity); @@ -205,6 +204,8 @@ int main(int ac, char *av[]) interval_updating_configuration += TickCount::now() - time_instance; } TickCount t2 = TickCount::now(); + vertical_stress.exec(); + accumulated_deviatoric_plastic_strain.exec(); body_states_recording.writeToFile(); TickCount t3 = TickCount::now(); interval += t3 - t2; diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index 92aa21433b..466b2bb589 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -146,10 +146,10 @@ int main(int ac, char *av[]) // outputs //----------------------------------------------------------------------------- BodyStatesRecordingToVtp write_beam_states(beam_body); - write_beam_states.addToWrite(beam_body, "VonMisesStress"); - write_beam_states.addToWrite(beam_body, "VonMisesStrain"); write_beam_states.addToWrite(beam_body, "Density"); write_beam_states.addToWrite(beam_body, "Pressure"); + SimpleDynamics beam_von_mises_stress(beam_body); + write_beam_states.addToWrite(beam_body, "VonMisesStress"); ObservedQuantityRecording write_beam_tip_displacement("Position", beam_observer_contact); RegressionTestDynamicTimeWarping> write_beam_kinetic_energy(beam_body); //---------------------------------------------------------------------- @@ -209,6 +209,7 @@ int main(int ac, char *av[]) beam_body_inner.updateConfiguration(); correction_matrix.exec(); } + beam_von_mises_stress.exec(); write_beam_tip_displacement.writeToFile(ite); write_beam_kinetic_energy.writeToFile(ite); TickCount t2 = TickCount::now(); diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp index 76beeaca52..505320bc43 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -144,10 +144,10 @@ int main(int ac, char *av[]) // outputs //----------------------------------------------------------------------------- BodyStatesRecordingToVtp write_beam_states(beam_body); - write_beam_states.addToWrite(beam_body, "VonMisesStress"); - write_beam_states.addToWrite(beam_body, "VonMisesStrain"); write_beam_states.addToWrite(beam_body, "Density"); write_beam_states.addToWrite(beam_body, "Pressure"); + SimpleDynamics beam_von_mises_stress(beam_body); + write_beam_states.addToWrite(beam_body, "VonMisesStress"); ObservedQuantityRecording write_beam_tip_displacement("Position", beam_observer_contact); RegressionTestDynamicTimeWarping> write_beam_kinetic_energy(beam_body); //---------------------------------------------------------------------- @@ -205,6 +205,7 @@ int main(int ac, char *av[]) beam_body_inner.updateConfiguration(); correction_matrix.exec(); } + beam_von_mises_stress.exec(); write_beam_tip_displacement.writeToFile(ite); write_beam_kinetic_energy.writeToFile(ite); TickCount t2 = TickCount::now(); diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp index 01a4153ded..e0b6da44d0 100644 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp +++ b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp @@ -160,10 +160,10 @@ int main(int ac, char *av[]) // and regression tests of the simulation. //---------------------------------------------------------------------- BodyStatesRecordingToVtp write_states(sph_system); - write_states.addToWrite(plate_body, "VonMisesStress"); - write_states.addToWrite(plate_body, "VonMisesStress"); write_states.addToWrite(plate_body, "Pressure"); write_states.addToWrite(plate_body, "Density"); + SimpleDynamics plate_von_mises_stress(plate_body); + write_states.addToWrite(plate_body, "VonMisesStress"); RestartIO restart_io(sph_system); ObservedQuantityRecording write_plate_displacement("Position", plate_observer_contact); RegressionTestDynamicTimeWarping> write_kinetic_energy(plate_body); @@ -233,6 +233,7 @@ int main(int ac, char *av[]) plate_body_inner.updateConfiguration(); corrected_configuration.exec(); } + plate_von_mises_stress.exec(); write_plate_displacement.writeToFile(number_of_iterations); write_kinetic_energy.writeToFile(number_of_iterations); TickCount t2 = TickCount::now(); diff --git a/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp b/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp index c80282c26c..dae594254e 100644 --- a/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp +++ b/tests/3d_examples/test_3d_repose_angle/repose_angle.cpp @@ -171,7 +171,9 @@ int main(int ac, char *av[]) BodyStatesRecordingToVtp body_states_recording(sph_system); body_states_recording.addToWrite(soil_block, "Density"); body_states_recording.addToWrite(soil_block, "Pressure"); + SimpleDynamics vertical_stress(soil_block); body_states_recording.addToWrite(soil_block, "VerticalStress"); + SimpleDynamics accumulated_deviatoric_plastic_strain(soil_block); body_states_recording.addToWrite(soil_block, "AccDeviatoricPlasticStrain"); RestartIO restart_io(sph_system); RegressionTestDynamicTimeWarping> write_soil_mechanical_energy(soil_block, gravity); @@ -263,6 +265,8 @@ int main(int ac, char *av[]) interval_updating_configuration += TickCount::now() - time_instance; } TickCount t2 = TickCount::now(); + vertical_stress.exec(); + accumulated_deviatoric_plastic_strain.exec(); body_states_recording.writeToFile(); TickCount t3 = TickCount::now(); interval += t3 - t2; diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 464f0a3535..40ece24af5 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -77,7 +77,7 @@ int main(int ac, char *av[]) Dynamics1Level column_shear_stress_relaxation(column_inner); Dynamics1Level column_pressure_relaxation(column_inner); Dynamics1Level column_density_relaxation(column_inner); - SimpleDynamics beam_volume_update(column); + SimpleDynamics column_volume_update(column); ReduceDynamics advection_time_step(column, U_max, 0.2); ReduceDynamics acoustic_time_step(column, 0.4); InteractionDynamics column_wall_contact_force(column_wall_contact); @@ -86,10 +86,10 @@ int main(int ac, char *av[]) //---------------------------------------------------------------------- /**define simple data file input and outputs functions. */ BodyStatesRecordingToVtp write_states(system); - write_states.addToWrite(column, "VonMisesStress"); - write_states.addToWrite(column, "VonMisesStrain"); write_states.addToWrite(column, "Pressure"); write_states.addToWrite(column, "Density"); + SimpleDynamics column_von_mises_stress(column); + write_states.addToWrite(column, "VonMisesStress"); RegressionTestDynamicTimeWarping> write_displacement("Position", my_observer_contact); //---------------------------------------------------------------------- // From here the time stepping begins. @@ -125,7 +125,7 @@ int main(int ac, char *av[]) { Real relaxation_time = 0.0; Real advection_dt = advection_time_step.exec(); - beam_volume_update.exec(); + column_volume_update.exec(); while (relaxation_time < advection_dt) { Real acoustic_dt = acoustic_time_step.exec(); @@ -153,6 +153,7 @@ int main(int ac, char *av[]) corrected_configuration.exec(); } TickCount t2 = TickCount::now(); + column_von_mises_stress.exec(); write_states.writeToFile(ite); write_displacement.writeToFile(ite); TickCount t3 = TickCount::now(); From 5d0e828d7d9a2ee05f7c171f3227986314f0205a Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 20:42:03 +0200 Subject: [PATCH 17/21] remove the neighbour searching in initialization and update --- .../continuum_integration.cpp | 45 ++++++++++++------- .../continuum_integration.h | 31 ++++++++----- .../continuum_integration.hpp | 2 +- .../oscillating_beam_ULGNOG.cpp | 9 ++-- .../test_3d_taylor_bar_UL/taylor_bar_UL.cpp | 6 ++- 5 files changed, 59 insertions(+), 34 deletions(-) diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 5921d97634..f36233a9ef 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -134,8 +134,8 @@ void StressDiffusion::interaction(size_t index_i, Real dt) stress_rate_3D_[index_i] = diffusion_stress_rate_; } //====================================================================================// -ShearStressRelaxationHourglassControl :: - ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi) +ShearStressRelaxationHourglassControl1stHalf :: + ShearStressRelaxationHourglassControl1stHalf(BaseInnerRelation &inner_relation, Real xi) : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), shear_stress_(particles_->registerStateVariable("ShearStress")), @@ -144,10 +144,7 @@ ShearStressRelaxationHourglassControl :: strain_tensor_(particles_->registerStateVariable("StrainTensor")), strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), - scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce", Real(1.0))), - acc_shear_(particles_->registerStateVariable("AccelerationByShear")), - acc_hourglass_(particles_->registerStateVariable("AccelerationHourglass")), - G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), xi_(xi) + scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce")), xi_(xi) { particles_->addVariableToSort("ShearStress"); particles_->addVariableToSort("ShearStressRate"); @@ -155,11 +152,9 @@ ShearStressRelaxationHourglassControl :: particles_->addVariableToSort("StrainTensor"); particles_->addVariableToSort("StrainTensorRate"); particles_->addVariableToSort("ScalePenaltyForce"); - particles_->addVariableToSort("AccelerationByShear"); - particles_->addVariableToSort("AccelerationHourglass"); } //====================================================================================// -void ShearStressRelaxationHourglassControl::initialization(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl1stHalf::interaction(size_t index_i, Real dt) { Matd velocity_gradient = Matd::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; @@ -176,15 +171,31 @@ void ShearStressRelaxationHourglassControl::initialization(size_t index_i, Real velocity_gradient_[index_i] = velocity_gradient; } //====================================================================================// -void ShearStressRelaxationHourglassControl::interaction(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl1stHalf::update(size_t index_i, Real dt) { + scale_penalty_force_[index_i] = xi_; shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt; } //====================================================================================// -void ShearStressRelaxationHourglassControl::update(size_t index_i, Real dt) +ShearStressRelaxationHourglassControl2ndHalf :: + ShearStressRelaxationHourglassControl2ndHalf(BaseInnerRelation &inner_relation) + : fluid_dynamics::BaseIntegration(inner_relation), + continuum_(DynamicCast(this, particles_->getBaseMaterial())), + shear_stress_(particles_->getVariableDataByName("ShearStress")), + velocity_gradient_(particles_->getVariableDataByName("VelocityGradient")), + acc_shear_(particles_->registerStateVariable("AccelerationByShear")), + acc_hourglass_(particles_->registerStateVariable("AccelerationHourglass")), + scale_penalty_force_(particles_->getVariableDataByName("ScalePenaltyForce")), + G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())) +{ + particles_->addVariableToSort("AccelerationByShear"); + particles_->addVariableToSort("AccelerationHourglass"); +} +//====================================================================================// +void ShearStressRelaxationHourglassControl2ndHalf::interaction(size_t index_i, Real dt) { Real rho_i = rho_[index_i]; Matd shear_stress_i = shear_stress_[index_i]; @@ -200,28 +211,28 @@ void ShearStressRelaxationHourglassControl::update(size_t index_i, Real dt) Vecd v_ij = vel_[index_i] - vel_[index_j]; Real r_ij = inner_neighborhood.r_ij_[n]; Vecd v_ij_correction = v_ij - 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_j]) * r_ij * e_ij; - acceleration_hourglass += 0.5 * xi_ * (scale_penalty_force_[index_i] + scale_penalty_force_[index_j]) * G_ * v_ij_correction * dW_ijV_j * dt / (rho_i * r_ij); + acceleration_hourglass += 0.5 * (scale_penalty_force_[index_i] + scale_penalty_force_[index_j]) * G_ * v_ij_correction * dW_ijV_j * dt / (rho_i * r_ij); } acc_hourglass_[index_i] += acceleration_hourglass; acc_shear_[index_i] = acceleration + acc_hourglass_[index_i]; } //====================================================================================// -ShearStressRelaxationHourglassControlJ2Plasticity :: - ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) - : ShearStressRelaxationHourglassControl(inner_relation, xi), +ShearStressRelaxationHourglassControl1stHalfJ2Plasticity :: + ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(BaseInnerRelation &inner_relation, Real xi) + : ShearStressRelaxationHourglassControl1stHalf(inner_relation, xi), J2_plasticity_(DynamicCast(this, particles_->getBaseMaterial())), hardening_factor_(particles_->registerStateVariable("HardeningFactor")) { particles_->addVariableToSort("HardeningFactor"); } //====================================================================================// -void ShearStressRelaxationHourglassControlJ2Plasticity::interaction(size_t index_i, Real dt) +void ShearStressRelaxationHourglassControl1stHalfJ2Plasticity::update(size_t index_i, Real dt) { shear_stress_rate_[index_i] = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]); Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate_[index_i] * dt; Real hardening_factor_increment = J2_plasticity_.HardeningFactorRate(shear_stress_try, hardening_factor_[index_i]); hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment; - scale_penalty_force_[index_i] = J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); + scale_penalty_force_[index_i] = xi_ * J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); shear_stress_[index_i] = J2_plasticity_.ReturnMappingShearStress(shear_stress_try, hardening_factor_[index_i]); Matd strain_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); strain_tensor_[index_i] += strain_rate * dt; diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index cc1464c839..78c2e4c219 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -214,30 +214,41 @@ class StressDiffusion : public BasePlasticIntegration Real smoothing_length_, sound_speed_; }; -class ShearStressRelaxationHourglassControl : public fluid_dynamics::BaseIntegration +class ShearStressRelaxationHourglassControl1stHalf : public fluid_dynamics::BaseIntegration { public: - explicit ShearStressRelaxationHourglassControl(BaseInnerRelation &inner_relation, Real xi_ = 4.0); - virtual ~ShearStressRelaxationHourglassControl(){}; - void initialization(size_t index_i, Real dt = 0.0); + explicit ShearStressRelaxationHourglassControl1stHalf(BaseInnerRelation &inner_relation, Real xi = 4.0); + virtual ~ShearStressRelaxationHourglassControl1stHalf(){}; void interaction(size_t index_i, Real dt = 0.0); void update(size_t index_i, Real dt = 0.0); protected: GeneralContinuum &continuum_; Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_; - Real *scale_penalty_force_; - Vecd *acc_shear_, *acc_hourglass_; - Real G_, xi_; + Real *scale_penalty_force_, xi_; }; -class ShearStressRelaxationHourglassControlJ2Plasticity : public ShearStressRelaxationHourglassControl +class ShearStressRelaxationHourglassControl2ndHalf : public fluid_dynamics::BaseIntegration { public: - explicit ShearStressRelaxationHourglassControlJ2Plasticity(BaseInnerRelation &inner_relation, Real xi_ = 0.2); - virtual ~ShearStressRelaxationHourglassControlJ2Plasticity(){}; + explicit ShearStressRelaxationHourglassControl2ndHalf(BaseInnerRelation &inner_relation); + virtual ~ShearStressRelaxationHourglassControl2ndHalf(){}; void interaction(size_t index_i, Real dt = 0.0); + protected: + GeneralContinuum &continuum_; + Matd *shear_stress_, *velocity_gradient_; + Vecd *acc_shear_, *acc_hourglass_; + Real *scale_penalty_force_, G_; +}; + +class ShearStressRelaxationHourglassControl1stHalfJ2Plasticity : public ShearStressRelaxationHourglassControl1stHalf +{ + public: + explicit ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(BaseInnerRelation &inner_relation, Real xi = 0.2); + virtual ~ShearStressRelaxationHourglassControl1stHalfJ2Plasticity(){}; + void update(size_t index_i, Real dt = 0.0); + protected: J2Plasticity &J2_plasticity_; Real *hardening_factor_; diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp index 87e50b6d69..ede9a43e39 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp @@ -11,7 +11,7 @@ namespace continuum_dynamics template BaseIntegration1stHalf::BaseIntegration1stHalf(BaseInnerRelation &inner_relation) : FluidDynamicsType(inner_relation), - acc_shear_(this->particles_->template getVariableDataByName("AccelerationByShear")) {} + acc_shear_(this->particles_->template registerStateVariable("AccelerationByShear")) {} //=================================================================================================// template void BaseIntegration1stHalf::update(size_t index_i, Real dt) diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp index 505320bc43..aca9b8b2e9 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp @@ -128,12 +128,12 @@ int main(int ac, char *av[]) // this section define all numerical methods will be used in this case //----------------------------------------------------------------------------- /** initial condition */ - InteractionWithUpdate correction_matrix(beam_body_inner); SimpleDynamics beam_initial_velocity(beam_body); - Dynamics1Level - beam_shear_acceleration_relaxation(beam_body_inner); + InteractionWithUpdate correction_matrix(beam_body_inner); Dynamics1Level beam_pressure_relaxation(beam_body_inner); Dynamics1Level beam_density_relaxation(beam_body_inner); + InteractionWithUpdate beam_shear_stress(beam_body_inner); + InteractionDynamics beam_shear_acceleration(beam_body_inner); SimpleDynamics beam_volume_update(beam_body); ReduceDynamics advection_time_step(beam_body, U_ref, 0.2); ReduceDynamics acoustic_time_step(beam_body, 0.4); @@ -187,7 +187,8 @@ int main(int ac, char *av[]) Real acoustic_dt = acoustic_time_step.exec(); beam_pressure_relaxation.exec(acoustic_dt); constraint_beam_base.exec(); - beam_shear_acceleration_relaxation.exec(acoustic_dt); + beam_shear_stress.exec(acoustic_dt); + beam_shear_acceleration.exec(acoustic_dt); beam_density_relaxation.exec(acoustic_dt); ite++; relaxation_time += acoustic_dt; diff --git a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp index 40ece24af5..943e53bc6f 100644 --- a/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp +++ b/tests/3d_examples/test_3d_taylor_bar_UL/taylor_bar_UL.cpp @@ -74,9 +74,10 @@ int main(int ac, char *av[]) SimpleDynamics initial_condition(column); SimpleDynamics wall_normal_direction(wall_boundary); InteractionWithUpdate corrected_configuration(column_inner); - Dynamics1Level column_shear_stress_relaxation(column_inner); Dynamics1Level column_pressure_relaxation(column_inner); Dynamics1Level column_density_relaxation(column_inner); + InteractionWithUpdate column_shear_stress(column_inner); + InteractionDynamics column_shear_acceleration(column_inner); SimpleDynamics column_volume_update(column); ReduceDynamics advection_time_step(column, U_max, 0.2); ReduceDynamics acoustic_time_step(column, 0.4); @@ -139,7 +140,8 @@ int main(int ac, char *av[]) column_wall_contact_force.exec(acoustic_dt); column_pressure_relaxation.exec(acoustic_dt); - column_shear_stress_relaxation.exec(acoustic_dt); + column_shear_stress.exec(acoustic_dt); + column_shear_acceleration.exec(acoustic_dt); column_density_relaxation.exec(acoustic_dt); ite++; From 36bc883c6aefed393f344338d982bfdcd4f57f67 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Tue, 17 Sep 2024 21:54:43 +0200 Subject: [PATCH 18/21] reduce the CI-test time for 2d_column_collapse --- .../continuum_integration.cpp | 2 +- .../continuum_integration.hpp | 2 - .../column_collapse.cpp | 4 +- .../GranularBody_TotalMechanicalEnergy.xml | 85 +++++++++---------- ...ody_TotalMechanicalEnergy_Run_0_result.xml | 4 +- ...dy_TotalMechanicalEnergy_Run_10_result.xml | 9 -- ...dy_TotalMechanicalEnergy_Run_11_result.xml | 9 -- ...dy_TotalMechanicalEnergy_Run_12_result.xml | 9 -- ...ody_TotalMechanicalEnergy_Run_1_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_2_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_3_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_4_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_5_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_6_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_7_result.xml | 4 +- ...ody_TotalMechanicalEnergy_Run_8_result.xml | 9 -- ...ody_TotalMechanicalEnergy_Run_9_result.xml | 9 -- ...Body_TotalMechanicalEnergy_dtwdistance.xml | 2 +- ...larBody_TotalMechanicalEnergy_runtimes.dat | 2 +- 19 files changed, 59 insertions(+), 115 deletions(-) delete mode 100644 tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml delete mode 100644 tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml delete mode 100644 tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml delete mode 100644 tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml delete mode 100644 tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index f36233a9ef..94b55ac354 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -173,9 +173,9 @@ void ShearStressRelaxationHourglassControl1stHalf::interaction(size_t index_i, R //====================================================================================// void ShearStressRelaxationHourglassControl1stHalf::update(size_t index_i, Real dt) { - scale_penalty_force_[index_i] = xi_; shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; + scale_penalty_force_[index_i] = xi_; strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt; } diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp index ede9a43e39..511c526c74 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.hpp @@ -1,7 +1,5 @@ #pragma once - #include "continuum_integration.h" - #include "base_particles.hpp" namespace SPH { diff --git a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp index 65f74ec306..a99589fc47 100644 --- a/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp +++ b/tests/2d_examples/test_2d_column_collapse/column_collapse.cpp @@ -134,8 +134,8 @@ int main(int ac, char *av[]) int screen_output_interval = 500; int observation_sample_interval = screen_output_interval * 2; int restart_output_interval = screen_output_interval * 10; - Real End_Time = 1.0; /**< End time. */ - Real D_Time = End_Time / 50; /**< Time stamps for output of body states. */ + Real End_Time = 0.8; /**< End time. */ + Real D_Time = End_Time / 40; /**< Time stamps for output of body states. */ Real Dt = 0.1 * D_Time; //---------------------------------------------------------------------- // Statistics for CPU time diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml index df2b1d25ba..461e092951 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy.xml @@ -1,51 +1,42 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml index 296f1f7722..b376707536 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_0_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml deleted file mode 100644 index 49c242b5e7..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_10_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml deleted file mode 100644 index fac5f2d438..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_11_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml deleted file mode 100644 index 8a1dc7daa1..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_12_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml index 9cacfb4cc1..f63c0bb733 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_1_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml index 20390edd95..a85d9d3ba4 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_2_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml index 9e018fae74..0e388bbb73 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_3_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml index fe62912d6f..e0347772fe 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_4_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml index ffb3d2a9ac..238c385f80 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_5_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml index fd67bd9b30..9e22ce14bf 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_6_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml index f71552ddb8..4facd994cc 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_7_result.xml @@ -1,9 +1,9 @@ - + - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml deleted file mode 100644 index 8297225a57..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_8_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml deleted file mode 100644 index f02d5b8586..0000000000 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_Run_9_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml index f3f38a290a..9f850afadb 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat index 19c6fe5aad..ffe6691bfd 100644 --- a/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat +++ b/tests/2d_examples/test_2d_column_collapse/regression_test_tool/GranularBody_TotalMechanicalEnergy_runtimes.dat @@ -1,3 +1,3 @@ true -13 +8 4 \ No newline at end of file From ab37a4c249edbc6e51a67acf2d32449e6e1b9d1a Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Wed, 18 Sep 2024 12:24:28 +0200 Subject: [PATCH 19/21] delete unnecassary registered variables --- .../continuum_integration.cpp | 19 +++++++------------ .../continuum_integration.h | 2 +- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 94b55ac354..9526869a24 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -139,18 +139,14 @@ ShearStressRelaxationHourglassControl1stHalf :: : fluid_dynamics::BaseIntegration(inner_relation), continuum_(DynamicCast(this, particles_->getBaseMaterial())), shear_stress_(particles_->registerStateVariable("ShearStress")), - shear_stress_rate_(particles_->registerStateVariable("ShearStressRate")), velocity_gradient_(particles_->registerStateVariable("VelocityGradient")), strain_tensor_(particles_->registerStateVariable("StrainTensor")), - strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), scale_penalty_force_(particles_->registerStateVariable("ScalePenaltyForce")), xi_(xi) { particles_->addVariableToSort("ShearStress"); - particles_->addVariableToSort("ShearStressRate"); particles_->addVariableToSort("VelocityGradient"); particles_->addVariableToSort("StrainTensor"); - particles_->addVariableToSort("StrainTensorRate"); particles_->addVariableToSort("ScalePenaltyForce"); } //====================================================================================// @@ -164,8 +160,7 @@ void ShearStressRelaxationHourglassControl1stHalf::interaction(size_t index_i, R size_t index_j = inner_neighborhood.j_[n]; Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; Vecd &e_ij = inner_neighborhood.e_ij_[n]; - Matd velocity_gradient_ij; - velocity_gradient_ij = -(vel_i - vel_[index_j]) * (B_[index_i] * e_ij * dW_ijV_j).transpose(); + Matd velocity_gradient_ij = -(vel_i - vel_[index_j]) * (B_[index_i] * e_ij * dW_ijV_j).transpose(); velocity_gradient += velocity_gradient_ij; } velocity_gradient_[index_i] = velocity_gradient; @@ -173,11 +168,11 @@ void ShearStressRelaxationHourglassControl1stHalf::interaction(size_t index_i, R //====================================================================================// void ShearStressRelaxationHourglassControl1stHalf::update(size_t index_i, Real dt) { - shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); - shear_stress_[index_i] += shear_stress_rate_[index_i] * dt; + Matd shear_stress_rate = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); + shear_stress_[index_i] += shear_stress_rate * dt; scale_penalty_force_[index_i] = xi_; - strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); - strain_tensor_[index_i] += strain_tensor_rate_[index_i] * dt; + Matd strain_tensor_rate = 0.5 * (velocity_gradient_[index_i] + velocity_gradient_[index_i].transpose()); + strain_tensor_[index_i] += strain_tensor_rate * dt; } //====================================================================================// ShearStressRelaxationHourglassControl2ndHalf :: @@ -228,8 +223,8 @@ ShearStressRelaxationHourglassControl1stHalfJ2Plasticity :: //====================================================================================// void ShearStressRelaxationHourglassControl1stHalfJ2Plasticity::update(size_t index_i, Real dt) { - shear_stress_rate_[index_i] = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]); - Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate_[index_i] * dt; + Matd shear_stress_rate = J2_plasticity_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i], hardening_factor_[index_i]); + Matd shear_stress_try = shear_stress_[index_i] + shear_stress_rate * dt; Real hardening_factor_increment = J2_plasticity_.HardeningFactorRate(shear_stress_try, hardening_factor_[index_i]); hardening_factor_[index_i] += sqrt(2.0 / 3.0) * hardening_factor_increment; scale_penalty_force_[index_i] = xi_ * J2_plasticity_.ScalePenaltyForce(shear_stress_try, hardening_factor_[index_i]); diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 78c2e4c219..76589e18a5 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -224,7 +224,7 @@ class ShearStressRelaxationHourglassControl1stHalf : public fluid_dynamics::Base protected: GeneralContinuum &continuum_; - Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_; + Matd *shear_stress_, *velocity_gradient_, *strain_tensor_, *B_; Real *scale_penalty_force_, xi_; }; From 29049ca2572b7a15a2d2b3451da94ba102925295 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Wed, 18 Sep 2024 12:26:34 +0200 Subject: [PATCH 20/21] samll changes --- .../continuum_integration.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 9526869a24..4fa8ca0b54 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -115,8 +115,8 @@ void StressDiffusion::interaction(size_t index_i, Real dt) Vecd acc_prior_i = force_prior_[index_i] / mass_[index_i]; Real gravity = abs(acc_prior_i(1, 0)); Real density = plastic_continuum_.getDensity(); - Mat3d diffusion_stress_rate_ = Mat3d::Zero(); - Mat3d diffusion_stress_ = Mat3d::Zero(); + Mat3d diffusion_stress_rate = Mat3d::Zero(); + Mat3d diffusion_stress = Mat3d::Zero(); Neighborhood &inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) { @@ -124,14 +124,14 @@ void StressDiffusion::interaction(size_t index_i, Real dt) Real r_ij = inner_neighborhood.r_ij_[n]; Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; Real y_ij = pos_[index_i](1, 0) - pos_[index_j](1, 0); - diffusion_stress_ = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; - diffusion_stress_(0, 0) -= (1 - sin(phi_)) * density * gravity * y_ij; - diffusion_stress_(1, 1) -= density * gravity * y_ij; - diffusion_stress_(2, 2) -= (1 - sin(phi_)) * density * gravity * y_ij; - diffusion_stress_rate_ += 2 * zeta_ * smoothing_length_ * sound_speed_ * - diffusion_stress_ * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); + diffusion_stress = stress_tensor_3D_[index_i] - stress_tensor_3D_[index_j]; + diffusion_stress(0, 0) -= (1 - sin(phi_)) * density * gravity * y_ij; + diffusion_stress(1, 1) -= density * gravity * y_ij; + diffusion_stress(2, 2) -= (1 - sin(phi_)) * density * gravity * y_ij; + diffusion_stress_rate += 2 * zeta_ * smoothing_length_ * sound_speed_ * + diffusion_stress * r_ij * dW_ijV_j / (r_ij * r_ij + 0.01 * smoothing_length_); } - stress_rate_3D_[index_i] = diffusion_stress_rate_; + stress_rate_3D_[index_i] = diffusion_stress_rate; } //====================================================================================// ShearStressRelaxationHourglassControl1stHalf :: From c302adcc296aaec8ba9b0f4b70841a7d01f31c97 Mon Sep 17 00:00:00 2001 From: Shuaihao-Zhang <1477856817@qq.com> Date: Wed, 18 Sep 2024 16:39:16 +0200 Subject: [PATCH 21/21] delete the SPH-ENOG method since the present SPH-GNOG is better --- .../continuum_integration.cpp | 74 ----- .../continuum_integration.h | 28 -- .../oscillating_beam_UL.cpp | 28 +- .../BeamBody_TotalKineticEnergy.xml | 0 ...amBody_TotalKineticEnergy_Run_0_result.xml | 2 +- ...amBody_TotalKineticEnergy_Run_1_result.xml | 0 ...amBody_TotalKineticEnergy_Run_2_result.xml | 0 ...amBody_TotalKineticEnergy_Run_3_result.xml | 2 +- ...amBody_TotalKineticEnergy_Run_4_result.xml | 0 ...amBody_TotalKineticEnergy_Run_5_result.xml | 2 +- ...eamBody_TotalKineticEnergy_dtwdistance.xml | 2 +- .../CMakeLists.txt | 21 -- .../oscillating_beam_ULGNOG.cpp | 231 ---------------- ...amBody_TotalKineticEnergy_Run_0_result.xml | 9 - ...amBody_TotalKineticEnergy_Run_3_result.xml | 9 - ...amBody_TotalKineticEnergy_Run_5_result.xml | 9 - ...eamBody_TotalKineticEnergy_dtwdistance.xml | 4 - .../BeamBody_TotalKineticEnergy_runtimes.dat | 3 - .../regression_test_tool.py | 37 --- .../CMakeLists.txt | 21 -- .../oscillating_plate_UL.cpp | 259 ------------------ ...teBody_TotalKineticEnergy_Run_0_result.xml | 9 - ...teBody_TotalKineticEnergy_Run_3_result.xml | 9 - ...teBody_TotalKineticEnergy_Run_5_result.xml | 9 - ...ateBody_TotalKineticEnergy_dtwdistance.xml | 4 - .../PlateBody_TotalKineticEnergy_runtimes.dat | 3 - .../regression_test_tool.py | 37 --- 27 files changed, 16 insertions(+), 796 deletions(-) rename tests/2d_examples/{test_2d_oscillating_beam_ULGNOG => test_2d_oscillating_beam_UL}/regression_test_tool/BeamBody_TotalKineticEnergy.xml (100%) rename tests/2d_examples/{test_2d_oscillating_beam_ULGNOG => test_2d_oscillating_beam_UL}/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml (100%) rename tests/2d_examples/{test_2d_oscillating_beam_ULGNOG => test_2d_oscillating_beam_UL}/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml (100%) rename tests/2d_examples/{test_2d_oscillating_beam_ULGNOG => test_2d_oscillating_beam_UL}/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml (100%) delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat delete mode 100644 tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat delete mode 100644 tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp index 4fa8ca0b54..aac2151fde 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.cpp @@ -29,80 +29,6 @@ Real AcousticTimeStep::outputResult(Real reduced_value) { return acousticCFL_ * smoothing_length_min_ / (reduced_value + TinyReal); } -//=================================================================================================// -ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation &inner_relation) - : fluid_dynamics::BaseIntegration(inner_relation), - continuum_(DynamicCast(this, particles_->getBaseMaterial())), - G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), - smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), - acc_shear_(particles_->registerStateVariable("AccelerationByShear")) -{ - particles_->addVariableToSort("AccelerationByShear"); -} -//=================================================================================================// -void ShearAccelerationRelaxation::interaction(size_t index_i, Real dt) -{ - Real rho_i = rho_[index_i]; - Vecd acceleration = Vecd::Zero(); - Neighborhood &inner_neighborhood = inner_configuration_[index_i]; - for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) - { - size_t index_j = inner_neighborhood.j_[n]; - Real r_ij = inner_neighborhood.r_ij_[n]; - Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j]; - Vecd &e_ij = inner_neighborhood.e_ij_[n]; - Real eta_ij = 2 * (0.7 * (Real)Dimensions + 2.1) * (vel_[index_i] - vel_[index_j]).dot(e_ij) / (r_ij + TinyReal); - acceleration += eta_ij * dW_ijV_j * e_ij; - } - acc_shear_[index_i] += G_ * acceleration * dt / rho_i; -} -//=================================================================================================// -ShearStressRelaxation::ShearStressRelaxation(BaseInnerRelation &inner_relation) - : fluid_dynamics::BaseIntegration(inner_relation), - continuum_(DynamicCast(this, particles_->getBaseMaterial())), - shear_stress_(particles_->registerStateVariable("ShearStress")), - shear_stress_rate_(particles_->registerStateVariable("ShearStressRate")), - velocity_gradient_(particles_->registerStateVariable("VelocityGradient")), - strain_tensor_(particles_->registerStateVariable("StrainTensor")), - strain_tensor_rate_(particles_->registerStateVariable("StrainTensorRate")), - B_(particles_->getVariableDataByName("LinearGradientCorrectionMatrix")), - Vol_(particles_->getVariableDataByName("VolumetricMeasure")) -{ - particles_->addVariableToSort("ShearStress"); - particles_->addVariableToSort("ShearStressRate"); - particles_->addVariableToSort("VelocityGradient"); - particles_->addVariableToSort("StrainTensor"); - particles_->addVariableToSort("StrainTensorRate"); -} -//====================================================================================// -void ShearStressRelaxation::initialization(size_t index_i, Real dt) -{ - strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; - shear_stress_[index_i] += shear_stress_rate_[index_i] * 0.5 * dt; -} -//====================================================================================// -void ShearStressRelaxation::interaction(size_t index_i, Real dt) -{ - Matd velocity_gradient = Matd::Zero(); - Neighborhood &inner_neighborhood = inner_configuration_[index_i]; - for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) - { - size_t index_j = inner_neighborhood.j_[n]; - Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_i]; - Vecd &e_ij = inner_neighborhood.e_ij_[n]; - Vecd v_ij = vel_[index_i] - vel_[index_j]; - velocity_gradient -= v_ij * (B_[index_i] * e_ij * dW_ijV_j).transpose(); - } - velocity_gradient_[index_i] = velocity_gradient; - strain_tensor_rate_[index_i] = 0.5 * (velocity_gradient + velocity_gradient.transpose()); - strain_tensor_[index_i] += strain_tensor_rate_[index_i] * 0.5 * dt; -} -//====================================================================================// -void ShearStressRelaxation::update(size_t index_i, Real dt) -{ - shear_stress_rate_[index_i] = continuum_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); - shear_stress_[index_i] += shear_stress_rate_[index_i] * dt * 0.5; -} //====================================================================================// StressDiffusion::StressDiffusion(BaseInnerRelation &inner_relation) : BasePlasticIntegration(inner_relation), diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 76589e18a5..85362e5bbc 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -79,34 +79,6 @@ class BaseIntegration1stHalf : public FluidDynamicsType using Integration1stHalf = BaseIntegration1stHalf; using Integration1stHalfRiemann = BaseIntegration1stHalf; -class ShearAccelerationRelaxation : public fluid_dynamics::BaseIntegration -{ - public: - explicit ShearAccelerationRelaxation(BaseInnerRelation &inner_relation); - virtual ~ShearAccelerationRelaxation(){}; - void interaction(size_t index_i, Real dt = 0.0); - - protected: - GeneralContinuum &continuum_; - Real G_, smoothing_length_; - Vecd *acc_shear_; -}; - -class ShearStressRelaxation : public fluid_dynamics::BaseIntegration -{ - public: - explicit ShearStressRelaxation(BaseInnerRelation &inner_relation); - virtual ~ShearStressRelaxation(){}; - void initialization(size_t index_i, Real dt = 0.0); - void interaction(size_t index_i, Real dt = 0.0); - void update(size_t index_i, Real dt = 0.0); - - protected: - GeneralContinuum &continuum_; - Matd *shear_stress_, *shear_stress_rate_, *velocity_gradient_, *strain_tensor_, *strain_tensor_rate_, *B_; - Real *Vol_; -}; - template class BasePlasticIntegration : public fluid_dynamics::BaseIntegration { diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index 466b2bb589..aca9b8b2e9 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -1,8 +1,9 @@ /* ----------------------------------------------------------------------------* - * SPHinXsys: 2D oscillation beam--updated Lagrangian method * + * SPHinXsys: 2D oscillation beam * * ----------------------------------------------------------------------------* * This is the one of the basic test cases for understanding SPH method for * - * solid simulation based on updated Lagrangian method * + * solid simulation based on updated Lagrangian method. * + * A generalized hourglass control method is used here. * * In this case, the constraint of the beam is implemented with * * internal constrained subregion. * * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu * @@ -107,7 +108,6 @@ int main(int ac, char *av[]) // Creating body, materials and particles. //---------------------------------------------------------------------- RealBody beam_body(sph_system, makeShared("BeamBody")); - beam_body.sph_adaptation_->resetKernel(); beam_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); beam_body.generateParticles(); @@ -128,17 +128,15 @@ int main(int ac, char *av[]) // this section define all numerical methods will be used in this case //----------------------------------------------------------------------------- /** initial condition */ - InteractionWithUpdate correction_matrix(beam_body_inner); SimpleDynamics beam_initial_velocity(beam_body); - - InteractionDynamics beam_shear_acceleration(beam_body_inner); - Dynamics1Level beam_shear_stress_relaxation(beam_body_inner); + InteractionWithUpdate correction_matrix(beam_body_inner); Dynamics1Level beam_pressure_relaxation(beam_body_inner); Dynamics1Level beam_density_relaxation(beam_body_inner); + InteractionWithUpdate beam_shear_stress(beam_body_inner); + InteractionDynamics beam_shear_acceleration(beam_body_inner); SimpleDynamics beam_volume_update(beam_body); - - ReduceDynamics fluid_advection_time_step(beam_body, U_ref, 0.2); - ReduceDynamics fluid_acoustic_time_step(beam_body, 0.4); + ReduceDynamics advection_time_step(beam_body, U_ref, 0.2); + ReduceDynamics acoustic_time_step(beam_body, 0.4); // clamping a solid body part. BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); SimpleDynamics constraint_beam_base(beam_base); @@ -182,17 +180,16 @@ int main(int ac, char *av[]) while (integration_time < D_Time) { Real relaxation_time = 0.0; - Real advection_dt = fluid_advection_time_step.exec(); + Real advection_dt = advection_time_step.exec(); beam_volume_update.exec(); - while (relaxation_time < advection_dt) { - Real acoustic_dt = fluid_acoustic_time_step.exec(); - beam_shear_stress_relaxation.exec(acoustic_dt); + Real acoustic_dt = acoustic_time_step.exec(); beam_pressure_relaxation.exec(acoustic_dt); constraint_beam_base.exec(); - beam_density_relaxation.exec(acoustic_dt); + beam_shear_stress.exec(acoustic_dt); beam_shear_acceleration.exec(acoustic_dt); + beam_density_relaxation.exec(acoustic_dt); ite++; relaxation_time += acoustic_dt; integration_time += acoustic_dt; @@ -230,6 +227,5 @@ int main(int ac, char *av[]) { write_beam_kinetic_energy.testResult(); } - return 0; } \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy.xml similarity index 100% rename from tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy.xml rename to tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy.xml diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml index 1f10003773..7a8cacca58 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml similarity index 100% rename from tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml rename to tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_1_result.xml diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml similarity index 100% rename from tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml rename to tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_2_result.xml diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml index bb2d859705..4f75fb4c37 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml similarity index 100% rename from tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml rename to tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_4_result.xml diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml index 5789bcd8a4..0313a9973c 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml @@ -4,6 +4,6 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml index 8be7b3f132..e59f0fa6a5 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml @@ -1,4 +1,4 @@ - + diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt deleted file mode 100644 index 36b6b1e42d..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) -PROJECT("${CURRENT_FOLDER}") - -SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) -SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") -SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") -SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") - -file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ - DESTINATION ${BUILD_INPUT_PATH}) - -add_executable(${PROJECT_NAME}) -aux_source_directory(. DIR_SRCS) -target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) -target_link_libraries(${PROJECT_NAME} sphinxsys_2d) -set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") - -add_test(NAME ${PROJECT_NAME} - COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp deleted file mode 100644 index aca9b8b2e9..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/oscillating_beam_ULGNOG.cpp +++ /dev/null @@ -1,231 +0,0 @@ -/* ----------------------------------------------------------------------------* - * SPHinXsys: 2D oscillation beam * - * ----------------------------------------------------------------------------* - * This is the one of the basic test cases for understanding SPH method for * - * solid simulation based on updated Lagrangian method. * - * A generalized hourglass control method is used here. * - * In this case, the constraint of the beam is implemented with * - * internal constrained subregion. * - * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu * - * ----------------------------------------------------------------------------*/ -#include "sphinxsys.h" -using namespace SPH; -//------------------------------------------------------------------------------ -// global parameters for the case -//------------------------------------------------------------------------------ -Real PL = 0.2; // beam length -Real PH = 0.02; // for thick plate -Real SL = 0.06; // depth of the insert -Real resolution_ref = PH / 10; -Real BW = resolution_ref * 4; // boundary width, at least three particles -/** Domain bounds of the system. */ -BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), - Vec2d(PL + 3.0 * BW, PL / 2.0)); -//---------------------------------------------------------------------- -// Material properties of the fluid. -//---------------------------------------------------------------------- -Real rho0_s = 1.0e3; // reference density -Real Youngs_modulus = 2.0e6; // reference Youngs modulus -Real poisson = 0.3975; // Poisson ratio -Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); -//---------------------------------------------------------------------- -// Parameters for initial condition on velocity -//---------------------------------------------------------------------- -Real kl = 1.875; -Real M = sin(kl) + sinh(kl); -Real N = cos(kl) + cosh(kl); -Real Q = 2.0 * (cos(kl) * sinh(kl) - sin(kl) * cosh(kl)); -Real vf = 0.05; -Real R = PL / (0.5 * Pi); -Real U_ref = vf * c0 * (M * (cos(kl) - cosh(kl)) - N * (sin(kl) - sinh(kl))) / Q; -//---------------------------------------------------------------------- -// Geometric shapes used in the system. -//---------------------------------------------------------------------- -// a beam base shape -std::vector beam_base_shape{ - Vecd(-SL - BW, -PH / 2 - BW), Vecd(-SL - BW, PH / 2 + BW), Vecd(0.0, PH / 2 + BW), - Vecd(0.0, -PH / 2 - BW), Vecd(-SL - BW, -PH / 2 - BW)}; -// a beam shape -std::vector beam_shape{ - Vecd(-SL, -PH / 2), Vecd(-SL, PH / 2), Vecd(PL, PH / 2), Vecd(PL, -PH / 2), Vecd(-SL, -PH / 2)}; -// Beam observer location -StdVec observation_location = {Vecd(PL, 0.0)}; -//---------------------------------------------------------------------- -// Define the beam body -//---------------------------------------------------------------------- -class Beam : public MultiPolygonShape -{ - public: - explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - multi_polygon_.addAPolygon(beam_base_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); - } -}; -//---------------------------------------------------------------------- -// application dependent initial condition -//---------------------------------------------------------------------- -class BeamInitialCondition - : public fluid_dynamics::FluidInitialCondition -{ - public: - explicit BeamInitialCondition(RealBody &beam_column) - : fluid_dynamics::FluidInitialCondition(beam_column){}; - - protected: - void update(size_t index_i, Real dt) - { - /** initial velocity profile */ - Real x = pos_[index_i][0] / PL; - if (x > 0.0) - { - vel_[index_i][1] = vf * c0 * - (M * (cos(kl * x) - cosh(kl * x)) - N * (sin(kl * x) - sinh(kl * x))) / Q; - } - }; -}; -//---------------------------------------------------------------------- -// define the beam base which will be constrained. -//---------------------------------------------------------------------- -MultiPolygon createBeamConstrainShape() -{ - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_base_shape, ShapeBooleanOps::add); - multi_polygon.addAPolygon(beam_shape, ShapeBooleanOps::sub); - return multi_polygon; -}; -//------------------------------------------------------------------------------ -// the main program -//------------------------------------------------------------------------------ -int main(int ac, char *av[]) -{ - //---------------------------------------------------------------------- - // Build up the environment of a SPHSystem with global controls. - //---------------------------------------------------------------------- - SPHSystem sph_system(system_domain_bounds, resolution_ref); - sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); - //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - RealBody beam_body(sph_system, makeShared("BeamBody")); - beam_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); - beam_body.generateParticles(); - - ObserverBody beam_observer(sph_system, "BeamObserver"); - beam_observer.sph_adaptation_->resetAdaptationRatios(1.15, 2.0); - beam_observer.generateParticles(observation_location); - //---------------------------------------------------------------------- - // Define body relation map. - // The contact map gives the topological connections between the bodies. - // Basically the the range of bodies to build neighbor particle lists. - // Generally, we first define all the inner relations, then the contact relations. - // At last, we define the complex relaxations by combining previous defined - // inner and contact relations. - //---------------------------------------------------------------------- - ContactRelation beam_observer_contact(beam_observer, {&beam_body}); - InnerRelation beam_body_inner(beam_body); - //----------------------------------------------------------------------------- - // this section define all numerical methods will be used in this case - //----------------------------------------------------------------------------- - /** initial condition */ - SimpleDynamics beam_initial_velocity(beam_body); - InteractionWithUpdate correction_matrix(beam_body_inner); - Dynamics1Level beam_pressure_relaxation(beam_body_inner); - Dynamics1Level beam_density_relaxation(beam_body_inner); - InteractionWithUpdate beam_shear_stress(beam_body_inner); - InteractionDynamics beam_shear_acceleration(beam_body_inner); - SimpleDynamics beam_volume_update(beam_body); - ReduceDynamics advection_time_step(beam_body, U_ref, 0.2); - ReduceDynamics acoustic_time_step(beam_body, 0.4); - // clamping a solid body part. - BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); - //----------------------------------------------------------------------------- - // outputs - //----------------------------------------------------------------------------- - BodyStatesRecordingToVtp write_beam_states(beam_body); - write_beam_states.addToWrite(beam_body, "Density"); - write_beam_states.addToWrite(beam_body, "Pressure"); - SimpleDynamics beam_von_mises_stress(beam_body); - write_beam_states.addToWrite(beam_body, "VonMisesStress"); - ObservedQuantityRecording write_beam_tip_displacement("Position", beam_observer_contact); - RegressionTestDynamicTimeWarping> write_beam_kinetic_energy(beam_body); - //---------------------------------------------------------------------- - // Setup computing and initial conditions. - //---------------------------------------------------------------------- - sph_system.initializeSystemCellLinkedLists(); - sph_system.initializeSystemConfigurations(); - beam_initial_velocity.exec(); - correction_matrix.exec(); - //---------------------------------------------------------------------- - // Setup computing time-step controls. - //---------------------------------------------------------------------- - Real &physical_time = *sph_system.getSystemVariableDataByName("PhysicalTime"); - int ite = 0; - Real T0 = 1.0; - Real End_Time = T0; - Real D_Time = End_Time / 100; /**< Time period for data observing */ - TickCount t1 = TickCount::now(); - TimeInterval interval; - //----------------------------------------------------------------------------- - // from here the time stepping begins - //----------------------------------------------------------------------------- - write_beam_states.writeToFile(0); - write_beam_tip_displacement.writeToFile(0); - write_beam_kinetic_energy.writeToFile(0); - // computation loop starts - while (physical_time < End_Time) - { - Real integration_time = 0.0; - while (integration_time < D_Time) - { - Real relaxation_time = 0.0; - Real advection_dt = advection_time_step.exec(); - beam_volume_update.exec(); - while (relaxation_time < advection_dt) - { - Real acoustic_dt = acoustic_time_step.exec(); - beam_pressure_relaxation.exec(acoustic_dt); - constraint_beam_base.exec(); - beam_shear_stress.exec(acoustic_dt); - beam_shear_acceleration.exec(acoustic_dt); - beam_density_relaxation.exec(acoustic_dt); - ite++; - relaxation_time += acoustic_dt; - integration_time += acoustic_dt; - physical_time += acoustic_dt; - if (ite % 500 == 0) - { - std::cout << "N=" << ite << " Time: " - << physical_time << " advection_dt: " - << advection_dt << " acoustic_dt: " - << acoustic_dt << "\n"; - } - } - beam_body.updateCellLinkedList(); - beam_body_inner.updateConfiguration(); - correction_matrix.exec(); - } - beam_von_mises_stress.exec(); - write_beam_tip_displacement.writeToFile(ite); - write_beam_kinetic_energy.writeToFile(ite); - TickCount t2 = TickCount::now(); - write_beam_states.writeToFile(); - TickCount t3 = TickCount::now(); - interval += t3 - t2; - } - TickCount t4 = TickCount::now(); - TimeInterval tt; - tt = t4 - t1 - interval; - std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; - - if (sph_system.GenerateRegressionData()) - { - write_beam_kinetic_energy.generateDataBase(1.0e-3); - } - else - { - write_beam_kinetic_energy.testResult(); - } - return 0; -} \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml deleted file mode 100644 index 7a8cacca58..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_0_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml deleted file mode 100644 index 4f75fb4c37..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_3_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml deleted file mode 100644 index 0313a9973c..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_Run_5_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml deleted file mode 100644 index e59f0fa6a5..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_dtwdistance.xml +++ /dev/null @@ -1,4 +0,0 @@ - - - - diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat deleted file mode 100644 index de07de18db..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/BeamBody_TotalKineticEnergy_runtimes.dat +++ /dev/null @@ -1,3 +0,0 @@ -true -6 -4 \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py b/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py deleted file mode 100644 index 4ad6ca7edc..0000000000 --- a/tests/2d_examples/test_2d_oscillating_beam_ULGNOG/regression_test_tool/regression_test_tool.py +++ /dev/null @@ -1,37 +0,0 @@ -# !/usr/bin/env python3 -import os -import sys - -path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') -sys.path.append(path) -from regression_test_base_tool import SphinxsysRegressionTest - -""" -case name: test_2d_oscillating_beam_ULGNOG -""" - -case_name = "test_2d_oscillating_beam_ULGNOG" -body_name = "BeamBody" -parameter_name = "TotalKineticEnergy" - -number_of_run_times = 0 -converged = 0 -sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) - - -while True: - print("Now start a new run......") - sphinxsys.run_case() - number_of_run_times += 1 - converged = sphinxsys.read_dat_file() - print("Please note: This is the", number_of_run_times, "run!") - if number_of_run_times <= 200: - if converged == "true": - print("The tested parameters of all variables are converged, and the run will stop here!") - break - elif converged != "true": - print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") - continue - else: - print("It's too many runs but still not converged, please try again!") - break diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt b/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt deleted file mode 100644 index 23ebe87350..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) -PROJECT("${CURRENT_FOLDER}") - -SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) -SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") -SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") -SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") - -file(MAKE_DIRECTORY ${BUILD_INPUT_PATH}) -file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/ - DESTINATION ${BUILD_INPUT_PATH}) - -add_executable(${PROJECT_NAME}) -aux_source_directory(. DIR_SRCS) -target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) -target_link_libraries(${PROJECT_NAME} sphinxsys_3d) -set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") - -add_test(NAME ${PROJECT_NAME} - COMMAND ${PROJECT_NAME} --state_recording=${TEST_STATE_RECORDING} - WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp deleted file mode 100644 index e0b6da44d0..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp +++ /dev/null @@ -1,259 +0,0 @@ -/** - * @file oscillating_plate_UL.cpp - * @brief This is the test case for the hourglass manuscript. - * @details We consider vibration deformation of a square plate under initial vertical velocity field. - * @author Shuaihao Zhang, Dong Wu and Xiangyu Hu - */ -#include "sphinxsys.h" -using namespace SPH; -//---------------------------------------------------------------------- -// Basic geometry parameters and numerical setup. -//---------------------------------------------------------------------- -Real PL = 0.4; /** Length of the square plate. */ -Real PH = 0.4; /** Width of the square plate. */ -Real PT = 0.01; /** Thickness of the square plate. */ -int particle_number = 3; /** Particle number in the direction of the thickness. */ -Real particle_spacing_ref = PT / (Real)particle_number; /** Initial reference particle spacing. */ -int particle_number_PL = PL / particle_spacing_ref; -int particle_number_PH = PH / particle_spacing_ref; -int BWD = 1; -Real BW = particle_spacing_ref * (Real)BWD; /** Boundary width, determined by specific layer of boundary particles. */ -BoundingBox system_domain_bounds(Vec3d(-BW, -BW, -PT / 2), Vec3d(PL + BW, PH + BW, PT / 2)); -StdVec observation_location = {Vecd(0.5 * PL, 0.5 * PH, 0.0), Vecd(-BW, -BW, 0.0)}; -//---------------------------------------------------------------------- -// Material parameters. -//---------------------------------------------------------------------- -Real rho0_s = 1000.0; /** Normalized density. */ -Real Youngs_modulus = 100.0e6; /** Normalized Youngs Modulus. */ -Real poisson = 0.3; /** Poisson ratio. */ -Real c0 = sqrt(Youngs_modulus / (3 * (1 - 2 * poisson) * rho0_s)); - -Real governing_vibration_integer_x = 2.0; -Real governing_vibration_integer_y = 2.0; -Real U_ref = 1.0; // Maximum velocity - -namespace SPH -{ -//---------------------------------------------------------------------- -// Complex shape for wall boundary, note that no partial overlap is allowed -// for the shapes in a complex shape. -//---------------------------------------------------------------------- -class Plate; -template <> -class ParticleGenerator : public ParticleGenerator -{ - public: - explicit ParticleGenerator(SPHBody &sph_body, BaseParticles &base_particles) - : ParticleGenerator(sph_body, base_particles){}; - virtual void prepareGeometricData() override - { - for (int k = 0; k < particle_number; k++) - { - for (int i = 0; i < particle_number_PL + 2 * BWD; i++) - { - for (int j = 0; j < particle_number_PH + 2 * BWD; j++) - { - Real x = particle_spacing_ref * i - BW + particle_spacing_ref * 0.5; - Real y = particle_spacing_ref * j - BW + particle_spacing_ref * 0.5; - Real z = particle_spacing_ref * (k - ((particle_number - 1.0) / 2.0)); - addPositionAndVolumetricMeasure( - Vecd(x, y, z), particle_spacing_ref * particle_spacing_ref * particle_spacing_ref); - } - } - } - } -}; - -/** Define the boundary geometry. */ -class BoundaryGeometry : public BodyPartByParticle -{ - public: - BoundaryGeometry(SPHBody &body, const std::string &body_part_name) - : BodyPartByParticle(body, body_part_name) - { - TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometry::tagManually, this, _1); - tagParticles(tagging_particle_method); - }; - virtual ~BoundaryGeometry(){}; - - private: - void tagManually(size_t index_i) - { - if ((base_particles_.ParticlePositions()[index_i][2] < 0.5 * particle_spacing_ref) && - (base_particles_.ParticlePositions()[index_i][2] > -0.5 * particle_spacing_ref) && - (base_particles_.ParticlePositions()[index_i][0] < 0.0 || - base_particles_.ParticlePositions()[index_i][0] > PL || - base_particles_.ParticlePositions()[index_i][1] < 0.0 || - base_particles_.ParticlePositions()[index_i][1] > PH)) - { - body_part_particles_.push_back(index_i); - } - }; -}; - -/** Define the initial condition. */ -class BeamInitialCondition - : public fluid_dynamics::FluidInitialCondition -{ - public: - explicit BeamInitialCondition(SPHBody &sph_body) - : fluid_dynamics::FluidInitialCondition(sph_body){}; - - void update(size_t index_i, Real dt) - { - /** initial velocity profile */ - vel_[index_i][2] = sin(governing_vibration_integer_x * Pi * pos_[index_i][0] / PL) * - sin(governing_vibration_integer_y * Pi * pos_[index_i][1] / PH); - }; -}; -} // namespace SPH - -//---------------------------------------------------------------------- -// Main program starts here. -//---------------------------------------------------------------------- -int main(int ac, char *av[]) -{ - //---------------------------------------------------------------------- - // Build up an SPHSystem and IO environment. - //---------------------------------------------------------------------- - SPHSystem sph_system(system_domain_bounds, particle_spacing_ref); - sph_system.handleCommandlineOptions(ac, av)->setIOEnvironment(); - //---------------------------------------------------------------------- - // Creating bodies with corresponding materials and particles. - //---------------------------------------------------------------------- - SolidBody plate_body(sph_system, makeShared("PlateBody")); - plate_body.defineMaterial(rho0_s, c0, Youngs_modulus, poisson); - plate_body.generateParticles(); - - ObserverBody plate_observer(sph_system, "PlateObserver"); - plate_observer.generateParticles(observation_location); - //---------------------------------------------------------------------- - // Define body relation map. - // The contact map gives the topological connections between the bodies. - // Basically the the range of bodies to build neighbor particle lists. - // Generally, we first define all the inner relations, then the contact relations. - //---------------------------------------------------------------------- - InnerRelation plate_body_inner(plate_body); - ContactRelation plate_observer_contact(plate_observer, {&plate_body}); - //---------------------------------------------------------------------- - // Define the numerical methods used in the simulation. - // Note that there may be data dependence on the sequence of constructions. - // Generally, the geometric models or simple objects without data dependencies, - // such as gravity, should be initiated first. - // Then the major physical particle dynamics model should be introduced. - // Finally, the auxillary models such as time step estimator, initial condition, - // boundary condition and other constraints should be defined. - //---------------------------------------------------------------------- - SimpleDynamics initial_velocity(plate_body); - InteractionWithUpdate corrected_configuration(plate_body_inner); - InteractionDynamics plate_shear_acceleration(plate_body_inner); - Dynamics1Level plate_pressure_relaxation(plate_body_inner); - Dynamics1Level plate_density_relaxation(plate_body_inner); - Dynamics1Level plate_shear_stress_relaxation(plate_body_inner); - ReduceDynamics fluid_advection_time_step(plate_body, U_ref, 0.2); - ReduceDynamics fluid_acoustic_time_step(plate_body, 0.4); - BoundaryGeometry boundary_geometry(plate_body, "BoundaryGeometry"); - SimpleDynamics constrain_holder(boundary_geometry, Vecd(1.0, 1.0, 0.0)); - SimpleDynamics constrain_mass_center(plate_body, Vecd(1.0, 1.0, 0.0)); - //---------------------------------------------------------------------- - // Define the methods for I/O operations, observations - // and regression tests of the simulation. - //---------------------------------------------------------------------- - BodyStatesRecordingToVtp write_states(sph_system); - write_states.addToWrite(plate_body, "Pressure"); - write_states.addToWrite(plate_body, "Density"); - SimpleDynamics plate_von_mises_stress(plate_body); - write_states.addToWrite(plate_body, "VonMisesStress"); - RestartIO restart_io(sph_system); - ObservedQuantityRecording write_plate_displacement("Position", plate_observer_contact); - RegressionTestDynamicTimeWarping> write_kinetic_energy(plate_body); - //---------------------------------------------------------------------- - // Prepare the simulation with cell linked list, configuration - // and case specified initial condition if necessary. - //---------------------------------------------------------------------- - sph_system.initializeSystemCellLinkedLists(); - sph_system.initializeSystemConfigurations(); - initial_velocity.exec(); - constrain_holder.exec(); - corrected_configuration.exec(); - //---------------------------------------------------------------------- - // Setup for time-stepping control - //---------------------------------------------------------------------- - Real &physical_time = *sph_system.getSystemVariableDataByName("PhysicalTime"); - size_t number_of_iterations = 0; - int screen_output_interval = 500; - int restart_output_interval = screen_output_interval * 10; - Real end_time = 0.02; - Real output_period = end_time / 50.0; - //---------------------------------------------------------------------- - // Statistics for CPU time - //---------------------------------------------------------------------- - TickCount t1 = TickCount::now(); - TimeInterval interval; - //---------------------------------------------------------------------- - // First output before the main loop. - //---------------------------------------------------------------------- - write_states.writeToFile(0); - write_plate_displacement.writeToFile(0); - write_kinetic_energy.writeToFile(0); - //---------------------------------------------------------------------- - // Main loop starts here. - //---------------------------------------------------------------------- - while (physical_time < end_time) - { - Real integration_time = 0.0; - while (integration_time < output_period) - { - Real relaxation_time = 0.0; - Real advection_dt = fluid_advection_time_step.exec(); - - while (relaxation_time < advection_dt) - { - Real acoustic_dt = fluid_acoustic_time_step.exec(); - plate_shear_stress_relaxation.exec(acoustic_dt); - plate_pressure_relaxation.exec(acoustic_dt); - constrain_holder.exec(acoustic_dt); - plate_density_relaxation.exec(acoustic_dt); - plate_shear_acceleration.exec(acoustic_dt); - number_of_iterations++; - relaxation_time += acoustic_dt; - integration_time += acoustic_dt; - physical_time += acoustic_dt; - if (number_of_iterations % screen_output_interval == 0) - { - std::cout << "N=" << number_of_iterations << " Time: " - << physical_time << " advection_dt: " - << advection_dt << " acoustic_dt: " - << acoustic_dt << "\n"; - if (number_of_iterations % restart_output_interval == 0 && number_of_iterations != sph_system.RestartStep()) - restart_io.writeToFile(Real(number_of_iterations)); - } - } - plate_body.updateCellLinkedList(); - plate_body_inner.updateConfiguration(); - corrected_configuration.exec(); - } - plate_von_mises_stress.exec(); - write_plate_displacement.writeToFile(number_of_iterations); - write_kinetic_energy.writeToFile(number_of_iterations); - TickCount t2 = TickCount::now(); - write_states.writeToFile(); - TickCount t3 = TickCount::now(); - interval += t3 - t2; - } - - TickCount t4 = TickCount::now(); - TickCount::interval_t tt; - tt = t4 - t1 - interval; - std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; - - if (sph_system.GenerateRegressionData()) - { - write_kinetic_energy.generateDataBase(1.0e-3); - } - else - { - write_kinetic_energy.testResult(); - } - return 0; -} diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml deleted file mode 100644 index 0595b41a58..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_0_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml deleted file mode 100644 index eedda0a708..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_3_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml deleted file mode 100644 index 89c7502971..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_Run_5_result.xml +++ /dev/null @@ -1,9 +0,0 @@ - - - - - - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml deleted file mode 100644 index 9497963594..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_dtwdistance.xml +++ /dev/null @@ -1,4 +0,0 @@ - - - - diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat deleted file mode 100644 index de07de18db..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateBody_TotalKineticEnergy_runtimes.dat +++ /dev/null @@ -1,3 +0,0 @@ -true -6 -4 \ No newline at end of file diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py b/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py deleted file mode 100644 index 52acb2cabc..0000000000 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py +++ /dev/null @@ -1,37 +0,0 @@ -# !/usr/bin/env python3 -import os -import sys - -path = os.path.abspath('../../../../../PythonScriptStore/RegressionTest') -sys.path.append(path) -from regression_test_base_tool import SphinxsysRegressionTest - -""" -case name: test_3d_oscillating_plate_UL -""" - -case_name = "test_3d_oscillating_plate_UL" -body_name = "PlateBody" -parameter_name = "TotalKineticEnergy" - -number_of_run_times = 0 -converged = 0 -sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name) - - -while True: - print("Now start a new run......") - sphinxsys.run_case() - number_of_run_times += 1 - converged = sphinxsys.read_dat_file() - print("Please note: This is the", number_of_run_times, "run!") - if number_of_run_times <= 200: - if converged == "true": - print("The tested parameters of all variables are converged, and the run will stop here!") - break - elif converged != "true": - print("The tested parameters of", sphinxsys.sphinxsys_parameter_name, "are not converged!") - continue - else: - print("It's too many runs but still not converged, please try again!") - break