Update method to find neighbors in grid

bfps/cpp/particles/p2p_distr_mpi.hpp

@@ -246,6 +246,12 @@ public:
                                                                               particles_positions[(idxPart)*size_particle_positions + IDX_Z]);
                     assert(my_down_z_cell_level <= get_cell_coord_z_from_index(particles_coord[idxPart]));
                     assert(get_cell_coord_z_from_index(particles_coord[idxPart]) <= my_top_z_cell_level);
+                    if(inout_index_particles[idxPart] == 3 || inout_index_particles[idxPart] == 52){// TODO
+                        printf("Coord index %ld (tree index %ld)\n", inout_index_particles[idxPart],particles_coord[idxPart]);
+                        printf(">> Box index %ld - %ld - %ld\n", get_cell_coord_x_from_index(particles_coord[idxPart]),
+                               get_cell_coord_y_from_index(particles_coord[idxPart]),
+                               get_cell_coord_z_from_index(particles_coord[idxPart]));
+                    }
                 }
             }
 
@@ -540,7 +546,8 @@ public:
                             }
 
                             const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
-                            const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors);
+                            long int neighbors_indexes[27];
+                            const int nbNeighbors = my_tree.getNeighbors(current_cell_idx, neighbors, neighbors_indexes, true);
 
                             // with other interval
                             for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
@@ -625,6 +632,24 @@ public:
                                                 &particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs],
                                                 dist_r2);
                         }
+                        if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
+                                && inout_index_particles[(intervals[idx_1].first+idx_p2)] == 52)
+                                || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
+                                    && inout_index_particles[(intervals[idx_1].first+idx_p2)] == 3)){// TODO
+                            printf("test interaction between :\n");
+                            printf("index %ld (%ld) pos %e %e %e\n",
+                                   (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
+                                    particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
+                                        particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
+                                    particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
+                            printf("index %ld (%ld) pos %e %e %e\n",
+                                   (intervals[idx_1].first+idx_p2),
+                                   inout_index_particles[(intervals[idx_1].first+idx_p2)],
+                                    particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_X],
+                                    particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Y],
+                                    particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Z]);
+                            printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
+                        }
                     }
                 }
 
@@ -646,6 +671,25 @@ public:
                                                     &particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs],
                                                     dist_r2);
                             }
+
+                            if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
+                                    && inout_index_particles[(intervals[idx_2].first+idx_p2)] == 52)
+                                    || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
+                                        && inout_index_particles[(intervals[idx_2].first+idx_p2)] == 3)){// TODO
+                                printf("interaction between :\n");
+                                printf("index %ld (%ld) pos %e %e %e\n",
+                                       (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
+                                        particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
+                                            particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
+                                        particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
+                                printf("index %ld (%ld) pos %e %e %e\n",
+                                       (intervals[idx_2].first+idx_p2),
+                                       inout_index_particles[(intervals[idx_2].first+idx_p2)],
+                                        particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
+                                        particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
+                                        particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
+                                printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
+                            }
                         }
                     }
                 }
@@ -654,27 +698,52 @@ public:
 
             const long int currenct_cell_idx = iter_cell.first;
             const std::vector<std::pair<partsize_t,partsize_t>>* neighbors[27];
-            const int nbNeighbors = my_tree.getNeighbors(currenct_cell_idx, neighbors);
+            long int neighbors_indexes[27];
+            const int nbNeighbors = my_tree.getNeighbors(currenct_cell_idx, neighbors, neighbors_indexes, false);
+
+            if(iter_cell.first == 3669){ // TODO
+                printf("Box %ld has %d neighbors\n", iter_cell.first, nbNeighbors);
+            }
 
             for(size_t idx_1 = 0 ; idx_1 < intervals.size() ; ++idx_1){
                 // with other interval
                 for(size_t idx_neighbor = 0 ; idx_neighbor < nbNeighbors ; ++idx_neighbor){
-                    for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
-                        for(partsize_t idx_p1 = 0 ; idx_p1 < intervals[idx_1].second ; ++idx_p1){
-                            for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
-                                const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
-                                                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
-                                                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z],
-                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
-                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
-                                                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
-                                if(dist_r2 < cutoff_radius*cutoff_radius){
-                                    in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
-                                                        &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
-                                                        &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
-                                                        &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
-                                                        &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-                                                        dist_r2);
+                    if(currenct_cell_idx < neighbors_indexes[idx_neighbor]){
+                        for(size_t idx_2 = 0 ; idx_2 < (*neighbors[idx_neighbor]).size() ; ++idx_2){
+                            for(partsize_t idx_p1 = 0 ; idx_p1 < intervals[idx_1].second ; ++idx_p1){
+                                for(partsize_t idx_p2 = 0 ; idx_p2 < (*neighbors[idx_neighbor])[idx_2].second ; ++idx_p2){
+                                    const real_number dist_r2 = compute_distance_r2(particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
+                                                                                    particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
+                                                                                    particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z],
+                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
+                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
+                                                                                    particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
+                                    if(dist_r2 < cutoff_radius*cutoff_radius){
+                                        in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
+                                                            &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
+                                                            &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
+                                                            &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
+                                                            &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
+                                                            dist_r2);
+                                    }
+                                    if((inout_index_particles[(intervals[idx_1].first+idx_p1)] == 3
+                                            && inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)] == 52)
+                                            || (inout_index_particles[(intervals[idx_1].first+idx_p1)] == 52
+                                                && inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)] == 3)){// TODO
+                                        printf("interaction between :\n");
+                                        printf("index %ld (%ld) pos %e %e %e\n",
+                                               (intervals[idx_1].first+idx_p1), inout_index_particles[(intervals[idx_1].first+idx_p1)],
+                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_X],
+                                                    particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Y],
+                                                particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions + IDX_Z]);
+                                        printf("index %ld (%ld) pos %e %e %e\n",
+                                               ((*neighbors[idx_neighbor])[idx_2].first+idx_p2),
+                                               inout_index_particles[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)],
+                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_X],
+                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Y],
+                                                particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z]);
+                                        printf("Radius = %e (%e)\n", sqrt(dist_r2), dist_r2);
+                                    }
                                 }
                             }
                         }

bfps/cpp/particles/p2p_tree.hpp

@@ -46,14 +46,16 @@ public:
         return emptyCell;
     }
 
-    int getNeighbors(const long int idx, const CellClass* output[27]) const{
+    int getNeighbors(const long int idx, const CellClass* output[27], long int output_indexes[27], const bool include_target) const{
         int nbNeighbors = 0;
 
+        std::fill_n(output, 27, nullptr);
+
         const long int idx_x = get_cell_coord_x_from_index(idx);
         const long int idx_y = get_cell_coord_y_from_index(idx);
         const long int idx_z = get_cell_coord_z_from_index(idx);
 
-        for(long int neigh_x = - 1 ; neigh_x <= 1 ; ++neigh_x){
+        for(long int neigh_x = -1 ; neigh_x <= 1 ; ++neigh_x){
             long int neigh_x_pbc = neigh_x+idx_x;
             if(neigh_x_pbc < 0){
                 neigh_x_pbc += nb_cell_levels[IDX_X];
@@ -62,7 +64,7 @@ public:
                 neigh_x_pbc -= nb_cell_levels[IDX_X];
             }
 
-            for(long int neigh_y = - 1 ; neigh_y <= 1 ; ++neigh_y){
+            for(long int neigh_y = -1 ; neigh_y <= 1 ; ++neigh_y){
                 long int neigh_y_pbc = neigh_y+idx_y;
                 if(neigh_y_pbc < 0){
                     neigh_y_pbc += nb_cell_levels[IDX_Y];
@@ -71,7 +73,7 @@ public:
                     neigh_y_pbc -= nb_cell_levels[IDX_Y];
                 }
 
-                for(long int neigh_z = - 1 ; neigh_z <= 1 ; ++neigh_z){
+                for(long int neigh_z = -1 ; neigh_z <= 1 ; ++neigh_z){
                     long int neigh_z_pbc = neigh_z+idx_z;
                     if(neigh_z_pbc < 0){
                         neigh_z_pbc += nb_cell_levels[IDX_Z];
@@ -81,13 +83,15 @@ public:
                     }
 
                     // Not the current cell
-                    if(neigh_x_pbc != idx_x || neigh_y_pbc != idx_y || neigh_z_pbc != idx_z){
+                    if(include_target || neigh_x_pbc != idx_x || neigh_y_pbc != idx_y || neigh_z_pbc != idx_z){
                         const long int idx_neigh = get_cell_idx(neigh_x_pbc,
                                                                   neigh_y_pbc,
                                                                   neigh_z_pbc);
                         const auto& iter = data.find(idx_neigh);
                         if(iter != data.end()){
                             output[nbNeighbors] = &(iter->second);
+                            output_indexes[nbNeighbors] = idx_neigh;
+                            nbNeighbors += 1;
                         }
                     }
                 }