diff --git a/src/common/aka_array.hh b/src/common/aka_array.hh
index 812266cd7..3b94a7d22 100644
--- a/src/common/aka_array.hh
+++ b/src/common/aka_array.hh
@@ -1,415 +1,417 @@
 /**
  * @file   aka_array.hh
  *
  * @author Till Junge <till.junge@epfl.ch>
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Sun Nov 22 2020
  *
  * @brief  Array container for Akantu This container differs from the
  * std::vector from the fact it as 2 dimensions a main dimension and the size
  * stored per entries
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  *
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_common.hh"
 #include "aka_types.hh"
 #include "aka_view_iterators.hh"
 /* -------------------------------------------------------------------------- */
 #include <typeinfo>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #ifndef AKANTU_ARRAY_HH_
 #define AKANTU_ARRAY_HH_
 
 namespace akantu {
 
 /// class that afford to store vectors in static memory
 // NOLINTNEXTLINE(cppcoreguidelines-special-member-functions)
 class ArrayBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   explicit ArrayBase(const ID & id = "") : id(id) {}
   ArrayBase(const ArrayBase & other, const ID & id = "") {
     this->id = (id.empty()) ? other.id : id;
   }
 
   ArrayBase(ArrayBase && other) = default;
   ArrayBase & operator=(const ArrayBase & other) = default;
   ArrayBase & operator=(ArrayBase && other) noexcept = default;
 
   virtual ~ArrayBase() = default;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// get the amount of space allocated in bytes
   virtual UInt getMemorySize() const = 0;
 
   // changed empty to match std::vector empty
   inline bool empty() const __attribute__((warn_unused_result)) {
     return size_ == 0;
   }
 
   /// function to print the containt of the class
   virtual void printself(std::ostream & stream, int indent = 0) const = 0;
 
   /* ------------------------------------------------------------------------ */
   /* Accessors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Get the Size of the Array
   UInt size() const { return size_; }
   /// Get the number of components
   AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
   /// Get the name of th array
   AKANTU_GET_MACRO(ID, id, const ID &);
   /// Set the name of th array
   AKANTU_SET_MACRO(ID, id, const ID &);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
   /// id of the vector
   ID id;
 
   /// the size used
   UInt size_{0};
 
   /// number of components
   UInt nb_component{1};
 };
 
 /* -------------------------------------------------------------------------- */
 /* Memory handling layer                                                      */
 /* -------------------------------------------------------------------------- */
 enum class ArrayAllocationType {
   _default,
   _pod,
 };
 
 template <typename T>
 struct ArrayAllocationTrait
     : public std::conditional_t<
           std::is_scalar<T>::value,
           std::integral_constant<ArrayAllocationType,
                                  ArrayAllocationType::_pod>,
           std::integral_constant<ArrayAllocationType,
                                  ArrayAllocationType::_default>> {};
 
 /* -------------------------------------------------------------------------- */
 template <typename T,
           ArrayAllocationType allocation_trait = ArrayAllocationTrait<T>::value>
 class ArrayDataLayer : public ArrayBase {
 public:
   using value_type = T;
+  using size_type = typename ArrayBase::size_type;
   using reference = value_type &;
   using pointer_type = value_type *;
   using const_reference = const value_type &;
 
 public:
   ~ArrayDataLayer() override = default;
 
   /// Allocation of a new vector
   explicit ArrayDataLayer(UInt size = 0, UInt nb_component = 1,
                           const ID & id = "");
 
   /// Allocation of a new vector with a default value
   ArrayDataLayer(UInt size, UInt nb_component, const_reference value,
                  const ID & id = "");
 
   /// Copy constructor (deep copy)
   ArrayDataLayer(const ArrayDataLayer & vect, const ID & id = "");
 
   /// Copy constructor (deep copy)
   explicit ArrayDataLayer(const std::vector<value_type> & vect);
 
   // copy operator
   ArrayDataLayer & operator=(const ArrayDataLayer & other);
 
   // move constructor
   ArrayDataLayer(ArrayDataLayer && other) noexcept = default;
 
   // move assign
   ArrayDataLayer & operator=(ArrayDataLayer && other) noexcept = default;
 
 protected:
   // deallocate the memory
   virtual void deallocate() {}
 
   // allocate the memory
   virtual void allocate(UInt size, UInt nb_component);
 
   // allocate and initialize the memory
   virtual void allocate(UInt size, UInt nb_component, const T & value);
 
 public:
   /// append a tuple of size nb_component containing value
   inline void push_back(const_reference value);
   /// append a vector
   // inline void push_back(const value_type new_elem[]);
 
   /// append a Vector or a Matrix
   template <template <typename> class C,
             typename = std::enable_if_t<aka::is_tensor<C<T>>::value>>
   inline void push_back(const C<T> & new_elem);
 
   /// changes the allocated size but not the size, if new_size = 0, the size is
   /// set to min(current_size and reserve size)
   virtual void reserve(UInt size, UInt new_size = UInt(-1));
 
   /// change the size of the Array
   virtual void resize(UInt size);
 
   /// change the size of the Array and initialize the values
   virtual void resize(UInt size, const T & val);
 
   /// get the amount of space allocated in bytes
   inline UInt getMemorySize() const override;
 
   /// Get the real size allocated in memory
   inline UInt getAllocatedSize() const;
 
   /// give the address of the memory allocated for this vector
   [[deprecated("use data instead to be stl compatible")]] T * storage() const {
     return values;
   };
 
   T * data() const { return values; };
 
 protected:
   /// allocation type agnostic  data access
   T * values{nullptr};
 
   /// data storage
   std::vector<T> data_storage;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Actual Array                                                               */
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal> class Array : public ArrayDataLayer<T> {
 private:
   using parent = ArrayDataLayer<T>;
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
   using value_type = typename parent::value_type;
+  using size_type = typename parent::size_type;
   using reference = typename parent::reference;
   using pointer_type = typename parent::pointer_type;
   using const_reference = typename parent::const_reference;
   using array_type = Array<T>;
 
   ~Array() override;
 
   Array() : Array(0){};
 
   /// Allocation of a new vector
   explicit Array(UInt size, UInt nb_component = 1, const ID & id = "");
 
   /// Allocation of a new vector with a default value
   explicit Array(UInt size, UInt nb_component, const_reference value,
                  const ID & id = "");
 
   /// Copy constructor
   Array(const Array & vect, const ID & id = "");
 
   /// Copy constructor (deep copy)
   explicit Array(const std::vector<T> & vect);
 
   // copy operator
   Array & operator=(const Array & other);
 
   // move constructor
   Array(Array && other) noexcept = default;
 
   // move assign
   Array & operator=(Array && other) noexcept = default;
 
   /* ------------------------------------------------------------------------ */
   /* Iterator                                                                 */
   /* ------------------------------------------------------------------------ */
   /// iterator for Array of nb_component = 1
   using scalar_iterator = view_iterator<T>;
   /// const_iterator for Array of nb_component = 1
   using const_scalar_iterator = const_view_iterator<T>;
 
   /// iterator returning Vectors of size n  on entries of Array with
   /// nb_component = n
   using vector_iterator = view_iterator<VectorProxy<T>>;
   /// const_iterator returning Vectors of n size on entries of Array with
   /// nb_component = n
   using const_vector_iterator = const_view_iterator<VectorProxy<const T>>;
 
   /// iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using matrix_iterator = view_iterator<MatrixProxy<T>>;
   /// const iterator returning Matrices of size (m, n) on entries of Array with
   /// nb_component = m*n
   using const_matrix_iterator = const_view_iterator<MatrixProxy<T>>;
 
   /// iterator returning Tensor3 of size (m, n, k) on entries of Array with
   /// nb_component = m*n*k
   using tensor3_iterator = view_iterator<Tensor3Proxy<T>>;
   /// const iterator returning Tensor3 of size (m, n, k) on entries of Array
   /// with nb_component = m*n*k
   using const_tensor3_iterator = const_view_iterator<Tensor3Proxy<T>>;
 
   /* ------------------------------------------------------------------------ */
   template <typename... Ns> inline decltype(auto) begin(Ns &&... n);
   template <typename... Ns> inline decltype(auto) end(Ns &&... n);
 
   template <typename... Ns> inline decltype(auto) begin(Ns &&... n) const;
   template <typename... Ns> inline decltype(auto) end(Ns &&... n) const;
 
   template <typename... Ns> inline decltype(auto) begin_reinterpret(Ns &&... n);
   template <typename... Ns> inline decltype(auto) end_reinterpret(Ns &&... n);
 
   template <typename... Ns>
   inline decltype(auto) begin_reinterpret(Ns &&... n) const;
   template <typename... Ns>
   inline decltype(auto) end_reinterpret(Ns &&... n) const;
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   /// search elem in the vector, return  the position of the first occurrence or
   /// -1 if not found
   UInt find(const_reference elem) const;
 
   /// @see Array::find(const_reference elem) const
   //  UInt find(T elem[]) const;
 
   /// append a value to the end of the Array
   inline void push_back(const_reference value) { parent::push_back(value); }
 
   /// append a Vector or a Matrix
   template <typename Derived>
   inline void push_back(const Eigen::MatrixBase<Derived> & new_elem) {
     parent::push_back(new_elem);
   }
 
   /// append the content of the iterator at the end of the Array
   template <typename Ret> inline void push_back(const view_iterator<Ret> & it) {
     push_back(*it);
   }
 
   /// erase the value at position i
   inline void erase(UInt i);
   /// ask Nico, clarify
   template <typename R> inline auto erase(const view_iterator<R> & it);
 
   /// @see Array::find(const_reference elem) const
   template <template <typename> class C,
             typename = std::enable_if_t<aka::is_tensor<C<T>>::value>>
   inline UInt find(const C<T> & elem);
 
   /// set all entries of the array to the value t
   /// @param t value to fill the array with
   inline void set(T t) {
     std::fill_n(this->values, this->size_ * this->nb_component, t);
   }
 
   /// set the array to T{}
   inline void zero() { this->set({}); }
 
   /// resize the array to 0
   inline void clear() { this->resize(0); }
 
   /// set all tuples of the array to a given vector or matrix
   /// @param vm Matrix or Vector to fill the array with
   template <template <typename> class C,
             typename = std::enable_if_t<aka::is_tensor<C<T>>::value>>
   inline void set(const C<T> & vm);
 
   /// Append the content of the other array to the current one
   void append(const Array<T> & other);
 
   /// copy another Array in the current Array, the no_sanity_check allows you to
   /// force the copy in cases where you know what you do with two non matching
   /// Arrays in terms of n
   void copy(const Array<T, is_scal> & other, bool no_sanity_check = false);
 
   /// function to print the containt of the class
   void printself(std::ostream & stream, int indent = 0) const override;
 
   /* ------------------------------------------------------------------------ */
   /* Operators                                                                */
   /* ------------------------------------------------------------------------ */
 public:
   /// substraction entry-wise
   Array<T, is_scal> & operator-=(const Array<T, is_scal> & other);
   /// addition entry-wise
   Array<T, is_scal> & operator+=(const Array<T, is_scal> & other);
   /// multiply evry entry by alpha
   Array<T, is_scal> & operator*=(const T & alpha);
 
   /// check if the array are identical entry-wise
   bool operator==(const Array<T, is_scal> & other) const;
   /// @see Array::operator==(const Array<T, is_scal> & other) const
   bool operator!=(const Array<T, is_scal> & other) const;
 
   /// return a reference to the j-th entry of the i-th tuple
   inline reference operator()(UInt i, UInt j = 0);
   /// return a const reference to the j-th entry of the i-th tuple
   inline const_reference operator()(UInt i, UInt j = 0) const;
 
   /// return a reference to the ith component of the 1D array
   inline reference operator[](UInt i);
   /// return a const reference to the ith component of the 1D array
   inline const_reference operator[](UInt i) const;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions Array<T, is_scal>                                         */
 /* -------------------------------------------------------------------------- */
 template <typename T, bool is_scal>
 inline std::ostream & operator<<(std::ostream & stream,
                                  const Array<T, is_scal> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Inline Functions ArrayBase                                                 */
 /* -------------------------------------------------------------------------- */
 inline std::ostream & operator<<(std::ostream & stream,
                                  const ArrayBase & _this) {
   _this.printself(stream);
   return stream;
 }
 
 } // namespace akantu
 
 #include "aka_array_tmpl.hh"
 
 #endif /* AKANTU_ARRAY_HH_ */
diff --git a/src/common/aka_types.hh b/src/common/aka_types.hh
index 8a733f9dc..03ead8ce5 100644
--- a/src/common/aka_types.hh
+++ b/src/common/aka_types.hh
@@ -1,700 +1,799 @@
 /**
  * @file   aka_types.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Thu Feb 17 2011
  * @date last modification: Wed Dec 09 2020
  *
  * @brief  description of the "simple" types
  *
  *
  * @section LICENSE
  *
  * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * Akantu is free software: you can redistribute it and/or modify it under the
  * terms of the GNU Lesser General Public License as published by the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  * 
  * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
  * details.
  * 
  * You should have received a copy of the GNU Lesser General Public License
  * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_compatibilty_with_cpp_standard.hh"
 #include "aka_error.hh"
 #include "aka_fwd.hh"
 #include "aka_math.hh"
 /* -------------------------------------------------------------------------- */
 #include <algorithm>
 #include <array>
 #include <initializer_list>
 #include <numeric>
 #include <type_traits>
 
 #ifndef __AKANTU_AKA_TYPES_HH__
 #define __AKANTU_AKA_TYPES_HH__
 
 /* -------------------------------------------------------------------------- */
 namespace aka {
 template <typename T> struct is_eigen_map : public std::false_type {};
 } // namespace aka
 
 /* -------------------------------------------------------------------------- */
 #define EIGEN_MATRIXBASE_PLUGIN "aka_types_eigen_matrix_base_plugin.hh"
 #define EIGEN_MATRIX_PLUGIN "aka_types_eigen_matrix_plugin.hh"
 #define EIGEN_PLAINOBJECTBASE_PLUGIN                                           \
   "aka_types_eigen_plain_object_base_plugin.hh"
 #include <Eigen/Dense>
 #include <Eigen/Eigenvalues>
 /* -------------------------------------------------------------------------- */
 
 namespace akantu {
 
 using Eigen::Ref;
 
 template <typename T, Eigen::Index n = Eigen::Dynamic>
 using Vector = Eigen::Matrix<T, n, 1>;
 
 template <typename T, Eigen::Index m = Eigen::Dynamic,
           Eigen::Index n = Eigen::Dynamic>
 using Matrix = Eigen::Matrix<T, m, n, Eigen::AutoAlign | Eigen::ColMajor>;
 
 template <typename T, Eigen::Index n = Eigen::Dynamic>
 using VectorProxy =
     Eigen::Map<std::conditional_t<std::is_const<T>::value,
                                   const Vector<std::remove_const_t<T>, n>,
                                   Vector<std::remove_const_t<T>, n>>>;
 
 template <typename T, Eigen::Index m = Eigen::Dynamic,
           Eigen::Index n = Eigen::Dynamic>
 using MatrixProxy =
     Eigen::Map<std::conditional_t<std::is_const<T>::value,
                                   const Matrix<std::remove_const_t<T>, m, n>,
                                   Matrix<std::remove_const_t<T>, m, n>>>;
 
 using VectorXr = Vector<Real>;
 using MatrixXr = Matrix<Real>;
 
 enum NormType : int8_t { L_1 = 1, L_2 = 2, L_inf = -1 };
 
 struct TensorTraitBase {};
 template <size_t n> struct TensorTrait : public TensorTraitBase {};
 
 } // namespace akantu
 
 namespace aka {
 template <typename Derived>
 using is_matrice = aka::bool_constant<not Derived::IsVectorAtCompileTime>;
 
 template <typename Derived>
 using is_vector = aka::bool_constant<Derived::IsVectorAtCompileTime>;
 
 template <typename... Ds>
 using are_vectors = aka::conjunction<is_vector<Ds>...>;
 
 template <typename... Ds>
 using are_matrices = aka::conjunction<is_matrice<Ds>...>;
 
 template <typename... Ds>
 using enable_if_matrices_t = std::enable_if_t<are_matrices<Ds...>::value>;
 
 // template <typename T> struct is_eigen_map : public std::false_type {};
 
 template <typename PlainObjectType, int MapOptions, typename StrideType>
 struct is_eigen_map<Eigen::Map<PlainObjectType, MapOptions, StrideType>>
     : public std::true_type {};
 } // namespace aka
 
 /* -------------------------------------------------------------------------- */
 namespace aka {
 template <typename T>
 using is_tensor = std::is_base_of<akantu::TensorTraitBase, T>;
 
 template <typename T, size_t n>
 using is_tensor_n = std::is_base_of<akantu::TensorTrait<n>, T>;
 
 //   template <typename T> using is_vector = is_tensor_n<T, 1>;
 
 //   template <typename T> using is_matrix = is_tensor_n<T, 2>;
 
 template <std::size_t n, typename T = void, typename... Ts>
 using enable_if_tensors_n = std::enable_if<
     aka::conjunction<
         is_tensor_n<Ts, n>...,
         std::is_same<
             std::common_type_t<std::decay_t<typename Ts::value_type>...>,
             std::decay_t<typename Ts::value_type>>...>::value,
     T>;
 
 template <std::size_t n, typename T = void, typename... Ts>
 using enable_if_tensors_n_t = typename enable_if_tensors_n<n, T, Ts...>::type;
 
 //   template <typename T = void, typename... Ms>
 //   using enable_if_vectors_t = typename enable_if_tensors_n<1, T,
 //   Ms...>::type; template <typename T = void, typename... Ms>
 
 //   using enable_if_matricies_t = typename enable_if_tensors_n<2, T,
 //   Ms...>::type;
 } // namespace aka
 
 namespace akantu {
 template <typename T, std::size_t ndim> class TensorBase;
 template <typename T, size_t ndim> class TensorProxy;
 template <typename T, size_t ndim> class Tensor;
 } // namespace akantu
 
 #include "aka_view_iterators.hh"
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 template <typename T, std::size_t ndim>
 class TensorBase : public TensorTrait<ndim> {
   using RetType = TensorBase<T, ndim>;
   static_assert(ndim > 2, "TensorBase cannot by used for dimensions < 3");
 
 protected:
   using idx_t = std::size_t;
 
   template <typename... Args>
   using valid_args_t = typename std::enable_if<
       aka::conjunction<aka::disjunction<std::is_integral<Args>,
                                         std::is_enum<Args>>...>::value and
           ndim == sizeof...(Args),
       int>::type;
 
 public:
   using proxy = TensorBase<T, ndim>;
 
   template <size_t _ndim = ndim,
             std::enable_if_t<_ndim == 1 or _ndim == 2, int> = 0>
   TensorBase() {
     n.fill(0);
   }
 
   TensorBase() { n.fill(0); }
 
   template <typename... Args, valid_args_t<Args...> = 0>
   constexpr TensorBase(Args... args)
       : n{idx_t(args)...}, _size(detail::product_all(args...)) {}
 
   constexpr TensorBase(const TensorBase & other)
       : n(other.n), _size(other._size), values(other.values) {}
 
   constexpr TensorBase(TensorBase && other)
       : n(std::move(other.n)), _size(std::exchange(other._size, 0)),
         values(std::exchange(other.values, nullptr)) {}
 
 protected:
   template <typename Array, std::size_t... I>
   constexpr auto check_indices(const Array & idx,
                                std::index_sequence<I...>) const {
     bool result = true;
     (void)std::initializer_list<int>{(result &= idx[I] < n[I], 0)...};
     return result;
   }
 
   template <typename... Args> constexpr auto compute_index(Args... args) const {
     std::array<idx_t, sizeof...(Args)> idx{idx_t(args)...};
     AKANTU_DEBUG_ASSERT(
         check_indices(idx, std::make_index_sequence<sizeof...(Args)>{}),
         "The indexes are out of bound");
 
     idx_t index = 0, i = (sizeof...(Args)) - 1;
     for (; i > 0; i--) {
       index += idx[i];
       if (i > 0)
         index *= n[i - 1];
     }
     return index + idx[0];
   }
 
   template <typename S, std::size_t... I>
   constexpr auto get_slice(idx_t s, std::index_sequence<I...>) {
     return S(this->values + s * detail::product_all(n[I]...), n[I]...);
   }
 
   template <typename S, std::size_t... I>
   constexpr auto get_slice(idx_t s, std::index_sequence<I...>) const {
     return S(this->values + s * detail::product_all(n[I]...), n[I]...);
   }
 
 public:
   template <typename... Args, valid_args_t<Args...> = 0>
   inline T & operator()(Args... args) {
     return *(this->values + compute_index(std::move(args)...));
   }
 
   template <typename... Args, valid_args_t<Args...> = 0>
   inline const T & operator()(Args... args) const {
     return *(this->values + compute_index(std::move(args)...));
   }
 
   template <idx_t _ndim = ndim, std::enable_if_t<(_ndim > 3), int> = 0>
   inline auto operator()(idx_t s) {
     return get_slice<TensorProxy<T, ndim - 1>>(
         std::move(s), std::make_index_sequence<ndim - 1>());
   }
 
   template <idx_t _ndim = ndim, std::enable_if_t<(_ndim > 3), int> = 0>
   inline auto operator()(idx_t s) const {
     return get_slice<TensorProxy<T, ndim - 1>>(
         std::move(s), std::make_index_sequence<ndim - 1>());
   }
 
   template <idx_t _ndim = ndim, std::enable_if_t<_ndim == 3, int> = 0>
   inline auto operator()(idx_t s) {
     return get_slice<
         Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>>(
         std::move(s), std::make_index_sequence<ndim - 1>());
   }
 
   template <idx_t _ndim = ndim, std::enable_if_t<_ndim == 3, int> = 0>
   inline auto operator()(idx_t s) const {
     return get_slice<
         Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>>(
         std::move(s), std::make_index_sequence<ndim - 1>());
   }
 
 protected:
   template <class Operator> RetType & transform(Operator && op) {
     std::transform(this->values, this->values + this->_size, this->values,
                    std::forward<Operator>(op));
     return *(static_cast<RetType *>(this));
   }
 
   template <class Other, class Operator>
   RetType & transform(Other && other, Operator && op) {
     AKANTU_DEBUG_ASSERT(_size == other.size(),
                         "The two tensors do not have the same size "
                             << this->_size << " != " << other._size);
 
     std::transform(this->values, this->values + this->_size, other.values,
                    this->values, std::forward<Operator>(op));
     return *(static_cast<RetType *>(this));
   }
 
   template <class Operator> T accumulate(T init, Operator && op) {
     return std::accumulate(this->values, this->values + this->_size,
                            std::move(init), std::forward<Operator>(op));
   }
 
   template <class Other, class Init, class Accumulate, class Operator>
   T transform_reduce(Other && other, T init, Accumulate && acc,
                      Operator && op) {
     return std::inner_product(
         this->values, this->values + this->_size, other.data(), std::move(init),
         std::forward<Accumulate>(acc), std::forward<Operator>(op));
   }
 
   // arithmetic operators ------------------------------------------------------
   /* ------------------------------------------------------------------------ */
 public:
   inline decltype(auto) operator+=(const TensorBase & other) {
     return transform(other, [](auto && a, auto && b) { return a + b; });
   }
 
   /* ------------------------------------------------------------------------ */
   inline TensorBase & operator-=(const TensorBase & other) {
     return transform(other, [](auto && a, auto && b) { return a - b; });
   }
 
   /* ------------------------------------------------------------------------ */
   inline TensorBase & operator+=(const T & x) {
     return transform([&x](auto && a) { return a + x; });
   }
 
   /* ------------------------------------------------------------------------ */
   inline TensorBase & operator-=(const T & x) {
     return transform([&x](auto && a) { return a - x; });
   }
 
   /* ------------------------------------------------------------------------ */
   inline TensorBase & operator*=(const T & x) {
     return transform([&x](auto && a) { return a * x; });
   }
 
   /* ---------------------------------------------------------------------- */
   inline TensorBase & operator/=(const T & x) {
     return transform([&x](auto && a) { return a / x; });
   }
 
   /// Y = \alpha X + Y
   inline TensorBase & aXplusY(const TensorBase & other, const T alpha = 1.) {
     return transform(other,
                      [&alpha](auto && a, auto && b) { return alpha * a + b; });
   }
 
   /* ------------------------------------------------------------------------ */
 
   template <typename TO, std::size_t ndim_o, bool is_proxy_o>
   friend class TensorBase;
 
   T * data() { return values; }
   const T * data() const { return values; }
 
   // clang-format off
   [[deprecated("use data instead to be stl compatible")]]
   T * storage() {
     return values;
   }
 
   [[deprecated("use data instead to be stl compatible")]]
   const T * storage() const {
     return values;
   }
   // clang-format on
 
   auto size() const { return _size; }
   auto size(idx_t i) const {
     AKANTU_DEBUG_ASSERT(i < ndim, "This tensor has only " << ndim
                                                           << " dimensions, not "
                                                           << (i + 1));
     return n[i];
   };
 
   inline void set(const T & t) { std::fill_n(values, _size, t); };
   inline void clear() { set(T()); };
 
 public:
   /// "Entrywise" norm norm<L_p> @f[ \|\boldsymbol{T}\|_p = \left(
   /// \sum_i^{n[0]}\sum_j^{n[1]}\sum_k^{n[2]} |T_{ijk}|^p \right)^{\frac{1}{p}}
   /// @f]
   template <NormType norm_type>
   auto norm() -> std::enable_if_t<norm_type == L_inf, T> const {
     return accumulate(
         T(), [](auto && init, auto && a) { return init + std::abs(a); });
   }
 
   template <NormType norm_type>
   auto norm() -> std::enable_if_t<norm_type == L_1, T> const {
     return accumulate(T(), [](auto && init, auto && a) {
       return std::max(init, std::abs(a));
     });
   }
 
   template <NormType norm_type>
   auto norm() -> std::enable_if_t<norm_type == L_2, T> const {
     return std::sqrt(
         accumulate(T(), [](auto && init, auto && a) { return init + a * a; }));
   }
 
   template <NormType norm_type>
   auto norm() -> std::enable_if_t<(norm_type > L_2), T> const {
     return std::pow(accumulate(T(),
                                [](auto && init, auto && a) {
                                  return init + std::pow(a, norm_type);
                                }),
                     1. / norm_type);
   }
 
   T norm() const { return norm<L_2>(); }
 
 protected:
   template <std::size_t N, typename... Args,
             std::enable_if_t<(sizeof...(Args) == ndim), int> = 0>
   void serialize(std::ostream & stream, Args... args) const {
     stream << this->operator()(std::move(args)...);
   }
 
   template <std::size_t N, typename... Args,
             std::enable_if_t<(sizeof...(Args) < ndim), int> = 0>
   void serialize(std::ostream & stream, Args... args) const {
     stream << "[";
     for (idx_t i = 0; i < n[N]; ++i) {
       if (i != 0)
         stream << ",";
       serialize<N + 1>(stream, std::move(args)..., i);
     }
     stream << "]";
   }
 
 public:
   void printself(std::ostream & stream) const { serialize<0>(stream); };
 
 protected:
   template <std::size_t... I>
   constexpr decltype(auto) begin(std::index_sequence<I...>) {
     return view_iterator<ViewIteratorHelper_t<sizeof...(I), T>>(values,
                                                                 n[I]...);
   }
 
   template <std::size_t... I>
   constexpr decltype(auto) end(std::index_sequence<I...>) {
     return view_iterator<ViewIteratorHelper_t<sizeof...(I), T>>(values + _size,
                                                                 n[I]...);
   }
 
 public:
   decltype(auto) begin() { return begin(std::make_index_sequence<ndim - 1>{}); }
   decltype(auto) end() { return end(std::make_index_sequence<ndim - 1>{}); }
 
 protected:
   // size per dimension
   std::array<idx_t, ndim> n;
 
   // total storage size
   idx_t _size{0};
 
   // actual data location
   T * values{nullptr};
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T, size_t ndim>
 class TensorProxy : public TensorBase<T, ndim> {
 private:
   using parent = TensorBase<T, ndim>;
 
 public:
   // proxy constructor
   template <typename... Args>
   constexpr TensorProxy(T * data, Args... args) : parent(args...) {
     this->values = data;
   }
 
   constexpr TensorProxy(const TensorProxy<T, ndim> & other) : parent(other) {
     this->values = other.values;
   }
 
   constexpr TensorProxy(const Tensor<T, ndim> & other) : parent(other) {
     this->values = other.values;
   }
 
   // move constructors ---------------------------------------------------------
   // proxy -> proxy
   TensorProxy(TensorProxy && other) : parent(other) {}
 
   TensorProxy & operator=(const TensorBase<T, ndim> & other) {
     AKANTU_DEBUG_ASSERT(
         other.size() == this->size(),
         "You are trying to copy too a tensors proxy with the wrong size "
             << this->_size << " != " << other._size);
 
     static_assert(std::is_trivially_copyable<T>{},
                   "Cannot copy a tensor on non trivial types");
 
     std::copy(other.values, other.values + this->_size, this->values);
     return *this;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T, size_t ndim> class Tensor : public TensorBase<T, ndim> {
 private:
   using parent = TensorBase<T, ndim>;
 
 public:
   template <typename... Args> constexpr Tensor(Args... args) : parent(args...) {
     static_assert(
         std::is_trivially_constructible<T>{},
         "Cannot create a tensor on non trivially constructible types");
     this->values = new T[this->_size];
   }
 
   /* ------------------------------------------------------------------------ */
   virtual ~Tensor() { delete[] this->values; }
 
   // copy constructors ---------------------------------------------------------
   constexpr Tensor(const Tensor & other) : parent(other) {
     this->values = new T[this->_size];
     std::copy(other.values, other.values + this->_size, this->values);
   }
 
   constexpr explicit Tensor(const TensorProxy<T, ndim> & other)
       : parent(other) {
     //    static_assert(false, "Copying data are you sure");
     this->values = new T[this->_size];
     std::copy(other.values, other.values + this->_size, this->values);
   }
 
   // move constructors ---------------------------------------------------------
   // proxy -> proxy, non proxy -> non proxy
   Tensor(Tensor && other) : parent(other) {}
 
   // copy operator -------------------------------------------------------------
   /// operator= copy-and-swap
   Tensor & operator=(const TensorBase<T, ndim> & other) {
     if (&other == this)
       return *this;
 
     std::cout << "Warning: operator= delete data" << std::endl;
     delete[] this->values;
     this->n = other.n;
     this->_size = other._size;
 
     static_assert(
         std::is_trivially_constructible<T>{},
         "Cannot create a tensor on non trivially constructible types");
 
     this->values = new T[this->_size];
 
     static_assert(std::is_trivially_copyable<T>{},
                   "Cannot copy a tensor on non trivial types");
 
     std::copy(other.values, other.values + this->_size, this->values);
 
     return *this;
   }
 };
 
 /* -------------------------------------------------------------------------- */
 template <typename T> using Tensor3 = Tensor<T, 3>;
 template <typename T> using Tensor3Proxy = TensorProxy<T, 3>;
 template <typename T> using Tensor3Base = TensorBase<T, 3>;
+
+class ArrayBase;
+
 /* -------------------------------------------------------------------------- */
 namespace details {
   template <typename T> struct MapPlainObjectType { using type = T; };
 
   template <typename PlainObjectType, int MapOptions, typename StrideType>
   struct MapPlainObjectType<
       Eigen::Map<PlainObjectType, MapOptions, StrideType>> {
     using type = PlainObjectType;
   };
 
   template <typename T>
   using MapPlainObjectType_t = typename MapPlainObjectType<T>::type;
 
+  template <typename Scalar, Idx...> struct EigenMatrixViewHelper {};
+
+  template <typename Scalar, Idx RowsAtCompileTime>
+  struct EigenMatrixViewHelper<Scalar, RowsAtCompileTime> {
+    using type = Eigen::Matrix<Scalar, RowsAtCompileTime, 1>;
+  };
+
+  template <typename Scalar, Idx RowsAtCompileTime, Idx ColsAtCompileTime>
+  struct EigenMatrixViewHelper<Scalar, RowsAtCompileTime, ColsAtCompileTime> {
+    using type = Eigen::Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime>;
+  };
+
+  template <typename Scalar, Idx... sizes>
+  using EigenMatrixViewHelper_t =
+      typename EigenMatrixViewHelper<Scalar, sizes...>::type;
+
+  template <typename Array, Idx... sizes> class EigenView {
+    static_assert(sizeof...(sizes) == 1 or sizeof...(sizes) == 2,
+                  "Eigen only supports Vector and Matrices");
+    using size_type = typename std::decay_t<Array>::size_type;
+    using value_type = typename std::decay_t<Array>::value_type;
+
+  public:
+    EigenView(Array && array, decltype(sizes)... sizes_)
+        : array(array), sizes_(sizes_...) {}
+
+    EigenView(Array && array) : array(array), sizes_(sizes...) {}
+
+    EigenView(const EigenView & other) = default;
+    EigenView(EigenView && other) = default;
+
+    template <typename T = std::remove_reference_t<
+                  decltype(*std::declval<Array>().data())>,
+              std::enable_if_t<not std::is_const<T>::value> * = nullptr>
+    decltype(auto) begin() {
+      return aka::make_from_tuple<::akantu::view_iterator<
+          Eigen::Map<EigenMatrixViewHelper_t<value_type, sizes...>>>>(
+          std::tuple_cat(std::make_tuple(array.get().data()), sizes_));
+    }
+
+    template <typename T = std::remove_reference_t<
+                  decltype(*std::declval<Array>().data())>,
+              std::enable_if_t<not std::is_const<T>::value> * = nullptr>
+    decltype(auto) end() {
+      return aka::make_from_tuple<::akantu::view_iterator<
+          Eigen::Map<EigenMatrixViewHelper_t<value_type, sizes...>>>>(
+          std::tuple_cat(std::make_tuple(array.get().data() + array_size()),
+                         sizes_));
+    }
+
+    decltype(auto) begin() const {
+      return aka::make_from_tuple<::akantu::view_iterator<
+          Eigen::Map<const EigenMatrixViewHelper_t<value_type, sizes...>>>>(
+          std::tuple_cat(std::make_tuple(array.get().data()), sizes_));
+    }
+    decltype(auto) end() const {
+      return aka::make_from_tuple<::akantu::view_iterator<
+          Eigen::Map<const EigenMatrixViewHelper_t<value_type, sizes...>>>>(
+          std::tuple_cat(std::make_tuple(array.get().data() + array_size()),
+                         sizes_));
+    }
+
+  private:
+    template <
+        class A = Array,
+        std::enable_if_t<std::is_base_of<ArrayBase, std::decay_t<A>>::value> * =
+            nullptr>
+    size_type array_size() {
+      return array.get().size() * array.get().getNbComponent();
+    }
+
+    template <class A = Array,
+              std::enable_if_t<not std::is_base_of<
+                  ArrayBase, std::decay_t<A>>::value> * = nullptr>
+    size_type array_size() {
+      return array.get().size();
+    }
+
+  private:
+    std::reference_wrapper<std::remove_reference_t<Array>> array;
+    std::tuple<decltype(sizes)...> sizes_;
+  };
+
 } // namespace details
+
+template <Idx RowsAtCompileTime, typename Array>
+decltype(auto) make_view(Array && array) {
+  return details::EigenView<Array, RowsAtCompileTime>(
+      std::forward<Array>(array), RowsAtCompileTime);
+}
+
+template <Idx RowsAtCompileTime, Idx ColsAtCompileTime, typename Array>
+decltype(auto) make_view(Array && array) {
+  return details::EigenView<Array, RowsAtCompileTime, ColsAtCompileTime>(
+      std::forward<Array>(array), RowsAtCompileTime, ColsAtCompileTime);
+}
+
 } // namespace akantu
 
 namespace Eigen {
 /* -------------------------------------------------------------------------- */
 /* Vector                                                                     */
 /* -------------------------------------------------------------------------- */
 template <typename Derived>
 template <typename ED, typename T,
           std::enable_if_t<not std::is_const<T>::value and
                            ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::begin() {
   return ::akantu::view_iterator<typename Derived::Scalar>(
       this->derived().data());
 }
 
 template <typename Derived>
 template <typename ED, typename T,
           std::enable_if_t<not std::is_const<T>::value and
                            ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::end() {
   return ::akantu::view_iterator<typename Derived::Scalar>(
       this->derived().data() + this->size());
 }
 
 template <typename Derived>
 template <typename ED, std::enable_if_t<ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::begin() const {
   using Scalar = typename Derived::Scalar;
   return ::akantu::const_view_iterator<Scalar>(this->derived().data());
 }
 
 template <typename Derived>
 template <typename ED, std::enable_if_t<ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::end() const {
   using Scalar = typename Derived::Scalar;
   return ::akantu::const_view_iterator<Scalar>(this->derived().data() +
                                                this->size());
 }
 
 /* -------------------------------------------------------------------------- */
 /* Matrix                                                                     */
 /* -------------------------------------------------------------------------- */
 template <typename Derived>
 template <typename ED, typename T,
           std::enable_if_t<not std::is_const<T>::value and
                            not ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::begin() {
   return ::akantu::view_iterator<
       Map<Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime, 1>>>(
       this->derived().data(), this->rows());
 }
 
 template <typename Derived>
 template <typename ED, typename T,
           std::enable_if_t<not std::is_const<T>::value and
                            not ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::end() {
   return ::akantu::view_iterator<
       Map<Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime, 1>>>(
       this->derived().data() + this->size(), this->rows());
 }
 
 template <typename Derived>
 template <typename ED, std::enable_if_t<not ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::begin() const {
   using Scalar = typename Derived::Scalar;
   return ::akantu::const_view_iterator<
       Map<const Matrix<Scalar, Derived::RowsAtCompileTime, 1>>>(
       const_cast<Scalar *>(this->derived().data()), this->rows());
 }
 
 template <typename Derived>
 template <typename ED, std::enable_if_t<not ED::IsVectorAtCompileTime> *>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE decltype(auto)
 MatrixBase<Derived>::end() const {
   using Scalar = typename Derived::Scalar;
   return ::akantu::const_view_iterator<
       Map<const Matrix<Scalar, Derived::RowsAtCompileTime, 1>>>(
       const_cast<Scalar *>(this->derived().data()) + this->size(),
       this->rows());
 }
 
 template <typename Derived>
 template <typename OtherDerived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void
 MatrixBase<Derived>::eig(MatrixBase<OtherDerived> & values) {
   EigenSolver<akantu::details::MapPlainObjectType_t<std::decay_t<Derived>>>
       solver(*this, false);
   values = solver.eigenvalues();
 }
 
 template <typename Derived>
 template <typename D1, typename D2>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void
 MatrixBase<Derived>::eig(MatrixBase<D1> & values, MatrixBase<D2> & vectors) {
   EigenSolver<akantu::details::MapPlainObjectType_t<std::decay_t<Derived>>>
       solver(*this, false);
   values = solver.eigenvalues();
   vectors = solver.eigenvectors();
 }
 
 template <typename Derived>
 template <typename OtherDerived>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void
 MatrixBase<Derived>::eigh(MatrixBase<OtherDerived> & values) {
   SelfAdjointEigenSolver<
       akantu::details::MapPlainObjectType_t<std::decay_t<Derived>>>
       solver(*this, false);
   values = solver.eigenvalues();
 }
 
 template <typename Derived>
 template <typename D1, typename D2>
 EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void
 MatrixBase<Derived>::eigh(MatrixBase<D1> & values, MatrixBase<D2> & vectors) {
   SelfAdjointEigenSolver<
       akantu::details::MapPlainObjectType_t<std::decay_t<Derived>>>
       solver(*this, false);
   values = solver.eigenvalues();
   vectors = solver.eigenvectors();
 }
 
 } // namespace Eigen
 
 #endif /* AKANTU_AKA_TYPES_HH_ */
diff --git a/src/common/aka_view_iterators.hh b/src/common/aka_view_iterators.hh
index 828c43e6a..94bbeb2ba 100644
--- a/src/common/aka_view_iterators.hh
+++ b/src/common/aka_view_iterators.hh
@@ -1,505 +1,515 @@
 /**
  * @file   aka_view_iterators.hh
  *
  * @author Nicolas Richart
  *
  * @date creation  Thu Nov 15 2018
  *
  * @brief View iterators
  *
  * @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_common.hh"
 #include "aka_types.hh"
 /* -------------------------------------------------------------------------- */
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 #ifndef __AKANTU_AKA_VIEW_ITERATORS_HH__
 #define __AKANTU_AKA_VIEW_ITERATORS_HH__
 
 namespace akantu {
 
 /* -------------------------------------------------------------------------- */
 /* Iterators                                                                  */
 /* -------------------------------------------------------------------------- */
 namespace detail {
   template <class R> struct IteratorHelper { static constexpr Int dim = 0; };
 
   template <class Derived> struct IteratorHelper<Eigen::MatrixBase<Derived>> {
   private:
     using T = typename Derived::Scalar;
     static constexpr Int m = Derived::RowsAtCompileTime;
     static constexpr Int n = Derived::ColsAtCompileTime;
 
   public:
     static constexpr Int dim =
         Eigen::MatrixBase<Derived>::IsVectorAtCompileTime ? 1 : 2;
     using pointer = T *;
     using proxy = Eigen::Map<Eigen::Matrix<T, m, n>>;
     using const_proxy = Eigen::Map<const Eigen::Matrix<T, m, n>>;
   };
 
   template <class Derived> struct IteratorHelper<Eigen::Map<Derived>> {
   private:
     using T = typename Derived::Scalar;
     static constexpr Int m = Derived::RowsAtCompileTime;
     static constexpr Int n = Derived::ColsAtCompileTime;
 
   public:
-    static constexpr Int dim = Derived::IsVectorAtCompileTime ? 1 : 2;
+    static constexpr Int dim =
+        Derived::IsVectorAtCompileTime and m != 1 ? 1 : 2;
     using pointer = T *;
     using proxy = Eigen::Map<Derived>;
     using const_proxy = Eigen::Map<const Derived>;
   };
 
   template <typename T, std::size_t _dim>
   struct IteratorHelper<Tensor<T, _dim>> {
     static constexpr Int dim = _dim;
     using pointer = T *;
     using proxy = TensorProxy<T, _dim>;
     using const_proxy = TensorProxy<const T, _dim>;
   };
 
   template <typename T, std::size_t _dim>
   struct IteratorHelper<TensorProxy<T, _dim>> {
     static constexpr Int dim = _dim;
     using pointer = T *;
     using proxy = TensorProxy<T, _dim>;
     using const_proxy = TensorProxy<const T, _dim>;
   };
 
   /* --------------------------------------------------------------------------
    */
   template <class R, class daughter, class IR = R,
             Int dim = IteratorHelper<std::decay_t<R>>::dim>
   class internal_view_iterator {
   protected:
     using helper = IteratorHelper<std::decay_t<R>>;
     using internal_value_type = IR;
     using internal_pointer = IR *;
     using scalar_pointer = typename helper::pointer;
     using proxy_t = typename helper::proxy;
     using const_proxy_t = typename helper::const_proxy;
 
   public:
     using pointer = proxy_t *;
     using value_type = proxy_t;
     using reference = proxy_t &;
     using const_reference = const_proxy_t &;
     using difference_type = Int;
     using iterator_category = std::random_access_iterator_tag;
 
   private:
     template <std::size_t... I>
     constexpr auto get_new_const_proxy(scalar_pointer data,
                                        std::index_sequence<I...>) const {
       return const_proxy_t(data, dims[I]...);
     }
 
     constexpr auto get_new_const_proxy(scalar_pointer data) const {
       return get_new_const_proxy(data, std::make_index_sequence<dim>());
     }
     template <std::size_t... I>
     constexpr auto get_new_proxy(scalar_pointer data,
                                  std::index_sequence<I...>) {
       return proxy_t(data, dims[I]...);
     }
 
     constexpr auto get_new_proxy(scalar_pointer data) {
       return get_new_proxy(data, std::make_index_sequence<dim>());
     }
 
     template <typename T, std::size_t... I>
     constexpr void reset_proxy(T & t, scalar_pointer data,
                                std::index_sequence<I...>) const {
       new (&t) T(data, dims[I]...);
     }
 
     template <typename T> constexpr auto reset_proxy(T & t) const {
       return reset_proxy(t, this->ret_ptr, std::make_index_sequence<dim>());
     }
 
   protected:
     template <typename OR, typename OD, typename OIR,
               std::enable_if_t<std::is_convertible<
                   decltype(std::declval<OIR>().data()),
                   decltype(std::declval<IR>().data())>::value> * = nullptr>
     explicit internal_view_iterator(
         internal_view_iterator<OR, OD, OIR, dim> & it)
         : dims(it.dims), _offset(it._offset), initial(it.initial),
           ret_ptr(it.ret_ptr), proxy(get_new_proxy(it.ret_ptr)),
           const_proxy(get_new_const_proxy(it.ret_ptr)) {}
 
     template <typename OR, typename OD, typename OIR, Int _dim>
     friend class internal_view_iterator;
 
     template <typename... Args>
-    using valid_args_t = typename std::enable_if<
+    using valid_args_t = std::enable_if_t<
         aka::conjunction<aka::disjunction<std::is_integral<Args>,
                                           std::is_enum<Args>>...>::value and
             dim == sizeof...(Args),
-        int>::type;
+        int>;
 
   public:
     template <typename... Ns, valid_args_t<Ns...> = 0>
     internal_view_iterator(scalar_pointer data, Ns... ns)
         : dims({Int(ns)...}),
           _offset(detail::product_all(std::forward<Ns>(ns)...)), initial(data),
           ret_ptr(data), proxy(data, ns...), const_proxy(data, ns...) {}
 
-    template<Int _dim = dim, std::enable_if_t<_dim == 1> * = nullptr>
+    // Static 1x1 matrix cannot be differenciated from vector in eigen
+    template <typename RD = std::decay_t<R>,
+              std::enable_if_t<RD::RowsAtCompileTime == 1 and
+                               RD::ColsAtCompileTime == 1> * = nullptr>
+    internal_view_iterator(scalar_pointer data, Idx rows)
+        : dims({rows, 1}), _offset(rows), initial(data), ret_ptr(data),
+          proxy(data, rows, 1), const_proxy(data, rows, 1) {}
+
+    template <Int _dim = dim, std::enable_if_t<_dim == 1> * = nullptr>
     internal_view_iterator() : proxy(nullptr, 0), const_proxy(nullptr, 0) {}
 
-    template<Int _dim = dim, std::enable_if_t<_dim == 2> * = nullptr>
-    internal_view_iterator() : proxy(nullptr, 0, 0), const_proxy(nullptr, 0, 0) {}
+    template <Int _dim = dim, std::enable_if_t<_dim == 2> * = nullptr>
+    internal_view_iterator()
+        : proxy(nullptr, 0, 0), const_proxy(nullptr, 0, 0) {}
 
     internal_view_iterator(const internal_view_iterator & it)
         : proxy(get_new_proxy(it.ret_ptr)),
           const_proxy(get_new_const_proxy(it.ret_ptr)) {
       if (this != &it) {
         this->dims = it.dims;
         this->_offset = it._offset;
         this->initial = it.initial;
         this->ret_ptr = it.ret_ptr;
       }
     }
 
     internal_view_iterator(internal_view_iterator && it) = default;
 
     virtual ~internal_view_iterator() = default;
 
     inline internal_view_iterator &
     operator=(const internal_view_iterator & it) {
       if (this != &it) {
         this->_offset = it._offset;
         this->initial = it.initial;
         this->ret_ptr = it.ret_ptr;
         reset_proxy(this->proxy);
         reset_proxy(this->const_proxy);
       }
       return *this;
     }
 
   public:
     UInt getCurrentIndex() {
       return (this->ret_ptr - this->initial) / this->_offset;
     }
 
     inline reference operator*() {
       reset_proxy(proxy);
       return proxy;
     }
 
     inline const_reference operator*() const {
       reset_proxy(const_proxy);
       return proxy;
     }
 
     inline pointer operator->() {
       reset_proxy(proxy);
       return &proxy;
     }
 
     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+=(UInt n) {
       ret_ptr += _offset * n;
       return static_cast<daughter &>(*this);
     }
 
     inline daughter & operator-=(UInt n) {
       ret_ptr -= _offset * n;
       return static_cast<daughter &>(*this);
     }
 
     inline auto operator[](UInt n) {
       return get_new_proxy(ret_ptr + n * _offset);
     }
 
     inline auto operator[](UInt n) const {
       return get_new_const_proxy(ret_ptr + n * _offset);
     }
 
     inline bool operator==(const internal_view_iterator & other) const {
       return this->ret_ptr == other.ret_ptr;
     }
 
     inline bool operator!=(const internal_view_iterator & other) const {
       return this->ret_ptr != other.ret_ptr;
     }
 
     inline bool operator<(const internal_view_iterator & other) const {
       return this->ret_ptr < other.ret_ptr;
     }
 
     inline bool operator<=(const internal_view_iterator & other) const {
       return this->ret_ptr <= other.ret_ptr;
     }
 
     inline bool operator>(const internal_view_iterator & other) const {
       return this->ret_ptr > other.ret_ptr;
     }
 
     inline bool operator>=(const internal_view_iterator & 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 internal_view_iterator & b) {
       return (this->ret_ptr - b.ret_ptr) / _offset;
     }
 
     inline scalar_pointer data() const { return ret_ptr; }
     inline difference_type offset() const { return _offset; }
 
   protected:
     std::array<Int, dim> dims;
     difference_type _offset{0};
     scalar_pointer initial{nullptr};
     scalar_pointer ret_ptr{nullptr};
     proxy_t proxy;
     const_proxy_t const_proxy;
   };
 
   /* ------------------------------------------------------------------------ */
   /**
    * Specialization for scalar types
    */
   template <class R, class daughter, class IR>
   class internal_view_iterator<R, daughter, IR, 0> {
   public:
     using value_type = R;
     using pointer = R *;
     using reference = R &;
     using const_reference = const R &;
     using difference_type = std::ptrdiff_t;
     using iterator_category = std::random_access_iterator_tag;
 
   protected:
     using internal_value_type = IR;
     using internal_pointer = IR *;
 
   public:
     internal_view_iterator(pointer data = nullptr) : ret(data), initial(data){};
     internal_view_iterator(const internal_view_iterator & it) = default;
     internal_view_iterator(internal_view_iterator && it) = default;
 
     virtual ~internal_view_iterator() = default;
 
     inline internal_view_iterator &
     operator=(const internal_view_iterator & 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 internal_view_iterator & other) const {
       return ret == other.ret;
     }
     inline bool operator!=(const internal_view_iterator & other) const {
       return ret != other.ret;
     }
     inline bool operator<(const internal_view_iterator & other) const {
       return ret < other.ret;
     }
     inline bool operator<=(const internal_view_iterator & other) const {
       return ret <= other.ret;
     }
     inline bool operator>(const internal_view_iterator & other) const {
       return ret > other.ret;
     }
     inline bool operator>=(const internal_view_iterator & 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 internal_view_iterator & b) {
       return ret - b.ret;
     }
 
     inline pointer data() const { return ret; }
 
   protected:
     pointer ret{nullptr};
     pointer initial{nullptr};
   };
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <typename R> class view_iterator;
 
 template <typename R>
 class const_view_iterator
     : public detail::internal_view_iterator<const R, const_view_iterator<R>,
                                             R> {
 public:
   using parent =
       detail::internal_view_iterator<const R, const_view_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_view_iterator() : parent(){};
   const_view_iterator(const const_view_iterator & it) = default;
   const_view_iterator(const_view_iterator && it) = default;
 
   template <typename P, typename... Ns>
   const_view_iterator(P * data, Ns... ns) : parent(data, ns...) {}
 
   const_view_iterator & operator=(const const_view_iterator & it) = default;
 
   template <typename OR,
             std::enable_if_t<not std::is_same<R, OR>::value> * = nullptr>
   const_view_iterator(const_view_iterator<OR> & it) : parent(it) {}
 };
 
 template <class R, bool is_tensor_ = aka::is_tensor<R>::value>
 struct ConstConverterIteratorHelper {
   static inline auto convert(const view_iterator<R> & it) {
     return const_view_iterator<R>(std::unique_ptr<R>(new R(*it, false)));
   }
 };
 
 template <class R> struct ConstConverterIteratorHelper<R, false> {
   static inline auto convert(const view_iterator<R> & it) {
     return const_view_iterator<R>(it.data());
   }
 };
 
 template <typename R>
 class view_iterator
     : public detail::internal_view_iterator<R, view_iterator<R>> {
 public:
   using parent = detail::internal_view_iterator<R, view_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:
   view_iterator() : parent(){};
   view_iterator(const view_iterator & it) = default;
   view_iterator(view_iterator && it) = default;
 
   template <typename P, typename... Ns>
   view_iterator(P * data, Ns... ns) : parent(data, ns...) {}
 
   view_iterator & operator=(const view_iterator & it) = default;
 
   operator const_view_iterator<R>() {
     return ConstConverterIteratorHelper<R>::convert(*this);
   }
 };
 
 namespace {
   template <std::size_t dim, typename T> struct ViewIteratorHelper {
     using type = TensorProxy<T, dim>;
   };
 
   template <typename T> struct ViewIteratorHelper<0, T> { using type = T; };
   template <typename T> struct ViewIteratorHelper<1, T> {
     using type = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>>;
   };
   template <typename T> struct ViewIteratorHelper<1, const T> {
     using type = Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>>;
   };
 
   template <typename T> struct ViewIteratorHelper<2, T> {
     using type = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
   };
 
   template <typename T> struct ViewIteratorHelper<2, const T> {
     using type =
         Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>;
   };
 
   template <std::size_t dim, typename T>
   using ViewIteratorHelper_t = typename ViewIteratorHelper<dim, T>::type;
 } // namespace
 
 // #include <iterator>
 
 // namespace std {
 // template <typename R>
 // struct iterator_traits<::akantu::types::details::vector_iterator<R>> {
 // protected:
 //   using iterator = ::akantu::types::details::vector_iterator<R>;
 
 // public:
 //   using iterator_category = typename iterator::iterator_category;
 //   using value_type = typename iterator::value_type;
 //   using difference_type = typename iterator::difference_type;
 //   using pointer = typename iterator::pointer;
 //   using reference = typename iterator::reference;
 // };
 // } // namespace std
 
 } // namespace akantu
 
 #endif /* !__AKANTU_AKA_VIEW_ITERATORS_HH__ */
diff --git a/third-party/akantu_iterators/include/aka_compatibilty_with_cpp_standard.hh b/third-party/akantu_iterators/include/aka_compatibilty_with_cpp_standard.hh
index 5861502d3..f3f13376b 100644
--- a/third-party/akantu_iterators/include/aka_compatibilty_with_cpp_standard.hh
+++ b/third-party/akantu_iterators/include/aka_compatibilty_with_cpp_standard.hh
@@ -1,193 +1,207 @@
 /**
  * @file   aka_compatibilty_with_cpp_standard.hh
  *
  * @author Nicolas Richart <nicolas.richart@epfl.ch>
  *
  * @date creation: Fri Jun 18 2010
  * @date last modification: Wed Jan 10 2018
  *
  * @brief  The content of this file is taken from the possible implementations
  * on
  * http://en.cppreference.com
  *
  *
  * Copyright (©)  2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
  * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  * akantu-iterators 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-iterators 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-iterators. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "aka_static_if.hh"
 /* -------------------------------------------------------------------------- */
 #include <iterator>
 #include <tuple>
 #include <type_traits>
 #include <utility>
 #if __cplusplus >= 201703L
 #include <functional>
 #endif
 /* -------------------------------------------------------------------------- */
 
 #ifndef AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH
 #define AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH
 
 namespace aka {
 
 /* -------------------------------------------------------------------------- */
 // Part taken from C++14
 #if __cplusplus < 201402L
 template <bool B, class T = void>
 using enable_if_t = typename enable_if<B, T>::type;
 #else
 template <bool B, class T = void> using enable_if_t = std::enable_if_t<B, T>;
 #endif
 
 /* -------------------------------------------------------------------------- */
 // Part taken from C++17
 #if __cplusplus < 201703L
 
 /* -------------------------------------------------------------------------- */
 // bool_constant
 template <bool B> using bool_constant = std::integral_constant<bool, B>;
 namespace {
   template <bool B> constexpr bool bool_constant_v = bool_constant<B>::value;
 }
 
 /* -------------------------------------------------------------------------- */
 // conjunction
 template <class...> struct conjunction : std::true_type {};
 template <class B1> struct conjunction<B1> : B1 {};
 template <class B1, class... Bn>
 struct conjunction<B1, Bn...>
     : std::conditional_t<bool(B1::value), conjunction<Bn...>, B1> {};
 
 /* -------------------------------------------------------------------------- */
 // disjunction
 template <class...> struct disjunction : std::false_type {};
 template <class B1> struct disjunction<B1> : B1 {};
 template <class B1, class... Bn>
 struct disjunction<B1, Bn...>
     : std::conditional_t<bool(B1::value), B1, disjunction<Bn...>> {};
 
 /* -------------------------------------------------------------------------- */
 // negations
 template <class B> struct negation : bool_constant<!bool(B::value)> {};
 
 /* -------------------------------------------------------------------------- */
 // invoke
 namespace detail {
   template <class T> struct is_reference_wrapper : std::false_type {};
   template <class U>
   struct is_reference_wrapper<std::reference_wrapper<U>> : std::true_type {};
 
   template <class T, class Type, class T1, class... Args>
   decltype(auto) INVOKE(Type T::*f, T1 && t1, Args &&... args) {
     static_if(std::is_member_function_pointer<decltype(f)>{})
         .then_([&](auto && f) {
           static_if(std::is_base_of<T, std::decay_t<T1>>{})
               .then_([&](auto && f) {
                 return (std::forward<T1>(t1).*f)(std::forward<Args>(args)...);
               })
               .elseif(is_reference_wrapper<std::decay_t<T1>>{})([&](auto && f) {
                 return (t1.get().*f)(std::forward<Args>(args)...);
               })
               .else_([&](auto && f) {
                 return ((*std::forward<T1>(t1)).*
                         f)(std::forward<Args>(args)...);
               })(std::forward<decltype(f)>(f));
         })
         .else_([&](auto && f) {
           static_assert(std::is_member_object_pointer<decltype(f)>::value,
                         "f is not a member object");
           static_assert(sizeof...(args) == 0, "f takes arguments");
           static_if(std::is_base_of<T, std::decay_t<T1>>{})
               .then_([&](auto && f) { return std::forward<T1>(t1).*f; })
               .elseif(std::is_base_of<T, std::decay_t<T1>>{})(
                   [&](auto && f) { return t1.get().*f; })
               .else_([&](auto && f) { return (*std::forward<T1>(t1)).*f; })(
                   std::forward<decltype(f)>(f));
         })(std::forward<decltype(f)>(f));
   }
 
   template <class F, class... Args>
   decltype(auto) INVOKE(F && f, Args &&... args) {
     return std::forward<F>(f)(std::forward<Args>(args)...);
   }
 } // namespace detail
 
 template <class F, class... Args>
 decltype(auto) invoke(F && f, Args &&... args) {
   return detail::INVOKE(std::forward<F>(f), std::forward<Args>(args)...);
 }
 
 /* -------------------------------------------------------------------------- */
 // apply
 namespace detail {
   template <class F, class Tuple, std::size_t... Is>
   constexpr decltype(auto) apply_impl(F && f, Tuple && t,
                                       std::index_sequence<Is...> /*unused*/) {
     return invoke(std::forward<F>(f), std::get<Is>(std::forward<Tuple>(t))...);
   }
 } // namespace detail
 
 /* -------------------------------------------------------------------------- */
 template <class F, class Tuple>
 constexpr decltype(auto) apply(F && f, Tuple && t) {
   return detail::apply_impl(
       std::forward<F>(f), std::forward<Tuple>(t),
       std::make_index_sequence<std::tuple_size<std::decay_t<Tuple>>::value>{});
 }
 
 /* -------------------------------------------------------------------------- */
 // count_if
 template <class InputIt, class UnaryPredicate>
 typename std::iterator_traits<InputIt>::difference_type
 count_if(InputIt first, InputIt last, UnaryPredicate p) {
   typename std::iterator_traits<InputIt>::difference_type ret = 0;
   for (; first != last; ++first) {
     if (p(*first)) {
       ret++;
     }
   }
   return ret;
 }
 
+namespace detail {
+  template <class T, class Tuple, std::size_t... I>
+  constexpr T make_from_tuple_impl(Tuple && t, std::index_sequence<I...>) {
+    return T(std::get<I>(std::forward<Tuple>(t))...);
+  }
+} // namespace detail
+
+template <class T, class Tuple> constexpr T make_from_tuple(Tuple && t) {
+  return detail::make_from_tuple_impl<T>(
+      std::forward<Tuple>(t),
+      std::make_index_sequence<
+          std::tuple_size<std::remove_reference_t<Tuple>>::value>{});
+}
+
 #else
 template <bool B> using bool_constant = std::bool_constant<B>;
 template <bool B> constexpr bool bool_constant_v = std::bool_constant<B>::value;
 
 template <class... Args> using conjunction = std::conjunction<Args...>;
 template <class... Args> using disjunction = std::disjunction<Args...>;
 template <class B> using negation = std::negation<B>;
 
 template <class F, class Tuple>
 constexpr decltype(auto) apply(F && f, Tuple && t) {
   return std::apply(std::forward<F>(f), std::forward<Tuple>(t));
 }
 
 template <class InputIt, class UnaryPredicate>
 decltype(auto) count_if(InputIt first, InputIt last, UnaryPredicate p) {
   return std::count_if(first, last, p);
 }
 #endif
 
 template <typename cat1, typename cat2>
 using is_iterator_category_at_least =
     std::is_same<std::common_type_t<cat1, cat2>, cat2>;
 
 } // namespace aka
 
 #endif /* AKANTU_AKA_COMPATIBILTY_WITH_CPP_STANDARD_HH */