From 5801ae50e762ea7b902d396fe0daa40b6dfef8f3 Mon Sep 17 00:00:00 2001 From: Tang-U-b Date: Fri, 1 Sep 2023 21:39:18 +0200 Subject: [PATCH] modify the comments and spell error checks --- src/shared/kernels/anisotropic_kernel.h | 4 +- src/shared/kernels/base_kernel.h | 4 +- .../solid_dynamics/elastic_dynamics.cpp | 2 +- .../test_2d_anisotropic_beam.cpp | 360 +++++++++--------- .../oscillating_beam.cpp | 2 +- 5 files changed, 185 insertions(+), 187 deletions(-) diff --git a/src/shared/kernels/anisotropic_kernel.h b/src/shared/kernels/anisotropic_kernel.h index 6820935465..43ea4112bc 100644 --- a/src/shared/kernels/anisotropic_kernel.h +++ b/src/shared/kernels/anisotropic_kernel.h @@ -26,7 +26,7 @@ * implemented in derived classes. The kernel function define the relevance * between two neighboring particles. Basically, the further the two * particles, the less relevance they have. - * @author Zhentong Wang and Xiaojing Tang + * @author Xiaojing Tang and Zhentong Wang * @version 0.1 * @version 0.3.0 */ @@ -88,9 +88,11 @@ class AnisotropicKernel : public KernelType this->factor_d2W_3D_ = this->factor_W_3D_; }; + /** Calculates the transform tensor form anisotropic space to isotropic space **/ Mat2d getCoordinateTransformationTensorG(Vec2d kernel_vector, Vec2d transform_vector); Mat3d getCoordinateTransformationTensorG(Vec3d kernel_vector, Vec3d transform_vector); + /** Calculates the unit vector between a pair of particles **/ virtual Vec2d e(const Real &distance, const Vec2d &displacement) const override; virtual Vec3d e(const Real &distance, const Vec3d &displacement) const override; diff --git a/src/shared/kernels/base_kernel.h b/src/shared/kernels/base_kernel.h index 3c9e40507a..ab9eb01f8b 100644 --- a/src/shared/kernels/base_kernel.h +++ b/src/shared/kernels/base_kernel.h @@ -96,7 +96,7 @@ class Kernel virtual Vec3d e(const Real &distance, const Vec3d &displacement) const { return displacement / (distance + TinyReal); }; /** - * check if particles are within cutOff radius + * check if particles are within cutoff radius */ virtual bool checkIfWithinCutOffRadius(Vec2d displacement) { @@ -216,4 +216,4 @@ class Kernel void reduceTwice(); /** reduce for linear structures or filaments */ }; } // namespace SPH -#endif // BASE_KERNELS_H \ No newline at end of file +#endif // BASE_KERNELS_H diff --git a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp index 37d66c2b1a..bf8182952b 100644 --- a/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/elastic_dynamics.cpp @@ -79,7 +79,7 @@ void Integration1stHalfPK2::initialization(size_t index_i, Real dt) F_[index_i] += dF_dt_[index_i] * dt * 0.5; rho_[index_i] = rho0_ / F_[index_i].determinant(); // obtain the first Piola-Kirchhoff stress from the second Piola-Kirchhoff stress - // it seems using reproducing correction here increases convergence rate near the free surface + // it seems using reproducing correction here increases convergence rate near the free surface, note that the correction matrix is in a form of transpose stress_PK1_B_[index_i] = elastic_solid_.StressPK1(F_[index_i], index_i) * B_[index_i].transpose(); } //=================================================================================================// diff --git a/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp index 045fc46666..39a510c19b 100644 --- a/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp +++ b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp @@ -1,65 +1,63 @@ /* ---------------------------------------------------------------------------* * SPHinXsys: 2D anisotropic beam example-one body version * * ----------------------------------------------------------------------------* - * This is the one of the basic test cases, also the first case for * - * understanding SPH method for solid simulation. * - * In this case, the constraint of the beam is implemented with * - * internal constrained subregion. * + * This is a test cases using ASPH method with anisotropic kerne for * + * simulating solid. Particle space is anisotropic in different directions in this beam. * + * @author Xiaojing Tang * ----------------------------------------------------------------------------*/ #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.02; // depth of the insert +Real PL = 0.2; // beam length +Real PH = 0.02; // beam width +Real SL = 0.02; // constrained length -int y_num = 10; -Real ratio_ = 4.0; -//reference particle spacing -Real resolution_ref = PH / y_num; -Real resolution_ref_large = ratio_* resolution_ref; +// particle spacings and particle numbers +int y_num = 10; // particle number in y direction +Real ratio_ = 4.0; // anisotropic ratio, also dp_x / dp_y +Real resolution_ref = PH / y_num; // particle spacing in y direction +Real resolution_ref_large = ratio_ * resolution_ref; // large particle spacing, also the particle spacing in x direction +Real Total_PL = PL + SL; // total length +int x_num = Total_PL / resolution_ref_large; // particle number in x direction +// anisotropic parameters +Vec2d scaling_vector = Vec2d(1.0, 1.0 / ratio_); // scaling_vector for defining the anisotropic kernel +Real scaling_factor = 1.0 / ratio_; // scaling factor to calculate the time step -Real Total_PL = PL + SL; //beam length -int x_num = Total_PL / resolution_ref_large; - -Vec2d scaling_vector = Vec2d (1.0, 1.0 / ratio_); -Real scaling_factor = 1.0 / ratio_; - -class NonisotropicParticleGenerator : public ParticleGenerator +//---------------------------------------------------------------------- +// particle generation considering the anisotropic resolution +//---------------------------------------------------------------------- +class AnisotropicParticleGenerator : public ParticleGenerator { -public: - NonisotropicParticleGenerator(SPHBody &sph_body) : ParticleGenerator(sph_body) {}; - - virtual void initializeGeometricVariables() override - { - // set particles directly - for (int i = 0; i < x_num; i++) - { - for (int j = 0; j < y_num; j++) - { - Real x = -SL + (i + 0.5) * resolution_ref_large; - Real y = -PH / 2+ (j + 0.5)* resolution_ref; - initializePositionAndVolumetricMeasure(Vecd(x, y), (resolution_ref * resolution_ref_large)); - } + public: + AnisotropicParticleGenerator(SPHBody &sph_body) : ParticleGenerator(sph_body){}; - } - } + virtual void initializeGeometricVariables() override + { + // set particles directly + for (int i = 0; i < x_num; i++) + { + for (int j = 0; j < y_num; j++) + { + Real x = -SL + (i + 0.5) * resolution_ref_large; + Real y = -PH / 2 + (j + 0.5) * resolution_ref; + initializePositionAndVolumetricMeasure(Vecd(x, y), (resolution_ref * resolution_ref_large)); + } + } + } }; -// reference particle spacing -//Real resolution_ref = PH / 4.0; 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)); + Vec2d(PL + 3.0 * BW, PL / 2.0)); //---------------------------------------------------------------------- -// Material properties of the fluid. +// Material properties of the solid. //---------------------------------------------------------------------- -Real rho0_s = 1.0e3; // reference density +Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus -Real poisson = 0.3975; // Poisson ratio +Real poisson = 0.3975; // Poisson ratio //---------------------------------------------------------------------- // Parameters for initial condition on velocity //---------------------------------------------------------------------- @@ -74,11 +72,11 @@ Real R = PL / (0.5 * Pi); //---------------------------------------------------------------------- // a beam base shape std::vector beam_base_shape{ - Vecd(-SL , -PH / 2 ), Vecd(-SL , PH / 2 ), Vecd(0.0, PH / 2 ), - Vecd(0.0, -PH / 2 ), Vecd(-SL , -PH / 2 )}; + Vecd(-SL, -PH / 2), Vecd(-SL, PH / 2), Vecd(0.0, PH / 2), + Vecd(0.0, -PH / 2), Vecd(-SL, -PH / 2)}; // a beam shape std::vector beam_shape{ - Vecd(0.0, -PH / 2), Vecd(0.0, PH / 2), Vecd(PL, PH / 2), Vecd(PL, -PH / 2), Vecd(0.0, -PH / 2)}; + Vecd(0.0, -PH / 2), Vecd(0.0, PH / 2), Vecd(PL, PH / 2), Vecd(PL, -PH / 2), Vecd(0.0, -PH / 2)}; // Beam observer location StdVec observation_location = {Vecd(PL, 0.0)}; //---------------------------------------------------------------------- @@ -86,15 +84,13 @@ StdVec observation_location = {Vecd(PL, 0.0)}; //---------------------------------------------------------------------- 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); - } + 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 //---------------------------------------------------------------------- @@ -126,58 +122,59 @@ MultiPolygon createBeamConstrainShape() multi_polygon.addAPolygon(beam_shape, ShapeBooleanOps::sub); return multi_polygon; }; +//---------------------------------------------------------------------- +// calculate correction matrix B to keep the accuracy +//---------------------------------------------------------------------- -class EulerCorrectConfiguration : - public LocalDynamics, public GeneralDataDelegateInner +class AnisotropicCorrectConfiguration : public LocalDynamics, public GeneralDataDelegateInner { -public: - EulerCorrectConfiguration(BaseInnerRelation &inner_relation, int beta = 0, Real alpha = Real(0)) - : LocalDynamics(inner_relation.getSPHBody()), - GeneralDataDelegateInner(inner_relation), - beta_(beta), alpha_(alpha), - B_(*particles_->getVariableByName("CorrectionMatrix")), - pos_(particles_->pos_) - { - particles_->registerVariable(show_neighbor_, "ShowingNeighbor", Real(0.0)); - } - virtual ~EulerCorrectConfiguration() {}; -protected: - protected: + public: + AnisotropicCorrectConfiguration(BaseInnerRelation &inner_relation, int beta = 0, Real alpha = Real(0)) + : LocalDynamics(inner_relation.getSPHBody()), + GeneralDataDelegateInner(inner_relation), + beta_(beta), alpha_(alpha), + B_(*particles_->getVariableByName("CorrectionMatrix")), + pos_(particles_->pos_) + { + particles_->registerVariable(show_neighbor_, "ShowingNeighbor", Real(0.0)); + } + virtual ~AnisotropicCorrectConfiguration(){}; + + protected: + protected: int beta_; Real alpha_; - StdLargeVec& B_; - StdLargeVec& pos_; - StdLargeVec show_neighbor_; - - void interaction(size_t index_i, Real dt = 0.0) - { - Matd local_configuration = Eps * Matd::Identity(); - const 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]; - if (index_i == 67) - { - show_neighbor_[index_j] = 1.0; - - }; + StdLargeVec &B_; + StdLargeVec &pos_; + StdLargeVec show_neighbor_; - Vecd gradW_ijV_j = dW_ijV_j * e_ij; - Vecd r_ji = pos_[index_i] - pos_[index_j]; - local_configuration -= r_ji * gradW_ijV_j.transpose(); + void interaction(size_t index_i, Real dt = 0.0) + { + Matd local_configuration = Eps * Matd::Identity(); + const 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]; + if (index_i == 67) + { + show_neighbor_[index_j] = 1.0; + }; - } - B_[index_i] = local_configuration; - }; + Vecd gradW_ijV_j = dW_ijV_j * e_ij; + Vecd r_ji = pos_[index_i] - pos_[index_j]; + local_configuration -= r_ji * gradW_ijV_j.transpose(); + } + B_[index_i] = local_configuration; + }; void update(size_t index_i, Real dt) - { - Real det_sqr = pow(B_[index_i].determinant(), beta_); - Matd inverse = B_[index_i].inverse(); - B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr); - } + { + Real det_sqr = pow(B_[index_i].determinant(), beta_); + Matd inverse = B_[index_i].inverse(); + B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr); + } }; //------------------------------------------------------------------------------ @@ -185,125 +182,124 @@ class EulerCorrectConfiguration : //------------------------------------------------------------------------------ int main(int ac, char *av[]) { - //---------------------------------------------------------------------- - // Build up the environment of a SPHSystem with global controls. - //---------------------------------------------------------------------- - SPHSystem system(system_domain_bounds, resolution_ref_large); + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem with global controls. + //---------------------------------------------------------------------- + SPHSystem system(system_domain_bounds, resolution_ref_large); #ifdef BOOST_AVAILABLE - // handle command line arguments - system.handleCommandlineOptions(ac, av); + // handle command line arguments + system.handleCommandlineOptions(ac, av); #endif //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - SolidBody beam_body(system, makeShared("BeamBody")); - beam_body.sph_adaptation_->resetKernel>(scaling_vector); - beam_body.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); - beam_body.generateParticles(); + // Creating body, materials and particles. + //---------------------------------------------------------------------- + SolidBody beam_body(system, makeShared("BeamBody")); + beam_body.sph_adaptation_->resetKernel>(scaling_vector); + beam_body.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); + beam_body.generateParticles(); - ObserverBody beam_observer(system, "BeamObserver"); - beam_observer.sph_adaptation_->resetKernel>(scaling_vector); - 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. - //---------------------------------------------------------------------- - InnerRelation beam_body_inner(beam_body); - ContactRelation beam_observer_contact(beam_observer, {&beam_body}); - //----------------------------------------------------------------------------- - // this section define all numerical methods will be used in this case - //----------------------------------------------------------------------------- - SimpleDynamics beam_initial_velocity(beam_body); - // corrected strong configuration - InteractionWithUpdate beam_corrected_configuration(beam_body_inner); - beam_body.addBodyStateForRecording("ShowingNeighbor"); - // time step size calculation - ReduceDynamics computing_time_step_size(beam_body); - // stress relaxation for the beam - Dynamics1Level stress_relaxation_first_half(beam_body_inner); - Dynamics1Level stress_relaxation_second_half(beam_body_inner); - // clamping a solid body part. This is softer than a direct constraint - BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); + ObserverBody beam_observer(system, "BeamObserver"); + beam_observer.sph_adaptation_->resetKernel>(scaling_vector); + 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. + //---------------------------------------------------------------------- + InnerRelation beam_body_inner(beam_body); + ContactRelation beam_observer_contact(beam_observer, {&beam_body}); + //----------------------------------------------------------------------------- + // this section define all numerical methods will be used in this case + //----------------------------------------------------------------------------- + SimpleDynamics beam_initial_velocity(beam_body); + // corrected strong configuration + InteractionWithUpdate beam_corrected_configuration(beam_body_inner); + beam_body.addBodyStateForRecording("ShowingNeighbor"); + // time step size calculation + ReduceDynamics computing_time_step_size(beam_body); + // stress relaxation for the beam + Dynamics1Level stress_relaxation_first_half(beam_body_inner); + Dynamics1Level stress_relaxation_second_half(beam_body_inner); + // clamping a solid body part. + BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- - // outputs - //----------------------------------------------------------------------------- - IOEnvironment io_environment(system); + // outputs + //----------------------------------------------------------------------------- + IOEnvironment io_environment(system); BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); - RegressionTestEnsembleAverage> + RegressionTestEnsembleAverage> write_beam_tip_displacement("Position", io_environment, beam_observer_contact); - //---------------------------------------------------------------------- - // Setup computing and initial conditions. - //---------------------------------------------------------------------- + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- system.initializeSystemCellLinkedLists(); system.initializeSystemConfigurations(); beam_initial_velocity.exec(); beam_corrected_configuration.exec(); - //---------------------------------------------------------------------- - // Setup computing time-step controls. - //---------------------------------------------------------------------- - int ite = 0; - Real T0 = 1.0; - Real end_time = T0; - // time step size for output file - Real output_interval = 0.01 * T0; - Real Dt = 0.1 * output_interval; /**< Time period for data observing */ - Real dt = 0.0; // default acoustic time step sizes + //---------------------------------------------------------------------- + // Setup computing time-step controls. + //---------------------------------------------------------------------- + int ite = 0; + Real T0 = 1.0; + Real end_time = T0; + // time step size for output file + Real output_interval = 0.01 * T0; + Real Dt = 0.1 * output_interval; /**< Time period for data observing */ + Real dt = 0.0; // default acoustic time step sizes // statistics for computing time TickCount t1 = TickCount::now(); TimeInterval interval; - //----------------------------------------------------------------------------- - // from here the time stepping begins - //----------------------------------------------------------------------------- - write_beam_states.writeToFile(0); - write_beam_tip_displacement.writeToFile(0); + //----------------------------------------------------------------------------- + // from here the time stepping begins + //----------------------------------------------------------------------------- + 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; - // integrate time (loop) until the next output time - while (integration_time < output_interval) - { + // computation loop starts + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integration_time = 0.0; + // integrate time (loop) until the next output time + while (integration_time < output_interval) + { - Real relaxation_time = 0.0; - while (relaxation_time < Dt) - { - stress_relaxation_first_half.exec(dt); + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + stress_relaxation_first_half.exec(dt); constraint_beam_base.exec(); stress_relaxation_second_half.exec(dt); ite++; - dt = scaling_factor * computing_time_step_size.exec(); + dt = scaling_factor * computing_time_step_size.exec(); relaxation_time += dt; integration_time += dt; - GlobalStaticVariables::physical_time_ += dt; + GlobalStaticVariables::physical_time_ += dt; - if (ite % 100 == 0) - { - std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ << " dt: " - << dt << "\n"; - } - } - } + if (ite % 100 == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " dt: " + << dt << "\n"; + } + } + } - write_beam_tip_displacement.writeToFile(ite); + 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(); + } + TickCount t4 = TickCount::now(); TimeInterval tt; tt = t4 - t1 - interval; std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; - - if (system.GenerateRegressionData()) + + if (system.GenerateRegressionData()) { - // The lift force at the cylinder is very small and not important in this case. write_beam_tip_displacement.generateDataBase(Vec2d(1.0e-2, 1.0e-2), Vec2d(1.0e-2, 1.0e-2)); } else @@ -311,5 +307,5 @@ int main(int ac, char *av[]) write_beam_tip_displacement.testResult(); } - return 0; + return 0; } diff --git a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp index c90ff2545f..653f454f64 100644 --- a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp @@ -21,7 +21,7 @@ Real BW = resolution_ref * 4; // boundary width, at least three particles BoundingBox system_domain_bounds(Vec2d(-SL - BW, -PL / 2.0), Vec2d(PL + 3.0 * BW, PL / 2.0)); //---------------------------------------------------------------------- -// Material properties of the fluid. +// Material properties of the solid. //---------------------------------------------------------------------- Real rho0_s = 1.0e3; // reference density Real Youngs_modulus = 2.0e6; // reference Youngs modulus