diff --git a/src/mesh/mesh.cc b/src/mesh/mesh.cc
index cd9cc9390..a5a2a0382 100644
--- a/src/mesh/mesh.cc
+++ b/src/mesh/mesh.cc
@@ -1,703 +1,710 @@
 /**
  * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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 "element_class.hh"
 #include "group_manager_inline_impl.hh"
 #include "mesh.hh"
 #include "mesh_global_data_updater.hh"
 #include "mesh_io.hh"
 #include "mesh_iterators.hh"
 #include "mesh_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include "communicator.hh"
 #include "element_synchronizer.hh"
 #include "facet_synchronizer.hh"
 #include "mesh_utils_distribution.hh"
 #include "node_synchronizer.hh"
 #include "periodic_node_synchronizer.hh"
 #if defined(AKANTU_COHESIVE_ELEMENT)
 #include "cohesive_element_inserter.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 /* -------------------------------------------------------------------------- */
 #include "dumper_field.hh"
 #include "dumper_internal_material_field.hh"
 /* -------------------------------------------------------------------------- */
 #include <limits>
 #include <sstream>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 Mesh::Mesh(Int spatial_dimension, const ID & id, Communicator & communicator)
     : GroupManager(*this, id + ":group_manager"), MeshData("mesh_data", id),
       id(id), connectivities("connectivities", id),
       ghosts_counters("ghosts_counters", id),
       spatial_dimension(spatial_dimension),
       size(Vector<double>::Zero(spatial_dimension)), bbox(spatial_dimension),
       bbox_local(spatial_dimension), communicator(&communicator) {
   AKANTU_DEBUG_IN();
   size.fill(0.);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Mesh::Mesh(Int spatial_dimension, Communicator & communicator, const ID & id)
     : Mesh(spatial_dimension, id, communicator) {
   AKANTU_DEBUG_IN();
 
   this->nodes =
       std::make_shared<Array<Real>>(0, spatial_dimension, id + ":coordinates");
   this->nodes_flags = std::make_shared<Array<NodeFlag>>(0, 1, NodeFlag::_normal,
                                                         id + ":nodes_flags");
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Mesh::Mesh(Int spatial_dimension, const ID & id)
     : Mesh(spatial_dimension, Communicator::getStaticCommunicator(), id) {}
 
 /* -------------------------------------------------------------------------- */
 Mesh::Mesh(Int spatial_dimension, const std::shared_ptr<Array<Real>> & nodes,
            const ID & id)
     : Mesh(spatial_dimension, id, Communicator::getStaticCommunicator()) {
   // NOLINTBEGIN(cppcoreguidelines-prefer-member-initializer)
   this->nodes = nodes;
   this->nb_global_nodes = this->nodes->size();
   // NOLINTEND(cppcoreguidelines-prefer-member-initializer)
   //
   this->nodes_to_elements.resize(nodes->size());
   for (auto & node_set : nodes_to_elements) {
     node_set = std::make_unique<std::set<Element>>();
   }
 
   this->computeBoundingBox();
 }
 
 /* -------------------------------------------------------------------------- */
 Mesh::~Mesh() = default;
 
+/* -------------------------------------------------------------------------- */
+template <>
+void Mesh::sendEvent<MeshIsDistributedEvent>(MeshIsDistributedEvent & event) {
+  //    if(event.getList().size() != 0)
+  EventHandlerManager<MeshEventHandler>::sendEvent(event);
+}
+
+/* -------------------------------------------------------------------------- */
+template <> void Mesh::sendEvent<NewNodesEvent>(NewNodesEvent & event) {
+  this->computeBoundingBox();
+  this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal);
+
+#if defined(AKANTU_COHESIVE_ELEMENT)
+  if (aka::is_of_type<CohesiveNewNodesEvent>(event)) {
+    // nodes might have changed in the connectivity
+    const auto & mesh_to_mesh_facet =
+        this->getData<Element>("mesh_to_mesh_facet");
+
+    for (auto ghost_type : ghost_types) {
+      for (auto type : connectivities.elementTypes(_spatial_dimension =
+                                                       spatial_dimension - 1,
+                       _ghost_type = ghost_type)) {
+        auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
+        if (not mesh_to_mesh_facet.exists(type, ghost_type)) {
+          continue;
+        }
+
+        const auto & mesh_to_mesh_facet_type =
+            mesh_to_mesh_facet(type, ghost_type);
+        auto && mesh_facet_conn_it =
+            make_view(mesh_facets->connectivities(type, ghost_type),
+                      nb_nodes_per_element)
+                .begin();
+
+        for (auto && [el, conn] : enumerate(make_view(
+                 connectivities(type, ghost_type), nb_nodes_per_element))) {
+          conn = mesh_facet_conn_it[mesh_to_mesh_facet_type(el).element];
+        }
+      }
+    }
+  }
+#endif
+
+  GroupManager::onNodesAdded(event.getList(), event);
+  EventHandlerManager<MeshEventHandler>::sendEvent(event);
+}
+
 /* -------------------------------------------------------------------------- */
 void Mesh::getBarycenters(Array<Real> & barycenter, ElementType type,
                           GhostType ghost_type) const {
   barycenter.resize(getNbElement(type, ghost_type));
   for (auto && data : enumerate(make_view(barycenter, spatial_dimension))) {
     getBarycenter(Element{type, Idx(std::get<0>(data)), ghost_type},
                   std::get<1>(data));
   }
 }
 
 class FacetGlobalConnectivityAccessor : public DataAccessor<Element> {
 public:
   FacetGlobalConnectivityAccessor(Mesh & mesh)
       : global_connectivity("global_connectivity",
                             "facet_connectivity_synchronizer") {
     global_connectivity.initialize(
         mesh, _spatial_dimension = _all_dimensions, _with_nb_element = true,
         _with_nb_nodes_per_element = true, _element_kind = _ek_regular);
     mesh.getGlobalConnectivity(global_connectivity);
   }
 
   [[nodiscard]] Int getNbData(const Array<Element> & elements,
                               const SynchronizationTag & tag) const override {
     Int size = 0;
     if (tag == SynchronizationTag::_smmc_facets_conn) {
       Int nb_nodes = Mesh::getNbNodesPerElementList(elements);
       size += nb_nodes * sizeof(Idx);
     }
     return size;
   }
 
   void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
                 const SynchronizationTag & tag) const override {
     if (tag == SynchronizationTag::_smmc_facets_conn) {
       for (const auto & element : elements) {
         const auto & conns =
             global_connectivity(element.type, element.ghost_type);
         for (auto n : arange(conns.getNbComponent())) {
           buffer << conns(element.element, n);
         }
       }
     }
   }
 
   void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
                   const SynchronizationTag & tag) override {
     if (tag == SynchronizationTag::_smmc_facets_conn) {
       for (const auto & element : elements) {
         auto & conns = global_connectivity(element.type, element.ghost_type);
         for (auto n : arange(conns.getNbComponent())) {
           buffer >> conns(element.element, n);
         }
       }
     }
   }
 
   AKANTU_GET_MACRO(GlobalConnectivity, (global_connectivity), decltype(auto));
 
 private:
   ElementTypeMapArray<Idx> global_connectivity;
 };
 
 /* -------------------------------------------------------------------------- */
 // const Array<Real> & Mesh::getNormals(ElementType element_type,
 //                                      GhostType ghost_type) {
 //   if (this->hasData<Real>("normals", element_type, ghost_type)) {
 //     return this->getData<Real>("normals", element_type, ghost_type);
 //   }
 
 //   auto & normals = getDataPointer<Real>("normals", element_type, ghost_type,
 //                                         spatial_dimension, true);
 //   for (auto && data [[gnu::unused]] :
 //        enumerate(make_view(normals, spatial_dimension))) {
 //     AKANTU_TO_IMPLEMENT();
 //   }
 
 //   AKANTU_TO_IMPLEMENT();
 // }
 
 /* -------------------------------------------------------------------------- */
 Mesh & Mesh::initMeshFacets(const ID & id) {
   AKANTU_DEBUG_IN();
 
   if (mesh_facets) {
     AKANTU_DEBUG_OUT();
     return *mesh_facets;
   }
 
   mesh_facets = std::make_unique<Mesh>(spatial_dimension, this->nodes,
                                        getID() + ":" + id);
   mesh_facets->mesh_parent = this;
   mesh_facets->is_mesh_facets = true;
   mesh_facets->nodes_flags = this->nodes_flags;
   mesh_facets->nodes_global_ids = this->nodes_global_ids;
 
   MeshUtils::buildAllFacets(*this, *mesh_facets, 0);
 
   if (mesh.isDistributed()) {
     mesh_facets->is_distributed = true;
     mesh_facets->element_synchronizer = std::make_unique<FacetSynchronizer>(
         *mesh_facets, mesh.getElementSynchronizer());
 
     FacetGlobalConnectivityAccessor data_accessor(*mesh_facets);
     /// communicate
     mesh_facets->element_synchronizer->synchronizeOnce(
         data_accessor, SynchronizationTag::_smmc_facets_conn);
 
     /// flip facets
     MeshUtils::flipFacets(*mesh_facets, data_accessor.getGlobalConnectivity(),
                           _ghost);
   }
 
   /// transfers the the mesh physical names to the mesh facets
   if (not this->hasData("physical_names")) {
     AKANTU_DEBUG_OUT();
     return *mesh_facets;
   }
 
   auto & mesh_phys_data = this->getData<std::string>("physical_names");
   auto & mesh_to_mesh_facet = this->getData<Element>("mesh_to_mesh_facet");
   mesh_to_mesh_facet.initialize(*this,
                                 _spatial_dimension = spatial_dimension - 1,
                                 _with_nb_element = true);
 
   auto & phys_data = mesh_facets->getData<std::string>("physical_names");
   phys_data.initialize(*mesh_facets, _spatial_dimension = spatial_dimension - 1,
                        _with_nb_element = true);
 
   ElementTypeMapArray<Real> barycenters(getID(), "temporary_barycenters");
   barycenters.initialize(*mesh_facets, _nb_component = spatial_dimension,
                          _spatial_dimension = spatial_dimension - 1,
                          _with_nb_element = true);
 
   for (auto && ghost_type : ghost_types) {
     for (auto && type :
          barycenters.elementTypes(spatial_dimension - 1, ghost_type)) {
       mesh_facets->getBarycenters(barycenters(type, ghost_type), type,
                                   ghost_type);
     }
   }
 
   for_each_element(
       mesh,
       [&](auto && element) {
         Vector<Real> barycenter(spatial_dimension);
         mesh.getBarycenter(element, barycenter);
         auto norm_barycenter = barycenter.norm();
         auto tolerance = Math::getTolerance();
         if (norm_barycenter > tolerance) {
           tolerance *= norm_barycenter;
         }
 
         Vector<Real> barycenter_facet(spatial_dimension);
 
         auto range = enumerate(make_view(
             barycenters(element.type, element.ghost_type), spatial_dimension));
 #ifndef AKANTU_NDEBUG
         auto min_dist = std::numeric_limits<Real>::max();
 #endif
         // this is a spacial search coded the most inefficient way.
         auto facet =
             std::find_if(range.begin(), range.end(), [&](auto && data) {
               auto norm_distance = barycenter.distance(std::get<1>(data));
 #ifndef AKANTU_NDEBUG
               min_dist = std::min(min_dist, norm_distance);
 #endif
               return (norm_distance < tolerance);
             });
 
         if (facet == range.end()) {
           AKANTU_DEBUG_INFO("The element "
                             << element
                             << " did not find its associated facet in the "
                                "mesh_facets! Try to decrease math tolerance. "
                                "The closest element was at a distance of "
                             << min_dist);
           return;
         }
 
         // set physical name
         auto && facet_element =
             Element{element.type, std::get<0>(*facet), element.ghost_type};
         phys_data(facet_element) = mesh_phys_data(element);
 
         mesh_to_mesh_facet(element) = facet_element;
       },
       _spatial_dimension = spatial_dimension - 1);
 
   mesh_facets->createGroupsFromMeshData<std::string>("physical_names");
 
   AKANTU_DEBUG_OUT();
   return *mesh_facets;
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::defineMeshParent(const Mesh & mesh) {
   AKANTU_DEBUG_IN();
 
   this->mesh_parent = &mesh;
   this->is_mesh_facets = true;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::read(const std::string & filename, const MeshIOType & mesh_io_type) {
 
   AKANTU_DEBUG_ASSERT(not is_distributed,
                       "You cannot read a mesh that is already distributed");
 
   MeshIO::read(filename, *this, mesh_io_type);
 
   auto types =
       this->elementTypes(spatial_dimension, _not_ghost, _ek_not_defined);
   auto it = types.begin();
   auto last = types.end();
   if (it == last) {
     AKANTU_DEBUG_WARNING(
         "The mesh contained in the file "
         << filename << " does not seem to be of the good dimension."
         << " No element of dimension " << spatial_dimension << " were read.");
   }
 
   this->makeReady();
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::write(const std::string & filename,
                  const MeshIOType & mesh_io_type) {
   MeshIO::write(filename, *this, mesh_io_type);
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::makeReady() {
   this->nb_global_nodes = this->nodes->size();
   this->computeBoundingBox();
   this->nodes_flags->resize(nodes->size(), NodeFlag::_normal);
   this->nodes_to_elements.resize(nodes->size());
   for (auto & node_set : nodes_to_elements) {
     node_set = std::make_unique<std::set<Element>>();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::printself(std::ostream & stream, int indent) const {
   std::string space(indent, AKANTU_INDENT);
 
   stream << space << "Mesh [" << std::endl;
   stream << space << " + id                : " << getID() << std::endl;
   stream << space << " + spatial dimension : " << this->spatial_dimension
          << std::endl;
   stream << space << " + nodes [" << std::endl;
   nodes->printself(stream, indent + 2);
   stream << space << " + connectivities [" << std::endl;
   connectivities.printself(stream, indent + 2);
   stream << space << " ]" << std::endl;
 
   GroupManager::printself(stream, indent + 1);
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::computeBoundingBox() {
   AKANTU_DEBUG_IN();
 
   bbox_local.reset();
 
   for (auto & pos : make_view(*nodes, spatial_dimension)) {
     //    if(!isPureGhostNode(i))
     bbox_local += pos;
   }
 
   if (this->is_distributed) {
     bbox = bbox_local.allSum(*communicator);
   } else {
     bbox = bbox_local;
   }
 
   size = bbox.size();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::getGlobalConnectivity(
     ElementTypeMapArray<Idx> & global_connectivity) {
   AKANTU_DEBUG_IN();
 
   for (auto && ghost_type : ghost_types) {
     for (auto type :
          global_connectivity.elementTypes(_spatial_dimension = _all_dimensions,
          _element_kind = _ek_not_defined, _ghost_type = ghost_type)) {
       if (not connectivities.exists(type, ghost_type)) {
         continue;
       }
 
       auto local_conn_view = make_view(connectivities(type, ghost_type));
       auto global_conn_view = make_view(global_connectivity(type, ghost_type));
 
       std::transform(local_conn_view.begin(), local_conn_view.end(),
                      global_conn_view.begin(),
                      [&](Idx l) -> Idx { return this->getNodeGlobalId(l); });
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 DumperIOHelper & Mesh::getGroupDumper(const std::string & dumper_name,
                                       const std::string & group_name) {
   if (group_name == "all") {
     return this->getDumper(dumper_name);
   }
   return element_groups[group_name]->getDumper(dumper_name);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 ElementTypeMap<Int> Mesh::getNbDataPerElem(ElementTypeMapArray<T> & arrays) {
   ElementTypeMap<Int> nb_data_per_elem;
 
   for (auto type : arrays.elementTypes(_element_kind = _ek_not_defined)) {
     auto nb_elements = this->getNbElement(type);
     auto & array = arrays(type);
 
     nb_data_per_elem(type) = array.getNbComponent() * array.size();
     nb_data_per_elem(type) /= nb_elements;
   }
 
   return nb_data_per_elem;
 }
 
 /* -------------------------------------------------------------------------- */
 template ElementTypeMap<Int>
 Mesh::getNbDataPerElem(ElementTypeMapArray<Real> & array);
 
 template ElementTypeMap<Int>
 Mesh::getNbDataPerElem(ElementTypeMapArray<Int> & array);
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 std::shared_ptr<dumpers::Field>
 Mesh::createFieldFromAttachedData(const std::string & field_id,
                                   const std::string & group_name,
                                   ElementKind element_kind) {
 
   std::shared_ptr<dumpers::Field> field;
   ElementTypeMapArray<T> * internal = nullptr;
   try {
     internal = &(this->getData<T>(field_id));
   } catch (...) {
     return nullptr;
   }
 
   auto && nb_data_per_elem = this->getNbDataPerElem(*internal);
 
   field = this->createElementalField<T, dumpers::InternalMaterialField>(
       *internal, group_name, this->spatial_dimension, element_kind,
       nb_data_per_elem);
 
   return field;
 }
 
 template std::shared_ptr<dumpers::Field>
 Mesh::createFieldFromAttachedData<Real>(const std::string & field_id,
                                         const std::string & group_name,
                                         ElementKind element_kind);
 
 template std::shared_ptr<dumpers::Field>
 Mesh::createFieldFromAttachedData<Int>(const std::string & field_id,
                                        const std::string & group_name,
                                        ElementKind element_kind);
 
 /* -------------------------------------------------------------------------- */
 void Mesh::distributeImpl(
     Communicator & communicator,
     const std::function<Int(const Element &, const Element &)> &
         edge_weight_function [[gnu::unused]],
     const std::function<Int(const Element &)> & vertex_weight_function
     [[gnu::unused]]) {
   AKANTU_DEBUG_ASSERT(is_distributed == false,
                       "This mesh is already distribute");
   this->communicator = &communicator;
 
   this->element_synchronizer = std::make_unique<ElementSynchronizer>(
       *this, this->getID() + ":element_synchronizer", true);
 
   this->node_synchronizer = std::make_unique<NodeSynchronizer>(
       *this, this->getID() + ":node_synchronizer", true);
 
   auto psize = this->communicator->getNbProc();
 
   if (psize > 1) {
 #ifdef AKANTU_USE_SCOTCH
     auto prank = this->communicator->whoAmI();
     if (prank == 0) {
       MeshPartitionScotch partition(*this, spatial_dimension);
       partition.partitionate(psize, edge_weight_function,
                              vertex_weight_function);
 
       MeshUtilsDistribution::distributeMeshCentralized(*this, 0, partition);
     } else {
       MeshUtilsDistribution::distributeMeshCentralized(*this, 0);
     }
 
 #else
     AKANTU_ERROR("Cannot distribute a mesh without a partitioning tool");
 #endif
   }
 
   // if (psize > 1)
   this->is_distributed = true;
 
   this->computeBoundingBox();
 
   MeshIsDistributedEvent event(AKANTU_CURRENT_FUNCTION);
   this->sendEvent(event);
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::getAssociatedElements(const Array<Idx> & node_list,
                                  Array<Element> & elements) {
   for (const auto & node : node_list) {
     for (const auto & element : *nodes_to_elements[node]) {
       elements.push_back(element);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::getAssociatedElements(const Idx & node,
                                  Array<Element> & elements) const {
   for (const auto & element : *nodes_to_elements[node]) {
     elements.push_back(element);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::fillNodesToElements(Int dimension) {
   Element element{ElementNull};
 
   auto nb_nodes = nodes->size();
   this->nodes_to_elements.resize(nb_nodes);
 
   for (Int n = 0; n < nb_nodes; ++n) {
     if (this->nodes_to_elements[n]) {
       this->nodes_to_elements[n]->clear();
     } else {
       this->nodes_to_elements[n] = std::make_unique<std::set<Element>>();
     }
   }
 
   for (auto ghost_type : ghost_types) {
     element.ghost_type = ghost_type;
     for (const auto & type :
          elementTypes(dimension, ghost_type, _ek_not_defined)) {
       element.type = type;
 
       auto nb_element = this->getNbElement(type, ghost_type);
       auto connectivity = connectivities(type, ghost_type);
       auto conn_it = connectivity.begin(connectivity.getNbComponent());
 
       for (Int el = 0; el < nb_element; ++el, ++conn_it) {
         element.element = el;
         const auto & conn = *conn_it;
         for (auto node : conn) {
           nodes_to_elements[node]->insert(element);
         }
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 std::tuple<Idx, Idx> Mesh::updateGlobalData(NewNodesEvent & nodes_event,
                                             NewElementsEvent & elements_event) {
   if (global_data_updater) {
     return this->global_data_updater->updateData(nodes_event, elements_event);
   }
   return std::make_tuple(nodes_event.getList().size(),
                          elements_event.getList().size());
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::registerGlobalDataUpdater(
     std::unique_ptr<MeshGlobalDataUpdater> && global_data_updater) {
   this->global_data_updater = std::move(global_data_updater);
 }
 
 /* -------------------------------------------------------------------------- */
 void Mesh::eraseElements(const Array<Element> & elements) {
   ElementTypeMap<Idx> last_element;
 
   RemovedElementsEvent event(*this, "new_numbering", AKANTU_CURRENT_FUNCTION);
   auto & remove_list = event.getList();
   auto & new_numbering = event.getNewNumbering();
 
   for (auto && el : elements) {
     if (el.ghost_type != _not_ghost) {
       auto & count = ghosts_counters(el);
       --count;
       AKANTU_DEBUG_ASSERT(count >= 0,
                           "Something went wrong in the ghost element counter");
 
       if (count > 0) {
         continue;
       }
     }
     remove_list.push_back(el);
     if (not new_numbering.exists(el.type, el.ghost_type)) {
       auto nb_element = mesh.getNbElement(el.type, el.ghost_type);
       auto & numbering =
           new_numbering.alloc(nb_element, 1, el.type, el.ghost_type);
       for (auto && [i, numb] : enumerate(numbering)) {
         numb = i;
       }
     }
 
     new_numbering(el) = -1;
   }
 
   auto find_last_not_deleted = [](auto && array, Int start) -> Int {
     do {
       --start;
     } while (start >= 0 and array[start] == -1);
 
     return start;
   };
 
   auto find_first_deleted = [](auto && array, Int start) -> Int {
     auto begin = array.begin();
     auto it = std::find_if(begin + start, array.end(),
                            [](auto & el) { return el == -1; });
     return Int(it - begin);
   };
 
   for (auto ghost_type : ghost_types) {
     for (auto type : new_numbering.elementTypes(_ghost_type = ghost_type)) {
       auto & numbering = new_numbering(type, ghost_type);
       auto last_not_delete = find_last_not_deleted(numbering, numbering.size());
 
       if (last_not_delete < 0) {
         continue;
       }
 
       auto pos = find_first_deleted(numbering, 0);
 
       while (pos < last_not_delete) {
         std::swap(numbering[pos], numbering[last_not_delete]);
         last_not_delete = find_last_not_deleted(numbering, last_not_delete);
         pos = find_first_deleted(numbering, pos + 1);
       }
     }
   }
 
   this->ghosts_counters.onElementsRemoved(new_numbering);
   this->sendEvent(event);
 }
 
-/* -------------------------------------------------------------------------- */
-template <> inline void Mesh::sendEvent<NewNodesEvent>(NewNodesEvent & event) {
-  this->computeBoundingBox();
-  this->nodes_flags->resize(this->nodes->size(), NodeFlag::_normal);
-
-#if defined(AKANTU_COHESIVE_ELEMENT)
-  if (aka::is_of_type<CohesiveNewNodesEvent>(event)) {
-    // nodes might have changed in the connectivity
-    const auto & mesh_to_mesh_facet =
-        this->getData<Element>("mesh_to_mesh_facet");
-
-    for (auto ghost_type : ghost_types) {
-      for (auto type : connectivities.elementTypes(_spatial_dimension =
-                                                       spatial_dimension - 1,
-                       _ghost_type = ghost_type)) {
-        auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
-        if (not mesh_to_mesh_facet.exists(type, ghost_type)) {
-          continue;
-        }
-
-        const auto & mesh_to_mesh_facet_type =
-            mesh_to_mesh_facet(type, ghost_type);
-        auto && mesh_facet_conn_it =
-            make_view(mesh_facets->connectivities(type, ghost_type),
-                      nb_nodes_per_element)
-                .begin();
-
-        for (auto && [el, conn] : enumerate(make_view(
-                 connectivities(type, ghost_type), nb_nodes_per_element))) {
-          conn = mesh_facet_conn_it[mesh_to_mesh_facet_type(el).element];
-        }
-      }
-    }
-  }
-#endif
-
-  GroupManager::onNodesAdded(event.getList(), event);
-  EventHandlerManager<MeshEventHandler>::sendEvent(event);
-}
-
 } // namespace akantu
diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh
index 0acf14980..27c0671c8 100644
--- a/src/mesh/mesh.hh
+++ b/src/mesh/mesh.hh
@@ -1,712 +1,709 @@
 /**
  * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * This file is part of Akantu
  *
  * 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/>.
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_MESH_HH_
 #define AKANTU_MESH_HH_
 
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 #include "aka_bbox.hh"
 #include "aka_event_handler_manager.hh"
 #include "communicator.hh"
 #include "dumpable.hh"
 #include "element.hh"
 #include "element_class.hh"
 #include "element_type_map.hh"
 #include "group_manager.hh"
 #include "mesh_data.hh"
 #include "mesh_events.hh"
 /* -------------------------------------------------------------------------- */
 #include <functional>
 #include <set>
 #include <unordered_map>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 class ElementSynchronizer;
 class NodeSynchronizer;
 class PeriodicNodeSynchronizer;
 class MeshGlobalDataUpdater;
 } // namespace akantu
 
 namespace akantu {
 
 namespace {
   DECLARE_NAMED_ARGUMENT(communicator);
   DECLARE_NAMED_ARGUMENT(edge_weight_function);
   DECLARE_NAMED_ARGUMENT(vertex_weight_function);
 } // namespace
 
 /* -------------------------------------------------------------------------- */
 /* Mesh                                                                       */
 /* -------------------------------------------------------------------------- */
 
 /**
  * @class  Mesh mesh.hh
  *
  * This class contaisn the coordinates of the nodes in the Mesh.nodes
  * akant::Array, and the connectivity. The connectivity are stored in by element
  * types.
  *
  * In order to loop on all element you have to loop on all types like this :
  * @code{.cpp}
  for(auto & type : mesh.elementTypes()) {
    Int nb_element  = mesh.getNbElement(type);
    const auto & conn = mesh.getConnectivity(type);
    for(Int e = 0; e < nb_element; ++e) {
      ...
    }
  }
 
  or
 
  for_each_element(mesh, [](Element & element) {
     std::cout << element << std::endl
   });
  @endcode
 */
 class Mesh : public EventHandlerManager<MeshEventHandler>,
              public GroupManager,
              public MeshData,
              public Dumpable {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 private:
   /// default constructor used for chaining, the last parameter is just to
   /// differentiate constructors
   Mesh(Int spatial_dimension, const ID & id, Communicator & communicator);
 
 public:
   /// constructor that create nodes coordinates array
   Mesh(Int spatial_dimension, const ID & id = "mesh");
 
   /// mesh not distributed and not using the default communicator
   Mesh(Int spatial_dimension, Communicator & communicator,
        const ID & id = "mesh");
 
   /**
    * constructor that use an existing nodes coordinates
    * array, by getting the vector of coordinates
    */
   Mesh(Int spatial_dimension, const std::shared_ptr<Array<Real>> & nodes,
        const ID & id = "mesh");
 
   ~Mesh() override;
 
   Mesh(const Mesh &) = delete;
   Mesh(Mesh &&) = delete;
   Mesh & operator=(const Mesh &) = delete;
   Mesh & operator=(Mesh &&) = delete;
 
   /// read the mesh from a file
   void read(const std::string & filename,
             const MeshIOType & mesh_io_type = _miot_auto);
   /// write the mesh to a file
   void write(const std::string & filename,
              const MeshIOType & mesh_io_type = _miot_auto);
 
 protected:
   void makeReady();
 
 private:
   /// initialize the connectivity to NULL and other stuff
   void init();
 
   /// function that computes the bounding box (fills xmin, xmax)
   void computeBoundingBox();
 
   /* ------------------------------------------------------------------------ */
   /* Distributed memory methods and accessors                                 */
   /* ------------------------------------------------------------------------ */
 public:
 protected:
   /// patitionate the mesh among the processors involved in their computation
   virtual void distributeImpl(
       Communicator & communicator,
       const std::function<Int(const Element &, const Element &)> &
           edge_weight_function,
       const std::function<Int(const Element &)> & vertex_weight_function);
 
 public:
   /// with the arguments to pass to the partitionner
   template <typename... pack>
   std::enable_if_t<are_named_argument<pack...>::value>
   distribute(pack &&... _pack) {
     distributeImpl(
         OPTIONAL_NAMED_ARG(communicator, Communicator::getStaticCommunicator()),
         OPTIONAL_NAMED_ARG(edge_weight_function,
                            [](auto &&, auto &&) { return 1; }),
         OPTIONAL_NAMED_ARG(vertex_weight_function, [](auto &&) { return 1; }));
   }
 
   /// defines is the mesh is distributed or not
   inline bool isDistributed() const { return this->is_distributed; }
 
   /* ------------------------------------------------------------------------ */
   /* Periodicity methods and accessors                                        */
   /* ------------------------------------------------------------------------ */
 public:
   /// set the periodicity in a given direction
   void makePeriodic(const SpatialDirection & direction);
   void makePeriodic(const SpatialDirection & direction, const ID & list_1,
                     const ID & list_2);
 
 protected:
   void makePeriodic(const SpatialDirection & direction,
                     const Array<Idx> & list_1, const Array<Idx> & list_2);
 
   /// Removes the face that the mesh is periodic
   void wipePeriodicInfo();
 
   inline void addPeriodicSlave(Idx slave, Idx master);
 
   template <typename T>
   void synchronizePeriodicSlaveDataWithMaster(Array<T> & data);
 
   // update the periodic synchronizer (creates it if it does not exists)
   void updatePeriodicSynchronizer();
 
 public:
   /// defines if the mesh is periodic or not
   inline bool isPeriodic() const { return this->is_periodic; }
 
   inline bool isPeriodic(const SpatialDirection & /*direction*/) const {
     return this->is_periodic;
   }
 
   class PeriodicSlaves;
 
   /// get the master node for a given slave nodes, except if node not a slave
   inline Idx getPeriodicMaster(Idx slave) const;
 
   /// get an iterable list of slaves for a given master node
   inline decltype(auto) getPeriodicSlaves(Idx master) const;
 
   /* ------------------------------------------------------------------------ */
   /* General Methods                                                          */
   /* ------------------------------------------------------------------------ */
 public:
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /// extract coordinates of nodes from an element
   template <typename T, class Derived1, class Derived2,
             std::enable_if_t<aka::is_vector_v<Derived2>> * = nullptr>
   inline void extractNodalValuesFromElement(
       const Array<T> & nodal_values,
       Eigen::MatrixBase<Derived1> & elemental_values,
       const Eigen::MatrixBase<Derived2> & connectivity) const;
 
   /// extract coordinates of nodes from an element
   template <typename T>
   inline decltype(auto)
   extractNodalValuesFromElement(const Array<T> & nodal_values,
                                 const Element & element) const;
 
   /// add a Array of connectivity for the given ElementType and GhostType .
   inline void addConnectivityType(ElementType type,
                                   GhostType ghost_type = _not_ghost);
 
   /* ------------------------------------------------------------------------ */
-  template <class Event> inline void sendEvent(Event & event) {
-    //    if(event.getList().size() != 0)
-    EventHandlerManager<MeshEventHandler>::sendEvent<Event>(event);
-  }
+  template <class Event> void sendEvent(Event & event);
 
   /// prepare the  event to remove the elements listed
   void eraseElements(const Array<Element> & elements);
 
   /* ------------------------------------------------------------------------ */
   template <typename T>
   inline void removeNodesFromArray(Array<T> & vect,
                                    const Array<Int> & new_numbering);
 
   /// init facets' mesh
   Mesh & initMeshFacets(const ID & id = "mesh_facets");
 
   /// define parent mesh
   void defineMeshParent(const Mesh & mesh);
 
   /// get global connectivity array
   void getGlobalConnectivity(ElementTypeMapArray<Int> & global_connectivity);
 
 public:
   void getAssociatedElements(const Array<Int> & node_list,
                              Array<Element> & elements);
 
   void getAssociatedElements(const Idx & node, Array<Element> & elements) const;
 
   inline decltype(auto) getAssociatedElements(const Idx & node) const;
 
 public:
   /// fills the nodes_to_elements for given dimension elements
   void fillNodesToElements(Int dimension = _all_dimensions);
 
 private:
   /// update the global ids, nodes type, ...
   std::tuple<Int, Int> updateGlobalData(NewNodesEvent & nodes_event,
                                         NewElementsEvent & elements_event);
 
   void registerGlobalDataUpdater(
       std::unique_ptr<MeshGlobalDataUpdater> && global_data_updater);
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the id of the mesh
   AKANTU_GET_MACRO(ID, id, const ID &);
 
   /// get the spatial dimension of the mesh = number of component of the
   /// coordinates
   AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, Int);
 
   /// get the nodes Array aka coordinates
   AKANTU_GET_MACRO(Nodes, *nodes, const Array<Real> &);
   AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array<Real> &);
 
   /// get the number of nodes
   auto getNbNodes() const { return nodes->size(); }
 
   /// get the Array of global ids of the nodes (only used in parallel)
   AKANTU_GET_MACRO_AUTO(GlobalNodesIds, *nodes_global_ids);
   // AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array<UInt>
   // &);
 
   /// get the global id of a node
   inline auto getNodeGlobalId(Idx local_id) const;
 
   /// get the global id of a node
   inline auto getNodeLocalId(Idx global_id) const;
 
   /// get the global number of nodes
   inline auto getNbGlobalNodes() const;
 
   /// get the nodes type Array
   AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array<NodeFlag> &);
 
 protected:
   AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array<NodeFlag> &);
 
 public:
   inline NodeFlag getNodeFlag(Idx local_id) const;
   inline auto getNodePrank(Idx local_id) const;
 
   /// say if a node is a pure ghost node
   inline bool isPureGhostNode(Idx n) const;
 
   /// say if a node is pur local or master node
   inline bool isLocalOrMasterNode(Idx n) const;
 
   inline bool isLocalNode(Idx n) const;
   inline bool isMasterNode(Idx n) const;
   inline bool isSlaveNode(Idx n) const;
 
   inline bool isPeriodicSlave(Idx n) const;
   inline bool isPeriodicMaster(Idx n) const;
 
   const Vector<Real> & getLowerBounds() const { return bbox.getLowerBounds(); }
   const Vector<Real> & getUpperBounds() const { return bbox.getUpperBounds(); }
   AKANTU_GET_MACRO(BBox, bbox, const BBox &);
 
   const Vector<Real> & getLocalLowerBounds() const {
     return bbox_local.getLowerBounds();
   }
   const Vector<Real> & getLocalUpperBounds() const {
     return bbox_local.getUpperBounds();
   }
   AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &);
 
   /// get the connectivity Array for a given type
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, Idx);
   AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, Idx);
   AKANTU_GET_MACRO(Connectivities, connectivities,
                    const ElementTypeMapArray<Idx> &);
 
   /// get the number of element of a type in the mesh
   inline auto getNbElement(ElementType type,
                            GhostType ghost_type = _not_ghost) const;
 
   /// get the number of element for a given ghost_type and a given dimension
   inline auto getNbElement(Int spatial_dimension = _all_dimensions,
                            GhostType ghost_type = _not_ghost,
                            ElementKind kind = _ek_not_defined) const;
 
   /// compute the barycenter of a given element
   template <class D, std::enable_if_t<aka::is_vector_v<D>> * = nullptr>
   inline void getBarycenter(const Element & element,
                             const Eigen::MatrixBase<D> & barycenter) const;
 
   inline Vector<Real> getBarycenter(const Element & element) const;
 
   void getBarycenters(Array<Real> & barycenter, ElementType type,
                       GhostType ghost_type) const;
 
   /// get the element connected to a subelement (element of lower dimension)
   decltype(auto) getElementToSubelement() const;
 
   /// get the element connected to a subelement
   const auto & getElementToSubelement(ElementType el_type,
                                       GhostType ghost_type = _not_ghost) const;
 
   /// get the elements connected to a subelement
   decltype(auto) getElementToSubelement(const Element & element) const;
 
   /// get the subelement (element of lower dimension) connected to a element
   decltype(auto) getSubelementToElement() const;
 
   /// get the subelement connected to an element
   const auto & getSubelementToElement(ElementType el_type,
                                       GhostType ghost_type = _not_ghost) const;
 
   /// get the subelement (element of lower dimension) connected to a element
   decltype(auto) getSubelementToElement(const Element & element) const;
 
   /// get connectivity of a given element
   inline decltype(auto) getConnectivity(const Element & element) const;
   inline decltype(auto)
   getConnectivityWithPeriodicity(const Element & element) const;
 
 protected:
   /// get the element connected to a subelement (element of lower dimension)
   auto & getElementToSubelementNC();
   auto & getSubelementToElementNC();
   inline auto & getElementToSubelementNC(const Element & element);
   inline decltype(auto) getSubelementToElementNC(const Element & element);
   /// get the element connected to a subelement
   auto & getElementToSubelementNC(ElementType el_type,
                                   GhostType ghost_type = _not_ghost);
   /// get the subelement connected to an element
   auto & getSubelementToElementNC(ElementType el_type,
                                   GhostType ghost_type = _not_ghost);
 
   inline decltype(auto) getConnectivityNC(const Element & element);
 
 public:
   /// get a name field associated to the mesh
   template <typename T>
   inline decltype(auto) getData(const ID & data_name, ElementType el_type,
                                 GhostType ghost_type = _not_ghost) const;
 
   /// get a name field associated to the mesh
   template <typename T>
   inline decltype(auto) getData(const ID & data_name, ElementType el_type,
                                 GhostType ghost_type = _not_ghost);
 
   /// get a name field associated to the mesh
   template <typename T>
   inline decltype(auto) getData(const ID & data_name) const;
 
   /// get a name field associated to the mesh
   template <typename T> inline decltype(auto) getData(const ID & data_name);
 
   template <typename T>
   auto getNbDataPerElem(ElementTypeMapArray<T> & array) -> ElementTypeMap<Int>;
 
   template <typename T>
   std::shared_ptr<dumpers::Field>
   createFieldFromAttachedData(const std::string & field_id,
                               const std::string & group_name,
                               ElementKind element_kind);
 
   /// templated getter returning the pointer to data in MeshData (modifiable)
   template <typename T>
   inline decltype(auto)
   getDataPointer(const std::string & data_name, ElementType el_type,
                  GhostType ghost_type = _not_ghost, Int nb_component = 1,
                  bool size_to_nb_element = true,
                  bool resize_with_parent = false);
 
   template <typename T>
   inline decltype(auto)
   getDataPointer(const ID & data_name, ElementType el_type,
                  GhostType ghost_type, Int nb_component,
                  bool size_to_nb_element, bool resize_with_parent,
                  const T & defaul_);
 
   /// Facets mesh accessor
   inline auto getMeshFacets() const -> const Mesh &;
   inline auto getMeshFacets() -> Mesh &;
 
   inline auto hasMeshFacets() const { return mesh_facets != nullptr; }
 
   /// Parent mesh accessor
   inline auto getMeshParent() const -> const Mesh &;
 
   inline auto isMeshFacets() const { return this->is_mesh_facets; }
 
   /// return the dumper from a group and and a dumper name
   auto getGroupDumper(const std::string & dumper_name,
                       const std::string & group_name) -> DumperIOHelper &;
 
   /* ------------------------------------------------------------------------ */
   /* Wrappers on ElementClass functions                                       */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the number of nodes per element for a given element type
   static inline constexpr auto getNbNodesPerElement(ElementType type) -> Int;
 
   /// get the number of nodes per element for a given element type considered as
   /// a first order element
   static inline constexpr auto getP1ElementType(ElementType type)
       -> ElementType;
 
   /// get the kind of the element type
   static inline constexpr auto getKind(ElementType type) -> ElementKind;
 
   /// get spatial dimension of a type of element
   static inline constexpr auto getSpatialDimension(ElementType type) -> Int;
 
   /// get the natural space dimension of a type of element
   static inline constexpr auto getNaturalSpaceDimension(ElementType type)
       -> Int;
 
   /// get number of facets of a given element type
   static inline constexpr auto getNbFacetsPerElement(ElementType type) -> Int;
 
   /// get number of facets of a given element type
   static inline constexpr auto getNbFacetsPerElement(ElementType type, Idx t)
       -> Int;
 
   /// get local connectivity of a facet for a given facet type
   static inline decltype(auto) getFacetLocalConnectivity(ElementType type,
                                                          Idx t = 0);
 
   /// get connectivity of facets for a given element
   inline auto getFacetConnectivity(const Element & element, Idx t = 0) const
       -> Matrix<Idx>;
 
   /// get the number of type of the surface element associated to a given
   /// element type
   static inline constexpr auto getNbFacetTypes(ElementType type, Idx t = 0)
       -> Int;
 
   /// get the type of the surface element associated to a given element
   static inline constexpr auto getFacetType(ElementType type, Idx t = 0)
       -> ElementType;
 
   /// get all the type of the surface element associated to a given element
   static inline decltype(auto) getAllFacetTypes(ElementType type);
 
   /// get the number of nodes in the given element list
   static inline auto getNbNodesPerElementList(const Array<Element> & elements)
       -> Int;
 
   /* ------------------------------------------------------------------------ */
   /* Element type Iterator                                                    */
   /* ------------------------------------------------------------------------ */
   using ElementTypesIteratorHelper =
       ElementTypeMapArray<Idx, ElementType>::ElementTypesIteratorHelper;
 
   template <typename... pack>
   auto elementTypes(pack &&... _pack) const -> ElementTypesIteratorHelper;
 
   AKANTU_GET_MACRO_DEREF_PTR(ElementSynchronizer, element_synchronizer);
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ElementSynchronizer,
                                        element_synchronizer);
   AKANTU_GET_MACRO_DEREF_PTR(NodeSynchronizer, node_synchronizer);
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(NodeSynchronizer, node_synchronizer);
   AKANTU_GET_MACRO_DEREF_PTR(PeriodicNodeSynchronizer,
                              periodic_node_synchronizer);
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(PeriodicNodeSynchronizer,
                                        periodic_node_synchronizer);
 
   // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator
   // &);
   AKANTU_GET_MACRO_DEREF_PTR(Communicator, communicator);
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Communicator, communicator);
   AKANTU_GET_MACRO_AUTO(PeriodicMasterSlaves, periodic_master_slave);
 
   /* ------------------------------------------------------------------------ */
   /* Private methods for friends                                              */
   /* ------------------------------------------------------------------------ */
 private:
   friend class MeshAccessor;
   friend class MeshUtils;
 
   AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(NodesPointer, nodes);
 
   /// get a pointer to the nodes_global_ids Array<UInt> and create it if
   /// necessary
   inline auto getNodesGlobalIdsPointer() -> Array<Idx> &;
 
   /// get a pointer to the nodes_type Array<Int> and create it if necessary
   inline auto getNodesFlagsPointer() -> Array<NodeFlag> &;
 
   /// get a pointer to the connectivity Array for the given type and create it
   /// if necessary
   inline auto getConnectivityPointer(ElementType type,
                                      GhostType ghost_type = _not_ghost)
       -> Array<Idx> &;
 
   /// get the ghost element counter
   inline auto getGhostsCounters(ElementType type, GhostType ghost_type = _ghost)
       -> Array<Idx> & {
     AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost,
                         "No ghost counter for _not_ghost elements");
     return ghosts_counters(type, ghost_type);
   }
 
   /// get a pointer to the element_to_subelement Array for the given type and
   /// create it if necessary
   inline decltype(auto)
   getElementToSubelementPointer(ElementType type,
                                 GhostType ghost_type = _not_ghost);
 
   /// get a pointer to the subelement_to_element Array for the given type and
   /// create it if necessary
   inline decltype(auto)
   getSubelementToElementPointer(ElementType type,
                                 GhostType ghost_type = _not_ghost);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
   ID id;
 
   /// array of the nodes coordinates
   std::shared_ptr<Array<Real>> nodes;
 
   /// global node ids
   std::shared_ptr<Array<Idx>> nodes_global_ids;
 
   /// node flags (shared/periodic/...)
   std::shared_ptr<Array<NodeFlag>> nodes_flags;
 
   /// processor handling the node when not local or master
   std::unordered_map<Idx, Int> nodes_prank;
 
   /// global number of nodes;
   Int nb_global_nodes{0};
 
   /// all class of elements present in this mesh (for heterogenous meshes)
   ElementTypeMapArray<Idx> connectivities;
 
   /// count the references on ghost elements
   ElementTypeMapArray<Idx> ghosts_counters;
 
   /// the spatial dimension of this mesh
   Idx spatial_dimension{0};
 
   /// size covered by the mesh on each direction
   Vector<Real> size;
 
   /// global bounding box
   BBox bbox;
 
   /// local bounding box
   BBox bbox_local;
 
   /// Extra data loaded from the mesh file
   // MeshData mesh_data;
 
   /// facets' mesh
   std::unique_ptr<Mesh> mesh_facets;
 
   /// parent mesh (this is set for mesh_facets meshes)
   const Mesh * mesh_parent{nullptr};
 
   /// defines if current mesh is mesh_facets or not
   bool is_mesh_facets{false};
 
   /// defines if the mesh is centralized or distributed
   bool is_distributed{false};
 
   /// defines if the mesh is periodic
   bool is_periodic{false};
 
   /// Communicator on which mesh is distributed
   Communicator * communicator;
 
   /// Element synchronizer
   std::unique_ptr<ElementSynchronizer> element_synchronizer;
 
   /// Node synchronizer
   std::unique_ptr<NodeSynchronizer> node_synchronizer;
 
   /// Node synchronizer for periodic nodes
   std::unique_ptr<PeriodicNodeSynchronizer> periodic_node_synchronizer;
 
   using NodesToElements = std::vector<std::unique_ptr<std::set<Element>>>;
 
   /// class to update global data using external knowledge
   std::unique_ptr<MeshGlobalDataUpdater> global_data_updater;
 
   /// This info is stored to simplify the dynamic changes
   NodesToElements nodes_to_elements;
 
   /// periodicity local info
   std::unordered_map<Idx, Idx> periodic_slave_master;
   std::unordered_multimap<Idx, Idx> periodic_master_slave;
 };
 
 /// standard output stream operator
 inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 inline constexpr auto Mesh::getNbNodesPerElement(ElementType type) -> Int {
   return tuple_dispatch_with_default<AllElementTypes>(
       [](auto && enum_type) {
         constexpr ElementType type = aka::decay_v<decltype(enum_type)>;
         return ElementClass<type>::getNbNodesPerElement();
       },
       type, [](auto && /*type*/) { return 0; });
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto Mesh::getNbElement(ElementType type, GhostType ghost_type) const {
   try {
     const auto & conn = connectivities(type, ghost_type);
     return conn.size();
   } catch (...) {
     return 0;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 inline auto Mesh::getNbElement(const Int spatial_dimension,
                                GhostType ghost_type, ElementKind kind) const {
   AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == Int(-1),
                       "spatial_dimension is " << spatial_dimension
                                               << " and is greater than 3 !");
   Int nb_element = 0;
 
   for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) {
     nb_element += getNbElement(type, ghost_type);
   }
 
   return nb_element;
 }
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
 
 /* -------------------------------------------------------------------------- */
 /* Inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 #include "element_type_map_tmpl.hh"
 #include "mesh_inline_impl.hh"
 
 #endif /* AKANTU_MESH_HH_ */