diff --git a/python/py_mesh.cc b/python/py_mesh.cc
index f30ba8505..e06154338 100644
--- a/python/py_mesh.cc
+++ b/python/py_mesh.cc
@@ -1,200 +1,207 @@
 /**
  * @file   py_mesh.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Philip Mueller <philip.paul.mueller@bluemail.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Jun 16 2019
  * @date last modification: Mon Mar 15 2021
  *
  * @brief  pybind11 interface to Mesh
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_config.hh"
 /* -------------------------------------------------------------------------- */
 #include "py_aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <mesh.hh>
 #include <mesh_accessor.hh>
 #include <mesh_utils.hh>
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 namespace py = pybind11;
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 namespace {
   /* ------------------------------------------------------------------------ */
   template <typename T>
   void register_element_type_map_array(py::module & mod,
                                        const std::string & name) {
     py::class_<ElementTypeMapArray<T>, std::shared_ptr<ElementTypeMapArray<T>>>(
         mod, ("ElementTypeMapArray" + name).c_str())
         .def(
             "__call__",
             [](ElementTypeMapArray<T> & self, ElementType type,
                GhostType ghost_type) -> decltype(auto) {
               return self(type, ghost_type);
             },
             py::arg("type"), py::arg("ghost_type") = _not_ghost,
             py::return_value_policy::reference)
         .def(
             "elementTypes",
             [](ElementTypeMapArray<T> & self, UInt _dim, GhostType _ghost_type,
                ElementKind _kind) -> std::vector<ElementType> {
               auto types = self.elementTypes(_dim, _ghost_type, _kind);
               std::vector<ElementType> _types;
               for (auto && t : types) {
                 _types.push_back(t);
               }
               return _types;
             },
             py::arg("dim") = _all_dimensions,
             py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular);
   }
 } // namespace
 /* -------------------------------------------------------------------------- */
 void register_mesh(py::module & mod) {
 
   register_element_type_map_array<Real>(mod, "Real");
   register_element_type_map_array<UInt>(mod, "UInt");
   // register_element_type_map_array<std::string>(mod, "String");
 
+  py::class_<Mesh::PeriodicSlaves>(mod, "PeriodicSlaves")
+      .def(
+          "__iter__",
+          [](Mesh::PeriodicSlaves & _this) {
+            py::make_iterator(_this.begin(), _this.end());
+          },
+          py::keep_alive<0, 1>());
+
   py::class_<MeshData>(mod, "MeshData")
       .def(
           "getElementalDataUInt",
           [](MeshData & _this, const ID & name) -> decltype(auto) {
             return _this.getElementalData<UInt>(name);
           },
           py::return_value_policy::reference)
       .def(
           "getElementalDataReal",
           [](MeshData & _this, const ID & name) -> decltype(auto) {
             return _this.getElementalData<Real>(name);
           },
           py::return_value_policy::reference);
 
   py::class_<Mesh, GroupManager, Dumpable, MeshData>(mod, "Mesh",
                                                      py::multiple_inheritance())
       .def(py::init<UInt, const ID &>(), py::arg("spatial_dimension"),
            py::arg("id") = "mesh")
       .def("read", &Mesh::read, py::arg("filename"),
            py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file")
       .def(
           "getNodes",
           [](Mesh & self) -> decltype(auto) { return self.getNodes(); },
           py::return_value_policy::reference)
       .def("getNbNodes", &Mesh::getNbNodes)
       .def(
           "getConnectivity",
           [](Mesh & self, ElementType type) -> decltype(auto) {
             return self.getConnectivity(type);
           },
           py::return_value_policy::reference)
       .def(
           "addConnectivityType",
           [](Mesh & self, ElementType type, GhostType ghost_type) -> void {
             self.addConnectivityType(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("distribute", [](Mesh & self) { self.distribute(); })
       .def("fillNodesToElements", &Mesh::fillNodesToElements,
            py::arg("dimension") = _all_dimensions)
       .def("getAssociatedElements",
            [](Mesh & self, const UInt & node, py::list list) {
              Array<Element> elements;
              self.getAssociatedElements(node, elements);
              for (auto && element : elements) {
                list.append(element);
              }
            })
       .def("makePeriodic",
            [](Mesh & self, const SpatialDirection & direction) {
              self.makePeriodic(direction);
            })
       .def(
           "getNbElement",
           [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type,
              ElementKind kind) {
             return self.getNbElement(spatial_dimension, ghost_type, kind);
           },
           py::arg("spatial_dimension") = _all_dimensions,
           py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined)
       .def(
           "getNbElement",
           [](Mesh & self, ElementType type, GhostType ghost_type) {
             return self.getNbElement(type, ghost_type);
           },
           py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def_static(
           "getSpatialDimension",
           [](ElementType & type) { return Mesh::getSpatialDimension(type); })
       .def(
           "getDataReal",
           [](Mesh & _this, const ID & name, ElementType type,
              GhostType ghost_type) -> decltype(auto) {
             return _this.getData<Real>(name, type, ghost_type);
           },
           py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost,
           py::return_value_policy::reference)
       .def(
           "hasDataReal",
           [](Mesh & _this, const ID & name, ElementType type,
              GhostType ghost_type) -> bool {
             return _this.hasData<Real>(name, type, ghost_type);
           },
           py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost)
       .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); })
       .def("getPeriodicMaster", &Mesh::getPeriodicMaster)
       .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves)
       .def("isPeriodicSlave", &Mesh::isPeriodicSlave)
-      .def("isPeriodicMaster", &Mesh::isPeriodicMaster)
-      .def("getPeriodicMasterSlaves", &Mesh::getPeriodicMasterSlaves);
+      .def("isPeriodicMaster", &Mesh::isPeriodicMaster);
 
   /* ------------------------------------------------------------------------ */
   py::class_<MeshUtils>(mod, "MeshUtils")
       .def_static("buildFacets", &MeshUtils::buildFacets);
 
   py::class_<MeshAccessor>(mod, "MeshAccessor")
       .def(py::init<Mesh &>(), py::arg("mesh"))
       .def(
           "resizeConnectivity",
           [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt)
               -> void { self.resizeConnectivity(new_size, type, gt); },
           py::arg("new_size"), py::arg("type"),
           py::arg("ghost_type") = _not_ghost)
       .def(
           "resizeNodes",
           [](MeshAccessor & self, UInt new_size) -> void {
             self.resizeNodes(new_size);
           },
           py::arg("new_size"))
       .def("makeReady", &MeshAccessor::makeReady);
 }
 } // namespace akantu
diff --git a/src/mesh/mesh_inline_impl.hh b/src/mesh/mesh_inline_impl.hh
index 7b5e3aff2..22d27cfd4 100644
--- a/src/mesh/mesh_inline_impl.hh
+++ b/src/mesh/mesh_inline_impl.hh
@@ -1,764 +1,765 @@
 /**
  * @file   mesh_inline_impl.hh
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Dana Christen <dana.christen@epfl.ch>
  * @author Mohit Pundir <mohit.pundir@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Thu Jul 15 2010
  * @date last modification: Fri Dec 11 2020
  *
  * @brief  Implementation of the inline functions of the mesh class
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_iterators.hh"
 #include "element_class.hh"
 #include "mesh.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_MESH_INLINE_IMPL_HH_
 #define AKANTU_MESH_INLINE_IMPL_HH_
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 inline ElementKind Element::kind() const { return Mesh::getKind(type); }
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 template <typename... pack>
 Mesh::ElementTypesIteratorHelper Mesh::elementTypes(pack &&... _pack) const {
   return connectivities.elementTypes(_pack...);
 }
 
 /* -------------------------------------------------------------------------- */
 inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh,
                                             const std::string & origin)
     : MeshEvent<UInt>(origin),
       new_numbering(mesh.getNbNodes(), 1, "new_numbering") {}
 
 /* -------------------------------------------------------------------------- */
 inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh,
                                                   const ID & new_numbering_id,
                                                   const std::string & origin)
     : MeshEvent<Element>(origin),
       new_numbering(new_numbering_id, mesh.getID()) {}
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void Mesh::sendEvent<NewElementsEvent>(NewElementsEvent & event) {
   this->fillNodesToElements();
   EventHandlerManager<MeshEventHandler>::sendEvent(event);
 }
 
 /* -------------------------------------------------------------------------- */
 template <> inline void Mesh::sendEvent<NewNodesEvent>(NewNodesEvent & event) {
   this->computeBoundingBox();
   this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal);
   EventHandlerManager<MeshEventHandler>::sendEvent(event);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void
 Mesh::sendEvent<RemovedElementsEvent>(RemovedElementsEvent & event) {
   this->connectivities.onElementsRemoved(event.getNewNumbering());
   this->fillNodesToElements();
   this->computeBoundingBox();
 
   EventHandlerManager<MeshEventHandler>::sendEvent(event);
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 inline void Mesh::sendEvent<RemovedNodesEvent>(RemovedNodesEvent & event) {
   const auto & new_numbering = event.getNewNumbering();
   this->removeNodesFromArray(*nodes, new_numbering);
   if (nodes_global_ids and not is_mesh_facets) {
     this->removeNodesFromArray(*nodes_global_ids, new_numbering);
   }
   if (not is_mesh_facets) {
     this->removeNodesFromArray(*nodes_flags, new_numbering);
   }
 
   if (not nodes_to_elements.empty()) {
     std::vector<std::unique_ptr<std::set<Element>>> tmp(
         nodes_to_elements.size());
     auto it = nodes_to_elements.begin();
 
     UInt new_nb_nodes = 0;
     for (auto new_i : new_numbering) {
       if (new_i != UInt(-1)) {
         tmp[new_i] = std::move(*it);
         ++new_nb_nodes;
       }
       ++it;
     }
 
     tmp.resize(new_nb_nodes);
     std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin());
   }
 
   computeBoundingBox();
 
   EventHandlerManager<MeshEventHandler>::sendEvent(event);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline void Mesh::removeNodesFromArray(Array<T> & vect,
                                        const Array<UInt> & new_numbering) {
   Array<T> tmp(vect.size(), vect.getNbComponent());
   UInt nb_component = vect.getNbComponent();
   UInt new_nb_nodes = 0;
   for (UInt i = 0; i < new_numbering.size(); ++i) {
     UInt new_i = new_numbering(i);
     if (new_i != UInt(-1)) {
       T * to_copy = vect.storage() + i * nb_component;
       std::uninitialized_copy(to_copy, to_copy + nb_component,
                               tmp.storage() + new_i * nb_component);
       ++new_nb_nodes;
     }
   }
 
   tmp.resize(new_nb_nodes);
   vect.copy(tmp);
 }
 
 /* -------------------------------------------------------------------------- */
 inline Array<UInt> & Mesh::getNodesGlobalIdsPointer() {
   AKANTU_DEBUG_IN();
   if (not nodes_global_ids) {
     nodes_global_ids = std::make_shared<Array<UInt>>(
         nodes->size(), 1, getID() + ":nodes_global_ids");
 
     for (auto && global_ids : enumerate(*nodes_global_ids)) {
       std::get<1>(global_ids) = std::get<0>(global_ids);
     }
   }
 
   AKANTU_DEBUG_OUT();
   return *nodes_global_ids;
 }
 
 /* -------------------------------------------------------------------------- */
 inline Array<UInt> & Mesh::getConnectivityPointer(ElementType type,
                                                   GhostType ghost_type) {
   if (connectivities.exists(type, ghost_type)) {
     return connectivities(type, ghost_type);
   }
 
   if (ghost_type != _not_ghost) {
     ghosts_counters.alloc(0, 1, type, ghost_type, 1);
   }
 
   AKANTU_DEBUG_INFO("The connectivity vector for the type " << type
                                                             << " created");
 
   UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
   return connectivities.alloc(0, nb_nodes_per_element, type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline Array<std::vector<Element>> &
 Mesh::getElementToSubelementPointer(ElementType type, GhostType ghost_type) {
   return getDataPointer<std::vector<Element>>("element_to_subelement", type,
                                               ghost_type, 1, true);
 }
 
 /* -------------------------------------------------------------------------- */
 inline Array<Element> &
 Mesh::getSubelementToElementPointer(ElementType type, GhostType ghost_type) {
   auto & array = getDataPointer<Element>(
       "subelement_to_element", type, ghost_type, getNbFacetsPerElement(type),
       false, is_mesh_facets, ElementNull);
   return array;
 }
 
 /* -------------------------------------------------------------------------- */
 inline const auto & Mesh::getElementToSubelement() const {
   return getData<std::vector<Element>>("element_to_subelement");
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto & Mesh::getElementToSubelementNC() {
   return getData<std::vector<Element>>("element_to_subelement");
 }
 
 /* -------------------------------------------------------------------------- */
 inline const auto & Mesh::getElementToSubelement(ElementType type,
                                                  GhostType ghost_type) const {
   return getData<std::vector<Element>>("element_to_subelement", type,
                                        ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto & Mesh::getElementToSubelementNC(ElementType type,
                                              GhostType ghost_type) {
   return getData<std::vector<Element>>("element_to_subelement", type,
                                        ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline const auto &
 Mesh::getElementToSubelement(const Element & element) const {
   return getData<std::vector<Element>>("element_to_subelement")(element, 0);
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto & Mesh::getElementToSubelementNC(const Element & element) {
   return getData<std::vector<Element>>("element_to_subelement")(element, 0);
 }
 
 /* -------------------------------------------------------------------------- */
 inline const auto & Mesh::getSubelementToElement() const {
   return getData<Element>("subelement_to_element");
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto & Mesh::getSubelementToElementNC() {
   return getData<Element>("subelement_to_element");
 }
 
 /* -------------------------------------------------------------------------- */
 inline const auto & Mesh::getSubelementToElement(ElementType type,
                                                  GhostType ghost_type) const {
   return getData<Element>("subelement_to_element", type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto & Mesh::getSubelementToElementNC(ElementType type,
                                              GhostType ghost_type) {
   return getData<Element>("subelement_to_element", type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline VectorProxy<Element>
 Mesh::getSubelementToElement(const Element & element) const {
   return this->getSubelementToElement().get(element);
 }
 
 /* -------------------------------------------------------------------------- */
 inline VectorProxy<Element>
 Mesh::getSubelementToElementNC(const Element & element) const {
   return this->getSubelementToElement().get(element);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline Array<T> &
 Mesh::getDataPointer(const ID & data_name, ElementType el_type,
                      GhostType ghost_type, UInt nb_component,
                      bool size_to_nb_element, bool resize_with_parent) {
   Array<T> & tmp = this->getElementalDataArrayAlloc<T>(
       data_name, el_type, ghost_type, nb_component);
 
   if (size_to_nb_element) {
     if (resize_with_parent) {
       tmp.resize(mesh_parent->getNbElement(el_type, ghost_type));
     } else {
       tmp.resize(this->getNbElement(el_type, ghost_type));
     }
   }
 
   return tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline Array<T> &
 Mesh::getDataPointer(const ID & data_name, ElementType el_type,
                      GhostType ghost_type, UInt nb_component,
                      bool size_to_nb_element, bool resize_with_parent,
                      const T & defaul_) {
   Array<T> & tmp = this->getElementalDataArrayAlloc<T>(
       data_name, el_type, ghost_type, nb_component);
 
   if (size_to_nb_element) {
     if (resize_with_parent) {
       tmp.resize(mesh_parent->getNbElement(el_type, ghost_type), defaul_);
     } else {
       tmp.resize(this->getNbElement(el_type, ghost_type), defaul_);
     }
   }
 
   return tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline const Array<T> & Mesh::getData(const ID & data_name, ElementType el_type,
                                       GhostType ghost_type) const {
   return this->getElementalDataArray<T>(data_name, el_type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline Array<T> & Mesh::getData(const ID & data_name, ElementType el_type,
                                 GhostType ghost_type) {
   return this->getElementalDataArray<T>(data_name, el_type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline const ElementTypeMapArray<T> &
 Mesh::getData(const ID & data_name) const {
   return this->getElementalData<T>(data_name);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline ElementTypeMapArray<T> & Mesh::getData(const ID & data_name) {
   return this->getElementalData<T>(data_name);
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbElement(ElementType type, GhostType ghost_type) const {
   try {
 
     const Array<UInt> & conn = connectivities(type, ghost_type);
     return conn.size();
   } catch (...) {
     return 0;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbElement(const UInt spatial_dimension,
                                GhostType ghost_type, ElementKind kind) const {
   AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1),
                       "spatial_dimension is " << spatial_dimension
                                               << " and is greater than 3 !");
   UInt nb_element = 0;
 
   for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) {
     nb_element += getNbElement(type, ghost_type);
   }
 
   return nb_element;
 }
 
 /* -------------------------------------------------------------------------- */
 inline void Mesh::getBarycenter(const Element & element,
                                 Vector<Real> & barycenter) const {
   Vector<UInt> conn = getConnectivity(element);
   Matrix<Real> local_coord(spatial_dimension, conn.size());
   auto node_begin = make_view(*nodes, spatial_dimension).begin();
 
   for (auto && node : enumerate(conn)) {
     local_coord(std::get<0>(node)) =
         Vector<Real>(node_begin[std::get<1>(node)]);
   }
 
   Math::barycenter(local_coord.storage(), conn.size(), spatial_dimension,
                    barycenter.storage());
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbNodesPerElement(ElementType type) {
   UInt nb_nodes_per_element = 0;
 #define GET_NB_NODES_PER_ELEMENT(type)                                         \
   nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement()
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT);
 #undef GET_NB_NODES_PER_ELEMENT
   return nb_nodes_per_element;
 }
 
 /* -------------------------------------------------------------------------- */
 inline ElementType Mesh::getP1ElementType(ElementType type) {
   ElementType p1_type = _not_defined;
 #define GET_P1_TYPE(type) p1_type = ElementClass<type>::getP1ElementType()
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE);
 #undef GET_P1_TYPE
   return p1_type;
 }
 
 /* -------------------------------------------------------------------------- */
 inline ElementKind Mesh::getKind(ElementType type) {
   ElementKind kind = _ek_not_defined;
 #define GET_KIND(type) kind = ElementClass<type>::getKind()
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND);
 #undef GET_KIND
   return kind;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getSpatialDimension(ElementType type) {
   UInt spatial_dimension = 0;
 #define GET_SPATIAL_DIMENSION(type)                                            \
   spatial_dimension = ElementClass<type>::getSpatialDimension()
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION);
 #undef GET_SPATIAL_DIMENSION
 
   return spatial_dimension;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNaturalSpaceDimension(const ElementType & type) {
   UInt natural_dimension = 0;
 #define GET_NATURAL_DIMENSION(type)                                            \
   natural_dimension = ElementClass<type>::getNaturalSpaceDimension()
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_DIMENSION);
 #undef GET_NATURAL_DIMENSION
 
   return natural_dimension;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbFacetTypes(ElementType type, UInt /*t*/) {
   UInt nb = 0;
 #define GET_NB_FACET_TYPE(type) nb = ElementClass<type>::getNbFacetTypes()
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE);
 #undef GET_NB_FACET_TYPE
   return nb;
 }
 
 /* -------------------------------------------------------------------------- */
 inline constexpr auto Mesh::getFacetType(ElementType type, UInt t) {
 #define GET_FACET_TYPE(type) return ElementClass<type>::getFacetType(t);
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE);
 
 #undef GET_FACET_TYPE
 
   return _not_defined;
 }
 
 /* -------------------------------------------------------------------------- */
 inline constexpr auto Mesh::getAllFacetTypes(ElementType type) {
 #define GET_FACET_TYPE(type) return ElementClass<type>::getFacetTypes();
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH_NO_DEFAULT(GET_FACET_TYPE);
 #undef GET_FACET_TYPE
 
   return ElementClass<_not_defined>::getFacetTypes();
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbFacetsPerElement(ElementType type) {
   AKANTU_DEBUG_IN();
 
   UInt n_facet = 0;
 #define GET_NB_FACET(type) n_facet = ElementClass<type>::getNbFacetsPerElement()
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET);
 #undef GET_NB_FACET
 
   AKANTU_DEBUG_OUT();
   return n_facet;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbFacetsPerElement(ElementType type, UInt t) {
   AKANTU_DEBUG_IN();
 
   UInt n_facet = 0;
 #define GET_NB_FACET(type)                                                     \
   n_facet = ElementClass<type>::getNbFacetsPerElement(t)
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET);
 #undef GET_NB_FACET
 
   AKANTU_DEBUG_OUT();
   return n_facet;
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto Mesh::getFacetLocalConnectivity(ElementType type, UInt t) {
   AKANTU_DEBUG_IN();
 
 #define GET_FACET_CON(type)                                                    \
   AKANTU_DEBUG_OUT();                                                          \
   return ElementClass<type>::getFacetLocalConnectivityPerElement(t)
 
   AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON);
 #undef GET_FACET_CON
 
   AKANTU_DEBUG_OUT();
   return ElementClass<_not_defined>::getFacetLocalConnectivityPerElement(0);
   // This avoid a compilation warning but will certainly
   // also cause a segfault if reached
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto Mesh::getFacetConnectivity(const Element & element, UInt t) const {
   AKANTU_DEBUG_IN();
 
   Matrix<const UInt> local_facets(getFacetLocalConnectivity(element.type, t));
   Matrix<UInt> facets(local_facets.rows(), local_facets.cols());
 
   const Array<UInt> & conn = connectivities(element.type, element.ghost_type);
 
   for (UInt f = 0; f < facets.rows(); ++f) {
     for (UInt n = 0; n < facets.cols(); ++n) {
       facets(f, n) = conn(element.element, local_facets(f, n));
     }
   }
 
   AKANTU_DEBUG_OUT();
   return facets;
 }
 
 /* -------------------------------------------------------------------------- */
 inline VectorProxy<UInt> Mesh::getConnectivity(const Element & element) const {
   return connectivities.get(element);
 }
 
 /* -------------------------------------------------------------------------- */
 inline VectorProxy<UInt> Mesh::getConnectivityNC(const Element & element) {
   return connectivities.get(element);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 inline void Mesh::extractNodalValuesFromElement(
     const Array<T> & nodal_values, T * local_coord, const UInt * connectivity,
     UInt n_nodes, UInt nb_degree_of_freedom) const {
   for (UInt n = 0; n < n_nodes; ++n) {
     memcpy(local_coord + n * nb_degree_of_freedom,
            nodal_values.storage() + connectivity[n] * nb_degree_of_freedom,
            nb_degree_of_freedom * sizeof(T));
   }
 }
 
 /* -------------------------------------------------------------------------- */
 inline void Mesh::addConnectivityType(ElementType type, GhostType ghost_type) {
   getConnectivityPointer(type, ghost_type);
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isPureGhostNode(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_pure_ghost;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isLocalOrMasterNode(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_local_master_mask) == NodeFlag::_normal;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isLocalNode(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_normal;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isMasterNode(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_master;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isSlaveNode(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_shared_mask) == NodeFlag::_slave;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isPeriodicSlave(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) ==
          NodeFlag::_periodic_slave;
 }
 
 /* -------------------------------------------------------------------------- */
 inline bool Mesh::isPeriodicMaster(UInt n) const {
   return ((*nodes_flags)(n)&NodeFlag::_periodic_mask) ==
          NodeFlag::_periodic_master;
 }
 
 /* -------------------------------------------------------------------------- */
 inline NodeFlag Mesh::getNodeFlag(UInt local_id) const {
   return (*nodes_flags)(local_id);
 }
 
 /* -------------------------------------------------------------------------- */
 inline Int Mesh::getNodePrank(UInt local_id) const {
   auto it = nodes_prank.find(local_id);
   return it == nodes_prank.end() ? -1 : it->second;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNodeGlobalId(UInt local_id) const {
   return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id;
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNodeLocalId(UInt global_id) const {
   if (nodes_global_ids == nullptr) {
     return global_id;
   }
   return nodes_global_ids->find(global_id);
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbGlobalNodes() const {
   return nodes_global_ids ? nb_global_nodes : nodes->size();
 }
 
 /* -------------------------------------------------------------------------- */
 inline UInt Mesh::getNbNodesPerElementList(const Array<Element> & elements) {
   UInt nb_nodes_per_element = 0;
   UInt nb_nodes = 0;
   ElementType current_element_type = _not_defined;
 
   for (const auto & el : elements) {
     if (el.type != current_element_type) {
       current_element_type = el.type;
       nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type);
     }
 
     nb_nodes += nb_nodes_per_element;
   }
 
   return nb_nodes;
 }
 
 /* -------------------------------------------------------------------------- */
 inline Mesh & Mesh::getMeshFacets() {
   if (this->mesh_facets == nullptr) {
     AKANTU_SILENT_EXCEPTION(
         "No facet mesh is defined yet! check the buildFacets functions");
   }
 
   return *this->mesh_facets;
 }
 
 /* -------------------------------------------------------------------------- */
 inline const Mesh & Mesh::getMeshFacets() const {
   if (this->mesh_facets == nullptr) {
     AKANTU_SILENT_EXCEPTION(
         "No facet mesh is defined yet! check the buildFacets functions");
   }
 
   return *this->mesh_facets;
 }
 /* -------------------------------------------------------------------------- */
 inline const Mesh & Mesh::getMeshParent() const {
   if (this->mesh_parent == nullptr) {
     AKANTU_SILENT_EXCEPTION(
         "No parent mesh is defined! This is only valid in a mesh_facets");
   }
 
   return *this->mesh_parent;
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::addPeriodicSlave(UInt slave, UInt master) {
   if (master == slave) {
     return;
   }
 
   // if pair already registered
   auto master_slaves = periodic_master_slave.equal_range(master);
   auto slave_it =
       std::find_if(master_slaves.first, master_slaves.second,
                    [&](auto & pair) { return pair.second == slave; });
   if (slave_it == master_slaves.second) {
     // no duplicates
     periodic_master_slave.insert(std::make_pair(master, slave));
     AKANTU_DEBUG_INFO("adding periodic slave, slave gid:"
                       << getNodeGlobalId(slave) << " [lid: " << slave << "]"
                       << ", master gid:" << getNodeGlobalId(master)
                       << " [lid: " << master << "]");
     // std::cout << "adding periodic slave, slave gid:" <<
     // getNodeGlobalId(slave)
     //           << " [lid: " << slave << "]"
     //           << ", master gid:" << getNodeGlobalId(master)
     //           << " [lid: " << master << "]" << std::endl;
   }
 
   periodic_slave_master[slave] = master;
 
   auto set_flag = [&](auto node, auto flag) {
     (*nodes_flags)[node] &= ~NodeFlag::_periodic_mask; // clean periodic flags
     (*nodes_flags)[node] |= flag;
   };
 
   set_flag(slave, NodeFlag::_periodic_slave);
   set_flag(master, NodeFlag::_periodic_master);
 }
 
 /* -------------------------------------------------------------------------- */
 UInt Mesh::getPeriodicMaster(UInt slave) const {
   return periodic_slave_master.at(slave);
 }
 
 /* -------------------------------------------------------------------------- */
 class Mesh::PeriodicSlaves {
   using internal_iterator = std::unordered_multimap<UInt, UInt>::const_iterator;
   std::pair<internal_iterator, internal_iterator> pair;
 
 public:
   PeriodicSlaves(const Mesh & mesh, UInt master)
       : pair(mesh.getPeriodicMasterSlaves().equal_range(master)) {}
 
   PeriodicSlaves(const PeriodicSlaves & other) = default;
   PeriodicSlaves(PeriodicSlaves && other) = default;
   PeriodicSlaves & operator=(const PeriodicSlaves & other) = default;
 
   class const_iterator {
     internal_iterator it;
 
   public:
     const_iterator(internal_iterator it) : it(it) {}
 
     const_iterator operator++() {
       ++it;
       return *this;
     }
     bool operator!=(const const_iterator & other) { return other.it != it; }
+    bool operator==(const const_iterator & other) { return other.it == it; }
     auto operator*() { return it->second; }
   };
 
   auto begin() const { return const_iterator(pair.first); }
   auto end() const { return const_iterator(pair.second); }
 };
 
 /* -------------------------------------------------------------------------- */
 inline decltype(auto) Mesh::getPeriodicSlaves(UInt master) const {
   return PeriodicSlaves(*this, master);
 }
 
 /* -------------------------------------------------------------------------- */
 inline Vector<UInt>
 Mesh::getConnectivityWithPeriodicity(const Element & element) const {
   Vector<UInt> conn = getConnectivity(element);
   if (not isPeriodic()) {
     return conn;
   }
 
   for (auto && node : conn) {
     if (isPeriodicSlave(node)) {
       node = getPeriodicMaster(node);
     }
   }
 
   return conn;
 }
 
 } // namespace akantu
 
 #endif /* AKANTU_MESH_INLINE_IMPL_HH_ */