Add rotation Patch and impact patch cases

src/shared/particle_dynamics/general_dynamics/general_interaction.cpp

@@ -9,6 +9,7 @@ CorrectedConfigurationInner::
     : LocalDynamics(inner_relation.getSPHBody()),
       GeneralDataDelegateInner(inner_relation),
       beta_(beta), alpha_(alpha),
+      smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()),
       B_(*particles_->getVariableByName<Matd>("CorrectionMatrix")) {}
 //=================================================================================================//
 void CorrectedConfigurationInner::interaction(size_t index_i, Real dt)
@@ -30,7 +31,9 @@ void CorrectedConfigurationInner::update(size_t index_i, Real dt)
 {
     Real det_sqr = pow(B_[index_i].determinant(), beta_);
     Matd inverse = B_[index_i].inverse();
-    B_[index_i] = (det_sqr * inverse + alpha_ * Matd::Identity()) / (alpha_ + det_sqr);
+    Real weight1_ = det_sqr / (alpha_ * pow(smoothing_length_, beta_) + det_sqr);
+    Real weight2_ = alpha_ * pow(smoothing_length_, beta_) / (alpha_ * pow(smoothing_length_, beta_) + det_sqr);
+    B_[index_i] = weight1_ * inverse + weight2_ * Matd::Identity();
 }
 //=================================================================================================//
 CorrectedConfigurationComplex::

src/shared/particle_dynamics/general_dynamics/general_interaction.h

@@ -42,6 +42,7 @@ class CorrectedConfigurationInner : public LocalDynamics, public GeneralDataDele
   protected:
     int beta_;
     Real alpha_;
+    Real smoothing_length_;
     StdLargeVec<Matd> &B_;
 
     void interaction(size_t index_i, Real dt = 0.0);

tests/user_examples/test_2d_impact_patch/CMakeLists.txt

@@ -0,0 +1,25 @@
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir
+
+set(CMAKE_VERBOSE_MAKEFILE on)
+
+STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} )
+PROJECT("${CURRENT_FOLDER}")
+
+SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/")
+SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input")
+SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload")
+
+file(MAKE_DIRECTORY ${BUILD_INPUT_PATH})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/
+     DESTINATION ${BUILD_INPUT_PATH})
+
+aux_source_directory(. DIR_SRCS)
+add_executable(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS} )
+set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}")
+target_link_libraries(${PROJECT_NAME} extra_sources_2d)
+
+add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME}
+                 WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+add_test(NAME ${PROJECT_NAME}_restart COMMAND ${PROJECT_NAME} --restart_step=4000
+                 WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})

tests/user_examples/test_2d_impact_patch/impact_patch.cpp

@@ -0,0 +1,270 @@
+/**
+ * @file	impactpatch.cpp
+ * @brief	2D impactpatch example.
+ * @author	Yaru Ren, Chi Zhang and Xiangyu Hu
+ */
+#include "sphinxsys.h" //SPHinXsys Library.
+using namespace SPH;   // Namespace cite here.
+#define PI 3.1415926
+//----------------------------------------------------------------------
+//	Basic geometry parameters and numerical setup.
+//----------------------------------------------------------------------
+Real LL = 0.667;                        /**< Liquid column length. */
+Real LH = 2.0;                         /**< Liquid column height. */
+Real particle_spacing_ref = LL / 150;    /**< Initial reference particle spacing. */
+Real BW = particle_spacing_ref * 4;  /**< Extending width for boundary conditions. */
+BoundingBox system_domain_bounds(Vec2d(-LH, -LH), Vec2d(LL + BW, LH + BW));
+//----------------------------------------------------------------------
+//	Material parameters.
+//----------------------------------------------------------------------
+Real rho0_f = 1.0;                    /**< Reference density of fluid. */
+Real U_max = 1.0;                    /**< Characteristic velocity. */
+Real c_f = 100.0;                   /**< Reference sound speed. */
+//----------------------------------------------------------------------
+//	Geometric shapes used in this case.
+//----------------------------------------------------------------------
+Vec2d DamP_lb(-LL / 2, -LH / 2);	 /**< Left bottom. */
+Vec2d DamP_lt(-LL / 2, LH / 2);	 /**< Left top. */
+Vec2d DamP_rt(LL / 2, LH / 2); /**< Right top. */
+Vec2d DamP_rb(LL / 2, -LH / 2);	 /**< Right bottom. */
+//----------------------------------------------------------------------
+//	Geometry definition.
+//----------------------------------------------------------------------
+std::vector<Vecd> createWaterBlockShape()
+{
+    std::vector<Vecd> water_block_shape;
+    water_block_shape.push_back(DamP_lb);
+    water_block_shape.push_back(DamP_lt);
+    water_block_shape.push_back(DamP_rt);
+    water_block_shape.push_back(DamP_rb);
+    water_block_shape.push_back(DamP_lb);
+
+    return water_block_shape;
+}
+class WaterBlock : public MultiPolygonShape
+{
+public:
+    explicit WaterBlock(const std::string& shape_name) : MultiPolygonShape(shape_name)
+    {
+        multi_polygon_.addAPolygon(createWaterBlockShape(), ShapeBooleanOps::add);
+    }
+};
+
+/**
+ * application dependent initial velocity
+ */
+class InitialVelocity
+    : public fluid_dynamics::FluidInitialCondition
+{
+public:
+    InitialVelocity(SPHBody& sph_body)
+        : fluid_dynamics::FluidInitialCondition(sph_body) {};
+
+    void update(size_t index_i, Real dt)
+    {
+        /** initial velocity profile */
+        if (pos_[index_i][1]>=0.0)
+        {
+            vel_[index_i][1] = -1.0;
+        }
+        else
+        {
+            vel_[index_i][1] =  1.0;
+        }        
+
+    }
+};
+//----------------------------------------------------------------------
+//	wave gauge
+//----------------------------------------------------------------------
+Real h = 1.3 * particle_spacing_ref;
+MultiPolygon createWaveProbeShape()
+{
+    std::vector<Vecd> pnts;
+    pnts.push_back(Vecd(1.0 - h, 0.0));
+    pnts.push_back(Vecd(1.0 - h, LH));
+    pnts.push_back(Vecd(1.0 + h, LH));
+    pnts.push_back(Vecd(1.0 + h, 0.0));
+    pnts.push_back(Vecd(1.0 - h, 0.0));
+
+    MultiPolygon multi_polygon;
+    multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add);
+    return multi_polygon;
+}
+//----------------------------------------------------------------------
+//	Main program starts here.
+//----------------------------------------------------------------------
+int main(int ac, char *av[])
+{
+    //----------------------------------------------------------------------
+    //	Build up an SPHSystem.
+    //----------------------------------------------------------------------
+    BoundingBox system_domain_bounds(Vec2d(-LH, -LH), Vec2d(LL + BW, LH + BW));
+    SPHSystem sph_system(system_domain_bounds, particle_spacing_ref);
+    /** Tag for run particle relaxation for the initial body fitted distribution. */
+    sph_system.setRunParticleRelaxation(false);
+    /** Tag for computation start with relaxed body fitted particles distribution. */
+    sph_system.setReloadParticles(false);
+    sph_system.handleCommandlineOptions(ac, av);
+    IOEnvironment io_environment(sph_system);
+    //----------------------------------------------------------------------
+    //	Creating bodies with corresponding materials and particles.
+    //----------------------------------------------------------------------
+    FluidBody water_block(sph_system, makeShared<WaterBlock>("WaterBody"));
+    water_block.defineAdaptation<SPHAdaptation>(1.3, 1.0);
+    water_block.defineParticlesAndMaterial<BaseParticles, WeaklyCompressibleFluid>(rho0_f, c_f);
+    // Using relaxed particle distribution if needed
+    (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
+        ? water_block.generateParticles<ParticleGeneratorReload>(io_environment, water_block.getName())
+        : water_block.generateParticles<ParticleGeneratorLattice>();
+    //----------------------------------------------------------------------
+    //	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 water_body_inner(water_block);
+    //----------------------------------------------------------------------
+    //----------------------------------------------------------------------
+    //	Define the numerical methods used in the simulation.
+    //	Note that there may be data dependence on the sequence of constructions.
+    //----------------------------------------------------------------------
+    SimpleDynamics<InitialVelocity> initial_condition(water_block);
+    /** time-space method to detect surface particles. */
+    InteractionWithUpdate<fluid_dynamics::SpatialTemporalFreeSurfaceIdentificationInner>
+        free_surface_indicator(water_body_inner);
+    /** Apply transport velocity formulation. */
+    InteractionDynamics<fluid_dynamics::TransportVelocityCorrectionInner<BulkParticles>> transport_velocity_correction(water_body_inner);
+
+    Dynamics1Level<fluid_dynamics::Integration1stHalfRiemannCorrect> fluid_pressure_relaxation_correct(water_body_inner);
+    Dynamics1Level<fluid_dynamics::Integration2ndHalfRiemann> fluid_density_relaxation(water_body_inner);
+    InteractionWithUpdate<CorrectedConfigurationInner> corrected_configuration_fluid(water_body_inner, 2, 1.0);
+    InteractionWithUpdate<fluid_dynamics::DensitySummationFreeSurfaceInner> fluid_density_by_summation(water_body_inner);
+    SimpleDynamics<TimeStepInitialization> fluid_step_initialization(water_block);
+    ReduceDynamics<fluid_dynamics::AdvectionTimeStepSize> fluid_advection_time_step(water_block, U_max);
+    ReduceDynamics<fluid_dynamics::AcousticTimeStepSize> fluid_acoustic_time_step(water_block);
+    /** We can output a method-specific particle data for debug */
+    water_block.addBodyStateForRecording<Real>("Pressure");
+    water_block.addBodyStateForRecording<int>("Indicator");
+    water_block.addBodyStateForRecording<Real>("PositionDivergence");
+    //----------------------------------------------------------------------
+    //	Define the methods for I/O operations, observations
+    //	and regression tests of the simulation.
+    //----------------------------------------------------------------------
+    BodyStatesRecordingToVtp body_states_recording(io_environment, sph_system.real_bodies_);
+    RestartIO restart_io(io_environment, sph_system.real_bodies_);
+    RegressionTestDynamicTimeWarping<ReducedQuantityRecording<ReduceDynamics<TotalMechanicalEnergy>>>
+        write_water_mechanical_energy(io_environment, water_block);
+    //----------------------------------------------------------------------
+    //	Prepare the simulation with cell linked list, configuration
+    //	and case specified initial condition if necessary.
+    //----------------------------------------------------------------------
+    sph_system.initializeSystemCellLinkedLists();
+    sph_system.initializeSystemConfigurations();
+    initial_condition.exec();
+    //----------------------------------------------------------------------
+    //	Load restart file if necessary.
+    //----------------------------------------------------------------------
+    if (sph_system.RestartStep() != 0)
+    {
+        GlobalStaticVariables::physical_time_ = restart_io.readRestartFiles(sph_system.RestartStep());
+        water_block.updateCellLinkedList();
+        water_body_inner.updateConfiguration();
+    }
+    //----------------------------------------------------------------------
+    //	Setup for time-stepping control
+    //----------------------------------------------------------------------
+    size_t number_of_iterations = sph_system.RestartStep();
+    int screen_output_interval = 100;
+    int observation_sample_interval = screen_output_interval * 2;
+    int restart_output_interval = screen_output_interval * 10;
+    Real end_time = 1.0;
+    Real output_interval = 0.005;
+    //----------------------------------------------------------------------
+    //	Statistics for CPU time
+    //----------------------------------------------------------------------
+    TickCount t1 = TickCount::now();
+    TimeInterval interval;
+    TimeInterval interval_computing_time_step;
+    TimeInterval interval_computing_fluid_pressure_relaxation;
+    TimeInterval interval_updating_configuration;
+    TickCount time_instance;
+    //----------------------------------------------------------------------
+    //	First output before the main loop.
+    //----------------------------------------------------------------------
+    body_states_recording.writeToFile();
+    write_water_mechanical_energy.writeToFile(number_of_iterations);
+    //----------------------------------------------------------------------
+    //	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 < output_interval)
+        {
+            /** outer loop for dual-time criteria time-stepping. */
+            time_instance = TickCount::now();
+            fluid_step_initialization.exec();
+            Real advection_dt = 0.3 *fluid_advection_time_step.exec();
+            free_surface_indicator.exec();
+            fluid_density_by_summation.exec();
+            corrected_configuration_fluid.exec();
+            transport_velocity_correction.exec();
+            interval_computing_time_step += TickCount::now() - time_instance;
+
+            time_instance = TickCount::now();
+            Real relaxation_time = 0.0;
+            Real acoustic_dt = 0.0;
+            while (relaxation_time < advection_dt)
+            {
+                /** inner loop for dual-time criteria time-stepping.  */
+                acoustic_dt = fluid_acoustic_time_step.exec();
+                fluid_pressure_relaxation_correct.exec(acoustic_dt);
+                fluid_density_relaxation.exec(acoustic_dt);
+                relaxation_time += acoustic_dt;
+                integration_time += acoustic_dt;
+                GlobalStaticVariables::physical_time_ += acoustic_dt;
+            }
+            interval_computing_fluid_pressure_relaxation += TickCount::now() - time_instance;
+
+            /** screen output, write body reduced values and restart files  */
+            if (number_of_iterations % screen_output_interval == 0)
+            {
+                std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << "	Time = "
+                          << GlobalStaticVariables::physical_time_
+                          << "	advection_dt = " << advection_dt << "	acoustic_dt = " << acoustic_dt << "\n";
+
+                if (number_of_iterations % restart_output_interval == 0)
+                    restart_io.writeToFile(number_of_iterations);
+            }
+            number_of_iterations++;
+
+            write_water_mechanical_energy.writeToFile(number_of_iterations);
+
+            /** Update cell linked list and configuration. */
+            time_instance = TickCount::now();
+            water_block.updateCellLinkedListWithParticleSort(100);
+            water_body_inner.updateConfiguration();
+            interval_updating_configuration += TickCount::now() - time_instance;
+        }
+
+        body_states_recording.writeToFile();
+        TickCount t2 = 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;
+    std::cout << std::fixed << std::setprecision(9) << "interval_computing_time_step ="
+              << interval_computing_time_step.seconds() << "\n";
+    std::cout << std::fixed << std::setprecision(9) << "interval_computing_fluid_pressure_relaxation = "
+              << interval_computing_fluid_pressure_relaxation.seconds() << "\n";
+    std::cout << std::fixed << std::setprecision(9) << "interval_updating_configuration = "
+              << interval_updating_configuration.seconds() << "\n";
+
+    return 0;
+};

tests/user_examples/test_2d_impact_patch/regression_test_tool/regression_test_tool.py

@@ -0,0 +1,44 @@
+# !/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_standing_wave
+"""
+
+case_name = "test_2d_standing_wave"
+body_name = "WaterBody"
+parameter_name = "TotalMechanicalEnergy"
+body_name_1 = "WaveProbe"
+parameter_name_1 = "FreeSurfaceHeight"
+
+number_of_run_times = 0
+converged = 0
+sphinxsys = SphinxsysRegressionTest(case_name, body_name, parameter_name)
+sphinxsys_1 = SphinxsysRegressionTest(case_name, body_name_1, parameter_name_1)
+
+
+while True:
+    print("Now start a new run......")
+    sphinxsys.run_case()
+    number_of_run_times += 1
+    converged = sphinxsys.read_dat_file()
+    converged_1 = sphinxsys_1.read_dat_file()
+    print("Please note: This is the", number_of_run_times, "run!")
+    if number_of_run_times <= 200:
+        if (converged == "true") and (converged_1 == "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
+        elif converged_1 != "true":
+            print("The tested parameters of", sphinxsys_1.sphinxsys_parameter_name, "are not converged!")
+            continue
+    else:
+        print("It's too many runs but still not converged, please try again!")
+        break

tests/user_examples/test_2d_ration_patch/CMakeLists.txt

@@ -0,0 +1,25 @@
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${CMAKE_SOURCE_DIR}/cmake) # main (top) cmake dir
+
+set(CMAKE_VERBOSE_MAKEFILE on)
+
+STRING( REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR} )
+PROJECT("${CURRENT_FOLDER}")
+
+SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/")
+SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input")
+SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload")
+
+file(MAKE_DIRECTORY ${BUILD_INPUT_PATH})
+file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/regression_test_tool/
+     DESTINATION ${BUILD_INPUT_PATH})
+
+aux_source_directory(. DIR_SRCS)
+add_executable(${PROJECT_NAME} ${EXECUTABLE_OUTPUT_PATH} ${DIR_SRCS} )
+set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}")
+target_link_libraries(${PROJECT_NAME} extra_sources_2d)
+
+add_test(NAME ${PROJECT_NAME} COMMAND ${PROJECT_NAME}
+                 WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})
+add_test(NAME ${PROJECT_NAME}_restart COMMAND ${PROJECT_NAME} --restart_step=4000
+                 WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH})