diff --git a/README.md b/README.md index 7f76a689b1..3c59fdda19 100644 --- a/README.md +++ b/README.md @@ -35,6 +35,8 @@ Here, we present several short examples in flow, solid dynamics, fluid structure + + ## Journal publications diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp index 748d4c6f0d..03cd7ae672 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp +++ b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp @@ -1,291 +1,292 @@ - -#include "sphinxsys.h" +/** + * @file 2d_flow_stream_around_fish.cpp + * @brief fish swimming driven by active muscles + * @author Yaru Ren and Xiangyu Hu + */ #include "2d_flow_stream_around_fish.h" #include "active_model.h" +#include "sphinxsys.h" using namespace SPH; int main(int ac, char *av[]) { - //---------------------------------------------------------------------- - // Build up the environment of a SPHSystem. - //---------------------------------------------------------------------- - SPHSystem system(system_domain_bounds, particle_spacing_ref); - /** Tag for run particle relaxation for the initial body fitted distribution. */ - system.setRunParticleRelaxation(true); - /** Tag for computation start with relaxed body fitted particles distribution. */ - system.setReloadParticles(false); - //handle command line arguments + //---------------------------------------------------------------------- + // Build up the environment of a SPHSystem. + //---------------------------------------------------------------------- + SPHSystem system(system_domain_bounds, particle_spacing_ref); + /** Tag for run particle relaxation for the initial body fitted distribution. */ + system.setRunParticleRelaxation(true); + /** Tag for computation start with relaxed body fitted particles distribution. */ + system.setReloadParticles(false); + // handle command line arguments #ifdef BOOST_AVAILABLE - system.handleCommandlineOptions(ac, av); - IOEnvironment io_environment(system); + system.handleCommandlineOptions(ac, av); + IOEnvironment io_environment(system); #endif - //---------------------------------------------------------------------- - // Creating body, materials and particles. - //---------------------------------------------------------------------- - //---------------------------------------------------------------------- - FluidBody water_block(system, makeShared("WaterBody")); - water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); - water_block.generateParticles(); - /** - * @brief Particles and body creation for fish. - */ - SolidBody fish_body(system, makeShared("FishBody")); - fish_body.defineAdaptationRatios(1.15, 2.0); - fish_body.defineBodyLevelSetShape()->writeLevelSet(io_environment); - fish_body.defineParticlesAndMaterial(); - //fish_body.defineParticlesAndMaterial>(rho0_s, bulk_modulus1, fiber_direction, sheet_direction, a01, b0); - // Using relaxed particle distribution if needed - (!system.RunParticleRelaxation() && system.ReloadParticles()) - ? fish_body.generateParticles(io_environment, fish_body.getName()) - : fish_body.generateParticles(); - //---------------------------------------------------------------------- - // 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 fish_inner(fish_body); - InnerRelation water_block_inner(water_block); - ComplexRelation water_block_complex(water_block, { &fish_body }); - ContactRelation fish_contact(fish_body, { &water_block }); - //---------------------------------------------------------------------- - // Run particle relaxation for body-fitted distribution if chosen. - //---------------------------------------------------------------------- - if (system.RunParticleRelaxation()) - { - /** - * @brief Methods used for particle relaxation. - */ - /** Random reset the insert body particle position. */ - SimpleDynamics random_fish_body_particles(fish_body); - /** Write the body state to Vtp file. */ - BodyStatesRecordingToVtp write_fish_body(io_environment, fish_body); - /** Write the particle reload files. */ - ReloadParticleIO write_particle_reload_files(io_environment, { &fish_body }); - - /** A Physics relaxation step. */ - relax_dynamics::RelaxationStepInner relaxation_step_inner(fish_inner); - /** - * @brief Particle relaxation starts here. - */ - random_fish_body_particles.exec(0.25); - relaxation_step_inner.SurfaceBounding().exec(); - write_fish_body.writeToFile(); - - /** relax particles of the insert body. */ - int ite_p = 0; - while (ite_p < 1000) - { - relaxation_step_inner.exec(); - ite_p += 1; - if (ite_p % 200 == 0) - { - std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the inserted body N = " << ite_p << "\n"; - write_fish_body.writeToFile(ite_p); - } - } - std::cout << "The physics relaxation process of inserted body finish !" << std::endl; + //---------------------------------------------------------------------- + // Creating body, materials and particles. + //---------------------------------------------------------------------- + //---------------------------------------------------------------------- + FluidBody water_block(system, makeShared("WaterBody")); + water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); + water_block.generateParticles(); + /** + * @brief Particles and body creation for fish. + */ + SolidBody fish_body(system, makeShared("FishBody")); + fish_body.defineAdaptationRatios(1.15, 2.0); + fish_body.defineBodyLevelSetShape()->writeLevelSet(io_environment); + fish_body.defineParticlesAndMaterial(); + // Using relaxed particle distribution if needed + (!system.RunParticleRelaxation() && system.ReloadParticles()) + ? fish_body.generateParticles(io_environment, fish_body.getName()) + : fish_body.generateParticles(); + //---------------------------------------------------------------------- + // 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 fish_inner(fish_body); + InnerRelation water_block_inner(water_block); + ComplexRelation water_block_complex(water_block, {&fish_body}); + ContactRelation fish_contact(fish_body, {&water_block}); + //---------------------------------------------------------------------- + // Run particle relaxation for body-fitted distribution if chosen. + //---------------------------------------------------------------------- + if (system.RunParticleRelaxation()) + { + /** + * @brief Methods used for particle relaxation. + */ + /** Random reset the insert body particle position. */ + SimpleDynamics random_fish_body_particles(fish_body); + /** Write the body state to Vtp file. */ + BodyStatesRecordingToVtp write_fish_body(io_environment, fish_body); + /** Write the particle reload files. */ + ReloadParticleIO write_particle_reload_files(io_environment, {&fish_body}); - /** Output results. */ - write_particle_reload_files.writeToFile(); - return 0; - } - //---------------------------------------------------------------------- - // Define the main numerical methods used in the simulation. - // Note that there may be data dependence on the constructors of these methods. - //---------------------------------------------------------------------- - SimpleDynamics initialize_a_fluid_step(water_block, makeShared(Vec2d::Zero())); - BodyAlignedBoxByParticle emitter( - water_block, makeShared(Transform(Vec2d(emitter_translation)), emitter_halfsize)); - SimpleDynamics emitter_inflow_injection(emitter, 10, 0); - /** Emitter buffer inflow condition. */ - BodyAlignedBoxByCell emitter_buffer( - water_block, makeShared(Transform(Vec2d(emitter_buffer_translation)), emitter_buffer_halfsize)); - SimpleDynamics> emitter_buffer_inflow_condition(emitter_buffer); - BodyAlignedBoxByCell disposer( - water_block, makeShared(Transform(Vec2d(disposer_translation)), disposer_halfsize)); - SimpleDynamics disposer_outflow_deletion(disposer, 0); - /** time-space method to detect surface particles. */ - InteractionWithUpdate - free_stream_surface_indicator(water_block_complex); - /** Evaluation of density by freestream approach. */ - InteractionWithUpdate update_fluid_density(water_block_complex); - /** We can output a method-specific particle data for debug */ - water_block.addBodyStateForRecording("Pressure"); - water_block.addBodyStateForRecording("SurfaceIndicator"); - //fish_body.addBodyStateForRecording("MaterailId"); - /** Time step size without considering sound wave speed. */ - ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); - /** Time step size with considering sound wave speed. */ - ReduceDynamics get_fluid_time_step_size(water_block); - /** modify the velocity of boundary particles with free-stream velocity. */ - SimpleDynamics> velocity_boundary_condition_constraint(water_block); - /** Pressure relaxation using verlet time stepping. */ - Dynamics1Level pressure_relaxation(water_block_complex); - /** correct the velocity of boundary particles with free-stream velocity through the post process of pressure relaxation. */ - pressure_relaxation.post_processes_.push_back(&velocity_boundary_condition_constraint); - Dynamics1Level density_relaxation(water_block_complex); - /** Computing viscous acceleration. */ - InteractionDynamics viscous_acceleration(water_block_complex); - /** Impose transport velocity formulation. */ - InteractionDynamics transport_velocity_correction(water_block_complex); - /** Computing vorticity in the flow. */ - InteractionDynamics compute_vorticity(water_block_inner); - //---------------------------------------------------------------------- - // Algorithms of FSI. - //---------------------------------------------------------------------- - SimpleDynamics fish_body_normal_direction(fish_body); - /** Corrected configuration for the elastic insert body. */ - InteractionWithUpdate fish_body_corrected_configuration(fish_inner); - /** Compute the force exerted on solid body due to fluid pressure and viscosity. */ - InteractionDynamics viscous_force_on_solid(fish_contact); - InteractionDynamics - fluid_force_on_solid_update(fish_contact, viscous_force_on_solid); - /** Compute the average velocity of the insert body. */ - solid_dynamics::AverageVelocityAndAcceleration average_velocity_and_acceleration(fish_body); - //---------------------------------------------------------------------- - // Algorithms of solid dynamics. - //---------------------------------------------------------------------- - SimpleDynamics CompositematerialId(fish_body); + /** A Physics relaxation step. */ + relax_dynamics::RelaxationStepInner relaxation_step_inner(fish_inner); + /** + * @brief Particle relaxation starts here. + */ + random_fish_body_particles.exec(0.25); + relaxation_step_inner.SurfaceBounding().exec(); + write_fish_body.writeToFile(); - ReduceDynamics fish_body_computing_time_step_size(fish_body); - /** Stress relaxation for the inserted body. */ - Dynamics1Level fish_body_stress_relaxation_first_half(fish_inner); - Dynamics1Level fish_body_stress_relaxation_second_half(fish_inner); - /** Update norm .*/ - SimpleDynamics fish_body_update_normal(fish_body); + /** relax particles of the insert body. */ + int ite_p = 0; + while (ite_p < 1000) + { + relaxation_step_inner.exec(); + ite_p += 1; + if (ite_p % 200 == 0) + { + std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the inserted body N = " << ite_p << "\n"; + write_fish_body.writeToFile(ite_p); + } + } + std::cout << "The physics relaxation process of inserted body finish !" << std::endl; - fish_body.addBodyStateForRecording("Density"); - fish_body.addBodyStateForRecording("ActiveStrain"); - fish_body.addBodyStateForRecording("MaterailId"); - //---------------------------------------------------------------------- - // Define the methods for I/O operations and observations of the simulation. - //---------------------------------------------------------------------- - BodyStatesRecordingToVtp write_real_body_states(io_environment, system.real_bodies_); - RestartIO restart_io(io_environment, system.real_bodies_); - RegressionTestTimeAverage>> - write_total_viscous_force_on_insert_body(io_environment, viscous_force_on_solid, "TotalViscousForceOnSolid"); - RegressionTestTimeAverage>> - write_total_force_on_insert_body(io_environment, fluid_force_on_solid_update, "TotalForceOnSolid"); - //---------------------------------------------------------------------- - // Prepare the simulation with cell linked list, configuration - // and case specified initial condition if necessary. - //---------------------------------------------------------------------- - /** initialize cell linked lists for all bodies. */ - system.initializeSystemCellLinkedLists(); - /** initialize configurations for all bodies. */ - system.initializeSystemConfigurations(); - /** computing surface normal direction for the insert body. */ - fish_body_normal_direction.exec(); - /** computing linear reproducing configuration for the insert body. */ - fish_body_corrected_configuration.exec(); - //CompositematerialID.parallel_exec(); - //---------------------------------------------------------------------- - // Setup computing and initial conditions. - //---------------------------------------------------------------------- - size_t number_of_iterations = 0; - int screen_output_interval = 100; - Real End_Time = 10.0; /**< End time. */ - Real D_Time = 0.01; /**< time stamps for output. */ - //---------------------------------------------------------------------- - // Statistics for CPU time - //---------------------------------------------------------------------- - TickCount t1 = TickCount::now(); - TimeInterval interval; - //---------------------------------------------------------------------- - // First output before the main loop. - //---------------------------------------------------------------------- - write_real_body_states.writeToFile(); - //---------------------------------------------------------------------- - // Main loop starts here. - //---------------------------------------------------------------------- - while (GlobalStaticVariables::physical_time_ < End_Time) - { - Real integration_time = 0.0; + /** Output results. */ + write_particle_reload_files.writeToFile(); + return 0; + } + //---------------------------------------------------------------------- + // Define the main numerical methods used in the simulation. + // Note that there may be data dependence on the constructors of these methods. + //---------------------------------------------------------------------- + SimpleDynamics initialize_a_fluid_step(water_block, makeShared(Vec2d::Zero())); + BodyAlignedBoxByParticle emitter( + water_block, makeShared(Transform(Vec2d(emitter_translation)), emitter_halfsize)); + SimpleDynamics emitter_inflow_injection(emitter, 10, 0); + /** Emitter buffer inflow condition. */ + BodyAlignedBoxByCell emitter_buffer( + water_block, makeShared(Transform(Vec2d(emitter_buffer_translation)), emitter_buffer_halfsize)); + SimpleDynamics> emitter_buffer_inflow_condition(emitter_buffer); + BodyAlignedBoxByCell disposer( + water_block, makeShared(Transform(Vec2d(disposer_translation)), disposer_halfsize)); + SimpleDynamics disposer_outflow_deletion(disposer, 0); + /** time-space method to detect surface particles. */ + InteractionWithUpdate + free_stream_surface_indicator(water_block_complex); + /** Evaluation of density by freestream approach. */ + InteractionWithUpdate update_fluid_density(water_block_complex); + /** We can output a method-specific particle data for debug */ + water_block.addBodyStateForRecording("Pressure"); + water_block.addBodyStateForRecording("SurfaceIndicator"); + /** Time step size without considering sound wave speed. */ + ReduceDynamics get_fluid_advection_time_step_size(water_block, U_f); + /** Time step size with considering sound wave speed. */ + ReduceDynamics get_fluid_time_step_size(water_block); + /** modify the velocity of boundary particles with free-stream velocity. */ + SimpleDynamics> velocity_boundary_condition_constraint(water_block); + /** Pressure relaxation using verlet time stepping. */ + Dynamics1Level pressure_relaxation(water_block_complex); + /** correct the velocity of boundary particles with free-stream velocity through the post process of pressure relaxation. */ + pressure_relaxation.post_processes_.push_back(&velocity_boundary_condition_constraint); + Dynamics1Level density_relaxation(water_block_complex); + /** Computing viscous acceleration. */ + InteractionDynamics viscous_acceleration(water_block_complex); + /** Impose transport velocity formulation. */ + InteractionDynamics transport_velocity_correction(water_block_complex); + /** Computing vorticity in the flow. */ + InteractionDynamics compute_vorticity(water_block_inner); + //---------------------------------------------------------------------- + // Algorithms of FSI. + //---------------------------------------------------------------------- + SimpleDynamics fish_body_normal_direction(fish_body); + /** Corrected configuration for the elastic insert body. */ + InteractionWithUpdate fish_body_corrected_configuration(fish_inner); + /** Compute the force exerted on solid body due to fluid pressure and viscosity. */ + InteractionDynamics viscous_force_on_solid(fish_contact); + InteractionDynamics + fluid_force_on_fish_update(fish_contact, viscous_force_on_solid); + /** Compute the average velocity of the insert body. */ + solid_dynamics::AverageVelocityAndAcceleration average_velocity_and_acceleration(fish_body); + //---------------------------------------------------------------------- + // Algorithms of solid dynamics. + //---------------------------------------------------------------------- + SimpleDynamics composite_material_id(fish_body); + SimpleDynamics imposing_active_strain(fish_body); + ReduceDynamics fish_body_computing_time_step_size(fish_body); + /** Stress relaxation for the inserted body. */ + Dynamics1Level fish_body_stress_relaxation_first_half(fish_inner); + Dynamics1Level fish_body_stress_relaxation_second_half(fish_inner); + /** Update norm .*/ + SimpleDynamics fish_body_update_normal(fish_body); + fish_body.addBodyStateForRecording("Density"); + fish_body.addBodyStateForRecording("MaterialID"); + fish_body.addBodyStateForRecording("ActiveStrain"); + //---------------------------------------------------------------------- + // Define the methods for I/O operations and observations of the simulation. + //---------------------------------------------------------------------- + BodyStatesRecordingToVtp write_real_body_states(io_environment, system.real_bodies_); + RestartIO restart_io(io_environment, system.real_bodies_); + RegressionTestTimeAverage>> + write_total_viscous_force_on_insert_body(io_environment, viscous_force_on_solid, "TotalViscousForceOnSolid"); + RegressionTestTimeAverage>> + write_total_force_on_insert_body(io_environment, fluid_force_on_fish_update, "TotalForceOnSolid"); + //---------------------------------------------------------------------- + // Prepare the simulation with cell linked list, configuration + // and case specified initial condition if necessary. + //---------------------------------------------------------------------- + /** initialize cell linked lists for all bodies. */ + system.initializeSystemCellLinkedLists(); + /** initialize configurations for all bodies. */ + system.initializeSystemConfigurations(); + /** computing surface normal direction for the fish. */ + fish_body_normal_direction.exec(); + /** computing linear reproducing configuration for the fish. */ + fish_body_corrected_configuration.exec(); + /** initialize material ids for the fish. */ + composite_material_id.exec(); + //---------------------------------------------------------------------- + // Setup computing and initial conditions. + //---------------------------------------------------------------------- + size_t number_of_iterations = 0; + int screen_output_interval = 100; + Real End_Time = 2.5; /**< End time. */ + Real D_Time = 0.01; /**< time stamps for output. */ + //---------------------------------------------------------------------- + // Statistics for CPU time + //---------------------------------------------------------------------- + TickCount t1 = TickCount::now(); + TimeInterval interval; + //---------------------------------------------------------------------- + // First output before the main loop. + //---------------------------------------------------------------------- + write_real_body_states.writeToFile(); + //---------------------------------------------------------------------- + // Main loop starts here. + //---------------------------------------------------------------------- + while (GlobalStaticVariables::physical_time_ < End_Time) + { + Real integration_time = 0.0; - /** Integrate time (loop) until the next output time. */ - while (integration_time < D_Time) - { - initialize_a_fluid_step.exec(); - Real Dt = get_fluid_advection_time_step_size.exec(); - free_stream_surface_indicator.exec(); - update_fluid_density.exec(); - viscous_acceleration.exec(); - transport_velocity_correction.exec(); + /** Integrate time (loop) until the next output time. */ + while (integration_time < D_Time) + { + initialize_a_fluid_step.exec(); + Real Dt = get_fluid_advection_time_step_size.exec(); + free_stream_surface_indicator.exec(); + update_fluid_density.exec(); + viscous_acceleration.exec(); + transport_velocity_correction.exec(); - /** FSI for viscous force. */ - viscous_force_on_solid.exec(); - /** Update normal direction on elastic body.*/ - fish_body_update_normal.exec(); - size_t inner_ite_dt = 0; - size_t inner_ite_dt_s = 0; - Real relaxation_time = 0.0; - while (relaxation_time < Dt) - { - Real dt = SMIN(get_fluid_time_step_size.exec(), Dt - relaxation_time); - /** Fluid pressure relaxation, first half. */ - pressure_relaxation.exec(dt); - /** FSI for pressure force. */ - fluid_force_on_solid_update.exec(); - /** Fluid pressure relaxation, second half. */ - density_relaxation.exec(dt); + /** FSI for viscous force. */ + viscous_force_on_solid.exec(); + /** Update normal direction on elastic body.*/ + fish_body_update_normal.exec(); + size_t inner_ite_dt = 0; + size_t inner_ite_dt_s = 0; + Real relaxation_time = 0.0; + while (relaxation_time < Dt) + { + Real dt = get_fluid_time_step_size.exec(); + /** Fluid pressure relaxation, first half. */ + pressure_relaxation.exec(dt); + /** FSI for fluid force on solid body. */ + fluid_force_on_fish_update.exec(); + /** Fluid pressure relaxation, second half. */ + density_relaxation.exec(dt); - /** Solid dynamics. */ - inner_ite_dt_s = 0; - Real dt_s_sum = 0.0; - average_velocity_and_acceleration.initialize_displacement_.exec(); - while (dt_s_sum < dt) - { - Real dt_s = SMIN(fish_body_computing_time_step_size.exec(), dt - dt_s_sum); - CompositematerialId.exec(dt_s); - fish_body_stress_relaxation_first_half.exec(dt_s); - fish_body_stress_relaxation_second_half.exec(dt_s); - dt_s_sum += dt_s; - inner_ite_dt_s++; - } - average_velocity_and_acceleration.update_averages_.exec(dt); + /** Solid dynamics. */ + inner_ite_dt_s = 0; + Real dt_s_sum = 0.0; + average_velocity_and_acceleration.initialize_displacement_.exec(); + while (dt_s_sum < dt) + { + Real dt_s = SMIN(fish_body_computing_time_step_size.exec(), dt - dt_s_sum); + imposing_active_strain.exec(); + fish_body_stress_relaxation_first_half.exec(dt_s); + fish_body_stress_relaxation_second_half.exec(dt_s); + dt_s_sum += dt_s; + inner_ite_dt_s++; + } + average_velocity_and_acceleration.update_averages_.exec(dt); - relaxation_time += dt; - integration_time += dt; - GlobalStaticVariables::physical_time_ += dt; - emitter_buffer_inflow_condition.exec(dt); - inner_ite_dt++; - } + relaxation_time += dt; + integration_time += dt; + GlobalStaticVariables::physical_time_ += dt; + emitter_buffer_inflow_condition.exec(dt); + inner_ite_dt++; + } - if (number_of_iterations % screen_output_interval == 0) - { - std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " - << GlobalStaticVariables::physical_time_ - << " Dt = " << Dt << " Dt / dt = " << inner_ite_dt << " dt / dt_s = " << inner_ite_dt_s << "\n"; - } - number_of_iterations++; + if (number_of_iterations % screen_output_interval == 0) + { + std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = " + << GlobalStaticVariables::physical_time_ + << " Dt = " << Dt << " Dt / dt = " << inner_ite_dt << " dt / dt_s = " << inner_ite_dt_s << "\n"; + } + number_of_iterations++; - /** Water block configuration and periodic condition. */ - emitter_inflow_injection.exec(); - disposer_outflow_deletion.exec(); + /** Water block configuration and periodic condition. */ + emitter_inflow_injection.exec(); + disposer_outflow_deletion.exec(); - water_block.updateCellLinkedListWithParticleSort(100); - fish_body.updateCellLinkedList(); - /** one need update configuration after periodic condition. */ - water_block_complex.updateConfiguration(); - /** one need update configuration after periodic condition. */ - fish_contact.updateConfiguration(); - - } + water_block.updateCellLinkedListWithParticleSort(100); + fish_body.updateCellLinkedList(); + /** one need update configuration after periodic condition. */ + water_block_complex.updateConfiguration(); + /** one need update configuration after periodic condition. */ + fish_contact.updateConfiguration(); + } - TickCount t2 = TickCount::now(); - /** write run-time observation into file */ - compute_vorticity.exec(); - write_real_body_states.writeToFile(); - write_total_viscous_force_on_insert_body.writeToFile(number_of_iterations); + TickCount t2 = TickCount::now(); + /** write run-time observation into file */ + compute_vorticity.exec(); + write_real_body_states.writeToFile(); + write_total_viscous_force_on_insert_body.writeToFile(number_of_iterations); - TickCount t3 = TickCount::now(); - interval += t3 - t2; - } - TickCount t4 = TickCount::now(); + 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; + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; - return 0; + return 0; } diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.h b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.h index 3595fe2bfa..c60e993ba6 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.h +++ b/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.h @@ -1,17 +1,17 @@ -#include "sphinxsys.h" #include "2d_fish_and_bones.h" #include "composite_material.h" +#include "sphinxsys.h" #define PI 3.1415926 using namespace SPH; //---------------------------------------------------------------------- // Basic geometry parameters and numerical setup. //---------------------------------------------------------------------- -Real DL = 0.8; /**< Channel length. */ -Real DH = 0.4; /**< Channel height. */ -Real particle_spacing_ref = 0.0025; /**< Initial reference particle spacing. */ +Real DL = 0.8; /**< Channel length. */ +Real DH = 0.4; /**< Channel height. */ +Real particle_spacing_ref = 0.0025; /**< Initial reference particle spacing. */ Real DL_sponge = particle_spacing_ref * 20.0; /**< Sponge region to impose inflow condition. */ -Real BW = particle_spacing_ref * 4.0; /**< Boundary width, determined by specific layer of boundary particles. */ +Real BW = particle_spacing_ref * 4.0; /**< Boundary width, determined by specific layer of boundary particles. */ /** Domain bounds of the system. */ BoundingBox system_domain_bounds(Vec2d(-DL_sponge - BW, -BW), Vec2d(DL + BW, DH + BW)); @@ -27,20 +27,20 @@ Vec2d disposer_translation = Vec2d(DL, DH + 0.25 * DH) - disposer_halfsize; //---------------------------------------------------------------------- // Material properties of the fluid. //---------------------------------------------------------------------- -Real rho0_f = 1000.0; /**< Density. */ -Real U_f = 1.0; /**< freestream velocity. */ -Real c_f = 10.0 * U_f; /**< Speed of sound. */ -Real Re = 30000.0; /**< Reynolds number. */ -Real mu_f = rho0_f * U_f * 0.3 / Re; /**< Dynamics viscosity. */ +Real rho0_f = 1000.0; /**< Density. */ +Real U_f = 1.0; /**< freestream velocity. */ +Real c_f = 10.0 * U_f; /**< Speed of sound. */ +Real Re = 30000.0; /**< Reynolds number. */ +Real mu_f = rho0_f * U_f * 0.3 / Re; /**< Dynamics viscosity. */ //---------------------------------------------------------------------- // Global parameters on the solid properties //---------------------------------------------------------------------- //---------------------------------------------------------------------- -Real cx = 0.3 * DL; /**< Center of fish in x direction. */ +Real cx = 0.3 * DL; /**< Center of fish in x direction. */ Real cy = DH / 2; /**< Center of fish in y direction. */ Real fish_length = 0.2; /**< Length of fish. */ Real fish_thickness = 0.03; /**< The maximum fish thickness. */ -Real muscel_thickness = 0.02; /**< The maximum fish thickness. */ +Real muscle_thickness = 0.02; /**< The maximum fish thickness. */ Real head_length = 0.03; /**< Length of fish bone. */ Real bone_thickness = 0.003; /**< Length of fish bone. */ Real fish_shape_resolution = particle_spacing_ref * 0.5; @@ -57,40 +57,40 @@ Real a3 = -15.73 * fish_thickness / pow(fish_length, 3); Real a4 = 21.87 * fish_thickness / pow(fish_length, 4); Real a5 = -10.55 * fish_thickness / pow(fish_length, 5); -Real b1 = 1.22 * muscel_thickness / fish_length; -Real b2 = 3.19 * muscel_thickness / fish_length / fish_length; -Real b3 = -15.73 * muscel_thickness / pow(fish_length, 3); -Real b4 = 21.87 * muscel_thickness / pow(fish_length, 4); -Real b5 = -10.55 * muscel_thickness / pow(fish_length, 5); +Real b1 = 1.22 * muscle_thickness / fish_length; +Real b2 = 3.19 * muscle_thickness / fish_length / fish_length; +Real b3 = -15.73 * muscle_thickness / pow(fish_length, 3); +Real b4 = 21.87 * muscle_thickness / pow(fish_length, 4); +Real b5 = -10.55 * muscle_thickness / pow(fish_length, 5); //---------------------------------------------------------------------- // SPH bodies with cases-dependent geometries (ComplexShape). //---------------------------------------------------------------------- /** create a water block shape */ std::vector createWaterBlockShape() { - //geometry - std::vector water_block_shape; - water_block_shape.push_back(Vecd(-DL_sponge, 0.0)); - water_block_shape.push_back(Vecd(-DL_sponge, DH)); - water_block_shape.push_back(Vecd(DL, DH)); - water_block_shape.push_back(Vecd(DL, 0.0)); - water_block_shape.push_back(Vecd(-DL_sponge, 0.0)); - - return water_block_shape; + // geometry + std::vector water_block_shape; + water_block_shape.push_back(Vecd(-DL_sponge, 0.0)); + water_block_shape.push_back(Vecd(-DL_sponge, DH)); + water_block_shape.push_back(Vecd(DL, DH)); + water_block_shape.push_back(Vecd(DL, 0.0)); + water_block_shape.push_back(Vecd(-DL_sponge, 0.0)); + + return water_block_shape; } /** -* Fish body with tethering constraint. -*/ + * Fish body with tethering constraint. + */ class FishBody : public MultiPolygonShape { -public: - explicit FishBody(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - std::vector fish_shape = CreatFishShape(cx, cy, fish_length, fish_shape_resolution); - multi_polygon_.addAPolygon(fish_shape, ShapeBooleanOps::add); - } + public: + explicit FishBody(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + std::vector fish_shape = CreatFishShape(cx, cy, fish_length, fish_shape_resolution); + multi_polygon_.addAPolygon(fish_shape, ShapeBooleanOps::add); + } }; //---------------------------------------------------------------------- // Define case dependent bodies material, constraint and boundary conditions. @@ -98,122 +98,137 @@ class FishBody : public MultiPolygonShape /** Fluid body definition */ class WaterBlock : public ComplexShape { -public: - explicit WaterBlock(const std::string& shape_name) : ComplexShape(shape_name) - { - /** Geomtry definition. */ - MultiPolygon outer_boundary(createWaterBlockShape()); - add(outer_boundary, "OuterBoundary"); - MultiPolygon fish(CreatFishShape(cx, cy, fish_length, fish_shape_resolution)); - subtract(fish); - } + public: + explicit WaterBlock(const std::string &shape_name) : ComplexShape(shape_name) + { + /** Geometry definition. */ + MultiPolygon outer_boundary(createWaterBlockShape()); + add(outer_boundary, "OuterBoundary"); + MultiPolygon fish(CreatFishShape(cx, cy, fish_length, fish_shape_resolution)); + subtract(fish); + } }; /** Case dependent inflow boundary condition. */ struct FreeStreamVelocity { - Real u_ref_, t_ref_; - - template - FreeStreamVelocity(BoundaryConditionType& boundary_condition) - : u_ref_(0.0), t_ref_(2.0) {} - - Vecd operator()(Vecd& position, Vecd& velocity) - { - Vecd target_velocity = Vecd::Zero(); - Real run_time = GlobalStaticVariables::physical_time_; - target_velocity[0] = run_time < t_ref_ ? 0.5 * u_ref_ * (1.0 - cos(Pi * run_time / t_ref_)) : u_ref_; - return target_velocity; - } + Real u_ref_, t_ref_; + + template + FreeStreamVelocity(BoundaryConditionType &boundary_condition) + : u_ref_(0.0), t_ref_(2.0) {} + + Vecd operator()(Vecd &position, Vecd &velocity) + { + Vecd target_velocity = Vecd::Zero(); + Real run_time = GlobalStaticVariables::physical_time_; + target_velocity[0] = run_time < t_ref_ ? 0.5 * u_ref_ * (1.0 - cos(Pi * run_time / t_ref_)) : u_ref_; + return target_velocity; + } }; - //---------------------------------------------------------------------- // Define time dependent acceleration in x-direction //---------------------------------------------------------------------- class TimeDependentAcceleration : public Gravity { - Real t_ref_, u_ref_, du_ave_dt_; + Real t_ref_, u_ref_, du_ave_dt_; -public: - explicit TimeDependentAcceleration(Vecd gravity_vector) - : Gravity(gravity_vector), t_ref_(2.0), u_ref_(0.00), du_ave_dt_(0) {} + public: + explicit TimeDependentAcceleration(Vecd gravity_vector) + : Gravity(gravity_vector), t_ref_(2.0), u_ref_(0.00), du_ave_dt_(0) {} - virtual Vecd InducedAcceleration(Vecd& position) override - { - Real run_time_ = GlobalStaticVariables::physical_time_; - du_ave_dt_ = 0.5 * u_ref_ * (Pi / t_ref_) * sin(Pi * run_time_ / t_ref_); + virtual Vecd InducedAcceleration(Vecd &position) override + { + Real run_time_ = GlobalStaticVariables::physical_time_; + du_ave_dt_ = 0.5 * u_ref_ * (Pi / t_ref_) * sin(Pi * run_time_ / t_ref_); - return run_time_ < t_ref_ ? Vecd(du_ave_dt_, 0.0) : global_acceleration_; - } + return run_time_ < t_ref_ ? Vecd(du_ave_dt_, 0.0) : global_acceleration_; + } }; - -//Material ID +// Material ID class SolidBodyMaterial : public CompositeMaterial { -public: - SolidBodyMaterial() : CompositeMaterial(rho0_s) - { - add(rho0_s, Youngs_modulus1, poisson); - add(rho0_s, Youngs_modulus2, poisson); - add(rho0_s, Youngs_modulus3, poisson); - }; + public: + SolidBodyMaterial() : CompositeMaterial(rho0_s) + { + add(rho0_s, Youngs_modulus1, poisson); + add(rho0_s, Youngs_modulus2, poisson); + add(rho0_s, Youngs_modulus3, poisson); + }; +}; +//---------------------------------------------------------------------- +// Case dependent initialization material ids +//---------------------------------------------------------------------- +class FishMaterialInitialization + : public MaterialIdInitialization +{ + public: + explicit FishMaterialInitialization(SolidBody &solid_body) + : MaterialIdInitialization(solid_body){}; + + void update(size_t index_i, Real dt = 0.0) + { + Real x = pos0_[index_i][0] - cx; + Real y = pos0_[index_i][1]; + + Real y1 = a1 * pow(x, 0 + 1) + a2 * pow(x, 1 + 1) + a3 * pow(x, 2 + 1) + a4 * pow(x, 3 + 1) + a5 * pow(x, 4 + 1); + if (x <= (fish_length - head_length) && y > (y1 - 0.004 + cy) && y > (cy + bone_thickness / 2)) + { + material_id_[index_i] = 0; // region for muscle + } + else if (x <= (fish_length - head_length) && y < (-y1 + 0.004 + cy) && y < (cy - bone_thickness / 2)) + { + material_id_[index_i] = 0; // region for muscle + } + else if ((x > (fish_length - head_length)) || ((y < (cy + bone_thickness / 2)) && (y > (cy - bone_thickness / 2)))) + { + material_id_[index_i] = 2; + } + else + { + material_id_[index_i] = 1; + } + }; }; -// Setup material ID -class MaterialId - : public solid_dynamics::ElasticDynamicsInitialCondition +// imposing active strain to fish muscle +class ImposingActiveStrain + : public solid_dynamics::ElasticDynamicsInitialCondition { -public: - explicit MaterialId(SolidBody& solid_body) - : solid_dynamics::ElasticDynamicsInitialCondition(solid_body), - solid_particles_(dynamic_cast(&solid_body.getBaseParticles())), - materail_id_(*solid_particles_->getVariableByName("MaterailId")), - pos0_(solid_particles_->pos0_) - { - solid_particles_->registerVariable(active_strain_, "ActiveStrain"); - }; - virtual void update(size_t index_i, Real dt = 0.0) - { - Real x = pos0_[index_i][0]-cx; - Real y = pos0_[index_i][1]; - Real y1(0); - - y1 = a1 * pow(x, 0 + 1) + a2 * pow(x, 1 + 1) + a3 * pow(x, 2 + 1) + a4 * pow(x, 3 + 1) + a5 * pow(x, 4 + 1); - - Real Am = 0.12; - Real frequency = 4.0; - Real w = 2 * PI * frequency; - Real lamda = 3.0 * fish_length; - Real wave_number = 2 * PI / lamda; - Real hx = -(pow(x, 2)- pow(fish_length, 2)) / pow(fish_length, 2); - Real ta = 0.2; - Real st = 1 - exp(-GlobalStaticVariables::physical_time_ / ta); - - active_strain_[index_i] = Matd::Zero(); - - if (x <=(fish_length - head_length) && y > (y1-0.004 + cy) && y > (cy + bone_thickness / 2)) - { - materail_id_[index_i] = 0; - active_strain_[index_i](0, 0) = -Am * hx * st * pow(sin(w * GlobalStaticVariables::physical_time_/2 + wave_number * x/2), 2); - } - else if (x <= (fish_length - head_length) && y < (-y1 + 0.004 + cy) && y <(cy - bone_thickness / 2)) - { - materail_id_[index_i] = 0; - active_strain_[index_i](0, 0) = -Am * hx * st * pow(sin(w * GlobalStaticVariables::physical_time_/2 + wave_number * x/2 + PI / 2), 2); - } - else if ((x > (fish_length - head_length)) || ((y < (cy + bone_thickness / 2)) && (y > (cy - bone_thickness / 2)))) - { - materail_id_[index_i] = 2; - } - else - { - materail_id_[index_i] = 1; - } - }; -protected: - SolidParticles* solid_particles_; - StdLargeVec& materail_id_; - StdLargeVec& pos0_; - StdLargeVec active_strain_; + public: + explicit ImposingActiveStrain(SolidBody &solid_body) + : solid_dynamics::ElasticDynamicsInitialCondition(solid_body), + material_id_(*particles_->getVariableByName("MaterialID")), + pos0_(*particles_->getVariableByName("InitialPosition")) + { + particles_->registerVariable(active_strain_, "ActiveStrain"); + }; + virtual void update(size_t index_i, Real dt = 0.0) + { + if (material_id_[index_i] == 0) + { + Real x = pos0_[index_i][0] - cx; + Real y = pos0_[index_i][1]; + + Real Am = 0.12; + Real frequency = 4.0; + Real w = 2 * Pi * frequency; + Real lambda = 3.0 * fish_length; + Real wave_number = 2 * Pi / lambda; + Real hx = -(pow(x, 2) - pow(fish_length, 2)) / pow(fish_length, 2); + Real start_time = 0.2; + Real current_time = GlobalStaticVariables::physical_time_; + Real strength = 1 - exp(-current_time / start_time); + + Real phase_shift = y > (cy + bone_thickness / 2) ? 0 : Pi / 2; + active_strain_[index_i](0, 0) = + -Am * hx * strength * pow(sin(w * current_time / 2 + wave_number * x / 2 + phase_shift), 2); + } + }; + + protected: + StdLargeVec &material_id_; + StdLargeVec &pos0_; + StdLargeVec active_strain_; }; diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/CMakeLists.txt b/tests/user_examples/test_2d_flow_stream_around_fish/CMakeLists.txt index 9c995da4c4..8fc0e129ec 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/CMakeLists.txt +++ b/tests/user_examples/test_2d_flow_stream_around_fish/CMakeLists.txt @@ -13,5 +13,17 @@ SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") aux_source_directory(. DIR_SRCS) ADD_EXECUTABLE(${PROJECT_NAME} ${DIR_SRCS}) - set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") - target_link_libraries(${PROJECT_NAME} sphinxsys_2d) +if(${CMAKE_SYSTEM_NAME} MATCHES "Windows") + add_test(NAME ${PROJECT_NAME}_particle_relaxation COMMAND ${PROJECT_NAME} --r=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME} --r=false --i=true + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +else() + file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/run_test.sh + DESTINATION ${EXECUTABLE_OUTPUT_PATH}) + add_test(NAME ${PROJECT_NAME} COMMAND bash ${EXECUTABLE_OUTPUT_PATH}/run_test.sh + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) +endif() + +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") +target_link_libraries(${PROJECT_NAME} sphinxsys_2d) diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/active_model.cpp b/tests/user_examples/test_2d_flow_stream_around_fish/active_model.cpp index 9b54977e31..f1bbc19a9f 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/active_model.cpp +++ b/tests/user_examples/test_2d_flow_stream_around_fish/active_model.cpp @@ -10,7 +10,7 @@ namespace SPH ActiveIntegration1stHalf:: ActiveIntegration1stHalf(BaseInnerRelation &inner_relation) : Integration1stHalfPK2(inner_relation), active_strain_(*particles_->getVariableByName("ActiveStrain")) - , materail_id_(*particles_->getVariableByName("MaterailId")) + , material_id_(*particles_->getVariableByName("MaterialID")) { particles_->registerVariable(F_0, "ActiveTensor"); particles_->registerVariable(E_e, "ElasticStrain"); @@ -22,14 +22,14 @@ namespace SPH // Calculate the active deformation gradient F_0[index_i] = Matd::Identity(); - F_0[index_i] = (materail_id_[index_i] == 0) ? ((2.0 * active_strain_[index_i] + Matd::Identity()).llt().matrixL()) : F_0[index_i]; + F_0[index_i] = (material_id_[index_i] == 0) ? ((2.0 * active_strain_[index_i] + Matd::Identity()).llt().matrixL()) : F_0[index_i]; // Calculate the elastic deformation Matd F_e = F_[index_i] * F_0[index_i].inverse(); Matd F_0_star = F_0[index_i].determinant() * (F_0[index_i].inverse()).transpose(); // Calculate the elastic strain - E_e[index_i] = (materail_id_[index_i] == 0) ? (0.5 * (F_[index_i].transpose() * F_[index_i] - Matd::Identity()) - active_strain_[index_i]) : F_[index_i]; + E_e[index_i] = (material_id_[index_i] == 0) ? (0.5 * (F_[index_i].transpose() * F_[index_i] - Matd::Identity()) - active_strain_[index_i]) : F_[index_i]; stress_PK1_B_[index_i] = F_e * elastic_solid_.StressPK2(E_e[index_i], index_i) * F_0_star * B_[index_i]; } diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/active_model.h b/tests/user_examples/test_2d_flow_stream_around_fish/active_model.h index a6d6db3527..3bcef409db 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/active_model.h +++ b/tests/user_examples/test_2d_flow_stream_around_fish/active_model.h @@ -20,7 +20,7 @@ namespace SPH protected: StdLargeVec &active_strain_,F_0, E_e; - StdLargeVec& materail_id_; + StdLargeVec& material_id_; }; } } diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.cpp b/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.cpp index ef753b3785..217bc9fdd9 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.cpp +++ b/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.cpp @@ -1,31 +1,30 @@ #include "composite_material.h" +#include "base_local_dynamics.h" #include "base_particles.hpp" -#include - namespace SPH { - //=================================================================================================// - void CompositeMaterial::initializeLocalParameters(BaseParticles* base_particles) - { - ElasticSolid::initializeLocalParameters(base_particles); - for (size_t i = 0; i < CompositeMaterails_.size(); ++i) - CompositeMaterails_[i]->initializeLocalParameters(base_particles); - - for (size_t j = 0; j < CompositeMaterails_.size(); ++j) - sound_speed_.push_back(CompositeMaterails_[j]->ReferenceSoundSpeed()); - - c0_ = *std::max_element(sound_speed_.begin(), sound_speed_.end()); - setContactStiffness(c0_); - - base_particles->registerVariable(materail_id_, "MaterailId"); - } +//=================================================================================================// +void CompositeMaterial::initializeLocalParameters(BaseParticles *base_particles) +{ + ElasticSolid::initializeLocalParameters(base_particles); + base_particles->registerVariable(material_id_, "MaterialID"); - //=================================================================================================// - Matd ActiveModelSolid::StressPK2(Matd& F, size_t particle_index_i) - { - return lambda0_ * F.trace() * Matd::Identity() + 2.0 * G0_ * F; - } - //=================================================================================================// + for (size_t i = 0; i < composite_materials_.size(); ++i) + { + composite_materials_[i]->initializeLocalParameters(base_particles); + } +} +//=================================================================================================// +MaterialIdInitialization::MaterialIdInitialization(SPHBody &sph_body) + : LocalDynamics(sph_body), GeneralDataDelegateSimple(sph_body), + material_id_(*particles_->getVariableByName("MaterialID")), + pos0_(*particles_->getVariableByName("InitialPosition")){}; +//=================================================================================================// +Matd ActiveModelSolid::StressPK2(Matd &F, size_t particle_index_i) +{ + return lambda0_ * F.trace() * Matd::Identity() + 2.0 * G0_ * F; } -//=================================================================================================// \ No newline at end of file +//=================================================================================================// +} // namespace SPH + //=================================================================================================// \ No newline at end of file diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.h b/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.h index 2ec4aeb42b..240654bae3 100644 --- a/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.h +++ b/tests/user_examples/test_2d_flow_stream_around_fish/composite_material.h @@ -3,73 +3,92 @@ #define COMPOSITE_MATERIAL_H #include "elastic_solid.h" -#include + +#include namespace SPH { - /**@class CompositeMaterial*/ - class CompositeMaterial : public ElasticSolid - { - protected: - UniquePtrsKeeper composite_ptr_keeper_; - StdVec CompositeMaterails_; - /** initialize the local properties, fiber and sheet direction. */ - std::vector sound_speed_; - public: - StdLargeVec materail_id_; - - explicit CompositeMaterial(Real rho0) : ElasticSolid(rho0) - { - material_type_name_ = "CompositeMaterial"; - }; - virtual ~CompositeMaterial() {}; - - virtual void initializeLocalParameters(BaseParticles* base_particles) override; - - virtual Matd StressPK2(Matd& deformation, size_t particle_index_i) override - { - return CompositeMaterails_[materail_id_[particle_index_i]]->StressPK2(deformation, particle_index_i); - }; - - Real CompositeDensity(size_t particle_index_i) - { - return CompositeMaterails_[materail_id_[particle_index_i]]->ReferenceDensity(); - }; - - virtual Matd StressCauchy(Matd& almansi_strain, Matd& F, size_t particle_index_i) override - { - return Matd::Identity(); - }; - - virtual Real VolumetricKirchhoff(Real J) override - { - return 0.0; - }; - - virtual std::string getRelevantStressMeasureName() override { return "PK2"; }; - - template - void add(ConstructorArgs &&...args) - { - CompositeMaterails_.push_back(composite_ptr_keeper_.createPtr(std::forward(args)...)); - }; - }; - - /** - * @class Activemodel - */ - class ActiveModelSolid : public SaintVenantKirchhoffSolid - { - public: - explicit ActiveModelSolid(Real rho0, Real youngs_modulus, Real poisson_ratio) - : SaintVenantKirchhoffSolid(rho0, youngs_modulus, poisson_ratio) - { - material_type_name_ = "ActiveModelSolid"; - }; - virtual ~ActiveModelSolid() {}; - - /** second Piola-Kirchhoff stress related with green-lagrangian deformation tensor */ - virtual Matd StressPK2(Matd& deformation, size_t particle_index_i); - }; -} -#endif //COMPOSITE_MATERAIL_H \ No newline at end of file +/**@class CompositeMaterial*/ +class CompositeMaterial : public ElasticSolid +{ + protected: + UniquePtrsKeeper composite_ptrs_keeper_; + StdVec composite_materials_; + + public: + StdLargeVec material_id_; + + explicit CompositeMaterial(Real rho0) : ElasticSolid(rho0) + { + material_type_name_ = "CompositeMaterial"; + }; + virtual ~CompositeMaterial(){}; + + virtual void initializeLocalParameters(BaseParticles *base_particles) override; + + virtual Matd StressPK2(Matd &deformation, size_t particle_index_i) override + { + return composite_materials_[material_id_[particle_index_i]]->StressPK2(deformation, particle_index_i); + }; + + Real CompositeDensity(size_t particle_index_i) + { + return composite_materials_[material_id_[particle_index_i]]->ReferenceDensity(); + }; + + virtual Matd StressCauchy(Matd &almansi_strain, Matd &F, size_t particle_index_i) override + { + return Matd::Identity(); + }; + + virtual Real VolumetricKirchhoff(Real J) override + { + return 0.0; + }; + + virtual std::string getRelevantStressMeasureName() override { return "PK2"; }; + + template + void add(Args &&...args) + { + ElasticSolid *added_material = + composite_ptrs_keeper_.createPtr(std::forward(args)...); + composite_materials_.push_back(added_material); + c0_ = SMAX(c0_, added_material->ReferenceSoundSpeed()); + setContactStiffness(c0_); + }; +}; + +/** + * @class MaterialIdInitialization + */ +class MaterialIdInitialization + : public LocalDynamics, + public GeneralDataDelegateSimple +{ + public: + explicit MaterialIdInitialization(SPHBody &sph_body); + + protected: + StdLargeVec &material_id_; + StdLargeVec &pos0_; +}; + +/** + * @class ActiveModelSolid + */ +class ActiveModelSolid : public SaintVenantKirchhoffSolid +{ + public: + explicit ActiveModelSolid(Real rho0, Real youngs_modulus, Real poisson_ratio) + : SaintVenantKirchhoffSolid(rho0, youngs_modulus, poisson_ratio) + { + material_type_name_ = "ActiveModelSolid"; + }; + virtual ~ActiveModelSolid(){}; + + /** second Piola-Kirchhoff stress related with green-lagrangian deformation tensor */ + virtual Matd StressPK2(Matd &deformation, size_t particle_index_i); +}; +} // namespace SPH +#endif // COMPOSITE_MATERIAL_H diff --git a/tests/user_examples/test_2d_flow_stream_around_fish/run_test.sh b/tests/user_examples/test_2d_flow_stream_around_fish/run_test.sh new file mode 100644 index 0000000000..5ad71f9c2e --- /dev/null +++ b/tests/user_examples/test_2d_flow_stream_around_fish/run_test.sh @@ -0,0 +1,2 @@ +./test_2d_flow_stream_around_fish --r=true +./test_2d_flow_stream_around_fish --r=false --i=true