diff --git a/python/cast.hh b/python/cast.hh
index 1c39758..5001747 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,211 +1,211 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __CAST_HH__
 #define __CAST_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid_base.hh"
 #include "numpy.hh"
 #include "surface.hh"
 /* -------------------------------------------------------------------------- */
 #include <boost/preprocessor/seq.hpp>
 #include <pybind11/cast.h>
 /* -------------------------------------------------------------------------- */
 
 namespace pybind11 {
 
 // Format descriptor necessary for correct wrap of tamaas complex type
 template <typename T>
 struct format_descriptor<
     tamaas::complex<T>, detail::enable_if_t<std::is_floating_point<T>::value>> {
   static constexpr const char c = format_descriptor<T>::c;
   static constexpr const char value[3] = {'Z', c, '\0'};
   static std::string format() { return std::string(value); }
 };
 
 #ifndef PYBIND11_CPP17
 
 template <typename T>
 constexpr const char format_descriptor<
     tamaas::complex<T>,
     detail::enable_if_t<std::is_floating_point<T>::value>>::value[3];
 
 #endif
 
 namespace detail {
 
 // declare tamaas complex as a complex type for pybind11
 template <typename T>
 struct is_complex<tamaas::complex<T>> : std::true_type {};
 
 template <typename T>
 struct is_fmt_numeric<tamaas::complex<T>,
                       detail::enable_if_t<std::is_floating_point<T>::value>>
     : std::true_type {
   static constexpr int index = is_fmt_numeric<T>::index + 3;
 };
 
 static inline handle policy_switch(return_value_policy policy, handle parent) {
   switch (policy) {
   case return_value_policy::copy:
   case return_value_policy::move:
     return handle();
   case return_value_policy::automatic_reference:  // happens in python-derived
                                                   // classes
   case return_value_policy::reference:
     return none();
   case return_value_policy::reference_internal:
     return parent;
   default:
     TAMAAS_EXCEPTION("Policy is not handled");
   }
 }
 
 template <class array, typename T, tamaas::UInt dim>
 handle grid_to_python(const tamaas::Grid<T, dim>& grid,
                       return_value_policy policy, handle parent) {
   parent = policy_switch(policy, parent);  // reusing variable
   std::vector<tamaas::UInt> sizes(dim);
   std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
 
   if (grid.getNbComponents() != 1)
     sizes.push_back(grid.getNbComponents());
 
   return array(sizes, grid.getInternalData(), parent).release();
 }
 
 /**
  * Type caster for grid classes
  * inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
  * and pybind11/eigen.h
  */
 template <template <typename, tamaas::UInt> class G, typename T,
           tamaas::UInt dim>
 struct type_caster<G<T, dim>> {
   using type = G<T, dim>;
   using array_type =
       array_t<typename type::value_type, 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 (!array_type::check_(src) or !convert)
       return false;
 
     auto buf = array_type::ensure(src);
     if (!buf) return false;
     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) {
     return grid_to_python<array_type>(src, policy, parent);
   }
 };
 
 /**
  * Type caster for GridBase classes
  */
 template <typename T>
 struct type_caster<tamaas::GridBase<T>> {
   using type = tamaas::GridBase<T>;
   using array_type =
       array_t<typename type::value_type, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("GridBaseWrap<T>"));
 
   bool load(handle src, bool convert) {
     if (!array_type::check_(src) or !convert)
       return false;
 
     auto buf = array_type::ensure(src);
     if (!buf) return false;
     value.move(tamaas::wrap::GridBaseNumpy<T>(buf));
 
     return true;
   }
 
   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 = dynamic_cast<decltype(conv)>(src);      \
+    const auto& conv = dynamic_cast<const tamaas::Grid<T, dim>&>(src);         \
     return grid_to_python<array_type>(conv, policy, parent);                   \
   }
 
     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>;
   using array_type =
       array_t<typename type::value_type, array::c_style | array::forcecast>;
 
 public:
   PYBIND11_TYPE_CASTER(type, _("SurfaceWrap<T, dim>"));
 
   bool load(handle src, bool convert) {
     if (!array_type::check_(src) or !convert)
       return false;
 
     auto buf = array_type::ensure(src);
     if (!buf) return false;
     value.move(tamaas::wrap::SurfaceNumpy<type>(buf));
 
     return true;
   }
 
   static handle cast(const type& src, return_value_policy policy,
                      handle parent) {
     return grid_to_python<array_type>(src, policy, parent);
   }
 };
 }  // namespace detail
 }  // namespace pybind11
 
 #endif  // __CAST_HH__
diff --git a/src/bem/bem_grid_polonski.hh b/src/bem/bem_grid_polonski.hh
index 7bfa662..c01163e 100644
--- a/src/bem/bem_grid_polonski.hh
+++ b/src/bem/bem_grid_polonski.hh
@@ -1,77 +1,77 @@
 /**
  *
  * @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_GRID_POLONSKI_HH__
 #define __BEM_GRID_POLONSKI_HH__
 /* -------------------------------------------------------------------------- */
 #include "bem_grid_kato.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 // typedef for swig
-typedef Grid<Real, 1> Grid1dReal;
+using Grid1dReal = Grid<Real, 1>;
 
 /**
  * @brief multi component implementation of Polonski algorithm
  * This method's main unknown is the pressure.
  */
 class BemGridPolonski : public BemGridKato {
 public:
   /// Constructor
   BemGridPolonski(Surface<Real>& surface);
   /// Destructor
   virtual ~BemGridPolonski();
 
 public:
   /* --------------------------------------------------------------------------
    */
   /* Actual computations */
   /* --------------------------------------------------------------------------
    */
   /// Compute equilibrium
   virtual void computeEquilibrium(Real tol, const Grid<Real, 1>& p_target);
   /// Compute gap functional
   Real computeGradientNorm();
   /// Compute optimal step size
   Real computeOptimalStep();
   /// Center field on contact area
   void centerOnContactArea(Grid<Real, 2>& field);
   /// Update search direction
   void updateSearchDirection(Real delta, Real G, Real G_old);
   /// Update tractions
   Real updateTractions(Real tau);
   /// Set pressure on all points
   void setPressure(const Grid<Real, 1>& p_target);
 
 protected:
   Grid<Real, 2> search_direction;
   Grid<Real, 2> projected_search_direction;
 };
 
 }  // namespace tamaas
 
 #endif  // __BEM_GRID_POLONSKI_HH__
diff --git a/src/core/fftransform.cpp b/src/core/fftransform.cpp
index 43bc827..9717c20 100644
--- a/src/core/fftransform.cpp
+++ b/src/core/fftransform.cpp
@@ -1,55 +1,55 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #include "fftransform.hh"
 
 namespace tamaas {
 
 template <typename T, UInt dim>
 FFTransform<T, dim>::FFTransform(Grid<T, dim>& real_,
                                  GridHermitian<T, dim>& spectral_)
     : real(), spectral() {
   this->real.wrap(real_);
   this->spectral.wrap(spectral_);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
-FFTransform<T, dim>::~FFTransform() {}
+FFTransform<T, dim>::~FFTransform() = default;
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void FFTransform<T, dim>::normalize() {
   real *= 1. / real.getNbPoints();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template class FFTransform<Real, 1>;
 template class FFTransform<Real, 2>;
 template class FFTransform<Real, 3>;
 
 }  // namespace tamaas
diff --git a/src/core/fftransform_fftw.cpp b/src/core/fftransform_fftw.cpp
index 2e8c58c..c69a16b 100644
--- a/src/core/fftransform_fftw.cpp
+++ b/src/core/fftransform_fftw.cpp
@@ -1,103 +1,103 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #include "fftransform_fftw.hh"
 
 #define TAMAAS_FFTW_FLAG FFTW_ESTIMATE
 
 namespace tamaas {
 
 template <typename T, UInt dim>
 FFTransformFFTW<T, dim>::FFTransformFFTW(Grid<T, dim>& real,
                                          GridHermitian<T, dim>& spectral)
     : FFTransform<T, dim>(real, spectral), forward_plan(nullptr),
       backward_plan(nullptr) {
   TAMAAS_ASSERT(real.getNbComponents() == spectral.getNbComponents(),
                 "components for FFTW do not match");
 
   // Data pointers
   Real* in = this->real.getInternalData();
-  complex_t* out =
+  auto* out =
       reinterpret_cast<complex_t*>(this->spectral.getInternalData());
   int components = real.getNbComponents();
 
   // fftw parameters
   int rank = dim;  // dimension of fft
   int n[dim] = {0};
   std::copy(real.sizes().begin(), real.sizes().end(),
             n);              // size of individual fft
   int howmany = components;  // how many fft to compute
   int idist = 1, odist = 1;  // components are next to each other in memory
   int istride = real.getStrides().back() * components,
       ostride = spectral.getStrides().back() * components;
   int *inembed = nullptr, *onembed = nullptr;  // row major
 
 #if TAMAAS_FFTW_FLAG != FFTW_ESTIMATE
   // backup the input
   Grid<Real, dim> backup(real);
 #endif
 
   forward_plan =
       fftw::plan_many_forward(rank, n, howmany, in, inembed, istride, idist,
                               out, onembed, ostride, odist, TAMAAS_FFTW_FLAG);
   backward_plan =
       fftw::plan_many_backward(rank, n, howmany, out, onembed, ostride, odist,
                                in, inembed, istride, idist, TAMAAS_FFTW_FLAG);
 
 #if TAMAAS_FFTW_FLAG != FFTW_ESTIMATE
   // restore the backup
   real = backup;
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 FFTransformFFTW<T, dim>::~FFTransformFFTW() {
   fftw::destroy(forward_plan);
   fftw::destroy(backward_plan);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void FFTransformFFTW<T, dim>::forwardTransform() {
   fftw::execute(forward_plan);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void FFTransformFFTW<T, dim>::backwardTransform() {
   fftw::execute(backward_plan);
   this->normalize();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template class FFTransformFFTW<Real, 1>;
 template class FFTransformFFTW<Real, 2>;
 template class FFTransformFFTW<Real, 3>;
 
 }  // namespace tamaas
diff --git a/src/core/fftransform_fftw.hh b/src/core/fftransform_fftw.hh
index 6313469..9755b6f 100644
--- a/src/core/fftransform_fftw.hh
+++ b/src/core/fftransform_fftw.hh
@@ -1,62 +1,62 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 
 #ifndef FFTRANSFORM_FFTW_HH
 #define FFTRANSFORM_FFTW_HH
 
 #include "fftransform.hh"
 #include "fftw_interface.hh"
 
 namespace tamaas {
 
 /**
  * @brief Class wrapping fftw's interface
  */
 template <typename T, UInt dim>
 class FFTransformFFTW : public FFTransform<T, dim> {
   using plan_t = typename fftw::helper<T>::plan;
   using complex_t = typename fftw::helper<T>::complex;
 
 public:
   /// Constructor
   FFTransformFFTW(Grid<T, dim>& real, GridHermitian<T, dim>& spectral);
 
   /// Destructor
-  virtual ~FFTransformFFTW();
+  ~FFTransformFFTW() override;
 
 public:
   /// Perform FFT
   void forwardTransform() override;
 
   /// Perform IFFT
   void backwardTransform() override;
 
 private:
   plan_t forward_plan;
   plan_t backward_plan;
 };
 
 }  // namespace tamaas
 
 #endif  // FFTRANSFORM_FFTW_HH
diff --git a/src/core/loop.hh b/src/core/loop.hh
index 33354c7..df91c9b 100644
--- a/src/core/loop.hh
+++ b/src/core/loop.hh
@@ -1,234 +1,235 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef LOOP_HH
 #define LOOP_HH
 /* -------------------------------------------------------------------------- */
 #include "loops/apply.hh"
 #include "loops/loop_utils.hh"
 #include "ranges.hh"
 #include "tamaas.hh"
 #include <thrust/execution_policy.h>
 #include <thrust/for_each.h>
 #include <thrust/iterator/counting_iterator.h>
 #include <thrust/iterator/zip_iterator.h>
 #include <thrust/transform_reduce.h>
 #include <thrust/tuple.h>
 #include <thrust/version.h>
 #include <type_traits>
 #include <utility>
 
 namespace tamaas {
 
 template <typename T>
 struct is_policy : std::false_type {};
 template <>
 struct is_policy<thrust::detail::host_t> : std::true_type {};
 template <>
 struct is_policy<thrust::detail::device_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::host_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::device_t> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::host_t&> : std::true_type {};
 template <>
 struct is_policy<const thrust::detail::device_t&> : std::true_type {};
 
 /**
  * @brief Singleton class for automated loops using lambdas
 
  * This class is sweet candy :) It provides abstraction of the paralelism
  * paradigm used in loops and allows simple and less error-prone loop syntax,
  * with minimum boiler plate. I love it <3
  */
 class Loop {
 public:
   /// Backends enumeration
   enum backend {
     omp,   ///< [OpenMP](http://www.openmp.org/specifications/) backend
     cuda,  ///< [Cuda](http://docs.nvidia.com/cuda/index.html) backend
   };
 
   /// Helper class to count iterations within lambda-loop
   template <typename T>
   class arange {
   public:
     using it_type = thrust::counting_iterator<T>;
     using reference = typename it_type::reference;
     arange(T start, T size) : start(start), range_size(size) {}
     it_type begin(UInt /*i*/ = 1) const { return it_type(T(start)); }
     it_type end(UInt /*i*/ = 1) const { return it_type(range_size); }
     UInt getNbComponents() const { return 1; }
 
   private:
     T start, range_size;
   };
 
   template <typename T>
   static arange<T> range(T size) {
     return arange<T>(0, size);
   }
 
   template <typename T, typename U>
   static arange<T> range(U start, T size) {
     return arange<T>(start, size);
   }
 
 private:
   /// Replacement for thrust::transform_iterator which copies values
   template <typename Iterator, typename Functor, typename Value>
   class transform_iterator
       : public thrust::iterator_adaptor<
             transform_iterator<Iterator, Functor, Value>, Iterator, Value,
             thrust::use_default, thrust::use_default, Value> {
     Functor func;
 
   public:
     using super_t =
         thrust::iterator_adaptor<transform_iterator<Iterator, Functor, Value>,
                                  Iterator, Value, thrust::use_default,
                                  thrust::use_default, Value>;
     __host__ __device__ transform_iterator(const Iterator& x,
                                            const Functor& func)
         : super_t(x), func(func) {}
 
     friend class thrust::iterator_core_access;
 
   private:
     __host__ __device__ auto dereference() const
         -> decltype(func(*this->base())) {
       return func(*this->base());
     }
   };
 
 public:
   /// Loop functor over ranges
   template <typename Functor, typename... Ranges>
   static auto loop(Functor&& func, Ranges&&... ranges) ->
       typename std::enable_if<not is_policy<Functor>::value, void>::type {
     loopImpl(thrust::device, std::forward<Functor>(func),
              std::forward<Ranges>(ranges)...);
   }
 
   /// Loop over ranges with non-default policy
   template <typename DerivedPolicy, typename Functor, typename... Ranges>
   static void loop(const thrust::execution_policy<DerivedPolicy>& policy,
                    Functor&& func, Ranges&&... ranges) {
     loopImpl(policy, std::forward<Functor>(func),
              std::forward<Ranges>(ranges)...);
   }
 
 private:
   template <typename T>
   using reference_type = typename std::decay<T>::type::reference;
 
 public:
   /// Reduce functor over ranges
   template <operation op, typename Functor, typename... Ranges>
   static auto reduce(Functor&& func, Ranges&&... ranges) ->
       typename std::enable_if<
           not is_policy<Functor>::value,
           decltype(func(std::declval<reference_type<Ranges>>()...))>::type {
     return reduceImpl<op>(thrust::device, std::forward<Functor>(func),
                           std::forward<Ranges>(ranges)...);
   }
 
   /// Reduce over ranges with non-default policy
   template <operation op, typename DerivedPolicy, typename Functor,
             typename... Ranges>
   static auto reduce(const thrust::execution_policy<DerivedPolicy>& policy,
                      Functor&& func, Ranges&&... ranges)
       -> decltype(func(std::declval<reference_type<Ranges>>()...)) {
     return reduceImpl<op>(policy, std::forward<Functor>(func),
                           std::forward<Ranges>(ranges)...);
   }
 
 private:
   /// Loop over ranges and apply functor
   template <typename DerivedPolicy, typename Functor, typename... Ranges>
   static void loopImpl(const thrust::execution_policy<DerivedPolicy>& policy,
                        Functor&& func, Ranges&&... ranges);
   /// Loop over ranges, apply functor and reduce result
   template <operation op, typename DerivedPolicy, typename Functor,
             typename... Ranges>
   static auto reduceImpl(const thrust::execution_policy<DerivedPolicy>& policy,
                          Functor&& func, Ranges&&... ranges)
       -> decltype(func(std::declval<reference_type<Ranges>>()...));
 
+public:
   /// Constructor
   Loop() = delete;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation                                                    */
 /* -------------------------------------------------------------------------- */
 template <typename DerivedPolicy, typename Functor, typename... Ranges>
 void Loop::loopImpl(const thrust::execution_policy<DerivedPolicy>& policy,
                     Functor&& func, Ranges&&... ranges) {
   auto begin = thrust::make_zip_iterator(thrust::make_tuple(ranges.begin()...));
   auto end = thrust::make_zip_iterator(thrust::make_tuple(ranges.end()...));
 
   checkLoopSize(ranges...);
 
   thrust::for_each(policy, begin, end,
                    detail::ApplyFunctor<Functor>(std::forward<Functor>(func)));
 
 #if (defined(USE_CUDA) && THRUST_VERSION < 100904)
   cudaDeviceSynchronize();
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 template <operation op, typename DerivedPolicy, typename Functor,
           typename... Ranges>
 auto Loop::reduceImpl(const thrust::execution_policy<DerivedPolicy>& policy,
                       Functor&& func, Ranges&&... ranges)
     -> decltype(func(std::declval<reference_type<Ranges>>()...)) {
   using return_type = decltype(func(std::declval<reference_type<Ranges>>()...));
   using apply_type = detail::ApplyFunctor<Functor, return_type>;
 
   checkLoopSize(ranges...);
 
   auto applier = apply_type(func);
   detail::reduction_helper<op, return_type> red_helper;
 
   auto begin_zip =
       thrust::make_zip_iterator(thrust::make_tuple(ranges.begin()...));
   auto end_zip = thrust::make_zip_iterator(thrust::make_tuple(ranges.end()...));
 
   transform_iterator<decltype(begin_zip), apply_type, return_type> begin(
       begin_zip, applier);
   transform_iterator<decltype(end_zip), apply_type, return_type> end(end_zip,
                                                                      applier);
 
   auto result =
       thrust::reduce(policy, begin, end, red_helper.init(), red_helper);
 #if (defined(USE_CUDA) && THRUST_VERSION < 100904)
   cudaDeviceSynchronize();
 #endif
   return result;
 }
 
 }  // namespace tamaas
 
 #endif  // LOOP_HH
diff --git a/src/core/ranges.hh b/src/core/ranges.hh
index c5817b4..1100e0e 100644
--- a/src/core/ranges.hh
+++ b/src/core/ranges.hh
@@ -1,114 +1,115 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef RANGES_HH
 #define RANGES_HH
 /* -------------------------------------------------------------------------- */
 #include "iterator.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /// Range class for complex iterators
 template <typename LocalType, typename ValueType, UInt local_size>
 class Range {
 
   class iterator : public iterator_::iterator<ValueType> {
     using parent = iterator_::iterator<ValueType>;
 
   public:
     using value_type = LocalType;
     using reference = value_type;
     iterator(const parent& o) : parent(o) {}
 
     inline reference operator*() { return LocalType(this->data); }
     inline reference operator*() const { return LocalType(this->data); }
   };
 
 public:
   // useful type trait
   template <typename Container>
   struct is_valid_container
       : std::is_same<
             std::remove_cv_t<typename std::decay_t<Container>::value_type>,
             std::remove_cv_t<ValueType>> {};
 
 public:
   using value_type = LocalType;
   using reference = value_type&&;
 
   /// Construct from a container
   template <class Container>
   Range(Container&& cont) : _begin(cont.begin()), _end(cont.end()) {
     if (cont.getNbComponents() != local_size)
       TAMAAS_EXCEPTION(
           "Number of components does not match local tensor type size ("
           << cont.getNbComponents() << ", expected " << local_size << ")");
     _begin.setStep(local_size);
     _end.setStep(local_size);
   }
 
   /// Construct from two iterators
-  Range(iterator _begin, iterator _end) : _begin(_begin), _end(_end) {}
+  Range(iterator _begin, iterator _end)
+      : _begin(std::move(_begin)), _end(std::move(_end)) {}
 
   iterator begin() { return _begin; }
   iterator end() { return _end; }
 
   Range headless() const {
     iterator shift(_begin);
     return Range{++shift, _end};
   }
 
 private:
   iterator _begin, _end;
 };
 
 /// Make range with proxy type
 template <typename LocalType, class Container>
 std::enable_if_t<
     is_proxy<LocalType>::value,
     Range<LocalType, typename LocalType::value_type, LocalType::size>>
 range(Container&& cont) {
   using RangeClass =
       Range<LocalType, typename LocalType::value_type, LocalType::size>;
   static_assert(
       RangeClass::template is_valid_container<Container>::value,
       "Mismatch between range's value type and container's value type");
   return RangeClass(std::forward<Container>(cont));
 }
 
 /// Make range with standard type
 template <typename LocalType, class Container>
 std::enable_if_t<not is_proxy<LocalType>::value, Range<LocalType, LocalType, 1>>
 range(Container&& cont) {
   using RangeClass = Range<LocalType, LocalType, 1>;
   static_assert(
       RangeClass::template is_valid_container<Container>::value,
       "Mismatch between range's value type and container's value type");
   return RangeClass(std::forward<Container>(cont));
 }
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 
 #endif
diff --git a/src/core/surface.cpp b/src/core/surface.cpp
index cb3158b..4d778d7 100644
--- a/src/core/surface.cpp
+++ b/src/core/surface.cpp
@@ -1,422 +1,421 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #include <fstream>
 #include <sstream>
 #ifdef USE_CUDA
 #include <thrust/random.h>
 #else
 #include <random>
 #endif
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef USING_IOHELPER
 #include "dumper_paraview.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 std::string fname = "surface_height.pdf";
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::expandByLinearInterpolation() {
   UInt n = this->size();
-  Surface* old_s = new Surface(n, 1.);
+  auto old_s = std::make_unique<Surface>(n, 1.);
   int new_n = n * 2;
   *old_s = *this;
   this->setGridSize(new_n);
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       (*this)(i, j) = (*old_s)(i, j);
     }
   }
-  delete old_s;
 }
 
 Real sinc(Real x) { return sin(M_PI * x) / (M_PI * x); }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByShannonInterpolation(int /*factor*/) {
 
   /* /!\ this->FFTTransform(); /!\ */
   TAMAAS_EXCEPTION("Shannon interpolation needs old FFTTransform interface");
 #if 0
   UInt n = this->size();
   UInt N = std::pow(2,factor)*this->size();
   std::cout << N << " " << factor << std::endl;
   Surface s = Surface(N,1.);
 
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = 0; j < n/2; ++j) {
       s(i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = 0; j < n/2; ++j) {
       s(N-n+i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = n/2; j < n; ++j) {
       s(N-n+i,N-n+j) = (*this)(i,j);
     }
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = n/2; j < n; ++j) {
       s(i,N-n+j) = (*this)(i,j);
     }
 
   //  s.dumpToTextFile<false,true>("temp.csv");
 
   /* /!\ s.FFTITransform(); /!\ */
 
   this->setGridSize(N);
   Real f = std::pow(2,factor);
   f = f*f;
 
   for (UInt i = 0; i < N*N; ++i)
     (*this)(i) = f*s(i);
 
 
 //   int new_n = 2*n;
 //   Surface s(*this);
 //   this->setGridSize(new_n);
 
 //   // put the original signal with zeros for to be interpolated points
 //   for (int i=0; i<n; i++)
 //     for (int j=0; j<n; j++)
 //       (*this)(2*i,j*2) = s(i,j);
 
 //   // do shannon interpolation
 //   int rcut = 64;
 
 //   //interpolate lines
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (i%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int i2 = i + k;
 // 	if (i2 < 0) i2 += new_n;
 // 	if (i2 >= new_n) i2 -= new_n;
 // 	if (i2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += s(i2,j)*sinc(d/2);
 //       }
 //     }
 //   }
 
 //   //interpolate columns
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (j%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int j2 = j + k;
 // 	if (j2 < 0) j2 += new_n;
 // 	if (j2 >= new_n) j2 -= new_n;
 // 	if (j2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2),
 // 				s[i+j2*new_n].imag()*sinc(d/2));
 //       }
 //     }
 //   }
 
 
 // //   for (int i=0; i<n; i++) {
 // //     for (int j=0; j<n; j++) {
 // //       s[2*i+j*2*new_n].real() = 0;
 // //       s[2*i+j*2*new_n].im = 0;
 // //     }
 // //   }
 #endif
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByMirroring() {
   UInt n = this->size();
   UInt new_n = n * 2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i = 0; i < new_n; i++) {
     for (UInt j = 0; j < new_n; j++) {
       UInt i1 = i;
       if (i1 > n)
         i1 = new_n - i1;
       UInt j1 = j;
       if (j1 > n)
         j1 = new_n - j1;
       if (i == n || j == n) {
         (*this)(i, j) = 0.;
         continue;
       }
       (*this)(i, j) = s(i1, j1);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByPuttingZeros() {
   UInt n = this->size();
   UInt new_n = n * 2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       (*this)(2 * i, j * 2) = s(i, j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::contractByTakingSubRegion() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n / 2);
 
   for (UInt i = 0; i < n / 2; i++) {
     for (UInt j = 0; j < n / 2; j++) {
       (*this)(i, j) = s(i, j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::contractByIgnoringMidPoints() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n / 2);
 
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       if (i % 2)
         continue;
       if (j % 2)
         continue;
       (*this)(i / 2, j / 2) = s(i, j);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 // void Surface<T>::dumpToTextFile(string filename, bool pr) {
 
 //   cout << "Writing surface file " << filename << endl;
 //   ofstream file(filename.c_str());
 //   file.precision(15);
 
 //   if (pr)
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         if (heights[i+j*n].real() < 1e-10)
 //           file << scientific << -1e10 << " " << 0;
 //         else
 //           file << scientific << heights[i+j*n].real() << " " <<
 //           heights[i+j*n].im;
 //         file << endl;
 //       }
 //       file << endl;
 //     }
 //   }
 //   else
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         file << scientific << heights[i+j*n].real() << " " <<
 //         heights[i+j*n].im; file << endl;
 //       }
 //     }
 //   }
 //   file.close();
 
 // }
 
 /* -------------------------------------------------------------------------- */
 
 // Surface * Surface<T>::loadFromTextFile(std::string filename) {
 
 // std::string line;
 // std::ifstream myfile (filename.c_str());
 
 // if (!myfile.is_open()){
 //   SURFACE_FATAL("cannot open file " << filename);
 // }
 
 // Surface * my_surface = nullptr;
 // UInt n = 0;
 // if (myfile.is_open())
 //   {
 //     while (! myfile.eof() )
 // 	{
 // 	  getline (myfile,line);
 // 	  if (line[0] == '#')
 // 	    continue;
 // 	  if (line.size() == 0)
 // 	    continue;
 
 // 	  std::stringstream s;
 // 	  s << line;
 // 	  int i,j;
 // 	  s >> i;
 // 	  s >> j;
 // 	  if (n < i) n = i;
 // 	  if (n < j) n = j;
 // 	}
 //     ++n;
 //     std::cout << "read " << n << " x " << n << " grid points" << std::endl;
 
 //     myfile.clear();              // forget we hit the end of file
 //     myfile.seekg(0, std::ios::beg);   // move to the start of the file
 
 //     my_surface = new Surface(n);
 
 //     while (! myfile.eof() )
 // 	{
 // 	  getline (myfile,line);
 // 	  if (line[0] == '#')
 // 	    continue;
 // 	  if (line.size() == 0)
 // 	    continue;
 
 // 	  std::stringstream s;
 // 	  s << line;
 // 	  int i,j;
 // 	  Real h,hIm;
 // 	  s >> i;
 // 	  s >> j;
 // 	  s >> h;
 // 	  s >> hIm;
 // 	  my_surface->heights[i+j*n] = complex(h,hIm);
 // 	}
 //     myfile.close();
 //   }
 //   return my_surface;
 // }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::dumpToParaview(std::string /*filename*/) const {
 #ifdef USING_IOHELPER
 
   UInt n = this->size();
   const Surface& data = *this;
   /* build 3D coordinates array to be usable with iohelper */
   //  cerr << "allocate coords for paraview" << endl;
   Real* coords = new Real[n * n * 3];
   // Real *  coords = new Real[(2*n-1)*(2*n-1)*3];
   // cerr << "allocate zcoords for paraview" << endl;
   Real* zcoords = new Real[n * n * 3];
   // Real *  zcoords = new Real[(2*n-1)*(2*n-1)*3];
   // cerr << "allocate zcoords2 for paraview" << endl;
   Real* zcoords2 = new Real[n * n * 3];
   //  Real *  zcoords2 = new Real[(2*n-1)*(2*n-1)*3];
   int index = 0;
   // for (int i=-n+1; i<n; i++) {
   //   for (int j=-n+1; j<n; j++) {
   for (int i = 0; i < n; i++) {
     for (int j = 0; j < n; j++) {
       coords[3 * index] = 1. * i / n;
       coords[3 * index + 1] = 1. * j / n;
       coords[3 * index + 2] = 0;
       zcoords[3 * index] = 0;
       zcoords[3 * index + 1] = 0;
       zcoords2[3 * index] = 0;
       zcoords2[3 * index + 1] = 0;
       int i1 = (i % n);
       int j1 = (j % n);
       if (i1 < 0)
         i1 += n;
       if (j1 < 0)
         j1 += n;
       zcoords[3 * index + 2] = data(i1, j1).real();
       zcoords2[3 * index + 2] = data(i1, j1).imag();
       ++index;
     }
   }
 
   iohelper::DumperParaview dumper;
   dumper.setPoints(coords, 3, index, filename.c_str());
   dumper.setConnectivity(nullptr, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE);
   dumper.addNodeDataField("data.re", zcoords, 3, index);
   dumper.addNodeDataField("data.im", zcoords2, 3, index);
   dumper.setPrefix("./paraview");
   dumper.init();
   dumper.dump();
 
   delete[] coords;
   delete[] zcoords;
   delete[] zcoords2;
 #endif  // USING_IOHELPER
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::recenterHeights() {
   UInt n = this->size();
   Real zmean = SurfaceStatistics::computeAverage(*this);
 
   for (UInt i = 0; i < n * n; ++i)
     (*this)(i) = (*this)(i)-zmean;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::generateWhiteNoiseSurface(Surface<Real>& sources,
                                            long random_seed) {
   // set the random gaussian generator
   random_engine gen(random_seed);
   normal_distribution<Real> gaussian_generator;
 
   UInt n = sources.size();
 
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       sources(i, j) = gaussian_generator(gen);
     }
 }
 /* -------------------------------------------------------------------------- */
 
 template class Surface<Real>;
 template class Surface<UInt>;
 template class Surface<Int>;
 template class Surface<bool>;
 
 }  // namespace tamaas
diff --git a/src/core/surface.hh b/src/core/surface.hh
index c95c163..3affae4 100644
--- a/src/core/surface.hh
+++ b/src/core/surface.hh
@@ -1,117 +1,117 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #ifndef __SURFACE_HH__
 #define __SURFACE_HH__
 /* -------------------------------------------------------------------------- */
 #include <complex>
 #include <iostream>
 #include <limits>
 #include <string>
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @deprecated should use Grid class instead
  * @brief Square 2D surface class
  */
 template <typename T>
 class Surface : public Grid<T, 2> {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 
 public:
   Surface() : Grid<T, 2>() {}
 
   Surface(UInt a, Real L) : Grid<T, 2>() {
     setGridSize(a);
     setL(L);
   }
 
-  virtual ~Surface() {}
+  ~Surface() override = default;
 
   Surface(const Surface& o) : Grid<T, 2>(o) {}
   Surface(Surface&& o) : Grid<T, 2>(std::forward<Surface>(o)) {}
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
   //! set the grid size
   void setGridSize(UInt n) { this->resize({n, n}); }
   //! get lateral number of discretization points
   inline UInt size() const { return this->n[0]; };
   //! get the lateral size of the sample
   inline Real getL() const { return 1.; };
   //! get the lateral size of the sample
   inline void setL(Real) {}
 
   Surface& operator=(T val) {
     Grid<T, 2>::operator=(val);
     return *this;
   }
   Surface& operator=(const Surface& o) {
     Grid<T, 2>::operator=(o);
     return *this;
   }
   Surface& operator=(Surface&& o) {
     Grid<T, 2>::operator=(std::forward<Surface>(o));
     return *this;
   }
 
   //! Write VTU file that is a set of points
   void dumpToParaview(std::string filename) const;
   //! create a new surface with linear interpolated points
   void expandByLinearInterpolation();
   //! create a new surface with linear interpolated points
   void expandByShannonInterpolation(int factor);
   //! create a new surface with mirroring
   void expandByMirroring();
   //! create a new surface with mirroring
   void expandByPuttingZeros();
   //! contract the surface by ignoring extra data
   void contractByIgnoringMidPoints();
   //! contract the surface by taking left bottom part of it
   void contractByTakingSubRegion();
   // recenter the heights to zero value
   void recenterHeights();
   // set all the members of the map to a given value
   static void generateWhiteNoiseSurface(Surface<Real>& surface,
                                         long random_seed);
   virtual const Surface<T>& getWrappedNumpy() { return *this; };
 };
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 
 #include "surface_complex.hh"
 
 #endif /* __SURFACE_HH__ */
diff --git a/src/core/surface_statistics.hh b/src/core/surface_statistics.hh
index 73dc73f..26e1a97 100644
--- a/src/core/surface_statistics.hh
+++ b/src/core/surface_statistics.hh
@@ -1,709 +1,702 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef __SURFACE_STATISTICS_HH__
 #define __SURFACE_STATISTICS_HH__
 /* -------------------------------------------------------------------------- */
 #include "fft_plan_manager.hh"
 #include "fftransform.hh"
 #include "surface.hh"
 #include <thrust/complex.h>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @deprecated should move code to Statistics class instead
  * @brief Statistics calculations on surfaces
  */
 class SurfaceStatistics {
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-public:
-  SurfaceStatistics(){};
-  virtual ~SurfaceStatistics(){};
-
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   inline static Real computeMaximum(const Surface<Real>& s);
   inline static Real computeMinimum(const Surface<Real>& s);
 
   inline static Real
   computeAverage(const Surface<Real>& s,
                  Real epsilon = std::numeric_limits<Real>::lowest()) {
     return computeAverage<Real>(s, epsilon);
   }
 
   template <typename T>
   inline static T computeAverage(const Surface<T>& s,
                                  T epsilon = std::numeric_limits<T>::lowest());
   inline static Real
   computeStdev(const Grid<Real, 2>& s,
                Real epsilon = std::numeric_limits<Real>::lowest());
 
   inline static Real computeSkewness(const Surface<Real>& surface, Real avg,
                                      Real std);
   inline static Real computeKurtosis(const Surface<Real>& surface, Real avg,
                                      Real std);
   inline static Real computeRMSSlope(const Surface<Real>& surface);
   inline static Real computeSpectralRMSSlope(Surface<Real>& surface);
   inline static Real computeSpectralMeanCurvature(Surface<Real>& surface);
   inline static Real computeSpectralStdev(Surface<Real>& surface);
   //  inline static void computeStatistics(Surface<Real> & surface);
   inline static std::vector<Real> computeMoments(Surface<Real>& surface);
   inline static Real computeAnalyticFractalMoment(UInt order, UInt k1, UInt k2,
                                                   Real hurst, Real rms, Real L);
   //! compute perimeter
   inline static Real computePerimeter(const Surface<Real>& surface,
                                       Real eps = 1e-10);
   //! get real contact area
   inline static Real computeContactArea(const Surface<Real>& surface);
   //! get real contact area
   inline static Real computeContactAreaRatio(const Surface<Real>& surface);
 
   template <bool consider_zeros = false>
   inline static std::vector<int>
   computeDistribution(const Surface<Real>& surface, unsigned int n_bins,
                       Real min_value = std::numeric_limits<Real>::max(),
                       Real max_value = std::numeric_limits<Real>::lowest());
 
   inline static std::vector<Complex>
   computeSpectralDistribution(const Surface<Real>& surface, unsigned int n_bins,
                               UInt cutoff = std::numeric_limits<UInt>::max());
 
   //! compute average amplitude
   inline static Real computeAverageAmplitude(const Surface<Real>& surface);
   inline static Complex
   computeAverageAmplitude(const SurfaceComplex<Real>& surface);
 
   //! sum all heights
   inline static Real computeSum(const Surface<Real>& surface);
   //! compute auto correlation
   inline static Surface<Real> computeAutocorrelation(Surface<Real>& surface);
   //! compute powerspectrum
   inline static SurfaceComplex<Real>
   computePowerSpectrum(Surface<Real>& surface);
 };
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 //#include "surface_statistics_inline_impl.cc"
 
 /* -------------------------------------------------------------------------- */
 
 template <bool consider_zeros>
 std::vector<int>
 SurfaceStatistics::computeDistribution(const Surface<Real>& surface,
                                        unsigned int n_bins, Real min_value,
                                        Real max_value) {
 
   if (n_bins == 0)
     SURFACE_FATAL("badly set number of bins");
 
   Real min_heights;
   Real max_heights;
 
   std::vector<int> bins;
   bins.resize(n_bins);
   bins.assign(n_bins, 0);
 
   min_heights = min_value;
   max_heights = max_value;
 
   if (max_value == -1. * std::numeric_limits<Real>::max())
     max_heights = computeMaximum(surface);
   if (min_value == std::numeric_limits<Real>::max())
     min_heights = computeMinimum(surface);
 
   int counter = 0;
   const int n = surface.size();
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     if (!consider_zeros && value < zero_threshold)
       continue;
     if (value > max_heights || value < min_heights)
       SURFACE_FATAL("out of bound value: "
                     << min_heights << " <? " << value << " <? " << max_heights
                     << " : bounds probably badly chosen");
     Real ind = (value - min_heights) / (max_heights - min_heights);
     unsigned int index = ind * n_bins;
     if (index == n_bins)
       --index;
     ++bins[index];
     ++counter;
   }
   if (consider_zeros && counter != n * n)
     SURFACE_FATAL("problem in counting points");
 
   //  Real offset = (max_heights - min_heights)/n_bins;
   // Real x = 0;
   //  Real f1 = 0, f2 = 0;
   //  for (unsigned int i = 0 ; i < n_bins ; ++i)
   //  {
   //    x  = i*offset+min_heights;
   // f1 = bins[i]/Real(n*n);
   // f2 = bins[i]/Real(counter);
   //}
   Real min_non_zero = std::numeric_limits<Real>::max();
 #pragma omp parallel for reduction(min : min_non_zero)
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     if (value < min_non_zero && value > 1e-14)
       min_non_zero = value;
   }
   return bins;
 }
 /* -------------------------------------------------------------------------- */
 std::vector<Complex> SurfaceStatistics::computeSpectralDistribution(
     const Surface<Real>& surface, unsigned int n_bins, UInt cutoff) {
 
   if (n_bins == 0)
     SURFACE_FATAL("badly set number of bins");
 
   const int n = surface.size();
 
   Real min_value = SurfaceStatistics::computeMinimum(surface);
   Real max_value = SurfaceStatistics::computeMaximum(surface);
 
   std::vector<Complex> distrib_fft(n_bins);
   const int N = n_bins;
   Real L = max_value - min_value;
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k)
     distrib_fft[k] = 0.;
 
 #pragma omp parallel for
   for (int k = 1; k < N / 2; ++k) {
     for (int i = 0; i < n * n; ++i) {
       Real value = surface(i);
       Real angle = -2. * M_PI * k * (value - min_value) / L;
       Complex phase = std::polar<Real>(1., angle);
       distrib_fft[k] += phase;
       distrib_fft[N - k] += phase;
     }
   }
 
   for (int i = 0; i < n * n; ++i) {
     Real value = surface(i);
     Real angle = -2. * M_PI * N / 2 * (value - min_value) / L;
     Complex phase = std::polar<Real>(1., angle);
     distrib_fft[N / 2] += phase;
   }
 
 #pragma omp parallel for
   for (int k = 1; k < N; ++k)
     distrib_fft[k] /= n * n;
 
   distrib_fft[0] = std::polar(1., 0.);
 
 #pragma omp parallel for
   for (int k = 1; k < N / 2; ++k) {
     if (k >= (int)cutoff) {
       distrib_fft[k] = 0.;
       distrib_fft[N - k] = 0.;
     }
   }
 
     // proceed to inverse fourier transform
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k) {
     distrib_fft[k] = conj(distrib_fft[k]);
   }
 
-  fftw_complex* in = (fftw_complex*)&distrib_fft[0];
+  auto* in = (fftw_complex*)&distrib_fft[0];
   fftw_plan p = fftw_plan_dft_1d(N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   std::vector<Complex> distrib;
   distrib.resize(N);
 
 #pragma omp parallel for
   for (int k = 0; k < N; ++k) {
     distrib[k] = conj(distrib_fft[k]);
   }
 
   return distrib;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAverageAmplitude(const Surface<Real>& surface) {
   return surface.mean();
 }
 /* -------------------------------------------------------------------------- */
 
 Complex SurfaceStatistics::computeAverageAmplitude(
     const SurfaceComplex<Real>& surface) {
   return surface.mean();
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeContactArea(const Surface<Real>& surface) {
   Real L = surface.getL();
   return SurfaceStatistics::computeContactAreaRatio(surface) * L * L;
 }
 
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeContactAreaRatio(const Surface<Real>& surface) {
   UInt area = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : area)
   for (UInt i = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       if (surface(i, j) > 1e-10)
         area++;
     }
   }
   return area / Real(n * n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 Surface<Real>
 SurfaceStatistics::computeAutocorrelation(Surface<Real>& surface) {
   SurfaceComplex<Real> correlation_heights = computePowerSpectrum(surface);
   Surface<Real> acf(correlation_heights.size(), correlation_heights.getL());
   FFTransform<Real, 2>& plan =
       FFTPlanManager::get().createPlan(acf, correlation_heights);
   plan.backwardTransform();
   FFTPlanManager::get().destroyPlan(acf, correlation_heights);
   return acf;
 }
 /* -------------------------------------------------------------------------- */
 
 SurfaceComplex<Real>
 SurfaceStatistics::computePowerSpectrum(Surface<Real>& surface) {
   SurfaceComplex<Real> power_spectrum(surface.size(), surface.getL());
   FFTransform<Real, 2>& plan =
       FFTPlanManager::get().createPlan(surface, power_spectrum);
   plan.forwardTransform();
 
   UInt n = power_spectrum.dataSize();
   UInt factor = power_spectrum.size() * power_spectrum.size();
 
 #pragma omp parallel for
   for (UInt i = 0; i < n; ++i)
     power_spectrum(i) /= factor;
 
   power_spectrum.makeItRealBySquare();
   FFTPlanManager::get().destroyPlan(surface, power_spectrum);
   return power_spectrum;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 T SurfaceStatistics::computeAverage(const Surface<T>& surface, T epsilon) {
 
   Real avg = 0.0;
   UInt cpt = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : avg, cpt)
   for (UInt i = 0; i < n * n; i++) {
     Real val = surface(i);
     if (std::abs(val) > epsilon) {
       avg += val;
       ++cpt;
     }
   }
   avg /= cpt;
   return avg;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSkewness(const Surface<Real>& surface, Real avg,
                                         Real std) {
   Real skw = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : skw)
   for (UInt i = 0; i < n * n; i++) {
     Real tmp = surface(i) - avg;
     skw += (tmp * tmp * tmp);
   }
   skw /= (n * std * std * std);
   return skw;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeKurtosis(const Surface<Real>& surface, Real avg,
                                         Real std) {
 
   Real kur = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : kur)
   for (UInt i = 0; i < n * n; i++) {
     Real tmp = surface.size() - avg;
     kur += (tmp * tmp * tmp * tmp);
   }
   kur /= (n * std * std * std * std);
   return kur - 3;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMaximum(const Surface<Real>& surface) {
   Real _max = -1. * std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(max : _max)
   for (UInt i = 0; i < n * n; i++) {
     _max = std::max(surface(i), _max);
   }
   return _max;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMinimum(const Surface<Real>& surface) {
   Real _min = std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(min : _min)
   for (UInt i = 0; i < n * n; i++) {
     _min = std::min(surface(i), _min);
   }
   return _min;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeRMSSlope(const Surface<Real>& surface) {
 
   UInt n = surface.size();
   Surface<Real> grad(n, 1.);
 
   for (UInt i = 0, m = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       UInt k = i + n * j;
 
       //      xl = ((i==0)?(k+n-1):(k-1));
       UInt xh = ((i == (n - 1)) ? (k - n + 1) : (k + 1));
       //      yl = ((j==0)?(k+n*n-n):(k-n));
       UInt yh = ((j == (n - 1)) ? (k - n * n + n) : (k + n));
 
       Real gx = surface(xh) - surface(k);
       Real gy = surface(yh) - surface(k);
       grad(m++) = gx * gx + gy * gy;
     }
   }
 
   Real avg = SurfaceStatistics::computeAverage(grad);
   Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
   Real scale_factor = surface.getL() / n;
   return slp / scale_factor;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralRMSSlope(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_slopes = 0;
   // /!\ commented code was causing race conditions with dummy variable in
   // GridHermitian
 #pragma omp parallel for collapse(2), reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i)
     for (UInt j = 1; j < n / 2; ++j) {
       rms_slopes += 2. * ((i * i) + (j * j)) * power(i, j).real();
       // rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
       rms_slopes += 2. * ((i * i) + (j * j)) * power(n - i, j).real();
       // rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
     }
 
       // external line (i or j == 0)
 #pragma omp parallel for reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i) {
     rms_slopes += 1. * (i * i) * power(i, 0).real();
     rms_slopes += 1. * (i * i) * power(n - i, 0).real();
     rms_slopes += 2. * (i * i) * power(0, i).real();
     // rms_slopes += 1.*(i*i)*power(0,n-i).real();
   }
 
     // internal line (i or j == n/2)
 #pragma omp parallel for reduction(+ : rms_slopes)
   for (UInt i = 1; i < n / 2; ++i) {
     rms_slopes += 0.5 * (i * i + n * n / 4) * power(i, n / 2).real();
     rms_slopes += 0.5 * (i * i + n * n / 4) * power(n - i, n / 2).real();
     rms_slopes += 1. * (i * i + n * n / 4) * power(n / 2, i).real();
     // rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
   }
 
   {
     // i == n/2 j == n/2
     rms_slopes += 0.25 * (n * n / 2) * power(n / 2, n / 2).real();
     // i == 0 j == n/2
     rms_slopes += 1. * (n * n / 4) * power(0, n / 2).real();
     // i == n/2 j == 0
     rms_slopes += 1. * (n * n / 4) * power(n / 2, 0).real();
   }
   rms_slopes = 2. * M_PI / n * sqrt(rms_slopes);
   Real scale_factor = surface.getL() / n;
   return rms_slopes / scale_factor;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeStdev(const Grid<Real, 2>& surface,
                                      Real epsilon) {
   UInt n = surface.dataSize();
   UInt cpt = 0;
   Real std = 0.0;
   Real avg = 0.0;
 
 #pragma omp parallel for reduction(+ : avg, std, cpt)
   for (UInt i = 0; i < n; i++) {
     Real val = surface(i);
     if (val > epsilon) {
       avg += val;
       std += val * val;
       ++cpt;
     }
   }
   avg /= cpt;
   return sqrt(std / cpt - avg * avg);
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralStdev(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real rms_heights = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : rms_heights)
   for (UInt i = 0; i < n; ++i)
     for (UInt j = 0; j < n; ++j) {
       rms_heights += power(i, j).real();
     }
   rms_heights = sqrt(rms_heights);
   return rms_heights;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralMeanCurvature(Surface<Real>& surface) {
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_curvatures = 0;
 #pragma omp parallel for collapse(2), reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i)
     for (UInt j = 1; j < n / 2; ++j) {
       Real coeff = ((i * i) + (j * j));
       coeff *= coeff;
       rms_curvatures += 1. * coeff * power(i, j).real();
       rms_curvatures += 1. * coeff * power(n - i, n - j).real();
       rms_curvatures += 1. * coeff * power(n - i, j).real();
       rms_curvatures += 1. * coeff * power(i, n - j).real();
     }
 
       // external line (i or j == 0)
 #pragma omp parallel for reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i) {
     Real coeff = i * i;
     coeff *= coeff;
     rms_curvatures += 4. / 3. * coeff * power(i, 0).real();
     rms_curvatures += 4. / 3. * coeff * power(n - i, 0).real();
     rms_curvatures += 4. / 3. * coeff * power(0, i).real();
     rms_curvatures += 4. / 3. * coeff * power(0, n - i).real();
   }
 
     // internal line (i or j == n/2)
 #pragma omp parallel for reduction(+ : rms_curvatures)
   for (UInt i = 1; i < n / 2; ++i) {
     Real coeff = (i * i + n * n / 4);
     coeff *= coeff;
     rms_curvatures += 0.5 * coeff * power(i, n / 2).real();
     rms_curvatures += 0.5 * coeff * power(n - i, n / 2).real();
     rms_curvatures += 0.5 * coeff * power(n / 2, i).real();
     rms_curvatures += 0.5 * coeff * power(n / 2, n - i).real();
   }
 
   {
     // i == n/2 j == n/2
     rms_curvatures +=
         0.25 * (n * n / 2) * (n * n / 2) * power(n / 2, n / 2).real();
     // i == 0 j == n/2
     rms_curvatures +=
         2. / 3. * (n * n / 4) * (n * n / 4) * power(0, n / 2).real();
     // i == n/2 j == 0
     rms_curvatures +=
         2. / 3. * (n * n / 4) * (n * n / 4) * power(n / 2, 0).real();
   }
 
   rms_curvatures *= 9. / 4.;
   rms_curvatures = M_PI / n / n * sqrt(rms_curvatures);
 
   return rms_curvatures;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computePerimeter(const Surface<Real>& surface,
                                          Real eps) {
   int perimeter = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+ : perimeter)
   for (UInt i = 0; i < n; ++i) {
     for (UInt j = 0; j < n; ++j) {
       if (j > 0 && ((surface(i, (j - 1)) < eps && surface(i, j) > eps) ||
                     (surface(i, (j - 1)) > eps && surface(i, j) < eps)))
         perimeter++;
       if (i > 0 && ((surface(j, (i - 1)) < eps && surface(j, i) > eps) ||
                     (surface(j, (i - 1)) > eps && surface(j, i) < eps)))
         perimeter++;
     }
   }
 
   return perimeter / Real(n);
 }
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeSum(const Surface<Real>& surface) {
   Real sum = 0.;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+ : sum)
   for (UInt i = 0; i < n * n; ++i)
     sum += surface(i);
 
   return sum;
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<Real> SurfaceStatistics::computeMoments(Surface<Real>& surface) {
 
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real L = surface.getL();
   UInt n = surface.size();
   Real k_factor = 2. * M_PI / L;
 
   std::vector<std::vector<Real>> moment;
   moment.resize(6);
   for (int i = 0; i < 5; i++)
     moment[i].resize(5);
 
   for (UInt i = 0; i < n / 2; i++) {
     Real i0 = pow(Real(i) * k_factor, 0);
     Real i2 = pow(Real(i) * k_factor, 2);
     Real i4 = pow(Real(i) * k_factor, 4);
     Real im0 = pow(-Real(i) * k_factor, 0);
     Real im2 = pow(-Real(i) * k_factor, 2);
     Real im4 = pow(-Real(i) * k_factor, 4);
 
     for (UInt j = 0; j < n / 2; j++) {
       Real Phi = power(i, j).real();
 
       Real Phi_10;
       Real Phi_01;
       Real Phi_11;
       if (i == 0 && j == 0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = 0;
       } else if (i == 0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = power(i, (n - j)).real();
       } else if (j == 0) {
         Phi_10 = power(n - i, j).real();
         Phi_11 = 0;
         Phi_01 = 0;
       } else {
         Phi_10 = power(n - i, j).real();
         Phi_11 = power(n - i, (n - j)).real();
         Phi_01 = power(i, (n - j)).real();
       }
       Real j0 = pow(Real(j) * k_factor, 0);
       Real j2 = pow(Real(j) * k_factor, 2);
       Real j4 = pow(Real(j) * k_factor, 4);
       Real jm0 = pow(-Real(j) * k_factor, 0);
       Real jm2 = pow(-Real(j) * k_factor, 2);
       Real jm4 = pow(-Real(j) * k_factor, 4);
 
       moment[0][0] += i0 * j0 * Phi + im0 * j0 * Phi_10 + i0 * jm0 * Phi_01 +
                       im0 * jm0 * Phi_11;
       moment[0][2] += i0 * j2 * Phi + im0 * j2 * Phi_10 + i0 * jm2 * Phi_01 +
                       im0 * jm2 * Phi_11;
       moment[2][0] += i2 * j0 * Phi + im2 * j0 * Phi_10 + i2 * jm0 * Phi_01 +
                       im2 * jm0 * Phi_11;
       moment[2][2] += i2 * j2 * Phi + im2 * j2 * Phi_10 + i2 * jm2 * Phi_01 +
                       im2 * jm2 * Phi_11;
       moment[0][4] += i0 * j4 * Phi + im0 * j4 * Phi_10 + i0 * jm4 * Phi_01 +
                       im0 * jm4 * Phi_11;
       moment[4][0] += i4 * j0 * Phi + im4 * j0 * Phi_10 + i4 * jm0 * Phi_01 +
                       im4 * jm0 * Phi_11;
     }
   }
 
   std::vector<Real> mean_moments;
   // m2
   Real mean_m2 = 0.5 * (moment[0][2] + moment[2][0]);
   mean_moments.push_back(mean_m2);
   // m4
   Real mean_m4 = (3 * moment[2][2] + moment[0][4] + moment[4][0]) / 3.;
   mean_moments.push_back(mean_m4);
   return mean_moments;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order, UInt k1,
                                                      UInt k2, Real hurst,
                                                      Real A, Real L) {
   Real xi = k2 / k1;
   Real exp = -2. * hurst;
   Real _k1 = 2. * M_PI / L * k1;
 
   Real Cq;
   switch (order) {
   case 0:
     Cq = 2. * M_PI;
     break;
   case 2:
     Cq = M_PI;
     break;
   case 4:
     Cq = 3. / 4. * M_PI;
     break;
   default:
     SURFACE_FATAL("invalid order");
   }
   exp = order - 2. * hurst;
   Real res = A * Cq * pow(_k1, exp) / exp * (pow(xi, exp) - 1);
 
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 }  // namespace tamaas
 
 #endif /* __SURFACE_STATISTICS_HH__ */
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index ad25f57..024bab9 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,171 +1,171 @@
 /**
  *  @mainpage Welcome to Tamaas !
  *
  *  @section Introduction
  *  Tamaas is a spectral-integral-equation based contact library. It is made
  *  with love to be fast and friendly!
  *
  *  @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *  @author Lucas Frérot <lucas.frerot@epfl.ch>
  *  @author Valentine Rey <valentine.rey@epfl.ch>
  *  @author Son Pham-Ba <son.phamba@epfl.ch>
  *
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __TAMAAS_HH__
 #define __TAMAAS_HH__
 /* -------------------------------------------------------------------------- */
 // Standard includes
 #include <exception>
 #include <iostream>
 #include <memory>
 #include <type_traits>
 #include <string>
 /* -------------------------------------------------------------------------- */
 // Special thrust includes
 #include <thrust/complex.h>
 #include <thrust/random.h>
 
 #ifdef USE_CUDA
 #include "unified_allocator.hh"
 #else
 #include "fftw_allocator.hh"
 #endif
 
 /// @section Cuda specific definitions
 #define CUDA_LAMBDA __device__ __host__
 
 #ifdef USE_CUDA
 #define DEFAULT_ALLOCATOR(T) UnifiedAllocator<T>
 #else
 #define DEFAULT_ALLOCATOR(T) FFTWAllocator<T>
 #endif
 
 /* -------------------------------------------------------------------------- */
 
 /// @section Standard definitions pulled from C++14
 namespace std {
 #if __cplusplus < 201402L
 /// exchange
 template <class T, class U = T>
 T exchange(T& obj, U&& new_value) {
   T old_value = move(obj);
   obj = forward<U>(new_value);
   return old_value;
 }
 /// make_unique
 template <typename T, typename... Args>
 unique_ptr<T> make_unique(Args&&... args) {
   return unique_ptr<T>(new T(forward<Args>(args)...));
 }
 #endif
 }  // namespace std
 
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /// @section Common types definitions
 using Real = REAL_TYPE;          ///< default floating point type
 using Int = INT_TYPE;            ///< default signed integer type
 using UInt = std::make_unsigned_t<Int>;  ///< default unsigned integer type
 template <typename T>
 using complex = thrust::complex<T>;  ///< template complex wrapper
 using Complex = complex<Real>;       ///< default floating point complex type
 static constexpr Real zero_threshold = 1e-14;
 
 /// @section Defining random toolbox
 using ::thrust::random::normal_distribution;
 using ::thrust::random::uniform_real_distribution;
 using random_engine = ::thrust::random::default_random_engine;
 
 namespace detail {
 template <bool acc, template <typename> class Trait, typename Head,
           typename... Tail>
 struct fold_trait_tail_rec
     : std::integral_constant<bool,
                              fold_trait_tail_rec<acc and Trait<Head>::value,
                                                  Trait, Tail...>::value> {};
 
 template <bool acc, template <typename> class Trait, typename Head>
 struct fold_trait_tail_rec<acc, Trait, Head>
     : std::integral_constant<bool, acc and Trait<Head>::value> {};
 }  // namespace detail
 
 template <template <typename> class Trait, typename... T>
 struct fold_trait : detail::fold_trait_tail_rec<true, Trait, T...> {};
 
 /* -------------------------------------------------------------------------- */
 /// initialize tamaas (0 threads => let OMP_NUM_THREADS decide)
 void initialize(UInt num_threads = 0);
 /// cleanup tamaas
 void finalize();
 /* -------------------------------------------------------------------------- */
 
 /// Generic exception class
 class Exception : public std::exception {
 public:
   /// Constructor
-  Exception(const std::string& mesg) : msg(mesg) {}
+  Exception(std::string mesg) : msg(std::move(mesg)) {}
 
-  virtual const char* what() const noexcept { return msg.c_str(); }
+  const char* what() const noexcept override { return msg.c_str(); }
 
-  virtual ~Exception() = default;
+  ~Exception() override = default;
 
 private:
   std::string msg;  ///< message of exception
 };
 /* -------------------------------------------------------------------------- */
 /// Enumeration of reduction operations
 enum class operation { plus, times, min, max };
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 
 /// @section Convenience macros
 #define TAMAAS_EXCEPTION(mesg)                                                 \
   {                                                                            \
     std::stringstream sstr;                                                    \
     sstr << __FILE__ << ":" << __LINE__ << ":FATAL: " << mesg << '\n';         \
     throw ::tamaas::Exception(sstr.str());                                     \
   }
 
 #define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
 
 #if defined(TAMAAS_DEBUG)
 #define TAMAAS_ASSERT(cond, reason)                                            \
   do {                                                                         \
     if (!(cond)) {                                                             \
       TAMAAS_EXCEPTION(#cond " assert failed: " << reason);                    \
     }                                                                          \
   } while (0)
 #define TAMAAS_DEBUG_EXCEPTION(reason) TAMAAS_EXCEPTION(reason)
 #else
 #define TAMAAS_ASSERT(cond, reason)
 #define TAMAAS_DEBUG_EXCEPTION(reason)
 #endif
 #define TAMAAS_ACCESSOR(var, type, name)                                       \
   type& get##name() { return var; }                                            \
   void set##name(const type& new_var) { var = new_var; }
 
 /* -------------------------------------------------------------------------- */
 
 #endif  // TAMAAS_HH
diff --git a/src/model/elasto_plastic/residual.cpp b/src/model/elasto_plastic/residual.cpp
index 89d497d..ab58269 100644
--- a/src/model/elasto_plastic/residual.cpp
+++ b/src/model/elasto_plastic/residual.cpp
@@ -1,169 +1,168 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "residual.hh"
 #include "grid_view.hh"
 #include "model_factory.hh"
 #include "model_type.hh"
 #include <list>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 template <UInt dim, UInt comp = dim>
 void symmetrize(Grid<Real, dim>& /*grad*/) {
 #if 0
   Loop::loop(
       [](MatrixProxy<Real, comp, comp> gradient) {
         // Explicit loop here to avoid temporaries
         for (UInt i = 0; i < dim; ++i) {
           for (UInt j = 0; j < dim; ++j) {
             gradient(i, j) = gradient(j, i) =
                 0.5 * (gradient(i, j) + gradient(j, i));
           }
         }
       },
       range<MatrixProxy<Real, comp, comp>>(grad));
 #endif
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 ResidualTemplate<type>::ResidualTemplate(Model* model, Real sigma_0, Real h)
     : Residual(model), hardening(model, sigma_0, h) {
 
   // Registering operators for residual and displacement calculation
   ModelFactory::registerVolumeOperators(*model);
 
   // Initializing grids
   for (auto grid : {&strain, &stress, &residual, &tmp}) {
     *grid = allocateGrid<type, false, Real>(model->getDiscretization(),
                                             trait::voigt);
   }
 
   // Registering fields
   model->registerField("stress", stress);
   model->registerField("strain", strain);
   model->registerField("residual", residual);
 
   plastic_filter = [this](UInt layer) {
 #if 0
     return true;
 #else
     return this->plastic_layers.find(layer) != this->plastic_layers.end();
 #endif
   };
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeResidual(GridBase<Real>& strain_increment) {
-  Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(strain_increment);
+  auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
   hardening.template computePlasticIncrement<false>(*residual, *strain, inc);
   this->updateFilter(*residual);
   this->model->applyElasticity(*residual, *residual);
   integralOperator("mindlin_gradient")
       .applyIf(*residual, *residual, plastic_filter);
   integralOperator("boussinesq_gradient")
       .apply(this->model->getTraction(), *tmp);
   *residual -= strain_increment;
   *residual += *tmp;
 
   symmetrize(*residual);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::applyTangent(
     GridBase<Real>& output, GridBase<Real>& input,
     GridBase<Real>& current_strain_increment) {
-  Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(current_strain_increment);
-  Grid<Real, dim>& out = dynamic_cast<decltype(out)>(output);
-  Grid<Real, dim>& in = dynamic_cast<decltype(in)>(input);
+  auto& inc = dynamic_cast<Grid<Real, dim>&>(current_strain_increment);
+  auto& out = dynamic_cast<Grid<Real, dim>&>(output);
+  auto& in = dynamic_cast<Grid<Real, dim>&>(input);
 
   hardening.applyTangentIncrement(out, in, *strain, inc);
   this->updateFilter(out);
   integralOperator("mindlin_gradient").applyIf(out, out, plastic_filter);
   out -= in;
 
   symmetrize(*residual);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeStress(GridBase<Real>& strain_increment) {
-  Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(strain_increment);
+  auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
   hardening.template computeStress<false>(*stress, *strain, inc);
   // integralOperator("boussinesq_gradient")
   //     .apply(this->model->getTraction(), tmp);
   // this->model->applyElasticity(tmp, tmp);
   // stress += tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::updateState(
     GridBase<Real>& converged_strain_increment) {
-  Grid<Real, dim>& inc =
-      dynamic_cast<decltype(inc)>(converged_strain_increment);
+  auto& inc = dynamic_cast<Grid<Real, dim>&>(converged_strain_increment);
   hardening.template computeStress<true>(*stress, *strain, inc);
   this->model->applyElasticity(*residual, hardening.getPlasticStrain());
   updateFilter(*residual);
   *strain += converged_strain_increment;
 
   // Computing displacements
   integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
                                       plastic_filter);
   Grid<Real, dim> disp_tmp(this->model->getDiscretization(), trait::components);
   integralOperator("boussinesq").apply(this->model->getTraction(), disp_tmp);
   this->model->getDisplacement() += disp_tmp;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeResidualDisplacement(
     GridBase<Real>& strain_increment) {
-  Grid<Real, dim>& inc = dynamic_cast<decltype(inc)>(strain_increment);
+  auto& inc = dynamic_cast<Grid<Real, dim>&>(strain_increment);
   hardening.template computePlasticIncrement<false>(*tmp, *strain, inc);
   updateFilter(*tmp);
   this->model->applyElasticity(*residual, *tmp);
 
   integralOperator("mindlin").applyIf(*residual, this->model->getDisplacement(),
                                       plastic_filter);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::updateFilter(
     Grid<Real, dim>& plastic_strain_increment) {
   // plastic_layers.clear();
   for (UInt i : Loop::range(plastic_strain_increment.sizes().front())) {
     auto slice = make_view(plastic_strain_increment, i);
     if (slice.dot(slice) / slice.dataSize() > 1e-14)
       plastic_layers.insert(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template class ResidualTemplate<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/model.cpp b/src/model/model.cpp
index f9e94d8..9df15c5 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,163 +1,163 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "model.hh"
 #include "be_engine.hh"
 #include "logger.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
-Model::Model(const std::vector<Real>& system_size,
-             const std::vector<UInt>& discretization)
-    : system_size(system_size), discretization(discretization) {}
+Model::Model(std::vector<Real> system_size, std::vector<UInt> discretization)
+    : system_size(std::move(system_size)),
+      discretization(std::move(discretization)) {}
 
 /* -------------------------------------------------------------------------- */
 
 Model::~Model() = default;
 
 /* -------------------------------------------------------------------------- */
 
 void Model::setElasticity(Real E, Real nu) {
   this->E = E;
   this->nu = nu;
   updateOperators();
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Model::getHertzModulus() const { return E / (1 - nu * nu); }
 
 /* -------------------------------------------------------------------------- */
 
 GridBase<Real>& Model::getTraction() { return getField("traction"); }
 
 const GridBase<Real>& Model::getTraction() const {
   return getField("traction");
 }
 
 /* -------------------------------------------------------------------------- */
 
 GridBase<Real>& Model::getDisplacement() { return getField("displacement"); }
 
 const GridBase<Real>& Model::getDisplacement() const {
   return getField("displacement");
 }
 
 /* -------------------------------------------------------------------------- */
 
 const std::vector<Real>& Model::getSystemSize() const { return system_size; }
 
 /* -------------------------------------------------------------------------- */
 
 const std::vector<UInt>& Model::getDiscretization() const {
   return discretization;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void Model::solveNeumann() {
   engine->registerNeumann();
   engine->solveNeumann(getField("traction"), getField("displacement"));
 }
 
 void Model::solveDirichlet() {
   engine->registerDirichlet();
   engine->solveDirichlet(getField("displacement"), getField("traction"));
 }
 
 /* -------------------------------------------------------------------------- */
 IntegralOperator* Model::getIntegralOperator(const std::string& name) {
   return operators[name].get();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::updateOperators() {
   for (auto& op : operators)
     op.second->updateFromModel();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumper(std::shared_ptr<ModelDumper> dumper) {
   this->dumpers.push_back(dumper);
 }
 
 void Model::dump() const {
   for (auto& dumper : dumpers)
     if (dumper)
       dumper->dump(*this);
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::registerField(const std::string& name,
                           std::shared_ptr<GridBase<Real>> field) {
   fields[name] = std::move(field);
 }
 
 const GridBase<Real>& Model::getField(const std::string& name) const try {
   return *fields.at(name);
 } catch (std::out_of_range& e) {
   Logger().get(LogLevel::warning)
       << "Field " << name << " not registered in model\n";
   throw e;
 }
 
 GridBase<Real>& Model::getField(const std::string& name) try {
   return *fields.at(name);
 } catch (std::out_of_range& e) {
   Logger().get(LogLevel::warning)
       << "Field " << name << " not registered in model\n";
   throw e;
 }
 
 std::vector<std::string> Model::getFields() const {
   std::vector<std::string> keys;
   keys.reserve(fields.size());
 
   std::transform(fields.begin(), fields.end(), std::back_inserter(keys),
                  [](auto&& pair) { return pair.first; });
   return keys;
 }
 
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this) {
   o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
     << ", nu = " << _this.getPoissonRatio() << ")\n";
 
   // Printing domain size
   o << "  - domain = [";
   std::for_each(_this.getSystemSize().begin(), _this.getSystemSize().end() - 1,
                 [&o](const Real& x) { o << x << ", "; });
   o << _this.getSystemSize().back() << "]\n";
 
   // Printing discretization
   o << "  - discretization = [";
   std::for_each(_this.getDiscretization().begin(),
                 _this.getDiscretization().end() - 1,
                 [&o](const UInt& x) { o << x << ", "; });
   o << _this.getDiscretization().back() << "]";
 
   return o;
 }
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/model.hh b/src/model/model.hh
index b258ff2..5048567 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,228 +1,228 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef MODEL_HH
 #define MODEL_HH
 /* -------------------------------------------------------------------------- */
 #include "be_engine.hh"
 #include "grid_base.hh"
 #include "integral_operator.hh"
 #include "model_dumper.hh"
 #include "model_type.hh"
 #include "tamaas.hh"
 #include <algorithm>
 #include <memory>
 #include <unordered_map>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Class for linear isotropic elasticity
  */
 namespace influence {
 template <UInt dim>
 struct ElasticHelper {
   ElasticHelper(Real mu, Real nu)
       : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
 
   template <typename DT, typename ST>
   Matrix<std::remove_cv_t<DT>, dim, dim>
   operator()(const StaticMatrix<DT, ST, dim, dim>& gradu) const {
     auto trace = gradu.trace();
     Matrix<std::remove_cv_t<DT>, dim, dim> sigma;
     for (UInt i = 0; i < dim; ++i)
       for (UInt j = 0; j < dim; ++j)
         sigma(i, j) =
             (i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
     return sigma;
   }
 
   template <typename DT, typename ST>
   SymMatrix<std::remove_cv_t<DT>, dim>
   operator()(const StaticSymMatrix<DT, ST, dim>& eps) const {
     SymMatrix<std::remove_cv_t<DT>, dim> sigma;
     auto trace = eps.trace();
     for (UInt i = 0; i < dim; ++i)
       sigma(i) = lambda * trace + 2 * mu * eps(i);
     for (UInt i = dim; i < voigt_size<dim>::value; ++i)
       sigma(i) = 2 * mu * eps(i);
     return sigma;
   }
 
   template <typename DT, typename ST>
   Matrix<std::remove_cv_t<DT>, dim, dim>
   inverse(const StaticMatrix<DT, ST, dim, dim>& sigma) const {
     Matrix<std::remove_cv_t<DT>, dim, dim> epsilon;
     auto trace = sigma.trace();
     auto c1 = 1. / (2 * mu);
     auto c2 = -lambda / (3 * lambda + 2 * mu);
 
     for (UInt i = 0; i < dim; ++i)
       for (UInt j = 0; j < dim; ++j)
         epsilon(i, j) = c1 * (sigma(i, j) + c2 * trace * (i == j));
   }
 
   const Real mu, nu, lambda;
 };
 }  // namespace influence
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Model containing pressure and displacement
  * This class is a container for the model fields. It is supposed to be
  * dimension agnostic, hence the GridBase members.
  */
 class Model {
 protected:
   /// Constructor
-  Model(const std::vector<Real>& system_size,
-        const std::vector<UInt>& discretization);
+  Model(std::vector<Real> system_size,
+        std::vector<UInt> discretization);
 
 public:
   /// Destructor
   virtual ~Model();
 
 public:
   /// Set elasticity parameters
   void setElasticity(Real E, Real nu);
   /// Get Hertz contact modulus
   Real getHertzModulus() const;
   /// Get Young's modulus
   Real getYoungModulus() const { return E; }
   /// Get Poisson's ratio
   Real getPoissonRatio() const { return nu; }
   /// Get shear modulus
   Real getShearModulus() const { return E / (2 * (1 + nu)); }
   /// Set Young's modulus
   void setYoungModulus(Real E_) { this->E = E_; }
   /// Set Poisson's ratio
   void setPoissonRatio(Real nu_) { this->nu = nu_; }
 
 public:
   virtual void applyElasticity(GridBase<Real>& stress,
                                const GridBase<Real>& strain) const = 0;
 
 public:
   /// Get pressure
   GridBase<Real>& getTraction();
   /// Get pressure
   const GridBase<Real>& getTraction() const;
   /// Get displacement
   GridBase<Real>& getDisplacement();
   /// Get displacement
   const GridBase<Real>& getDisplacement() const;
   /// Get model type
   virtual model_type getType() const = 0;
   /// Get system physical size
   const std::vector<Real>& getSystemSize() const;
   /// Get boundary system physical size
   virtual std::vector<Real> getBoundarySystemSize() const = 0;
   /// Get discretization
   const std::vector<UInt>& getDiscretization() const;
   /// Get boundary discretization
   virtual std::vector<UInt> getBoundaryDiscretization() const = 0;
   /// Get boundary element engine
   BEEngine& getBEEngine() {
     TAMAAS_ASSERT(engine, "BEEngine was not initialized");
     return *engine;
   }
 
   /// Solve Neumann problem using default neumann operator
   void solveNeumann();
   /// Solve Dirichlet problem using default dirichlet operator
   void solveDirichlet();
 
 public:
   /// Register a new integral operator
   template <typename Operator>
   IntegralOperator* registerIntegralOperator(const std::string& name) {
     operators[name] = std::make_unique<Operator>(this);
     return operators[name].get();
   }
 
   /// Get a registerd integral operator
   IntegralOperator* getIntegralOperator(const std::string& name);
   /// Tell operators to update their cache
   void updateOperators();
 
 public:
   /// Register a field
   void registerField(const std::string& name,
                      std::shared_ptr<GridBase<Real>> field);
   /// Get a field
   const GridBase<Real>& getField(const std::string& name) const;
   /// Get a non-const field
   GridBase<Real>& getField(const std::string& name);
 
   /// Get fields
   std::vector<std::string> getFields() const;
 
 public:
   /// Set the dumper object
   void addDumper(std::shared_ptr<ModelDumper> dumper);
   /// Dump the model
   void dump() const;
 
 protected:
   Real E = 1, nu = 0;
   std::vector<Real> system_size;
   std::vector<UInt> discretization;
   std::unique_ptr<BEEngine> engine = nullptr;
   std::unordered_map<std::string, std::shared_ptr<IntegralOperator>> operators;
   std::unordered_map<std::string, std::shared_ptr<GridBase<Real>>> fields;
   std::vector<std::shared_ptr<ModelDumper>> dumpers;
 };
 
 /* --------------------------------------------------------------------------
  */
 /* Output model to stream */
 /* --------------------------------------------------------------------------
  */
 std::ostream& operator<<(std::ostream& o, const Model& _this);
 
 /* --------------------------------------------------------------------------
  */
 /* Simpler grid allocation */
 /* --------------------------------------------------------------------------
  */
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model) {
   return allocateGrid<boundary, T>(
       model.getType(), (boundary) ? model.getBoundaryDiscretization()
                                   : model.getDiscretization());
 }
 
 template <bool boundary, typename T>
 std::unique_ptr<GridBase<T>> allocateGrid(Model& model, UInt nc) {
   return allocateGrid<boundary, T>(model.getType(),
                                    (boundary)
                                        ? model.getBoundaryDiscretization()
                                        : model.getDiscretization(),
                                    nc);
 }
 
 }  // namespace tamaas
 
 #endif  // MODEL_HH
diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh
index 2ec6760..32807a1 100644
--- a/src/model/volume_potential.hh
+++ b/src/model/volume_potential.hh
@@ -1,134 +1,134 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef VOLUME_POTENTIAL_HH
 #define VOLUME_POTENTIAL_HH
 /* -------------------------------------------------------------------------- */
 #include "fft_plan_manager.hh"
 #include "grid_hermitian.hh"
 #include "grid_view.hh"
 #include "integral_operator.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 #include <functional>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 /**
  * @brief Volume potential operator class. Applies the operators for computation
  * of displacements and strains due to residual/eigen strains
  */
 template <model_type type>
 class VolumePotential : public IntegralOperator {
   using trait = model_type_traits<type>;
 
 protected:
   using filter_t = std::function<bool(UInt)>;
 
 public:
   VolumePotential(Model* model);
 
   /// Update from model (does nothing)
   void updateFromModel() override {}
 
   /// Kind
   IntegralOperator::kind getKind() const override {
     return IntegralOperator::neumann;
   }
 
   /// Type
   model_type getType() const override;
 
   /// Applying to all source
   void apply(GridBase<Real>& input, GridBase<Real>& output) const override {
     applyIf(input, output, [](UInt) { return true; });
   }
 
 protected:
   /// Function to handle layer-by-layer Fourier treatment of operators
   template <typename Func>
   void fourierApply(Func func, GridBase<Real>& in, GridBase<Real>& out) const {
     fourierApplyIf(func, in, out, [](UInt) { return true; });
   }
 
   /// Function to handle layer-by-layer Fourier treatment of operators
   template <typename Func>
   void fourierApplyIf(Func func, GridBase<Real>& in, GridBase<Real>& out,
                       filter_t pred) const;
 
   /// Initialize fourier buffers
   void initialize(UInt source_components, UInt out_components);
 
 protected:
   Grid<Real, trait::boundary_dimension> wavevectors;
   using BufferType = GridHermitian<Real, trait::boundary_dimension>;
   mutable std::vector<BufferType> source_buffer;
   mutable std::vector<BufferType> disp_buffer;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation */
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 template <typename Func>
 void VolumePotential<type>::fourierApplyIf(Func func, GridBase<Real>& in,
                                            GridBase<Real>& out,
                                            filter_t pred) const {
   constexpr UInt dim = trait::dimension;
-  Grid<Real, dim>& i = dynamic_cast<decltype(i)>(in);
-  Grid<Real, dim>& o = dynamic_cast<decltype(o)>(out);
+  auto& i = dynamic_cast<Grid<Real, dim>&>(in);
+  auto& o = dynamic_cast<Grid<Real, dim>&>(out);
 
   // TAMAAS_ASSERT(i.sizes().front() == o.sizes().front(),
   //               "Number of layers does not match");
   // ^^^^ not applicable for Boussinesq
 
   // Transforming source
   for (UInt layer : Loop::range(i.sizes().front())) {
     if (not pred(layer))
       continue;
 
     auto in_layer = make_view(i, layer);
     FFTPlanManager::get()
         .createPlan(in_layer, source_buffer[layer])
         .forwardTransform();
   }
 
   // Setting output buffer to zero
   for (auto&& grid : disp_buffer)
     grid = 0;
 
   // Applying operator
   func(source_buffer, disp_buffer);
 
   // Transforming output
   for (UInt layer : Loop::range(o.sizes().front())) {
     auto out_layer = make_view(o, layer);
     FFTPlanManager::get()
         .createPlan(out_layer, disp_buffer[layer])
         .backwardTransform();
   }
 }
 
 }  // namespace tamaas
 
 #endif  // VOLUME_POTENTIAL_HH
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index e6a79ab..7a328ab 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,271 +1,271 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fft_plan_manager.hh"
 #include "grid_view.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "static_types.hh"
 
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <model_type mtype, IntegralOperator::kind otype>
 Westergaard<mtype, otype>::Westergaard(Model* model) : IntegralOperator(model) {
   // Copy sizes
   auto bdisc = model->getBoundaryDiscretization();
 
   auto boundary_hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           bdisc);
 
   constexpr UInt nb_components = trait::components;
 
   buffer.setNbComponents(nb_components);
   buffer.resize(boundary_hermitian_sizes);
   influence.setNbComponents(nb_components * nb_components);
   influence.resize(boundary_hermitian_sizes);
   initInfluence();
 }
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <UInt m, UInt j>
 struct boundary_fft_helper {
   template <typename Buffer, typename Out>
   static void backwardTransform(Buffer&& buffer, Out&& out) {
     auto view = make_view(out, 0);
     FFTPlanManager::get().createPlan(view, buffer).backwardTransform();
   }
 };
 
 template <UInt m>
 struct boundary_fft_helper<m, m> {
   template <typename Buffer, typename Out>
   static void backwardTransform(Buffer&& buffer, Out&& out) {
     FFTPlanManager::get().createPlan(out, buffer).backwardTransform();
   }
 };
 }  // namespace detail
 
 template <model_type mtype, IntegralOperator::kind otype>
 template <typename Functor>
 void Westergaard<mtype, otype>::fourierApply(Functor func, GridBase<Real>& in,
                                              GridBase<Real>& out) const try {
-  Grid<Real, bdim>& i = dynamic_cast<decltype(i)>(in);
-  Grid<Real, dim>& full_out = dynamic_cast<decltype(full_out)>(out);
+  auto& i = dynamic_cast<Grid<Real, bdim>&>(in);
+  auto& full_out = dynamic_cast<Grid<Real, dim>&>(out);
 
   FFTPlanManager::get().createPlan(i, buffer).forwardTransform();
 
   /// Applying influence
   func(buffer, influence);
 
   /// Backward fourier transform on boundary only
   detail::boundary_fft_helper<bdim, dim>::backwardTransform(buffer, full_out);
 
 } catch (const std::bad_cast& e) {
   TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type mtype, IntegralOperator::kind otype>
 template <typename Functor>
 void Westergaard<mtype, otype>::initFromFunctor(Functor func) {
   // Compute wavevectors for influence
   auto wavevectors = FFTransform<Real, bdim>::template computeFrequencies<true>(
       influence.sizes());
   // Get boundary physical size
   auto system_size = this->model->getBoundarySystemSize();
   VectorProxy<const Real, bdim> domain(system_size[0]);
 
   // Normalize wavevectors
   wavevectors *= 2 * M_PI;
   wavevectors /= domain;
 
   // Apply functor
   Loop::loop(func, range<VectorProxy<Real, bdim>>(wavevectors),
              range<MatrixProxy<Complex, comp, comp>>(influence));
 
   // Set fundamental mode to zero
   MatrixProxy<Complex, comp, comp> mat(influence(0));
   mat = 0;
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define NEUMANN_BASIC(type)                                                    \
   template <>                                                                  \
   void Westergaard<type, IntegralOperator::neumann>::initInfluence() {         \
     auto E_star = model->getHertzModulus();                                    \
     auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> q,               \
                                       MatrixProxy<Complex, comp, comp> k) {    \
       k(0, 0) = 2. / (E_star * q.l2norm());                                    \
     };                                                                         \
     initFromFunctor(basic);                                                    \
   }
 #define DIRICHLET_BASIC(type)                                                  \
   template <>                                                                  \
   void Westergaard<type, IntegralOperator::dirichlet>::initInfluence() {       \
     auto E_star = model->getHertzModulus();                                    \
     auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> q,               \
                                       MatrixProxy<Complex, comp, comp> k) {    \
       k(0, 0) = E_star * q.l2norm() / 2;                                       \
     };                                                                         \
     initFromFunctor(basic);                                                    \
   }
 
 NEUMANN_BASIC(model_type::basic_1d);
 NEUMANN_BASIC(model_type::basic_2d);
 DIRICHLET_BASIC(model_type::basic_1d);
 DIRICHLET_BASIC(model_type::basic_2d);
 
 #undef NEUMANN_BASIC
 #undef DIRICHLET_BASIC
 /* -------------------------------------------------------------------------- */
 
 template <>
 void Westergaard<model_type::surface_1d,
                  IntegralOperator::neumann>::initInfluence() {
   auto E = model->getYoungModulus();
   auto nu = model->getPoissonRatio();
   const Complex I(0, 1);
   auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> q_,
                                         MatrixProxy<Complex, comp, comp> F) {
     Real q = q_(0);
     q *= 1. / q_.l2norm();
 
     F(0, 0) = 2 * (1 + nu) * (1 - nu * q * q);
     F(1, 1) = 2 * (1 - nu * nu);
 
     F(0, 1) = I * q * (1 + nu) * (1 - 2 * nu);
     F(1, 0) = -F(0, 1);
 
     F *= 1. / (E * q_.l2norm());
   };
 
   initFromFunctor(surface);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <>
 void Westergaard<model_type::surface_2d,
                  IntegralOperator::neumann>::initInfluence() {
   auto E = model->getYoungModulus();
   auto nu = model->getPoissonRatio();
   const Complex I(0, 1);
   auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> q_,
                                         MatrixProxy<Complex, comp, comp> F) {
     Vector<Real, bdim> q(q_);
     q *= 1 / q_.l2norm();
 
     F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0));
     F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1));
     F(2, 2) = 2 * (1 - nu * nu);
 
     F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu);
 
     F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu);
     F(1, 2) = I * q(1) * (1 + nu) * (1. - 2. * nu);
 
     F(2, 0) = -F(0, 2);
     F(2, 1) = -F(1, 2);
 
     F *= 1. / (E * q_.l2norm());
   };
 
   initFromFunctor(surface);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type mtype, IntegralOperator::kind otype>
 void Westergaard<mtype, otype>::initInfluence() {}
 
 /* ------------------------------------------------------------------------ */
 template <model_type mtype, IntegralOperator::kind otype>
 void Westergaard<mtype, otype>::apply(GridBase<Real>& input,
                                       GridBase<Real>& output) const {
   auto apply = [](decltype(buffer)& buffer,
                   const decltype(influence)& influence) {
     Loop::loop(
         [] CUDA_LAMBDA(VectorProxy<Complex, comp> i,
                        MatrixProxy<const Complex, comp, comp> F) { i = F * i; },
         range<VectorProxy<Complex, comp>>(buffer),
         range<MatrixProxy<const Complex, comp, comp>>(influence));
   };
   fourierApply(apply, input, output);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, IntegralOperator::kind otype>
 Real Westergaard<type, otype>::getOperatorNorm() {
   constexpr UInt comp = model_type_traits<type>::components;
   Real _norm = 0.0;
 
   _norm = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(MatrixProxy<Complex, comp, comp> m) {
         Real n = thrust::norm(m.l2squared());
         return std::isnan(n) ? 0 : n;
       },
       range<MatrixProxy<Complex, comp, comp>>(influence));
 
   auto size = model->getSystemSize();
   auto disc = model->getDiscretization();
   auto dim = model_type_traits<type>::dimension;
 
   /// TODO: why?
   switch (dim) {
   case 1:
     _norm /= (size[0] / disc[0]) * (size[0] / disc[0]);
     break;
   default:
     for (UInt i = 0; i < dim; i++) {
       _norm /= size[i] / disc[i];
     }
   }
 
   return std::sqrt(_norm);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation */
 /* -------------------------------------------------------------------------- */
 template class Westergaard<model_type::basic_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::basic_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::basic_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::basic_2d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::surface_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::surface_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::surface_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::surface_2d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::volume_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::volume_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::volume_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::volume_2d, IntegralOperator::dirichlet>;
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
diff --git a/src/surface/isopowerlaw.hh b/src/surface/isopowerlaw.hh
index d2ef5b8..97c1e4f 100644
--- a/src/surface/isopowerlaw.hh
+++ b/src/surface/isopowerlaw.hh
@@ -1,95 +1,95 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef ISOPOWERLAW_H
 #define ISOPOWERLAW_H
 /* -------------------------------------------------------------------------- */
 #include "filter.hh"
 #include "grid_hermitian.hh"
 #include "static_types.hh"
 #include "tamaas.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Class representing an isotropic power law spectrum
 template <UInt dim>
 class Isopowerlaw : public Filter<dim> {
 public:
   /// Compute filter coefficients
-  virtual void
-  computeFilter(GridHermitian<Real, dim>& filter_coefficients) const;
+  void
+  computeFilter(GridHermitian<Real, dim>& filter_coefficients) const override;
 
   /// Compute a point of the PSD
   CUDA_LAMBDA inline Real operator()(const VectorProxy<Real, dim>& q_vec) const;
 
   /// Analytical rms of heights
   Real rmsHeights() const;
 
   /// Analytical moments
   std::vector<Real> moments() const;
 
   /// Analytical Nayak's parameter
   Real alpha() const;
 
   /// Analytical RMS slopes
   Real rmsSlopes() const;
 
 public:
   TAMAAS_ACCESSOR(q0, UInt, Q0);
   TAMAAS_ACCESSOR(q1, UInt, Q1);
   TAMAAS_ACCESSOR(q2, UInt, Q2);
   TAMAAS_ACCESSOR(hurst, Real, Hurst);
 
 protected:
   UInt q0, q1, q2;
   Real hurst;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline implementations                                                     */
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 CUDA_LAMBDA inline Real Isopowerlaw<dim>::
 operator()(const VectorProxy<Real, dim>& q_vec) const {
   const Real C = 1.;
   const Real q = q_vec.l2norm();
 
   Real val;
 
   if (q < q0)
     val = 0;
   else if (q > q2)
     val = 0;
   else if (q < q1)
     val = C;
   else
     val = C * std::pow((q / q1), -2 * (hurst + 1));
 
   return std::sqrt(val);
 }
 
 }  // namespace tamaas
 
 #endif  // ISOPOWERLAW_H
diff --git a/src/surface/regularized_powerlaw.hh b/src/surface/regularized_powerlaw.hh
index c4e5764..b45f2cd 100644
--- a/src/surface/regularized_powerlaw.hh
+++ b/src/surface/regularized_powerlaw.hh
@@ -1,78 +1,78 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef REGULARIZED_POWERLAW_H
 #define REGULARIZED_POWERLAW_H
 /* -------------------------------------------------------------------------- */
 #include "filter.hh"
 #include "grid_hermitian.hh"
 #include "static_types.hh"
 #include "tamaas.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Class representing an isotropic power law spectrum
 template <UInt dim>
 class RegularizedPowerlaw : public Filter<dim> {
 public:
   /// Compute filter coefficients
-  virtual void
-  computeFilter(GridHermitian<Real, dim>& filter_coefficients) const;
+  void
+  computeFilter(GridHermitian<Real, dim>& filter_coefficients) const override;
 
   /// Compute a point of the PSD
   CUDA_LAMBDA inline Real operator()(const VectorProxy<Real, dim>& q_vec) const;
 
 public:
   TAMAAS_ACCESSOR(q1, UInt, Q1);
   TAMAAS_ACCESSOR(q2, UInt, Q2);
   TAMAAS_ACCESSOR(hurst, Real, Hurst);
 
 protected:
   UInt q1, q2;
   Real hurst;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline implementations                                                     */
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 CUDA_LAMBDA inline Real RegularizedPowerlaw<dim>::
 operator()(const VectorProxy<Real, dim>& q_vec) const {
   const Real C = 1.;
   const Real q = q_vec.l2norm();
 
   Real val;
 
   if (q > q2)
     val = 0;
   else
     val = C * std::pow(1 + (q / q1) * (q / q1), -(hurst + 1));
 
   return std::sqrt(val);
 }
 
 }  // namespace tamaas
 
 #endif  // REGULARIZED_POWERLAW_H
diff --git a/src/surface/surface_generator_filter_fft.cpp b/src/surface/surface_generator_filter_fft.cpp
index 25abed9..dc973f0 100644
--- a/src/surface/surface_generator_filter_fft.cpp
+++ b/src/surface/surface_generator_filter_fft.cpp
@@ -1,57 +1,55 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #include "surface_generator_filter_fft.hh"
 #include "isopowerlaw.hh"
 #include "logger.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 extern template class Grid<Real, 2>;
 
 SurfaceGeneratorFilterFFT::SurfaceGeneratorFilterFFT() {
   pl_filter = std::make_shared<Isopowerlaw<2>>();
   generator.setFilter(pl_filter);
 }
 
-SurfaceGeneratorFilterFFT::~SurfaceGeneratorFilterFFT() {}
-
 Real SurfaceGeneratorFilterFFT::constrainRMS() {
   Logger().get(LogLevel::warning)
       << "Statistical bias in using the generated RMS as "
          "normalizing factor"
       << std::endl;
   Real hrms_grid = SurfaceStatistics::computeStdev(generator.grid);
   generator.grid *= rms / hrms_grid;
   return rms;
 }
 
 Grid<Real, 2>& SurfaceGeneratorFilterFFT::buildSurface() {
   generator.setSizes({{grid_size, grid_size}});
   return generator.buildSurface();
 }
 
 }  // namespace tamaas
diff --git a/src/surface/surface_generator_filter_fft.hh b/src/surface/surface_generator_filter_fft.hh
index 32d0871..fc6cfa2 100644
--- a/src/surface/surface_generator_filter_fft.hh
+++ b/src/surface/surface_generator_filter_fft.hh
@@ -1,75 +1,73 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef SURFACE_GENERATOR_FILTER_FFT_H
 #define SURFACE_GENERATOR_FILTER_FFT_H
 /* -------------------------------------------------------------------------- */
 #include "isopowerlaw.hh"
 #include "surface_generator_filter.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Legacy class to maintain interface
 class SurfaceGeneratorFilterFFT {
 public:
   //! can construct a generator with default values
   SurfaceGeneratorFilterFFT();
-  //! destructor
-  virtual ~SurfaceGeneratorFilterFFT();
 
 public:
   //! get high frequency cut
   UInt& getQ2() { return pl_filter->getQ2(); }
   //! get low frequency cut
   UInt& getQ0() { return pl_filter->getQ0(); }
   //! get roll frequency cut
   UInt& getQ1() { return pl_filter->getQ1(); }
   //! get hurst
   Real& getHurst() { return pl_filter->getHurst(); }
   /// compute analytical RMS
   Real analyticalRMS() const { return pl_filter->rmsHeights(); }
   /// Swig does not wrap the numpy otherwise
   Grid<Real, 2>& buildSurface();
   /// Get gridsize
   UInt& getGridSize() { return grid_size; }
   /// Get random seed
   long& getRandomSeed() { return generator.getRandomSeed(); }
   /// Get RMS
   Real& getRMS() { return rms; }
   /// Contrain RMS
   Real constrainRMS();
   /// Init
   void Init() {}
 
 protected:
   SurfaceGeneratorFilter<2> generator;
   std::shared_ptr<Isopowerlaw<2>> pl_filter = nullptr;
 
   UInt grid_size;
   Real rms;
 };
 
 }  // namespace tamaas
 
 #endif
diff --git a/tests/test_loop.cpp b/tests/test_loop.cpp
index 2f27dee..7968df2 100644
--- a/tests/test_loop.cpp
+++ b/tests/test_loop.cpp
@@ -1,362 +1,362 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "grid_view.hh"
 #include "loop.hh"
 #include "static_types.hh"
 #include "test.hh"
 /* -------------------------------------------------------------------------- */
 /* WARNING: here we cannot use lambdas for tests because GoogleTest declares  */
 /* test functions as private members of classes which is incompatible with    */
 /* cuda's extended lambdas. I know... it's f*cking stupid                     */
 /* -------------------------------------------------------------------------- */
 
 using namespace tamaas;
 
 template <typename T>
 struct AddOneInplace {
   CUDA_LAMBDA void operator()(T& x) { x += 1; }
 };
 
 // Testing loops on one grid
 TEST(TestLoops, OneArgument) {
   Grid<Real, 1> grid({20}, 1);
   Grid<Real, 1> solution({20}, 1);
 
   auto add_one = [](Real& x) { return x + 1; };
 
   std::iota(grid.begin(), grid.end(), 1);
 
   // Makeing solution
   std::transform(grid.begin(), grid.end(), solution.begin(), add_one);
 
   auto add_one_inplace = AddOneInplace<Real>();
 
   Loop::loop(add_one_inplace, grid);
 
   ASSERT_TRUE(compare(grid, solution, AreFloatEqual()))
       << "One argument loop failed";
 }
 
 struct PrimalTest {
   CUDA_LAMBDA void operator()(Int& primal, Int& val) {
     val = (primal > 0) ? -1 : 1;
   }
 };
 
 // Testing loops on two grids
 TEST(TestLoops, TwoArguments) {
   // Why no ints?
   Grid<Int, 2> grid({20, 20}, 1);
   Grid<Int, 2> primal({20, 20}, 1);
   Grid<Int, 2> solution({20, 20}, 1);
 
   primal(0, 0) = 1;
   primal(0, 1) = 1;
   primal(1, 0) = 1;
   primal(1, 1) = 1;
 
   std::transform(primal.begin(), primal.end(), solution.begin(),
                  [](Int& primal) { return (primal > 0) ? -1 : 1; });
 
   auto primal_test = PrimalTest();
 
   Loop::loop(primal_test, primal, grid);
 
   ASSERT_TRUE(compare(solution, grid)) << "Two argument loop failed";
 }
 
 struct AssignUInt {
   CUDA_LAMBDA void operator()(UInt& x, UInt i) { x = i; }
 };
 
 // Testing an enumeration
 TEST(TestLoops, Enumeration) {
   Grid<UInt, 1> grid({100}, 1);
   Grid<UInt, 1> solution({100}, 1);
   std::iota(solution.begin(), solution.end(), 0);
 
   auto assign_uint = AssignUInt();
 
   Loop::loop(assign_uint, grid, Loop::range(100));
 
   ASSERT_TRUE(compare(solution, grid)) << "Enumeration loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 struct Identity {
   CUDA_LAMBDA UInt operator()(UInt& x) const { return x; }
 };
 
 // Testing one grid reductions
 TEST(TestReductions, OneArgument) {
   Grid<UInt, 1> grid({6}, 1);
   std::iota(grid.begin(), grid.end(), 1);
 
   const auto id = Identity();
 
   // Sum reduction
-  UInt sol = std::accumulate(grid.begin(), grid.end(), 0, std::plus<UInt>());
+  UInt sol = std::accumulate(grid.begin(), grid.end(), 0, std::plus<>());
   UInt red = Loop::reduce<operation::plus>(id, grid);
   ASSERT_TRUE(sol == red) << "Addition reduction failed on one argument";
 
   // Product reduction
-  sol = std::accumulate(grid.begin(), grid.end(), 1, std::multiplies<UInt>());
+  sol = std::accumulate(grid.begin(), grid.end(), 1, std::multiplies<>());
   red = Loop::reduce<operation::times>(id, grid);
   ASSERT_TRUE(sol == red) << "Multiplication reduction failed on one argument";
 
   // Min reduction
   sol = *std::min_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::min>(id, grid);
   ASSERT_TRUE(sol == red) << "Min reduction failed on one argument";
 
   // Max reduction
   sol = *std::max_element(grid.begin(), grid.end());
   red = Loop::reduce<operation::max>(id, grid);
   ASSERT_TRUE(sol == red) << "Max reduction failed on one argument";
 }
 
 TEST(TestReductions, ReduceAndTransform) {
   UInt n = 20;
   Grid<UInt, 1> grid({n}, 1), solution({n}, 1);
   std::iota(solution.begin(), solution.end(), 0);
   UInt sum_value = (n - 1) * n / 2;
 
   UInt res = Loop::reduce<operation::plus>(
       [](UInt& x, UInt i) {
         x = i;
         return x;
       },
       grid, Loop::range(n));
 
   EXPECT_EQ(res, sum_value) << "Reduction failed";
   EXPECT_TRUE(compare(grid, solution)) << "Assign failed";
 }
 
 struct PrimalReduce {
   CUDA_LAMBDA UInt operator()(UInt& p, UInt& val) { return (p > 0) ? val : 0; }
 };
 
 TEST(TestReductions, TwoArguments) {
   Grid<UInt, 1> grid({20}, 1);
   Grid<UInt, 1> primal({20}, 1);
   grid = 1;
   primal(0) = 1;
   primal(1) = 1;
 
   auto primal_reduce = PrimalReduce();
 
   // Reduce on values where primal > 0
   UInt red = Loop::reduce<operation::plus>(primal_reduce, primal, grid);
   ASSERT_TRUE(red == 2) << "Two args reduction failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestRange, type_trait) {
   Grid<UInt, 1> grid({1}, 1);
   auto gridrange = range<VectorProxy<UInt, 1>>(grid);
 
   static_assert(decltype(gridrange)::is_valid_container<Grid<UInt, 1>>::value,
                 "is_valid_container Type trait is wrong");
 
   static_assert(
       not decltype(gridrange)::is_valid_container<Grid<Real, 1>&>::value,
       "is_valid_container Type trait is wrong");
 
   static_assert(not Range<VectorProxy<Real, 1>, Real,
                           1>::is_valid_container<decltype(grid)>::value,
                 "is_valid_container Type trait is wrong");
 }
 
 TEST(TestRange, headless) {
   Grid<UInt, 1> grid({10}, 1), solution({10}, 1);
   std::fill(++solution.begin(), solution.end(), 1);
   auto gridrange = range<VectorProxy<UInt, 1>>(grid).headless();
   Loop::loop([](auto x) { x = 1; }, gridrange);
   ASSERT_TRUE(compare(grid, solution)) << "Headless fail";
 }
 
 template <typename T>
 using WrapVector = VectorProxy<T, 2>;
 
 struct AddOneVector {
   CUDA_LAMBDA void operator()(WrapVector<UInt> x) { x(0) += 1; }
 };
 
 TEST(TestStridedLoops, VectorStride) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 1);
 
   Grid<UInt, 2> solution({10, 10}, 2);
   solution = grid;
 
   std::for_each(solution.begin(), solution.end(), [](UInt& x) {
     if (x % 2 == 1)
       x += 1;
   });
 
   auto add_one_inplace = AddOneVector();
 
   Loop::loop(add_one_inplace, range<WrapVector<UInt>>(grid));
 
   ASSERT_TRUE(compare(solution, grid)) << "Static vector strided loop failed";
 }
 
 template <typename T>
 using WrapMatrix = MatrixProxy<T, 2, 2>;
 
 struct SetOneMatrix {
   CUDA_LAMBDA void operator()(WrapMatrix<UInt> x) {
     x(0, 0) = 1;
     x(1, 1) = 1;
   }
 };
 
 TEST(TestStridedLoops, MatrixStride) {
   Grid<UInt, 2> grid({10, 10}, 4);
   Grid<UInt, 2> solution({10, 10}, 4);
   std::iota(solution.begin(), solution.end(), 0);
 
   std::for_each(solution.begin(), solution.end(), [](UInt& x) {
     if (x % 4 == 0 || x % 4 == 3)
       x = 1;
     else
       x = 0;
   });
 
   auto set_one = SetOneMatrix();
 
   Loop::loop(set_one, range<WrapMatrix<UInt>>(grid));
 
   ASSERT_TRUE(compare(solution, grid)) << "Static matrix strided loop failed";
 }
 
 struct VectorReduction {
   CUDA_LAMBDA Vector<UInt, 3> operator()(VectorProxy<UInt, 3> v) const {
     return v;
   }
 };
 
 struct BroadcastSet123 {
   CUDA_LAMBDA inline void operator()(VectorProxy<UInt, 3> v) const {
     v(0) = 1;
     v(1) = 2;
     v(2) = 3;
   }
 };
 
 TEST(TestStridedReduction, VectorReduce) {
   Grid<UInt, 2> grid({10, 10}, 3);
 
   Loop::loop(BroadcastSet123(), range<VectorProxy<UInt, 3>>(grid));
 
   auto res = Loop::reduce<operation::plus>(VectorReduction(),
                                            range<VectorProxy<UInt, 3>>(grid));
 
   ASSERT_EQ(res(0), 100);
   ASSERT_EQ(res(1), 200);
   ASSERT_EQ(res(2), 300);
 }
 
 TEST(TestViewReduction, ScalarReduce) {
   Grid<UInt, 2> grid({10, 10}, 3);
   Loop::loop(BroadcastSet123(), range<VectorProxy<UInt, 3>>(grid));
 
   auto view = make_component_view(grid, 2);
 
   UInt res = Loop::reduce<operation::plus>([](UInt& x) { return x; }, view);
   EXPECT_EQ(res, 300) << "Reduce on component view fail";
 }
 
 TEST(TestViewReduction, VectorReduce) {
   Grid<UInt, 2> grid({10, 10}, 3);
   auto view2 = make_view(grid, 0);
 
   Loop::loop(BroadcastSet123(), range<VectorProxy<UInt, 3>>(view2));
 
   auto res2 = Loop::reduce<operation::plus>(VectorReduction(),
                                             range<VectorProxy<UInt, 3>>(view2));
   EXPECT_EQ(res2(0), 10);
   EXPECT_EQ(res2(1), 20);
   EXPECT_EQ(res2(2), 30);
 }
 
 TEST(TestViewLoop, ScalarLoop) {
   Grid<UInt, 2> grid({10, 10}, 3), solution({10, 10}, 3);
   auto view = make_component_view(grid, 2);
 
   Loop::loop([](auto& x) { x = 1; }, view);
 
   Loop::loop([](auto v) { v(2) = 1; }, range<VectorProxy<UInt, 3>>(solution));
 
   ASSERT_TRUE(compare(grid, solution)) << "View loop fail";
 }
 
 TEST(TestLoopChecks, Components) {
   Grid<UInt, 2> grid({10, 10}, 3);
 
   EXPECT_THROW(
       Loop::loop([](auto v) { v = 0; }, range<VectorProxy<UInt, 2>>(grid)),
       Exception)
       << "Broken check on number of components";
 }
 
 TEST(TestLoopChecks, LoopSize) {
   Grid<UInt, 1> grid({10}, 2), other({10}, 1);
 
   EXPECT_THROW(Loop::loop([](auto& x, auto& y) { x = y; }, grid, other),
                Exception)
       << "Check on loop size without ranges fail";
   other.resize({11});
   EXPECT_THROW(Loop::loop([](auto x, auto y) { x(0) = y(0); },
                           range<VectorProxy<UInt, 2>>(grid),
                           range<VectorProxy<UInt, 1>>(other)),
                Exception)
       << "Check on loop size with ranges fail";
 
   Grid<UInt, 2> twod({10, 11}, 2);
   auto view = make_view(twod, 0);
   EXPECT_THROW(Loop::loop([](auto& x, auto& y) { x = y; }, grid, view),
                Exception)
       << "Check on loop size with view fail";
 }
 
 TEST(TestReductions, ReduceAndTransformVector) {
   UInt n = 20;
   Grid<UInt, 1> grid({n}, 2), solution({n}, 2);
   std::iota(solution.begin(), solution.end(), 1);
   std::iota(grid.begin(), grid.end(), 0);
   UInt sum_value = (2 * n + 1) * 2 * n / 2;
 
   UInt res = Loop::reduce<operation::plus>(
       [](auto x, UInt i) {
         x += 1;
         return x(0) + x(1);
       },
       range<VectorProxy<UInt, 2>>(grid), Loop::range(n));
 
   EXPECT_EQ(res, sum_value) << "Reduction failed";
   EXPECT_TRUE(compare(grid, solution)) << "Assign failed";
 }