diff --git a/src/core/grid.cpp b/src/core/grid.cpp
index f45eb26..7098849 100644
--- a/src/core/grid.cpp
+++ b/src/core/grid.cpp
@@ -1,156 +1,155 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <cstring>
 #include <complex>
 #include <algorithm>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid():
   GridBase<T>() {
   this->n.fill(0);
   this->strides.fill(1);
   this->nb_components = 1;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const std::array<UInt, dim> & n,
 		   UInt nb_components):
   GridBase<T>() {
   this->init(n, nb_components);
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const std::vector<UInt> & n,
 		   UInt nb_components):
   GridBase<T>() {
   if (n.size() != dim)
     TAMAAS_EXCEPTION("Provided sizes for grid do not match dimension");
   this->init(n, nb_components);
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const std::initializer_list<UInt> & n,
 		   UInt nb_components):
   GridBase<T>() {
   if (n.size() != dim)
     TAMAAS_EXCEPTION("Provided sizes for grid do not match dimension");
   this->init(n, nb_components);
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const Grid<T, dim> & o):
   GridBase<T>(o),
   n(o.n),
   strides(o.strides)
 {}
   
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(Grid<T, dim>&& o):
   GridBase<T>(o),
   n(std::move(o.n)),
   strides(std::move(o.strides))
 {}
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const std::array<UInt, dim> & n) {
   this->n = n;
   UInt size = this->computeSize();
   this->data.resize(size);
   this->data.assign(size, T(0.));
   this->computeStrides();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::computeStrides() {
   std::copy(n.begin()+1, n.end(), strides.rbegin()+2);
   strides[dim] = 1;
   strides[dim-1] = this->nb_components;
   std::partial_sum(strides.rbegin(), strides.rend(),
 		   strides.rbegin(), std::multiplies<UInt>());
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::printself(std::ostream & str) const {
   str << "Grid(" << dim << ", " << this->nb_components << ") {";
-  const auto size = this->dataSize();
-  for (UInt i = 0 ; i < size ; i++) {
-    str << (*this)(i) << ", ";
+  for (auto & val : *this) {
+    str << val << ", ";
   }
   str << "\b\b}";
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define GRID_SCALAR_OPERATOR_IMPL(op)			\
   template <typename T, UInt dim>			\
   inline Grid<T, dim> & Grid<T, dim>::operator op(const T & e) {	\
     const auto size = this->dataSize();					\
     _Pragma("omp parallel for")						\
     for (UInt i = 0 ; i < size ; i++) {					\
       (*this)(i) op e;							\
     }									\
     return *this;							\
   }									\
 
 GRID_SCALAR_OPERATOR_IMPL(+=);
 GRID_SCALAR_OPERATOR_IMPL(*=);
 GRID_SCALAR_OPERATOR_IMPL(-=);
 GRID_SCALAR_OPERATOR_IMPL(/=);
 GRID_SCALAR_OPERATOR_IMPL(=);
 
 #undef GRID_SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 /// Class instanciation
 #define GRID_INSTANCIATE_TYPE(type) template class Grid<type, 1>;	\
   template class Grid<type, 2>;						\
   template class Grid<type, 3>
 
 GRID_INSTANCIATE_TYPE(Real);
 GRID_INSTANCIATE_TYPE(UInt);
 GRID_INSTANCIATE_TYPE(Complex);
 GRID_INSTANCIATE_TYPE(int);
 GRID_INSTANCIATE_TYPE(bool);
 GRID_INSTANCIATE_TYPE(unsigned long);
 #undef GRID_INSTANCIATE_TYPE
 
 __END_TAMAAS__
diff --git a/src/core/grid.hh b/src/core/grid.hh
index 4e30eec..6cf7dba 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,212 +1,230 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_HH__
 #define __GRID_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid_base.hh"
 #include <array>
 #include <vector>
 #include <numeric>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Multi-dimensional & multi-component array class
  *
  * This class is a container for multi-component data stored on a multi-
  * dimensional grid.
  *
  * The access function is the parenthesis operator. For a grid of dimension d,
  * the operator takes d+1 arguments: the first d arguments are the position on
  * the grid and the last one is the component of the stored data.
  *
  * It is also possible to use the access operator with only one argument, it is
  * then considering the grid as a flat array, accessing the given cell of the
  * array.
  */
 template <typename T, UInt dim>
 class Grid: public GridBase<T> {
 public:
 /* -------------------------------------------------------------------------- */
 /* Types                                                                      */
 /* -------------------------------------------------------------------------- */
   typedef T value_type;
-  typedef typename GridBase<T>::iterator iterator;
   static constexpr UInt dimension = dim;
 
 /* -------------------------------------------------------------------------- */
 /* Constructors                                                               */
 /* -------------------------------------------------------------------------- */
 public:
   /// Constructor by default (empty array)
   Grid();
   /// Constructor
   Grid(const std::array<UInt, dim> & n,
        UInt nb_components);
   /// Constructor with vectors
   Grid(const std::vector<UInt> & n,
        UInt nb_components);
   /// Constructor with initializer list
   Grid(const std::initializer_list<UInt> & n,
        UInt nb_components);
   /// Copy constructor
   Grid(const Grid & o);
   /// Move constructor (transfers data ownership)
   Grid(Grid&& o);
 
 private:
   /// Init from standard container
   template <typename Container>
   void init(const Container & n, UInt nb_components);
 
 public:
 /* -------------------------------------------------------------------------- */
 /* Common operations                                                          */
 /* -------------------------------------------------------------------------- */
   /// Resize array
   void resize(const std::array<UInt, dim> & n);
   /// Compute size
   inline UInt computeSize() const;
   /// Compute strides
   virtual void computeStrides();
   /// Print
   virtual void printself(std::ostream & str) const;
   /// Get sizes
   const std::array<UInt, dim> & sizes() const {return n;}
   /// Get strides
   const std::array<UInt, dim+1> & getStrides() const {return this->strides;}
 
 /* -------------------------------------------------------------------------- */
 /* Access operators                                                           */
 /* -------------------------------------------------------------------------- */
   /// Variadic access operator (non-const)
   template <typename... T1>
   inline T & operator()(T1... args);
   /// Variadic access operator
   template <typename... T1>
   inline const T & operator()(T1... args) const;
 
   /// Tuple index access operator
   inline T & operator()(std::array<UInt, dim+1> tuple);
   inline const T & operator()(std::array<UInt, dim+1> tuple) const;
 
 private:
   /// Unpacking the arguments of variadic ()
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos,
 			   UInt index, T1... rest) const;
   /// End case for recursion
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos,
 			   UInt index) const;
   /// Computing offset for a tuple index
   inline UInt computeOffset(std::array<UInt, dim+1> tuple) const;
 
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators                                                       */
 /* -------------------------------------------------------------------------- */
 public:
 #define GRID_VEC_OPERATOR(op)					\
   template <typename T1>					\
   Grid & operator op(const Grid<T1, dim> & other)		\
 
   GRID_VEC_OPERATOR(+=);
   GRID_VEC_OPERATOR(*=);
   GRID_VEC_OPERATOR(-=);
   GRID_VEC_OPERATOR(/=);
 #undef GRID_VEC_OPERATOR
 
 #define GRID_SCALAR_OPERATOR(op)			\
   Grid & operator op(const T & e)			\
 
   GRID_SCALAR_OPERATOR(+=);
   GRID_SCALAR_OPERATOR(*=);
   GRID_SCALAR_OPERATOR(-=);
   GRID_SCALAR_OPERATOR(/=);
   GRID_SCALAR_OPERATOR(=);
 #undef GRID_SCALAR_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 /* Move/Copy operators                                                        */
 /* -------------------------------------------------------------------------- */
 
   // = operator
   Grid & operator=(const Grid & other);
   // = operator (not const input, otherwise gcc is confused)
   Grid & operator=(Grid & other);
   // = operator (move)
   Grid & operator=(Grid && other);
 
   // Copy data from another grid
   template <typename T1>
   void copy(const Grid<T1, dim> & other);
   // Move data from another grid
   template <typename T1>
   void move(Grid<T1, dim> && other);
 
 /* -------------------------------------------------------------------------- */
 /* Iterators                                                                  */
 /* -------------------------------------------------------------------------- */
 public:
-  virtual iterator begin(UInt n = 1) {
+  using iterator = typename GridBase<T>::iterator;
+  using const_iterator = typename GridBase<T>::const_iterator;
+
+public:
+  iterator begin(UInt n = 1) override {
     std::vector<UInt> strides(this->strides.begin(), this->strides.end());
     std::vector<UInt> sizes(this->n.begin(), this->n.end());
     sizes.push_back(this->nb_components);
     return iterator(this->getInternalData()+this->offset, 0, n, strides, sizes);
   }
 
-  virtual iterator end(UInt n = 1) {
+  iterator end(UInt n = 1) override {
     std::vector<UInt> strides(this->strides.begin(), this->strides.end());
     std::vector<UInt> sizes(this->n.begin(), this->n.end());
     sizes.push_back(this->nb_components);
     return iterator(this->getInternalData()+this->offset,
 		    this->dataSize(), n, strides, sizes);
   }
 
+  const_iterator begin(UInt n = 1) const override {
+    std::vector<UInt> strides(this->strides.begin(), this->strides.end());
+    std::vector<UInt> sizes(this->n.begin(), this->n.end());
+    sizes.push_back(this->nb_components);
+    return const_iterator(this->getInternalData()+this->offset, 0, n, strides, sizes);
+  }
+
+  const_iterator end(UInt n = 1) const override {
+    std::vector<UInt> strides(this->strides.begin(), this->strides.end());
+    std::vector<UInt> sizes(this->n.begin(), this->n.end());
+    sizes.push_back(this->nb_components);
+    return const_iterator(this->getInternalData()+this->offset,
+			  this->dataSize(), n, strides, sizes);
+  }
+
 /* -------------------------------------------------------------------------- */
 /* Member variables                                                           */
 /* -------------------------------------------------------------------------- */
 protected:
   std::array<UInt, dim> n;
   std::array<UInt, dim+1> strides;
 };
 
 __END_TAMAAS__
 
 
 /* -------------------------------------------------------------------------- */
 /* Inline/template function definitions                                       */
 /* -------------------------------------------------------------------------- */
 
 #include "grid_tmpl.hh"
 
 /* -------------------------------------------------------------------------- */
 
 
 #endif // __GRID_HH__
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 9b77ec6..c2a482c 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,198 +1,151 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_BASE_HH__
 #define __GRID_BASE_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "array.hh"
+#include "iterator.hh"
 #include <cstddef>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /// Namespace for view index classes
 namespace view {
   /// Index class
   struct index {
     index(Int i = -1):i(i){}
     Int i;
   };
 
   /// Blocked index class
   struct blocked : public index {
     blocked(Int i):index(i){}
   };
 
   /// Free index class
   struct free : public index {
     free():index(-1){}
   };
 }
 
 template <typename T>
 class GridBase {
 public:
   /// Constructor by default
   GridBase() = default;
   /// Copy constructor
   GridBase(const GridBase & o): data(o.data),
 				nb_components(o.nb_components),
 				offset(o.offset) {}
   /// Move constructor (transfers data ownership)
   GridBase(GridBase&& o): data(std::move(o.data)),
 			  nb_components(std::exchange(o.nb_components, 1)),
 			  offset(std::exchange(o.offset, 0)) {}
 
   /// Destructor
   virtual ~GridBase() = default;
 
 /* -------------------------------------------------------------------------- */
 /* Iterator class                                                             */
 /* -------------------------------------------------------------------------- */
-protected:
-  struct iterator {
-    /// constructor
-    iterator(T * data,
-	     std::size_t start,
-	     ptrdiff_t offset,
-	     std::vector<UInt> strides,
-	     std::vector<UInt> sizes):
-      data(data),
-      index(start),
-      offset(offset),
-      strides(strides),
-      n(sizes) {}
-    /// destructor
-    ~iterator() {}
-    /// pre-increment
-    inline iterator & operator++() { index += offset; return *this; }
-    /// comparison
-    inline bool operator<(const iterator & a) { return index < a.index; }
-    /// inequality
-    inline bool operator!=(const iterator & a) { return index != a.index; }
-    /// increment with given offset
-    inline iterator & operator+=(ptrdiff_t a) { index += a*offset; return *this;}
-    /// dereference iterator
-    inline T & operator*() {
-      std::vector<UInt> tuple(strides.size());
-      UInt index_copy = index;
-
-      tuple.back() = index_copy % n.back();
-      index_copy -= tuple.back();
-      index_copy /= n.back();
-
-      /// Computing tuple from index
-      for (UInt d = tuple.size()-2 ; d > 0 ; d--) {
-	tuple[d] = index_copy % this->n[d];
-	index_copy -= tuple[d];
-	index_copy /= this->n[d];
-      }
-
-      tuple[0] = index_copy;
-
-      UInt access_offset = 0;
-      for (UInt i = 0 ; i < tuple.size() ; ++i)
-	access_offset += tuple[i] * strides[i];
-      return data[access_offset];
-    }
-    /// needed for OpenMP range calculations
-    inline Int operator-(const iterator & a) const { return (index - a.index)/offset; }
-  protected:
-    T * data;
-    std::size_t index;
-    ptrdiff_t offset;
-    std::vector<UInt> strides;
-    std::vector<UInt> n;
-  };
+public:
+  using iterator = iterator_::iterator<T>;
+  using const_iterator = iterator_::iterator<const T>;
 
 public:
   /// Compute size
   virtual UInt computeSize() const = 0;
   /// Get internal data pointer (const)
   const T * getInternalData() const {return this->data.getPtr();}
   /// Get internal data pointer (non-const)
   T * getInternalData() {return this->data.getPtr();}
   /// Get number of components
   UInt getNbComponents() const {return nb_components;}
   /// Set number of components
   void setNbComponents(UInt n) {nb_components = n;}
   /// Get offset
   UInt getOffset() const {return offset;}
   /// Get total size
   virtual UInt dataSize() const {return this->data.size(); }
   /// Get number of points
   UInt getNbPoints() const {return this->computeSize() / this->getNbComponents(); }
   /// Set components
   void uniformSetComponents(const GridBase<T> & vec);
 
 /* -------------------------------------------------------------------------- */
 /* Iterator methods                                                           */
 /* -------------------------------------------------------------------------- */
   virtual iterator begin(UInt n = 1) = 0;
   virtual iterator end(UInt n = 1) = 0;
+  virtual const_iterator begin(UInt n = 1) const = 0;
+  virtual const_iterator end(UInt n = 1) const = 0;
 
 public:
   template <typename T1>
   void copy(const GridBase<T1> & other) {
     data = other.data;
     nb_components = other.nb_components;
     offset = other.offset;
   }
 
   template <typename T1>
   void move(GridBase<T1> && other) {
     data = std::move(other.data);
     nb_components = std::exchange(other.nb_components, 1);
     offset = std::exchange(other.offset, 0);
   }
 
 protected:
   Array<T> data;
   UInt nb_components = 1;
   UInt offset = 0;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void GridBase<T>::uniformSetComponents(const GridBase<T> &vec) {
   if (!(vec.dataSize() == this->nb_components))
       TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
 
-  const auto size = this->dataSize();
+  auto begin_it = begin(), end_it = end();
+
 #pragma omp parallel for
-  for (UInt i = 0 ; i < size ; i++) {
-    this->data[i] = vec.data[i % this->nb_components];
+  for (auto it = begin() ; it < end_it ; ++it) {
+    UInt i = it - begin_it;
+    *it = vec.data[i % this->nb_components];
   }
 }
 
 
 __END_TAMAAS__
 
 
 #endif // __GRID_BASE_HH__
diff --git a/src/core/grid_tmpl.hh b/src/core/grid_tmpl.hh
index a9a953f..c27c14b 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,195 +1,200 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_TMPL_HH__
 #define __GRID_TMPL_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 
 __BEGIN_TAMAAS__
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::init(const T1 & n, UInt nb_components) {
   std::copy(n.begin(), n.end(), this->n.begin());
   this->nb_components = nb_components;
   this->resize(this->n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeSize() const {
   UInt size = 1;
   for (UInt i = 0 ; i < dim ; i++) size *= n[i];
   size *= this->nb_components;
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos,
 				       UInt index, T1... rest) const {
   offset += index * strides[index_pos];
   return unpackOffset(offset, index_pos+1, rest...); // tail-rec bb
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos,
 				       UInt index) const {
   return offset + index * strides[index_pos];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline T & Grid<T, dim>::operator()(T1... args) {
   /// Checking access integrity
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs == dim + 1 ||
 		nargs == 1 ||
 		nargs == dim,
 		"number of arguments in operator() does not match dimension");
 
   constexpr UInt start = (nargs == 1) ? dim : 0;
   UInt offset = unpackOffset(this->offset, start, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const T & Grid<T, dim>::operator()(T1... args) const {
   /// Checking access integrity
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs == dim + 1 ||
 		nargs == 1 ||
 		nargs == dim,
 		"number of arguments in operator() does not match dimension");
 
   constexpr UInt start = (nargs == 1) ? dim : 0;
   UInt offset = unpackOffset(this->offset, start, args...);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeOffset(std::array<UInt, dim+1> tuple) const {
   UInt offset = this->offset;
   for (UInt i = 0 ; i < dim+1 ; i++) {
     offset += tuple[i] * strides[i];
   }
 
   return offset;
 }
 
 template <typename T, UInt dim>
 inline T & Grid<T, dim>::operator()(std::array<UInt, dim+1> tuple) {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline const T & Grid<T, dim>::operator()(std::array<UInt, dim+1> tuple) const {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define GRID_VEC_OPERATOR_IMPL(op)					\
   template <typename T, UInt dim>					\
   template <typename T1>						\
   inline Grid<T, dim> & Grid<T, dim>::operator op(const Grid<T1, dim> & other) { \
-    TAMAAS_ASSERT(other.computeSize() == this->data.size(), "surface size does not match");\
+    TAMAAS_ASSERT(other.computeSize() == this->computeSize(),		\
+		  "surface size does not match");			\
+    auto begin_it = this->begin(), end_it = this->end();		\
+    auto other_begin = other.begin();					\
     _Pragma("omp parallel for")						\
-    for (UInt i = 0 ; i < this->data.size() ; i++) {			\
-      this->data[i] op other(i);					\
+    for (auto it = begin() ; it < end_it ; ++it) {			\
+      auto other_it = other_begin;					\
+      other_it += it - begin_it;					\
+      *it op *other_it;							\
     }									\
     return *this;							\
   }									\
 
 GRID_VEC_OPERATOR_IMPL(+=);
 GRID_VEC_OPERATOR_IMPL(-=);
 GRID_VEC_OPERATOR_IMPL(*=);
 GRID_VEC_OPERATOR_IMPL(/=);
 #undef GRID_VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim> & Grid<T, dim>::operator=(const Grid<T, dim> & other) {
   this->copy(other);
   return *this;
 }
 
 template <typename T, UInt dim>
 Grid<T, dim> & Grid<T, dim>::operator=(Grid<T, dim> & other) {
   this->copy(other);
   return *this;
 }
 
 template <typename T, UInt dim>
 Grid<T, dim> & Grid<T, dim>::operator=(Grid<T, dim> && other) {
-  this->move(std::move(other));
+  this->move(std::forward<Grid<T, dim>>(other));
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::copy(const Grid<T1, dim> & other) {
   GridBase<T>::copy(other);
   this->n = other.n;
   this->strides = other.strides;
 }
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::move(Grid<T1, dim> && other) {
   GridBase<T>::move(std::forward<Grid<T1, dim>>(other));
   this->n = std::move(other.n);
   this->strides = std::move(other.strides);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Stream operator                                                            */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline std::ostream & operator<<(std::ostream & stream,
 				 const Grid<T, dim> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif // __GRID_TMPL_HH__
diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh
index bc5843a..75aa587 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,215 +1,201 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas 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.
  *
  * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_VIEW_HH__
 #define __GRID_VIEW_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
 class GridView : public Base<T, dim> {
 public:
   typedef typename GridBase<T>::iterator iterator;
 public:
   /// Constructor
   GridView(GridBase<T> & grid_base,
 	   const std::vector<view::index> & multi_index);
   /// Destructor
   virtual ~GridView() = default;
 
 public:
   /// Compute strides and offset
   virtual void computeStrides();
   /// Return view size
   virtual UInt dataSize() const { return this->computeSize(); }
   /// Delete resize
   void resize(const std::array<UInt, dim> &) = delete;
-  /// Print to stream
-  virtual void printself(std::ostream & str) const;
 
 public:
   /// Access operator
   template <typename... T1>
   T & operator()(T1... args);
   /// Const access operator
   template <typename... T1>
   const T & operator()(T1... args) const;
 
 protected:
   Base<T, base_dim> * grid;
   std::vector<view::index> indexes;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implentation                                                      */
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
 GridView<Base, T, base_dim, dim>::GridView(GridBase<T> & grid_base,
   const std::vector<view::index> & multi_index):
   Base<T, dim>(),
   grid(nullptr),
   indexes(multi_index)
 {
   static_assert(base_dim >= dim,
 		"view dimension must be >= than the base class");
 
   // Checking view type
   grid = dynamic_cast<Base<T, base_dim> *>(&grid_base);
 
   if (!grid) TAMAAS_EXCEPTION("given base type is incompatible with view");
   
   // Checking view integrity
   if (multi_index.size() == base_dim) {
     this->nb_components = grid->getNbComponents();
   } else if (multi_index.size() == base_dim+1) {
     this->nb_components = 1;
   } else {
     TAMAAS_EXCEPTION("Multi-index has wrong number of components");
   }
 
   UInt view_dim = 0;
   for (auto & it : multi_index) {
     if (it.i == -1) view_dim += 1;
   }
 
   if (view_dim != dim)
     TAMAAS_EXCEPTION("Multi-index has the wrong number of free components");
 
   if (indexes.size() == base_dim
       && this->nb_components == 1) indexes.push_back(view::blocked(0));
 
   this->data.wrapMemory(grid->getInternalData(), grid->dataSize());
   computeStrides();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
 void GridView<Base, T, base_dim, dim>::computeStrides() {
   /// Copying sizes with number of components in last position
   auto stride_it = this->strides.begin();
   auto size_it = this->n.begin();
   auto grid_stride_it = grid->getStrides().begin();
   auto grid_size_it = grid->sizes().begin();
 
   for (auto & index : indexes) {
     if (index.i == -1) {
       *stride_it = *grid_stride_it;
       *size_it = *grid_size_it;
       ++stride_it, ++size_it;
     } else
       this->offset += index.i * *grid_stride_it;
     ++grid_stride_it, ++grid_size_it;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
 template <typename... T1>
 T & GridView<Base, T, base_dim, dim>::operator()(T1... args) {
   constexpr UInt nargs = sizeof...(args);
 
   static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
 		"number of arguments in operator() does not match dimension");
 
   if (nargs == 1) {
     UInt index = UInt(args...);
     std::array<UInt, dim+1> tuple = {0};
 
     tuple[dim] = index % this->nb_components;
     index -= tuple[dim];
     index /= this->nb_components;
 
     /// Computing tuple from index
     for (UInt d = dim-1 ; d > 0 ; d--) {
       tuple[d] = index % this->n[d];
       index -= tuple[d];
       index /= this->n[d];
     }
 
     tuple[0] = index;
     return Base<T, dim>::operator()(tuple);
   }
   else {
     return Base<T, dim>::operator()(args...);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
 template <typename... T1>
 const T & GridView<Base, T, base_dim, dim>::operator()(T1... args) const {
   constexpr UInt nargs = sizeof...(args);
 
   static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
 		"number of arguments in operator() does not match dimension");
 
   if (nargs == 1) {
     UInt index = UInt(args...);
     std::array<UInt, dim+1> tuple = {0};
 
     tuple[dim] = index % this->nb_components;
     index -= tuple[dim];
     index /= this->nb_components;
 
     /// Computing tuple from index
     for (UInt d = dim-1 ; d > 0 ; d--) {
       tuple[d] = index % this->n[d];
       index -= tuple[d];
       index /= this->n[d];
     }
 
     tuple[0] = index;
     return Base<T, dim>::operator()(tuple);
   }
   else {
     return Base<T, dim>::operator()(args...);
   }
 }
 
-/* -------------------------------------------------------------------------- */
-
-template <template <typename, UInt> class Base, typename T, UInt base_dim, UInt dim>
-void GridView<Base, T, base_dim, dim>::printself(std::ostream & str) const {
-  str << "GridView(" << dim << ", " << this->nb_components << ") {";
-  const auto size = this->dataSize();
-  for (UInt i = 0 ; i < size ; i++) {
-    str << (*this)(i) << ", ";
-  }
-  str << "\b\b}";
-}
-
 
 __END_TAMAAS__
 
 #endif // __GRID_VIEW_HH__
diff --git a/src/core/iterator.hh b/src/core/iterator.hh
new file mode 100644
index 0000000..6e13911
--- /dev/null
+++ b/src/core/iterator.hh
@@ -0,0 +1,141 @@
+/**
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
+ * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+ * Solides)
+ *
+ * Tamaas 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.
+ *
+ * Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef __ITERATOR_HH__
+#define __ITERATOR_HH__
+/* -------------------------------------------------------------------------- */
+#include "tamaas.hh"
+#include <vector>
+#include <utility>
+#include <cstddef>
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+
+namespace iterator_ {
+
+template <typename T>
+class iterator {
+public:
+  typedef T value_type;
+  typedef std::ptrdiff_t difference_type;
+  typedef T * pointer;
+  typedef T & reference;
+  
+/* -------------------------------------------------------------------------- */
+/* Constructors/Destructors                                                   */
+/* -------------------------------------------------------------------------- */
+  /// constructor
+  iterator(pointer data,
+	   std::size_t start,
+	   ptrdiff_t offset,
+	   std::vector<UInt> strides,
+	   std::vector<UInt> sizes):
+    data(data),
+    index(start),
+    offset(offset),
+    strides(strides),
+    n(sizes) {}
+  /// copy constructor
+  iterator(const iterator & o):
+    data(o.data),
+    index(o.index),
+    offset(o.offset),
+    strides(o.strides),
+    n(o.n) {}
+  /// destructor
+  ~iterator() {}
+  /// copy operator
+  iterator & operator=(const iterator & o) {
+    data = o.data;
+    index = o.index;
+    offset = o.offset;
+    strides = o.strides;
+    n = o.n;
+  }
+
+  /// dereference iterator
+  inline reference operator*() {
+    std::vector<UInt> tuple(strides.size());
+    UInt index_copy = index;
+
+    tuple.back() = index_copy % n.back();
+    index_copy -= tuple.back();
+    index_copy /= n.back();
+
+    /// Computing tuple from index
+    for (UInt d = tuple.size()-2 ; d > 0 ; d--) {
+      tuple[d] = index_copy % this->n[d];
+      index_copy -= tuple[d];
+      index_copy /= this->n[d];
+    }
+
+    tuple[0] = index_copy;
+
+    UInt access_offset = 0;
+    for (UInt i = 0 ; i < tuple.size() ; ++i)
+      access_offset += tuple[i] * strides[i];
+    return data[access_offset];
+  }
+
+  /// increment with given offset
+  inline iterator & operator+=(difference_type a) {
+    index += a*offset;
+    return *this;
+  }
+
+  /// increment iterator
+  iterator & operator++() {
+    *this += 1;
+    return *this;
+  }
+
+  /// comparison
+  inline bool operator<(const iterator & a) {
+    return index < a.index;
+  }
+
+  /// inequality
+  inline bool operator!=(const iterator & a) {
+    return index != a.index;
+  }
+
+  /// needed for OpenMP range calculations
+  inline Int operator-(const iterator & a) const {
+    return (index - a.index)/offset;
+  }
+
+protected:
+  T * data;
+  std::size_t index;
+  ptrdiff_t offset;
+  std::vector<UInt> strides;
+  std::vector<UInt> n;
+};
+
+} // end namespace iterator
+
+__END_TAMAAS__
+/* -------------------------------------------------------------------------- */
+#endif // __ITERATOR_HH__