fixes storing of colliding particle pairs

cpp/particles/interpolation/particles_field_computer.hpp

@@ -123,15 +123,15 @@ public:
 
     real_number get_norm_pos_in_cell(const real_number in_pos, const int idx_pos) const {
         const real_number shifted_pos = in_pos - spatial_box_offset[idx_pos];
-        DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, shifted_pos = %g\n", shifted_pos);
+        //DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, shifted_pos = %g\n", shifted_pos);
         const int nb_box_repeat = int(floor(shifted_pos/spatial_box_width[idx_pos]));
-        DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, nb_box_repeat = %d\n", nb_box_repeat);
+        //DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, nb_box_repeat = %d\n", nb_box_repeat);
         const real_number pos_in_box = shifted_pos - nb_box_repeat*spatial_box_width[idx_pos];
-        DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, pos_in_box = %g\n", pos_in_box);
+        //DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, pos_in_box = %g\n", pos_in_box);
         const int cell_idx = int(floor(pos_in_box/box_step_width[idx_pos]));
-        DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, cell_idx = %d\n", cell_idx);
+        //DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, cell_idx = %d\n", cell_idx);
         const real_number pos_in_cell = (pos_in_box - cell_idx*box_step_width[idx_pos]) / box_step_width[idx_pos];
-        DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, pos_in_cell = %g\n", pos_in_cell);
+        //DEBUG_MSG("particles_field_computer::get_norm_pos_in_cell, pos_in_cell = %g\n", pos_in_cell);
         assert(0 <= pos_in_cell && pos_in_cell < 1);
         return pos_in_cell;
     }

cpp/particles/p2p/p2p_computer.hpp

@@ -69,6 +69,9 @@ public:
         }
     }
 
+    /**
+     * NOTE: this is called only ONCE for each pair of interacting particles.
+     */
     template <int size_particle_positions, int size_particle_rhs>
     void compute_interaction(const partsize_t /*idx_part1*/,
                              const real_number pos_part1[],

cpp/particles/p2p/p2p_computer_empty.hpp

@@ -39,6 +39,9 @@ public:
     void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
     }
 
+    /**
+     * NOTE: this is called only ONCE for each pair of interacting particles.
+     */
     template <int size_particle_positions, int size_particle_rhs>
     void compute_interaction(const partsize_t /*idx_part1*/,
                              const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],

cpp/particles/p2p/p2p_ghost_collisions.hpp

@@ -57,18 +57,17 @@ class p2p_ghost_collisions
         bool isActive;
         void add_colliding_pair(partsize_t idx_part1, partsize_t idx_part2)
         {
-            //DEBUG_MSG("thread %d, compute_interaction %ld, computer %ld, with collision_pairs_local size %d\n",
-            //        omp_get_thread_num(),
-            //        &this->collision_pairs_local,
-            //        this,
-            //        this->collision_pairs_local.size());
-            if (idx_part1 < idx_part2) // only store unique pairs of colliding particles
+            // store colliding particle ids in order, to be able to identify pairs more easily
+            assert(idx_part1!=idx_part2);
+            partsize_t idx_part_small = idx_part1;
+            partsize_t idx_part_big = idx_part2;
+            if (idx_part1 > idx_part2)
             {
-                this->collision_pairs_local.push_back(idx_part1);
-                this->collision_pairs_local.push_back(idx_part2);
+                idx_part_small = idx_part2;
+                idx_part_big = idx_part1;
             }
-            //DEBUG_MSG("inside compute interaction idx_part1 = %ld and idx_part2 = %ld\n", idx_part1, idx_part2);
-            assert(idx_part1!=idx_part2);
+            this->collision_pairs_local.push_back(idx_part_small);
+            this->collision_pairs_local.push_back(idx_part_big);
         }
     protected:
         bool synchronisation;
@@ -185,6 +184,9 @@ public:
     }
 
 
+    /**
+     * NOTE: this is called only ONCE for each pair of interacting particles.
+     */
     template <int size_particle_positions, int size_particle_rhs>
     void compute_interaction(const partsize_t idx_part1,
                              const real_number pos_part1[],