diff --git a/bfps/cpp/particles/p2p_distr_mpi.hpp b/bfps/cpp/particles/p2p_distr_mpi.hpp
index 7f610be18c749077ef73730ebec365b21e212569..a32c33756ccd88df08d5834bfcc33f7bbf7bbf40 100644
--- a/bfps/cpp/particles/p2p_distr_mpi.hpp
+++ b/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);
+ }
}
}
}
diff --git a/bfps/cpp/particles/p2p_tree.hpp b/bfps/cpp/particles/p2p_tree.hpp
index b11cd826f74e59b8014d4563098069873cf7ded8..4715dfad7fecc1c564a43972504ddd1e621493b2 100644
--- a/bfps/cpp/particles/p2p_tree.hpp
+++ b/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;
}
}
}