Fix zero pos particles (#31)

* change pressure BC to be define by body part

* SpringNormalOnSurfaceParticles with body part

* option to use inner or outer surface for spring

* zero position of particles fixed

SPHINXsys/src/shared/particle_dynamics/base_particle_dynamics.h

@@ -174,6 +174,7 @@ namespace SPH
 			sorted_id_(body_->base_particles_->sorted_id_),
 			unsorted_id_(body_->base_particles_->unsorted_id_) {};
 		virtual ~DataDelegateSimple() {};
+		ParticlesType* GetParticles(){ return particles_; };
 	protected:
 		BodyType* body_;
 		ParticlesType* particles_;

SPHINXsys/src/shared/particle_dynamics/relax_dynamics/relax_dynamics.cpp

@@ -231,6 +231,9 @@ namespace SPH
 			Real dt_square = get_time_step_square_.parallel_exec();
 			update_solid_particle_position_.parallel_exec(dt_square);
 			surface_bounding_.parallel_exec();
+
+			// copy the updated position at the end of the relaxation step to avoid bugs
+			std::copy(GetParticles()->pos_n_.begin(), GetParticles()->pos_n_.end(), GetParticles()->pos_0_.begin());
 		}
 		//=================================================================================================//
 	}

SPHINXsys/src/shared/particle_dynamics/relax_dynamics/relax_dynamics.h

@@ -241,11 +241,11 @@ namespace SPH
 		* @class SolidRelaxationStepInner
 		* @brief carry out particle relaxation step of particles within the body
 		*/
-		class SolidRelaxationStepInner : public RelaxationStepInner
+		class SolidRelaxationStepInner : public RelaxationStepInner, public solid_dynamics::SolidDataSimple
 		{
 		public:
 			explicit SolidRelaxationStepInner(BaseBodyRelationInner* body_inner_relation, bool level_set_correction = false) :
-				RelaxationStepInner(body_inner_relation, level_set_correction),
+				RelaxationStepInner(body_inner_relation, level_set_correction), solid_dynamics::SolidDataSimple(body_inner_relation->sph_body_),
 				update_solid_particle_position_(real_body_) {};
 			virtual ~SolidRelaxationStepInner() {};
 
@@ -254,7 +254,6 @@ namespace SPH
 			virtual void exec(Real dt = 0.0) override;
 			virtual void parallel_exec(Real dt = 0.0) override;
 		};
-
 	}
 }
 #endif //RELAX_DYNAMICS_H

SPHINXsys/src/shared/particle_dynamics/solid_dynamics/solid_dynamics.cpp

@@ -618,8 +618,8 @@ namespace SPH
 		}
 		//=================================================================================================//
 		SpringNormalOnSurfaceParticles::
-			SpringNormalOnSurfaceParticles(SolidBody *body, Vecd source_point, Real stiffness, Real damping_ratio)
-			: ParticleDynamicsSimple(body), SolidDataSimple(body),
+			SpringNormalOnSurfaceParticles(SolidBody *body, BodyPartByParticle* body_part, bool outer_surface, Vecd source_point, Real stiffness, Real damping_ratio)
+			: PartSimpleDynamicsByParticle(body, body_part), SolidDataSimple(body),
 			  pos_n_(particles_->pos_n_),
 			  pos_0_(particles_->pos_0_),
 			  n_(particles_->n_),
@@ -649,8 +649,16 @@ namespace SPH
 				// get the cos of the angle between the vector and the normal
 				Real cos_teta = getAngleBetweenTwoVectors (vector_to_particle, normal);
 
+				// if outer surface, the normals close an angle greater than 90°
+				// if the angle is greater than 90°, we apply the spring force to the surface particle
+				Real epsilon = 1e-6; // to ignore exactly perpendicular surfaces
+				if (outer_surface && cos_teta < -epsilon)
+				{
+					apply_spring_force_to_particle_[particle_i] = true;
+				}
+				// if not outer surface, it's inner surface, meaning the normals close an angle smaller than 90°
 				// if the angle is less than 90°, we apply the spring force to the surface particle
-				if (cos_teta > 1e-3)
+				if (!outer_surface && cos_teta > epsilon)
 				{
 					apply_spring_force_to_particle_[particle_i] = true;
 				}
@@ -777,8 +785,8 @@ namespace SPH
 		}
 		//=================================================================================================//
 		SurfacePressureFromSource::
-			SurfacePressureFromSource(SPHBody* body, Vecd source_point, StdVec<array<Real, 2>> pressure_over_time)
-			: ParticleDynamicsSimple(body), SolidDataSimple(body),
+			SurfacePressureFromSource(SPHBody* body, BodyPartByParticle* body_part, Vecd source_point, StdVec<array<Real, 2>> pressure_over_time)
+			: PartSimpleDynamicsByParticle(body, body_part), SolidDataSimple(body),
 			  pos_0_(particles_->pos_0_),
 			  n_(particles_->n_),
 			  dvel_dt_prior_(particles_->dvel_dt_prior_),

SPHINXsys/src/shared/particle_dynamics/solid_dynamics/solid_dynamics.h

@@ -411,10 +411,10 @@ namespace SPH
 		* Only for 3D applications
 		*/
 		class SpringNormalOnSurfaceParticles 
-			: public ParticleDynamicsSimple, public SolidDataSimple
+			: public PartSimpleDynamicsByParticle, public SolidDataSimple
 		{
 		public:
-			SpringNormalOnSurfaceParticles(SolidBody* body, Vecd source_point, Real stiffness, Real damping_ratio = 0.05);
+			SpringNormalOnSurfaceParticles(SolidBody* body, BodyPartByParticle* body_part, bool outer_surface, Vecd source_point, Real stiffness, Real damping_ratio = 0.05);
 			~SpringNormalOnSurfaceParticles();
 
 			StdLargeVec<bool>& GetApplySpringForceToParticle(){ return apply_spring_force_to_particle_; }
@@ -472,10 +472,10 @@ namespace SPH
 		* @brief SurfacePressureFromSource, applies pressure on the surface particles coming from a source point
 		*/
 		class SurfacePressureFromSource :
-			public ParticleDynamicsSimple, public SolidDataSimple
+			public PartSimpleDynamicsByParticle, public SolidDataSimple
 		{
 		public:
-			SurfacePressureFromSource(SPHBody* body, Vecd source_point, StdVec<array<Real, 2>> pressure_over_time);
+			SurfacePressureFromSource(SPHBody* body, BodyPartByParticle* body_part, Vecd source_point, StdVec<array<Real, 2>> pressure_over_time);
 			virtual ~SurfacePressureFromSource() {};
 
 			StdLargeVec<bool>& GetApplyPressureToParticle(){ return apply_pressure_to_particle_; }