Optimize contact

apps/twop/example/input_0.yaml

@@ -124,5 +124,5 @@ Output:
   Perform_FE_Out: false
   Perform_Out: true
   Test_Output_Interval: 4000
-  Debug: 1
+  Debug: 3
   Tag_PP: 0

apps/twop/main.cpp

@@ -78,8 +78,8 @@ class Model : public model::DEMModel {
     // check if init() successfully created quadrature data which we need for postprocessing
     {
       size_t totalQuadPoints = 0;
-      for (auto &p: d_particles) {
-        const auto *particle_mesh_p = p->d_rp_p->getMeshP();
+      for (auto &p: d_particlesListTypeAll) {
+        const auto &particle_mesh_p = p->d_rp_p->getMeshP();
 
         // get Quadrature
         fe::BaseElem *elem;

examples/PeriDEM/attrition_tests/sim1_multi_particle_circ_tri_drum_with_rotating_cylinder_with_protrusion/problem_setup.py

@@ -903,6 +903,11 @@ def create_input_file(inp_dir, pp_tag):
   ## wall rotation rate
   rotation_rate = -20. * np.pi
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
+
   ### ---------------------------------------------------------------- ###
   # generate mesh and particle location data
   ### ---------------------------------------------------------------- ###
@@ -1119,8 +1124,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

examples/PeriDEM/attrition_tests/sim1_multi_particle_circ_tri_drum_with_rotating_cylinder_with_protrusion_increased_Kn_factor/problem_setup.py

@@ -903,6 +903,11 @@ def create_input_file(inp_dir, pp_tag):
   ## wall rotation rate
   rotation_rate = -20. * np.pi
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
+
   ### ---------------------------------------------------------------- ###
   # generate mesh and particle location data
   ### ---------------------------------------------------------------- ###
@@ -1119,8 +1124,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

examples/PeriDEM/attrition_tests/sim2_multi_particle_circ_tri_drum_hex_with_rotating_cylinder_with_protrusion/problem_setup.py

@@ -895,6 +895,11 @@ def create_input_file(inp_dir, pp_tag):
   ## wall rotation rate
   rotation_rate = -20. * np.pi
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
+
   ### ---------------------------------------------------------------- ###
   # generate mesh and particle location data
   ### ---------------------------------------------------------------- ###
@@ -1080,8 +1085,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

examples/PeriDEM/attrition_tests/sim2_multi_particle_circ_tri_drum_hex_with_rotating_cylinder_with_protrusion_thin_container_and_change_rotation_rate/problem_setup.py

@@ -894,6 +894,11 @@ def create_input_file(inp_dir, pp_tag):
   wall_rotation_rate = -40. * np.pi
   wall_rotation_center = [-0.2*R_in, 0.2*R_in, 0.]
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
+
   ### ---------------------------------------------------------------- ###
   # generate mesh and particle location data
   ### ---------------------------------------------------------------- ###
@@ -1081,8 +1086,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

examples/PeriDEM/attrition_tests/sim3_multi_particle_circ_tri_drum_hex_with_rotating_rectangle_container_with_protrusion_and_opening/problem_setup.py

@@ -1103,6 +1103,11 @@ def create_input_file(inp_dir, pp_tag):
   wall_rotation_rate = -40. * np.pi
   wall_rotation_center = [0.15*Lin, 0.15*Win, 0.]
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
+
   ### ---------------------------------------------------------------- ###
   # generate mesh and particle location data
   ### ---------------------------------------------------------------- ###
@@ -1327,8 +1332,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

examples/PeriDEM/attrition_tests/sim4_multi_particle_circ_tri_drum_hex_with_rotating_rectangle_container_with_protrusion_and_opening_within_bigger_rectangle_container/problem_setup.py

@@ -392,6 +392,7 @@ def particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R
 
     pcount = 0
     count = 0
+    select_count = 0
     while pcount < N1 and count < 100*N1:
       if count%N1 == 0:
         print('small particle iter = ', count)
@@ -400,9 +401,10 @@ def particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R
       r = R1 + np.random.uniform(-0.1 * R1, 0.1 * R1)
       x = center[0] + np.random.uniform(-0.5*Lin + R1 + padding, 0.5*Lin - R1- padding)
       y = np.random.uniform(-0.5*Lin + R1 + padding, max_y - R1 - padding)
-      p_zone = np.random.choice(id_choices1, size=1)[0]
-      if np.random.uniform(0., 1.) > 0.25:
-        p_zone = np.random.choice(id_choices1, size=1)[0]
+      # p_zone = np.random.choice(id_choices1, size=1)[0]
+      # if np.random.uniform(0., 1.) > 0.25:
+      #   p_zone = np.random.choice(id_choices1, size=1)[0]
+      p_zone = id_choices1[select_count % len(id_choices1)]
       p = [p_zone, x,y,center[2], r]
 
       # check if it collides of existing particles
@@ -411,6 +413,8 @@ def particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R1, R
       if pintersect == False:
         particles.append(p)
         pcount += 1
+        select_count += 1
+
       count +=1
 
   inpf = open(inp_dir + 'particle_locations_' + str(pp_tag) + '.csv','w')
@@ -1025,17 +1029,19 @@ def create_input_file(inp_dir, pp_tag):
 
   L = 0.04
   W = 0.05
-  Lin = L - horizon
-  Win = W - horizon
+  Lin = L - 3*horizon
+  Win = W - 3*horizon
 
   ## wall rotation rate
-  wall_rotation_rate = -20. * np.pi
+  wall_rotation_rate = -30. * np.pi
   wall_rotation_center = [0.15*Lin, 0.15*Win, 0.]
 
   ## rotation of outer rectangle will be within following rectangle
   wall_rot_max_dist_from_rot_center = L - (0.5*L - wall_rotation_center[0])
   if wall_rot_max_dist_from_rot_center < W - (0.5*W - wall_rotation_center[1]):
     wall_rot_max_dist_from_rot_center = W - (0.5*W - wall_rotation_center[1])
+
+  wall_rot_max_dist_from_rot_center = wall_rot_max_dist_from_rot_center + 0.05*L
 
   # two corner points
   wall_rot_working_rect_relative_corner_point_1 = [-wall_rot_max_dist_from_rot_center, -wall_rot_max_dist_from_rot_center, 0.]
@@ -1054,8 +1060,8 @@ def create_input_file(inp_dir, pp_tag):
   # right-top corner
   base_rect_inner_corner_point_2 = [wall_rot_working_rect_corner_point_2[0] + 0.2*L, wall_rot_working_rect_corner_point_2[1] + 0.2*L, wall_rot_working_rect_corner_point_2[2]]
 
-  base_rect_x_thickness = 0.5*horizon
-  base_rect_y_thickness = 0.5*horizon
+  base_rect_x_thickness = 3*horizon
+  base_rect_y_thickness = 3*horizon
 
   # annulus rectangle requires 12 parameters (first 6 are for outer rectangle and rest are inner/void rectangle). Six parameters are coordinates of two corner points
   base_rect_innr_pts = get_ref_rect_points(get_center(base_rect_inner_corner_point_1, base_rect_inner_corner_point_1), base_rect_inner_corner_point_2[0] - base_rect_inner_corner_point_1[0], base_rect_inner_corner_point_2[1] - base_rect_inner_corner_point_1[1])
@@ -1090,9 +1096,9 @@ def create_input_file(inp_dir, pp_tag):
 
   ## time 
   final_time = 0.1
-  num_steps = 4000000
-  # final_time = 0.00002
-  # num_steps = 2
+  num_steps = 1000000
+  # final_time = 0.1
+  # num_steps = 40000
   num_outputs = 400
   dt_out_n = num_steps / num_outputs
   test_dt_out_n = dt_out_n / 100
@@ -1102,7 +1108,7 @@ def create_input_file(inp_dir, pp_tag):
   poisson = 0.25
 
   rho_wall = 1200.
-  K_wall = 1.e+6
+  K_wall = 1.e+4
   E_wall = get_E(K_wall, poisson)
   G_wall = get_G(E_wall, poisson)
   Gc_wall = 100.
@@ -1114,7 +1120,7 @@ def create_input_file(inp_dir, pp_tag):
   Gc_large = Gc_wall
 
   rho_small = 1200.
-  K_small = 2.e+5
+  K_small = 2.e+3
   E_small = get_E(K_small, poisson)
   G_small = get_G(E_small, poisson)
   Gc_small = 50.
@@ -1151,6 +1157,10 @@ def create_input_file(inp_dir, pp_tag):
   gravity_active = True
   gravity = [0., -10., 0.]
 
+  ## neighbor search details
+  neigh_search_factor = 10.
+  neigh_search_interval = 40
+  neigh_search_criteria = "simple_all"
 
 
   ### ---------------------------------------------------------------- ###
@@ -1161,7 +1171,7 @@ def create_input_file(inp_dir, pp_tag):
   padding = 1.1 * R_contact_factor * mesh_size
   max_y = 0.5*Lin - 3*mesh_size
   # number of particles of small and large sizes
-  N1, N2 = 30, 8
+  N1, N2 = 40, 12
   id_choices1 = [0, 3, 1, 2]
   id_choices2 = [4, 5, 6, 7]
   num_particles_zones_1_to_8, particles_zones_1_to_8 = particle_locations(inp_dir, pp_tag, center, padding, max_y, mesh_size, R_small, R_large, id_choices1, id_choices2, N1, N2, Lin, Win, z_coord = 0., add_orient = True)
@@ -1378,8 +1388,9 @@ def create_input_file(inp_dir, pp_tag):
 
   # Neighbor info
   inpf.write("Neighbor:\n")
-  inpf.write("  Update_Criteria: simple_all\n")
-  inpf.write("  Search_Factor: 5.0\n")
+  inpf.write("  Update_Criteria: %s\n" % (neigh_search_criteria))
+  inpf.write("  Search_Factor: %4.e\n" % (neigh_search_factor))
+  inpf.write("  Search_Interval: %d\n" % (neigh_search_interval))
 
   # Material info
   inpf.write("Material:\n")

src/fe/mesh.cpp

@@ -263,6 +263,7 @@ void fe::Mesh::computeVol() {
                  "element mesh as the element-node connectivity data is "
                  "invalid."
               << std::endl;
+    exit(EXIT_FAILURE);
   }
 
   if (false) {

src/geometry/geometryUtil.cpp

@@ -12,7 +12,7 @@
 #include "nsearch/nsearch.h"
 #include <fmt/format.h>
 
-typedef nsearch::NFlannSearchKd NSearch;
+typedef nsearch::NFlannSearchKd<3> NSearch;
 
 void geometry::computeNonlocalNeighborhood(const std::vector<util::Point> &nodes,
                                       double horizon,

src/inp/input.cpp

@@ -366,6 +366,10 @@ void inp::Input::setParticleDeck() {
       d_particleDeck_p->d_pNeighDeck.d_sFactor =
           ne["Search_Factor"].as<double>();
 
+    if (ne["Search_Interval"])
+      d_particleDeck_p->d_pNeighDeck.d_neighUpdateInterval =
+              ne["Search_Interval"].as<size_t>();
+
     for (auto p : d_particleDeck_p->d_particleZones) {
 
       // if particle is actually a wall then we do it differently
@@ -675,7 +679,7 @@ void inp::Input::setContactDeck() {
     for (size_t zz = z; zz <= n; zz++) {
       std::string read_zone = "Zone_" + std::to_string(z) + std::to_string(zz);
 
-      if (d_outputDeck_p->d_debug > 1)
+      if (d_outputDeck_p->d_debug > 3)
         std::cout << "Processing pair (z, zz) = (" << z << ", " << zz << ").\n";
 
       auto e = config["Contact"][read_zone];
@@ -708,7 +712,7 @@ void inp::Input::setContactDeck() {
         inp::ContactPairDeck cd =
             d_contactDeck_p->d_data[zone_pair[0] - 1][zone_pair[1] - 1];
 
-        if (d_outputDeck_p->d_debug > 1) {
+        if (d_outputDeck_p->d_debug > 3) {
           std::cout << "Processing pair (z, zz) = (" << z << ", " << zz << "). "
                     << "Copying contact data from (n, m) = (" << zone_pair[0]
                     << ", " << zone_pair[1] << ").\n";

src/inp/pdecks/pNeighborDeck.h

@@ -27,14 +27,22 @@ struct PNeighborDeck {
    * radius of particle) */
   double d_sFactor;
 
-  /*! @brief Neighbor search tolerance */
+  /*!
+   * @brief Neighbor search tolerance
+   */
   double d_sTol;
 
+  /*! @brief Neighbor update time interval (for contact) */
+  size_t d_neighUpdateInterval;
+
   /*!
    * @brief Constructor
    */
   PNeighborDeck()
-      : d_updateCriteria("max_distance_travel"), d_sFactor(5.), d_sTol(0.){};
+      : d_updateCriteria("simple_all"),
+        d_sFactor(1.),
+        d_sTol(0.),
+        d_neighUpdateInterval(1) {};
 
   /*!
    * @brief Returns the string containing printable information about the object
@@ -51,6 +59,7 @@ struct PNeighborDeck {
     oss << tabS << "Update criteria  = " << d_updateCriteria << std::endl;
     oss << tabS << "Search factor = " << d_sFactor << std::endl;
     oss << tabS << "Search tolerance = " << d_sTol << std::endl;
+    oss << tabS << "Search update interval = " << d_neighUpdateInterval << std::endl;
     oss << tabS << std::endl;
 
     return oss.str();

src/loading/particleFLoading.cpp

@@ -149,7 +149,7 @@ void loading::ParticleFLoading::apply(const double &time,
 
                   auto force_i = util::Point();
                   for (auto d : bc.d_direction) {
-                    force_i[d] = fmax;
+                    force_i[d-1] = fmax;
                   }
 
                   // add force

src/loading/particleULoading.cpp

@@ -182,8 +182,8 @@ void loading::ParticleULoading::apply(const double &time,
                   auto u_i = util::Point();
                   auto v_i = util::Point();
                   for (auto d : bc.d_direction) {
-                    u_i[d] = du;
-                    v_i[d] = dv;
+                    u_i[d-1] = du;
+                    v_i[d-1] = dv;
                   }
 
                   if (bc.d_timeFnType == "rotation") {