diff --git a/cmake/libakantu/v3/printers.py b/cmake/libakantu/v3/printers.py
index de6356340..a8ca32bee 100755
--- a/cmake/libakantu/v3/printers.py
+++ b/cmake/libakantu/v3/printers.py
@@ -1,234 +1,339 @@
 #!/usr/bin/env python3
 # encoding: utf-8
 #
 # Inspired from boost's pretty printers from
 # Rüdiger Sonderfeld <ruediger@c-plusplus.de>
 # and from Pretty-printers for libstc++ from Free Software Foundation, Inc.
 #
-
 import gdb
 import re
-#import sys
-# import libstdcxx.v6.printers as std
+import sys
 
 __use_gdb_pp__ = True
 try:
     import gdb.printing
 except ImportError:
     __use_gdb_pp__ = False
 
 
 class AkantuPrinter(object):
     regex = None
 
     @classmethod
     def supports(cls, typename):
         #print('{0} ~= {1}'.format(typename, cls.regex), file=sys.stderr)
         return cls.regex.search(typename)
 
     @classmethod
     def get_basic_type(cls, value):
         """ Determines the type associated to a value"""
 
         _type = value.type
         # If it points to a reference, get the reference.
         if _type.code == gdb.TYPE_CODE_REF:
             _type = _type.target()
 
         # Get the unqualified type, stripped of typedefs.
         _type = _type.unqualified().strip_typedefs()
 
         return _type.tag
 
 if __use_gdb_pp__:
     __akantu_pretty_printers__ = \
         gdb.printing.RegexpCollectionPrettyPrinter("libakantu-v3")
 else:
     class AkantuPrettyPrinters(object):
         def __init__(self, name):
             super(AkantuPrettyPrinters, self).__init__()
             self.name = name
             self.printers = {}
 
         @property
         def enabled(self):
             return True
 
         def add_printer(self, name, regex, printer):
             self.printers[name] = printer
 
         def __call__(self, val):
             typename = AkantuPrinter.get_basic_type(val)
             if not typename:
                 return None
 
             for name, printer in self.printers.items():
                 if(printer.supports(typename)):
                     return printer
 
             return None
 
     __akantu_pretty_printers__ = AkantuPrettyPrinters("libakantu-v3")
 
 
 def register_pretty_printer(pretty_printer):
     "Registers a Pretty Printer"
 
     __akantu_pretty_printers__.add_printer(pretty_printer.name,
                                            pretty_printer.regex,
                                            pretty_printer)
 
     return pretty_printer
 
 @register_pretty_printer
 class AkaArrayPrinter(AkantuPrinter):
     """Pretty printer for akantu::Array<T>"""
     regex = re.compile('^akantu::Array<(.*?), (true|false)>$')
     name = 'akantu::Array'
 
     def __init__(self, value):
         self.typename = self.get_basic_type(value)
         self.value = value
         self.ptr = self.value['values']
         self.size = int(self.value['size_'])
         self.nb_component = int(self.value['nb_component'])
 
     def display_hint(self):
         return 'array'
 
     def to_string(self):
         m = self.regex.search(self.typename)
         return 'Array<{0}>({1}, {2}) stored at {3}'.format(
             m.group(1), self.size, self.nb_component, self.ptr)
 
     def children(self):
         _ptr = self.ptr
         for i in range(self.size):
             _values = ["{0}".format((_ptr + j).dereference())
                        for j in range(self.nb_component)]
             _ptr = _ptr + self.nb_component
             yield ('[{0}]'.format(i),
                    ('{0}' if self.nb_component == 1 else '[{0}]').format(
                        ', '.join(_values)))
 
-# @register_pretty_printer
-# class AkaArrayIteratorPrinter(AkantuPrinter):
-#     """Pretty printer for akantu::Array<T>"""
-#     regex = re.compile('^akantu::Array<(.*?), (true|false)>::internal_iterator<(.*?), (.*?), (false|true)>$')
-#     name = 'akantu::Array::iterator'
-
-#     def __init__(self, value):
-#         self.typename = self.get_basic_type(value)
-#         self.value = value
-#         self.ret = self.value['ret'].dereference()
-
-#     def to_string(self):
-#         m = self.regex.search(self.typename)
-#         return 'Array<{0}>::iterator<{3}>'.format(
-#             m.group(1), m.group(1))
-
-#     def children(self):
-#         yield ('[data]', self.ret)
+if sys.version_info[0] > 2:
+    ### Python 3 stuff
+    Iterator = object
+else:
+    ### Python 2 stuff
+    class Iterator:
+        """Compatibility mixin for iterators
+
+        Instead of writing next() methods for iterators, write
+        __next__() methods and use this mixin to make them work in
+        Python 2 as well as Python 3.
+
+        Idea stolen from the "six" documentation:
+        <http://pythonhosted.org/six/#six.Iterator>
+        """
+        def next(self):
+            return self.__next__()
+
+
+class RbtreeIterator(Iterator):
+    """
+    Turn an RB-tree-based container (std::map, std::set etc.) into
+    a Python iterable object.
+    """
+
+    def __init__(self, rbtree):
+        self.size = rbtree['_M_t']['_M_impl']['_M_node_count']
+        self.node = rbtree['_M_t']['_M_impl']['_M_header']['_M_left']
+        self.count = 0
+
+    def __iter__(self):
+        return self
+
+    def __len__(self):
+        return int (self.size)
+
+    def __next__(self):
+        if self.count == self.size:
+            raise StopIteration
+        result = self.node
+        self.count = self.count + 1
+        if self.count < self.size:
+            # Compute the next node.
+            node = self.node
+            if node.dereference()['_M_right']:
+                node = node.dereference()['_M_right']
+                while node.dereference()['_M_left']:
+                    node = node.dereference()['_M_left']
+            else:
+                parent = node.dereference()['_M_parent']
+                while node == parent.dereference()['_M_right']:
+                    node = parent
+                    parent = parent.dereference()['_M_parent']
+                if node.dereference()['_M_right'] != parent:
+                    node = parent
+                    self.node = node
+        return result
+
+def get_value_from_Rb_tree_node(node):
+    """Returns the value held in an _Rb_tree_node<_Val>"""
+    try:
+        member = node.type.fields()[1].name
+        if member == '_M_value_field':
+            # C++03 implementation, node contains the value as a member
+            return node['_M_value_field']
+        elif member == '_M_storage':
+            # C++11 implementation, node stores value in __aligned_membuf
+            valtype = node.type.template_argument(0)
+            return get_value_from_aligned_membuf(node['_M_storage'], valtype)
+    except:
+        pass
+    raise ValueError("Unsupported implementation for %s" % str(node.type))
+
+@register_pretty_printer
+class AkaElementTypeMapArrayPrinter(AkantuPrinter):
+    """Pretty printer for akantu::ElementTypeMap<Array<T>>"""
+    regex = re.compile('^akantu::ElementTypeMap<akantu::Array<(.*?), (true|false)>\*, akantu::(.*?)>$')
+    name = 'akantu::ElementTypeMapArray'
+
+    # Turn an RbtreeIterator into a pretty-print iterator.
+    class _iter(Iterator):
+        def __init__(self, rbiter, type):
+            self.rbiter = rbiter
+            self.count = 0
+            self.type = type
+
+        def __iter__(self):
+            return self
+
+        def __next__(self):
+            if self.count % 2 == 0:
+                n = next(self.rbiter)
+                n = n.cast(self.type).dereference()
+                n = get_value_from_Rb_tree_node(n)
+                self.pair = n
+                item = n['first']
+            else:
+                item = self.pair['second']
+            result = ('[{0}]'.format(self.count), item.dereference())
+            self.count = self.count + 1
+            return result
+
+    def __init__(self, value):
+        self.typename = self.get_basic_type(value)
+        self.value = value
+        self.data = self.value['data']
+        self.ghost_data = self.value['ghost_data']
+
+    def to_string(self):
+        m = self.regex.search(self.typename)
+        return 'ElementTypMapArray<{0}> with {1} _not_ghost and {2} _ghost'.format(
+            m.group(1), len(RbtreeIterator(self.data)), len(RbtreeIterator(self.ghost_data)))
+
+    def children(self):
+        m = self.regex.search(self.typename)
+        try:
+            _type = gdb.lookup_type("akantu::Array<{0}, {1}>".format(
+                m.group(1), m.group(2))).strip_typedefs().pointer()
+            yield ('[_not_ghost]', self._iter(RbtreeIterator(self.data), _type))
+            yield ('[_ghost]', self._iter(RbtreeIterator(self.ghost_data), _type))
+        except RuntimeError:
+            pass
+
+    def display_hint(self):
+        return 'map'
 
 
 # @register_pretty_printer
 class AkaTensorPrinter(AkantuPrinter):
     """Pretty printer for akantu::Tensor<T>"""
     regex = re.compile('^akantu::Tensor<(.*), +(.*), +(.*)>$')
     name = 'akantu::Tensor'
 
     value = None
     typename = ""
     ptr = None
     dims = []
     ndims = 0
 
     def pretty_print(self):
         def ij2str(i, j, m):
             return "{0}".format((self.ptr+m*j + i).dereference())
 
         def line(i, m, n):
             return "[{0}]".format(", ".join((ij2str(i, j, m) for j in
                                              range(n))))
 
         m = int(self.dims[0])
         if (self.ndims == 1):
             n = 1
         else:
             n = int(self.dims[1])
 
         return "[{0}]".format(", ".join(line(i, m, n) for i in range(m)))
 
     def __init__(self, value):
         self.typename = self.get_basic_type(value)
         self.value = value
         self.ptr = self.value['values']
         self.dims = self.value['n']
 
     def children(self):
         yield ('values', self.pretty_print())
         yield ('wrapped', self.value['wrapped'])
 
 @register_pretty_printer
 class AkaVectorPrinter(AkaTensorPrinter):
     """Pretty printer for akantu::Vector<T>"""
     regex = re.compile('^akantu::Vector<(.*)>$')
     name = 'akantu::Vector'
     n = 0
     ptr = 0x0
 
     def __init__(self, value):
         super(AkaVectorPrinter, self).__init__(value)
         self.ndims = 1
 
     def to_string(self):
         m = self.regex.search(self.typename)
         return 'Vector<{0}>({1}) [{2}]'.format(m.group(1), int(self.dims[0]),
                                                str(self.ptr))
 
 
 @register_pretty_printer
 class AkaMatrixPrinter(AkaTensorPrinter):
     """Pretty printer for akantu::Matrix<T>"""
     regex = re.compile('^akantu::Matrix<(.*)>$')
     name = 'akantu::Matrix'
 
     def __init__(self, value):
         super(AkaMatrixPrinter, self).__init__(value)
         self.ndims = 2
 
     def to_string(self):
         m = self.regex.search(self.typename)
         return 'Matrix<%s>(%d, %d) [%s]' % (m.group(1), int(self.dims[0]),
                                             int(self.dims[1]),
                                             str(self.ptr))
 
 
 @register_pretty_printer
 class AkaElementPrinter(AkantuPrinter):
     """Pretty printer for akantu::Element"""
     regex = re.compile('^akantu::Element$')
     name = 'akantu::Element'
 
     def __init__(self, value):
         self.typename = self.get_basic_type(value)
         self.value = value
 
         self.element = self.value['element']
         self.eltype = self.value['type']
         self.ghost_type = self.value['ghost_type']
 
     def to_string(self):
         return 'Element({0}, {1}, {2})'.format(self.element, self.eltype,
                                                self.ghost_type)
 
 
 def register_akantu_printers(obj):
     "Register Akantu Pretty Printers."
 
     if __use_gdb_pp__:
         gdb.printing.register_pretty_printer(obj, __akantu_pretty_printers__)
     else:
         if obj is None:
             obj = gdb
         obj.pretty_printers.append(__akantu_pretty_printers__)
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index a5c46d45c..ac1b92fd9 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1260 +1,1263 @@
 /**
  * @file   aka_array_tmpl.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Jul 15 2010
  * @date last modification: Fri Jan 22 2016
  *
  * @brief  Inline functions of the classes Array<T> and ArrayBase
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions Array<T>                                                 */
 /* -------------------------------------------------------------------------- */
 #include "aka_array.hh"
 /* -------------------------------------------------------------------------- */
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
 #define __AKANTU_AKA_ARRAY_TMPL_HH__
 
 namespace akantu {
 
 namespace debug {
   struct ArrayException : public Exception {};
 } // namespace debug
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline T & Array<T, is_scal>::operator()(UInt i, UInt j) {
   AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
                       "The value at position ["
                           << i << "," << j << "] is out of range in array \""
                           << id << "\"");
   return values[i * nb_component + j];
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline const T & Array<T, is_scal>::operator()(UInt i, UInt j) const {
   AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
                       "The value at position ["
                           << i << "," << j << "] is out of range in array \""
                           << id << "\"");
   return values[i * nb_component + j];
 }
 
 template <class T, bool is_scal>
 inline T & Array<T, is_scal>::operator[](UInt i) {
   AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
                       "The value at position ["
                           << i << "] is out of range in array \"" << id
                           << "\"");
   return values[i];
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 inline const T & Array<T, is_scal>::operator[](UInt i) const {
   AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
   AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
                       "The value at position ["
                           << i << "] is out of range in array \"" << id
                           << "\"");
   return values[i];
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * append a tuple to the array with the value value for all components
  * @param value the new last tuple or the array will contain nb_component copies
  * of value
  */
 template <class T, bool is_scal>
 inline void Array<T, is_scal>::push_back(const T & value) {
   resizeUnitialized(size_ + 1, true, value);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * append a tuple to the array
  * @param new_elem a C-array containing the values to be copied to the end of
  * the array */
 // template <class T, bool is_scal>
 // inline void Array<T, is_scal>::push_back(const T new_elem[]) {
 //   UInt pos = size_;
 
 //   resizeUnitialized(size_ + 1, false);
 
 //   T * tmp = values + nb_component * pos;
 //   std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
 // }
 
 /* -------------------------------------------------------------------------- */
 #ifndef SWIG
 /**
  * append a matrix or a vector to the array
  * @param new_elem a reference to a Matrix<T> or Vector<T> */
 template <class T, bool is_scal>
 template <template <typename> class C, typename>
 inline void Array<T, is_scal>::push_back(const C<T> & new_elem) {
   AKANTU_DEBUG_ASSERT(
       nb_component == new_elem.size(),
       "The vector("
           << new_elem.size()
           << ") as not a size compatible with the Array (nb_component="
           << nb_component << ").");
   UInt pos = size_;
   resizeUnitialized(size_ + 1, false);
 
   T * tmp = values + nb_component * pos;
   std::uninitialized_copy(new_elem.storage(), new_elem.storage() + nb_component,
                           tmp);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * append a tuple to the array
  * @param it an iterator to the tuple to be copied to the end of the array */
 template <class T, bool is_scal>
 template <class Ret>
 inline void
 Array<T, is_scal>::push_back(const Array<T, is_scal>::iterator<Ret> & it) {
   UInt pos = size_;
 
   resizeUnitialized(size_ + 1, false);
 
   T * tmp = values + nb_component * pos;
   T * new_elem = it.data();
   std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
 }
 #endif
 /* -------------------------------------------------------------------------- */
 /**
  * erase an element. If the erased element is not the last of the array, the
  * last element is moved into the hole in order to maintain contiguity. This
  * may invalidate existing iterators (For instance an iterator obtained by
  * Array::end() is no longer correct) and will change the order of the
  * elements.
  * @param i index of element to erase
  */
 template <class T, bool is_scal> inline void Array<T, is_scal>::erase(UInt i) {
   AKANTU_DEBUG_IN();
   AKANTU_DEBUG_ASSERT((size_ > 0), "The array is empty");
   AKANTU_DEBUG_ASSERT((i < size_),
                       "The element at position [" << i << "] is out of range ("
                                                   << i << ">=" << size_ << ")");
 
   if (i != (size_ - 1)) {
     for (UInt j = 0; j < nb_component; ++j) {
       values[i * nb_component + j] = values[(size_ - 1) * nb_component + j];
     }
   }
 
   resize(size_ - 1);
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Subtract another array entry by entry from this array in place. Both arrays
  * must
  * have the same size and nb_component. If the arrays have different shapes,
  * code compiled in debug mode will throw an expeption and optimised code
  * will behave in an unpredicted manner
  * @param other array to subtract from this
  * @return reference to modified this
  */
 template <class T, bool is_scal>
 Array<T, is_scal> & Array<T, is_scal>::
 operator-=(const Array<T, is_scal> & vect) {
   AKANTU_DEBUG_ASSERT((size_ == vect.size_) &&
                           (nb_component == vect.nb_component),
                       "The too array don't have the same sizes");
 
   T * a = values;
   T * b = vect.storage();
   for (UInt i = 0; i < size_ * nb_component; ++i) {
     *a -= *b;
     ++a;
     ++b;
   }
 
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Add another array entry by entry to this array in place. Both arrays must
  * have the same size and nb_component. If the arrays have different shapes,
  * code compiled in debug mode will throw an expeption and optimised code
  * will behave in an unpredicted manner
  * @param other array to add to this
  * @return reference to modified this
  */
 template <class T, bool is_scal>
 Array<T, is_scal> & Array<T, is_scal>::
 operator+=(const Array<T, is_scal> & vect) {
   AKANTU_DEBUG_ASSERT((size_ == vect.size) &&
                           (nb_component == vect.nb_component),
                       "The too array don't have the same sizes");
 
   T * a = values;
   T * b = vect.storage();
   for (UInt i = 0; i < size_ * nb_component; ++i) {
     *a++ += *b++;
   }
 
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Multiply all entries of this array by a scalar in place
  * @param alpha scalar multiplicant
  * @return reference to modified this
  */
 template <class T, bool is_scal>
 Array<T, is_scal> & Array<T, is_scal>::operator*=(const T & alpha) {
   T * a = values;
   for (UInt i = 0; i < size_ * nb_component; ++i) {
     *a++ *= alpha;
   }
 
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Compare this array element by element to another.
  * @param other array to compare to
  * @return true it all element are equal and arrays have the same shape, else
  * false
  */
 template <class T, bool is_scal>
 bool Array<T, is_scal>::operator==(const Array<T, is_scal> & array) const {
   bool equal = nb_component == array.nb_component && size_ == array.size_ &&
                id == array.id;
   if (!equal)
     return false;
 
   if (values == array.storage())
     return true;
   else
     return std::equal(values, values + size_ * nb_component, array.storage());
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & array) const {
   return !operator==(array);
 }
 
 /* -------------------------------------------------------------------------- */
 #ifndef SWIG
 /**
  * set all tuples of the array to a given vector or matrix
  * @param vm Matrix or Vector to fill the array with
  */
 template <class T, bool is_scal>
 template <template <typename> class C, typename>
 inline void Array<T, is_scal>::set(const C<T> & vm) {
   AKANTU_DEBUG_ASSERT(
       nb_component == vm.size(),
       "The size of the object does not match the number of components");
   for (T * it = values; it < values + nb_component * size_;
        it += nb_component) {
     std::copy(vm.storage(), vm.storage() + nb_component, it);
   }
 }
 #endif
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 void Array<T, is_scal>::append(const Array<T> & other) {
   AKANTU_DEBUG_ASSERT(
       nb_component == other.nb_component,
       "Cannot append an array with a different number of component");
   UInt old_size = this->size_;
   this->resizeUnitialized(this->size_ + other.size(), false);
 
   T * tmp = values + nb_component * old_size;
   std::uninitialized_copy(other.storage(),
                           other.storage() + other.size() * nb_component, tmp);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Functions Array<T, is_scal>                                                */
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
     : ArrayBase(id), values(nullptr) {
   AKANTU_DEBUG_IN();
   allocate(size, nb_component);
 
   if (!is_scal) {
     T val = T();
     std::uninitialized_fill(values, values + size * nb_component, val);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(UInt size, UInt nb_component, const T def_values[],
                          const ID & id)
     : ArrayBase(id), values(NULL) {
   AKANTU_DEBUG_IN();
   allocate(size, nb_component);
 
   T * tmp = values;
 
   for (UInt i = 0; i < size; ++i) {
     tmp = values + nb_component * i;
     std::uninitialized_copy(def_values, def_values + nb_component, tmp);
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(UInt size, UInt nb_component, const T & value,
                          const ID & id)
     : ArrayBase(id), values(nullptr) {
   AKANTU_DEBUG_IN();
   allocate(size, nb_component);
 
   std::uninitialized_fill_n(values, size * nb_component, value);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(const Array<T, is_scal> & vect, bool deep,
                          const ID & id)
     : ArrayBase(vect) {
   AKANTU_DEBUG_IN();
   this->id = (id == "") ? vect.id : id;
 
   if (deep) {
     allocate(vect.size_, vect.nb_component);
     T * tmp = values;
     std::uninitialized_copy(vect.storage(),
                             vect.storage() + size_ * nb_component, tmp);
   } else {
     this->values = vect.storage();
     this->size_ = vect.size_;
     this->nb_component = vect.nb_component;
     this->allocated_size = vect.allocated_size;
     this->size_of_type = vect.size_of_type;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 #ifndef SWIG
 template <class T, bool is_scal>
 Array<T, is_scal>::Array(const std::vector<T> & vect) {
   AKANTU_DEBUG_IN();
   this->id = "";
 
   allocate(vect.size(), 1);
   T * tmp = values;
   std::uninitialized_copy(&(vect[0]), &(vect[size_ - 1]), tmp);
 
   AKANTU_DEBUG_OUT();
 }
 #endif
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal> Array<T, is_scal>::~Array() {
   AKANTU_DEBUG_IN();
   AKANTU_DEBUG(dblAccessory,
                "Freeing " << printMemorySize<T>(allocated_size * nb_component)
                           << " (" << id << ")");
 
   if (values) {
     if (!is_scal)
       for (UInt i = 0; i < size_ * nb_component; ++i) {
         T * obj = values + i;
         obj->~T();
       }
     free(values);
   }
   size_ = allocated_size = 0;
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * perform the allocation for the constructors
  * @param size is the size of the array
  * @param nb_component is the number of component of the array
  */
 template <class T, bool is_scal>
 void Array<T, is_scal>::allocate(UInt size, UInt nb_component) {
   AKANTU_DEBUG_IN();
   if (size == 0) {
     values = nullptr;
   } else {
     values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
     AKANTU_DEBUG_ASSERT(values != nullptr,
                         "Cannot allocate "
                             << printMemorySize<T>(size * nb_component) << " ("
                             << id << ")");
   }
 
   if (values == nullptr) {
     this->size_ = this->allocated_size = 0;
   } else {
     AKANTU_DEBUG(dblAccessory,
                  "Allocated " << printMemorySize<T>(size * nb_component) << " ("
                               << id << ")");
     this->size_ = this->allocated_size = size;
   }
 
   this->size_of_type = sizeof(T);
   this->nb_component = nb_component;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 void Array<T, is_scal>::reserve(UInt new_size) {
   UInt tmp_size = this->size_;
   resizeUnitialized(new_size, false);
   this->size_ = tmp_size;
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * change the size of the array and allocate or free memory if needed. If the
  * size increases, the new tuples are filled with zeros
  * @param new_size new number of tuples contained in the array */
 template <class T, bool is_scal> void Array<T, is_scal>::resize(UInt new_size) {
   resizeUnitialized(new_size, !is_scal);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * change the size of the array and allocate or free memory if needed. If the
  * size increases, the new tuples are filled with zeros
  * @param new_size new number of tuples contained in the array */
 template <class T, bool is_scal>
 void Array<T, is_scal>::resize(UInt new_size, const T & val) {
   this->resizeUnitialized(new_size, true, val);
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * change the size of the array and allocate or free memory if needed.
  * @param new_size new number of tuples contained in the array */
 template <class T, bool is_scal>
 void Array<T, is_scal>::resizeUnitialized(UInt new_size, bool fill,
                                           const T & val) {
   //  AKANTU_DEBUG_IN();
   // free some memory
   if (new_size <= allocated_size) {
     if (!is_scal) {
       T * old_values = values;
       if (new_size < size_) {
         for (UInt i = new_size * nb_component; i < size_ * nb_component; ++i) {
           T * obj = old_values + i;
           obj->~T();
         }
       }
     }
 
     if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
       AKANTU_DEBUG(dblAccessory,
                    "Freeing " << printMemorySize<T>((allocated_size - size_) *
                                                     nb_component)
                               << " (" << id << ")");
 
       // Normally there are no allocation problem when reducing an array
       if (new_size == 0) {
         free(values);
         values = nullptr;
       } else {
         auto * tmp_ptr = static_cast<T *>(
             realloc(values, new_size * nb_component * sizeof(T)));
 
         if (tmp_ptr == nullptr) {
           AKANTU_EXCEPTION("Cannot free data ("
                            << id << ")"
                            << " [current allocated size : " << allocated_size
                            << " | "
                            << "requested size : " << new_size << "]");
         }
         values = tmp_ptr;
       }
       allocated_size = new_size;
     }
   } else {
     // allocate more memory
     UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
                              ? allocated_size + AKANTU_MIN_ALLOCATION
                              : new_size;
 
     auto * tmp_ptr = static_cast<T *>(
         realloc(values, size_to_alloc * nb_component * sizeof(T)));
     AKANTU_DEBUG_ASSERT(
         tmp_ptr != nullptr,
         "Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
     if (tmp_ptr == nullptr) {
       AKANTU_DEBUG_ERROR("Cannot allocate more data ("
                          << id << ")"
                          << " [current allocated size : " << allocated_size
                          << " | "
                          << "requested size : " << new_size << "]");
     }
 
     AKANTU_DEBUG(dblAccessory,
                  "Allocating " << printMemorySize<T>(
                      (size_to_alloc - allocated_size) * nb_component));
 
     allocated_size = size_to_alloc;
     values = tmp_ptr;
   }
 
   if (fill && this->size_ < new_size) {
     std::uninitialized_fill(values + size_ * nb_component,
                             values + new_size * nb_component, val);
   }
 
   size_ = new_size;
   //  AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * change the number of components by interlacing data
  * @param multiplicator number of interlaced components add
  * @param block_size blocks of data in the array
  * Examaple for block_size = 2, multiplicator = 2
  * array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn */
 template <class T, bool is_scal>
 void Array<T, is_scal>::extendComponentsInterlaced(UInt multiplicator,
                                                    UInt block_size) {
   AKANTU_DEBUG_IN();
 
   if (multiplicator == 1)
     return;
 
   AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator");
   AKANTU_DEBUG_ASSERT(nb_component % block_size == 0,
                       "stride must divide actual number of components");
 
   values = static_cast<T *>(
       realloc(values, nb_component * multiplicator * size_ * sizeof(T)));
 
   UInt new_component = nb_component / block_size * multiplicator;
 
   for (UInt i = 0, k = size_ - 1; i < size_; ++i, --k) {
     for (UInt j = 0; j < new_component; ++j) {
       UInt m = new_component - j - 1;
       UInt n = m / multiplicator;
       for (UInt l = 0, p = block_size - 1; l < block_size; ++l, --p) {
         values[k * nb_component * multiplicator + m * block_size + p] =
             values[k * nb_component + n * block_size + p];
       }
     }
   }
 
   nb_component = nb_component * multiplicator;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * search elem in the array, return  the position of the first occurrence or
  * -1 if not found
  *  @param elem the element to look for
  *  @return index of the first occurrence of elem or -1 if elem is not present
  */
 template <class T, bool is_scal>
 UInt Array<T, is_scal>::find(const T & elem) const {
   AKANTU_DEBUG_IN();
 
   auto begin = this->begin();
   auto end = this->end();
   auto it = std::find(begin, end, elem);
 
   AKANTU_DEBUG_OUT();
   return (it != end) ? it - begin : UInt(-1);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal> UInt Array<T, is_scal>::find(T elem[]) const {
   AKANTU_DEBUG_IN();
   T * it = values;
   UInt i = 0;
   for (; i < size_; ++i) {
     if (*it == elem[0]) {
       T * cit = it;
       UInt c = 0;
       for (; (c < nb_component) && (*cit == elem[c]); ++c, ++cit)
         ;
       if (c == nb_component) {
         AKANTU_DEBUG_OUT();
         return i;
       }
     }
     it += nb_component;
   }
   return UInt(-1);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <template <typename> class C, typename>
 inline UInt Array<T, is_scal>::find(const C<T> & elem) {
   AKANTU_DEBUG_ASSERT(elem.size() == nb_component,
                       "Cannot find an element with a wrong size ("
                           << elem.size() << ") != " << nb_component);
   return this->find(elem.storage());
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * copy the content of another array. This overwrites the current content.
  * @param other Array to copy into this array. It has to have the same
  * nb_component as this. If compiled in debug mode, an incorrect other will
  * result in an exception being thrown. Optimised code may result in
  * unpredicted behaviour.
  */
 template <class T, bool is_scal>
 void Array<T, is_scal>::copy(const Array<T, is_scal> & vect,
                              bool no_sanity_check) {
   AKANTU_DEBUG_IN();
 
   if (!no_sanity_check)
     if (vect.nb_component != nb_component)
       AKANTU_DEBUG_ERROR(
           "The two arrays do not have the same number of components");
 
   resize((vect.size_ * vect.nb_component) / nb_component);
 
   T * tmp = values;
   std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
                           tmp);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <bool is_scal> class ArrayPrintHelper {
 public:
   template <typename T>
   static void print_content(const Array<T> & vect, std::ostream & stream,
                             int indent) {
     if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
       std::string space;
       for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
         ;
 
       stream << space << " + values         : {";
       for (UInt i = 0; i < vect.size(); ++i) {
         stream << "{";
         for (UInt j = 0; j < vect.getNbComponent(); ++j) {
           stream << vect(i, j);
           if (j != vect.getNbComponent() - 1)
             stream << ", ";
         }
         stream << "}";
         if (i != vect.size() - 1)
           stream << ", ";
       }
       stream << "}" << std::endl;
     }
   }
 };
 
 template <> class ArrayPrintHelper<false> {
 public:
   template <typename T>
   static void print_content(__attribute__((unused)) const Array<T> & vect,
                             __attribute__((unused)) std::ostream & stream,
                             __attribute__((unused)) int indent) {}
 };
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
   std::string space;
   for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
     ;
 
   std::streamsize prec = stream.precision();
   std::ios_base::fmtflags ff = stream.flags();
 
   stream.setf(std::ios_base::showbase);
   stream.precision(2);
 
   stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
          << std::endl;
   stream << space << " + id             : " << this->id << std::endl;
   stream << space << " + size           : " << this->size_ << std::endl;
   stream << space << " + nb_component   : " << this->nb_component << std::endl;
   stream << space << " + allocated size : " << this->allocated_size
          << std::endl;
   stream << space << " + memory size    : "
          << printMemorySize<T>(allocated_size * nb_component) << std::endl;
   if (!AKANTU_DEBUG_LEVEL_IS_TEST())
     stream << space << " + address        : " << std::hex << this->values
            << std::dec << std::endl;
 
   stream.precision(prec);
   stream.flags(ff);
 
   ArrayPrintHelper<is_scal>::print_content(*this, stream, indent);
 
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions ArrayBase                                                */
 /* -------------------------------------------------------------------------- */
 
 inline UInt ArrayBase::getMemorySize() const {
   return allocated_size * nb_component * size_of_type;
 }
 
 inline void ArrayBase::empty() { size_ = 0; }
 
 #ifndef SWIG
 /* -------------------------------------------------------------------------- */
 /* Iterators                                                                  */
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <class R, class daughter, class IR, bool is_tensor>
 class Array<T, is_scal>::iterator_internal {
 public:
   using value_type = R;
   using pointer = R *;
   using reference = R &;
   using proxy = typename R::proxy;
   using const_proxy = const typename R::proxy;
   using const_reference = const R &;
   using internal_value_type = IR;
   using internal_pointer = IR *;
   using difference_type = std::ptrdiff_t;
   using iterator_category = std::random_access_iterator_tag;
 
 public:
   iterator_internal() = default;
 
   iterator_internal(pointer_type data, UInt _offset)
       : _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
     AKANTU_DEBUG_ERROR(
         "The constructor should never be called it is just an ugly trick...");
   }
 
   iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
       : _offset(wrapped->size()), initial(wrapped->storage()),
         ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
 
   iterator_internal(const iterator_internal & it) {
     if (this != &it) {
       this->_offset = it._offset;
       this->initial = it.initial;
       this->ret_ptr = it.ret_ptr;
       this->ret = std::make_unique<internal_value_type>(*it.ret, false);
     }
   }
 
   iterator_internal(iterator_internal && it) = default;
 
   virtual ~iterator_internal() = default;
 
   inline iterator_internal & operator=(const iterator_internal & it) {
     if (this != &it) {
       this->_offset = it._offset;
       this->initial = it.initial;
       this->ret_ptr = it.ret_ptr;
       if (this->ret)
         this->ret->shallowCopy(*it.ret);
       else
         this->ret = std::make_unique<internal_value_type>(*it.ret, false);
     }
     return *this;
   }
 
   UInt getCurrentIndex() {
     return (this->ret_ptr - this->initial) / this->_offset;
   };
 
   inline reference operator*() {
     ret->values = ret_ptr;
     return *ret;
   };
   inline const_reference operator*() const {
     ret->values = ret_ptr;
     return *ret;
   };
   inline pointer operator->() {
     ret->values = ret_ptr;
     return ret.get();
   };
   inline daughter & operator++() {
     ret_ptr += _offset;
     return static_cast<daughter &>(*this);
   };
   inline daughter & operator--() {
     ret_ptr -= _offset;
     return static_cast<daughter &>(*this);
   };
 
   inline daughter & operator+=(const UInt n) {
     ret_ptr += _offset * n;
     return static_cast<daughter &>(*this);
   }
   inline daughter & operator-=(const UInt n) {
     ret_ptr -= _offset * n;
     return static_cast<daughter &>(*this);
   }
 
   inline proxy operator[](const UInt n) {
     ret->values = ret_ptr + n * _offset;
     return proxy(*ret);
   }
   inline const_proxy operator[](const UInt n) const {
     ret->values = ret_ptr + n * _offset;
     return const_proxy(*ret);
   }
 
   inline bool operator==(const iterator_internal & other) const {
     return this->ret_ptr == other.ret_ptr;
   }
   inline bool operator!=(const iterator_internal & other) const {
     return this->ret_ptr != other.ret_ptr;
   }
   inline bool operator<(const iterator_internal & other) const {
     return this->ret_ptr < other.ret_ptr;
   }
   inline bool operator<=(const iterator_internal & other) const {
     return this->ret_ptr <= other.ret_ptr;
   }
   inline bool operator>(const iterator_internal & other) const {
     return this->ret_ptr > other.ret_ptr;
   }
   inline bool operator>=(const iterator_internal & other) const {
     return this->ret_ptr >= other.ret_ptr;
   }
 
   inline daughter operator+(difference_type n) {
     daughter tmp(static_cast<daughter &>(*this));
     tmp += n;
     return tmp;
   }
   inline daughter operator-(difference_type n) {
     daughter tmp(static_cast<daughter &>(*this));
     tmp -= n;
     return tmp;
   }
 
   inline difference_type operator-(const iterator_internal & b) {
     return (this->ret_ptr - b.ret_ptr) / _offset;
   }
 
   inline pointer_type data() const { return ret_ptr; }
   inline difference_type offset() const { return _offset; }
 
 protected:
   UInt _offset{0};
   pointer_type initial{nullptr};
   std::unique_ptr<internal_value_type> ret{nullptr};
   pointer_type ret_ptr{nullptr};
 };
 
 /* -------------------------------------------------------------------------- */
 /**
  * Specialization for scalar types
  */
 template <class T, bool is_scal>
 template <class R, class daughter, class IR>
 class Array<T, is_scal>::iterator_internal<R, daughter, IR, false> {
 public:
   using value_type = R;
   using pointer = R *;
   using reference = R &;
   using const_reference = const R &;
   using internal_value_type = IR;
   using internal_pointer = IR *;
   using difference_type = std::ptrdiff_t;
   using iterator_category = std::random_access_iterator_tag;
 
 public:
   iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
   iterator_internal(const iterator_internal & it) = default;
   iterator_internal(iterator_internal && it) = default;
 
   virtual ~iterator_internal() = default;
 
   inline iterator_internal & operator=(const iterator_internal & it) = default;
 
   UInt getCurrentIndex() { return (this->ret - this->initial); };
 
   inline reference operator*() { return *ret; };
   inline const_reference operator*() const { return *ret; };
   inline pointer operator->() { return ret; };
   inline daughter & operator++() {
     ++ret;
     return static_cast<daughter &>(*this);
   };
   inline daughter & operator--() {
     --ret;
     return static_cast<daughter &>(*this);
   };
 
   inline daughter & operator+=(const UInt n) {
     ret += n;
     return static_cast<daughter &>(*this);
   }
   inline daughter & operator-=(const UInt n) {
     ret -= n;
     return static_cast<daughter &>(*this);
   }
 
   inline reference operator[](const UInt n) { return ret[n]; }
 
   inline bool operator==(const iterator_internal & other) const {
     return ret == other.ret;
   }
   inline bool operator!=(const iterator_internal & other) const {
     return ret != other.ret;
   }
   inline bool operator<(const iterator_internal & other) const {
     return ret < other.ret;
   }
   inline bool operator<=(const iterator_internal & other) const {
     return ret <= other.ret;
   }
   inline bool operator>(const iterator_internal & other) const {
     return ret > other.ret;
   }
   inline bool operator>=(const iterator_internal & other) const {
     return ret >= other.ret;
   }
 
   inline daughter operator-(difference_type n) { return daughter(ret - n); }
   inline daughter operator+(difference_type n) { return daughter(ret + n); }
 
   inline difference_type operator-(const iterator_internal & b) {
     return ret - b.ret;
   }
 
   inline pointer data() const { return ret; }
 
 protected:
   pointer ret{nullptr};
   pointer initial{nullptr};
 };
 
 /* -------------------------------------------------------------------------- */
 /* Begin/End functions implementation                                         */
 /* -------------------------------------------------------------------------- */
 namespace detail {
   template <class Tuple, size_t... Is>
   constexpr auto take_front_impl(Tuple && t, std::index_sequence<Is...>) {
     return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
   }
 
   template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
     return take_front_impl(std::forward<Tuple>(t),
                            std::make_index_sequence<N>{});
   }
 
   template <typename... V>
   constexpr auto product_all(V &&... v) ->
       typename std::common_type<V...>::type {
     typename std::common_type<V...>::type result = 1;
     (void)std::initializer_list<int>{(result *= v, 0)...};
     return result;
   }
 
   template <typename... T> std::string to_string_all(T &&... t) {
     if (sizeof...(T) == 0)
       return "";
 
     std::stringstream ss;
     bool noComma = true;
     ss << "(";
     (void)std::initializer_list<bool>{
         (ss << (noComma ? "" : ", ") << t, noComma = false)...};
     ss << ")";
     return ss.str();
   }
 
   template <std::size_t N> struct InstantiationHelper {
     template <typename type, typename T, typename... Ns>
     static auto instantiate(T && data, Ns... ns) {
       return std::make_unique<type>(data, ns...);
     }
   };
 
   template <> struct InstantiationHelper<0> {
     template <typename type, typename T> static auto instantiate(T && data) {
       return data;
     }
   };
 
   template <typename Arr, typename T, typename... Ns>
   decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
     using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
     using array_type = std::decay_t<Arr>;
     using iterator =
-        std::conditional_t<std::is_const<Arr>::value,
+        std::conditional_t<std::is_const<std::remove_reference_t<Arr>>::value,
                            typename array_type::template const_iterator<type>,
                            typename array_type::template iterator<type>>;
+    static_assert(sizeof...(Ns), "You should provide a least one size");
+
     if (array.getNbComponent() * array.size() !=
         product_all(std::forward<Ns>(ns)...)) {
       AKANTU_CUSTOM_EXCEPTION_INFO(
           debug::ArrayException(),
-          "The iterator on Array "
+          "The iterator on " << debug::demangle(typeid(Arr).name())
               << to_string_all(array.size(), array.getNbComponent())
               << "is not compatible with the type "
               << debug::demangle(typeid(type).name()) << to_string_all(ns...));
     }
 
     auto && wrapped = aka::apply(
         [&](auto... n) {
           return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
               data, n...);
         },
         take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
 
     return iterator(std::move(wrapped));
   }
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
   return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
   return detail::get_iterator(*this, values + nb_component * size_,
                               std::forward<Ns>(ns)..., size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
   return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
   return detail::get_iterator(*this, values + nb_component * size_,
                               std::forward<Ns>(ns)..., size_);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
   return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
   return detail::get_iterator(
       *this, values + detail::product_all(std::forward<Ns>(ns)...),
       std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
   return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
 }
 
 template <class T, bool is_scal>
 template <typename... Ns>
 inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
   return detail::get_iterator(
       *this, values + detail::product_all(std::forward<Ns>(ns)...),
       std::forward<Ns>(ns)...);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Views                                                                      */
 /* -------------------------------------------------------------------------- */
 namespace detail {
   template <typename Array, typename... Ns> class ArrayView {
     using tuple = std::tuple<Ns...>;
 
   public:
     ArrayView(Array && array, Ns... ns)
         : array(std::forward<Array>(array)), sizes(std::move(ns)...){};
 
     decltype(auto) begin() {
       return aka::apply(
           [&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
     }
 
     decltype(auto) end() {
       return aka::apply(
           [&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
     }
 
     decltype(auto) size() {
       return std::get<std::tuple_size<tuple>::value - 1>(sizes);
     }
 
     decltype(auto) dims() {
       return std::tuple_size<tuple>::value - 1;
     }
   private:
     Array array;
     tuple sizes;
   };
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <typename Array, typename... Ns>
 decltype(auto) make_view(Array && array, Ns... ns) {
+  static_assert(sizeof...(Ns), "You should provide a least one dimension");
   auto size = std::forward<decltype(array)>(array).size() *
               std::forward<decltype(array)>(array).getNbComponent() /
               detail::product_all(ns...);
 
   return detail::ArrayView<Array, std::common_type_t<size_t, Ns>...,
                            std::common_type_t<size_t, decltype(size)>>(
       std::forward<Array>(array), std::move(ns)..., size);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename R>
 class Array<T, is_scal>::const_iterator
     : public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
                                R> {
 public:
   using parent = iterator_internal<const R, const_iterator, R>;
   using value_type = typename parent::value_type;
   using pointer = typename parent::pointer;
   using reference = typename parent::reference;
   using difference_type = typename parent::difference_type;
   using iterator_category = typename parent::iterator_category;
 
 public:
   const_iterator() : parent(){};
   // const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
   // const_iterator(pointer warped) : parent(warped) {}
   // const_iterator(const parent & it) : parent(it) {}
 
   const_iterator(const const_iterator & it) = default;
   const_iterator(const_iterator && it) = default;
 
   template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
   const_iterator(P * data) : parent(data) {}
 
   template <typename UP_P,
             typename =
                 std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
   const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
 
   const_iterator & operator=(const const_iterator & it) = default;
 };
 
 template <class T, class R, bool is_tensor_ = is_tensor<R>::value>
 struct ConstConverterIteratorHelper {
   using const_iterator = typename Array<T>::template const_iterator<R>;
   using iterator = typename Array<T>::template iterator<R>;
 
   static inline const_iterator convert(const iterator & it) {
     return const_iterator(std::unique_ptr<R>(new R(*it, false)));
   }
 };
 
 template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
   using const_iterator = typename Array<T>::template const_iterator<R>;
   using iterator = typename Array<T>::template iterator<R>;
   static inline const_iterator convert(const iterator & it) {
     return const_iterator(it.data());
   }
 };
 
 template <class T, bool is_scal>
 template <typename R>
 class Array<T, is_scal>::iterator
     : public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
 public:
   using parent = iterator_internal<R, iterator>;
   using value_type = typename parent::value_type;
   using pointer = typename parent::pointer;
   using reference = typename parent::reference;
   using difference_type = typename parent::difference_type;
   using iterator_category = typename parent::iterator_category;
 
 public:
   iterator() : parent(){};
   iterator(const iterator & it) = default;
   iterator(iterator && it) = default;
 
   template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
   iterator(P * data) : parent(data) {}
 
   template <typename UP_P,
             typename =
                 std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
   iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
 
   iterator & operator=(const iterator & it) = default;
 
   operator const_iterator<R>() {
     return ConstConverterIteratorHelper<T, R>::convert(*this);
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <class T, bool is_scal>
 template <typename R>
 inline typename Array<T, is_scal>::template iterator<R>
 Array<T, is_scal>::erase(const iterator<R> & it) {
   T * curr = it.data();
   UInt pos = (curr - values) / nb_component;
   erase(pos);
   iterator<R> rit = it;
   return --rit;
 }
 #endif
 
 } // namespace akantu
 
 #endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in
index 8f7e5057c..ecbec6274 100644
--- a/src/common/aka_config.hh.in
+++ b/src/common/aka_config.hh.in
@@ -1,103 +1,103 @@
 /**
  * @file   aka_config.hh.in
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Sun Sep 26 2010
  * @date last modification: Thu Jan 21 2016
  *
  * @brief  Compilation time configuration of Akantu
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_AKA_CONFIG_HH__
 #define __AKANTU_AKA_CONFIG_HH__
 
 #define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@
 #define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@
 #define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@
-#define AKANTU_VERSION (AKANTU_VERSION_MAJOR * 100000 \
-                        + AKANTU_VERSION_MINOR * 1000 \
+#define AKANTU_VERSION (AKANTU_VERSION_MAJOR * 10000 \
+                        + AKANTU_VERSION_MINOR * 100 \
                         + AKANTU_VERSION_PATCH)
 
 @AKANTU_TYPES_EXTRA_INCLUDES@
 namespace akantu {
 using Real = @AKANTU_FLOAT_TYPE@;
 using Int = @AKANTU_SIGNED_INTEGER_TYPE@;
 using UInt = @AKANTU_UNSIGNED_INTEGER_TYPE@;
 } // akantu
 
 #define AKANTU_INTEGER_SIZE @AKANTU_INTEGER_SIZE@
 #define AKANTU_FLOAT_SIZE @AKANTU_FLOAT_SIZE@
 
 #cmakedefine AKANTU_HAS_STRDUP
 
 #cmakedefine AKANTU_USE_BLAS
 #cmakedefine AKANTU_USE_LAPACK
 
 #cmakedefine AKANTU_PARALLEL
 #cmakedefine AKANTU_USE_MPI
 
 #cmakedefine AKANTU_USE_SCOTCH
 #cmakedefine AKANTU_USE_PTSCOTCH
 #cmakedefine AKANTU_SCOTCH_NO_EXTERN
 
 #cmakedefine AKANTU_IMPLICIT
 #cmakedefine AKANTU_USE_MUMPS
 #cmakedefine AKANTU_USE_PETSC
 
 #cmakedefine AKANTU_USE_IOHELPER
 #cmakedefine AKANTU_USE_QVIEW
 #cmakedefine AKANTU_USE_BLACKDYNAMITE
 
 #cmakedefine AKANTU_USE_PYBIND11
 
 #cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY
 
 #cmakedefine AKANTU_EXTRA_MATERIALS
 #cmakedefine AKANTU_STUDENTS_EXTRA_PACKAGE
 #cmakedefine AKANTU_DAMAGE_NON_LOCAL
 
 #cmakedefine AKANTU_SOLID_MECHANICS
 #cmakedefine AKANTU_STRUCTURAL_MECHANICS
 #cmakedefine AKANTU_HEAT_TRANSFER
 #cmakedefine AKANTU_PYTHON_INTERFACE
 
 #cmakedefine AKANTU_COHESIVE_ELEMENT
 #cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT
 
 #cmakedefine AKANTU_IGFEM
 
 #cmakedefine AKANTU_USE_CGAL
 #cmakedefine AKANTU_EMBEDDED
 
 // Debug tools
 //#cmakedefine AKANTU_NDEBUG
 #cmakedefine AKANTU_DEBUG_TOOLS
 #cmakedefine READLINK_COMMAND @READLINK_COMMAND@
 #cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@
 
 #define __aka_inline__ inline
 
 #endif /* __AKANTU_AKA_CONFIG_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index fcc686edf..3693383f0 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1194 +1,1195 @@
 /**
  * @file   solid_mechanics_model.cc
  *
  * @author Ramin Aghababaei <ramin.aghababaei@epfl.ch>
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author David Simon Kammer <david.kammer@epfl.ch>
  * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Clement Roux <clement.roux@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Tue Jul 27 2010
  * @date last modification: Tue Jan 19 2016
  *
  * @brief  Implementation of the SolidMechanicsModel class
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 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 "solid_mechanics_model.hh"
 #include "integrator_gauss.hh"
 #include "shape_lagrange.hh"
 #include "solid_mechanics_model_tmpl.hh"
 
 #include "communicator.hh"
 #include "element_synchronizer.hh"
 #include "sparse_matrix.hh"
 #include "synchronizer_registry.hh"
 
 #include "dumpable_inline_impl.hh"
 #ifdef AKANTU_USE_IOHELPER
 #include "dumper_iohelper_paraview.hh"
 #endif
 
 #include "material_non_local.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /**
  * A solid mechanics model need a mesh  and a dimension to be created. the model
  * by it  self can not  do a lot,  the good init  functions should be  called in
  * order to configure the model depending on what we want to do.
  *
  * @param  mesh mesh  representing  the model  we  want to  simulate
  * @param dim spatial  dimension of the problem, if dim =  0 (default value) the
  * dimension of the problem is assumed to be the on of the mesh
  * @param id an id to identify the model
  */
 SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id,
                                          const MemoryID & memory_id, const ModelType model_type)
     : Model(mesh, model_type, dim, id, memory_id),
       BoundaryCondition<SolidMechanicsModel>(), f_m2a(1.0),
       displacement(nullptr), previous_displacement(nullptr),
       displacement_increment(nullptr), mass(nullptr), velocity(nullptr),
       acceleration(nullptr), external_force(nullptr), internal_force(nullptr),
       blocked_dofs(nullptr), current_position(nullptr),
       material_index("material index", id, memory_id),
       material_local_numbering("material local numbering", id, memory_id),
       increment_flag(false),
       are_materials_instantiated(false) { //, pbc_synch(nullptr) {
   AKANTU_DEBUG_IN();
 
   this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh,
                                                Model::spatial_dimension);
 
 #if defined(AKANTU_USE_IOHELPER)
   this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
   this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
                          _ek_regular);
 #endif
 
   material_selector = std::make_shared<DefaultMaterialSelector>(material_index),
 
   this->initDOFManager();
 
   this->registerDataAccessor(*this);
 
   if (this->mesh.isDistributed()) {
     auto & synchronizer = this->mesh.getElementSynchronizer();
     this->registerSynchronizer(synchronizer, _gst_material_id);
     this->registerSynchronizer(synchronizer, _gst_smm_mass);
     this->registerSynchronizer(synchronizer, _gst_smm_stress);
     this->registerSynchronizer(synchronizer, _gst_smm_boundary);
     this->registerSynchronizer(synchronizer, _gst_for_dump);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 SolidMechanicsModel::~SolidMechanicsModel() {
   AKANTU_DEBUG_IN();
 
   for (auto & internal : this->registered_internals) {
     delete internal.second;
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) {
   Model::setTimeStep(time_step, solver_id);
 
 #if defined(AKANTU_USE_IOHELPER)
   this->mesh.getDumper().setTimeStep(time_step);
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 /* Initialization                                                             */
 /* -------------------------------------------------------------------------- */
 /**
  * This function groups  many of the initialization in on  function. For most of
  * basics  case the  function should  be  enough. The  functions initialize  the
  * model, the internal  vectors, set them to 0, and  depending on the parameters
  * it also initialize the explicit or implicit solver.
  *
  * @param material_file the  file containing the materials to  use
  * @param method the analysis method wanted.  See the akantu::AnalysisMethod for
  * the different possibilities
  */
 void SolidMechanicsModel::initFullImpl(const ModelOptions & options) {
   material_index.initialize(mesh, _element_kind = _ek_not_defined,
                             _default_value = UInt(-1), _with_nb_element = true);
   material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
                                       _with_nb_element = true);
 
   Model::initFullImpl(options);
 
   // initialize pbc
   if (this->pbc_pair.size() != 0)
     this->initPBC();
 
   // initialize the materials
   if (this->parser.getLastParsedFile() != "") {
     this->instantiateMaterials();
   }
 
   this->initMaterials();
 
   this->initBC(*this, *displacement, *displacement_increment, *external_force);
 }
 
 /* -------------------------------------------------------------------------- */
 TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const {
   return _tsst_dynamic_lumped;
 }
 
 /* -------------------------------------------------------------------------- */
 ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions(
     const TimeStepSolverType & type) const {
   ModelSolverOptions options;
 
   switch (type) {
   case _tsst_dynamic_lumped: {
     options.non_linear_solver_type = _nls_lumped;
     options.integration_scheme_type["displacement"] = _ist_central_difference;
     options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     break;
   }
   case _tsst_static: {
     options.non_linear_solver_type = _nls_newton_raphson;
     options.integration_scheme_type["displacement"] = _ist_pseudo_time;
     options.solution_type["displacement"] = IntegrationScheme::_not_defined;
     break;
   }
   case _tsst_dynamic: {
     if (this->method == _explicit_consistent_mass) {
       options.non_linear_solver_type = _nls_newton_raphson;
       options.integration_scheme_type["displacement"] = _ist_central_difference;
       options.solution_type["displacement"] = IntegrationScheme::_acceleration;
     } else {
       options.non_linear_solver_type = _nls_newton_raphson;
       options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2;
       options.solution_type["displacement"] = IntegrationScheme::_displacement;
     }
     break;
   }
   default:
     AKANTU_EXCEPTION(type << " is not a valid time step solver type");
   }
 
   return options;
 }
 
 /* -------------------------------------------------------------------------- */
 std::tuple<ID, TimeStepSolverType>
 SolidMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) {
   switch (method) {
   case _explicit_lumped_mass: {
     return std::make_tuple("explicit_lumped", _tsst_dynamic_lumped);
   }
   case _explicit_consistent_mass: {
     return std::make_tuple("explicit", _tsst_dynamic);
   }
   case _static: {
     return std::make_tuple("static", _tsst_static);
   }
   case _implicit_dynamic: {
     return std::make_tuple("implicit", _tsst_dynamic);
   }
   default:
     return std::make_tuple("unknown", _tsst_not_defined);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type,
                                      NonLinearSolverType) {
   auto & dof_manager = this->getDOFManager();
 
   /* ------------------------------------------------------------------------ */
   // for alloc type of solvers
   this->allocNodalField(this->displacement, spatial_dimension, "displacement");
   this->allocNodalField(this->previous_displacement, spatial_dimension,
                         "previous_displacement");
   this->allocNodalField(this->displacement_increment, spatial_dimension,
                         "displacement_increment");
   this->allocNodalField(this->internal_force, spatial_dimension,
                         "internal_force");
   this->allocNodalField(this->external_force, spatial_dimension,
                         "external_force");
   this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs");
   this->allocNodalField(this->current_position, spatial_dimension,
                         "current_position");
 
   // initialize the current positions
   this->current_position->copy(this->mesh.getNodes());
 
   /* ------------------------------------------------------------------------ */
   if (!dof_manager.hasDOFs("displacement")) {
     dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal);
     dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs);
     dof_manager.registerDOFsIncrement("displacement",
                                       *this->displacement_increment);
     dof_manager.registerDOFsPrevious("displacement",
                                      *this->previous_displacement);
   }
 
   /* ------------------------------------------------------------------------ */
   // for dynamic
   if (time_step_solver_type == _tsst_dynamic ||
       time_step_solver_type == _tsst_dynamic_lumped) {
     this->allocNodalField(this->velocity, spatial_dimension, "velocity");
     this->allocNodalField(this->acceleration, spatial_dimension,
                           "acceleration");
 
     if (!dof_manager.hasDOFsDerivatives("displacement", 1)) {
       dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity);
       dof_manager.registerDOFsDerivative("displacement", 2,
                                          *this->acceleration);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /**
  * Initialize the model,basically it  pre-compute the shapes, shapes derivatives
  * and jacobian
  */
 void SolidMechanicsModel::initModel() {
   /// \todo add  the current position  as a parameter to  initShapeFunctions for
   /// large deformation
   getFEEngine().initShapeFunctions(_not_ghost);
   getFEEngine().initShapeFunctions(_ghost);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleResidual() {
   AKANTU_DEBUG_IN();
 
   /* ------------------------------------------------------------------------ */
   // computes the internal forces
   this->assembleInternalForces();
 
   /* ------------------------------------------------------------------------ */
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->external_force, 1);
   this->getDOFManager().assembleToResidual("displacement",
                                            *this->internal_force, 1);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) {
   // \TODO check the materials to know what is the correct answer
   if (matrix_id == "C")
     return _mt_not_defined;
 
   return _symmetric;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) {
   if (matrix_id == "K") {
     this->assembleStiffnessMatrix();
   } else if (matrix_id == "M") {
     this->assembleMass();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) {
   if (matrix_id == "M") {
     this->assembleMassLumped();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::beforeSolveStep() {
   for (auto & material : materials)
     material->beforeSolveStep();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::afterSolveStep() {
   for (auto & material : materials)
     material->afterSolveStep();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::predictor() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::corrector() { ++displacement_release; }
 
 /* -------------------------------------------------------------------------- */
 /**
  * This function computes the internal forces as F_{int} = \int_{\Omega} N
  * \sigma d\Omega@f$
  */
 void SolidMechanicsModel::assembleInternalForces() {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the internal forces");
 
   this->internal_force->clear();
 
   // compute the stresses of local elements
   AKANTU_DEBUG_INFO("Compute local stresses");
   for (auto & material : materials) {
     material->computeAllStresses(_not_ghost);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Computation of the non local part */
   if (this->non_local_manager)
     this->non_local_manager->computeAllNonLocalStresses();
 
   // communicate the stresses
   AKANTU_DEBUG_INFO("Send data for residual assembly");
   this->asynchronousSynchronize(_gst_smm_stress);
 
   // assemble the forces due to local stresses
   AKANTU_DEBUG_INFO("Assemble residual for local elements");
   for (auto & material : materials) {
     material->assembleInternalForces(_not_ghost);
   }
 
   // finalize communications
   AKANTU_DEBUG_INFO("Wait distant stresses");
   this->waitEndSynchronize(_gst_smm_stress);
 
   // assemble the stresses due to ghost elements
   AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
   for (auto & material : materials) {
     material->assembleInternalForces(_ghost);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::assembleStiffnessMatrix() {
   AKANTU_DEBUG_IN();
 
   AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
 
   // Check if materials need to recompute the matrix
   bool need_to_reassemble = false;
 
   for (auto & material : materials) {
     need_to_reassemble |= material->hasStiffnessMatrixChanged();
   }
 
   if (need_to_reassemble) {
     this->getDOFManager().getMatrix("K").clear();
 
     // call compute stiffness matrix on each local elements
     for (auto & material : materials) {
       material->assembleStiffnessMatrix(_not_ghost);
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateCurrentPosition() {
   if (this->current_position_release == this->displacement_release)
     return;
 
   this->current_position->copy(this->mesh.getNodes());
 
   auto cpos_it = this->current_position->begin(Model::spatial_dimension);
   auto cpos_end = this->current_position->end(Model::spatial_dimension);
   auto disp_it = this->displacement->begin(Model::spatial_dimension);
 
   for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) {
     *cpos_it += *disp_it;
   }
 
   this->current_position_release = this->displacement_release;
 }
 
 /* -------------------------------------------------------------------------- */
 const Array<Real> & SolidMechanicsModel::getCurrentPosition() {
   this->updateCurrentPosition();
   return *this->current_position;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateDataForNonLocalCriterion(
     ElementTypeMapReal & criterion) {
   const ID field_name = criterion.getName();
   for (auto & material : materials) {
     if (!material->isInternal<Real>(field_name, _ek_regular))
       continue;
 
     for (auto ghost_type : ghost_types) {
       material->flattenInternal(field_name, criterion, ghost_type, _ek_regular);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /* Information                                                                */
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep() {
   AKANTU_DEBUG_IN();
 
   Real min_dt = getStableTimeStep(_not_ghost);
 
   /// reduction min over all processors
   mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min);
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
   AKANTU_DEBUG_IN();
 
   Real min_dt = std::numeric_limits<Real>::max();
 
   this->updateCurrentPosition();
 
   Element elem;
   elem.ghost_type = ghost_type;
 
   for (auto type :
        mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) {
     elem.type = type;
     UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);
     UInt nb_element = mesh.getNbElement(type);
 
     auto mat_indexes = material_index(type, ghost_type).begin();
     auto mat_loc_num = material_local_numbering(type, ghost_type).begin();
 
     Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension);
     FEEngine::extractNodalToElementField(mesh, *current_position, X, type,
                                          _not_ghost);
 
     auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element);
 
     for (UInt el = 0; el < nb_element;
          ++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
       elem.element = *mat_loc_num;
       Real el_h = getFEEngine().getElementInradius(*X_el, type);
       Real el_c = this->materials[*mat_indexes]->getCelerity(elem);
       Real el_dt = el_h / el_c;
 
       min_dt = std::min(min_dt, el_dt);
     }
   }
 
   AKANTU_DEBUG_OUT();
   return min_dt;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy() {
   AKANTU_DEBUG_IN();
 
   Real ekin = 0.;
   UInt nb_nodes = mesh.getNbNodes();
 
   if (this->getDOFManager().hasLumpedMatrix("M")) {
     auto m_it = this->mass->begin(Model::spatial_dimension);
     auto m_end = this->mass->end(Model::spatial_dimension);
     auto v_it = this->velocity->begin(Model::spatial_dimension);
 
     for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) {
       const Vector<Real> & v = *v_it;
       const Vector<Real> & m = *m_it;
 
       Real mv2 = 0;
       bool is_local_node = mesh.isLocalOrMasterNode(n);
       // bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
       bool count_node = is_local_node; // && is_not_pbc_slave_node;
       if (count_node) {
         for (UInt i = 0; i < Model::spatial_dimension; ++i) {
           if (m(i) > std::numeric_limits<Real>::epsilon())
             mv2 += v(i) * v(i) * m(i);
         }
       }
 
       ekin += mv2;
     }
   } else if (this->getDOFManager().hasMatrix("M")) {
     Array<Real> Mv(nb_nodes, Model::spatial_dimension);
     this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv);
 
     auto mv_it = Mv.begin(Model::spatial_dimension);
     auto mv_end = Mv.end(Model::spatial_dimension);
     auto v_it = this->velocity->begin(Model::spatial_dimension);
 
     for (; mv_it != mv_end; ++mv_it, ++v_it) {
       ekin += v_it->dot(*mv_it);
     }
   } else {
     AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
   }
 
   mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return ekin * .5;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getKineticEnergy(const ElementType & type,
                                            UInt index) {
   AKANTU_DEBUG_IN();
 
   UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type);
 
   Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension);
   Array<UInt> filter_element(1, 1, index);
 
   getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad,
                                                Model::spatial_dimension, type,
                                                _not_ghost, filter_element);
 
   Array<Real>::vector_iterator vit =
       vel_on_quad.begin(Model::spatial_dimension);
   Array<Real>::vector_iterator vend = vel_on_quad.end(Model::spatial_dimension);
 
   Vector<Real> rho_v2(nb_quadrature_points);
 
   Real rho = materials[material_index(type)(index)]->getRho();
 
   for (UInt q = 0; vit != vend; ++vit, ++q) {
     rho_v2(q) = rho * vit->dot(*vit);
   }
 
   AKANTU_DEBUG_OUT();
 
   return .5 * getFEEngine().integrate(rho_v2, type, index);
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getExternalWork() {
   AKANTU_DEBUG_IN();
 
   auto ext_force_it = external_force->begin(Model::spatial_dimension);
   auto int_force_it = internal_force->begin(Model::spatial_dimension);
   auto boun_it = blocked_dofs->begin(Model::spatial_dimension);
 
   decltype(ext_force_it) incr_or_velo_it;
   if (this->method == _static) {
     incr_or_velo_it =
         this->displacement_increment->begin(Model::spatial_dimension);
   } else {
     incr_or_velo_it = this->velocity->begin(Model::spatial_dimension);
   }
 
   Real work = 0.;
 
   UInt nb_nodes = this->mesh.getNbNodes();
 
   for (UInt n = 0; n < nb_nodes;
        ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) {
     const auto & int_force = *int_force_it;
     const auto & ext_force = *ext_force_it;
     const auto & boun = *boun_it;
     const auto & incr_or_velo = *incr_or_velo_it;
 
     bool is_local_node = this->mesh.isLocalOrMasterNode(n);
     // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n);
     bool count_node = is_local_node; // && is_not_pbc_slave_node;
 
     if (count_node) {
       for (UInt i = 0; i < Model::spatial_dimension; ++i) {
         if (boun(i))
           work -= int_force(i) * incr_or_velo(i);
         else
           work += ext_force(i) * incr_or_velo(i);
       }
     }
   }
 
   mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum);
 
   if (this->method != _static)
     work *= this->getTimeStep();
 
   AKANTU_DEBUG_OUT();
   return work;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy();
   } else if (energy_id == "external work") {
     return getExternalWork();
   }
 
   Real energy = 0.;
   for (auto & material : materials)
     energy += material->getEnergy(energy_id);
 
   /// reduction sum over all processors
   mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 /* -------------------------------------------------------------------------- */
 Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
                                     const ElementType & type, UInt index) {
   AKANTU_DEBUG_IN();
 
   if (energy_id == "kinetic") {
     return getKineticEnergy(type, index);
   }
 
   UInt mat_index = this->material_index(type, _not_ghost)(index);
   UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
   Real energy =
       this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
 
   AKANTU_DEBUG_OUT();
   return energy;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
                                           const NewElementsEvent & event) {
   AKANTU_DEBUG_IN();
 
   this->material_index.initialize(mesh, _element_kind = _ek_not_defined,
                                   _with_nb_element = true,
                                   _default_value = UInt(-1));
   this->material_local_numbering.initialize(
       mesh, _element_kind = _ek_not_defined, _with_nb_element = true,
       _default_value = UInt(-1));
 
   ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(),
                                    this->getMemoryID());
 
   for (auto & elem : element_list) {
     if (!filter.exists(elem.type, elem.ghost_type))
       filter.alloc(0, 1, elem.type, elem.ghost_type);
     filter(elem.type, elem.ghost_type).push_back(elem.element);
   }
 
   this->assignMaterialToElements(&filter);
 
   for (auto & material : materials)
     material->onElementsAdded(element_list, event);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onElementsRemoved(
     const Array<Element> & element_list,
     const ElementTypeMapArray<UInt> & new_numbering,
     const RemovedElementsEvent & event) {
   for (auto & material : materials) {
     material->onElementsRemoved(element_list, new_numbering, event);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
                                        const NewNodesEvent & event) {
   AKANTU_DEBUG_IN();
   UInt nb_nodes = mesh.getNbNodes();
 
   if (displacement) {
     displacement->resize(nb_nodes, 0.);
     ++displacement_release;
   }
   if (mass)
     mass->resize(nb_nodes, 0.);
   if (velocity)
     velocity->resize(nb_nodes, 0.);
   if (acceleration)
     acceleration->resize(nb_nodes, 0.);
   if (external_force)
     external_force->resize(nb_nodes, 0.);
   if (internal_force)
     internal_force->resize(nb_nodes, 0.);
   if (blocked_dofs)
     blocked_dofs->resize(nb_nodes, 0.);
   if (current_position)
     current_position->resize(nb_nodes, 0.);
 
   if (previous_displacement)
     previous_displacement->resize(nb_nodes, 0.);
   if (displacement_increment)
     displacement_increment->resize(nb_nodes, 0.);
 
   for (auto & material : materials) {
     material->onNodesAdded(nodes_list, event);
   }
 
   need_to_reassemble_lumped_mass = true;
   need_to_reassemble_mass = true;
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::onNodesRemoved(const Array<UInt> & /*element_list*/,
                                          const Array<UInt> & new_numbering,
                                          const RemovedNodesEvent & /*event*/) {
   if (displacement) {
     mesh.removeNodesFromArray(*displacement, new_numbering);
     ++displacement_release;
   }
   if (mass)
     mesh.removeNodesFromArray(*mass, new_numbering);
   if (velocity)
     mesh.removeNodesFromArray(*velocity, new_numbering);
   if (acceleration)
     mesh.removeNodesFromArray(*acceleration, new_numbering);
   if (internal_force)
     mesh.removeNodesFromArray(*internal_force, new_numbering);
   if (external_force)
     mesh.removeNodesFromArray(*external_force, new_numbering);
   if (blocked_dofs)
     mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
 
   // if (increment_acceleration)
   //   mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
   if (displacement_increment)
     mesh.removeNodesFromArray(*displacement_increment, new_numbering);
 
   if (previous_displacement)
     mesh.removeNodesFromArray(*previous_displacement, new_numbering);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
   std::string space;
   for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
     ;
 
   stream << space << "Solid Mechanics Model [" << std::endl;
   stream << space << " + id                : " << id << std::endl;
   stream << space << " + spatial dimension : " << Model::spatial_dimension
          << std::endl;
   stream << space << " + fem [" << std::endl;
   getFEEngine().printself(stream, indent + 2);
   stream << space << AKANTU_INDENT << "]" << std::endl;
   stream << space << " + nodals information [" << std::endl;
   displacement->printself(stream, indent + 2);
   if (velocity)
     velocity->printself(stream, indent + 2);
   if (acceleration)
     acceleration->printself(stream, indent + 2);
   if (mass)
     mass->printself(stream, indent + 2);
   external_force->printself(stream, indent + 2);
   internal_force->printself(stream, indent + 2);
   blocked_dofs->printself(stream, indent + 2);
   stream << space << AKANTU_INDENT << "]" << std::endl;
 
   stream << space << " + material information [" << std::endl;
   material_index.printself(stream, indent + 2);
   stream << space << AKANTU_INDENT << "]" << std::endl;
 
   stream << space << " + materials [" << std::endl;
   for (auto & material : materials) {
     material->printself(stream, indent + 1);
   }
   stream << space << AKANTU_INDENT << "]" << std::endl;
 
   stream << space << "]" << std::endl;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::initializeNonLocal() {
   this->non_local_manager->synchronize(*this, _gst_material_id);
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods(
     const GhostType & ghost_type) {
   for (auto & mat : materials) {
     MaterialNonLocalInterface * mat_non_local;
     if ((mat_non_local =
              dynamic_cast<MaterialNonLocalInterface *>(mat.get())) == nullptr)
       continue;
 
     ElementTypeMapArray<Real> quadrature_points_coordinates(
         "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id);
     quadrature_points_coordinates.initialize(this->getFEEngine(),
                                              _nb_component = spatial_dimension,
                                              _ghost_type = ghost_type);
 
     for (auto & type : quadrature_points_coordinates.elementTypes(
              Model::spatial_dimension, ghost_type)) {
       this->getFEEngine().computeIntegrationPointsCoordinates(
           quadrature_points_coordinates(type, ghost_type), type, ghost_type);
     }
 
     mat_non_local->initMaterialNonLocal();
 
     mat_non_local->insertIntegrationPointsInNeighborhoods(
         ghost_type, quadrature_points_coordinates);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::computeNonLocalStresses(
     const GhostType & ghost_type) {
   for (auto & mat : materials) {
     if (dynamic_cast<MaterialNonLocalInterface *>(mat.get()) == nullptr)
       continue;
 
     auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
     mat_non_local.computeNonLocalStresses(ghost_type);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateLocalInternal(
     ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
     const ElementKind & kind) {
   const ID field_name = internal_flat.getName();
   for (auto & material : materials) {
     if (material->isInternal<Real>(field_name, kind))
       material->flattenInternal(field_name, internal_flat, ghost_type, kind);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::updateNonLocalInternal(
     ElementTypeMapReal & internal_flat, const GhostType & ghost_type,
     const ElementKind & kind) {
 
   const ID field_name = internal_flat.getName();
 
   for (auto & mat : materials) {
     if (dynamic_cast<MaterialNonLocalInterface *>(mat.get()) == nullptr)
       continue;
 
     auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat);
     mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type,
                                           kind);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) {
   return dynamic_cast<FEEngine &>(
       getFEEngineClassBoundary<MyFEEngineType>(name));
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::splitElementByMaterial(
     const Array<Element> & elements,
     std::vector<Array<Element>> & elements_per_mat) const {
   ElementType current_element_type = _not_defined;
   GhostType current_ghost_type = _casper;
   const Array<UInt> * mat_indexes = nullptr;
   const Array<UInt> * mat_loc_num = nullptr;
 
-  for (auto && el : elements) {
+  for (const auto & el : elements) {
     if (el.type != current_element_type ||
         el.ghost_type != current_ghost_type) {
       current_element_type = el.type;
       current_ghost_type = el.ghost_type;
       mat_indexes = &(this->material_index(el.type, el.ghost_type));
       mat_loc_num = &(this->material_local_numbering(el.type, el.ghost_type));
     }
 
-    UInt old_id = el.element;
-    el.element = (*mat_loc_num)(old_id);
-    elements_per_mat[(*mat_indexes)(old_id)].push_back(el);
+    Element mat_el = el;
+
+    mat_el.element = (*mat_loc_num)(el.element);
+    elements_per_mat[(*mat_indexes)(el.element)].push_back(mat_el);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<Element> & elements,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   UInt nb_nodes_per_element = 0;
 
   for (const Element & el : elements) {
     nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
   }
 
   switch (tag) {
   case _gst_material_id: {
     size += elements.size() * sizeof(UInt);
     break;
   }
   case _gst_smm_mass: {
     size += nb_nodes_per_element * sizeof(Real) *
             Model::spatial_dimension; // mass vector
     break;
   }
   case _gst_smm_for_gradu: {
     size += nb_nodes_per_element * Model::spatial_dimension *
             sizeof(Real); // displacement
     break;
   }
   case _gst_smm_boundary: {
     // force, displacement, boundary
     size += nb_nodes_per_element * Model::spatial_dimension *
             (2 * sizeof(Real) + sizeof(bool));
     break;
   }
   case _gst_for_dump: {
     // displacement, velocity, acceleration, residual, force
     size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5;
     break;
   }
   default: {}
   }
 
   if (tag != _gst_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       size += mat.getNbData(elements, tag);
     });
   }
 
   AKANTU_DEBUG_OUT();
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<Element> & elements,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case _gst_material_id: {
     this->packElementalDataHelper(material_index, buffer, elements, false,
                                   getFEEngine());
     break;
   }
   case _gst_smm_mass: {
     packNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case _gst_smm_for_gradu: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case _gst_for_dump: {
     packNodalDataHelper(*displacement, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*acceleration, buffer, elements, mesh);
     packNodalDataHelper(*internal_force, buffer, elements, mesh);
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case _gst_smm_boundary: {
     packNodalDataHelper(*external_force, buffer, elements, mesh);
     packNodalDataHelper(*velocity, buffer, elements, mesh);
     packNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {}
   }
 
   if (tag != _gst_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       mat.packData(buffer, elements, tag);
     });
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<Element> & elements,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case _gst_material_id: {
     for (auto && element : elements) {
       UInt recv_mat_index;
       buffer >> recv_mat_index;
       UInt & mat_index = material_index(element);
       if (mat_index != UInt(-1))
         continue;
 
       // add ghosts element to the correct material
       mat_index = recv_mat_index;
       UInt index = materials[mat_index]->addElement(element);
       material_local_numbering(element) = index;
     }
     break;
   }
   case _gst_smm_mass: {
     unpackNodalDataHelper(*mass, buffer, elements, mesh);
     break;
   }
   case _gst_smm_for_gradu: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     break;
   }
   case _gst_for_dump: {
     unpackNodalDataHelper(*displacement, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*acceleration, buffer, elements, mesh);
     unpackNodalDataHelper(*internal_force, buffer, elements, mesh);
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     break;
   }
   case _gst_smm_boundary: {
     unpackNodalDataHelper(*external_force, buffer, elements, mesh);
     unpackNodalDataHelper(*velocity, buffer, elements, mesh);
     unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh);
     break;
   }
   default: {}
   }
 
   if (tag != _gst_material_id) {
     splitByMaterial(elements, [&](auto && mat, auto && elements) {
       mat.unpackData(buffer, elements, tag);
     });
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs,
                                     const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   UInt size = 0;
   //  UInt nb_nodes = mesh.getNbNodes();
 
   switch (tag) {
   case _gst_smm_uv: {
     size += sizeof(Real) * Model::spatial_dimension * 2;
     break;
   }
   case _gst_smm_res: {
     size += sizeof(Real) * Model::spatial_dimension;
     break;
   }
   case _gst_smm_mass: {
     size += sizeof(Real) * Model::spatial_dimension;
     break;
   }
   case _gst_for_dump: {
     size += sizeof(Real) * Model::spatial_dimension * 5;
     break;
   }
   default: {
     AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
   return size * dofs.size();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
                                    const Array<UInt> & dofs,
                                    const SynchronizationTag & tag) const {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case _gst_smm_uv: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case _gst_smm_res: {
     packDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case _gst_smm_mass: {
     packDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case _gst_for_dump: {
     packDOFDataHelper(*displacement, buffer, dofs);
     packDOFDataHelper(*velocity, buffer, dofs);
     packDOFDataHelper(*acceleration, buffer, dofs);
     packDOFDataHelper(*internal_force, buffer, dofs);
     packDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
                                      const Array<UInt> & dofs,
                                      const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   switch (tag) {
   case _gst_smm_uv: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     break;
   }
   case _gst_smm_res: {
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     break;
   }
   case _gst_smm_mass: {
     unpackDOFDataHelper(*mass, buffer, dofs);
     break;
   }
   case _gst_for_dump: {
     unpackDOFDataHelper(*displacement, buffer, dofs);
     unpackDOFDataHelper(*velocity, buffer, dofs);
     unpackDOFDataHelper(*acceleration, buffer, dofs);
     unpackDOFDataHelper(*internal_force, buffer, dofs);
     unpackDOFDataHelper(*external_force, buffer, dofs);
     break;
   }
   default: {
     AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
   }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 } // namespace akantu
diff --git a/src/synchronizer/communication_descriptor.hh b/src/synchronizer/communication_descriptor.hh
index e3f24ad99..609218b86 100644
--- a/src/synchronizer/communication_descriptor.hh
+++ b/src/synchronizer/communication_descriptor.hh
@@ -1,153 +1,153 @@
 /**
  * @file   synchronizer_tmpl.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date   Wed Aug  3 13:49:36 2016
  *
  * @brief  Implementation of the helper classes for the synchronizer
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2011 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_array.hh"
 #include "communication_tag.hh"
 #include "data_accessor.hh"
 #include "communication_request.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_COMMUNICATION_DESCRIPTOR_HH__
 #define __AKANTU_COMMUNICATION_DESCRIPTOR_HH__
 
 namespace akantu {
 
 /* ------------------------------------------------------------------------ */
 enum CommunicationSendRecv  { _send, _recv, _csr_not_defined };
 
 /* -------------------------------------------------------------------------- */
 struct CommunicationSRType {
   using type = CommunicationSendRecv;
   static const type _begin_ = _send;
   static const type _end_ = _csr_not_defined;
 };
 
 using send_recv_t = safe_enum<CommunicationSRType>;
 
 namespace {
   send_recv_t iterate_send_recv{};
 }
 
 /* ------------------------------------------------------------------------ */
 class Communication {
 public:
   explicit Communication(const CommunicationSendRecv & type = _csr_not_defined)
       : _type(type) {}
 
   Communication(const Communication &) = delete;
   Communication & operator=(const Communication &) = delete;
 
   void resize(UInt size) {
     this->_size = size;
     this->_buffer.resize(size);
   }
 
   inline const CommunicationSendRecv & type() const { return this->_type; }
   inline const UInt & size() const { return this->_size; }
 
   inline const CommunicationRequest & request() const { return this->_request; }
   inline CommunicationRequest & request() { return this->_request; }
 
   inline const CommunicationBuffer & buffer() const { return this->_buffer; }
   inline CommunicationBuffer & buffer() { return this->_buffer; }
 
 private:
   UInt _size{0};
   CommunicationBuffer _buffer;
   CommunicationRequest _request;
   CommunicationSendRecv _type;
 };
 
 template <class Entity> class Communications;
 
 /* ------------------------------------------------------------------------ */
 template <class Entity> class CommunicationDescriptor {
 public:
   CommunicationDescriptor(Communication & communication, Array<Entity> & scheme,
                           Communications<Entity> & communications,
                           const SynchronizationTag & tag, UInt proc);
 
   CommunicationDescriptor(const CommunicationDescriptor &) = default;
 
   CommunicationDescriptor &
   operator=(const CommunicationDescriptor &) = default;
 
   /// get the quantity of data in the buffer
   UInt getNbData() { return communication.size(); }
   /// set the quantity of data in the buffer
   void setNbData(UInt size) { communication.resize(size); }
 
   /// get the corresponding tag
   const SynchronizationTag & getTag() const { return tag; }
   /// get the data buffer
   CommunicationBuffer & getBuffer();
 
   /// get the corresponding request
   CommunicationRequest & getRequest();
 
   /// get the communication scheme
-  Array<Entity> & getScheme();
+  const Array<Entity> & getScheme();
 
   /// reset the buffer before pack or after unpack
   void resetBuffer();
 
   /// pack data for entities in the buffer
   void packData(const DataAccessor<Entity> & accessor);
 
   /// unpack data for entities from the buffer
   void unpackData(DataAccessor<Entity> & accessor);
 
   /// posts asynchronous send requests
   void postSend(int hash_id);
 
   /// posts asynchronous receive requests
   void postRecv(int hash_id);
 
   /// free the request
   void freeRequest();
 
   UInt getProc() { return proc; }
 
 protected:
   Communication & communication;
-  Array<Entity> & scheme;
+  const Array<Entity> & scheme;
   Communications<Entity> & communications;
   const SynchronizationTag & tag;
   UInt proc;
   UInt rank;
   UInt counter;
 };
 
 /* -------------------------------------------------------------------------- */
 } // namespace akantu
 
 #include "communication_descriptor_tmpl.hh"
 
 #endif /* __AKANTU_COMMUNICATION_DESCRIPTOR_HH__ */
diff --git a/src/synchronizer/communication_descriptor_tmpl.hh b/src/synchronizer/communication_descriptor_tmpl.hh
index 45c57628f..8e6982fec 100644
--- a/src/synchronizer/communication_descriptor_tmpl.hh
+++ b/src/synchronizer/communication_descriptor_tmpl.hh
@@ -1,152 +1,152 @@
 /**
  * @file   communication_descriptor_tmpl.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date   Tue Sep  6 14:03:16 2016
  *
  * @brief  implementation of CommunicationDescriptor
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2011 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 "communication_descriptor.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__
 #define __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* Implementations                                                            */
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 CommunicationDescriptor<Entity>::CommunicationDescriptor(
     Communication & communication, Array<Entity> & scheme,
     Communications<Entity> & communications, const SynchronizationTag & tag,
     UInt proc)
     : communication(communication), scheme(scheme),
       communications(communications), tag(tag), proc(proc),
       rank(communications.getCommunicator().whoAmI()) {
   counter = communications.getCounter(tag);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 CommunicationBuffer & CommunicationDescriptor<Entity>::getBuffer() {
   return communication.buffer();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 CommunicationRequest & CommunicationDescriptor<Entity>::getRequest() {
   return communication.request();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity> void CommunicationDescriptor<Entity>::freeRequest() {
   const auto & comm = communications.getCommunicator();
 
   // comm.test(communication.request());
   comm.freeCommunicationRequest(communication.request());
 
   communications.decrementPending(tag, communication.type());
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
-Array<Entity> & CommunicationDescriptor<Entity>::getScheme() {
+const Array<Entity> & CommunicationDescriptor<Entity>::getScheme() {
   return scheme;
 }
 
 template <class Entity>
 void CommunicationDescriptor<Entity>::resetBuffer() {
   this->communication.buffer().reset();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void CommunicationDescriptor<Entity>::packData(
     const DataAccessor<Entity> & accessor) {
   AKANTU_DEBUG_ASSERT(
       communication.type() == _send,
       "Cannot pack data on communication that is not of type _send");
   accessor.packData(communication.buffer(), scheme, tag);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void CommunicationDescriptor<Entity>::unpackData(
     DataAccessor<Entity> & accessor) {
   AKANTU_DEBUG_ASSERT(
       communication.type() == _recv,
       "Cannot unpack data from communication that is not of type _recv");
   accessor.unpackData(communication.buffer(), scheme, tag);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void CommunicationDescriptor<Entity>::postSend(int hash_id) {
   AKANTU_DEBUG_ASSERT(communication.type() == _send,
                       "Cannot send a communication that is not of type _send");
   Tag comm_tag = Tag::genTag(rank, counter, tag, hash_id);
   AKANTU_DEBUG_ASSERT(communication.buffer().getPackedSize() ==
                           communication.size(),
                       "a problem have been introduced with "
                           << "false sent sizes declaration "
                           << communication.buffer().getPackedSize()
                           << " != " << communication.size());
 
   AKANTU_DEBUG_INFO("Posting send to proc " << proc << " (tag: " << tag << " - "
                                             << communication.size()
                                             << " data to send) "
                                             << " [ " << comm_tag << " ]");
 
   CommunicationRequest & request = communication.request();
   request = communications.getCommunicator().asyncSend(communication.buffer(),
                                                        proc, comm_tag);
   communications.incrementPending(tag, communication.type());
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void CommunicationDescriptor<Entity>::postRecv(int hash_id) {
   AKANTU_DEBUG_ASSERT(communication.type() == _recv,
                       "Cannot receive data for communication ("
                           << communication.type()
                           << ")that is not of type _recv");
 
   Tag comm_tag = Tag::genTag(proc, counter, tag, hash_id);
   AKANTU_DEBUG_INFO("Posting receive from proc "
                     << proc << " (tag: " << tag << " - " << communication.size()
                     << " data to receive) "
                     << " [ " << comm_tag << " ]");
 
   CommunicationRequest & request = communication.request();
   request = communications.getCommunicator().asyncReceive(
       communication.buffer(), proc, comm_tag);
   communications.incrementPending(tag, communication.type());
 }
 
 } // akantu
 
 #endif /* __AKANTU_COMMUNICATION_DESCRIPTOR_TMPL_HH__ */
diff --git a/src/synchronizer/element_synchronizer.cc b/src/synchronizer/element_synchronizer.cc
index 76a26d983..25a7f77d1 100644
--- a/src/synchronizer/element_synchronizer.cc
+++ b/src/synchronizer/element_synchronizer.cc
@@ -1,320 +1,349 @@
 /**
  * @file   element_synchronizer.cc
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Dana Christen <dana.christen@epfl.ch>
  * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  * @author Marco Vocialta <marco.vocialta@epfl.ch>
  *
  * @date creation: Wed Sep 01 2010
  * @date last modification: Fri Jan 22 2016
  *
  * @brief  implementation of a  communicator using a static_communicator for
  * real
  * send/receive
  *
  * @section LICENSE
  *
  * Copyright (©)  2010-2012, 2014,  2015 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 "element_synchronizer.hh"
 #include "aka_common.hh"
 #include "mesh.hh"
 #include "mesh_utils.hh"
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 #include <iostream>
 #include <map>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 ElementSynchronizer::ElementSynchronizer(Mesh & mesh, const ID & id,
                                          MemoryID memory_id,
                                          bool register_to_event_manager,
                                          EventHandlerPriority event_priority)
     : SynchronizerImpl<Element>(mesh.getCommunicator(), id, memory_id),
       mesh(mesh), element_to_prank("element_to_prank", id, memory_id) {
 
   AKANTU_DEBUG_IN();
 
   if (register_to_event_manager)
     this->mesh.registerEventHandler(*this, event_priority);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 ElementSynchronizer::~ElementSynchronizer() = default;
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::substituteElements(
     const std::map<Element, Element> & old_to_new_elements) {
-  auto found_element_end =
-    old_to_new_elements.end();
+  auto found_element_end = old_to_new_elements.end();
 
   // substitute old elements with new ones
   for (auto && sr : iterate_send_recv) {
     for (auto && scheme_pair : communications.iterateSchemes(sr)) {
       auto & list = scheme_pair.second;
       for (auto & el : list) {
         auto found_element_it = old_to_new_elements.find(el);
         if (found_element_it != found_element_end)
           el = found_element_it->second;
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::onElementsChanged(
     const Array<Element> & old_elements_list,
     const Array<Element> & new_elements_list, const ElementTypeMapArray<UInt> &,
     const ChangedElementsEvent &) {
   // create a map to link old elements to new ones
   std::map<Element, Element> old_to_new_elements;
 
   for (UInt el = 0; el < old_elements_list.size(); ++el) {
     AKANTU_DEBUG_ASSERT(old_to_new_elements.find(old_elements_list(el)) ==
                             old_to_new_elements.end(),
                         "The same element cannot appear twice in the list");
 
     old_to_new_elements[old_elements_list(el)] = new_elements_list(el);
   }
 
   substituteElements(old_to_new_elements);
 
   communications.invalidateSizes();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::onElementsRemoved(
     const Array<Element> & element_to_remove,
     const ElementTypeMapArray<UInt> & new_numbering,
     const RemovedElementsEvent &) {
   AKANTU_DEBUG_IN();
 
   this->removeElements(element_to_remove);
   this->renumberElements(new_numbering);
 
   communications.invalidateSizes();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::buildElementToPrank() {
   AKANTU_DEBUG_IN();
 
   UInt spatial_dimension = mesh.getSpatialDimension();
   element_to_prank.initialize(mesh, _spatial_dimension = spatial_dimension,
                               _element_kind = _ek_not_defined,
                               _with_nb_element = true, _default_value = rank);
 
   /// assign prank to all ghost elements
   for (auto && scheme : communications.iterateSchemes(_recv)) {
     auto & recv = scheme.second;
     auto proc = scheme.first;
 
     for (auto & element : recv) {
       element_to_prank(element) = proc;
     }
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 Int ElementSynchronizer::getRank(const Element & element) const {
-  if(not element_to_prank.exists(element.type, element.ghost_type)) {
+  if (not element_to_prank.exists(element.type, element.ghost_type)) {
     // Nicolas: Ok This is nasty I know....
     const_cast<ElementSynchronizer *>(this)->buildElementToPrank();
   }
 
   return element_to_prank(element);
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::reset() {
   AKANTU_DEBUG_IN();
 
   communications.resetSchemes();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::removeElements(
     const Array<Element> & element_to_remove) {
   AKANTU_DEBUG_IN();
 
   std::vector<CommunicationRequest> send_requests;
   std::map<UInt, Array<UInt>> list_of_elements_per_proc;
 
   auto filter_list = [](const Array<UInt> & keep, Array<Element> & list) {
     Array<Element> new_list;
     for (UInt e = 0; e < keep.size() - 1; ++e) {
       Element & el = list(keep(e));
       new_list.push_back(el);
     }
     list.copy(new_list);
   };
 
   // Handling ghost elements
   for (auto && scheme_pair : communications.iterateSchemes(_recv)) {
     auto proc = scheme_pair.first;
     auto & recv = scheme_pair.second;
     Array<UInt> & keep_list = list_of_elements_per_proc[proc];
 
     auto rem_it = element_to_remove.begin();
     auto rem_end = element_to_remove.end();
 
     for (auto && pair : enumerate(recv)) {
       const auto & el = std::get<1>(pair);
       auto pos = std::find(rem_it, rem_end, el);
 
       if (pos == rem_end) {
         keep_list.push_back(std::get<0>(pair));
       }
     }
 
     keep_list.push_back(UInt(-1)); // To no send empty arrays
 
     auto && tag = Tag::genTag(proc, 0, Tag::_MODIFY_SCHEME, this->hash_id);
     AKANTU_DEBUG_INFO("Sending a message of size "
                       << keep_list.size() << " to proc " << proc << " TAG("
                       << tag << ")");
     send_requests.push_back(this->communicator.asyncSend(keep_list, proc, tag));
 
     auto old_size = recv.size();
     filter_list(keep_list, recv);
 
     AKANTU_DEBUG_INFO("I had " << old_size << " elements to recv from proc "
                                << proc << " and " << keep_list.size()
                                << " elements to keep. I have " << recv.size()
                                << " elements left.");
   }
 
   for (auto && scheme_pair : communications.iterateSchemes(_send)) {
     auto proc = scheme_pair.first;
     auto & send = scheme_pair.second;
 
     CommunicationStatus status;
     auto && tag = Tag::genTag(rank, 0, Tag::_MODIFY_SCHEME, this->hash_id);
     AKANTU_DEBUG_INFO("Getting number of elements of proc "
                       << proc << " to keep - TAG(" << tag << ")");
     this->communicator.probe<UInt>(proc, tag, status);
     Array<UInt> keep_list(status.size());
 
     AKANTU_DEBUG_INFO("Receiving list of elements ("
                       << keep_list.size() << " elements) to keep for proc "
                       << proc << " TAG(" << tag << ")");
     this->communicator.receive(keep_list, proc, tag);
 
     auto old_size = send.size();
     filter_list(keep_list, send);
 
     AKANTU_DEBUG_INFO("I had " << old_size << " elements to send to proc "
                                << proc << " and " << keep_list.size()
                                << " elements to keep. I have " << send.size()
                                << " elements left.");
   }
 
   this->communicator.waitAll(send_requests);
   this->communicator.freeCommunicationRequest(send_requests);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::renumberElements(
     const ElementTypeMapArray<UInt> & new_numbering) {
-  for(auto && sr : iterate_send_recv) {
+  for (auto && sr : iterate_send_recv) {
     for (auto && scheme_pair : communications.iterateSchemes(sr)) {
       auto & list = scheme_pair.second;
       for (auto && el : list) {
-        if(new_numbering.exists(el.type, el.ghost_type))
+        if (new_numbering.exists(el.type, el.ghost_type))
           el.element = new_numbering(el);
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 UInt ElementSynchronizer::sanityCheckDataSize(
     const Array<Element> & elements, const SynchronizationTag &) const {
-  return (elements.size() * mesh.getSpatialDimension() * sizeof(Real) +
-          sizeof(SynchronizationTag));
+  UInt size = 0;
+  size += sizeof(SynchronizationTag); // tag
+  size += sizeof(UInt);               // comm_desc.getNbData();
+  size += sizeof(UInt);               // comm_desc.getProc();
+  size += sizeof(UInt);               // mesh.getCommunicator().whoAmI();
+
+  // barycenters
+  size += (elements.size() * mesh.getSpatialDimension() * sizeof(Real));
+  return size;
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::packSanityCheckData(
     CommunicationDescriptor<Element> & comm_desc) const {
   auto & buffer = comm_desc.getBuffer();
   buffer << comm_desc.getTag();
+  buffer << comm_desc.getNbData();
+  buffer << comm_desc.getProc();
+  buffer << mesh.getCommunicator().whoAmI();
 
   auto & send_element = comm_desc.getScheme();
 
   /// pack barycenters in debug mode
-  for(auto && element : send_element) {
+  for (auto && element : send_element) {
     Vector<Real> barycenter(mesh.getSpatialDimension());
     mesh.getBarycenter(element, barycenter);
     buffer << barycenter;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void ElementSynchronizer::unpackSanityCheckData(
     CommunicationDescriptor<Element> & comm_desc) const {
   auto & buffer = comm_desc.getBuffer();
   const auto & tag = comm_desc.getTag();
 
+  auto nb_data = comm_desc.getNbData();
+  auto proc = comm_desc.getProc();
+  auto rank = mesh.getCommunicator().whoAmI();
+
+  decltype(nb_data) recv_nb_data;
+  decltype(proc) recv_proc;
+  decltype(rank) recv_rank;
+
   SynchronizationTag t;
   buffer >> t;
+  buffer >> recv_nb_data;
+  buffer >> recv_proc;
+  buffer >> recv_rank;
 
   AKANTU_DEBUG_ASSERT(
       t == tag, "The tag received does not correspond to the tag expected");
 
+  AKANTU_DEBUG_ASSERT(
+      nb_data == recv_nb_data,
+      "The nb_data received does not correspond to the nb_data expected");
+
+  AKANTU_DEBUG_ASSERT(UInt(recv_rank) == proc,
+                      "The rank received does not correspond to the proc");
+
+  AKANTU_DEBUG_ASSERT(recv_proc == UInt(rank),
+                      "The proc received does not correspond to the rank");
+
   auto & recv_element = comm_desc.getScheme();
   auto spatial_dimension = mesh.getSpatialDimension();
 
-  for(auto && element : recv_element) {
+  for (auto && element : recv_element) {
     Vector<Real> barycenter_loc(spatial_dimension);
     mesh.getBarycenter(element, barycenter_loc);
 
     Vector<Real> barycenter(spatial_dimension);
     buffer >> barycenter;
-    for (UInt i = 0; i < spatial_dimension; ++i) {
-      if (!Math::are_float_equal(barycenter_loc(i), barycenter(i)))
-        AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: "
-                           << element << "(barycenter[" << i
-                           << "] = " << barycenter_loc(i) << " and buffer[" << i
-                           << "] = " << barycenter(i) << ") ["
-                           << std::abs(barycenter(i) - barycenter_loc(i))
-                           << "] - tag: " << tag);
-    }
+
+    auto dist = barycenter_loc.distance(barycenter);
+    if (not Math::are_float_equal(dist, 0.))
+      AKANTU_EXCEPTION("Unpacking an unknown value for the element "
+                       << element << "(barycenter " << barycenter_loc
+                       << " != buffer " << barycenter << ") [" << dist
+                       << "] - tag: " << tag << " comm from " << proc << " to "
+                       << rank);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 } // namespace akantu
diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc
index ff2fced01..8a4759d6d 100644
--- a/src/synchronizer/master_element_info_per_processor.cc
+++ b/src/synchronizer/master_element_info_per_processor.cc
@@ -1,462 +1,461 @@
 /**
  * @file   master_element_info_per_processor.cc
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date   Fri Mar 11 14:57:13 2016
  *
  * @brief
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2011 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 "communicator.hh"
 #include "element_group.hh"
 #include "element_info_per_processor.hh"
 #include "element_synchronizer.hh"
 #include "mesh_iterators.hh"
 #include "mesh_utils.hh"
-#include "communicator.hh"
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 #include <iostream>
 #include <map>
 #include <tuple>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 MasterElementInfoPerProc::MasterElementInfoPerProc(
     ElementSynchronizer & synchronizer, UInt message_cnt, UInt root,
     ElementType type, const MeshPartition & partition)
     : ElementInfoPerProc(synchronizer, message_cnt, root, type),
       partition(partition), all_nb_local_element(nb_proc, 0),
       all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) {
   Vector<UInt> size(5);
   size(0) = (UInt)type;
 
   if (type != _not_defined) {
     nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
     nb_element = mesh.getNbElement(type);
     const auto & partition_num =
         this->partition.getPartition(this->type, _not_ghost);
     const auto & ghost_partition =
         this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
 
     for (UInt el = 0; el < nb_element; ++el) {
       this->all_nb_local_element[partition_num(el)]++;
       for (auto part = ghost_partition.begin(el);
            part != ghost_partition.end(el); ++part) {
         this->all_nb_ghost_element[*part]++;
       }
       this->all_nb_element_to_send[partition_num(el)] +=
           ghost_partition.getNbCols(el) + 1;
     }
 
     /// tag info
     std::vector<std::string> tag_names;
     this->getMeshData().getTagNames(tag_names, type);
     this->nb_tags = tag_names.size();
 
     size(4) = nb_tags;
 
     for (UInt p = 0; p < nb_proc; ++p) {
       if (p != root) {
         size(1) = this->all_nb_local_element[p];
         size(2) = this->all_nb_ghost_element[p];
         size(3) = this->all_nb_element_to_send[p];
         AKANTU_DEBUG_INFO(
             "Sending connectivities informations to proc "
             << p << " TAG("
             << Tag::genTag(this->rank, this->message_count, Tag::_SIZES)
             << ")");
         comm.send(size, p,
                   Tag::genTag(this->rank, this->message_count, Tag::_SIZES));
       } else {
         this->nb_local_element = this->all_nb_local_element[p];
         this->nb_ghost_element = this->all_nb_ghost_element[p];
       }
     }
   } else {
     for (UInt p = 0; p < this->nb_proc; ++p) {
       if (p != this->root) {
         AKANTU_DEBUG_INFO(
             "Sending empty connectivities informations to proc "
             << p << " TAG("
             << Tag::genTag(this->rank, this->message_count, Tag::_SIZES)
             << ")");
         comm.send(size, p,
                   Tag::genTag(this->rank, this->message_count, Tag::_SIZES));
       }
     }
   }
 }
 
 /* ------------------------------------------------------------------------ */
 void MasterElementInfoPerProc::synchronizeConnectivities() {
   const auto & partition_num =
       this->partition.getPartition(this->type, _not_ghost);
   const auto & ghost_partition =
       this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
 
   std::vector<Array<UInt>> buffers(this->nb_proc);
 
-  auto conn_it = this->mesh.getConnectivity(this->type, _not_ghost)
-                     .begin(this->nb_nodes_per_element);
-  auto conn_end = this->mesh.getConnectivity(this->type, _not_ghost)
-                      .end(this->nb_nodes_per_element);
+  const auto & connectivities =
+      this->mesh.getConnectivity(this->type, _not_ghost);
 
   /// copying the local connectivity
-  auto part_it = partition_num.begin();
-  for (; conn_it != conn_end; ++conn_it, ++part_it) {
-    const auto & conn = *conn_it;
+  for (auto && part_conn :
+       zip(partition_num,
+           make_view(connectivities, this->nb_nodes_per_element))) {
+    auto && part = std::get<0>(part_conn);
+    auto && conn = std::get<1>(part_conn);
     for (UInt i = 0; i < conn.size(); ++i) {
-      buffers[*part_it].push_back(conn[i]);
+      buffers[part].push_back(conn[i]);
     }
   }
 
   /// copying the connectivity of ghost element
-  conn_it = this->mesh.getConnectivity(this->type, _not_ghost)
-                .begin(this->nb_nodes_per_element);
-  for (UInt el = 0; conn_it != conn_end; ++el, ++conn_it) {
+  for (auto && tuple :
+           enumerate(make_view(connectivities, this->nb_nodes_per_element))) {
+    auto && el = std::get<0>(tuple);
+    auto && conn = std::get<1>(tuple);
     for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
          ++part) {
       UInt proc = *part;
-      const Vector<UInt> & conn = *conn_it;
       for (UInt i = 0; i < conn.size(); ++i) {
         buffers[proc].push_back(conn[i]);
       }
     }
   }
 
 #ifndef AKANTU_NDEBUG
-  for (UInt p = 0; p < this->nb_proc; ++p) {
+  for (auto p : arange(this->nb_proc)) {
     UInt size = this->nb_nodes_per_element *
                 (this->all_nb_local_element[p] + this->all_nb_ghost_element[p]);
     AKANTU_DEBUG_ASSERT(
         buffers[p].size() == size,
         "The connectivity data packed in the buffer are not correct");
   }
 #endif
 
   /// send all connectivity and ghost information to all processors
   std::vector<CommunicationRequest> requests;
-  for (UInt p = 0; p < this->nb_proc; ++p) {
+  for (auto p : arange(this->nb_proc)) {
     if (p == this->root)
       continue;
     auto && tag =
         Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY);
     AKANTU_DEBUG_INFO("Sending connectivities to proc " << p << " TAG(" << tag
                                                         << ")");
     requests.push_back(comm.asyncSend(buffers[p], p, tag));
   }
 
   Array<UInt> & old_nodes = this->getNodesGlobalIds();
 
   /// create the renumbered connectivity
   AKANTU_DEBUG_INFO("Renumbering local connectivities");
   MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root],
                                all_nb_ghost_element[root], type, old_nodes);
 
-
   comm.waitAll(requests);
   comm.freeCommunicationRequest(requests);
 }
 
 /* ------------------------------------------------------------------------ */
 void MasterElementInfoPerProc::synchronizePartitions() {
   const auto & partition_num =
       this->partition.getPartition(this->type, _not_ghost);
   const auto & ghost_partition =
       this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
 
   std::vector<Array<UInt>> buffers(this->partition.getNbPartition());
 
   /// splitting the partition information to send them to processors
   Vector<UInt> count_by_proc(nb_proc, 0);
   for (UInt el = 0; el < nb_element; ++el) {
     UInt proc = partition_num(el);
     buffers[proc].push_back(ghost_partition.getNbCols(el));
 
     UInt i(0);
     for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
          ++part, ++i) {
       buffers[proc].push_back(*part);
     }
   }
 
   for (UInt el = 0; el < nb_element; ++el) {
     UInt i(0);
     for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el);
          ++part, ++i) {
       buffers[*part].push_back(partition_num(el));
     }
   }
 
 #ifndef AKANTU_NDEBUG
   for (UInt p = 0; p < this->nb_proc; ++p) {
     AKANTU_DEBUG_ASSERT(buffers[p].size() == (this->all_nb_ghost_element[p] +
                                               this->all_nb_element_to_send[p]),
                         "Data stored in the buffer are most probably wrong");
   }
 #endif
 
   std::vector<CommunicationRequest> requests;
   /// last data to compute the communication scheme
   for (UInt p = 0; p < this->nb_proc; ++p) {
     if (p == this->root)
       continue;
 
     auto && tag =
         Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS);
     AKANTU_DEBUG_INFO("Sending partition informations to proc " << p << " TAG("
                                                                 << tag << ")");
     requests.push_back(comm.asyncSend(buffers[p], p, tag));
   }
 
   if (Mesh::getSpatialDimension(this->type) ==
       this->mesh.getSpatialDimension()) {
     AKANTU_DEBUG_INFO("Creating communications scheme");
     this->fillCommunicationScheme(buffers[this->rank]);
   }
 
   comm.waitAll(requests);
   comm.freeCommunicationRequest(requests);
 }
 
 /* -------------------------------------------------------------------------- */
 void MasterElementInfoPerProc::synchronizeTags() {
   AKANTU_DEBUG_IN();
 
   if (this->nb_tags == 0) {
     AKANTU_DEBUG_OUT();
     return;
   }
 
   UInt mesh_data_sizes_buffer_length;
   auto & mesh_data = this->getMeshData();
 
   /// tag info
   std::vector<std::string> tag_names;
   mesh_data.getTagNames(tag_names, type);
   // Make sure the tags are sorted (or at least not in random order),
   // because they come from a map !!
   std::sort(tag_names.begin(), tag_names.end());
 
   // Sending information about the tags in mesh_data: name, data type and
   // number of components of the underlying array associated to the current
   // type
   DynamicCommunicationBuffer mesh_data_sizes_buffer;
   for (auto && tag_name : tag_names) {
     mesh_data_sizes_buffer << tag_name;
     mesh_data_sizes_buffer << mesh_data.getTypeCode(tag_name);
     mesh_data_sizes_buffer << mesh_data.getNbComponent(tag_name, type);
   }
 
   mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.size();
   AKANTU_DEBUG_INFO(
       "Broadcasting the size of the information about the mesh data tags: ("
       << mesh_data_sizes_buffer_length << ").");
   comm.broadcast(mesh_data_sizes_buffer_length, root);
   AKANTU_DEBUG_INFO(
       "Broadcasting the information about the mesh data tags, addr "
       << (void *)mesh_data_sizes_buffer.storage());
 
   if (mesh_data_sizes_buffer_length != 0)
     comm.broadcast(mesh_data_sizes_buffer, root);
 
   if (mesh_data_sizes_buffer_length != 0) {
     // Sending the actual data to each processor
     std::vector<DynamicCommunicationBuffer> buffers(nb_proc);
     // Loop over each tag for the current type
     for (auto && tag_name : tag_names) {
       // Type code of the current tag (i.e. the tag named *names_it)
       this->fillTagBuffer(buffers, tag_name);
     }
 
     std::vector<CommunicationRequest> requests;
     for (UInt p = 0; p < nb_proc; ++p) {
       if (p == root)
         continue;
 
       auto && tag =
           Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA);
       AKANTU_DEBUG_INFO("Sending " << buffers[p].size()
                                    << " bytes of mesh data to proc " << p
                                    << " TAG(" << tag << ")");
 
       requests.push_back(comm.asyncSend(buffers[p], p, tag));
     }
 
     // Loop over each tag for the current type
     for (auto && tag_name : tag_names) {
       // Reinitializing the mesh data on the master
       this->fillMeshData(buffers[root], tag_name,
                          mesh_data.getTypeCode(tag_name),
                          mesh_data.getNbComponent(tag_name, type));
     }
 
     comm.waitAll(requests);
     comm.freeCommunicationRequest(requests);
     requests.clear();
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void MasterElementInfoPerProc::fillTagBufferTemplated(
     std::vector<DynamicCommunicationBuffer> & buffers,
     const std::string & tag_name) {
   MeshData & mesh_data = this->getMeshData();
 
   const auto & data = mesh_data.getElementalDataArray<T>(tag_name, type);
   const auto & partition_num =
       this->partition.getPartition(this->type, _not_ghost);
   const auto & ghost_partition =
       this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
 
   // Not possible to use the iterator because it potentially triggers the
   // creation of complex
   // type templates (such as akantu::Vector< std::vector<Element> > which don't
   // implement the right interface
   // (e.g. operator<< in that case).
   // typename Array<T>::template const_iterator< Vector<T> > data_it  =
   // data.begin(data.getNbComponent());
   // typename Array<T>::template const_iterator< Vector<T> > data_end =
   // data.end(data.getNbComponent());
 
   const T * data_it = data.storage();
   const T * data_end = data.storage() + data.size() * data.getNbComponent();
   const UInt * part = partition_num.storage();
 
   /// copying the data, element by element
   for (; data_it != data_end; ++part) {
     for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) {
       buffers[*part] << *data_it;
     }
   }
 
   data_it = data.storage();
   /// copying the data for the ghost element
   for (UInt el(0); data_it != data_end;
        data_it += data.getNbComponent(), ++el) {
     auto it = ghost_partition.begin(el);
     auto end = ghost_partition.end(el);
     for (; it != end; ++it) {
       UInt proc = *it;
       for (UInt j(0); j < data.getNbComponent(); ++j) {
         buffers[proc] << data_it[j];
       }
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 void MasterElementInfoPerProc::fillTagBuffer(
     std::vector<DynamicCommunicationBuffer> & buffers,
     const std::string & tag_name) {
   MeshData & mesh_data = this->getMeshData();
 
 #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem)          \
   case BOOST_PP_TUPLE_ELEM(2, 0, elem): {                                      \
     this->fillTagBufferTemplated<BOOST_PP_TUPLE_ELEM(2, 1, elem)>(buffers,     \
                                                                   tag_name);   \
     break;                                                                     \
   }
 
   MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name);
   switch (data_type_code) {
     BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, ,
                           AKANTU_MESH_DATA_TYPES)
   default:
     AKANTU_DEBUG_ERROR("Could not obtain the type of tag" << tag_name << "!");
     break;
   }
 #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA
 }
 
 /* -------------------------------------------------------------------------- */
 void MasterElementInfoPerProc::synchronizeGroups() {
   AKANTU_DEBUG_IN();
 
   std::vector<DynamicCommunicationBuffer> buffers(nb_proc);
 
   using ElementToGroup = std::vector<std::vector<std::string>>;
   ElementToGroup element_to_group(nb_element);
 
   for (auto & eg : ElementGroupsIterable(mesh)) {
     const auto & name = eg.getName();
 
     for (const auto & element : eg.getElements(type, _not_ghost)) {
       element_to_group[element].push_back(name);
     }
 
     auto eit = eg.begin(type, _not_ghost);
     if (eit != eg.end(type, _not_ghost))
       const_cast<Array<UInt> &>(eg.getElements(type)).empty();
   }
 
   const auto & partition_num =
       this->partition.getPartition(this->type, _not_ghost);
   const auto & ghost_partition =
       this->partition.getGhostPartitionCSR()(this->type, _not_ghost);
 
   /// copying the data, element by element
   for (auto && pair : zip(partition_num, element_to_group)) {
     buffers[std::get<0>(pair)] << std::get<1>(pair);
   }
 
   /// copying the data for the ghost element
   for (auto && pair : enumerate(element_to_group)) {
     auto && el = std::get<0>(pair);
     auto it = ghost_partition.begin(el);
     auto end = ghost_partition.end(el);
     for (; it != end; ++it) {
       UInt proc = *it;
       buffers[proc] << std::get<1>(pair);
     }
   }
 
   std::vector<CommunicationRequest> requests;
   for (UInt p = 0; p < this->nb_proc; ++p) {
     if (p == this->rank)
       continue;
 
     auto && tag = Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP);
     AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << tag
                                                         << ")");
     requests.push_back(comm.asyncSend(buffers[p], p, tag));
   }
 
   this->fillElementGroupsFromBuffer(buffers[this->rank]);
 
   comm.waitAll(requests);
   comm.freeCommunicationRequest(requests);
   requests.clear();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 
 } // namespace akantu
diff --git a/src/synchronizer/synchronizer_impl_tmpl.hh b/src/synchronizer/synchronizer_impl_tmpl.hh
index 7f3037156..8cd7a2cdf 100644
--- a/src/synchronizer/synchronizer_impl_tmpl.hh
+++ b/src/synchronizer/synchronizer_impl_tmpl.hh
@@ -1,305 +1,305 @@
 /**
  * @file   synchronizer_impl_tmpl.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date   Wed Aug 24 13:29:47 2016
  *
  * @brief  Implementation of the SynchronizerImpl
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2011 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 "synchronizer_impl.hh"
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__
 #define __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 SynchronizerImpl<Entity>::SynchronizerImpl(const Communicator & comm, const ID & id, MemoryID memory_id)
 
     : Synchronizer(comm, id, memory_id), communications(comm) {}
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::synchronizeOnceImpl(
     DataAccessor<Entity> & data_accessor,
     const SynchronizationTag & tag) const {
   // no need to synchronize
   if (this->nb_proc == 1)
     return;
 
   using CommunicationRequests = std::vector<CommunicationRequest>;
   using CommunicationBuffers = std::map<UInt, CommunicationBuffer>;
 
   CommunicationRequests send_requests, recv_requests;
   CommunicationBuffers send_buffers, recv_buffers;
 
   auto postComm = [&](const CommunicationSendRecv & sr,
                       CommunicationBuffers & buffers,
                       CommunicationRequests & requests) -> void {
     for (auto && pair : communications.iterateSchemes(sr)) {
       auto & proc = pair.first;
       const auto & scheme = pair.second;
 
       auto & buffer = buffers[proc];
       UInt buffer_size = data_accessor.getNbData(scheme, tag);
       buffer.resize(buffer_size);
 
       if (sr == _recv) {
         requests.push_back(communicator.asyncReceive(
             buffer, proc,
             Tag::genTag(this->rank, 0, Tag::_SYNCHRONIZE, this->hash_id)));
       } else {
         data_accessor.packData(buffer, scheme, tag);
         send_requests.push_back(communicator.asyncSend(
             buffer, proc,
             Tag::genTag(proc, 0, Tag::_SYNCHRONIZE, this->hash_id)));
       }
     }
   };
 
   // post the receive requests
   postComm(_recv, recv_buffers, recv_requests);
 
   // post the send data requests
   postComm(_send, send_buffers, send_requests);
 
   // treat the receive requests
   UInt request_ready;
   while ((request_ready = communicator.waitAny(recv_requests)) != UInt(-1)) {
     CommunicationRequest & req = recv_requests[request_ready];
     UInt proc = req.getSource();
 
     CommunicationBuffer & buffer = recv_buffers[proc];
     const auto & scheme = this->communications.getScheme(proc, _recv);
 
     data_accessor.unpackData(buffer, scheme, tag);
 
     req.free();
     recv_requests.erase(recv_requests.begin() + request_ready);
   }
 
   communicator.waitAll(send_requests);
   communicator.freeCommunicationRequest(send_requests);
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::asynchronousSynchronizeImpl(
     const DataAccessor<Entity> & data_accessor,
     const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   if (not this->communications.hasCommunicationSize(tag))
     this->computeBufferSize(data_accessor, tag);
 
   this->communications.incrementCounter(tag);
 
   // Posting the receive -------------------------------------------------------
   if (this->communications.hasPendingRecv(tag)) {
     AKANTU_CUSTOM_EXCEPTION_INFO(
         debug::CommunicationException(),
         "There must still be some pending receive communications."
             << " Tag is " << tag << " Cannot start new ones");
   }
 
   for (auto && comm_desc : this->communications.iterateRecv(tag)) {
     comm_desc.postRecv(this->hash_id);
   }
 
   // Posting the sends -------------------------------------------------------
   if (communications.hasPendingSend(tag)) {
     AKANTU_CUSTOM_EXCEPTION_INFO(
         debug::CommunicationException(),
         "There must be some pending sending communications."
             << " Tag is " << tag);
   }
 
   for (auto && comm_desc : this->communications.iterateSend(tag)) {
     comm_desc.resetBuffer();
 
 #ifndef AKANTU_NDEBUG
     this->packSanityCheckData(comm_desc);
 #endif
 
     comm_desc.packData(data_accessor);
     comm_desc.postSend(this->hash_id);
   }
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::waitEndSynchronizeImpl(
     DataAccessor<Entity> & data_accessor, const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
 #ifndef AKANTU_NDEBUG
   if (this->communications.begin(tag, _recv) !=
       this->communications.end(tag, _recv) &&
       !this->communications.hasPendingRecv(tag))
     AKANTU_CUSTOM_EXCEPTION_INFO(debug::CommunicationException(),
                                  "No pending communication with the tag \""
                                      << tag);
 #endif
 
   auto recv_end = this->communications.end(tag, _recv);
   decltype(recv_end) recv_it;
 
   while ((recv_it = this->communications.waitAnyRecv(tag)) != recv_end) {
     auto && comm_desc = *recv_it;
 #ifndef AKANTU_NDEBUG
     this->unpackSanityCheckData(comm_desc);
 #endif
 
     comm_desc.unpackData(data_accessor);
     comm_desc.resetBuffer();
     comm_desc.freeRequest();
   }
 
   this->communications.waitAllSend(tag);
   this->communications.freeSendRequests(tag);
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::computeAllBufferSizes(
     const DataAccessor<Entity> & data_accessor) {
   for (auto && tag : this->communications.iterateTags()) {
     this->computeBufferSize(data_accessor, tag);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::computeBufferSizeImpl(
     const DataAccessor<Entity> & data_accessor,
     const SynchronizationTag & tag) {
   AKANTU_DEBUG_IN();
 
   if (not this->communications.hasCommunication(tag)) {
     this->communications.initializeCommunications(tag);
     AKANTU_DEBUG_ASSERT(communications.hasCommunication(tag) == true,
                         "Communications where not properly initialized");
   }
 
   for (auto sr : iterate_send_recv) {
     for (auto && pair : this->communications.iterateSchemes(sr)) {
       auto proc = pair.first;
-      auto & scheme = pair.second;
+      const auto & scheme = pair.second;
       UInt size = 0;
 #ifndef AKANTU_NDEBUG
       size += this->sanityCheckDataSize(scheme, tag);
 #endif
       size += data_accessor.getNbData(scheme, tag);
       AKANTU_DEBUG_INFO("I have "
                         << size << "(" << printMemorySize<char>(size) << " - "
                         << scheme.size() << " element(s)) data to "
                         << std::string(sr == _recv ? "receive from" : "send to")
                         << proc << " for tag " << tag);
 
       this->communications.setCommunicationSize(tag, proc, size, sr);
     }
   }
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename Entity>
 template <typename Pred>
 void SynchronizerImpl<Entity>::split(SynchronizerImpl<Entity> & in_synchronizer,
                                      Pred && pred) {
   AKANTU_DEBUG_IN();
 
   auto filter_list = [&](auto & list, auto & new_list) {
     auto copy = list;
     list.resize(0);
     new_list.resize(0);
 
     for (auto && entity : copy) {
       if (std::forward<Pred>(pred)(entity)) {
         new_list.push_back(entity);
       } else {
         list.push_back(entity);
       }
     }
   };
 
   for (auto sr : iterate_send_recv) {
     for (auto & scheme_pair :
          in_synchronizer.communications.iterateSchemes(sr)) {
       auto proc = scheme_pair.first;
       auto & scheme = scheme_pair.second;
       auto & new_scheme = communications.createScheme(proc, sr);
       filter_list(new_scheme, scheme);
     }
   }
 
   in_synchronizer.communications.invalidateSizes();
   communications.invalidateSizes();
 
   AKANTU_DEBUG_OUT();
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename Entity>
 template <typename Updater>
 void SynchronizerImpl<Entity>::updateSchemes(Updater && scheme_updater) {
   for (auto sr : iterate_send_recv) {
     for (auto & scheme_pair : communications.iterateSchemes(sr)) {
       auto proc = scheme_pair.first;
       auto & scheme = scheme_pair.second;
       std::forward<Updater>(scheme_updater)(scheme, proc, sr);
     }
   }
 
   communications.invalidateSizes();
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 UInt SynchronizerImpl<Entity>::sanityCheckDataSize(
     const Array<Entity> &, const SynchronizationTag &) const {
   return 0;
 }
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::packSanityCheckData(
     CommunicationDescriptor<Entity> &) const {}
 
 /* -------------------------------------------------------------------------- */
 template <class Entity>
 void SynchronizerImpl<Entity>::unpackSanityCheckData(
     CommunicationDescriptor<Entity> &) const {}
 
 /* -------------------------------------------------------------------------- */
 } // namespace akantu
 
 #endif /* __AKANTU_SYNCHRONIZER_IMPL_TMPL_HH__ */