diff --git a/.clang-tidy b/.clang-tidy
index 15ca48f..95f00f1 100644
--- a/.clang-tidy
+++ b/.clang-tidy
@@ -1,3 +1,3 @@
-Checks:          '-*, modernize-use-*, -modernize-use-trailing-return-type*, performance-*, readability-*, -readability-braces-around-statements, -readability-implicit-bool-conversion, clang-analyzer-*, -clang-analyzer-core.NonNullParamChecker'
+Checks:          '-*, modernize-use-*, -modernize-use-trailing-return-type*, performance-*, readability-*, -readability-magic-numbers, -readability-isolate-declaration, -readability-braces-around-statements, -readability-implicit-bool-conversion, clang-analyzer-*, -clang-analyzer-core.NonNullParamChecker'
 AnalyzeTemporaryDtors: false
 FormatStyle:     file
diff --git a/src/core/array.hh b/src/core/array.hh
index dd648cc..9c6c7d3 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,175 +1,171 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 ARRAY_HH
 #define ARRAY_HH
 /* -------------------------------------------------------------------------- */
 #include "logger.hh"
 #include "tamaas.hh"
 #include <memory>
 #include <thrust/copy.h>
 #include <thrust/fill.h>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
-/**
- * @brief Generic storage class with wrapping capacities
- */
+/// Generic storage class with wrapping capacities
 template <typename T>
-class Array final {
-public:
+struct Array final {
   /// Default
   Array() = default;
   /// Empty array of given size
   Array(UInt size) : Array() { this->resize(size); }
   /// Copy constructor (deep)
   Array(const Array& v) : Array() {
     this->resize(v._size);
     thrust::copy_n(v._data, v._size, this->_data);
   }
   /// Move constructor (transfers data ownership)
   Array(Array&& v) noexcept : Array() {
     _data = std::exchange(v._data, nullptr);
     _size = std::exchange(v._size, 0);
     _reserved = std::exchange(v._reserved, 0);
     wrapped = std::exchange(v.wrapped, false);
   }
   /// Wrap array on data
   Array(T* data, UInt size) : Array() { this->wrapMemory(data, size); }
   /// Destructor
   ~Array() {
     if (not wrapped)
       _alloc.deallocate(_data, _reserved);
   }
 
-public:
   /// Copy operator
   Array& operator=(const Array& v) {
     this->wrapped = false;
     if (v.size() != this->size())
       this->resize(v.size());
     thrust::copy_n(v._data, _size, _data);
     return *this;
   }
 
   /// Move operator
   Array& operator=(Array&& v) noexcept {
     if (this != &v) {
       if (not wrapped)
         _alloc.deallocate(_data, _reserved);
       _data = std::exchange(v._data, nullptr);
       _size = std::exchange(v._size, 0);
       _reserved = std::exchange(v._reserved, 0);
       wrapped = std::exchange(v.wrapped, false);
     }
     return *this;
   }
 
   void wrap(Array& other) { this->wrapMemory(other._data, other._size); }
 
   void wrap(const Array& other) { this->wrapMemory(other._data, other._size); }
 
   /// Wrap a memory pointer
   void wrapMemory(T* data, UInt size) {
     _data = data;
     _size = size;
     wrapped = true;
     _reserved = 0;
   }
 
   /// Assign a sigle value to the array
   void assign(UInt size, const T& value) {
     this->resize(size);
     thrust::fill(_data, _data + size, value);
   }
 
   void setWrapped(bool w) { wrapped = w; }
 
   /// Data pointer access (const)
   const T* data() const { return _data; }
   /// Data pointer access (non-const)
   T* data() { return _data; }
   /// Data pointer setter
   void setData(T* new_ptr) { _data = new_ptr; }
   /// Resize array
   void resize(UInt size) {
     if (wrapped)
       TAMAAS_EXCEPTION("cannot resize wrapped array");
 
     // Erase array
     if (size == 0) {
       _alloc.deallocate(_data, _reserved);
       _data = nullptr;
       _size = 0;
       _reserved = 0;
       return;
     }
 
     // Do nothing
     if (size == this->_size)
       return;
 
     // Allocate new data
     _alloc.deallocate(_data, _reserved);
 
     _data = _alloc.allocate(size);
     _size = size;
     _reserved = size;
   }
   /// Reserve storage space
   void reserve(UInt size) {
     if (_reserved >= size)
       return;
 
     auto new_data = _alloc.allocate(size);
     if (new_data != _data) {
       thrust::copy_n(_data, _size, new_data);
       _alloc.deallocate(_data, _reserved);
       _data = new_data;
       _reserved = size;
     }
   }
   /// Access operator
   inline T& operator[](UInt i) {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Access operator (const)
   inline const T& operator[](UInt i) const {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Get size of array
   inline UInt size() const { return _size; }
 
 private:
   T* _data = nullptr;
   UInt _size = 0;
   std::size_t _reserved = 0;
   bool wrapped = false;
   Allocator<T> _alloc;
 };
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif /* ARRAY_HH */
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index ec9ef2d..83f5baf 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,344 +1,341 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 GRID_BASE_HH
 #define GRID_BASE_HH
 /* -------------------------------------------------------------------------- */
 #include "array.hh"
 #include "iterator.hh"
 #include "loop.hh"
 #include "mpi_interface.hh"
 #include "static_types.hh"
 #include "tamaas.hh"
 #include <cstddef>
 #include <limits>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @brief Dimension agnostic grid with components stored per points
  */
 template <typename T>
 class GridBase {
 public:
   /// Constructor by default
   GridBase() = default;
   /// Copy constructor
   GridBase(const GridBase& o) { this->copy(o); }
   /// Move constructor (transfers data ownership)
   GridBase(GridBase&& o) noexcept
       : data(std::move(o.data)),
         nb_components(std::exchange(o.nb_components, 1)) {}
   /// Initialize with size
   explicit GridBase(UInt nb_points, UInt nb_components = 1) {
     setNbComponents(nb_components);
     this->resize(nb_points * nb_components);
   }
 
   /// Destructor
   virtual ~GridBase() = default;
 
   /* ------------------------------------------------------------------------ */
   /* Iterator class */
   /* ------------------------------------------------------------------------ */
-public:
   using iterator = iterator_::iterator<T>;
   using const_iterator = iterator_::iterator<const T>;
   using value_type = T;
   using reference = value_type&;
 
-public:
   /// Get grid dimension
   virtual UInt getDimension() const { return 1; }
   /// Get internal data pointer (const)
   const T* getInternalData() const { return this->data.data(); }
   /// Get internal data pointer (non-const)
   T* getInternalData() { return this->data.data(); }
   /// Get number of components
   UInt getNbComponents() const { return nb_components; }
   /// Set number of components
   void setNbComponents(UInt n) { nb_components = n; }
   /// Get total size
   virtual UInt dataSize() const { return this->data.size(); }
   /// Get global data size
   UInt globalDataSize() const {
     return mpi::allreduce<operation::plus>(dataSize());
   }
   /// Get number of points
   UInt getNbPoints() const {
     return this->dataSize() / this->getNbComponents();
   }
   /// Get global number of points
   UInt getGlobalNbPoints() const {
     return this->globalDataSize() / this->getNbComponents();
   }
   /// Set components
   void uniformSetComponents(const GridBase<T>& vec);
   /// Resize
   void resize(UInt size) { this->data.assign(size, T(0.)); }
   /// Reserve storage w/o changing data logic
   void reserve(UInt size) { this->data.reserve(size); }
 
   /* ------------------------------------------------------------------------ */
   /* Iterator methods */
   /* ------------------------------------------------------------------------ */
   virtual iterator begin(UInt n = 1) {
     return iterator(this->getInternalData(), n);
   }
 
   virtual iterator end(UInt n = 1) {
     return iterator(this->getInternalData() + this->dataSize(), n);
   }
 
   virtual const_iterator begin(UInt n = 1) const {
     return const_iterator(this->getInternalData(), n);
   }
 
   virtual const_iterator end(UInt n = 1) const {
     return const_iterator(this->getInternalData() + this->dataSize(), n);
   }
 
   /* ------------------------------------------------------------------------ */
   /* Operators */
   /* ------------------------------------------------------------------------ */
   inline T& operator()(UInt i) { return this->data[i]; }
 
   inline const T& operator()(UInt i) const { return this->data[i]; }
 
 #define VEC_OPERATOR(op)                                                       \
   template <typename T1>                                                       \
   GridBase& operator op(const GridBase<T1>& other)
 
   VEC_OPERATOR(+=);
   VEC_OPERATOR(*=);
   VEC_OPERATOR(-=);
   VEC_OPERATOR(/=);
 #undef VEC_OPERATOR
 
   template <typename T1>
   T dot(const GridBase<T1>& other) const;
 
 #define SCALAR_OPERATOR(op) GridBase& operator op(const T& e)
 
   SCALAR_OPERATOR(+=);
   SCALAR_OPERATOR(*=);
   SCALAR_OPERATOR(-=);
   SCALAR_OPERATOR(/=);
   SCALAR_OPERATOR(=);
 #undef SCALAR_OPERATOR
 
 /// Broadcast operators
 #define BROADCAST_OPERATOR(op)                                                 \
   template <typename DT, typename ST, UInt size>                               \
   GridBase& operator op(const StaticArray<DT, ST, size>& b)
 
   BROADCAST_OPERATOR(+=);
   BROADCAST_OPERATOR(-=);
   BROADCAST_OPERATOR(*=);
   BROADCAST_OPERATOR(/=);
   BROADCAST_OPERATOR(=);
 #undef BROADCAST_OPERATOR
 
   /* ------------------------------------------------------------------------ */
   /* Min/max */
   /* ------------------------------------------------------------------------ */
   T min() const;
   T max() const;
   T sum() const;
   T mean() const { return sum() / static_cast<T>(globalDataSize()); }
   T var() const;
   T l2norm() const { return std::sqrt(dot(*this)); }
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy */
   /* ------------------------------------------------------------------------ */
 
-public:
   GridBase& operator=(const GridBase& o) {
     this->copy(o);
     return *this;
   }
 
   GridBase& operator=(GridBase& o) {
     this->copy(o);
     return *this;
   }
 
   GridBase& operator=(GridBase&& o) noexcept {
     this->move(std::move(o));
     return *this;
   }
 
   template <typename T1>
   void copy(const GridBase<T1>& other) {
     if (other.dataSize() != this->dataSize())
       this->resize(other.dataSize());
     thrust::copy(other.begin(), other.end(), this->begin());
     nb_components = other.nb_components;
   }
 
   template <typename T1>
   void move(GridBase<T1>&& other) noexcept {
     data = std::move(other.data);
     nb_components = std::exchange(other.nb_components, 1);
   }
 
   GridBase& wrap(GridBase& o) {
     data.wrap(o.data);
     nb_components = o.nb_components;
     return *this;
   }
 
   GridBase& wrap(const GridBase& o) {
     data.wrap(o.data);
     nb_components = o.nb_components;
     return *this;
   }
 
 protected:
   Array<T> data;
   UInt nb_components = 1;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void GridBase<T>::uniformSetComponents(const GridBase<T>& vec) {
   if (vec.dataSize() != this->nb_components)
     TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
 
   auto begin_it{begin()};
   auto end_it{end()};  // silencing nvcc warnings
 
   for (auto it = begin(); it < end_it; ++it) {
     UInt i = it - begin_it;
     *it = vec(i % this->nb_components);
   }
 }
 
 /* --------------------------------------------------------------------------
  */
 
 #define SCALAR_OPERATOR_IMPL(op)                                               \
   template <typename T>                                                        \
   inline GridBase<T>& GridBase<T>::operator op(const T& e) {                   \
     Loop::loop([e] CUDA_LAMBDA(T& val) { val op e; }, *this);                  \
     return *this;                                                              \
   }
 
 SCALAR_OPERATOR_IMPL(+=)
 SCALAR_OPERATOR_IMPL(*=)
 SCALAR_OPERATOR_IMPL(-=)
 SCALAR_OPERATOR_IMPL(/=)
 SCALAR_OPERATOR_IMPL(=)
 
 #undef SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 #define VEC_OPERATOR_IMPL(op)                                                  \
   template <typename T>                                                        \
   template <typename T1>                                                       \
   inline GridBase<T>& GridBase<T>::operator op(const GridBase<T1>& other) {    \
     TAMAAS_ASSERT(other.dataSize() == this->dataSize(),                        \
                   "surface size does not match");                              \
     Loop::loop([] CUDA_LAMBDA(T& x, const T1& y) { x op y; }, *this, other);   \
     return *this;                                                              \
   }
 
 VEC_OPERATOR_IMPL(+=)
 VEC_OPERATOR_IMPL(-=)
 VEC_OPERATOR_IMPL(*=)
 VEC_OPERATOR_IMPL(/=)
 #undef VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 #define BROADCAST_OPERATOR_IMPL(op)                                            \
   template <typename T>                                                        \
   template <typename DT, typename ST, UInt size>                               \
   GridBase<T>& GridBase<T>::operator op(const StaticArray<DT, ST, size>& b) {  \
     TAMAAS_ASSERT(this->dataSize() % size == 0,                                \
                   "Broadcast operator cannot broadcast"                        \
                       << size << " on array of size " << this->dataSize());    \
     Tensor<StaticArray, typename std::remove_const<DT>::type, size> a(b);      \
     Loop::loop([a] CUDA_LAMBDA(VectorProxy<T, size> g) { g op a; },            \
                range<VectorProxy<T, size>>(*this));                            \
     return *this;                                                              \
   }
 
 BROADCAST_OPERATOR_IMPL(+=)
 BROADCAST_OPERATOR_IMPL(-=)
 BROADCAST_OPERATOR_IMPL(*=)
 BROADCAST_OPERATOR_IMPL(/=)
 BROADCAST_OPERATOR_IMPL(=)
 #undef BROADCAST_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::min() const {
   // const auto id = [] CUDA_LAMBDA (const T&x){return x;};
   return Loop::reduce<operation::min>([] CUDA_LAMBDA(const T& x) { return x; },
                                       *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::max() const {
   return Loop::reduce<operation::max>([] CUDA_LAMBDA(const T& x) { return x; },
                                       *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::sum() const {
   return Loop::reduce<operation::plus>([] CUDA_LAMBDA(const T& x) { return x; },
                                        *this);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::var() const {
   const T mu = mean();
   const T var = Loop::reduce<operation::plus>(
       [mu] CUDA_LAMBDA(const T& x) { return (x - mu) * (x - mu); }, *this);
   return var / (globalDataSize() - 1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <typename T1>
 T GridBase<T>::dot(const GridBase<T1>& other) const {
   return Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(const T& x, const T1& y) { return x * y; }, *this, other);
 }
 
 }  // namespace tamaas
 
 #endif  // GRID_BASE_HH
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index d84c453..739dd07 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,211 +1,211 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 "statistics.hh"
 #include "fft_engine.hh"
 #include "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <UInt dim>
 Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
   return std::sqrt(surface.var());
 }
 
 template <UInt dim>
 Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
   wavevectors *= 2 * M_PI;  // need q for slopes
   const auto psd = computePowerSpectrum(surface);
 
   const Real rms_slope_mean = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(VectorProxy<const Real, dim> q, const Complex& psd_val) {
         // Checking if we're in the zone that does not have hermitian symmetry
         if (std::abs(q.back()) < 1e-15)
           return q.l2squared() * psd_val.real();
         return 2 * q.l2squared() * psd_val.real();
       },
       range<VectorProxy<const Real, dim>>(wavevectors), psd);
 
   return std::sqrt(rms_slope_mean);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 GridHermitian<Real, dim>
 Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   GridHermitian<Real, dim> psd(h_size, surface.getNbComponents());
 
   FFTEngine::makeEngine()->forward(surface, psd);
 
   Real factor = 1. / surface.getGlobalNbPoints();
 
   // Squaring the fourier transform of surface and normalizing
   Loop::loop(
       [factor] CUDA_LAMBDA(Complex & c) {
         c *= factor;
         c *= conj(c);
       },
       psd);
 
   return psd;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Grid<Real, dim>
 Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
   Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
   auto psd = computePowerSpectrum(surface);
 
   FFTEngine::makeEngine()->backward(acf, psd);
   acf *= acf.getGlobalNbPoints();
   return acf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 Real Statistics<dim>::contact(const Grid<Real, dim>& tractions,
                               UInt perimeter) {
   Real points = 0;
   UInt nc = tractions.getNbComponents();
 
   switch (nc) {
   case 1:
     points = Loop::reduce(
         [] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions);
     break;
   case 2:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 2> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 2>>(tractions));
     break;
   case 3:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 3> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 3>>(tractions));
     break;
   default:
     TAMAAS_EXCEPTION("Invalid number of components in traction");
   }
 
   auto area = points / tractions.getGlobalNbPoints();
 
   if (dim == 1)
     perimeter = 0;
 
   // Correction from Yastrebov et al. (Trib. Intl., 2017)
   // 10.1016/j.triboint.2017.04.023
   return area - (M_PI - 1 + std::log(2)) /
                     (24. * tractions.getGlobalNbPoints()) * perimeter;
 }
 
 /* -------------------------------------------------------------------------- */
 namespace {
 template <UInt dim>
 class moment_helper {
 public:
   moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
 
   CUDA_LAMBDA Complex operator()(VectorProxy<Real, dim> q,
                                  const Complex& phi) const {
     Real mul = 1;
     for (UInt i = 0; i < dim; ++i)
       mul *= std::pow(q(i), exponent[i]);
 
     // Do not duplicate everything from hermitian symmetry
     if (std::abs(q.back()) < 1e-15)
       return mul * phi;
     return 2 * mul * phi;
   }
 
 private:
   std::array<UInt, dim> exponent;
 };
 }  // namespace
 
 template <>
 std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
   constexpr UInt dim = 1;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   return moments;
 }
 
 template <>
 std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
   constexpr UInt dim = 2;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
-  moments[1] = 0.5 * (m02 + m20);
+  moments[1] = (m02 + m20) / 2;
   auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
-  moments[2] = (3 * m22 + m40 + m04) / 3.;
+  moments[2] = (3 * m22 + m40 + m04) / 3;
   return moments;
 }
 
 template struct Statistics<1>;
 template struct Statistics<2>;
 
 }  // namespace tamaas
diff --git a/src/model/elastic_functional.cpp b/src/model/elastic_functional.cpp
index 7e4a3eb..12dcd34 100644
--- a/src/model/elastic_functional.cpp
+++ b/src/model/elastic_functional.cpp
@@ -1,71 +1,71 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 "elastic_functional.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 namespace functional {
 
 Real ElasticFunctionalPressure::computeF(GridBase<Real>& pressure,
                                          GridBase<Real>& dual) const {
   // This doesn't work if dual is multi-component
   *buffer = dual;
   // this can get large
   const unsigned long long n = dual.getGlobalNbPoints();
-  return (0.5 * pressure.dot(*buffer) - pressure.dot(surface)) / n / n;
+  return (pressure.dot(*buffer) / 2 - pressure.dot(surface)) / n / n;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void ElasticFunctionalPressure::computeGradF(GridBase<Real>& pressure,
                                              GridBase<Real>& gradient) const {
   this->op.apply(pressure, *buffer);
   // This doesn't work if dual is multi-component
   *buffer -= surface;
   gradient += *buffer;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real ElasticFunctionalGap::computeF(GridBase<Real>& gap,
                                     GridBase<Real>& dual) const {
   // This doesn't work if dual is multi-component
   *buffer = gap;
   *buffer += surface;
   const UInt n = dual.getGlobalNbPoints();
   const Real F = 0.5 * buffer->dot(dual);
   return F / n / n;  // avoid integer overflow
 }
 
 /* -------------------------------------------------------------------------- */
 
 void ElasticFunctionalGap::computeGradF(GridBase<Real>& gap,
                                         GridBase<Real>& gradient) const {
   *buffer = gap;
   *buffer += surface;  // Now we have displacements
   this->op.apply(*buffer, *buffer);
   gradient += *buffer;
 }
 
 }  // namespace functional
 
 }  // namespace tamaas
diff --git a/src/model/elasto_plastic/residual.cpp b/src/model/elasto_plastic/residual.cpp
index 4fab6ae..0d710e5 100644
--- a/src/model/elasto_plastic/residual.cpp
+++ b/src/model/elasto_plastic/residual.cpp
@@ -1,159 +1,159 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 <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}) {
+  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(strain->sizes(), strain->getNbComponents(),
                       strain_increment.getInternalData());
   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;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::applyTangent(
     GridBase<Real>& output, GridBase<Real>& input,
     GridBase<Real>& current_strain_increment) {
   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;
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::computeStress(GridBase<Real>& 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) {
   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) {
   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.globalDataSize() > 1e-14)
       plastic_layers.insert(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 void ResidualTemplate<type>::setIntegrationMethod(integration_method method,
                                                   Real cutoff) {
   auto& mindlin = dynamic_cast<Mindlin<type, 1>&>(
       *this->model->getIntegralOperator("mindlin"));
   auto& mindlin_g = dynamic_cast<Mindlin<type, 2>&>(
       *this->model->getIntegralOperator("mindlin_gradient"));
 
   mindlin.setIntegrationMethod(method, cutoff);
   mindlin_g.setIntegrationMethod(method, cutoff);
 }
 
 /* -------------------------------------------------------------------------- */
 template class ResidualTemplate<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/meta_functional.cpp b/src/model/meta_functional.cpp
index b2bfba7..69d2175 100644
--- a/src/model/meta_functional.cpp
+++ b/src/model/meta_functional.cpp
@@ -1,54 +1,54 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 "meta_functional.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 namespace functional {
 
 Real MetaFunctional::computeF(GridBase<Real>& variable,
                               GridBase<Real>& dual) const {
   Real F = 0;
-  for (auto& f : functionals)
+  for (const auto& f : functionals)
     F += f->computeF(variable, dual);
   return F;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void MetaFunctional::computeGradF(GridBase<Real>& variable,
                                   GridBase<Real>& gradient) const {
   gradient = 0;
-  for (auto& f : functionals)
+  for (const auto& f : functionals)
     f->computeGradF(variable, gradient);
 }
 
 /* -------------------------------------------------------------------------- */
 
 void MetaFunctional::addFunctionalTerm(std::shared_ptr<Functional> functional) {
   functionals.push_back(std::move(functional));
 }
 
 }  // namespace functional
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
diff --git a/src/model/model.cpp b/src/model/model.cpp
index d82d2b8..778c5e9 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,183 +1,183 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 {
 /* -------------------------------------------------------------------------- */
 void Model::setElasticity(Real E_, Real nu_) {
   setYoungModulus(E_);
   setPoissonRatio(nu_);
   updateOperators();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::applyElasticity(GridBase<Real>& stress,
                             const GridBase<Real>& strain) const {
   operators.at("hooke")->apply(const_cast<GridBase<Real>&>(strain), stress);
 }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getTraction() { return getField("traction"); }
 const GridBase<Real>& Model::getTraction() const { return (*this)["traction"]; }
 
 /* -------------------------------------------------------------------------- */
 GridBase<Real>& Model::getDisplacement() { return getField("displacement"); }
 const GridBase<Real>& Model::getDisplacement() const {
   return (*this)["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(getTraction(), getDisplacement());
 }
 
 void Model::solveDirichlet() {
   engine->registerDirichlet();
   engine->solveDirichlet(getDisplacement(), getTraction());
 }
 
 /* -------------------------------------------------------------------------- */
 IntegralOperator* Model::getIntegralOperator(const std::string& name) const {
   return operators.at(name).get();
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<std::string> Model::getIntegralOperators() const {
   std::vector<std::string> keys;
   keys.reserve(operators.size());
 
   std::transform(operators.begin(), operators.end(), std::back_inserter(keys),
                  [](auto&& pair) { return pair.first; });
   return keys;
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::updateOperators() {
   for (auto& op : operators)
     op.second->updateFromModel();
 }
 
 /* -------------------------------------------------------------------------- */
 void Model::addDumper(std::shared_ptr<ModelDumper> dumper) {
   this->dumpers.push_back(std::move(dumper));
 }
 
 void Model::dump() const {
-  for (auto& dumper : dumpers)
+  for (const 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;
 }
 
 const GridBase<Real>& Model::operator[](const std::string& name) const {
   return getField(name);
 }
 
 GridBase<Real>& Model::operator[](const std::string& name) {
   return getField(name);
 }
 
 /* -------------------------------------------------------------------------- */
 std::ostream& operator<<(std::ostream& o, const Model& _this) {
   o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
     << ", nu = " << _this.getPoissonRatio() << ")\n";
 
   auto out_collec = [&o](auto&& collec) {
     std::for_each(collec.begin(), collec.end() - 1,
                   [&o](const auto& x) { o << x << ", "; });
     o << collec.back();
   };
 
   // Printing domain size
   o << "  - domain = [";
   out_collec(_this.getSystemSize());
   o << "]\n";
 
   // Printing discretization
   o << "  - discretization = [";
   out_collec(_this.getDiscretization());
   o << "]\n";
 
   if (mpi::size() > 1) {
     o << "  - global discretization = [";
     out_collec(_this.getGlobalDiscretization());
     o << "]\n";
   }
 
   // Print fields
   o << "  - registered fields = [";
   out_collec(_this.getFields());
   o << "]\n";
 
   o << "  - registered operators = [";
   out_collec(_this.getIntegralOperators());
   o << "]";
 
   if (not _this.dumpers.empty())
     o << "\n  - " << _this.dumpers.size() << " registered dumpers";
 
   return o;
 }
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/solvers/dfsane_solver.hh b/src/solvers/dfsane_solver.hh
index 456b993..8bdd049 100644
--- a/src/solvers/dfsane_solver.hh
+++ b/src/solvers/dfsane_solver.hh
@@ -1,69 +1,69 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 DFSANE_SOLVER_HH
 #define DFSANE_SOLVER_HH
 /* -------------------------------------------------------------------------- */
 #include "ep_solver.hh"
 #include "grid_base.hh"
 #include "residual.hh"
 #include <deque>
 #include <functional>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /**
  @brief Derivative-free non-linear solver
 
  This algorithm is based on W. La Cruz, J. Martínez, and M. Raydan, “Spectral
  residual method without gradient information for solving large-scale nonlinear
  systems of equations,” Math. Comp., vol. 75, no. 255, pp. 1429–1448, 2006, doi:
  10.1090/S0025-5718-06-01840-0.
 
  The same algorithm is available in scipy.optimize, but this version is robustly
  parallel by default (i.e. does not depend on BLAS's parallelism and is
  future-proof for MPI parallelism).
 */
 class DFSANESolver : public EPSolver {
   // Public interface
 public:
   DFSANESolver(Residual& residual);
   void solve() override;
 
   // Algorithm functions
 protected:
   Real computeSpectralCoeff(const std::pair<Real, Real>& bounds);
   void computeSearchDirection(Real sigma);
   void lineSearch(Real eta_k);
-  Real computeAlpha(Real alpha, Real f, Real fk,
-                    const std::pair<Real, Real>& bounds);
+  static Real computeAlpha(Real alpha, Real f, Real fk,
+                           const std::pair<Real, Real>& bounds);
 
 protected:
   GridBase<Real> search_direction, previous_residual, current_x, delta_x,
       delta_residual;
   std::deque<Real> previous_merits;
   std::function<Real(UInt)> eta;
 };
 
 }  // namespace tamaas
 #endif
diff --git a/src/surface/isopowerlaw.cpp b/src/surface/isopowerlaw.cpp
index 0378f48..4d65109 100644
--- a/src/surface/isopowerlaw.cpp
+++ b/src/surface/isopowerlaw.cpp
@@ -1,100 +1,100 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 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 "isopowerlaw.hh"
 #include <map>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 template <UInt dim>
 void Isopowerlaw<dim>::computeFilter(
     GridHermitian<Real, dim>& filter_coefficients) const {
   Filter<dim>::computeFilter(
       [this] CUDA_LAMBDA(Complex & coeff, VectorProxy<Real, dim> q) {
         coeff = (*this)(q);
       },
       filter_coefficients);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::rmsHeights() const {
   return std::sqrt(moments()[0]);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Analytical moments, cf. Yastrebov et al. (2015)
  * "From infinitesimal to full contact between rough surfaces: Evolution
  * of the contact area", appendix A
  */
 template <>
 std::vector<Real> Isopowerlaw<2>::moments() const {
   std::map<UInt, Real> T;
   T[0] = 2 * M_PI;
   T[2] = M_PI;
-  T[4] = 3 * M_PI / 4.;
+  T[4] = (3 * M_PI) / 4;
   Real q1 = static_cast<Real>(this->q1);  // shadowing this->q1
   Real xi = q0 / q1;
   Real zeta = q2 / q1;
   Real A = std::pow(q1, 2 * (hurst + 1));
 
   std::vector<Real> moments;
   moments.reserve(3);
   for (UInt q : {0, 2, 4}) {
-    Real m = A * T[q] * std::pow(q1, q - 2. * hurst) *
-             ((1 - std::pow(xi, q + 2)) / (q + 2.) +
-              (std::pow(zeta, q - 2. * hurst) - 1) / (q - 2. * hurst));
+    Real m = A * T[q] * std::pow(q1, q - 2 * hurst) *
+             ((1 - std::pow(xi, q + 2)) / (q + 2) +
+              (std::pow(zeta, q - 2 * hurst) - 1) / (q - 2 * hurst));
     moments.push_back(m);
   }
 
   return moments;
 }
 
 template <>
 std::vector<Real> Isopowerlaw<1>::moments() const {
   TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::alpha() const {
   const auto m = moments();
   return m[0] * m[2] / (m[1] * m[1]);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Real Isopowerlaw<dim>::rmsSlopes() const {
   // 2 pi because moments are computed with k instead of q
   // and slopes should be computed with q
   return 2 * M_PI * std::sqrt(2 * moments()[1]);
 }
 
 template class Isopowerlaw<1>;
 template class Isopowerlaw<2>;
 
 }  // namespace tamaas