diff --git a/include/samurai/algorithm.hpp b/include/samurai/algorithm.hpp
index d07483d53..bd6c21617 100644
--- a/include/samurai/algorithm.hpp
+++ b/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);
diff --git a/include/samurai/field.hpp b/include/samurai/field.hpp
index be02e9162..b4b3de870 100644
--- a/include/samurai/field.hpp
+++ b/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>
diff --git a/include/samurai/mesh.hpp b/include/samurai/mesh.hpp
index cf8f74e24..d53f60feb 100644
--- a/include/samurai/mesh.hpp
+++ b/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()
     {
diff --git a/include/samurai/mr/mesh.hpp b/include/samurai/mr/mesh.hpp
index b847cad98..c4160e06c 100644
--- a/include/samurai/mr/mesh.hpp
+++ b/include/samurai/mr/mesh.hpp
@@ -90,6 +90,17 @@ namespace samurai
 
         void update_sub_mesh_impl();
 
+        template <class MeshT>
+        void construct_ghost_cells(MeshT& mesh, cl_type& cell_list);
+        template <class MeshT>
+        void construct_mra_cells(MeshT& mesh, cl_type& cell_list);
+        template <class MeshT>
+        void construct_periodic_cells(MeshT& mesh, cl_type& cell_list);
+        template <class MeshT>
+        void construct_projection_cells(MeshT& mesh, cl_type& cell_list);
+        template <class MeshT>
+        void construct_neighbour_cells(MeshT& mesh, cl_type& cell_list);
+
         template <typename... T>
         xt::xtensor<bool, 1> exists(mesh_id_t type, std::size_t level, interval_t interval, T... index) const;
     };
@@ -140,38 +151,11 @@ namespace samurai
     }
 
     template <class Config>
-    inline void MRMesh<Config>::update_sub_mesh_impl()
+    template <class MeshT>
+    void MRMesh<Config>::construct_ghost_cells(MeshT& mesh, cl_type& cell_list)
     {
-#ifdef SAMURAI_WITH_MPI
-        mpi::communicator world;
-        // cppcheck-suppress redundantInitialization
-        auto max_level = mpi::all_reduce(world, this->cells()[mesh_id_t::cells].max_level(), mpi::maximum<std::size_t>());
-        // cppcheck-suppress redundantInitialization
-        auto min_level = mpi::all_reduce(world, this->cells()[mesh_id_t::cells].min_level(), mpi::minimum<std::size_t>());
-        cl_type cell_list;
-#else
-        // cppcheck-suppress redundantInitialization
-        auto max_level = this->cells()[mesh_id_t::cells].max_level();
-        // cppcheck-suppress redundantInitialization
-        auto min_level = this->cells()[mesh_id_t::cells].min_level();
-        cl_type cell_list;
-#endif
-        // Construction of ghost cells
-        // ===========================
-        //
-        // Example with max_stencil_width = 1
-        //
-        // level 2                       |.|-|-|-|-|.|                   |-|
-        // cells
-        //                                                               |.|
-        //                                                               ghost
-        //                                                               cells
-        // level 1             |...|---|---|...|...|---|---|...|
-        //
-        // level 0 |.......|-------|.......|       |.......|-------|.......|
-        //
         for_each_interval(
-            this->cells()[mesh_id_t::cells],
+            mesh[mesh_id_t::cells],
             [&](std::size_t level, const auto& interval, const auto& index_yz)
             {
                 lcl_type& lcl = cell_list[level];
@@ -182,168 +166,279 @@ namespace samurai
                         lcl[index].add_interval({interval.start - config::max_stencil_width, interval.end + config::max_stencil_width});
                     });
             });
-        this->cells()[mesh_id_t::cells_and_ghosts] = {cell_list, false};
+        mesh[mesh_id_t::cells_and_ghosts] = {cell_list, false};
+    }
 
-        // Add cells for the MRA
-        if (this->max_level() != this->min_level())
+    template <class Config>
+    template <class MeshT>
+    void MRMesh<Config>::construct_mra_cells(MeshT& mesh, cl_type& cell_list)
+    {
+        // cppcheck-suppress redundantInitialization
+        auto max_level = mesh.cells()[mesh_id_t::cells].max_level();
+        // cppcheck-suppress redundantInitialization
+        auto min_level = mesh.cells()[mesh_id_t::cells].min_level();
+
+        for (std::size_t level = max_level; level >= ((min_level == 0) ? 1 : min_level); --level)
         {
-            for (std::size_t level = max_level; level >= ((min_level == 0) ? 1 : min_level); --level)
-            {
-                auto expr = difference(this->cells()[mesh_id_t::cells_and_ghosts][level], this->get_union()[level]).on(level);
+            auto expr = difference(mesh.cells()[mesh_id_t::cells_and_ghosts][level], mesh.get_union()[level]).on(level);
+            expr(
+                [&](const auto& interval, const auto& index_yz)
+                {
+                    lcl_type& lcl = cell_list[level - 1];
+                    static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                        [&](auto stencil)
+                        {
+                            auto new_interval = interval >> 1;
+                            lcl[(index_yz >> 1) + stencil].add_interval(
+                                {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
+                        });
+                });
 
-                expr(
-                    [&](const auto& interval, const auto& index_yz)
+            auto expr_2 = intersection(mesh.cells()[mesh_id_t::cells][level], mesh.cells()[mesh_id_t::cells][level]);
+
+            expr_2(
+                [&](const auto& interval, const auto& index_yz)
+                {
+                    if (level - 1 > 0)
                     {
-                        lcl_type& lcl = cell_list[level - 1];
+                        lcl_type& lcl = cell_list[level - 2];
 
                         static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
                             [&](auto stencil)
                             {
-                                auto new_interval = interval >> 1;
-                                lcl[(index_yz >> 1) + stencil].add_interval(
+                                auto new_interval = interval >> 2;
+                                lcl[(index_yz >> 2) + stencil].add_interval(
                                     {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
                             });
-                    });
+                    }
+                });
+        }
+        mesh[mesh_id_t::all_cells] = {cell_list, false};
+    }
 
-                auto expr_2 = intersection(this->cells()[mesh_id_t::cells][level], this->cells()[mesh_id_t::cells][level]);
+    template <class Config>
+    template <class MeshT>
+    void MRMesh<Config>::construct_periodic_cells(MeshT& mesh, cl_type& cell_list)
+    {
+        // cppcheck-suppress redundantInitialization
+        auto max_level = mesh.cells()[mesh_id_t::cells].max_level();
+        // cppcheck-suppress redundantInitialization
+        auto min_level = mesh.cells()[mesh_id_t::cells].min_level();
 
-                expr_2(
-                    [&](const auto& interval, const auto& index_yz)
-                    {
-                        if (level - 1 > 0)
-                        {
-                            lcl_type& lcl = cell_list[level - 2];
+        xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> stencil;
+        xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> min_corner;
+        xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> max_corner;
 
-                            static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
-                                [&](auto stencil)
-                                {
-                                    auto new_interval = interval >> 2;
-                                    lcl[(index_yz >> 2) + stencil].add_interval(
-                                        {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
-                                });
-                        }
-                    });
-            }
-            this->cells()[mesh_id_t::all_cells] = {cell_list, false};
+        auto& domain     = mesh.domain();
+        auto min_indices = domain.min_indices();
+        auto max_indices = domain.max_indices();
 
-            this->update_mesh_neighbour();
+        // cppcheck-suppress constStatement
+        for (std::size_t level = min_level; level <= max_level; ++level)
+        {
+            std::size_t delta_l = domain.level() - level;
+            lcl_type& lcl       = cell_list[level];
 
-            for (auto& neighbour : this->mpi_neighbourhood())
+            for (std::size_t d = 0; d < dim; ++d)
             {
-                for (std::size_t level = 0; level <= max_level; ++level)
+                if (mesh.is_periodic(d))
                 {
-                    auto expr = intersection(this->subdomain(), neighbour.mesh[mesh_id_t::reference][level]).on(level);
-                    expr(
-                        [&](const auto& interval, const auto& index_yz)
+                    stencil.fill(0);
+                    stencil[d] = (max_indices[d] - min_indices[d]) >> delta_l;
+
+                    min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
+                    max_corner[d] = (min_indices[d] >> delta_l) + config::ghost_width;
+                    for (std::size_t dd = 0; dd < dim; ++dd)
+                    {
+                        if (dd != d)
                         {
-                            lcl_type& lcl = cell_list[level];
-                            lcl[index_yz].add_interval(interval);
+                            min_corner[dd] = (min_indices[dd] >> delta_l) - config::ghost_width;
+                            max_corner[dd] = (max_indices[dd] >> delta_l) + config::ghost_width;
+                        }
+                    }
 
-                            if (level > neighbour.mesh[mesh_id_t::reference].min_level())
-                            {
-                                lcl_type& lclm1 = cell_list[level - 1];
-
-                                static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
-                                    [&](auto stencil)
-                                    {
-                                        auto new_interval = interval >> 1;
-                                        lclm1[(index_yz >> 1) + stencil].add_interval(
-                                            {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
-                                    });
-                            }
+                    lca_type lca1{
+                        level,
+                        Box<typename interval_t::value_t, dim>{min_corner, max_corner}
+                    };
+
+                    auto set1 = intersection(mesh.cells()[mesh_id_t::reference][level], lca1);
+                    set1(
+                        [&](const auto& i, const auto& index_yz)
+                        {
+                            lcl[index_yz + xt::view(stencil, xt::range(1, _))].add_interval(i + stencil[0]);
                         });
+
+                    min_corner[d] = (max_indices[d] >> delta_l) - config::ghost_width;
+                    max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
+                    lca_type lca2{
+                        level,
+                        Box<typename interval_t::value_t, dim>{min_corner, max_corner}
+                    };
+
+                    auto set2 = intersection(mesh.cells()[mesh_id_t::reference][level], lca2);
+                    set2(
+                        [&](const auto& i, const auto& index_yz)
+                        {
+                            lcl[index_yz - xt::view(stencil, xt::range(1, _))].add_interval(i - stencil[0]);
+                        });
+                    mesh.cells()[mesh_id_t::all_cells][level] = {lcl};
                 }
             }
-            this->cells()[mesh_id_t::all_cells] = {cell_list, false};
-
-            // add ghosts for periodicity
-            xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> stencil;
-            xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> min_corner;
-            xt::xtensor_fixed<typename interval_t::value_t, xt::xshape<dim>> max_corner;
+        }
+    }
 
-            auto& domain     = this->domain();
-            auto min_indices = domain.min_indices();
-            auto max_indices = domain.max_indices();
+    template <class Config>
+    template <class MeshT>
+    void MRMesh<Config>::construct_projection_cells(MeshT& mesh, cl_type& cell_list)
+    {
+        // cppcheck-suppress redundantInitialization
+        auto max_level = mesh.cells()[mesh_id_t::cells].max_level();
 
-            for (std::size_t level = this->cells()[mesh_id_t::reference].min_level();
-                 level <= this->cells()[mesh_id_t::reference].max_level();
-                 ++level)
-            {
-                std::size_t delta_l = domain.level() - level;
-                lcl_type& lcl       = cell_list[level];
+        for (std::size_t level = 0; level < max_level; ++level)
+        {
+            lcl_type& lcl = cell_list[level + 1];
+            lcl_type lcl_proj{level};
+            auto expr = intersection(mesh.cells()[mesh_id_t::all_cells][level], mesh.get_union()[level]);
 
-                for (std::size_t d = 0; d < dim; ++d)
+            expr(
+                [&](const auto& interval, const auto& index_yz)
                 {
-                    if (this->is_periodic(d))
-                    {
-                        stencil.fill(0);
-                        stencil[d] = (max_indices[d] - min_indices[d]) >> delta_l;
-
-                        min_corner[d] = (min_indices[d] >> delta_l) - config::ghost_width;
-                        max_corner[d] = (min_indices[d] >> delta_l) + config::ghost_width;
-                        for (std::size_t dd = 0; dd < dim; ++dd)
+                    static_nested_loop<dim - 1, 0, 2>(
+                        [&](auto s)
                         {
-                            if (dd != d)
-                            {
-                                min_corner[dd] = (min_indices[dd] >> delta_l) - config::ghost_width;
-                                max_corner[dd] = (max_indices[dd] >> delta_l) + config::ghost_width;
-                            }
-                        }
-
-                        lca_type lca1{
-                            level,
-                            Box<typename interval_t::value_t, dim>{min_corner, max_corner}
-                        };
-
-                        auto set1 = intersection(this->cells()[mesh_id_t::reference][level], lca1);
-                        set1(
-                            [&](const auto& i, const auto& index_yz)
-                            {
-                                lcl[index_yz + xt::view(stencil, xt::range(1, _))].add_interval(i + stencil[0]);
-                            });
-
-                        min_corner[d] = (max_indices[d] >> delta_l) - config::ghost_width;
-                        max_corner[d] = (max_indices[d] >> delta_l) + config::ghost_width;
-                        lca_type lca2{
-                            level,
-                            Box<typename interval_t::value_t, dim>{min_corner, max_corner}
-                        };
+                            lcl[(index_yz << 1) + s].add_interval(interval << 1);
+                        });
+                    lcl_proj[index_yz].add_interval(interval);
+                });
+            mesh.cells()[mesh_id_t::all_cells][level + 1] = lcl;
+            mesh.cells()[mesh_id_t::proj_cells][level]    = lcl_proj;
+        }
+    }
 
-                        auto set2 = intersection(this->cells()[mesh_id_t::reference][level], lca2);
-                        set2(
-                            [&](const auto& i, const auto& index_yz)
-                            {
-                                lcl[index_yz - xt::view(stencil, xt::range(1, _))].add_interval(i - stencil[0]);
-                            });
-                        this->cells()[mesh_id_t::all_cells][level] = {lcl};
-                    }
-                }
-            }
+    template <class Config>
+    template <class MeshT>
+    void MRMesh<Config>::construct_neighbour_cells(MeshT& mesh, cl_type& cell_list)
+    {
+        // cppcheck-suppress redundantInitialization
+        auto max_level = mesh.cells()[mesh_id_t::cells].max_level();
 
-            for (std::size_t level = 0; level < max_level; ++level)
+        for (auto& neighbour : mesh.mpi_neighbourhood())
+        {
+            for (std::size_t level = 0; level <= max_level; ++level)
             {
-                lcl_type& lcl = cell_list[level + 1];
-                lcl_type lcl_proj{level};
-                auto expr = intersection(this->cells()[mesh_id_t::all_cells][level], this->get_union()[level]);
-
+                auto expr = intersection(mesh.subdomain(), neighbour.mesh[mesh_id_t::reference][level]).on(level);
                 expr(
                     [&](const auto& interval, const auto& index_yz)
                     {
-                        static_nested_loop<dim - 1, 0, 2>(
-                            [&](auto s)
-                            {
-                                lcl[(index_yz << 1) + s].add_interval(interval << 1);
-                            });
-                        lcl_proj[index_yz].add_interval(interval);
+                        lcl_type& lcl = cell_list[level];
+                        lcl[index_yz].add_interval(interval);
+
+                        if (level > neighbour.mesh[mesh_id_t::reference].min_level())
+                        {
+                            lcl_type& lclm1 = cell_list[level - 1];
+
+                            static_nested_loop<dim - 1, -config::prediction_order, config::prediction_order + 1>(
+                                [&](auto stencil)
+                                {
+                                    auto new_interval = interval >> 1;
+                                    lclm1[(index_yz >> 1) + stencil].add_interval(
+                                        {new_interval.start - config::prediction_order, new_interval.end + config::prediction_order});
+                                });
+                        }
                     });
-                this->cells()[mesh_id_t::all_cells][level + 1] = lcl;
-                this->cells()[mesh_id_t::proj_cells][level]    = lcl_proj;
             }
-            this->update_mesh_neighbour();
+        }
+        mesh.cells()[mesh_id_t::all_cells] = {cell_list, false};
+    }
+
+    template <class Config>
+    inline void MRMesh<Config>::update_sub_mesh_impl()
+    {
+        cl_type cell_list;
+#ifdef SAMURAI_WITH_MPI
+        std::vector<cl_type> neighbour_cell_list(this->mpi_neighbourhood().size());
+#endif
+        // Construction of ghost cells
+        // ===========================
+        //
+        // Example with max_stencil_width = 1
+        //
+        // level 2                       |.|-|-|-|-|.|                   |-|
+        // cells
+        //                                                               |.|
+        //                                                               ghost
+        //                                                               cells
+        // level 1             |...|---|---|...|...|---|---|...|
+        //
+        // level 0 |.......|-------|.......|       |.......|-------|.......|
+        //
+        //
+        //
+
+        construct_ghost_cells(*this, cell_list);
+        // reconstruct ghost cells for neighbourhood
+#ifdef SAMURAI_WITH_MPI
+        for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+        {
+            auto& neighbour = this->mpi_neighbourhood()[i];
+            construct_ghost_cells((neighbour.mesh), neighbour_cell_list[i]);
+        }
+#endif
+        // Add cells for the MRA
+        if (this->max_level() != this->min_level())
+        {
+            construct_mra_cells(*this, cell_list);
+            // this->update_neighbour_reference();
+#ifdef SAMURAI_WITH_MPI
+            for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+            {
+                auto& neighbour = this->mpi_neighbourhood()[i];
+                construct_mra_cells((neighbour.mesh), neighbour_cell_list[i]);
+            }
+#endif
+            construct_neighbour_cells(*this, cell_list);
+#ifdef SAMURAI_WITH_MPI
+            for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+            {
+                auto& neighbour = this->mpi_neighbourhood()[i];
+                construct_neighbour_cells(neighbour.mesh, neighbour_cell_list[i]);
+            }
+#endif
+            this->cells()[mesh_id_t::all_cells] = {cell_list, false};
+#ifdef SAMURAI_WITH_MPI
+            for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+            {
+                auto& neighbour                              = this->mpi_neighbourhood()[i];
+                neighbour.mesh.cells()[mesh_id_t::all_cells] = {neighbour_cell_list[i], false};
+            }
+#endif
+            // For now, this leads to a segfault.
+            construct_periodic_cells(*this, cell_list);
+#ifdef SAMURAI_WITH_MPI
+            //          for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+            //          {
+            //              auto& neighbour = this->mpi_neighbourhood()[i];
+            //              construct_periodic_cells(neighbour.mesh, neighbour_cell_list[i]);
+            //          }
+#endif
+            construct_projection_cells(*this, cell_list);
+#ifdef SAMURAI_WITH_MPI
+            for (std::size_t i = 0; i < this->mpi_neighbourhood().size(); i++)
+            {
+                auto& neighbour = this->mpi_neighbourhood()[i];
+                construct_projection_cells(neighbour.mesh, neighbour_cell_list[i]);
+            }
+#endif
+            //
+            // this->update_mesh_neighbour();
+            this->update_neighbour_subdomain();
+            //            this->update_neighbour_all_cells();
         }
         else
         {
             this->cells()[mesh_id_t::all_cells] = {cell_list, false};
+            // TODO : what to send here ?
             this->update_mesh_neighbour();
         }
     }