diff --git a/tests/user_examples/test_2d_membrane/CMakeLists.txt b/tests/user_examples/test_2d_membrane/CMakeLists.txt index b4b2b11b8b..532de66c3e 100644 --- a/tests/user_examples/test_2d_membrane/CMakeLists.txt +++ b/tests/user_examples/test_2d_membrane/CMakeLists.txt @@ -7,6 +7,9 @@ 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) diff --git a/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_ensemble_averaged_mean_variance.xml b/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_ensemble_averaged_mean_variance.xml new file mode 100644 index 0000000000..2ab1860b4b --- /dev/null +++ b/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_ensemble_averaged_mean_variance.xml @@ -0,0 +1,257 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_runtimes.dat b/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_runtimes.dat new file mode 100644 index 0000000000..0320826438 --- /dev/null +++ b/tests/user_examples/test_2d_membrane/regression_test_tool/MembraneObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +6 +5 \ No newline at end of file diff --git a/tests/user_examples/test_2d_membrane/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_membrane/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..206155021f --- /dev/null +++ b/tests/user_examples/test_2d_membrane/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_membrane +""" + +case_name = "test_2d_membrane" +body_name = "MembraneObserver" +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 diff --git a/tests/user_examples/test_2d_membrane/2d_membrane.cpp b/tests/user_examples/test_2d_membrane/test_2d_membrane.cpp similarity index 56% rename from tests/user_examples/test_2d_membrane/2d_membrane.cpp rename to tests/user_examples/test_2d_membrane/test_2d_membrane.cpp index f8c01223d1..aedb4e2b36 100644 --- a/tests/user_examples/test_2d_membrane/2d_membrane.cpp +++ b/tests/user_examples/test_2d_membrane/test_2d_membrane.cpp @@ -1,12 +1,10 @@ -/* ---------------------------------------------------------------------------* - * SPHinXsys: 2D menmbrane 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. * - * ----------------------------------------------------------------------------*/ -#include "particle_momentum_dissipation.h" +/** + * @file membrane.cpp + * @brief fluid diffusion coupled with porous structure deformation + * @details This is the one of the test cases using fluid diffusion inside the porous elastic media model. * + In this case, the multi-time step algorithm is used. + * @author Xiaojing Tang and Xiangyu Hu + */ #include "particle_momentum_dissipation.hpp" #include "porous_media_dynamics.h" #include "porous_media_solid.h" @@ -17,89 +15,94 @@ using namespace SPH; // global parameters for the case //------------------------------------------------------------------------------ Real PL = 10.0; // membrane length -Real PH = 0.125; // membrane thickenss -Real BC = PL * 0.15; +Real PH = 0.125; // membrane thickness +Real BC = PL * 0.15; // half length of the saturation -int num = 8; +int y_num = 8; // reference particle spacing -Real resolution_ref = PH / num; +Real resolution_ref = PH / y_num; /** Domain bounds of the system. */ BoundingBox system_domain_bounds(Vec2d(-PL, -PL), Vec2d(2.0 * PL, PL)); //---------------------------------------------------------------------- -// Material properties of the fluid. +// Material properties of the solid. //---------------------------------------------------------------------- -Real rho_0 = 2.0; // // reference density non-dimensionlaize -Real poisson = 0.26316; -Real Youngs_modulus = 8.242e6; +Real rho_0 = 2.0; // reference solid density non-dimensional +Real poisson = 0.26316; /**< Poisson ratio. */ +Real Youngs_modulus = 8.242e6; /**< Youngs modulus. */ Real physical_viscosity = 5000.0; - -Real diffusivity_constant_ = 1.0e-4; -Real fulid_initial_density_ = 1.0; +Real saturation = 0.4; /**< solid porosity. */ +//---------------------------------------------------------------------- +// Material properties of the fluid. +//---------------------------------------------------------------------- +Real diffusivity_constant_ = 1.0e-4; // reference diffusion constant non-dimensional +Real fluid_initial_density_ = 1.0; // reference fluid density non-dimensional Real water_pressure_constant_ = 3.0e6; -Real saturation = 0.4; +// the criterion to represent the static state Real refer_density_energy = 0.5 * water_pressure_constant_; //---------------------------------------------------------------------- // Geometric shapes used in the system. //---------------------------------------------------------------------- // a membrane base shape - -std::vector beam_base_shape{ +std::vector membrane_base_shape{ Vecd(-resolution_ref * 3.0, -PH / 2.0), Vecd(-resolution_ref * 3.0, PH / 2.0), Vecd(0.0, PH / 2.0), Vecd(0.0, -PH / 2.0), Vecd(-resolution_ref * 3.0, -PH / 2.0)}; // a membrane shape -std::vector beam_shape{Vecd(0.0, -PH / 2.0), Vecd(0.0, PH / 2.0), +std::vector membrane_shape{Vecd(0.0, -PH / 2.0), Vecd(0.0, PH / 2.0), Vecd(PL, PH / 2.0), Vecd(PL, -PH / 2.0), Vecd(0.0, -PH / 2.0)}; // a membrane end shape -std::vector beam_end_shape{ +std::vector membrane_end_shape{ Vecd(PL, -PH / 2.0), Vecd(PL, PH / 2.0), Vecd(PL + 4.0 * resolution_ref, PH / 2.0), Vecd(PL + 4.0 * resolution_ref, -PH / 2.0), Vecd(PL, -PH / 2.0)}; // a membrane saturation shape -std::vector beam_saturation_shape{ +std::vector membrane_saturation_shape{ Vecd(PL / 2.0 - BC, 0.0), Vecd(PL / 2.0 - BC, PH / 2.0), Vecd(PL / 2.0 + BC, PH / 2.0), Vecd(PL / 2.0 + BC, 0.0), Vecd(PL / 2.0 - BC, 0.0)}; -// Beam observer location -StdVec observation_location = {Vecd(PL / 4.0, 0.0)}; +// membrane observer location +StdVec observation_location = {Vecd(PL / 2.0, 0.0)}; //---------------------------------------------------------------------- -// Define the beam body +// Define the membrane body //---------------------------------------------------------------------- -class Beam : public MultiPolygonShape +class Membrane : public MultiPolygonShape { public: - explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) + explicit Membrane(const std::string &shape_name) : MultiPolygonShape(shape_name) { - multi_polygon_.addAPolygon(beam_base_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(beam_end_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(membrane_base_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(membrane_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(membrane_end_shape, ShapeBooleanOps::add); } }; //---------------------------------------------------------------------- -// define the beam base which will be constrained. +// define the membrane base which will be constrained. //---------------------------------------------------------------------- -MultiPolygon createBeamConstrainShape() +MultiPolygon createMembraneConstrainShape() { MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_base_shape, ShapeBooleanOps::add); - multi_polygon.addAPolygon(beam_shape, ShapeBooleanOps::sub); - multi_polygon.addAPolygon(beam_end_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(membrane_base_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(membrane_shape, ShapeBooleanOps::sub); + multi_polygon.addAPolygon(membrane_end_shape, ShapeBooleanOps::add); return multi_polygon; }; +//---------------------------------------------------------------------- +// define the membrane part which the saturation condition is applied. +//---------------------------------------------------------------------- MultiPolygon createSaturationConstrainShape() { MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_saturation_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(membrane_saturation_shape, ShapeBooleanOps::add); return multi_polygon; }; @@ -116,7 +119,7 @@ class SaturationInitialCondition : public multi_species_continuum::PorousMediaSa void update(size_t index_i, Real dt = 0.0) { fluid_saturation_[index_i] = saturation; - fluid_mass_[index_i] = saturation * fulid_initial_density_ * Vol_update_[index_i]; + fluid_mass_[index_i] = saturation * fluid_initial_density_ * Vol_update_[index_i]; total_mass_[index_i] = rho_n_[index_i] * Vol_update_[index_i] + fluid_mass_[index_i]; }; }; @@ -136,14 +139,14 @@ int main(int ac, char *av[]) #endif //---------------------------------------------------------------------- // Creating body, materials and particles. //---------------------------------------------------------------------- - SolidBody beam_body(sph_system, makeShared("2dMembrane")); - beam_body.defineParticlesAndMaterial( - rho_0, Youngs_modulus, poisson, diffusivity_constant_, fulid_initial_density_, water_pressure_constant_); - beam_body.generateParticles(); - - ObserverBody beam_observer(sph_system, "MembraneObserver"); - beam_observer.defineAdaptationRatios(1.15, 2.0); - beam_observer.generateParticles(observation_location); + SolidBody membrane(sph_system, makeShared("MembraneBody")); + membrane.defineParticlesAndMaterial( + rho_0, Youngs_modulus, poisson, diffusivity_constant_, fluid_initial_density_, water_pressure_constant_); + membrane.generateParticles(); + + ObserverBody membrane_observer(sph_system, "MembraneObserver"); + membrane_observer.defineAdaptationRatios(1.15, 2.0); + membrane_observer.generateParticles(observation_location); //---------------------------------------------------------------------- // Define body relation map. // The contact map gives the topological connections between the bodies. @@ -152,66 +155,65 @@ int main(int ac, char *av[]) // At last, we define the complex relaxations by combining previous defined // inner and contact relations. //---------------------------------------------------------------------- - InnerRelation beam_body_inner(beam_body); - ContactRelation beam_observer_contact(beam_observer, {&beam_body}); + InnerRelation membrane_body_inner(membrane); + ContactRelation membrane_observer_contact(membrane_observer, {&membrane}); //----------------------------------------------------------------------------- // this section define all numerical methods will be used in this case //----------------------------------------------------------------------------- - // corrected strong configuration - InteractionWithUpdate beam_corrected_configuration(beam_body_inner); - // time step size calculation - ReduceDynamics computing_time_step_size(beam_body); - ReduceDynamics saturation_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); - Dynamics1Level saturation_relaxation(beam_body_inner); - - // clamping a solid body part. This is softer than a direct constraint - BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics clamp_constrain_beam_base(beam_base); - - BodyRegionByParticle beam_saturation(beam_body, makeShared(createSaturationConstrainShape())); - SimpleDynamics constrain_beam_saturation(beam_saturation); - - // need to be done - ReduceDynamics get_kinetic_energy(beam_body); - - /** Damping */ + InteractionWithUpdate membrane_corrected_configuration(membrane_body_inner); + // time step size calculation for solid stress relaxation + ReduceDynamics computing_time_step_size(membrane); + // time step size calculation for fluid diffusion relaxation + ReduceDynamics saturation_time_step_size(membrane); + + // stress relaxation for the membrane + Dynamics1Level stress_relaxation_first_half(membrane_body_inner); + Dynamics1Level stress_relaxation_second_half(membrane_body_inner); + // fluid diffusion relaxation inside the membrane + Dynamics1Level saturation_relaxation(membrane_body_inner); + + // clamping a solid body part using momentum constraint + BodyRegionByParticle membrane_base(membrane, makeShared(createMembraneConstrainShape())); + SimpleDynamics clamp_constrain_membrane_base(membrane_base); + + // applying the saturation boundary condition. + BodyRegionByParticle membrane_saturation(membrane, makeShared(createSaturationConstrainShape())); + SimpleDynamics constrain_membrane_saturation(membrane_saturation); + + //total mechanical energy to check the static state is achieved or not + ReduceDynamics get_kinetic_energy(membrane); + + /** Damping applied on momentum */ DampingWithRandomChoice>> - beam_damping(0.5, beam_body_inner, "TotalMomentum", physical_viscosity); + membrane_damping(0.5, membrane_body_inner, "TotalMomentum", physical_viscosity); + //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- IOEnvironment io_environment(sph_system); - BodyStatesRecordingToVtp write_beam_states(io_environment, sph_system.real_bodies_); - // note there is a line observation - - ObservedQuantityRecording - write_beam_tip_position("Position", io_environment, beam_observer_contact); - + BodyStatesRecordingToVtp write_membrane_states(io_environment, sph_system.real_bodies_); + RegressionTestEnsembleAverage> + write_membrane_tip_position("Position", io_environment, membrane_observer_contact); + //---------------------------------------------------------------------- // Setup computing and initial conditions. //---------------------------------------------------------------------- sph_system.initializeSystemCellLinkedLists(); sph_system.initializeSystemConfigurations(); - constrain_beam_saturation.exec(); - beam_corrected_configuration.exec(); + constrain_membrane_saturation.exec(); + membrane_corrected_configuration.exec(); + //---------------------------------------------------------------------- + // Setup computing time-step controls. + //---------------------------------------------------------------------- int ite = 0; int total_ite = 0; - GlobalStaticVariables::physical_time_ = 0.0; - //---------------------------------------------------------------------- - // Setup computing time-step controls. - //---------------------------------------------------------------------- - Real End_Time = 100.0; + // time for applying saturation condition Real setup_saturation_time_ = End_Time * 0.1; - // time step size for output file Real D_Time = End_Time / 100.0; Real dt = 0.0; // default acoustic time step sizes @@ -222,43 +224,44 @@ int main(int ac, char *av[]) //----------------------------------------------------------------------------- // from here the time stepping begins //----------------------------------------------------------------------------- - write_beam_states.writeToFile(0); - write_beam_tip_position.writeToFile(0); - + write_membrane_states.writeToFile(0); + write_membrane_tip_position.writeToFile(0); // computation loop starts while (GlobalStaticVariables::physical_time_ < End_Time) { Real integration_time = 0.0; - // integrate time (loop) until the next output time + // the outer loop for fluid diffusion, also integrate time (loop) until the next output time while (integration_time < D_Time) { Real Dt = saturation_time_step_size.exec(); if (GlobalStaticVariables::physical_time_ < setup_saturation_time_) { - constrain_beam_saturation.exec(); + constrain_membrane_saturation.exec(); } saturation_relaxation.exec(Dt); - int stress_ite = 0; + int stress_ite = 0; // solid stress relaxation times Real relaxation_time = 0.0; Real total_kinetic_energy = 1000.0; while (relaxation_time < Dt) { - if (total_kinetic_energy > (5e-9 * refer_density_energy)) // this is because we change the total mehanical energy calculation + //the inner loop for solid stress relaxation + if (total_kinetic_energy > (5e-9 * refer_density_energy)) { stress_relaxation_first_half.exec(dt); - clamp_constrain_beam_base.exec(); - beam_damping.exec(dt); - clamp_constrain_beam_base.exec(); + clamp_constrain_membrane_base.exec(); + membrane_damping.exec(dt); + clamp_constrain_membrane_base.exec(); stress_relaxation_second_half.exec(dt); + //calculate the total kinetic energy to check the static state total_kinetic_energy = get_kinetic_energy.exec(); ite++; stress_ite++; dt = SMIN(computing_time_step_size.exec(), Dt); - if (ite % 1000 == 0) + if (ite % 10000 == 0) { std::cout << "N=" << ite << " Time: " << GlobalStaticVariables::physical_time_ << " Dt:" << Dt << " dt: " @@ -271,15 +274,14 @@ int main(int ac, char *av[]) integration_time += dt; GlobalStaticVariables::physical_time_ += dt; } - + write_membrane_tip_position.writeToFile(ite); std::cout << "One Diffusion finishes " - << "total_kinetic_energy = " << total_kinetic_energy - << " stress_ite = " << stress_ite << std::endl; + << "total_kinetic_energy = " << total_kinetic_energy + << " stress_ite = " << stress_ite << std::endl; } TickCount t2 = TickCount::now(); - write_beam_states.writeToFile(ite); - write_beam_tip_position.writeToFile(ite); + write_membrane_states.writeToFile(ite); TickCount t3 = TickCount::now(); interval += t3 - t2; @@ -293,6 +295,13 @@ int main(int ac, char *av[]) std::cout << "Total iterations computation: " << GlobalStaticVariables::physical_time_ / dt << " Total iterations: " << total_ite << std::endl; - + if (sph_system.GenerateRegressionData()) + { + write_membrane_tip_position.generateDataBase(Vec2d(1.0e-2, 1.0e-2), Vec2d(1.0e-2, 1.0e-2)); + } + else + { + write_membrane_tip_position.testResult(); + } return 0; } diff --git a/tests/user_examples/test_2d_stretching/CMakeLists.txt b/tests/user_examples/test_2d_stretching/CMakeLists.txt index 63dd6a3ac8..6c3edab17e 100644 --- a/tests/user_examples/test_2d_stretching/CMakeLists.txt +++ b/tests/user_examples/test_2d_stretching/CMakeLists.txt @@ -11,14 +11,14 @@ 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}) aux_source_directory(. DIR_SRCS) ADD_EXECUTABLE(${PROJECT_NAME} ${DIR_SRCS}) - target_link_libraries(${PROJECT_NAME} sphinxsys_2d GTest::gtest GTest::gtest_main) target_link_libraries(${PROJECT_NAME} extra_sources_2d) - set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") diff --git a/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_dtwdistance.xml b/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_dtwdistance.xml new file mode 100644 index 0000000000..24b979334b --- /dev/null +++ b/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_dtwdistance.xml @@ -0,0 +1,4 @@ + + + + diff --git a/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_runtimes.dat b/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_runtimes.dat new file mode 100644 index 0000000000..34de7daedb --- /dev/null +++ b/tests/user_examples/test_2d_stretching/regression_test_tool/BeamObserver_Position_runtimes.dat @@ -0,0 +1,3 @@ +true +7 +4 \ No newline at end of file diff --git a/tests/user_examples/test_2d_stretching/regression_test_tool/regression_test_tool.py b/tests/user_examples/test_2d_stretching/regression_test_tool/regression_test_tool.py new file mode 100644 index 0000000000..161bf99136 --- /dev/null +++ b/tests/user_examples/test_2d_stretching/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_2d_stretching +""" + +case_name = "test_2d_stretching" +body_name = "BeamObserver" +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/user_examples/test_2d_stretching/stretching.cpp b/tests/user_examples/test_2d_stretching/stretching.cpp deleted file mode 100644 index fa398942ea..0000000000 --- a/tests/user_examples/test_2d_stretching/stretching.cpp +++ /dev/null @@ -1,370 +0,0 @@ -/* ---------------------------------------------------------------------------* - * SPHinXsys: 2D oscillation 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. * - * ----------------------------------------------------------------------------*/ -#include "inelastic_solid_hardening.h" -#include "sphinxsys.h" - -using namespace SPH; -//------------------------------------------------------------------------------ -// global parameters for the case -//------------------------------------------------------------------------------ -Real PL = 0.05334; // beam length -Real PH = 0.012826; // for thick plate; - -// reference particle spacing -Real resolution_ref = PH / 30; -Real BW = resolution_ref * 4.0; // boundary width, at least three particles - -/** Domain bounds of the system. */ -BoundingBox system_domain_bounds(Vec2d(-PL / 2.0, -PL / 2.0), - Vec2d(2.0 * PL, PL / 2.0)); -// two dimensional should be circle smooth between two parts. -//---------------------------------------------------------------------- -Real rho0_s = 7850.0; /**< Reference density. */ -Real Shear_modulus = 80.1938e9; /**< Poisson ratio. */ -Real Bulk_modulus = 164.21e9; - -Real poisson = (3.0 * Bulk_modulus - 2.0 * Shear_modulus) / (6.0 * Bulk_modulus + 2.0 * Shear_modulus); /**< Poisson ratio. */ -Real Youngs_modulus = (9.0 * Shear_modulus * Bulk_modulus) / (3.0 * Bulk_modulus + Shear_modulus); - -Real yield_stress = 0.45e9; -Real hardening_modulus = 1.2924e8; -Real saturation_flow_stress = 7.15e8; -Real saturation_exponent = 16.93; - -Real physical_viscosity = 1.0e4; -Real refer_energy = 0.5 * 8000 * 0.01; // 40 - -Vecd norm_(1.0, 0.0); -Vecd upper_face_point_(0.02 + 3.0 * resolution_ref, 0.0); -Vecd lower_face_point_(0.02, 0.0); - -Vecd norm_4(1.0, 0.0); -Vecd upper_face_point_4(0.04 + 3.0 * resolution_ref, 0.0); -Vecd lower_face_point_4(0.04, 0.0); - -// Beam observer location -StdVec observation_location = {Vecd(PL / 2.0, PH / 2.0 - PH * 0.01)}; - -//---------------------------------------------------------------------- -// Geometric shapes used in the system. -//---------------------------------------------------------------------- - -std::vector beam_left_stretch_shape{ - Vecd(-BW, -PH / 2), Vecd(-BW, PH / 2), Vecd(0.0, PH / 2), - Vecd(0.0, -PH / 2), Vecd(-BW, -PH / 2)}; - -// a beam shape -std::vector beam_shape{ - Vecd(0.0, -PH / 2), Vecd(0.0, PH / 2), - Vecd(PL / 2.0, PH / 2 - PH * 0.01), - Vecd(PL, PH / 2), Vecd(PL, -PH / 2), - Vecd(PL / 2.0, -PH / 2 + PH * 0.01), - Vecd(0.0, -PH / 2)}; - -std::vector beam_right_stretch_shape{ - Vecd(PL, -PH / 2), Vecd(PL, PH / 2), Vecd(PL + BW, PH / 2), - Vecd(PL + BW, -PH / 2), Vecd(PL, -PH / 2)}; - -//---------------------------------------------------------------------- -// Define the beam body -//---------------------------------------------------------------------- -class Beam : public MultiPolygonShape -{ - public: - explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - multi_polygon_.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); - } -}; - -class LeftStretchSolidBodyRegion : public solid_dynamics::BaseMotionConstraint -{ - public: - // TODO: use only body part as argment since body can be referred from it already - LeftStretchSolidBodyRegion(BodyPartByParticle &body_part) - : solid_dynamics::BaseMotionConstraint(body_part), - vel_(particles_->vel_), pos_(particles_->pos_){}; - - virtual ~LeftStretchSolidBodyRegion(){}; - - protected: - StdLargeVec &vel_; - StdLargeVec &pos_; - virtual void update(size_t index_i, Real Dt = 0.0) - { - pos_[index_i][0] -= 0.5e-4 * Dt; - }; -}; - -class RightStretchSolidBodyRegion : public solid_dynamics::BaseMotionConstraint -{ - public: - // TODO: use only body part as argment since body can be referred from it already - RightStretchSolidBodyRegion(BodyPartByParticle &body_part) - : solid_dynamics::BaseMotionConstraint(body_part), - vel_(particles_->vel_), pos_(particles_->pos_){}; - - virtual ~RightStretchSolidBodyRegion(){}; - - protected: - StdLargeVec &vel_; - StdLargeVec &pos_; - virtual void update(size_t index_i, Real Dt = 0.0) - { - pos_[index_i][0] += 0.5e-4 * Dt; - }; -}; - -MultiPolygon createBeamRightStretchShape() -{ - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); - return multi_polygon; -}; - -MultiPolygon createBeamLeftStretchShape() -{ - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); - return multi_polygon; -}; - -MultiPolygon createConstrainBeamShape() -{ - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); - multi_polygon.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); - - return multi_polygon; -}; - -class ConstrainXVelocity : public solid_dynamics::BaseMotionConstraint -{ - public: - // TODO: use only body part as argment since body can be referred from it already - ConstrainXVelocity(BodyPartByParticle &body_part) - : solid_dynamics::BaseMotionConstraint(body_part), - vel_(particles_->vel_), pos_(particles_->pos_){}; - - virtual ~ConstrainXVelocity(){}; - - protected: - StdLargeVec &vel_; - StdLargeVec &pos_; - virtual void update(size_t index_i, Real dt = 0.0) - { - vel_[index_i] = Vecd(0.0, vel_[index_i][1]); - }; -}; - -//------------------------------------------------------------------------------ -// 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); - /** Tag for running particle relaxation for the initially body-fitted distribution */ - system.setRunParticleRelaxation(false); - /** Tag for starting with relaxed body-fitted particles distribution */ - system.setReloadParticles(false); - system.handleCommandlineOptions(ac, av); - IOEnvironment io_environment(system); - - //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - SolidBody beam_body(system, makeShared("StretchingBody")); - beam_body.defineBodyLevelSetShape(); - beam_body.defineParticlesAndMaterial( - rho0_s, Youngs_modulus, poisson, yield_stress, hardening_modulus, saturation_flow_stress, saturation_exponent); - - (!system.RunParticleRelaxation() && system.ReloadParticles()) - ? beam_body.generateParticles(io_environment, beam_body.getName()) - : beam_body.generateParticles(); - - ObserverBody beam_observer(system, "BeamObserver"); - 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. - //---------------------------------------------------------------------- - if (system.RunParticleRelaxation()) - { - //---------------------------------------------------------------------- - // Define body relation map used for particle relaxation. - //---------------------------------------------------------------------- - InnerRelation beam_body_inner(beam_body); - //---------------------------------------------------------------------- - // Define the methods for particle relaxation. - //---------------------------------------------------------------------- - SimpleDynamics beam_body_random_particles(beam_body); - relax_dynamics::RelaxationStepInner beam_body_relaxation_step_inner(beam_body_inner); - //---------------------------------------------------------------------- - // Output for particle relaxation. - //---------------------------------------------------------------------- - BodyStatesRecordingToVtp write_ball_state(io_environment, system.real_bodies_); - ReloadParticleIO write_particle_reload_files(io_environment, {&beam_body}); - //---------------------------------------------------------------------- - // Particle relaxation starts here. - //---------------------------------------------------------------------- - beam_body_random_particles.exec(0.25); - write_ball_state.writeToFile(0); - //---------------------------------------------------------------------- - // From here iteration for particle relaxation begins. - //---------------------------------------------------------------------- - int ite = 0; - int relax_step = 1000; - while (ite < relax_step) - { - beam_body_relaxation_step_inner.exec(); - ite += 1; - if (ite % 100 == 0) - { - std::cout << std::fixed << std::setprecision(9) << "Relaxation steps N = " << ite << "\n"; - write_ball_state.writeToFile(ite); - } - } - std::cout << "The physics relaxation process of particles finish !" << std::endl; - write_particle_reload_files.writeToFile(0); - return 0; - } - - 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 - //----------------------------------------------------------------------------- - - // corrected strong configuration - InteractionWithUpdate beam_corrected_configuration(beam_body_inner); - - // 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); - ReduceDynamics get_kinetic_energy(beam_body); - - BodyRegionByParticle beam_left_stretch(beam_body, makeShared(createBeamLeftStretchShape())); - SimpleDynamics stretch_beam_left_end(beam_left_stretch); - BodyRegionByParticle beam_right_stretch(beam_body, makeShared(createBeamRightStretchShape())); - SimpleDynamics stretch_beam_right_end(beam_right_stretch); - BodyRegionByParticle beam_constrain(beam_body, makeShared(createConstrainBeamShape())); - SimpleDynamics constrain_beam_end(beam_constrain); - - InteractionDynamics beam_deformation_gradient_tensor(beam_body_inner); - DampingWithRandomChoice>> - damping(0.5, beam_body_inner, "Velocity", physical_viscosity); - - //----------------------------------------------------------------------------- - // outputs - //----------------------------------------------------------------------------- - - BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); - ReducedQuantityRecording> - write_total_mechanical_energy(io_environment, beam_body); - //---------------------------------------------------------------------- - // Setup computing and initial conditions. - //---------------------------------------------------------------------- - system.initializeSystemCellLinkedLists(); - system.initializeSystemConfigurations(); - beam_corrected_configuration.exec(); - - //---------------------------------------------------------------------- - // Setup computing time-step controls. - //---------------------------------------------------------------------- - int ite = 0; - int Dt_ite = 0; - Real End_Time = 100.0; - - Real D_Time = End_Time / 100.0; // time step size for output file - Real Dt = End_Time / 10000.0; /**< Time period for stretching */ - 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_total_mechanical_energy.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 < D_Time) - { - Real relaxation_time = 0.0; - stretch_beam_left_end.exec(Dt); - stretch_beam_right_end.exec(Dt); - beam_deformation_gradient_tensor.exec(Dt); - Dt_ite++; - - int stress_ite = 0; - Real refer_total_kinetic_energy = 10000.0; - while (relaxation_time < Dt) - { - if (refer_total_kinetic_energy > 0.005) - { - stress_relaxation_first_half.exec(dt); - constrain_beam_end.exec(dt); - damping.exec(Dt); - constrain_beam_end.exec(dt); - stress_relaxation_second_half.exec(dt); - - refer_total_kinetic_energy = get_kinetic_energy.exec() / refer_energy; - - ite++; - stress_ite++; - - dt = computing_time_step_size.exec(); - if (ite % 500 == 0) - { - std::cout << "N=" << ite << " Time: " - << GlobalStaticVariables::physical_time_ - << " Dt: " << Dt << " dt: " << dt - << " Dt:dt = " << Dt / dt << "\n"; - } - } - relaxation_time += dt; - integration_time += dt; - GlobalStaticVariables::physical_time_ += dt; - } - - std::cout << "refer_total_kinetic_energy " << refer_total_kinetic_energy - << " stress_ite " << stress_ite << std::endl; - } - - 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." - << " Iterations: " << ite << std::endl; - std::cout << "Total iterations computation: " << GlobalStaticVariables::physical_time_ / dt - << std::endl; - return 0; -} diff --git a/tests/user_examples/test_2d_stretching/test_2d_stretching.cpp b/tests/user_examples/test_2d_stretching/test_2d_stretching.cpp new file mode 100644 index 0000000000..8605aa2bd3 --- /dev/null +++ b/tests/user_examples/test_2d_stretching/test_2d_stretching.cpp @@ -0,0 +1,391 @@ +/** + * @file stretching.cpp + * @brief Stretching a tensile from the two surfaces with displacement boundary condition. + * @details This is the one of the test cases using nonlinear hardening plastic model. * + In this case, the multi-time step algorithm is used. + * @author Xiaojing Tang and Xiangyu Hu + */ + +#include "sphinxsys.h" +#include "inelastic_solid_hardening.h" +using namespace SPH; +//------------------------------------------------------------------------------ +//global parameters for the case +//------------------------------------------------------------------------------ +Real PL = 0.05334; // beam length +Real PH = 0.012826; // beam width + +//reference particle spacing +Real resolution_ref = PH / 30; +Real BW = resolution_ref * 4.0 ; //boundary width, at least three particles + +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec2d(-PL / 2.0, -PL / 2.0), + Vec2d(2.0*PL , PL / 2.0)); + +//---------------------------------------------------------------------- +// Material properties of the solid. +//---------------------------------------------------------------------- +Real rho0_s = 7850.0; /**< Reference density. */ +Real Shear_modulus = 80.1938e9; /**< Shear modulus. */ +Real Bulk_modulus = 164.21e9; /**< Bulk modulus. */ + +Real poisson = (3.0 * Bulk_modulus -2.0 * Shear_modulus)/(6.0 * Bulk_modulus + 2.0 * Shear_modulus); /**< Poisson ratio. */ +Real Youngs_modulus = (9.0 * Shear_modulus * Bulk_modulus)/(3.0* Bulk_modulus + Shear_modulus); /**< Youngs modulus. */ +Real physical_viscosity = 1.0e4; + +/** plastic parameters. */ +Real yield_stress = 0.45e9; +Real hardening_modulus = 1.2924e8; +Real saturation_flow_stress = 7.15e8; +Real saturation_exponent = 16.93; + +// the criterion to represent the static state +Real refer_energy = 0.5*8000*0.01; +//Beam observer location +StdVec observation_location = { Vecd( PL / 2.0 , PH / 2.0 - PH * 0.01) }; + +//---------------------------------------------------------------------- +// Geometric shapes used in the system. +//---------------------------------------------------------------------- + +//a beam left stretch shape +std::vector beam_left_stretch_shape{ + Vecd(-BW , -PH / 2 ), Vecd(-BW, PH / 2 ), Vecd(0.0, PH / 2 ), + Vecd(0.0, -PH / 2 ), Vecd(-BW, -PH / 2 )}; + +//a beam shape +std::vector beam_shape{ + Vecd(0.0, -PH / 2), Vecd(0.0, PH / 2), + Vecd(PL / 2.0 , PH / 2 - PH * 0.01), + Vecd(PL, PH / 2), Vecd(PL , -PH / 2), + Vecd(PL / 2.0 , -PH / 2 + PH * 0.01), + Vecd(0.0, -PH / 2) }; + +//a beam right stretch shape +std::vector beam_right_stretch_shape{ + Vecd(PL, -PH / 2 ), Vecd(PL, PH / 2 ), Vecd(PL +BW , PH / 2 ), + Vecd(PL+BW, -PH / 2 ), Vecd(PL, -PH / 2 ) }; + +//---------------------------------------------------------------------- +// Define the beam body +//---------------------------------------------------------------------- +class Beam : public MultiPolygonShape +{ +public: + explicit Beam(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + multi_polygon_.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(beam_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); + } +}; + +//---------------------------------------------------------------------- +// application boundary displacement condition +//---------------------------------------------------------------------- +class LeftStretchSolidBodyRegion : public solid_dynamics::BaseMotionConstraint +{ +public: + // TODO: use only body part as argment since body can be referred from it already + LeftStretchSolidBodyRegion(BodyPartByParticle &body_part) + :solid_dynamics::BaseMotionConstraint(body_part), + vel_(particles_->vel_), pos_(particles_->pos_){}; + + virtual ~LeftStretchSolidBodyRegion() {}; +protected: + StdLargeVec &vel_; + StdLargeVec &pos_; + virtual void update(size_t index_i, Real Dt = 0.0) + { + pos_[index_i][0] -= 0.5e-4 * Dt; + }; +}; + +//---------------------------------------------------------------------- +// application boundary displacement condition +//---------------------------------------------------------------------- +class RightStretchSolidBodyRegion : public solid_dynamics::BaseMotionConstraint +{ +public: + // TODO: use only body part as argment since body can be referred from it already + RightStretchSolidBodyRegion(BodyPartByParticle &body_part) + :solid_dynamics::BaseMotionConstraint(body_part), + vel_(particles_->vel_), pos_(particles_->pos_){}; + + virtual ~RightStretchSolidBodyRegion() {}; + +protected: + StdLargeVec &vel_; + StdLargeVec &pos_; + virtual void update(size_t index_i, Real Dt = 0.0) + { + pos_[index_i][0] += 0.5e-4 * Dt; + }; +}; + +//---------------------------------------------------------------------- +// define the beam part which will be stretched +//---------------------------------------------------------------------- +MultiPolygon createBeamRightStretchShape() +{ + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); + return multi_polygon; +}; + +//---------------------------------------------------------------------- +// define the beam part which will be stretched +//---------------------------------------------------------------------- +MultiPolygon createBeamLeftStretchShape() +{ + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); + return multi_polygon; +}; + + +MultiPolygon createConstrainBeamShape() +{ + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(beam_left_stretch_shape, ShapeBooleanOps::add); + multi_polygon.addAPolygon(beam_right_stretch_shape, ShapeBooleanOps::add); + return multi_polygon; +}; + +//---------------------------------------------------------------------- +// application constrain condition +//---------------------------------------------------------------------- +class ConstrainXVelocity : public solid_dynamics::BaseMotionConstraint +{ +public: + // TODO: use only body part as argment since body can be referred from it already + ConstrainXVelocity(BodyPartByParticle &body_part) + :solid_dynamics::BaseMotionConstraint(body_part), + vel_(particles_->vel_), pos_(particles_->pos_){}; + + virtual ~ConstrainXVelocity() {}; +protected: + StdLargeVec &vel_; + StdLargeVec &pos_; + virtual void update(size_t index_i, Real dt = 0.0) + { + vel_[index_i] = Vecd(0.0, vel_[index_i][1]); + }; +}; + +//------------------------------------------------------------------------------ +//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); + /** Tag for running particle relaxation for the initially body-fitted distribution */ + system.setRunParticleRelaxation(false); + /** Tag for starting with relaxed body-fitted particles distribution */ + system.setReloadParticles(true); + system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(system); + + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + SolidBody beam_body(system, makeShared("StretchingBody")); + beam_body.defineBodyLevelSetShape(); + beam_body.defineParticlesAndMaterial( + rho0_s, Youngs_modulus, poisson, yield_stress, hardening_modulus, saturation_flow_stress, saturation_exponent); + (!system.RunParticleRelaxation() && system.ReloadParticles()) + ? beam_body.generateParticles(io_environment, beam_body.getName()) + : beam_body.generateParticles(); + + ObserverBody beam_observer(system, "BeamObserver"); + beam_observer.generateParticles(observation_location); + + if (system.RunParticleRelaxation()) + { + //---------------------------------------------------------------------- + // Define body relation map used for particle relaxation. + //---------------------------------------------------------------------- + InnerRelation beam_body_inner(beam_body); + //---------------------------------------------------------------------- + // Define the methods for particle relaxation. + //---------------------------------------------------------------------- + SimpleDynamics beam_body_random_particles(beam_body); + relax_dynamics::RelaxationStepInner beam_body_relaxation_step_inner(beam_body_inner); + //---------------------------------------------------------------------- + // Output for particle relaxation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp write_ball_state(io_environment, system.real_bodies_); + ReloadParticleIO write_particle_reload_files(io_environment, { &beam_body }); + //---------------------------------------------------------------------- + // Particle relaxation starts here. + //---------------------------------------------------------------------- + beam_body_random_particles.exec(0.25); + write_ball_state.writeToFile(0); + //---------------------------------------------------------------------- + // From here iteration for particle relaxation begins. + //---------------------------------------------------------------------- + int ite = 0; + int relax_step = 1000; + while (ite < relax_step) + { + beam_body_relaxation_step_inner.exec(); + ite += 1; + if (ite % 100 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps N = " << ite << "\n"; + write_ball_state.writeToFile(ite); + } + } + std::cout << "The physics relaxation process of particles finish !" << std::endl; + write_particle_reload_files.writeToFile(0); + return 0; + } + + //---------------------------------------------------------------------- + // 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 + //----------------------------------------------------------------------------- + //corrected strong configuration + InteractionWithUpdate beam_corrected_configuration(beam_body_inner); + //time step size calculation + ReduceDynamics computing_time_step_size(beam_body); + //stress relaxation for the body using plastic integration + Dynamics1Level stress_relaxation_first_half(beam_body_inner); + Dynamics1Level stress_relaxation_second_half(beam_body_inner); + //total mechanical energy to check the static state is achieved or not + ReduceDynamics get_kinetic_energy(beam_body); + + //boundary conditions + BodyRegionByParticle beam_left_stretch(beam_body, makeShared(createBeamLeftStretchShape())); + SimpleDynamics stretch_beam_left_end(beam_left_stretch); + BodyRegionByParticle beam_right_stretch(beam_body, makeShared(createBeamRightStretchShape())); + SimpleDynamics stretch_beam_right_end(beam_right_stretch); + BodyRegionByParticle beam_constrain(beam_body, makeShared(createConstrainBeamShape())); + SimpleDynamics constrain_beam_end(beam_constrain); + + InteractionDynamics beam_deformation_gradient_tensor(beam_body_inner); + //damping term + DampingWithRandomChoice>> + damping(0.5, beam_body_inner, "Velocity", physical_viscosity); + + //----------------------------------------------------------------------------- + //outputs + //----------------------------------------------------------------------------- + BodyStatesRecordingToVtp write_beam_states(io_environment, system.real_bodies_); + ReducedQuantityRecording> + write_total_mechanical_energy(io_environment, beam_body); + RegressionTestDynamicTimeWarping> + write_beam_tip_displacement("Position", io_environment, beam_observer_contact); + + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + beam_corrected_configuration.exec(); + + //---------------------------------------------------------------------- + // Setup computing time-step controls. + //---------------------------------------------------------------------- + int ite = 0; + int Dt_ite = 0; // stretching times + Real End_Time = 100.0; + + Real D_Time = End_Time /100.0; //time step size for output file + Real Dt = End_Time /10000.0; /**< Time period for stretching */ + 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); + //computation loop starts + while (GlobalStaticVariables::physical_time_ < End_Time) + { + Real integration_time = 0.0; + //the outer loop for stretching, also integrate time (loop) until the next output time + while (integration_time < D_Time) + { + Real relaxation_time = 0.0; + stretch_beam_left_end.exec(Dt); + stretch_beam_right_end.exec(Dt); + beam_deformation_gradient_tensor.exec(Dt); + Dt_ite ++; + + + int stress_ite = 0; // solid stress relaxation times + Real refer_total_kinetic_energy = 1.0; + while (relaxation_time < Dt) + { + //the inner loop for solid stress relaxation + if (refer_total_kinetic_energy > 0.005) + { + stress_relaxation_first_half.exec(dt); + constrain_beam_end.exec(dt); + damping.exec(Dt); + constrain_beam_end.exec(dt); + stress_relaxation_second_half.exec(dt); + //calculate the total kinetic energy to check the static state + refer_total_kinetic_energy = get_kinetic_energy.exec() / refer_energy; + + ite++; + stress_ite++; + + dt = computing_time_step_size.exec(); + if (ite % 1000 == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ + << " Dt: " << Dt << " dt: " << dt + << " Dt:dt = " << Dt / dt << + "\n"; + } + } + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + } + + + } + + 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." + << " Iterations: " << ite << std::endl; + + + if (system.GenerateRegressionData()) + { + write_beam_tip_displacement.generateDataBase(1.0e-3); + } + else + { + write_beam_tip_displacement.testResult(); + } + + return 0; +}