diff --git a/python/SConscript b/python/SConscript
index 6998e2e..1870dc4 100644
--- a/python/SConscript
+++ b/python/SConscript
@@ -1,91 +1,93 @@
 from __future__ import print_function
 
 from os.path import abspath, basename
 from distutils.version import StrictVersion
 from subprocess import check_output, STDOUT
 import re
 
 Import('main_env')
 
 env_swig = main_env.Clone(
     tools=['swig'],
     SHLIBPREFIX='',
 )
 
 if env_swig['SWIGVERSION'] is None:
     print("Could not find swig version")
     Exit(1)
 
 if StrictVersion(env_swig['SWIGVERSION']) < StrictVersion("3.0"):
     print("Swig version {} is not supported by tamaas".format(env_swig['SWIGVERSION']))
     Exit(1)
 
 # Check user's prefered version
 version_out = check_output([main_env['py_exec'], "--version"], stderr=STDOUT, universal_newlines=True).split(' ')
 
 # Matching version
 version_re = re.compile(r"^\d(\.\d+)*")
 version_string = version_re.match(version_out[-1]).group()
 # Verifying version
 version = StrictVersion(version_string)
 python3 = version >= StrictVersion("3.0")
 print("Building wrapper for python version {}".format(version))
 
 # Run code below to get versions from user preference and not scons
 versions_script = """
 from __future__ import print_function
 import distutils.sysconfig as dsys
 from numpy.distutils.misc_util import get_numpy_include_dirs as numpy_dirs
 
 print(dsys.get_python_inc())
 for d in numpy_dirs():
     print(d)"""
 
 includes = check_output([main_env['py_exec'], "-c", versions_script],
                         universal_newlines=True).split('\n')
 
 env_swig.AppendUnique(CPPPATH=includes)
 env_swig.AppendUnique(SWIGPATH=['$CPPPATH'])
 env_swig.AppendUnique(SWIGFLAGS=['-python', '-c++'])
 
 if python3:
     env_swig.AppendUnique(SWIGFLAGS=['-py3'])
 
 verbose = main_env['verbose']
 if not verbose:
     env_swig.AppendUnique(SWIGFLAGS=['-w325', '-w315'])
 if basename(env_swig['CXX']) != "clang++":  # for some warning in wrapper code
     env_swig.AppendUnique(CXXFLAGS=['-Wno-maybe-uninitialized'])
 else:
     env_swig.AppendUnique(CXXFLAGS=['-Wno-uninitialized',
                                     '-Wno-dynamic-class-memaccess'])
 
 # env_swig.SharedLibrary(
 #     target='_tamaas.so',
 #     source=['tamaas.i'],
 #     LIBS=['Tamaas'],
 #     RPATH=[abspath('../src')]
 # )
 
 # Pybind11 wrapper
 env_pybind = main_env.Clone(SHLIBPREFIX='')
 
 includes += [Dir('#third-party/pybind11/include')]
 
 extension = check_output([main_env['py_exec'] + '-config', "--extension-suffix"],
                          universal_newlines=True).split('\n')[0]
 env_pybind.AppendUnique(CPPPATH=includes)
 
 pybind_sources = Split("""
 tamaas_module.cpp
 wrap/core.cpp
 wrap/percolation.cpp
 wrap/surface.cpp
+wrap/model.cpp
+wrap/solvers.cpp
 """)
 
 env_pybind.SharedLibrary(
     target='tamaas' + extension,
     source=pybind_sources,
     LIBS=['Tamaas'],
     RPATH=[abspath('../src')]
 )
diff --git a/python/wrap.hh b/python/cast.hh
similarity index 64%
copy from python/wrap.hh
copy to python/cast.hh
index a137b7a..fa1ddcb 100644
--- a/python/wrap.hh
+++ b/python/cast.hh
@@ -1,207 +1,195 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
-#ifndef __WRAP_HH__
-#define __WRAP_HH__
+#ifndef __CAST_HH__
+#define __CAST_HH__
 /* -------------------------------------------------------------------------- */
-#include "grid.hh"
+#include "grid_base.hh"
 #include "surface.hh"
-#include "tamaas.hh"
+#include "numpy.hh"
 /* -------------------------------------------------------------------------- */
-#include <pybind11/pybind11.h>
-// keep order
 #include <pybind11/cast.h>
-#include <pybind11/numpy.h>
-/* -------------------------------------------------------------------------- */
-
-namespace tamaas {
-namespace wrap {
-namespace py = pybind11;
-
-template <typename T>
-using numpy = py::array_t<T, py::array::c_style | py::array::forcecast>;
-
-/// Class wrapping a numpy array
-template <class Parent>
-class GridNumpy : public Parent {
-public:
-  GridNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
-    auto* array_shape = buffer.shape();
-
-    if (buffer.ndim() > Parent::dimension + 1 ||
-        buffer.ndim() < Parent::dimension)
-      TAMAAS_EXCEPTION(
-          "Numpy array dimension do not match expected grid dimensions");
-
-    if (buffer.ndim() == Parent::dimension + 1)
-      this->nb_components = array_shape[buffer.ndim() - 1];
-
-    std::copy_n(array_shape, Parent::dimension, this->n.begin());
-    this->computeStrides();
-    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
-  }
-};
-
-/// Class wrapping a numpy array
-template <class Parent>
-class SurfaceNumpy : public Parent {
-public:
-  SurfaceNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
-    auto* array_shape = buffer.shape();
-
-    if (buffer.ndim() != 2)
-      TAMAAS_EXCEPTION(
-          "Numpy array dimension do not match expected surface dimensions");
-
-    std::copy_n(array_shape, Parent::dimension, this->n.begin());
-    this->computeStrides();
-    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
-  }
-};
-
-}  // namespace wrap
-}  // namespace tamaas
-
+#include <boost/preprocessor/seq.hpp>
 /* -------------------------------------------------------------------------- */
 
 namespace pybind11 {
 namespace detail {
 
 /**
  * Type caster for grid classes
  * inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
  * and pybind11/eigen.h
  */
 template <template <typename, unsigned int> class G, typename T,
           unsigned int dim>
 struct type_caster<G<T, dim>> {
   using type = G<T, dim>;
   template <typename U>
   using array_type = array_t<U, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("GridWrap<T, dim>"));
 
   /**
    * Conversion part 1 (Python->C++): convert a PyObject into a grid
    * instance or return false upon failure. The second argument
    * indicates whether implicit conversions should be applied.
    */
   bool load(handle src, bool convert) {
     if (!convert && !array_type<T>::check_(src))
       return false;
 
     auto buf = array_type<T>::ensure(src);
     value.move(tamaas::wrap::GridNumpy<G<T, dim>>(buf));
 
     return true;
   }
 
   /**
    * Conversion part 2 (C++ -> Python): convert a grid instance into
    * a Python object. The second and third arguments are used to
    * indicate the return value policy and parent object (for
    * ``return_value_policy::reference_internal``) and are generally
    * ignored by implicit casters.
    */
   static handle cast(const type& src, return_value_policy /* policy */,
                      handle /*parent*/) {
     if (src.getNbComponents() == 1) {
       auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
       return a.release();
     } else {
       std::array<unsigned int, dim> sizes;
       std::copy_n(src.sizes().begin(), dim, sizes.begin());
       sizes.back() = src.getNbComponents();
       auto a = array_type<T>(sizes, src.getInternalData(), none());
       return a.release();
     }
   }
 };
 
+/**
+ * Type caster for GridBase classes
+ */
+template <typename T>
+struct type_caster<tamaas::GridBase<T>> {
+  using type = tamaas::GridBase<T>;
+  template <typename U>
+  using array_type = array_t<U, array::c_style | array::forcecast>;
+
+public:
+  PYBIND11_TYPE_CASTER(type, _("GridBaseWrap<T>"));
+
+  /**
+   * Conversion part 1 (Python->C++): convert a PyObject into a grid
+   * instance or return false upon failure. The second argument
+   * indicates whether implicit conversions should be applied.
+   */
+  bool load(handle src, bool convert) {
+    if (!convert && !array_type<T>::check_(src))
+      return false;
+
+    auto buf = array_type<T>::ensure(src);
+    value.move(tamaas::wrap::GridBaseNumpy<T>(buf));
+
+    return true;
+  }
+
+  /**
+   * Conversion part 2 (C++ -> Python): convert a grid instance into
+   * a Python object. The second and third arguments are used to
+   * indicate the return value policy and parent object (for
+   * ``return_value_policy::reference_internal``) and are generally
+   * ignored by implicit casters.
+   */
+  static handle cast(const type& src, return_value_policy /* policy */,
+                     handle /*parent*/) {
+#define GRID_BASE_CASE(unused1, unused2, dim)                                  \
+  case dim: {                                                                  \
+    const tamaas::Grid<T, dim>& conv = static_cast<decltype(conv)>(src);       \
+    if (src.getNbComponents() != 1) {                                          \
+      std::array<unsigned int, dim + 1> sizes;                                 \
+      std::copy_n(conv.sizes().begin(), dim, sizes.begin());                   \
+      sizes.back() = src.getNbComponents();                                    \
+      return array_type<T>(sizes, src.getInternalData(), none()).release();    \
+    } else {                                                                   \
+      return array_type<T>(conv.sizes(), src.getInternalData(), none())        \
+          .release();                                                          \
+    }                                                                          \
+  }
+
+    switch (src.getDimension()) {
+      BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
+    default:
+      return nullptr;
+    }
+#undef GRID_BASE_CASE
+  }
+};
+
+/**
+ * Type caster for surface class
+ */
 template <typename T>
 struct type_caster<tamaas::Surface<T>> {
   using type = tamaas::Surface<T>;
   template <typename U>
   using array_type = array_t<U, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("SurfaceWrap<T, dim>"));
 
   /**
    * Conversion part 1 (Python->C++): convert a PyObject into a grid
    * instance or return false upon failure. The second argument
    * indicates whether implicit conversions should be applied.
    */
   bool load(handle src, bool convert) {
     if (!convert && !array_type<T>::check_(src))
       return false;
 
     auto buf = array_type<T>::ensure(src);
     value.move(tamaas::wrap::SurfaceNumpy<type>(buf));
 
     return true;
   }
 
   /**
    * Conversion part 2 (C++ -> Python): convert a grid instance into
    * a Python object. The second and third arguments are used to
    * indicate the return value policy and parent object (for
    * ``return_value_policy::reference_internal``) and are generally
    * ignored by implicit casters.
    */
   static handle cast(const type& src, return_value_policy /* policy */,
                      handle /*parent*/) {
     auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
     return a.release();
   }
 };
 }  // namespace detail
 }  // namespace pybind11
 
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-namespace wrap {
-
-inline std::string makeDimensionName(const std::string& name, UInt dim) {
-  std::stringstream str;
-  str << name << dim << "D";
-  return str.str();
-}
-
-/* -------------------------------------------------------------------------- */
-void wrapCore(py::module& mod);
-void wrapPercolation(py::module& mod);
-void wrapSurface(py::module& mod);
-/* -------------------------------------------------------------------------- */
-
-}  // namespace wrap
-
-__END_TAMAAS__
-
-#endif  // __WRAP_HH__
+#endif // __CAST_HH__
diff --git a/python/numpy.hh b/python/numpy.hh
new file mode 100644
index 0000000..1a889cc
--- /dev/null
+++ b/python/numpy.hh
@@ -0,0 +1,101 @@
+/**
+ * @file
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  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/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef __NUMPY_HH__
+#define __NUMPY_HH__
+/* -------------------------------------------------------------------------- */
+#include "tamaas.hh"
+#include "grid_base.hh"
+/* -------------------------------------------------------------------------- */
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+__BEGIN_TAMAAS__
+
+namespace wrap {
+namespace py = pybind11;
+
+template <typename T>
+using numpy = py::array_t<T, py::array::c_style | py::array::forcecast>;
+
+template <typename T>
+class GridBaseNumpy : public GridBase<T> {
+public:
+  GridBaseNumpy(numpy<T>& buffer) : GridBase<T>() {
+    this->nb_components = 1;
+
+    if (buffer.ndim() == 2)
+      this->nb_components = buffer.shape()[1];
+    else if (buffer.ndim() > 2)
+      TAMAAS_EXCEPTION("Dimension of numpy array not compatible with GridBase");
+
+    this->data.wrapMemory(buffer.mutable_data(), buffer.size());
+  }
+};
+
+/// Class wrapping a numpy array
+template <class Parent>
+class GridNumpy : public Parent {
+public:
+  GridNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
+    auto* array_shape = buffer.shape();
+
+    if (buffer.ndim() > Parent::dimension + 1 ||
+        buffer.ndim() < Parent::dimension)
+      TAMAAS_EXCEPTION(
+          "Numpy array dimension do not match expected grid dimensions");
+
+    if (buffer.ndim() == Parent::dimension + 1)
+      this->nb_components = array_shape[buffer.ndim() - 1];
+
+    std::copy_n(array_shape, Parent::dimension, this->n.begin());
+    this->computeStrides();
+    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
+  }
+};
+
+/// Class wrapping a numpy array
+template <class Parent>
+class SurfaceNumpy : public Parent {
+public:
+  SurfaceNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
+    auto* array_shape = buffer.shape();
+
+    if (buffer.ndim() != 2)
+      TAMAAS_EXCEPTION(
+          "Numpy array dimension do not match expected surface dimensions");
+
+    std::copy_n(array_shape, Parent::dimension, this->n.begin());
+    this->computeStrides();
+    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
+  }
+};
+
+}  // namespace wrap
+
+__END_TAMAAS__
+
+#endif // __NUMPY_HH__
diff --git a/python/tamaas_module.cpp b/python/tamaas_module.cpp
index b5aec9f..c56b698 100644
--- a/python/tamaas_module.cpp
+++ b/python/tamaas_module.cpp
@@ -1,51 +1,53 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  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 "tamaas.hh"
 #include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace py = pybind11;
 
 PYBIND11_MODULE(tamaas, mod) {
   mod.doc() = "tamaas module for python";
 
   // Wrapping the base methods
   mod.def("initialize", &initialize, py::arg("num_threads") = 0,
           "Initialize tamaas with desired number of threads");
   mod.def("finalize", &finalize, "Final cleanup");
 
   // Wrapping tamaas components
   wrap::wrapCore(mod);
   wrap::wrapPercolation(mod);
   wrap::wrapSurface(mod);
+  wrap::wrapModel(mod);
+  wrap::wrapSolvers(mod);
 }
 
 __END_TAMAAS__
diff --git a/python/wrap.hh b/python/wrap.hh
index a137b7a..9aad699 100644
--- a/python/wrap.hh
+++ b/python/wrap.hh
@@ -1,207 +1,61 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __WRAP_HH__
 #define __WRAP_HH__
 /* -------------------------------------------------------------------------- */
-#include "grid.hh"
-#include "surface.hh"
 #include "tamaas.hh"
+#include "numpy.hh"
+#include "cast.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/pybind11.h>
-// keep order
-#include <pybind11/cast.h>
-#include <pybind11/numpy.h>
-/* -------------------------------------------------------------------------- */
-
-namespace tamaas {
-namespace wrap {
-namespace py = pybind11;
-
-template <typename T>
-using numpy = py::array_t<T, py::array::c_style | py::array::forcecast>;
-
-/// Class wrapping a numpy array
-template <class Parent>
-class GridNumpy : public Parent {
-public:
-  GridNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
-    auto* array_shape = buffer.shape();
-
-    if (buffer.ndim() > Parent::dimension + 1 ||
-        buffer.ndim() < Parent::dimension)
-      TAMAAS_EXCEPTION(
-          "Numpy array dimension do not match expected grid dimensions");
-
-    if (buffer.ndim() == Parent::dimension + 1)
-      this->nb_components = array_shape[buffer.ndim() - 1];
-
-    std::copy_n(array_shape, Parent::dimension, this->n.begin());
-    this->computeStrides();
-    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
-  }
-};
-
-/// Class wrapping a numpy array
-template <class Parent>
-class SurfaceNumpy : public Parent {
-public:
-  SurfaceNumpy(numpy<typename Parent::value_type>& buffer) : Parent() {
-    auto* array_shape = buffer.shape();
-
-    if (buffer.ndim() != 2)
-      TAMAAS_EXCEPTION(
-          "Numpy array dimension do not match expected surface dimensions");
-
-    std::copy_n(array_shape, Parent::dimension, this->n.begin());
-    this->computeStrides();
-    this->data.wrapMemory(buffer.mutable_data(), this->computeSize());
-  }
-};
-
-}  // namespace wrap
-}  // namespace tamaas
-
-/* -------------------------------------------------------------------------- */
-
-namespace pybind11 {
-namespace detail {
-
-/**
- * Type caster for grid classes
- * inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
- * and pybind11/eigen.h
- */
-template <template <typename, unsigned int> class G, typename T,
-          unsigned int dim>
-struct type_caster<G<T, dim>> {
-  using type = G<T, dim>;
-  template <typename U>
-  using array_type = array_t<U, array::c_style | array::forcecast>;
-
-public:
-  PYBIND11_TYPE_CASTER(type, _("GridWrap<T, dim>"));
-
-  /**
-   * Conversion part 1 (Python->C++): convert a PyObject into a grid
-   * instance or return false upon failure. The second argument
-   * indicates whether implicit conversions should be applied.
-   */
-  bool load(handle src, bool convert) {
-    if (!convert && !array_type<T>::check_(src))
-      return false;
-
-    auto buf = array_type<T>::ensure(src);
-    value.move(tamaas::wrap::GridNumpy<G<T, dim>>(buf));
-
-    return true;
-  }
-
-  /**
-   * Conversion part 2 (C++ -> Python): convert a grid instance into
-   * a Python object. The second and third arguments are used to
-   * indicate the return value policy and parent object (for
-   * ``return_value_policy::reference_internal``) and are generally
-   * ignored by implicit casters.
-   */
-  static handle cast(const type& src, return_value_policy /* policy */,
-                     handle /*parent*/) {
-    if (src.getNbComponents() == 1) {
-      auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
-      return a.release();
-    } else {
-      std::array<unsigned int, dim> sizes;
-      std::copy_n(src.sizes().begin(), dim, sizes.begin());
-      sizes.back() = src.getNbComponents();
-      auto a = array_type<T>(sizes, src.getInternalData(), none());
-      return a.release();
-    }
-  }
-};
-
-template <typename T>
-struct type_caster<tamaas::Surface<T>> {
-  using type = tamaas::Surface<T>;
-  template <typename U>
-  using array_type = array_t<U, array::c_style | array::forcecast>;
-
-public:
-  PYBIND11_TYPE_CASTER(type, _("SurfaceWrap<T, dim>"));
-
-  /**
-   * Conversion part 1 (Python->C++): convert a PyObject into a grid
-   * instance or return false upon failure. The second argument
-   * indicates whether implicit conversions should be applied.
-   */
-  bool load(handle src, bool convert) {
-    if (!convert && !array_type<T>::check_(src))
-      return false;
-
-    auto buf = array_type<T>::ensure(src);
-    value.move(tamaas::wrap::SurfaceNumpy<type>(buf));
-
-    return true;
-  }
-
-  /**
-   * Conversion part 2 (C++ -> Python): convert a grid instance into
-   * a Python object. The second and third arguments are used to
-   * indicate the return value policy and parent object (for
-   * ``return_value_policy::reference_internal``) and are generally
-   * ignored by implicit casters.
-   */
-  static handle cast(const type& src, return_value_policy /* policy */,
-                     handle /*parent*/) {
-    auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
-    return a.release();
-  }
-};
-}  // namespace detail
-}  // namespace pybind11
-
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace wrap {
 
+namespace py = pybind11;
+
 inline std::string makeDimensionName(const std::string& name, UInt dim) {
   std::stringstream str;
   str << name << dim << "D";
   return str.str();
 }
 
 /* -------------------------------------------------------------------------- */
 void wrapCore(py::module& mod);
 void wrapPercolation(py::module& mod);
 void wrapSurface(py::module& mod);
+void wrapModel(py::module& mod);
+void wrapSolvers(py::module& mod);
 /* -------------------------------------------------------------------------- */
 
 }  // namespace wrap
 
 __END_TAMAAS__
 
 #endif  // __WRAP_HH__
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
new file mode 100644
index 0000000..8fd97eb
--- /dev/null
+++ b/python/wrap/model.cpp
@@ -0,0 +1,76 @@
+/**
+ * @file
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  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 "wrap.hh"
+#include "model.hh"
+#include "model_factory.hh"
+#include <pybind11/stl.h>
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_TAMAAS__
+
+namespace wrap {
+
+using namespace py::literals;
+
+void wrapModelClass(py::module& mod) {
+  py::enum_<model_type>(mod, "model_type")
+    .value("basic_1d", model_type::basic_1d)
+    .value("basic_2d", model_type::basic_1d)
+    .value("surface_1d", model_type::surface_1d)
+    .value("surface_2d", model_type::surface_1d)
+    .value("volume_1d", model_type::volume_1d)
+    .value("volume_2d", model_type::volume_1d)
+    .export_values();
+
+  py::class_<Model>(mod, "Model")
+      .def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
+      .def("getHertzModulus", &Model::getHertzModulus)
+      .def("getYoungModulus", &Model::getYoungModulus)
+      .def("getPoissonRatio", &Model::getPoissonRatio)
+      .def("getTraction", (GridBase<Real> & (Model::*)()) & Model::getTraction)
+      .def("getDisplacement",
+           (GridBase<Real> & (Model::*)()) & Model::getDisplacement)
+      .def("getSystemSize", &Model::getSystemSize)
+      .def("getDiscretization", &Model::getDiscretization)
+      .def("solveNeumann", &Model::solveNeumann)
+      .def("solveDirichlet", &Model::solveDirichlet);
+}
+
+void wrapModelFactory(py::module& mod) {
+  py::class_<ModelFactory>(mod, "ModelFactory")
+      .def_static("createModel", &ModelFactory::createModel, "model_type"_a,
+                  "system_size"_a, "discretization"_a);
+}
+
+void wrapModel(py::module& mod) {
+  wrapModelClass(mod);
+  wrapModelFactory(mod);
+}
+
+}  // namespace wrap
+
+__END_TAMAAS__
diff --git a/python/tamaas_module.cpp b/python/wrap/solvers.cpp
similarity index 51%
copy from python/tamaas_module.cpp
copy to python/wrap/solvers.cpp
index b5aec9f..75ae669 100644
--- a/python/tamaas_module.cpp
+++ b/python/wrap/solvers.cpp
@@ -1,51 +1,63 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  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 "tamaas.hh"
 #include "wrap.hh"
-/* -------------------------------------------------------------------------- */
-#include <pybind11/pybind11.h>
+#include "contact_solver.hh"
+#include "polonsky_keer_rey.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
-namespace py = pybind11;
+namespace wrap {
 
-PYBIND11_MODULE(tamaas, mod) {
-  mod.doc() = "tamaas module for python";
+/* -------------------------------------------------------------------------- */
+using namespace py::literals;
 
-  // Wrapping the base methods
-  mod.def("initialize", &initialize, py::arg("num_threads") = 0,
-          "Initialize tamaas with desired number of threads");
-  mod.def("finalize", &finalize, "Final cleanup");
+/* -------------------------------------------------------------------------- */
+void wrapSolvers(py::module& mod) {
+  py::class_<ContactSolver>(mod, "ContactSolver")
+      .def(py::init<Model&, const GridBase<Real>&, Real>(), "model"_a,
+           "surface"_a, "tolerance"_a)
+      .def("setMaxIterations", &ContactSolver::setMaxIterations, "max_iter"_a)
+      .def("setDumpFrequency", &ContactSolver::setDumpFrequency, "dump_freq"_a);
 
-  // Wrapping tamaas components
-  wrap::wrapCore(mod);
-  wrap::wrapPercolation(mod);
-  wrap::wrapSurface(mod);
+  py::class_<PolonskyKeerRey, ContactSolver> pkr(mod, "PolonskyKeerRey");
+  pkr.def(py::init<Model&, const GridBase<Real>&, Real, PolonskyKeerRey::type,
+                   PolonskyKeerRey::type>(),
+          "model"_a, "surface"_a, "tolerance"_a, "primal_type"_a,
+          "constraint_type"_a)
+      .def("solve", &PolonskyKeerRey::solve, "target"_a)
+      .def("computeError", &PolonskyKeerRey::computeError);
+
+  py::enum_<PolonskyKeerRey::type>(pkr, "type")
+      .value("gap", PolonskyKeerRey::gap)
+      .value("pressure", PolonskyKeerRey::pressure)
+      .export_values();
 }
 
+}  // namespace wrap
+
 __END_TAMAAS__
diff --git a/python/wrap/surface.cpp b/python/wrap/surface.cpp
index a4419a5..d1fb3c3 100644
--- a/python/wrap/surface.cpp
+++ b/python/wrap/surface.cpp
@@ -1,93 +1,94 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  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 "wrap.hh"
 #include "isopowerlaw.hh"
 #include "surface_generator_filter.hh"
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 namespace wrap {
 
 using namespace py::literals;
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapFilter(py::module& mod) {
   auto name = makeDimensionName("Filter", dim);
   py::class_<Filter<dim>>(mod, name.c_str());
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapIsopowerlaw(py::module& mod) {
   std::string name = makeDimensionName("Isopowerlaw", dim);
   py::class_<Isopowerlaw<dim>, Filter<dim>>(mod, name.c_str())
       .def(py::init<>())
       .def_property("q0", &Isopowerlaw<dim>::getQ0, &Isopowerlaw<dim>::setQ0)
       .def_property("q1", &Isopowerlaw<dim>::getQ1, &Isopowerlaw<dim>::setQ1)
       .def_property("q2", &Isopowerlaw<dim>::getQ2, &Isopowerlaw<dim>::setQ2)
       .def_property("hurst", &Isopowerlaw<dim>::getHurst,
                     &Isopowerlaw<dim>::setHurst)
       .def("rmsHeights", &Isopowerlaw<dim>::rmsHeights)
       .def("moments", &Isopowerlaw<dim>::moments)
       .def("alpha", &Isopowerlaw<dim>::alpha)
       .def("rmsSlopes", &Isopowerlaw<dim>::rmsSlopes);
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapSurfaceGeneratorFilter(py::module& mod) {
   std::string name = makeDimensionName("SurfaceGeneratorFilter", dim);
   py::class_<SurfaceGeneratorFilter<dim>>(mod, name.c_str())
       .def(py::init<>())
       .def("buildSurface", &SurfaceGeneratorFilter<dim>::buildSurface,
            "Build rough surface")
       .def("setFilter", &SurfaceGeneratorFilter<dim>::setFilter,
            "Set PSD filter", "filter"_a)
       .def("setSizes",
            (void (SurfaceGeneratorFilter<dim>::*)(
                const std::array<UInt, dim>&)) &
-               SurfaceGeneratorFilter<dim>::setSizes)
+               SurfaceGeneratorFilter<dim>::setSizes,
+           "surface_sizes"_a)
       .def_property("random_seed", &SurfaceGeneratorFilter<dim>::getRandomSeed,
                     &SurfaceGeneratorFilter<dim>::setRandomSeed);
 }
 
 /* -------------------------------------------------------------------------- */
 void wrapSurface(py::module& mod) {
   wrapFilter<1>(mod);
   wrapFilter<2>(mod);
   wrapIsopowerlaw<1>(mod);
   wrapIsopowerlaw<2>(mod);
   wrapSurfaceGeneratorFilter<1>(mod);
   wrapSurfaceGeneratorFilter<2>(mod);
 }
 
 }  // namespace wrap
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
diff --git a/src/core/grid.hh b/src/core/grid.hh
index 4a022a6..26403af 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,198 +1,176 @@
 /**
  * @file
  *
  * @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 "grid_base.hh"
 #include "tamaas.hh"
 #include <array>
 #include <numeric>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 __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 */
   /* ------------------------------------------------------------------------ */
   using value_type = T;
   static constexpr UInt dimension = dim;
 
   /* ------------------------------------------------------------------------ */
   /* Constructors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Constructor by default (empty array)
   Grid();
   /// Constructor
   Grid(const std::array<UInt, dim>& n, UInt nb_components);
   /// Wrap on given data
   Grid(const std::array<UInt, dim>& n, UInt nb_components, value_type * data);
   /// 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);
   /// Destructor
   ~Grid() override = default;
 
 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 override;
+  inline UInt computeSize() const;
+  /// Get grid dimension
+  inline UInt getDimension() const override { return dim; }
   /// 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;
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy operators */
   /* ------------------------------------------------------------------------ */
 public:
   // = operator
   Grid& operator=(const Grid& other);
   // = operator (not const input, otherwise gcc is confused)
   Grid& operator=(Grid& other);
   // = operator (move)
   Grid& operator=(Grid&& other);
   // = operator for scalar
   Grid& operator=(T val) {
     GridBase<T>::operator=(val);
     return *this;
   }
 
   // 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;            ///< shape of grid: size per dimension
   std::array<UInt, dim + 1> strides;  ///< strides for access
 };
 
 __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 52ed98d..e692c0c 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,290 +1,303 @@
 /**
  * @file
  *
  * @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 "array.hh"
 #include "iterator.hh"
 #include "loop.hh"
 #include "tamaas.hh"
 #include <cstddef>
 #include <limits>
 #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) {}
 };
 }  // namespace view
 
 /**
  * @brief Dimension agnostic grid with components stored per points
  */
 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 */
   /* ------------------------------------------------------------------------ */
 public:
   using iterator = iterator_::iterator<T>;
   using const_iterator = iterator_::iterator<const T>;
 
 public:
-  /// Compute size
-  virtual UInt computeSize() const = 0;
+  /// Get grid dimension
+  virtual UInt getDimension() const { return 1; }
   /// Get internal data pointer (const)
   const T* getInternalData() const { return this->data.data(); }
   /// Get internal data pointer (non-const)
   T* getInternalData() { return this->data.data(); }
   /// 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();
+    return this->dataSize() / 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;
+  virtual iterator begin(UInt n = 1) {
+    return iterator(this->getInternalData() + this->offset, 0, n);
+  }
+
+  virtual iterator end(UInt n = 1) {
+    return iterator(this->getInternalData() + this->offset, this->dataSize(),
+                    n);
+  }
+
+  virtual const_iterator begin(UInt n = 1) const {
+    return const_iterator(this->getInternalData() + this->offset, 0, n);
+  }
+
+  virtual const_iterator end(UInt n = 1) const {
+    return const_iterator(this->getInternalData() + this->offset,
+                          this->dataSize(), n);
+  }
 
   /* ------------------------------------------------------------------------ */
   /* Operators */
   /* ------------------------------------------------------------------------ */
   inline T& operator()(UInt i) { return this->data[i]; }
 
   inline const T& operator()(UInt i) const { return this->data[i]; }
 
 #define VEC_OPERATOR(op)                                                       \
   template <typename T1>                                                       \
   GridBase& operator op(const GridBase<T1>& other)
 
   VEC_OPERATOR(+=);
   VEC_OPERATOR(*=);
   VEC_OPERATOR(-=);
   VEC_OPERATOR(/=);
 #undef VEC_OPERATOR
 
   template <typename T1>
   T dot(const GridBase<T1>& other) const;
 
 #define SCALAR_OPERATOR(op) GridBase& operator op(const T& e)
 
   SCALAR_OPERATOR(+=);
   SCALAR_OPERATOR(*=);
   SCALAR_OPERATOR(-=);
   SCALAR_OPERATOR(/=);
   SCALAR_OPERATOR(=);
 #undef SCALAR_OPERATOR
 
   GridBase& operator=(const GridBase& o);
 
   /* ------------------------------------------------------------------------ */
   /* Min/max */
   /* ------------------------------------------------------------------------ */
   T min() const;
   T max() const;
   T sum() const;
   T mean() const { return sum() / static_cast<T>(dataSize()); }
   T var() const;
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy */
   /* ------------------------------------------------------------------------ */
 
 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");
 
   auto begin_it(begin());
   auto end_it(end());  // silencing nvcc warnings
 
 #pragma omp parallel for
   for (auto it = begin(); it < end_it; ++it) {
     UInt i = it - begin_it;
     *it = vec.data[i % this->nb_components];
   }
 }
 
   /* --------------------------------------------------------------------------
    */
 
 #define SCALAR_OPERATOR_IMPL(op)                                               \
   template <typename T>                                                        \
   inline GridBase<T>& GridBase<T>::operator op(const T& e) {                   \
     Loop::loop([e] CUDA_LAMBDA(T& val) { val op e; }, *this);                  \
     return *this;                                                              \
   }
 
 SCALAR_OPERATOR_IMPL(+=);
 SCALAR_OPERATOR_IMPL(*=);
 SCALAR_OPERATOR_IMPL(-=);
 SCALAR_OPERATOR_IMPL(/=);
 SCALAR_OPERATOR_IMPL(=);
 
 #undef SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 #define VEC_OPERATOR_IMPL(op)                                                  \
   template <typename T>                                                        \
   template <typename T1>                                                       \
   inline GridBase<T>& GridBase<T>::operator op(const GridBase<T1>& other) {    \
-    TAMAAS_ASSERT(other.computeSize() == this->computeSize(),                  \
+    TAMAAS_ASSERT(other.dataSize() == this->dataSize(),                        \
                   "surface size does not match");                              \
     Loop::loop([] CUDA_LAMBDA(T& x, const T1& y) { x op y; }, *this, other);   \
     return *this;                                                              \
   }
 
 VEC_OPERATOR_IMPL(+=);
 VEC_OPERATOR_IMPL(-=);
 VEC_OPERATOR_IMPL(*=);
 VEC_OPERATOR_IMPL(/=);
 #undef VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 GridBase<T>& GridBase<T>::operator=(const GridBase& o) {
   this->copy(o);
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 using op = operation;
 
 template <typename T>
 inline T GridBase<T>::min() const {
   // const auto id = [] CUDA_LAMBDA (const T&x){return x;};
   return Loop::reduce<op::min>([] CUDA_LAMBDA(const T& x) { return x; }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::max() const {
   return Loop::reduce<op::max>([] CUDA_LAMBDA(const T& x) { return x; }, *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::sum() const {
   return Loop::reduce<op::plus>([] CUDA_LAMBDA(const T& x) { return x; },
                                 *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::var() const {
   const T mu = mean();
   const T var = Loop::reduce<op::plus>(
       [mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
   return var / (dataSize() - 1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <typename T1>
 T GridBase<T>::dot(const GridBase<T1>& other) const {
   return Loop::reduce<op::plus>(
       [] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
 }
 
 __END_TAMAAS__
 
 #endif  // __GRID_BASE_HH__