feat: [WIP] Reconstruct neighbour mesh locally 2

include/samurai/algorithm.hpp

@@ -59,7 +59,7 @@ namespace samurai
     }
 
     template <class Mesh, class Func>
-    inline void for_each_level(Mesh& mesh, Func&& f, bool include_empty_levels = false)
+    inline void for_each_level(const Mesh& mesh, Func&& f, bool include_empty_levels = false)
     {
         using mesh_id_t = typename Mesh::mesh_id_t;
         for_each_level(mesh[mesh_id_t::cells], std::forward<Func>(f), include_empty_levels);

include/samurai/field.hpp

@@ -1188,14 +1188,15 @@ namespace samurai
     inline auto ScalarField<mesh_t, value_t>::begin() -> iterator
     {
         using mesh_id_t = typename mesh_t::mesh_id_t;
-        return iterator(this, this->mesh()[mesh_id_t::cells].begin());
+        // add comment : the CI is buggy ahahaha
+        return iterator(this, std::as_const(this->mesh()[mesh_id_t::cells]).begin());
     }
 
     template <class mesh_t, class value_t>
     inline auto ScalarField<mesh_t, value_t>::end() -> iterator
     {
         using mesh_id_t = typename mesh_t::mesh_id_t;
-        return iterator(this, this->mesh()[mesh_id_t::cells].end());
+        return iterator(this, std::as_const(this->mesh()[mesh_id_t::cells]).end());
     }
 
     template <class mesh_t, class value_t>

include/samurai/mesh.hpp

@@ -89,6 +89,7 @@ namespace samurai
         std::size_t nb_cells(std::size_t level, mesh_id_t mesh_id = mesh_id_t::reference) const;
 
         const ca_type& operator[](mesh_id_t mesh_id) const;
+        ca_type& operator[](mesh_id_t mesh_id);
 
         std::size_t max_level() const;
         std::size_t& max_level();
@@ -133,6 +134,16 @@ namespace samurai
         cell_t get_cell(std::size_t level, const xt::xexpression<E>& coord) const;
 
         void update_mesh_neighbour();
+
+        void update_neighbour_subdomain();
+        template <mesh_id_t MeshID>
+        void update_meshid_neighbour();
+
+        void update_neighbour_cells();
+        void update_neighbour_cells_and_ghosts();
+        void update_neighbour_all_cells();
+        void update_neighbour_reference();
+
         void to_stream(std::ostream& os) const;
 
       protected:
@@ -403,6 +414,12 @@ namespace samurai
         return m_cells[mesh_id];
     }
 
+    template <class D, class Config>
+    inline auto Mesh_base<D, Config>::operator[](mesh_id_t mesh_id) -> ca_type&
+    {
+        return m_cells[mesh_id];
+    }
+
     template <class D, class Config>
     inline std::size_t Mesh_base<D, Config>::max_level() const
     {
@@ -644,6 +661,94 @@ namespace samurai
 #endif
     }
 
+    // TODO : find a clever way to factorize the two next functions. For new, I have to duplicate the code 2 times.
+
+    // This function is to only send m_subdomain instead of the whole mesh data
+    template <class D, class Config>
+    inline void Mesh_base<D, Config>::update_neighbour_subdomain()
+    {
+#ifdef SAMURAI_WITH_MPI
+        // send/recv the meshes of the neighbouring subdomains
+        mpi::communicator world;
+        std::vector<mpi::request> req;
+
+        boost::mpi::packed_oarchive::buffer_type buffer;
+        boost::mpi::packed_oarchive oa(world, buffer);
+        oa << derived_cast().m_subdomain;
+
+        std::transform(m_mpi_neighbourhood.cbegin(),
+                       m_mpi_neighbourhood.cend(),
+                       std::back_inserter(req),
+                       [&](const auto& neighbour)
+                       {
+                           return world.isend(neighbour.rank, neighbour.rank, buffer);
+                       });
+
+        for (auto& neighbour : m_mpi_neighbourhood)
+        {
+            world.recv(neighbour.rank, world.rank(), neighbour.mesh.m_subdomain);
+        }
+
+        mpi::wait_all(req.begin(), req.end());
+#endif
+    }
+
+    // This function is to only send cells[i] instead of the whole mesh
+    template <class D, class Config>
+    template <typename Mesh_base<D, Config>::mesh_id_t MeshID>
+    inline void Mesh_base<D, Config>::update_meshid_neighbour()
+    {
+#ifdef SAMURAI_WITH_MPI
+        // send/recv the meshes of the neighbouring subdomains
+        mpi::communicator world;
+        std::vector<mpi::request> req;
+
+        boost::mpi::packed_oarchive::buffer_type buffer;
+        boost::mpi::packed_oarchive oa(world, buffer);
+        oa << derived_cast()[MeshID];
+
+        std::transform(m_mpi_neighbourhood.cbegin(),
+                       m_mpi_neighbourhood.cend(),
+                       std::back_inserter(req),
+                       [&](const auto& neighbour)
+                       {
+                           return world.isend(neighbour.rank, neighbour.rank, buffer);
+                       });
+
+        for (auto& neighbour : m_mpi_neighbourhood)
+        {
+            world.recv(neighbour.rank, world.rank(), neighbour.mesh[MeshID]);
+        }
+
+        mpi::wait_all(req.begin(), req.end());
+#endif
+    }
+
+    // These are helper functions to send a specific cells[i]
+    template <class D, class Config>
+    inline void Mesh_base<D, Config>::update_neighbour_cells()
+    {
+        update_meshid_neighbour<mesh_id_t::cells>();
+    }
+
+    template <class D, class Config>
+    inline void Mesh_base<D, Config>::update_neighbour_cells_and_ghosts()
+    {
+        update_meshid_neighbour<mesh_id_t::cells_and_ghosts>();
+    }
+
+    template <class D, class Config>
+    inline void Mesh_base<D, Config>::update_neighbour_all_cells()
+    {
+        update_meshid_neighbour<mesh_id_t::all_cells>();
+    }
+
+    template <class D, class Config>
+    inline void Mesh_base<D, Config>::update_neighbour_reference()
+    {
+        update_meshid_neighbour<mesh_id_t::reference>();
+    }
+
     template <class D, class Config>
     inline void Mesh_base<D, Config>::construct_subdomain()
     {