diff --git a/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h b/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h index e604de1ef3..258ea393d6 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h +++ b/tests/user_examples/extra_src/shared/general_continuum/all_continuum.h @@ -27,7 +27,7 @@ * @details The fluid dynamics algorithms begin for fluid bulk without boundary condition, * then algorithm interacting with wall is defined, further algorithms for multiphase flow interaction * built upon these basic algorithms. - * @author Chi Zhang and Xiangyu Hu + * @author Chi ZHang and Xiangyu Hu */ #pragma once diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/all_continuum_dynamics.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/all_continuum_dynamics.h index 57e8d46448..78eb5db7c0 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/all_continuum_dynamics.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/all_continuum_dynamics.h @@ -27,7 +27,7 @@ * @details The fluid dynamics algorithms begin for fluid bulk without boundary condition, * then algorithm interacting with wall is defined, further algorithms for multiphase flow interaction * built upon these basic algorithms. - * @author Chi Zhang and Xiangyu Hu + * @author Chi ZHang and Xiangyu Hu */ #pragma once diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.cpp index 3f6e41d3d9..b8520c28ca 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.cpp @@ -3,7 +3,7 @@ namespace SPH { //=================================================================================================// - namespace continuum_dyannmics + namespace continuum_dynamics { } diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h index eac408d174..6db8a9da24 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.h @@ -8,7 +8,7 @@ namespace SPH { -namespace continuum_dyannmics +namespace continuum_dynamics { typedef DataDelegateContact @@ -51,7 +51,7 @@ class BaseShearStressRelaxation2ndHalfWithWall : public fluid_dynamics::Interact }; using ShearStressRelaxation2ndHalfWithWall = BaseShearStressRelaxation2ndHalfWithWall; -} // namespace continuum_dyannmics +} // namespace continuum_dynamics } // namespace SPH #endif // CONTINUUM_DYNAMICS_COMPLEX_H diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp index c45c65f9b5..37e7286f35 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_complex.hpp @@ -5,7 +5,7 @@ namespace SPH { //=================================================================================================// -namespace continuum_dyannmics +namespace continuum_dynamics { //=================================================================================================// //==========================BaseShearStressRelaxation1stHalfWithWall================================// @@ -16,7 +16,7 @@ void BaseShearStressRelaxation1stHalfWithWallshear_stress_[index_i]; Real rho_i = this->rho_[index_i]; - Real rho_in_wall = this->granular_material_.getDensity(); + Real rho_in_wall = this->continuum_.getDensity(); Vecd acceleration = Vecd::Zero(); for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) // There may be several wall bodies. { @@ -48,7 +48,7 @@ void BaseShearStressRelaxation2ndHalfWithWallvel_[index_i]; Matd velocity_gradient = Matd::Zero(); - Real rho_in_wall = this->granular_material_.getDensity(); + Real rho_in_wall = this->continuum_.getDensity(); for (size_t k = 0; k < fluid_dynamics::FluidWallData::contact_configuration_.size(); ++k) { @@ -70,5 +70,5 @@ void BaseShearStressRelaxation2ndHalfWithWallvelocity_gradient_[index_i] += velocity_gradient; } -} // namespace continuum_dyannmics +} // namespace continuum_dynamics } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp index 7ecf69500e..be9da45338 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.cpp @@ -2,7 +2,7 @@ namespace SPH { //=================================================================================================// - namespace continuum_dyannmics + namespace continuum_dynamics { ContinuumInitialCondition::ContinuumInitialCondition(SPHBody &sph_body) : LocalDynamics(sph_body), ContinuumDataSimple(sph_body), @@ -23,7 +23,7 @@ namespace SPH //=================================================================================================// BaseRelaxation::BaseRelaxation(BaseInnerRelation &inner_relation) : LocalDynamics(inner_relation.getSPHBody()), ContinuumDataInner(inner_relation), - granular_material_(particles_->granular_material_), rho_(particles_->rho_), + continuum_(particles_->continuum_), rho_(particles_->rho_), p_(*particles_->getVariableByName("Pressure")), drho_dt_(*particles_->getVariableByName("DensityChangeRate")), pos_(particles_->pos_), vel_(particles_->vel_), acc_(particles_->acc_), acc_prior_(particles_->acc_prior_) {} @@ -118,6 +118,7 @@ namespace SPH acceleration += rho_[index_j] * (shear_stress_i / (rho_i * rho_i) + shear_stress_[index_j] / (rho_[index_j] * rho_[index_j])) * nablaW_ijV_j; } acc_shear_[index_i] += acceleration; // for with artificial stress + //acc_shear_[index_i] = acceleration; // for without artificial stress } void ShearStressRelaxation1stHalf::update(size_t index_i, Real dt) @@ -151,10 +152,138 @@ namespace SPH void ShearStressRelaxation2ndHalf::update(size_t index_i, Real dt) { - shear_stress_rate_[index_i] = granular_material_.ConstitutiveRelationShearStress(velocity_gradient_[index_i], shear_stress_[index_i]); + 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); } - } // namespace continuum_dyannmics + + //=================================================================================================// + //===================================Non-hourglass formulation=====================================// + //=================================================================================================// + ShearAccelerationRelaxation::ShearAccelerationRelaxation(BaseInnerRelation& inner_relation) + : BaseRelaxation(inner_relation), + G_(continuum_.getShearModulus(continuum_.getYoungsModulus(), continuum_.getPoissonRatio())), + smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), shear_stress_(particles_->shear_stress_), + B_(*this->particles_->template registerSharedVariable("CorrectionMatrix", Matd::Identity())), acc_shear_(particles_->acc_shear_) {} + //original + void ShearAccelerationRelaxation::interaction(size_t index_i, Real dt) + { + Real rho_i = rho_[index_i]; + Vecd acceleration = Vecd::Zero(); + Vecd vel_i = vel_[index_i]; + Neighborhood& inner_neighborhood = inner_configuration_[index_i]; + //Matd B_i = B_[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_ijV_j_[n]; + + Vecd v_ij = (vel_i - vel_[index_j]); + acceleration += 2 * v_ij * dW_ijV_j / (r_ij + 0.01 * smoothing_length_); + } + acc_shear_[index_i] += G_ * acceleration * dt / rho_i; + } + //=================================================================================================// + void AngularConservativeShearAccelerationRelaxation::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_ijV_j_[n]; + 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) + : BaseRelaxation(inner_relation), + shear_stress_(particles_->shear_stress_), shear_stress_rate_(particles_->shear_stress_rate_), + velocity_gradient_(particles_->velocity_gradient_), strain_tensor_(particles_->strain_tensor_), + strain_tensor_rate_(particles_->strain_tensor_rate_), von_mises_stress_(particles_->von_mises_stress_), + von_mises_strain_(particles_->von_mises_strain_), Vol_(particles_->Vol_), + B_(*particles_->getVariableByName("CorrectionMatrix")) {} + 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_ijV_j_[n]; + Vecd& e_ij = inner_neighborhood.e_ij_[n]; + + //Matd velocity_gradient_ij = - (vel_[index_i] - vel_[index_j]) * nablaW_ijV_j.transpose(); + Vecd v_ij = vel_[index_i] - vel_[index_j]; + Matd velocity_gradient_ij = -v_ij * (B_[index_i] * e_ij * dW_ijV_j_).transpose(); + velocity_gradient += velocity_gradient_ij; + } + 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_[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); + } + 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; + + //VonMises Stress + Matd stress_tensor_i = shear_stress_[index_i] - p_[index_i] * Matd::Identity(); + von_mises_stress_[index_i] = getVonMisesStressFromMatrix(stress_tensor_i); + } + + //=================================================================================================// + FixedInAxisDirection::FixedInAxisDirection(BodyPartByParticle& body_part, Vecd constrained_axises) + : BaseMotionConstraint(body_part), constrain_matrix_(Matd::Identity()) + { + for (int k = 0; k != Dimensions; ++k) + constrain_matrix_(k, k) = constrained_axises[k]; + }; + //=================================================================================================// + void FixedInAxisDirection::update(size_t index_i, Real dt) + { + vel_[index_i] = constrain_matrix_ * vel_[index_i]; + }; + + //=================================================================================================// + ConstrainSolidBodyMassCenter:: + ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction) + : LocalDynamics(sph_body), ContinuumDataSimple(sph_body), + correction_matrix_(Matd::Identity()), vel_(particles_->vel_), + compute_total_momentum_(sph_body, "Velocity") + { + for (int i = 0; i != Dimensions; ++i) + correction_matrix_(i, i) = constrain_direction[i]; + ReduceDynamics> compute_total_mass_(sph_body, "MassiveMeasure"); + total_mass_ = compute_total_mass_.exec(); + } + //=================================================================================================// + void ConstrainSolidBodyMassCenter::setupDynamics(Real dt) + { + velocity_correction_ = + correction_matrix_ * compute_total_momentum_.exec(dt) / total_mass_; + } + //=================================================================================================// + void ConstrainSolidBodyMassCenter::update(size_t index_i, Real dt) + { + vel_[index_i] -= velocity_correction_; + } + } // namespace continuum_dynamics } // namespace SPH diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h index a6ce9a2b80..9f9919a62f 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.h @@ -7,7 +7,7 @@ namespace SPH { -namespace continuum_dyannmics +namespace continuum_dynamics { typedef DataDelegateSimple ContinuumDataSimple; typedef DataDelegateInner ContinuumDataInner; @@ -42,7 +42,7 @@ class BaseRelaxation : public LocalDynamics, public ContinuumDataInner virtual ~BaseRelaxation(){}; protected: - GeneralContinuum &granular_material_; + GeneralContinuum &continuum_; StdLargeVec &rho_, &p_, &drho_dt_; StdLargeVec &pos_, &vel_, &acc_, &acc_prior_; }; @@ -110,6 +110,77 @@ class ShearStressRelaxation2ndHalf : public BaseRelaxation StdLargeVec &von_mises_stress_; }; +//=================================================================================================// +//===================================Non-hourglass formulation=====================================// +//=================================================================================================// +/** + * @class Integration1stHalf + * @brief Pressure relaxation scheme with the mostly used Riemann solver. +*/ +template +class BaseIntegration1stHalf : public FluidDynamicsType +{ +public: + explicit BaseIntegration1stHalf(BaseInnerRelation& inner_relation); + virtual ~BaseIntegration1stHalf() {}; + 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: + StdLargeVec& acc_shear_; +}; +using Integration1stHalf = BaseIntegration1stHalf; +using Integration1stHalfRiemann = BaseIntegration1stHalf; + +/** +* @class ShearAccelerationRelaxation +*/ +class ShearAccelerationRelaxation : public BaseRelaxation +{ +public: + explicit ShearAccelerationRelaxation(BaseInnerRelation& inner_relation); + virtual ~ShearAccelerationRelaxation() {}; + void interaction(size_t index_i, Real dt = 0.0); +protected: + Real G_, smoothing_length_; + StdLargeVec& shear_stress_, & B_; + StdLargeVec& acc_shear_; +}; + +/** + * @class AngularConservativeShearAccelerationRelaxation + */ +class AngularConservativeShearAccelerationRelaxation : public ShearAccelerationRelaxation +{ +public: + explicit AngularConservativeShearAccelerationRelaxation(BaseInnerRelation& inner_relation) + : ShearAccelerationRelaxation(inner_relation) {}; + virtual ~AngularConservativeShearAccelerationRelaxation() {}; + + void interaction(size_t index_i, Real dt = 0.0); + //void update(size_t index_i, Real dt = 0.0); + +}; + +/** +* @class ShearStressRelaxation +*/ +class ShearStressRelaxation : public BaseRelaxation +{ +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: + StdLargeVec& shear_stress_, & shear_stress_rate_, & velocity_gradient_, & strain_tensor_, & strain_tensor_rate_; + StdLargeVec& von_mises_stress_, & von_mises_strain_, & Vol_; + StdLargeVec& B_; +}; + /** * @class BaseMotionConstraint */ @@ -145,6 +216,44 @@ class FixConstraint : public BaseMotionConstraint }; using FixBodyConstraint = FixConstraint; using FixBodyPartConstraint = FixConstraint; -} // namespace continuum_dyannmics + +/** + * @class FixedInAxisDirection + * @brief Constrain the velocity of a solid body part. + */ +class FixedInAxisDirection : public BaseMotionConstraint +{ +public: + FixedInAxisDirection(BodyPartByParticle& body_part, Vecd constrained_axises = Vecd::Zero()); + virtual ~FixedInAxisDirection() {}; + void update(size_t index_i, Real dt = 0.0); + +protected: + Matd constrain_matrix_; +}; + +/** + * @class ConstrainSolidBodyMassCenter + * @brief Constrain the mass center of a solid body. + */ +class ConstrainSolidBodyMassCenter : public LocalDynamics, public ContinuumDataSimple +{ +private: + Real total_mass_; + Matd correction_matrix_; + Vecd velocity_correction_; + StdLargeVec& vel_; + ReduceDynamics> compute_total_momentum_; + +protected: + virtual void setupDynamics(Real dt = 0.0) override; + +public: + explicit ConstrainSolidBodyMassCenter(SPHBody& sph_body, Vecd constrain_direction = Vecd::Ones()); + virtual ~ConstrainSolidBodyMassCenter() {}; + + void update(size_t index_i, Real dt = 0.0); +}; +} // namespace continuum_dynamics } // namespace SPH #endif // CONTINUUM_DYNAMICS_INNER_H \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp index b55251dc87..a653c918c6 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_dynamics/continuum_dynamics_inner.hpp @@ -2,8 +2,27 @@ #include "continuum_dynamics_inner.h" namespace SPH { - namespace continuum_dyannmics + namespace continuum_dynamics { + template + BaseIntegration1stHalf::BaseIntegration1stHalf(BaseInnerRelation& inner_relation) + : FluidDynamicsType(inner_relation), + acc_shear_(*this->particles_->template getVariableByName("AccelerationByShear")) {} + template + void BaseIntegration1stHalf::initialization(size_t index_i, Real dt) + { + FluidDynamicsType::initialization(index_i, dt); + } + template + void BaseIntegration1stHalf::interaction(size_t index_i, Real dt) + { + FluidDynamicsType::interaction(index_i, dt); + } + template + void BaseIntegration1stHalf::update(size_t index_i, Real dt) + { + this->vel_[index_i] += (this->acc_prior_[index_i] + this->acc_[index_i] + this->acc_shear_[index_i]) * dt; + } } } // namespace SPH \ No newline at end of file diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp index 55f6f150a4..c195c1f653 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.cpp @@ -3,8 +3,8 @@ namespace SPH { ContinuumParticles:: - ContinuumParticles(SPHBody &sph_body, GeneralContinuum *granular_material) - : BaseParticles(sph_body, granular_material), granular_material_(*granular_material) {} + ContinuumParticles(SPHBody &sph_body, GeneralContinuum *continuum) + : BaseParticles(sph_body, continuum), continuum_(*continuum) {} //=================================================================================================// void ContinuumParticles::initializeOtherVariables() { @@ -16,6 +16,7 @@ namespace SPH registerVariable(shear_stress_, "ShearStress"); registerVariable(shear_stress_rate_, "ShearStressRate"); registerVariable(von_mises_stress_, "VonMisesStress"); + registerVariable(von_mises_strain_, "VonMisesStrain"); registerVariable(velocity_gradient_, "VelocityGradient"); registerVariable(strain_tensor_, "StrainTensor"); registerVariable(strain_tensor_rate_, "StrainTensorRate"); @@ -29,6 +30,7 @@ namespace SPH registerSortableVariable("ShearStress"); registerSortableVariable("ShearStressRate"); registerSortableVariable("VonMisesStress"); + registerSortableVariable("VonMisesStrain"); registerSortableVariable("VelocityGradient"); registerSortableVariable("StrainTensor"); registerSortableVariable("StrainTensorRate"); diff --git a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h index 1f87c71399..6f5ed620a8 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h +++ b/tests/user_examples/extra_src/shared/general_continuum/continuum_particles.h @@ -21,14 +21,15 @@ namespace SPH StdLargeVec velocity_gradient_; StdLargeVec von_mises_stress_; + StdLargeVec von_mises_strain_; StdLargeVec pos0_; /**< initial position */ StdLargeVec n_; /**< current normal direction */ StdLargeVec n0_; /**< initial normal direction */ - GeneralContinuum &granular_material_; + GeneralContinuum &continuum_; - ContinuumParticles(SPHBody &sph_body, GeneralContinuum *granular_material); + ContinuumParticles(SPHBody &sph_body, GeneralContinuum *continuum); virtual ~ContinuumParticles(){}; virtual void initializeOtherVariables() override; diff --git a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h index d91953deb0..0e016ba79f 100644 --- a/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h +++ b/tests/user_examples/extra_src/shared/general_continuum/general_continuum.h @@ -12,10 +12,11 @@ class GeneralContinuum : public WeaklyCompressibleFluid Real G_; /*< shearmodules */ Real K_; /*< bulkmodules */ Real nu_; /*< Poisson ratio */ + 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) + : WeaklyCompressibleFluid(rho0, c0), E_(0.0), G_(0.0), K_(0.0), nu_(0.0), contact_stiffness_(c0* c0) { material_type_name_ = "GeneralContinuum"; E_ = youngs_modulus; @@ -34,6 +35,8 @@ class GeneralContinuum : public WeaklyCompressibleFluid Real getShearModulus(Real youngs_modulus, Real poisson_ratio); Real getLambda(Real youngs_modulus, Real poisson_ratio); + Real ContactStiffness() { return contact_stiffness_; }; + virtual Matd ConstitutiveRelationShearStress(Matd &velocity_gradient, Matd &shear_stress); virtual GeneralContinuum *ThisObjectPtr() override { return this; }; diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/user_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp new file mode 100644 index 0000000000..e8c5f26b96 --- /dev/null +++ b/tests/user_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -0,0 +1,240 @@ +/* ---------------------------------------------------------------------------* +* SPHinXsys: 2D oscillation beam example-update Lagrange * +* ----------------------------------------------------------------------------* +* This is the one of the basic test cases for understanding SPH method for * +* solid simulation based on update Lagrange method * +* In this case, the constraint of the beam is implemented with * +* internal constrained subregion. * +* ----------------------------------------------------------------------------*/ +#include "sphinxsys.h" +#include "all_continuum.h" +using namespace SPH; +//------------------------------------------------------------------------------ +//global parameters for the case +//------------------------------------------------------------------------------ +Real PL = 0.2; //beam length +Real PH = 0.02; //for thick plate; =0.01 for thin 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 poisson = 0.4; //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); +//for dual time-step +Real U_max = 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 system(system_domain_bounds, resolution_ref); +#ifdef BOOST_AVAILABLE + // handle command line arguments + system.handleCommandlineOptions(ac, av); +#endif + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + RealBody beam_body(system, makeShared("BeamBody")); + beam_body.sph_adaptation_->resetKernel(); + beam_body.defineParticlesAndMaterial(rho0_s, c0, Youngs_modulus, poisson); + beam_body.generateParticles(); + beam_body.addBodyStateForRecording("VonMisesStress"); + beam_body.addBodyStateForRecording("VonMisesStrain"); + beam_body.addBodyStateForRecording("Pressure"); + beam_body.addBodyStateForRecording("Density"); + + ObserverBody beam_observer(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. + //---------------------------------------------------------------------- + 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); + SharedPtr gravity_ptr = makeShared(Vecd(0.0, -gravity_g)); + Dynamics1Level beam_pressure_relaxation(beam_body_inner); + Dynamics1Level beam_density_relaxation(beam_body_inner); + InteractionDynamics + beam_shear_acceleration_angular_conservative(beam_body_inner); + InteractionWithUpdate correcttion_matrix(beam_body_inner); + Dynamics1Level beam_shear_stress_relaxation(beam_body_inner); + //for dula timestep + ReduceDynamics fluid_advection_time_step(beam_body, U_max, 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 + //----------------------------------------------------------------------------- + IOEnvironment io_environment(system); + BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); + RegressionTestEnsembleAverage> + write_beam_tip_displacement("Position", io_environment, beam_observer_contact); + RegressionTestDynamicTimeWarping>> + write_water_mechanical_energy(io_environment, beam_body, gravity_ptr); + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + beam_initial_velocity.exec(); + correcttion_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 begines + //----------------------------------------------------------------------------- + write_beam_states.writeToFile(0); + write_beam_tip_displacement.writeToFile(0); + write_water_mechanical_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(); + + while (relaxation_time < advection_dt) + { + Real acoustic_dt = fluid_acoustic_time_step.exec(); + beam_shear_stress_relaxation.exec(acoustic_dt); + beam_pressure_relaxation.exec(acoustic_dt); + constraint_beam_base.exec(); + beam_density_relaxation.exec(acoustic_dt); + //shear acceleration with angular conservative + beam_shear_acceleration_angular_conservative.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(); + correcttion_matrix.exec(); + } + write_beam_tip_displacement.writeToFile(ite); + write_water_mechanical_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; + + //system.generate_regression_data_ = true; + if (system.generate_regression_data_) + { + write_beam_tip_displacement.generateDataBase(Vec2d(1.0e-2, 1.0e-2), Vec2d(1.0e-2, 1.0e-2)); + } + else + { + write_beam_tip_displacement.testResult(); + } + return 0; +} \ No newline at end of file diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalMechanicalEnergy_runtimes.dat b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalMechanicalEnergy_runtimes.dat new file mode 100644 index 0000000000..7b9618df5d --- /dev/null +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamBody_TotalMechanicalEnergy_runtimes.dat @@ -0,0 +1,3 @@ +false +7 +0 \ No newline at end of file diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position.xml b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position.xml new file mode 100644 index 0000000000..aa8121c730 --- /dev/null +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position.xml @@ -0,0 +1,103 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml index d2ebf5ca38..f47aa858b4 100644 --- a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_ensemble_averaged_mean_variance.xml @@ -1,209 +1,207 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_result.xml b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_result.xml new file mode 100644 index 0000000000..b7abab0b6c --- /dev/null +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_result.xml @@ -0,0 +1,615 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_runtimes.dat b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_runtimes.dat index 1284659881..0320826438 100644 --- a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_runtimes.dat +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/BeamObserver_Position_runtimes.dat @@ -1,3 +1,3 @@ true -10 +6 5 \ No newline at end of file diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/regression_test_tool.py index e09a211fe3..eb8359898f 100644 --- a/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/regression_test_tool.py +++ b/tests/user_examples/test_2d_oscillating_beam_UL/regression_test_tool/regression_test_tool.py @@ -10,7 +10,7 @@ case name: test_2d_oscillating_beam_UL """ -case_name = "test_2d_oscillating_beam_UL" +case_name = "zsh_2d_oscillating_beam_UL_non-hourglass" body_name = "BeamObserver" parameter_name = "Position" diff --git a/tests/user_examples/test_2d_oscillating_beam_UL/test_2d_oscillating_beam_UL.cpp b/tests/user_examples/test_2d_oscillating_beam_UL/test_2d_oscillating_beam_UL.cpp deleted file mode 100644 index a3271df045..0000000000 --- a/tests/user_examples/test_2d_oscillating_beam_UL/test_2d_oscillating_beam_UL.cpp +++ /dev/null @@ -1,226 +0,0 @@ -/* ---------------------------------------------------------------------------* - * SPHinXsys: 2D oscillation beam example-update Lagrange * - * ----------------------------------------------------------------------------* - * This is the one of the basic test cases for understanding SPH method for * - * solid simulation based on update Lagrange method * - * In this case, the constraint of the beam is implemented with * - * internal constrained subregion. * - * ----------------------------------------------------------------------------*/ -#include "all_continuum.h" -#include "sphinxsys.h" -using namespace SPH; -//------------------------------------------------------------------------------ -// global parameters for the case -//------------------------------------------------------------------------------ -Real PL = 0.2; // beam length -Real PH = 0.02; // for thick plate; =0.01 for thin 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); -//---------------------------------------------------------------------- -// 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 &granular_column) - : fluid_dynamics::FluidInitialCondition(granular_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 system(system_domain_bounds, resolution_ref); -#ifdef BOOST_AVAILABLE - // handle command line arguments - system.handleCommandlineOptions(ac, av); -#endif //---------------------------------------------------------------------- - //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - RealBody beam_body(system, makeShared("BeamBody")); - beam_body.sph_adaptation_->resetKernel(); - beam_body.defineParticlesAndMaterial(rho0_s, c0, Youngs_modulus, poisson); - beam_body.generateParticles(); - beam_body.addBodyStateForRecording("VonMisesStress"); - - ObserverBody beam_observer(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. - //---------------------------------------------------------------------- - 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); - Dynamics1Level granular_pressure_relaxation(beam_body_inner); - Dynamics1Level granular_density_relaxation(beam_body_inner); - Dynamics1Level granular_shear_stress_relaxation(beam_body_inner); - InteractionWithUpdate granular_shear_stress_rate_relaxation(beam_body_inner); - InteractionDynamics granular_artificial_normal_shear_stress_relaxation(beam_body_inner); - ReduceDynamics computing_time_step_size(beam_body); - // clamping a solid body part. - BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics> constraint_beam_base(beam_base); - //----------------------------------------------------------------------------- - // outputs - //----------------------------------------------------------------------------- - IOEnvironment io_environment(system); - BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); - RegressionTestEnsembleAverage> - write_beam_tip_displacement("Position", io_environment, beam_observer_contact); - //---------------------------------------------------------------------- - // Setup computing and initial conditions. - //---------------------------------------------------------------------- - system.initializeSystemCellLinkedLists(); - system.initializeSystemConfigurations(); - beam_initial_velocity.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 */ - Real Dt = 0.1 * D_Time; - Real dt = 0.0; // default acoustic time step sizes - TickCount t1 = TickCount::now(); - TimeInterval interval; - //----------------------------------------------------------------------------- - // from here the time stepping begines - //----------------------------------------------------------------------------- - write_beam_states.writeToFile(0); - write_beam_tip_displacement.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; - while (relaxation_time < Dt) - { - granular_pressure_relaxation.exec(dt); - granular_artificial_normal_shear_stress_relaxation.exec(); - granular_shear_stress_relaxation.exec(dt); - constraint_beam_base.exec(); - granular_shear_stress_rate_relaxation.exec(dt); - granular_density_relaxation.exec(dt); - - ite++; - dt = computing_time_step_size.exec(); - relaxation_time += dt; - integration_time += dt; - GlobalStaticVariables::physical_time_ += dt; - if (ite % 100 == 0) - { - std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ << " dt: " - << dt << "\n"; - } - beam_body.updateCellLinkedList(); - beam_body_inner.updateConfiguration(); - } - } - write_beam_tip_displacement.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; - - // system.generate_regression_data_ = true; - if (system.generate_regression_data_) - { - write_beam_tip_displacement.generateDataBase(Vec2d(3.0e-2, 3.0e-2), Vec2d(1.0e-2, 1.0e-2)); - } - else - { - write_beam_tip_displacement.testResult(); - } - return 0; -} \ No newline at end of file diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/CMakeLists.txt b/tests/user_examples/test_3d_oscillating_plate_UL/CMakeLists.txt new file mode 100644 index 0000000000..0ee37144e2 --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/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} extra_sources_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME} + COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) \ No newline at end of file diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/user_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp new file mode 100644 index 0000000000..38790ca431 --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp @@ -0,0 +1,255 @@ +/** +* @file dw_3d_volum_oscillating_plate.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 Dong Wu, Chi Zhang and Xiangyu Hu +* @version 0.1 + */ +#include "sphinxsys.h" +#include "all_continuum.h" +using namespace SPH; +/** + * @brief 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_nuber_PL = PL / particle_spacing_ref; +int particle_nuber_PH = PH / particle_spacing_ref; +int BWD = 1; +Real BW = particle_spacing_ref * (Real)BWD; /** Boundary width, determined by specific layer of boundary particles. */ +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec3d(-BW, -BW, -PT / 2), Vec3d(PL + BW, PH + BW, PT / 2)); +// Observer location +StdVec observation_location = { Vecd(0.5 * PL, 0.5 * PH, 0.0), Vecd(-BW, -BW, 0.0) }; +/** For material properties of the solid. */ +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 = 1.0; +Real governing_vibration_integer_y = 1.0; +Real U_max = 1.0; //Maximum velocity +/** Define application dependent particle generator for thin structure. */ +class PlateParticleGenerator : public ParticleGenerator +{ +public: + explicit PlateParticleGenerator(SPHBody& sph_body) : ParticleGenerator(sph_body) {}; + virtual void initializeGeometricVariables() override + { + for (int k = 0; k < particle_number; k++) + { + for (int i = 0; i < particle_nuber_PL + 2 * BWD; i++) + { + for (int j = 0; j < particle_nuber_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)); + initializePositionAndVolumetricMeasure(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_.pos_[index_i][2] < 0.5 * particle_spacing_ref) + & (base_particles_.pos_[index_i][2] > -0.5 * particle_spacing_ref) + & (base_particles_.pos_[index_i][0] < 0.0 || base_particles_.pos_[index_i][0] > PL + || base_particles_.pos_[index_i][1] < 0.0 || base_particles_.pos_[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); + }; +}; + +/** + * The main program + */ +int main() +{ + /** Setup the system. */ + SPHSystem system(system_domain_bounds, particle_spacing_ref); + /** Tag for computation from restart files. 0: start with initial condition. */ + system.setRestartStep(0); + + /** create a plate body. */ + SolidBody plate_body(system, makeShared("PlateBody")); + plate_body.defineParticlesAndMaterial(rho0_s, c0, Youngs_modulus, poisson); + plate_body.generateParticles(); + //plate_body.addBodyStateForRecording("VolumetricStress"); + plate_body.addBodyStateForRecording("VonMisesStress"); + plate_body.addBodyStateForRecording("VonMisesStrain"); + plate_body.addBodyStateForRecording("Pressure"); + plate_body.addBodyStateForRecording("Density"); + + /** Define Observer. */ + ObserverBody plate_observer(system, "PlateObserver"); + plate_observer.defineParticlesAndMaterial(); + plate_observer.generateParticles(observation_location); + + /** Set body contact map + * The contact map gives the data conntections between the bodies + * basically the the range of bodies to build neighbor particle lists + */ + InnerRelation plate_body_inner(plate_body); + ContactRelation plate_observer_contact(plate_observer, { &plate_body }); + /** + * This section define all numerical methods will be used in this case. + */ + SimpleDynamics initial_velocity(plate_body); + /** Time step size calculation. */ + ReduceDynamics fluid_advection_time_step(plate_body, U_max, 0.2); + ReduceDynamics fluid_acoustic_time_step(plate_body, 0.4); + /** stress relaxation. */ + Dynamics1Level plate_pressure_relaxation(plate_body_inner); + Dynamics1Level plate_density_relaxation(plate_body_inner); + InteractionDynamics + plate_shear_acceleration_angular_conservative(plate_body_inner); + /** Corrected configuration. */ + InteractionWithUpdate corrected_configuration(plate_body_inner); + Dynamics1Level plate_shear_stress_relaxation(plate_body_inner); + /** Constrain the Boundary. */ + 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)); + /** Output */ + IOEnvironment io_environment(system); + BodyStatesRecordingToVtp write_states(io_environment, system.real_bodies_); + RestartIO restart_io(io_environment, system.real_bodies_); + // ObservedQuantityRecording + // write_plate_displacement("Position", io_environment, plate_observer_contact); + RegressionTestEnsembleAverage> + write_plate_displacement("Position", io_environment, plate_observer_contact); + ReducedQuantityRecording> + write_kinetic_energy(io_environment, plate_body); + + /** Apply initial condition. */ + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + initial_velocity.exec(); + constrain_holder.exec(); + corrected_configuration.exec(); + + /** + * @brief The time stepping starts here. + */ + if (system.RestartStep() != 0) + { + GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(system.RestartStep()); + } + GlobalStaticVariables::physical_time_ = 0.0; + /** first output. */ + write_states.writeToFile(0); + write_plate_displacement.writeToFile(0); + write_kinetic_energy.writeToFile(0); + //write_elastic_strain_energy.writeToFile(0); + + /** Setup physical parameters. */ + size_t number_of_iterations = system.RestartStep(); + int screen_output_interval = 500; + int restart_output_interval = screen_output_interval * 10; + Real end_time = 0.1; + Real output_period = end_time / 100.0; + /** Statistics for computing time. */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + /** + * Main 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 = 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); + //shear acceleration with angular conservative + plate_shear_acceleration_angular_conservative.exec(acoustic_dt); + number_of_iterations++; + relaxation_time += acoustic_dt; + integration_time += acoustic_dt; + GlobalStaticVariables::physical_time_ += acoustic_dt; + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << "N=" << number_of_iterations << " Time: " + << GlobalStaticVariables::physical_time_ << " advection_dt: " + << advection_dt << " acoustic_dt: " + << acoustic_dt << "\n"; + if (number_of_iterations % restart_output_interval == 0 && number_of_iterations != system.RestartStep()) + restart_io.writeToFile(Real(number_of_iterations)); + } + } + plate_body.updateCellLinkedList(); + plate_body_inner.updateConfiguration(); + corrected_configuration.exec(); //put correct matrix in the end + } + 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; + + //system.generate_regression_data_ = true; + if (system.generate_regression_data_) + { + write_plate_displacement.generateDataBase(Vecd(1.0e-2, 1.0e-2, 1.0e-2), Vecd(1.0e-2, 1.0e-2, 1.0e-2)); + } + else + { + write_plate_displacement.testResult(); + } + return 0; +} \ No newline at end of file diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position.xml b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position.xml new file mode 100644 index 0000000000..d9b57e3ea2 --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position.xml @@ -0,0 +1,70 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_ensemble_averaged_mean_variance.xml b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_ensemble_averaged_mean_variance.xml new file mode 100644 index 0000000000..fd184fceaa --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_ensemble_averaged_mean_variance.xml @@ -0,0 +1,141 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_result.xml b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_result.xml new file mode 100644 index 0000000000..2ce2e2e3f1 --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_result.xml @@ -0,0 +1,417 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_runtimes.dat b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_runtimes.dat new file mode 100644 index 0000000000..0320826438 --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/PlateObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +5 \ No newline at end of file diff --git a/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_3d_oscillating_plate_UL/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..57e3caa3bc --- /dev/null +++ b/tests/user_examples/test_3d_oscillating_plate_UL/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 +""" + +case_name = "zsh_3d_oscillating_plate" +body_name = "PlateObserver" +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_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