diff --git a/src/core/fftransform.hh b/src/core/fftransform.hh
index 586d0f7..e9d6148 100644
--- a/src/core/fftransform.hh
+++ b/src/core/fftransform.hh
@@ -1,115 +1,119 @@
 /**
  *
  * @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 FFTRANSFORM_HH
 #define FFTRANSFORM_HH
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "types.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 class FFTransform {
 public:
   /// Constructor
   FFTransform(Grid<T, dim> & real,
               GridHermitian<T, dim> & spectral);
 
   /// Destructor
   virtual ~FFTransform();
 
 public:
   /// Perform FFT
   virtual void forwardTransform() = 0;
 
   /// Perform IFFT
   virtual void backwardTransform() = 0;
 
   /// Normalize the real surface after backward transform
   virtual void normalize();
 
 public:
   /// Fill a grid with modevectors: boolean if grid has hermitian dimensions
   template <bool hermitian>
-  static void computeFrequencies(Grid<T, dim> & freq);
+  static Grid<T, dim> computeFrequencies(const std::array<UInt, dim> & sizes);
 
 protected:
   Grid<T, dim> & real;
   GridHermitian<T, dim> & spectral;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template<typename T, UInt dim>
 template<bool hermitian>
-void FFTransform<T, dim>::computeFrequencies(Grid<T, dim> & freq) {
+Grid<T, dim>
+FFTransform<T, dim>::computeFrequencies(const std::array<UInt, dim> & sizes) {
   // If hermitian is true, we suppose the dimensions of freq are
   // reduced based on hermitian symetry and that it has dim components
 
-  auto n = freq.sizes();
+  auto & n = sizes;
+  Grid<T, dim> freq(n, dim);
 
 #pragma omp parallel for
   for (auto it = freq.begin(dim) ; it < freq.end(dim) ; ++it) {
     VectorProxy<T> wavevector(&(*it), dim);
     UInt index = freq.getNbPoints() - (freq.end(dim) - it);
-    std::array<UInt, dim> tuple = {0};
+    std::array<UInt, dim> tuple = {{0}};
 
     /// Computing tuple from index
     for (UInt d = dim-1 ; d > 0 ; d--) {
       tuple[d] = index % n[d];
       index -= tuple[d];
       index /= n[d];
     }
 
     tuple[0] = index;
 
     UInt dmax = dim;
 
     if (hermitian) {
       dmax = dim - 1;
       wavevector(dim-1) = tuple[dim-1];
     }
 
     for (UInt d = 0 ; d < dmax ; d++) {
       T td = tuple[d];
       T nd = n[d];
       T q = (tuple[d] <= n[d]/2) ? td : td-nd;
       wavevector(d) = q;
     }
 
   }
+
+  return freq;
 }
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 #endif // FFTRANSFORM_HH
diff --git a/src/core/grid.hh b/src/core/grid.hh
index 869ab05..fc26beb 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,218 +1,218 @@
 /**
  *
  * @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;
   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;
+  inline UInt computeSize() const override;
   /// 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:
   using iterator = typename GridBase<T>::iterator;
   using const_iterator = typename GridBase<T>::const_iterator;
 
 public:
   iterator begin(UInt n = 1) override {
     return iterator(this->getInternalData()+this->offset, 0, n);
   }
 
   iterator end(UInt n = 1) override {
     return iterator(this->getInternalData()+this->offset,
 		    this->dataSize(), n);
   }
 
   const_iterator begin(UInt n = 1) const override {
     return const_iterator(this->getInternalData()+this->offset, 0, n);
   }
 
   const_iterator end(UInt n = 1) const override {
     return const_iterator(this->getInternalData()+this->offset,
 			  this->dataSize(), n);
   }
 
 /* -------------------------------------------------------------------------- */
 /* 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_hermitian.hh b/src/core/grid_hermitian.hh
index 0d6752e..e78d0c2 100644
--- a/src/core/grid_hermitian.hh
+++ b/src/core/grid_hermitian.hh
@@ -1,180 +1,180 @@
 /**
  *
  * @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_HERMITIAN_HH__
 #define __GRID_HERMITIAN_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <vector>
 #include <complex>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 static std::complex<Real> dummy;
 
 /**
  * @brief Multi-dimensional, multi-component herimitian array
  *
  * This class represents an array of hermitian data, meaning it has dimensions
  * of: n1 * n2 * n3 * ... * (nx / 2 + 1)
  * 
  * However, it acts as a fully dimensioned array, returning a dummy reference
  * for data outside its real dimension, allowing one to write full (and
  * inefficient) loops without really worrying about the reduced dimension.
  *
  * It would however be better to just use the true dimensions of the surface
  * for all intents and purposes, as it saves computation time.
  */
 template <typename T, UInt dim>
 class GridHermitian : public Grid<std::complex<T>, dim> {
 public:
   typedef std::complex<T> value_type;
 /* -------------------------------------------------------------------------- */
 /* Constructors                                                               */
 /* -------------------------------------------------------------------------- */
 public:
   GridHermitian() = default;
   GridHermitian(const GridHermitian & o) = default;
   GridHermitian(GridHermitian && o) = default;
   using Grid<value_type, dim>::Grid;
   using Grid<std::complex<T>, dim>::operator=;
 /* -------------------------------------------------------------------------- */
 /* Access operators                                                           */
 /* -------------------------------------------------------------------------- */
   template <typename... T1>
   inline std::complex<T> & operator()(T1... args);
   template <typename... T1>
   inline const std::complex<T> & operator()(T1... args) const;
 
 
   static std::array<UInt, dim> hermitianDimensions(const std::array<UInt, dim> & n) {
     std::array<UInt, dim> hn(n);
     hn[dim-1] /= 2;
     hn[dim-1] += 1;
     return hn;
   }
 
 private:
   template <typename... T1>
   inline void packTuple(UInt * t, UInt i, T1... args) const;
   template <typename... T1>
   inline void packTuple(UInt * t, UInt i) const;
 
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline function definitions                                                */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt * t, UInt i, T1... args) const {
   *t = i;
   packTuple(t+1, args...);
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt * t, UInt i) const {
   *t = i;
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline std::complex<T> & GridHermitian<T, dim>::operator()(T1... args) {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   if (sizeof...(T1) == 1 && dim != 1) {
     UInt tuple[1] = {0};
     packTuple(tuple, args...);
     if (tuple[0] >= this->data.size()) {
       TAMAAS_DEBUG_EXCEPTION("out of bonds access on hermitian surface");
       dummy = std::complex<T>(0,0);
       return dummy;
     }
     else return this->Grid<std::complex<T>, dim>::operator()(args...);
   }
 
-  std::array<UInt, dim+1> tuple = {0};
+  std::array<UInt, dim+1> tuple = {{0}};
   packTuple(tuple.data(), args...);
 
   // reverse tuple if necessary
   if (tuple[dim-1] >= this->n[dim-1]) {
     for (UInt i = 0 ; i < dim ; i++) {
       if (tuple[i] && i != dim-1)
 	tuple[i] = this->n[i]-tuple[i];
       else if (tuple[i] && i == dim-1)
 	tuple[i] = 2*(this->n[i]-1) - tuple[i];
     }
 
     dummy = std::conj(this->Grid<std::complex<T>, dim>::operator()(tuple));
     return dummy;
   } else {
     return this->Grid<std::complex<T>, dim>::operator()(tuple);
   }
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const std::complex<T> & GridHermitian<T, dim>::operator()(T1... args) const {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   if (sizeof...(T1) == 1 && dim != 1) {
     UInt tuple[1] = {0};
     packTuple(tuple, args...);
     if (tuple[0] >= this->data.size()) {
       TAMAAS_DEBUG_EXCEPTION("out of bonds access on hermitian surface");
       dummy = std::complex<T>(0,0);
       return dummy;
     }
     else return this->Grid<std::complex<T>, dim>::operator()(args...);
   }
 
   std::array<UInt, dim+1> tuple = {0};
   packTuple(tuple.data(), args...);
 
   // reverse tuple if necessary
   if (tuple[dim-1] >= this->n[dim-1]) {
     for (UInt i = 0 ; i < dim ; i++) {
       if (tuple[i] && i != dim-1)
 	tuple[i] = this->n[i]-tuple[i];
       else if (tuple[i] && i == dim-1)
 	tuple[i] = 2*(this->n[i]-1) - tuple[i];
     }
 
     dummy = std::conj(this->Grid<std::complex<T>, dim>::operator()(tuple));
     return dummy;
   } else {
     return this->Grid<std::complex<T>, dim>::operator()(tuple);
   }
 }
 
 __END_TAMAAS__
 
 
 
 #endif // __GRID_HERMITIAN_HH__
diff --git a/src/core/grid_tmpl.hh b/src/core/grid_tmpl.hh
index c27c14b..3e09b20 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,200 +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->computeSize(),		\
 		  "surface size does not match");			\
     auto begin_it = this->begin(), end_it = this->end();		\
-    auto other_begin = other.begin();					\
-    _Pragma("omp parallel for")						\
+    auto other_begin = other.begin(), other_it = other.begin();		\
+    _Pragma("omp parallel for firstprivate(other_it)")			\
     for (auto it = begin() ; it < end_it ; ++it) {			\
-      auto other_it = other_begin;					\
+      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::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 702a6e6..7d80aaa 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,240 +1,265 @@
 /**
  *
  * @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__
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief View type on grid
  * /!\ iterators do not work on this type of grid
+ *
  */
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 class GridView : public Base<T, dim> {
 public:
   /// Constructor
   GridView(GridBase<T> & grid_base,
 	   const std::vector<view::index> & multi_index);
   /// Move constructor
   GridView(GridView&& o) :
     Base<T, dim>(std::forward<Base<T, dim>>(o)),
     grid(std::exchange(o.grid, nullptr)),
     indexes(std::move(o.indexes)) {}
 
   /// Destructor
   virtual ~GridView() = default;
 
 public:
   /// Compute strides and offset
   void computeStrides() override;
   /// Return view size
   UInt dataSize() const override { return this->computeSize(); }
   /// Delete resize
   void resize(const std::array<UInt, dim> &) = delete;
   /// Get a view only for the last component
   GridView lastComponentView() const;
 
+public:
+  using iterator = typename GridBase<T>::iterator;
+  using const_iterator = typename GridBase<T>::const_iterator;
+
+  iterator begin(UInt n [[gnu::unused]] = 1) override {
+    TAMAAS_EXCEPTION("iterators on views do not work as expected");
+    return iterator(nullptr, 0, 0);
+  }
+
+  iterator end(UInt n [[gnu::unused]] = 1) override {
+    TAMAAS_EXCEPTION("iterators on views do not work as expected");
+    return iterator(nullptr, 0, 0);
+  }
+
+  const_iterator begin(UInt n [[gnu::unused]] = 1) const override {
+    TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
+    return const_iterator(nullptr, 0, 0);
+  }
+
+  const_iterator end(UInt n [[gnu::unused]] = 1) const override {
+    TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
+    return const_iterator(nullptr, 0, 0);
+  }
+
 public:
   /// Access operator (1D access is extremely slow)
   template <typename... T1>
   T & operator()(T1... args);
   /// Const access operator (1D access is extremely slow)
   template <typename... T1>
   const T & operator()(T1... args) const;
 
 protected:
   Base<T, base_dim> * grid;
   std::vector<view::index> indexes;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation                                                    */
 /* -------------------------------------------------------------------------- */
 
 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>
 GridView<Base, T, base_dim, dim>
 GridView<Base, T, base_dim, dim>::lastComponentView() const {
   GridView<Base, T, base_dim, dim> tmp(*grid, indexes);
   tmp.indexes.back() = view::blocked(this->nb_components-1);
   tmp.computeStrides();
   return tmp;
 }
 
 
 __END_TAMAAS__
 
 #endif // __GRID_VIEW_HH__
diff --git a/src/core/iterator.hh b/src/core/iterator.hh
index 55e1caa..7ba034d 100644
--- a/src/core/iterator.hh
+++ b/src/core/iterator.hh
@@ -1,184 +1,185 @@
 /**
  *
  * @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):
     data(data),
     index(start),
     step(offset),
     access_offset(0) {}
   /// copy constructor
   iterator(const iterator & o):
     data(o.data),
     index(o.index),
     step(o.step),
     access_offset(o.access_offset) {}
   /// destructor
   ~iterator() {}
   /// copy operator
   iterator & operator=(const iterator & o) {
     data = o.data;
     index = o.index;
     step = o.step;
     access_offset = o.access_offset;
+    return *this;
   }
 
 protected:
   /// Compute access offset for data from index
   inline void computeAccessOffset() {
     access_offset = index;
   }
 
 public:
   /// dereference iterator
   inline reference operator*() {
     return data[access_offset];
   }
 
   /// increment with given offset
   inline iterator & operator+=(difference_type a) {
     index += a*step;
     computeAccessOffset();
     return *this;
   }
 
   /// increment iterator
   inline 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 difference_type operator-(const iterator & a) const {
     return (index - a.index)/step;
   }
 
 protected:
   T * data;
   std::size_t index;
   difference_type step;
   difference_type access_offset;
 };
 
 
 template <typename T>
 class iterator_view : iterator<T> {
 public:
   typedef T value_type;
   typedef std::ptrdiff_t difference_type;
   typedef T * pointer;
   typedef T & reference;
 
   /// Constructor
   iterator_view(pointer data,
 		std::size_t start,
 		ptrdiff_t offset,
 		std::vector<UInt> strides,
 		std::vector<UInt> sizes):
     iterator<T>(data, start, offset),
     strides(strides),
     n(sizes),
     tuple(strides.size()) {}
 
   iterator_view(const iterator_view & o):
     iterator<T>(o),
     strides(o.strides),
     n(o.n),
     tuple(o.tuple) {}
 
   iterator_view & operator=(const iterator_view & o) {
     iterator<T>::operator=(o);
     strides = o.strides;
     n = o.n;
     tuple = o.tuple;
   }
 
 protected:
   void computeAccessOffset() {
     UInt index_copy = this->index;
 
     this->tuple.back() = index_copy % this->n.back();
     index_copy -= this->tuple.back();
     index_copy /= this->n.back();
 
     /// Computing tuple from index
     for (UInt d = this->tuple.size()-2 ; d > 0 ; d--) {
       this->tuple[d] = index_copy % this->n[d];
       index_copy -= this->tuple[d];
       index_copy /= this->n[d];
     }
 
     this->tuple[0] = index_copy;
 
     this->access_offset = 0;
     for (UInt i = 0 ; i < this->tuple.size() ; ++i)
       this->access_offset += this->tuple[i] * this->strides[i];
   }
 
   std::vector<UInt> strides;
   std::vector<UInt> n;
   std::vector<UInt> tuple;
 };
 
 } // end namespace iterator_
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif // __ITERATOR_HH__
diff --git a/src/model/bem.hh b/src/model/bem.hh
index 922be8c..7f90f4c 100644
--- a/src/model/bem.hh
+++ b/src/model/bem.hh
@@ -1,61 +1,60 @@
 /**
  *
  * @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 __BEM_HH__
 #define __BEM_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "model.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 class Bem {
 public:
   /// Constructor
-  Bem(Model * model, model_type type): model(model), type(type) {}
+  Bem(Model * model): model(model) {}
 
 public:
   /// Destructor
   virtual ~Bem() {}
 
 public:
   /// Initialize influence coefficients
   virtual void initInfluence() = 0;
   /// Solve Neumann problem
   virtual void solveNeumann(const GridBase<Real> * neumann, GridBase<Real> * dirichlet) const = 0;
   /// Solve Dirichlet problem
   virtual void solveDirichlet(const GridBase<Real> * dirichlet, GridBase<Real> * neumann) const = 0;
   
 protected:
   Model * model;
-  model_type type;
 };
 
 
 __END_TAMAAS__
 
 #endif // __BEM_HH__
diff --git a/src/model/model.cpp b/src/model/model.cpp
index 2bf5ee1..a199903 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,68 +1,76 @@
 /**
  *
  * @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 "model.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
-Model::Model(const std::vector<Real> & system_size):
+Model::Model(const std::vector<Real> & system_size,
+	     const std::vector<UInt> & discretization):
   system_size(system_size),
+  discretization(discretization),
   pressure(nullptr),
   displacement(nullptr)
 {}
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setElasticity(Real E, Real nu) {
   this->E = E; this->nu = nu;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Model::getHertzModulus() const {
   return E / (1-nu*nu);
 }
 
 /* -------------------------------------------------------------------------- */
 
 GridBase<Real> & Model::getPressure() {
   return *pressure;
 }
 
 /* -------------------------------------------------------------------------- */
 
 GridBase<Real> & Model::getDisplacement() {
   return *pressure;
 }
 
 /* -------------------------------------------------------------------------- */
 
 const std::vector<Real> & Model::getSystemSize() const {
   return system_size;
 }
 
+/* -------------------------------------------------------------------------- */
+
+const std::vector<UInt> & Model::getDiscretization() const {
+  return discretization;
+}
+
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
diff --git a/src/model/model.hh b/src/model/model.hh
index 8637dd1..be08830 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,67 +1,70 @@
 /**
  *
  * @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 __MODEL_HH__
 #define __MODEL_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid_base.hh"
 #include <memory>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
-/// Boundary model: information here is only representing the boundary
+/// Model containing pressure and displacement
 class Model {
 protected:
   /// Constructor
-  Model(const std::vector<Real> & system_size);
+  Model(const std::vector<Real> & system_size,
+	const std::vector<UInt> & discretization);
 public:
   /// Destructor
   virtual ~Model() = default;
 
 public:
   /// Set elasticity parameters
   void setElasticity(Real E, Real nu);
   /// Get Hertz contact modulus
   Real getHertzModulus() const;
   
 public:
   GridBase<Real> & getPressure();
   GridBase<Real> & getDisplacement();
   virtual GridBase<Real> & getBoundaryDisplacement() = 0;
   const std::vector<Real> & getSystemSize() const;
+  const std::vector<UInt> & getDiscretization() const;
 
 protected:
   Real E = 1, nu = 0;
   std::vector<Real> system_size;
+  std::vector<UInt> discretization;
   std::unique_ptr<GridBase<Real>> pressure;
   std::unique_ptr<GridBase<Real>> displacement;
 };
 
 __END_TAMAAS__
 
 #endif // __MODEL_HH__
diff --git a/src/model/model_template.cpp b/src/model/model_template.cpp
index 2b65a84..347775b 100644
--- a/src/model/model_template.cpp
+++ b/src/model/model_template.cpp
@@ -1,69 +1,69 @@
 /**
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2017 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 "model_template.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 ModelTemplate<type>::ModelTemplate(const std::vector<Real> & system_size,
 				   const std::vector<UInt> & discretization):
-  Model(system_size),
+  Model(system_size, discretization),
   boundary_displacement(nullptr)
 {
   constexpr UInt dim = trait::dimension;
   constexpr UInt nb_components = trait::components;
 
   if (system_size.size() != dim)
     TAMAAS_EXCEPTION("System size does not match model type");
   if (discretization.size() != dim)
     TAMAAS_EXCEPTION("Discretization size does not match model type");
 
   std::array<UInt, dim-1> pressure_size;
   std::copy(discretization.begin()+1, discretization.end(),
 	    pressure_size.begin());
   pressure = new Grid<Real, dim-1>(pressure_size, nb_components);
   displacement = new Grid<Real, dim>(discretization, nb_components);
   boundary_displacement = new ViewType(*displacement, trait::indices);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 GridBase<Real> & ModelTemplate<type>::getBoundaryDisplacement() {
   return *boundary_displacement;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template<> class ModelTemplate<model_type::basic_1d>;
 template<> class ModelTemplate<model_type::basic_2d>;
 template<> class ModelTemplate<model_type::surface_1d>;
 template<> class ModelTemplate<model_type::surface_2d>;
 template<> class ModelTemplate<model_type::volume_1d>;
 template<> class ModelTemplate<model_type::volume_2d>;
 
 __END_TAMAAS__
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index f71b7c6..217529e 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,81 +1,107 @@
 /**
  *
  * @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 "westergaard.hh"
+#include "fftransform.hh"
+#include "types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
-template <UInt dim>
-Westergaard<dim>::Westergaard(Model * model,
-			      model_type type,
-			      const std::array<UInt, dim> & sizes):
-  Bem(model, type)
+template <model_type type>
+Westergaard<type>::Westergaard(Model * model):
+  Bem(model)
 {
-  auto hermitian_sizes = GridHermitian<Real, dim>::hermitianDimensions(sizes);
-  UInt nb_components = 0;
+  auto sizes = model->getDiscretization();
+  auto hermitian_sizes =
+    GridHermitian<Real, trait::dimension>::hermitianDimensions(sizes);
+  UInt nb_components = trait::components;
 
   // Cases where we want stuff below the surface in normal contact
   nb_components *= nb_components; // square for matrix
 
   influence.setNbComponents(nb_components);
   influence.resize(sizes);
 }
 
 /* -------------------------------------------------------------------------- */
 
-template <UInt dim>
-Westergaard<dim>::~Westergaard()
-{}
-
-/* -------------------------------------------------------------------------- */
+template <model_type type>
+void Westergaard<type>::initInfluence() {
+  TAMAAS_EXCEPTION("Not implemented");
+}
 
-template <UInt dim>
-void Westergaard<dim>::initInfluence() {
+template <>
+void Westergaard<model_type::basic_2d>::initInfluence() {
+  auto wavevectors =
+    FFTransform<Real, 2>::computeFrequencies<true>(influence.sizes());
+  auto wave_begin = wavevectors.begin(), wave_it = wavevectors.begin();
 
+  VectorProxy<Real> domain(const_cast<Real*>(model->getSystemSize().data()), 2);
+  Real E_star = model->getHertzModulus();
+#pragma omp parallel for firstprivate(wave_it)
+  for (auto inf_it = influence.begin(1) ;
+       inf_it < influence.end(1) ;
+       ++inf_it) {
+    wave_it = wave_begin;
+    wave_it += inf_it - influence.end(1);
+    VectorProxy<Real> q_vec(&(*wave_it), 2);
+    q_vec *= 2 * M_PI;
+    q_vec /= domain;
+    Real q = q_vec.l2norm();
+    *inf_it = 2. / (q*E_star);
+  }
+       
 }
 
 /* -------------------------------------------------------------------------- */
 
-template <UInt dim>
-void Westergaard<dim>::solveNeumann(const GridBase<Real> * neumann,
-				    GridBase<Real> * dirichlet) const {
+template <model_type type>
+void Westergaard<type>::solveNeumann(const GridBase<Real> * neumann,
+				     GridBase<Real> * dirichlet) const {
   
 }
 
 /* -------------------------------------------------------------------------- */
 
-template <UInt dim>
-void Westergaard<dim>::solveDirichlet(const GridBase<Real> * dirichlet,
-				      GridBase<Real> * neumann) const {
+template <model_type type>
+void Westergaard<type>::solveDirichlet(const GridBase<Real> * dirichlet,
+				       GridBase<Real> * neumann) const {
 
 }
 
 /* -------------------------------------------------------------------------- */
-  
+/* Template instanciation                                                     */
+/* -------------------------------------------------------------------------- */
+template<> class Westergaard<model_type::basic_1d>;
+//template<> class Westergaard<model_type::basic_2d>;
+template<> class Westergaard<model_type::surface_1d>;
+template<> class Westergaard<model_type::surface_2d>;
+template<> class Westergaard<model_type::volume_1d>;
+template<> class Westergaard<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/src/model/westergaard.hh b/src/model/westergaard.hh
index 27b49fc..fc49425 100644
--- a/src/model/westergaard.hh
+++ b/src/model/westergaard.hh
@@ -1,60 +1,61 @@
 /**
  *
  * @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 __WESTERGAARD_HH__
 #define __WESTERGAARD_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "bem.hh"
 #include "grid_hermitian.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
-template <UInt dim>
+template <model_type type>
 class Westergaard : public Bem {
+  using trait = model_type_traits<type>;
 public:
   /// Constructor
-  Westergaard(Model * model,
-	      model_type type,
-	      const std::array<UInt, dim> & sizes);
+  Westergaard(Model * model);
   /// Destructor
-  virtual ~Westergaard();
+  virtual ~Westergaard() = default;
 
 public:
   /// Initialize influence coefficients
-  virtual void initInfluence();
+  void initInfluence() override;
   /// Solve Neumann problem
-  virtual void solveNeumann(const GridBase<Real> * neumann, GridBase<Real> * dirichlet) const;
+  void solveNeumann(const GridBase<Real> * neumann,
+		    GridBase<Real> * dirichlet) const override;
   /// Solve Dirichlet problem
-  virtual void solveDirichlet(const GridBase<Real> * dirichlet, GridBase<Real> * neumann) const;
+  void solveDirichlet(const GridBase<Real> * dirichlet,
+		      GridBase<Real> * neumann) const override;
 
 public:
-  GridHermitian<Real, dim> influence;
+  GridHermitian<Real, trait::dimension> influence;
 };
 
 __END_TAMAAS__
 
 #endif // __WESTERGAARD_HH__
diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp
index 8a6f2ad..174a714 100644
--- a/src/surface/isopowerlaw.cpp
+++ b/src/surface/isopowerlaw.cpp
@@ -1,108 +1,108 @@
 /**
  *
  * @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 "isopowerlaw.hh"
 #include "fftransform.hh"
 #include <map>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <UInt dim>
 void Isopowerlaw<dim>::computeFilter(GridHermitian<Real, dim> & filter_coefficients) const {
-  Grid<Real, dim> wavevectors(filter_coefficients.sizes(), dim);
-  FFTransform<Real, dim>::template computeFrequencies<true>(wavevectors);
+  auto wavevectors =
+    FFTransform<Real, dim>::template computeFrequencies<true>(filter_coefficients.sizes());
   auto end = wavevectors.end(dim);
 
 #pragma omp parallel for
   for (auto it = wavevectors.begin(dim) ; it < end ; ++it) {
     VectorProxy<Real> q(&(*it), dim);
     UInt i = wavevectors.getNbPoints() - (end - it);
     filter_coefficients(i) = (*this)(q);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::rmsHeights() const {
   return std::sqrt(M_PI * ( (hurst + 1)/hurst * q1*q1
 			    - 1./hurst * std::pow(q1, 2*(hurst+1)) * std::pow(q2, -2*hurst)
 			    - q0*q0 ));
 }
 
 /* -------------------------------------------------------------------------- */
 
 /** 
  * Analytical moments, cf. Yastrebov et al. (2015)
  * "From infinitesimal to full contact between rough surfaces: Evolution
  * of the contact area", appendix A
  */
 template <>
 std::vector<Real> Isopowerlaw<2>::moments() const {
   std::map<UInt, Real> T;
   T[0] = 2*M_PI; T[2] = M_PI; T[4] = 3 * M_PI / 4.;
   Real xi = q0 / q1;
   Real zeta = q2 / q1;
   
   std::vector<Real> moments;
 
   using namespace std;
   for (int q : {0, 2, 4}) {
     Real m = T[q] * pow(q1, q-2*hurst) * ((1-pow(xi, q+2))/(q+2) +
 					  (pow(zeta, q-2*hurst)-1)/(q-2*hurst));
     moments.push_back(m);
   }
 
   return moments;
 }
 
 template <>
 std::vector<Real> Isopowerlaw<1>::moments() const {
   TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces");
   return std::vector<Real>();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::alpha() const {
   std::vector<Real> m = moments();
   return m[0] * m[2] / (m[1]*m[1]);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::rmsSlopes() const {
   std::vector<Real> m = moments();
   return std::sqrt(M_PI * m[1] / 2);
 }
 
 
 template class Isopowerlaw<1>;
 template class Isopowerlaw<2>;
 
 __END_TAMAAS__