feat: [WIP] Find neighbourhood in MPI

include/samurai/mesh.hpp

@@ -170,6 +170,7 @@ namespace samurai
         void construct_union();
         void update_sub_mesh();
         void renumbering();
+        void find_neighbour();
         void partition_mesh(std::size_t start_level, const Box<double, dim>& global_box);
         void load_balancing();
         void load_transfer(const std::vector<double>& load_fluxes);
@@ -714,6 +715,137 @@ namespace samurai
         }
     }
 
+    template <std::size_t Dim>
+    struct StencilProvider
+    {
+        static_assert(Dim == 2 || Dim == 3, "Dimension non supportée pour les stencils");
+    };
+
+    // Spécialisation pour dim = 2
+    template <>
+    struct StencilProvider<2>
+    {
+        static constexpr std::size_t stencil_size = 8;
+        inline static const xt::xtensor_fixed<int, xt::xshape<stencil_size, 2>> stencils{
+            {-1, 0 },
+            {1,  0 },
+            {0,  -1},
+            {0,  1 },
+            {-1, 1 },
+            {1,  1 },
+            {-1, -1},
+            {1,  -1}
+        };
+    };
+
+    // Spécialisation pour dim = 3
+    template <>
+    struct StencilProvider<3>
+    {
+        static constexpr std::size_t stencil_size = 27;
+        inline static const xt::xtensor_fixed<int, xt::xshape<stencil_size, 3>> stencils{
+            {-1, -1, -1},
+            {-1, -1, 0 },
+            {-1, -1, 1 },
+            {-1, 0,  -1},
+            {-1, 0,  0 },
+            {-1, 0,  1 },
+            {-1, 1,  -1},
+            {-1, 1,  0 },
+            {-1, 1,  1 },
+            {0,  -1, -1},
+            {0,  -1, 0 },
+            {0,  -1, 1 },
+            {0,  0,  -1},
+            {0,  0,  1 },
+            {0,  1,  -1},
+            {0,  1,  0 },
+            {0,  1,  1 },
+            {1,  -1, -1},
+            {1,  -1, 0 },
+            {1,  -1, 1 },
+            {1,  0,  -1},
+            {1,  0,  0 },
+            {1,  0,  1 },
+            {1,  1,  -1},
+            {1,  1,  0 },
+            {1,  1,  1 }
+        };
+    };
+
+    template <class D, class Config>
+    void Mesh_base<D, Config>::find_neighbour()
+    {
+        mpi::communicator world;
+        auto rank = world.rank();
+        auto size = world.size();
+
+        m_mpi_neighbourhood.reserve(static_cast<std::size_t>(size) - 1);
+        // This code find the neighbourhood of each mpi rank
+        // 1) we have to exchange the local mesh with all the neighbourhood. It is costly and not scalable but is it a easy way to do it
+        // 2) we find intersection between each subdomain cells_and_ghosts.
+        // 3) be careful : for periodic conditions, how to treat neighbourhood ?? maybe in another mpi neighbourhood.
+
+        for (int ir = 0; ir < size; ++ir)
+        {
+            if (ir != rank)
+            {
+                m_mpi_neighbourhood.push_back(ir);
+                ;
+            }
+        }
+
+        // send
+
+        // Please think about if they are all necessary
+        construct_subdomain();
+        construct_union();
+        update_sub_mesh();
+        renumbering();
+        update_mesh_neighbour();
+
+        std::vector<mpi_subdomain_t> m_mpi_neighbourhood_temp;
+
+        const auto& stencils = StencilProvider<dim>::stencils;
+        for (auto& neighbour : m_mpi_neighbourhood)
+        {
+            int est_voisin_global = 0;
+            {
+                // I use cells and ghosts becose cell only is not enough
+                // in contrary, cells and ghost are overkill
+                // maybe we can do cells for mesh 1 and cells_and_ghosts for mesh 2
+
+                // We need a kind of ""expand"" operator. but for now we use bruteforce translation
+
+                for (std::size_t is = 0; is < stencils.shape()[0]; ++is)
+                {
+                    auto s           = xt::view(stencils, is);
+                    auto translation = translate(this->m_subdomain, s);
+                    auto overlap     = intersection(translation, neighbour.mesh.m_subdomain);
+                    // tester si voisin
+                    int est_voisin_niveau = 0;
+                    // maybe it exist a method like "overlap.is_empty() so we should use it"
+                    overlap(
+                        [&](const auto& i, const auto& index)
+                        {
+                            est_voisin_niveau = 1;
+                        });
+                    if (est_voisin_niveau == 1)
+                    {
+                        est_voisin_global = 1;
+                        break;
+                    }
+                }
+            }
+            if (est_voisin_global)
+            {
+                m_mpi_neighbourhood_temp.push_back(neighbour);
+            }
+        }
+
+        std::swap(m_mpi_neighbourhood_temp, m_mpi_neighbourhood);
+    }
+
     template <class D, class Config>
     void Mesh_base<D, Config>::partition_mesh([[maybe_unused]] std::size_t start_level, [[maybe_unused]] const Box<double, dim>& global_box)
     {
@@ -787,7 +919,7 @@ namespace samurai
         std::size_t subdomain_end   = 0;
         lcl_type subdomain_cells(start_level, m_domain.origin_point(), m_domain.scaling_factor());
         // in 1D MPI, we need a specific partitioning
-        if (dim == 1)
+        if constexpr (dim == 1)
         {
             std::size_t n_cells               = m_domain.nb_cells();
             std::size_t n_cells_per_subdomain = n_cells / static_cast<std::size_t>(size);
@@ -810,7 +942,7 @@ namespace samurai
                                       }
                                   });
         }
-        else if (dim >= 2)
+        else if constexpr (dim >= 2)
         {
             auto subdomain_nb_intervals = m_domain.nb_intervals() / static_cast<std::size_t>(size);
             subdomain_start             = static_cast<std::size_t>(rank) * subdomain_nb_intervals;
@@ -833,14 +965,8 @@ namespace samurai
 
         this->m_cells[mesh_id_t::cells][start_level] = subdomain_cells;
 
-        m_mpi_neighbourhood.reserve(static_cast<std::size_t>(size) - 1);
-        for (int ir = 0; ir < size; ++ir)
-        {
-            if (ir != rank)
-            {
-                m_mpi_neighbourhood.push_back(ir);
-            }
-        }
+        // in theory, at this step, we can construct the neighbourhood by construction.
+        find_neighbour();
 
         // // Neighbours
         // m_mpi_neighbourhood.reserve(static_cast<std::size_t>(pow(3, dim) - 1));