diff --git a/src/core/array.hh b/src/core/array.hh
index 56d43eb..c803920 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,163 +1,156 @@
 /**
  * @file
  *
  * @author Guillaume Anciaux <guillaume.anciaux@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 __ARRAY__HH__
 #define __ARRAY__HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
+#include <thrust/fill.h>
+#include <thrust/copy.h>
 #include <memory>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Generic storage class with wrapping capacities
  */
 template <typename T>
 class Array {
 public:
   /// Default
   Array() = default;
   /// Empty array of given size
   Array(UInt size) : Array() { this->resize(size); }
   /// Copy constructor (deep)
   Array(const Array& v) : Array() {
     this->resize(v._size);
     std::copy(v._data, v._data + v._size, this->_data);
   }
   /// Move constructor (transfers data ownership)
   Array(Array&& v) : Array() {
     _data = std::exchange(v._data, nullptr);
     _size = std::exchange(v._size, 0);
     wrapped = std::exchange(v.wrapped, false);
   }
   /// Wrap array on data
   Array(T* data, UInt size) : Array() { this->wrapMemory(data, size); }
   /// Destructor
   virtual ~Array() {
     if (wrapped == false)
       _alloc.deallocate(_data, _size);
   }
 
 public:
   /// Copy operator
   Array& operator=(const Array& v) {
     this->wrapped = false;
     if (v.size() != this->size())
       this->resize(v.size());
-#pragma omp parallel for
-    for (UInt i = 0; i < this->_size; ++i) {
-      _data[i] = v._data[i];
-    }
+    thrust::copy(v._data, v._data + _size, _data);
     return *this;
   }
 
   /// Move operator
   Array& operator=(Array&& v) {
     if (this != &v) {
       if (!wrapped)
         _alloc.deallocate(_data, _size);
-      _data = v._data;
-      _size = v._size;
-      wrapped = v.wrapped;
-      v._data = nullptr;
-      v._size = 0;
-      v.wrapped = false;
+      _data = std::exchange(v._data, nullptr);
+      _size = std::exchange(v._size, 0);
+      wrapped = std::exchange(v.wrapped, false);
     }
     return *this;
   }
 
   /// Wrap a memory pointer
   void wrapMemory(T* data, UInt size) {
     this->_data = data;
     this->_size = size;
     this->wrapped = true;
   }
 
   /// Assign a sigle value to the array
   void assign(UInt size, const T& value) {
     this->resize(size);
-#pragma omp parallel for
-    for (UInt i = 0; i < size; i++) {
-      _data[i] = value;
-    }
+    thrust::fill(_data, _data + size, value);
   }
 
   void setWrapped(bool w) { wrapped = w; }
 
   /// Data pointer access (const)
   const T* data() const { return _data; }
   /// Data pointer access (non-const)
   T* data() { return _data; }
   /// Data pointer setter
   void setData(T* new_ptr) { _data = new_ptr; }
   /// Resize array
   void resize(UInt size) {
     if (wrapped == true)
       TAMAAS_EXCEPTION("cannot resize wrapped array");
 
     // Erase array
     if (size == 0) {
       _alloc.deallocate(_data, _size);
       this->_data = nullptr;
       this->_size = 0;
       return;
     }
 
     // Do nothing
     if (size == this->_size)
       return;
 
     // Allocate new data
     _alloc.deallocate(_data, _size);
 
     this->_data = _alloc.allocate(size);
     this->_size = size;
   }
   /// Access operator
   inline T& operator[](UInt i) {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Access operator (const)
   inline const T& operator[](UInt i) const {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Get size of array
   inline UInt size() const { return _size; }
 
 private:
   T* _data = nullptr;
   UInt _size = 0;
   bool wrapped = false;
   DEFAULT_ALLOCATOR _alloc = DEFAULT_ALLOCATOR();
 };
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif /* __ARRAY__HH__ */
diff --git a/src/core/grid.hh b/src/core/grid.hh
index 153f00c..8460aa0 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,194 +1,196 @@
 /**
  * @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);
   /// 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;
   /// 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/types.hh b/src/core/types.hh
index 40ee71e..d2480ab 100644
--- a/src/core/types.hh
+++ b/src/core/types.hh
@@ -1,120 +1,120 @@
 /**
  * @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 __TYPES_HH__
 #define __TYPES_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 #include <cstddef>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T, UInt dim>
 class Proxy : public Grid<T, dim> {
 public:
   /// Default constructor
   Proxy() : Grid<T, dim>() {}
   /// Copy constructor
   Proxy(Proxy& o) : Grid<T, dim>() {
     this->n = o.n;
     this->nb_components = o.nb_components;
     this->offset = o.offset;
     this->strides = o.strides;
     this->data.wrapMemory(o.getInternalData(), o.dataSize());
   }
   /// Destructor
-  virtual ~Proxy() {}
+  ~Proxy() override = default;
 
 public:
   /// Set pointer
   void setPointer(T* ptr) { this->data.setData(ptr); }
   /// Get increment
   ptrdiff_t increment() const { return this->data.size(); }
   ///
   Proxy& operator=(const Proxy& o) {
     this->copy(o);
     return *this;
   }
   Proxy& operator=(T val) {
     Grid<T, dim>::operator=(val);
     return *this;
   }
   template <typename T1>
   Proxy& operator=(const Grid<T1, dim>& o) {
     this->copy(o);
     return *this;
   }
 
   template <typename T1>
   void copy(const Grid<T1, dim>& o) {
     Grid<T, dim>::copy(o);
     this->data.setWrapped(true);
   }
 };
 
 // Forward declarations
 template <typename T>
 class MatrixProxy;
 
 template <typename T>
 class VectorProxy : public Proxy<T, 1> {
 public:
   /// Wrap on memory
   VectorProxy(T* mem, UInt n);
-  VectorProxy(VectorProxy& o) : Proxy<T, 1>(o) {}
-  virtual ~VectorProxy() {}
+  using Proxy<T, 1>::Proxy;
+  ~VectorProxy() override = default;
 
 public:
   /// = operator
   using Proxy<T, 1>::operator=;
   /// Matrix-vector multiplication
   void mul(const MatrixProxy<T>& mat, const VectorProxy<T>& vec);
   /// Dot product
   T dot(const VectorProxy<T>& vec);
   /// L2 norm
   T l2norm() const;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 class MatrixProxy : public Proxy<T, 2> {
 public:
   /// Wrap on memory
   MatrixProxy(T* mem, UInt n, UInt m);
-  MatrixProxy(MatrixProxy& o) : Proxy<T, 2>(o) {}
-  virtual ~MatrixProxy() {}
+  using Proxy<T, 2>::Proxy;
+  ~MatrixProxy() override = default;
 
 public:
   /// = operator
   using Proxy<T, 2>::operator=;
 };
 
 __END_TAMAAS__
 
 #endif  // __TYPES_HH__
diff --git a/tests/test_rough_surface.cpp b/tests/test_rough_surface.cpp
index 3dd884e..7976846 100644
--- a/tests/test_rough_surface.cpp
+++ b/tests/test_rough_surface.cpp
@@ -1,90 +1,91 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "isopowerlaw.hh"
 #include "surface_generator_filter.hh"
 #include "polonsky_keer_rey.hh"
 #include "model_factory.hh"
 /* -------------------------------------------------------------------------- */
 
 using namespace tamaas;
 
 int main(int argc, char * argv[]) {
   /// Surface parameters
   UInt grid_size = 512;
   const UInt k0 = 4, k1 = 4, k2 = 32;
   const Real hurst = 0.8;
 
   /// Initialize Tamaas
   initialize();
 
   /// Surface generation (2D)
   SurfaceGeneratorFilter<2> sg;
   sg.setSizes({{grid_size, grid_size}});
   sg.setRandomSeed(0);
 
   Isopowerlaw<2> iso;
   iso.setQ0(k0);
   iso.setQ1(k1);
   iso.setQ2(k2);
   iso.setHurst(hurst);
 
   sg.setFilter(&iso);
 
   std::cout << "Building rough surface with properties\n"
 	    << "    Q0 = " << k0 << '\n'
 	    << "    Q1 = " << k1 << '\n'
 	    << "    Q2 = " << k2 << '\n'
 	    << "    Hurst = " << hurst << '\n'
 	    << "    size = " << grid_size << 'x' << grid_size << std::endl;
 
   /// Making surface
   GridBase<Real> & grid = sg.buildSurface();
 
   std::cout << "Creating model" << std::endl;
 
   /// Model initialization
-  Model* model = ModelFactory::createModel(model_type::basic_2d, {1., 1.},
-                                          {grid_size, grid_size});
+  std::unique_ptr<Model> model{ModelFactory::createModel(
+      model_type::basic_2d, {1., 1.}, {grid_size, grid_size})};
+
   Real load = 0.1;
   Real precision = 1e-12;
 
   std::cout << "Solving contact problem" << std::endl;
 
   /// Solve normal contact problem
   PolonskyKeerRey solver{*model, grid, precision, PolonskyKeerRey::pressure,
                          PolonskyKeerRey::pressure};
 
   solver.solve(load);
 
   std::cout << "Finished" << std::endl;
 
   /// Cleanup Tamaas
   finalize();
 
   return EXIT_SUCCESS;
 }