diff --git a/src/grid.hh b/src/grid.hh
index 422403d..4f4a23c 100644
--- a/src/grid.hh
+++ b/src/grid.hh
@@ -1,298 +1,298 @@
 /**
  *
  * @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 "array.hh"
 #include <cstddef>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Multi-dimensional & multi-component array class
  */
 template <typename T, UInt dim>
 class Grid {
 /* -------------------------------------------------------------------------- */
 /* Constructors                                                               */
 /* -------------------------------------------------------------------------- */
 public:
   /// Constructor by default (empty array)
   Grid();
   /// Constructor
   Grid(const UInt n[dim], const Real L[dim], UInt nb_components);
   /// Destructor
   virtual ~Grid();
 
 /* -------------------------------------------------------------------------- */
 /* Iterator class                                                             */
 /* -------------------------------------------------------------------------- */
   struct iterator {
     /// constructor
     iterator(T * start, ptrdiff_t offset):ptr(start), offset(offset) {}
     /// destructor
     ~iterator() {}
     /// pre-increment
     iterator & operator++() { ptr += offset; return *this; }
     /// comparison
     bool operator<(const iterator & a) { return ptr < a.ptr; }
     /// increment with given offset
-    iterator & operator+=(ptrdiff_t a) { ptr += a; return *this;}
+    iterator & operator+=(ptrdiff_t a) { ptr += a*offset; return *this;}
     /// dereference iterator
     T * operator*() { return ptr; }
     /// needed for OpenMP range calculations
     Int operator-(const iterator & a) const { return (ptr - a.ptr)/offset; }
   private:
     T * ptr;
     ptrdiff_t offset;
   };
 
 public:
   iterator begin(UInt n = 1) { return iterator(data.getPtr(), n); }
   iterator end(UInt n = 1) { return iterator(data.getPtr()+computeSize(), n); }
 
 public:
 /* -------------------------------------------------------------------------- */
 /* Common operations                                                          */
 /* -------------------------------------------------------------------------- */
   /// Resize array
   void resize(const UInt n[dim]);
   /// Compute size
   inline UInt computeSize() const;
   /// Print
   void printself(std::ostream & str) const;
   /// Get sizes
   const UInt * sizes() const {return n;}
   /// Get total size
   UInt dataSize() const {return data.size(); }
   /// Get number of points
   UInt getNbPoints() const {return dataSize() / getNbComponents(); }
   /// Get lengths
   const Real * lengths() const {return L;}
   /// Get number of components
   UInt getNbComponents() const {return nb_components;}
   /// Set number of components
   void setNbComponents(UInt n) {nb_components = n;}
   /// Get internal data pointer (const)
   const T * getInternalData() const {return data.getPtr();}
   /// Get internal data pointer (non-const)
   T * getInternalData() {return data.getPtr();}
   /// Set components
   void uniformSetComponents(const Grid<T, 1> & vec);
 
 /* -------------------------------------------------------------------------- */
 /* Access operators                                                           */
 /* -------------------------------------------------------------------------- */
   template <typename... T1>
   inline T & operator()(T1... args);
   template <typename... T1>
   inline const T & operator()(T1... args) const;
 
   inline T & operator()(UInt tuple[dim+1]);
   inline const T & operator()(UInt tuple[dim+1]) const;
 
 private:
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index, T1... rest) const;
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index) const;
   inline UInt computeOffset(UInt tuple[dim+1]) const;
 
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators                                                       */
 /* -------------------------------------------------------------------------- */
 public:
 #define GRID_VEC_OPERATOR(op)					\
   template <typename T1>					\
   void 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)		\
   void operator op(const T & e)			\
 
   GRID_SCALAR_OPERATOR(+=);
   GRID_SCALAR_OPERATOR(*=);
   GRID_SCALAR_OPERATOR(-=);
   GRID_SCALAR_OPERATOR(/=);
   GRID_SCALAR_OPERATOR(=);
 #undef GRID_SCALAR_OPERATOR
 
   // = operator
   void operator=(const Grid<T, dim> & other);
   // = operator (not const input, otherwise gcc is confused)
   void operator=(Grid<T, dim> & other);
 
   // Copy data from another grid
   template <typename T1>
   void copy(const Grid<T1, dim> & other);
 
 /* -------------------------------------------------------------------------- */
 /* Member variables                                                           */
 /* -------------------------------------------------------------------------- */
 protected:
   UInt n[dim];
   Real L[dim];
   Array<T> data;
   UInt nb_components;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline function definitions                                                */
 /* -------------------------------------------------------------------------- */
 
 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 *= nb_components;
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index, T1... rest) const {
   UInt size = sizeof...(T1);
   offset += index;
   offset *= (size == 1) ? this->nb_components : this->n[dim-size+1];
   return unpackOffset(offset, rest...); // tail-rec bb
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index) const {
   return offset + index;
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline T & Grid<T, dim>::operator()(T1... args) {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   UInt offset = unpackOffset(0, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const T & Grid<T, dim>::operator()(T1... args) const {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   UInt offset = unpackOffset(0, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeOffset(UInt tuple[dim+1]) const {
   UInt offset = 0;
   for (UInt i = 0 ; i < dim-1 ; i++) {
     offset += tuple[i];
     offset *= n[i+1];
   }
 
   offset += tuple[dim-1];
   offset *= nb_components;
   return offset + tuple[dim];
 }
 
 template <typename T, UInt dim>
 inline T & Grid<T, dim>::operator()(UInt tuple[dim+1]) {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline const T & Grid<T, dim>::operator()(UInt tuple[dim+1]) const {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define GRID_VEC_OPERATOR_IMPL(op)					\
   template <typename T, UInt dim>					\
   template <typename T1>						\
   inline void Grid<T, dim>::operator op(const Grid<T1, dim> & other) {	\
     TAMAAS_ASSERT(other.computeSize() == data.size(), "surface size does not match");\
     _Pragma("omp parallel for")						\
     for (UInt i = 0 ; i < this->data.size() ; i++) {			\
       data[i] op other(i);						\
     }									\
   }									\
 
 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>
 void Grid<T, dim>::operator=(const Grid<T, dim> & other) {
   this->copy(other);
 }
 
 template <typename T, UInt dim>
 void Grid<T, dim>::operator=(Grid<T, dim> & other) {
   this->copy(other);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::copy(const Grid<T1, dim> & other) {
   resize(other.sizes());
 #pragma omp parallel for
   for (UInt i = 0 ; i < data.size() ; i++)
     data[i] = other(i);
 }
 
 /* -------------------------------------------------------------------------- */
 /* 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_HH__