Merge pull request #420 from DongWuTUM/shell_particle_factor_scaling

Scale the factor of ShellMidSurfaceBounding

src/for_2D_build/particle_generator/particle_generator_lattice_supplementary.cpp

@@ -54,6 +54,11 @@ void ThickSurfaceParticleGeneratorLattice::initializeGeometricVariables()
     Real interval = (Real)planned_number_of_particles_ / (Real)all_cells_;
     if (interval <= 0)
         interval = 1; // It has to be lager than 0.
+
+    // initialize a uniform distribution between 0 (inclusive) and 1 (exclusive)
+    std::mt19937_64 rng;
+    std::uniform_real_distribution<Real> unif(0, 1);
+
     // Add a particle in each interval, randomly. We will skip the last intervals if we already reach the number of particles
     for (int i = 0; i < number_of_lattices[0]; ++i)
         for (int j = 0; j < number_of_lattices[1]; ++j)
@@ -63,7 +68,7 @@ void ThickSurfaceParticleGeneratorLattice::initializeGeometricVariables()
             {
                 if (body_shape_.checkContain(particle_position))
                 {
-                    Real random_real = (Real)rand() / (RAND_MAX);
+                    Real random_real = unif(rng);
                     // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles
                     if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_)
                     {

src/for_3D_build/particle_generator/particle_generator_lattice_supplementary.cpp

@@ -52,10 +52,15 @@ void ThickSurfaceParticleGeneratorLattice::initializeGeometricVariables()
     Real number_of_particles = total_volume_ / avg_particle_volume_ + 0.5;
     planned_number_of_particles_ = int(number_of_particles);
 
+    // initialize a uniform distribution between 0 (inclusive) and 1 (exclusive)
+    std::mt19937_64 rng;
+    std::uniform_real_distribution<Real> unif(0, 1);
+
     // Calculate the interval based on the number of particles.
     Real interval = planned_number_of_particles_ / (all_cells_ + TinyReal);
     if (interval <= 0)
         interval = 1; // It has to be lager than 0.
+
     // Add a particle in each interval, randomly. We will skip the last intervals if we already reach the number of particles.
     for (int i = 0; i < number_of_lattices[0]; ++i)
         for (int j = 0; j < number_of_lattices[1]; ++j)
@@ -66,7 +71,7 @@ void ThickSurfaceParticleGeneratorLattice::initializeGeometricVariables()
                 {
                     if (body_shape_.checkContain(particle_position))
                     {
-                        Real random_real = (Real)rand() / (RAND_MAX);
+                        Real random_real = unif(rng);
                         // If the random_real is smaller than the interval, add a particle, only if we haven't reached the max. number of particles.
                         if (random_real <= interval && base_particles_.total_real_particles_ < planned_number_of_particles_)
                         {

src/shared/particle_dynamics/relax_dynamics/relax_dynamics.cpp

@@ -159,19 +159,17 @@ void RelaxationStepComplex::exec(Real dt)
 }
 //=================================================================================================//
 ShellMidSurfaceBounding::
-    ShellMidSurfaceBounding(NearShapeSurface &body_part, BaseInnerRelation &inner_relation,
-                            Real thickness, Real level_set_refinement_ratio)
+    ShellMidSurfaceBounding(NearShapeSurface &body_part, BaseInnerRelation &inner_relation)
     : BaseLocalDynamics<BodyPartByCell>(body_part), RelaxDataDelegateInner(inner_relation),
       pos_(particles_->pos_), constrained_distance_(0.5 * sph_body_.sph_adaptation_->MinimumSpacing()),
       particle_spacing_ref_(sph_body_.sph_adaptation_->MinimumSpacing()),
-      thickness_(thickness), level_set_refinement_ratio_(level_set_refinement_ratio),
       level_set_shape_(DynamicCast<LevelSetShape>(this, sph_body_.body_shape_)) {}
 //=================================================================================================//
 void ShellMidSurfaceBounding::update(size_t index_i, Real dt)
 {
     Vecd none_normalized_normal = level_set_shape_->findLevelSetGradient(pos_[index_i]);
     Vecd normal = level_set_shape_->findNormalDirection(pos_[index_i]);
-    Real factor = none_normalized_normal.squaredNorm() / level_set_refinement_ratio_;
+    Real factor = 0.2 * none_normalized_normal.norm();
     pos_[index_i] -= factor * constrained_distance_ * normal;
 }
 //=================================================================================================//
@@ -292,12 +290,10 @@ void ShellNormalDirectionPrediction::SmoothingNormal::update(size_t index_i, Rea
 }
 //=================================================================================================//
 ShellRelaxationStepInner::
-    ShellRelaxationStepInner(BaseInnerRelation &inner_relation, Real thickness,
-                             Real level_set_refinement_ratio, bool level_set_correction)
+    ShellRelaxationStepInner(BaseInnerRelation &inner_relation, bool level_set_correction)
     : RelaxationStepInner(inner_relation, level_set_correction),
       update_shell_particle_position_(*real_body_),
-      mid_surface_bounding_(near_shape_surface_, inner_relation,
-                            thickness, level_set_refinement_ratio) {}
+      mid_surface_bounding_(near_shape_surface_, inner_relation) {}
 //=================================================================================================//
 void ShellRelaxationStepInner::exec(Real ite_p)
 {

src/shared/particle_dynamics/relax_dynamics/relax_dynamics.h

@@ -302,15 +302,14 @@ class ShellMidSurfaceBounding : public BaseLocalDynamics<BodyPartByCell>,
                                 public RelaxDataDelegateInner
 {
   public:
-    ShellMidSurfaceBounding(NearShapeSurface &body_part, BaseInnerRelation &inner_relation,
-                            Real thickness, Real level_set_refinement_ratio);
+    ShellMidSurfaceBounding(NearShapeSurface &body_part, BaseInnerRelation &inner_relation);
     virtual ~ShellMidSurfaceBounding(){};
     void update(size_t index_i, Real dt = 0.0);
 
   protected:
     StdLargeVec<Vecd> &pos_;
     Real constrained_distance_;
-    Real particle_spacing_ref_, thickness_, level_set_refinement_ratio_;
+    Real particle_spacing_ref_;
     LevelSetShape *level_set_shape_;
 };
 
@@ -440,8 +439,7 @@ class ShellNormalDirectionPrediction : public BaseDynamics<void>
 class ShellRelaxationStepInner : public RelaxationStepInner
 {
   public:
-    explicit ShellRelaxationStepInner(BaseInnerRelation &inner_relation, Real thickness,
-                                      Real level_set_refinement_ratio, bool level_set_correction = false);
+    explicit ShellRelaxationStepInner(BaseInnerRelation &inner_relation, bool level_set_correction = false);
     virtual ~ShellRelaxationStepInner(){};
 
     SimpleDynamics<UpdateParticlePosition> update_shell_particle_position_;

tests/2d_examples/test_2d_ball_shell_collision/ball_shell_collision.cpp

@@ -112,8 +112,7 @@ int main(int ac, char *av[])
         //	Define the methods for particle relaxation for wall boundary.
         //----------------------------------------------------------------------
         SimpleDynamics<RandomizeParticlePosition> wall_boundary_random_particles(wall_boundary);
-        relax_dynamics::ShellRelaxationStepInner
-            relaxation_step_wall_boundary_inner(wall_boundary_inner, thickness, level_set_refinement_ratio);
+        relax_dynamics::ShellRelaxationStepInner relaxation_step_wall_boundary_inner(wall_boundary_inner);
         relax_dynamics::ShellNormalDirectionPrediction shell_normal_prediction(wall_boundary_inner, thickness, cos(Pi / 3.75));
         wall_boundary.addBodyStateForRecording<int>("UpdatedIndicator");
         //----------------------------------------------------------------------

tests/2d_examples/test_2d_shell_beam_collision/src/shell_beam_collision.cpp

@@ -85,7 +85,7 @@ int main(int ac, char *av[])
     /** Tag for running particle relaxation for the initially body-fitted distribution */
     sph_system.setRunParticleRelaxation(false);
     /** Tag for starting with relaxed body-fitted particles distribution */
-    sph_system.setReloadParticles(false);
+    sph_system.setReloadParticles(true);
     sph_system.handleCommandlineOptions(ac, av);
     IOEnvironment io_environment(sph_system);
     //----------------------------------------------------------------------
@@ -138,8 +138,7 @@ int main(int ac, char *av[])
         //	Define the methods for particle relaxation for wall boundary.
         //----------------------------------------------------------------------
         SimpleDynamics<RandomizeParticlePosition> shell_random_particles(shell);
-        relax_dynamics::ShellRelaxationStepInner
-            relaxation_step_shell_inner(shell_inner, thickness, level_set_refinement_ratio);
+        relax_dynamics::ShellRelaxationStepInner relaxation_step_shell_inner(shell_inner);
         relax_dynamics::ShellNormalDirectionPrediction shell_normal_prediction(shell_inner, thickness, cos(Pi / 3.75));
         shell.addBodyStateForRecording<int>("UpdatedIndicator");
         //----------------------------------------------------------------------

tests/2d_examples/test_2d_shell_particle_relaxation/test_2d_shell_particle_relaxation.cpp

@@ -63,8 +63,7 @@ int main(int ac, char *av[])
     /** Random reset the particle position. */
     SimpleDynamics<RandomizeParticlePosition> random_pipe_body_particles(pipe_body);
     /** A  Physics relaxation step. */
-    relax_dynamics::ShellRelaxationStepInner
-        relaxation_step_pipe_body_inner(pipe_body_inner, thickness, level_set_refinement_ratio);
+    relax_dynamics::ShellRelaxationStepInner relaxation_step_pipe_body_inner(pipe_body_inner);
     relax_dynamics::ShellNormalDirectionPrediction shell_normal_prediction(pipe_body_inner, thickness);
     pipe_body.addBodyStateForRecording<int>("UpdatedIndicator");
     /**

tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp

@@ -14,13 +14,13 @@ using namespace SPH;
 
 Real to_rad(Real angle) { return angle * Pi / 180; }
 
-void relax_shell(RealBody &plate_body, Real thickness, Real level_set_refinement_ratio)
+void relax_shell(RealBody &plate_body, Real thickness)
 {
     // BUG: apparently only works if dp > thickness, otherwise ShellNormalDirectionPrediction::correctNormalDirection() throws error
 
     InnerRelation imported_model_inner(plate_body);
     SimpleDynamics<RandomizeParticlePosition> random_imported_model_particles(plate_body);
-    relax_dynamics::ShellRelaxationStepInner relaxation_step_inner(imported_model_inner, thickness, level_set_refinement_ratio);
+    relax_dynamics::ShellRelaxationStepInner relaxation_step_inner(imported_model_inner);
     relax_dynamics::ShellNormalDirectionPrediction shell_normal_prediction(imported_model_inner, thickness);
     //----------------------------------------------------------------------
     //	Particle relaxation starts here.

tests/3d_examples/test_3d_shell_particle_relaxation/test_3d_shell_particle_relaxation.cpp

@@ -73,7 +73,7 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     SimpleDynamics<RandomizeParticlePosition> random_imported_model_particles(imported_model);
     /** A  Physics relaxation step. */
-    relax_dynamics::ShellRelaxationStepInner relaxation_step_inner(imported_model_inner, thickness, level_set_refinement_ratio);
+    relax_dynamics::ShellRelaxationStepInner relaxation_step_inner(imported_model_inner);
     relax_dynamics::ShellNormalDirectionPrediction shell_normal_prediction(imported_model_inner, thickness);
     //----------------------------------------------------------------------
     //	Particle relaxation starts here.