Xiangyu/check fish swiming case

README.md

@@ -35,6 +35,8 @@ Here, we present several short examples in flow, solid dynamics, fluid structure
 <img src="https://github.com/Xiangyu-Hu/SPHinXsys-public-files/blob/master/videos/half-sphere.gif" height="192px"></a>
 <a href="https://github.com/Xiangyu-Hu/SPHinXsys/blob/master/tests/3d_examples/test_3d_twisting_column/twisting_column.cpp"> 
 <img src="https://github.com/Xiangyu-Hu/SPHinXsys-public-files/blob/master/videos/twisting.gif" height="168px"></a>
+<a href="https://github.com/Xiangyu-Hu/SPHinXsys/blob/master/tests/user_examples/test_2d_flow_stream_around_fish/2d_flow_stream_around_fish.cpp"> 
+<img src="https://github.com/Xiangyu-Hu/SPHinXsys-public-files/blob/master/videos/fish-swimming.gif" height="168px"></a>
 
 ## Journal publications
 

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,163 +57,178 @@ 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<Vecd> createWaterBlockShape()
 {
-	//geometry
-	std::vector<Vecd> 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<Vecd> 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<Vecd> 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<Vecd> 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.
 //----------------------------------------------------------------------
 /** 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<MultiPolygonShape>(outer_boundary, "OuterBoundary");
-		MultiPolygon fish(CreatFishShape(cx, cy, fish_length, fish_shape_resolution));
-		subtract<MultiPolygonShape>(fish);
-	}
+  public:
+    explicit WaterBlock(const std::string &shape_name) : ComplexShape(shape_name)
+    {
+        /** Geometry definition. */
+        MultiPolygon outer_boundary(createWaterBlockShape());
+        add<MultiPolygonShape>(outer_boundary, "OuterBoundary");
+        MultiPolygon fish(CreatFishShape(cx, cy, fish_length, fish_shape_resolution));
+        subtract<MultiPolygonShape>(fish);
+    }
 };
 
 /** Case dependent inflow boundary condition. */
 struct FreeStreamVelocity
 {
-	Real u_ref_, t_ref_;
-
-	template <class BoundaryConditionType>
-	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 <class BoundaryConditionType>
+    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<ActiveModelSolid>(rho0_s, Youngs_modulus1, poisson);
-		add<SaintVenantKirchhoffSolid>(rho0_s, Youngs_modulus2, poisson);
-		add<SaintVenantKirchhoffSolid>(rho0_s, Youngs_modulus3, poisson);
-	};
+  public:
+    SolidBodyMaterial() : CompositeMaterial(rho0_s)
+    {
+        add<ActiveModelSolid>(rho0_s, Youngs_modulus1, poisson);
+        add<SaintVenantKirchhoffSolid>(rho0_s, Youngs_modulus2, poisson);
+        add<SaintVenantKirchhoffSolid>(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<SolidParticles*>(&solid_body.getBaseParticles())),
-		materail_id_(*solid_particles_->getVariableByName<int>("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<int>& materail_id_;
-	StdLargeVec<Vecd>& pos0_;
-	StdLargeVec<Matd> active_strain_;
+  public:
+    explicit ImposingActiveStrain(SolidBody &solid_body)
+        : solid_dynamics::ElasticDynamicsInitialCondition(solid_body),
+          material_id_(*particles_->getVariableByName<int>("MaterialID")),
+          pos0_(*particles_->getVariableByName<Vecd>("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<int> &material_id_;
+    StdLargeVec<Vecd> &pos0_;
+    StdLargeVec<Matd> active_strain_;
 };

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)