diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9d75005c..fb6dbd0f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,421 +1,422 @@
 # Changelog
 
 All notable changes to this project will be documented in this file.
 
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html) for final versions and [PEP440](https://www.python.org/dev/peps/pep-0440/) in case intermediate versions need to be released (e.g. development version `2.2.3.dev1` or release candidates `2.2.3rc1`), or individual commits are packaged.
 
 ## Unreleased
 
 ### Added
 
 - Added a `Material` class tree to decouple material behavior from residuals.
   The `Material` class can be derived in Python for custom material laws.
 - `Model` can set the integration method for volume operators with the
   `setIntegrationMethod` method.
 - `EPICSolver` can be used in `utils.load_path`
 - `residual.material` gives a reference to the material object used in the residual
 - Added a DC-FFT/Boussinesq operator for non-periodic elastic solutions
 - Can now choose integral operator in `PolonskyKeerRey`, so non-periodic contact
   is possible!
 - Python sub-classes of `IntegralOperator` are possible
 - `TAMAAS_MSG` variadic macro which concatenates arguments with file, line no
   and function info
 - `assertion_error`, `not_implemented_error`, `nan_error` and `model_type_error`
   classes for exceptions, which allow for finer exceptions for debugging
 - SCons targets `clang-format`, `clang-tidy`, `flake8`, `mypy` for linting, and
   target alias `lint` to run all available lint tools
 
 ### Changed
 
 - Removed `model_type` template for `Residual`: they are only defined on
   the `model_type::volume_2d` type.
 - Simplified interface for `Residual`: it now takes a `Material` object at
   construction which handles the computation of eigenstresses involved in the
   residual equation
 - Added function name to `TAMAAS_MSG`
 - `TAMAAS_ASSERT` raises `assertion_error` and is not disabled in release
   builds: it should be used to check class invariants
 - `TAMAAS_EXCEPTION` calls have been replaced by either `throw <exception>` or
   `TAMAAS_ASSERT` calls, so that specific more error types can be used instead
   of the catch-all `tamaas::Exception` type
   
 ### Fixed
 
 - Random phases set to zero where the C2R transform expects real inputs
+- Issue where large HDF5 files could not be read in parallel
 
 ### Deprecated
 
 - `getStress()`, `setIntegrationMethod()` and `getVector()` in `Residual` are
   deprecated, use `model["stress"]`, `model.setIntegrationMethod()` and
   `residual.vector` respectively instead.
 - `ModelFactor.createResidual` is deprecated. Instead, instanciate a `Material`
   object and use the `Residual` constructor to create the residual object.
   
 ### Removed
 
 - `tamaas::Exception` class
 - `TAMAAS_EXCEPTION`, `TAMAAS_DEBUG_MSG`, `TAMAAS_DEBUG_EXCEPTION` macros
 
 ## v2.6.0 -- 2022-12-23
 
 ### Added
 
 - Added `matvec` interface for the following operators:
   - `Hooke`
   - `Kelvin`
   - `Westergaard`
 - Documented error handling with MPI
 - Added `span<T>` class
 - `FieldContainer` class that defines an interface to access fields registered
   with an object
 - Operators now expose their fields, which can be modified from Python
 - Added `HookeField` operator for heterogeneous elastic properties
 - Documented proper unit choice
 - Added `Statistics::computeFDRMSSlope()`, finite difference approximation of
   RMS of surface slopes
 - Kato saturated uses L2 norm on residual displacement instead of L∞
 - Request system to avoid erasing model data on restart
 - Boundary field names are dumped
 
 ### Changed
 
 - Resizing a grid no longer sets values to zero
 - Wheels built for `manylinux2014` instead of `manylinux2010`
 - Cleaned up `Array` interface
 - Shared-memory compile options disabled by default
 - Unknown build variables raise a warning
 - Removed conditional branches from PKR
 - NetCDF dumps in `NETCDF4` format instead of `NETCDF4_CLASSIC`
 - VTK no longer writes FieldData, this was causing issues reading converted files
 
 ### Fixed
 
 - Fixed flood fill algorithm to properly account for PBCs
 
 ## v2.5.0 -- 2022-07-05
 
 ### Added
 
 - Added `tamaas.utils.publications` to print a list of publications relevant to
   the parts of Tamaas used in a script. Call at the end of a script for an
   exhaustive list.
 - Added `tamaas.mpi.gather` and `tamaas.mpi.scatter` to help handling of 2D data
   in MPI contexts.
 - Added `getBoundaryFields()/boundary_fields` function/property to `Model`.
 - Added Python constructor to `Model`, which is an alias to
   `ModelFactory.createModel`.
 - Added convenience `tamaas.dumpers._helper.netCDFtoParaview` function to
   convert model sequence to Paraview data file.
 - Added dump of surface fields as `FieldData` in VTK files
 - Added `scipy.sparse.linalg.LinearOperator` interface to `Westergaard`. A
   `LinearOperator` can be created by calling
   `scipy.sparse.linalg.aslinearoperator` on the operator object.
 - Allowed sub-classing of `ContactSolver` in Python.
 - Exposed elastic contact functionals in Python to help with custom solver
   implementation.
 - Added an example of custom solver with penalty method using Scipy.
 - Added `graphArea` function in statistics to compute area of function graph
 - Added `tamaas.utils.radial_average` to radialy average a 2D field
 - Added option to not make a dump directory (e.g. `netcdf`, `hdf5` directories
   for dump files). This allows to specify an absolute path for dump files, and
   is useful if simulation workflow is managed with file targets, as do systems
   like GNU Make and Snakemake.
 
 ### Changed
 
 - Filter for boundary fields in dumper should be smarter.
 - `read` functions in dumpers are now `@classmethod`
 - Changed sign of plastic correction in KatoSaturated solver to match the
   EPICSolver sign.
 - Plastic correction of KatoSaturated is registered in the model as
   `KatoSaturated::residual_displacement`.
 - Changed to modern Python packaging (PEP517) with `setup.cfg` and
   `pyproject.toml`. A `setup.py` script is still provided for editable
   installations, see [gh#7953](https://github.com/pypa/pip/issues/7953). Once
   setuptools has stable support for
   [PEP621](https://peps.python.org/pep-0621/), `setup.cfg` will be merged into
   `pyproject.toml`. Note that source distributions that can be built with
   [`build`](https://pypi.org/project/build/) are not true source distributions,
   since they contain the compiled Python module of Tamaas.
 - Updated copyright to reflect Tamaas' development history since leaving EPFL
 
 ### Fixed
 
 - Fixed `tamaas.dumpers.NetCDFDumper` to dump in MPI context.
 
 ### Deprecated
 
 - SCons versions < 3 and Python < 3.6 are explicitely no longer supported.
 
 ## v2.4.0 -- 2022-03-22
 
 ### Added
 
 - Added a `tamaas.utils` module.
 - Added `tamaas.utils.load_path` generator, which yields model objects for a
   sequence of applied loads.
 - Added `tamaas.utils.seeded_surfaces` generator, which yields surfaces for a
   sequence of random seeds.
 - Added `tamaas.utils.hertz_surface` which generates a parabolic (Hertzian)
   surface.
 - Added `tamaas.utils.log_context` context manager which locally sets the log
   level for Tamaas' logger.
 - Added `tamaas.compute.from_voigt` to convert Voigt/Mendel fields to dense
   tensor fields.
 - Added a deep-copy function to `ModelFactory` to copy `Model` objects. Use
   `model_copy = copy.deepcopy(model)` in Python to make a copy. Currently only
   copies registered fields, same as dumpers/readers.
 - Added a `read_sequence` method for dumpers to read model frames.
 - Automatic draft release on Zenodo.
 
 ### Changed
 
 - `*_ROOT` variables for vendored dependencies GoogleTest and pybind11 now
   default to empty strings, so that the vendored trees in `third-party/` are not
   selected by default. This is so system packages are given priority. Vendored
   dependency submodules will eventually be deprecated.
 
 ### Fixed
 
 - Fixed an issue with `scons -c` when GTest was missing.
 
 ## v2.3.1 -- 2021-11-08
 
 ### Added
 
 - Now using `clang-format`, `clang-tidy` and `flake8` for linting
 - Pre-built Tamaas images can be pulled with `docker pull
   registry.gitlab.com/tamaas/tamaas`
 - Added a `--version` flag to the `tamaas` command
 
 ### Changed
 
 - The root `Dockerfile` now compiles Tamaas, so using Tamaas in Docker is easier.
 - The model attributes dumped to Numpy files are now written in a JSON-formatted
   string to avoid unsafe loading/unpickling of objects.
 - Removed the `build_doc` build option: now the doc targets are automatically
   added if the dependencies are met, and built if `scons doc` is called.
 - Removed the `use_googletest` build option: if tests are built and gtest is
   present, the corresponding tests will be built
 
 ### Fixed
 
 - The command `tamaas plot` gracefully fails if Matplotlib is missing
 - Better heuristic to guess which fields are defined on boundary in dumpers
   (still not perfect and may give false negatives)
 
 ## v2.3.0 -- 2021-06-15
 
 ### Added
 
 - Added `read()` method to dumpers to create a model from a dump file
 - `getClusters()` can be called in MPI contact with partial contact maps
 - Added a JSON encoder class for models and a JSON dumper
 - CUDA compatibility is re-established, but has not been tested
 - Docstrings in the Python bindings for many classes/methods
 
 ### Changed
 
 - Tamaas version numbers are now managed by
   [versioneer](https://github.com/python-versioneer/python-versioneer). This
   means that Git tags prefixed with `v` (e.g. `v2.2.3`) carry meaning and
   determine the version. When no tag is set, versioneer uses the last tag,
   specifies the commit short hash and the distance to the last tag (e.g.
   `2.2.2+33.ge314b0e`). This version string is used in the compiled library, the
   `setup.py` script and the `__version__` variable in the python module.
 - Tamaas migrated to [GitLab](https://gitlab.com/tamaas/tamaas)
 - Continuous delivery has been implemented:
 
     - the `master` branch will now automatically build and publish Python wheels
       to `https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple`. These
       "nightly" builds can be installed with:
 
           pip install \
             --extra-index-url https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple \
             tamaas
 
     - version tags pushed to `master` will automatically publish the wheels to
       [PyPI](https://pypi.org/project/tamaas/)
 
 ### Deprecated
 
 - The `finalize()` function is now deprecated, since it is automatically called
   when the process terminates
 - Python versions 3.5 and below are not supported anymore
 
 ### Fixed
 
 - Fixed a host of dump read/write issues when model type was not `volume_*d`.
   Dumper tests are now streamlined and systematic.
 - Fixed a bug where `Model::solveDirichlet` would not compute correctly
 - Fixed a bug where `Statistics::contact` would not normalize by the global
   number of surface points
 
 ## v2.2.2 -- 2021-04-02
 
 ### Added
 
 - Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try
   executing `tamaas surface -h` from the command-line!
 
 ### Changed
 
 - `CXXFLAGS` are now passed to the linker
 - Added this changelog
 - Using absolute paths for environmental variables when running `scons test`
 - Reorganized documentation layout
 - Gave the build system a facelift (docs are now generated directly with SCons
   instead of a Makefile)
 
 ### Deprecated
 
 - Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with
   a Python 2 wheel.
 - The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command
 
 ### Fixed
 
 - `UVWDumper` no longer imports `mpi4py` in sequential
 - Compiling with different Thrust/FFTW backends
 
 
 ## v2.2.1 -- 2021-03-02
 
 ### Added
 
 - Output registered fields and dumpers in `print(model)`
 - Added `operator[]` to the C++ model class (for fields)
 - Added `traction` and `displacement` properties to Python model bindings
 - Added `operators` property to Python model bindings, which provides a
   dict-like access to registered operators
 - Added `shape` and `spectrum` to properties to Python surface generator
   bindings
 - Surface generator constructor accepts surface global shape as argument
 - Choice of FFTW thread model
 
 ### Changed
 
 - Tests use `/tmp` for temporary files
 - Updated dependency versions (Thrust, Pybind11)
 
 ### Deprecated
 
 - Most `get___()` and `set___()` in Python bindings have been deprecated. They
   will generate a `DeprecationWarning`.
 
 ### Removed
 
 - All legacy code
 
 
 ## v2.2.0 -- 2020-12-31
 
 ### Added
 
 - More accurate function for computation of contact area
 - Function to compute deviatoric of tensor fields
 - MPI implementation
 - Convenience `hdf5toVTK` function
 - Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give
   shape information
 
 ### Changed
 
 - Preprocessor defined macros are prefixed with `TAMAAS_`
 - Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt`
 
 ### Fixed
 
 - Warning about deprecated constructors with recent GCC versions
 - Wrong computation of grid strides
 - Wrong computation of grid sizes in views
 
 
 ## v2.1.4 -- 2020-08-07
 
 ### Added
 
 - Possibility to generate a static `libTamaas`
 - C++ implementation of DFSANE solver
 - Allowing compilation without OpenMP
 
 ### Changed
 
 - NetCDF dumper writes frames to a single file
 
 ### Fixed
 
 - Compatibility with SCons+Python 3
 
 ## v2.1.3 -- 2020-07-27
 
 ### Added
 
 - Version number to `TamaasInfo`
 
 ### Changed
 
 - Prepending root directory when generating archive
 
 
 ## v2.1.2 -- 2020-07-24
 
 This release changes some core internals related to discrete Fourier transforms
 for future MPI support.
 
 ### Added
 
 - Caching `CXXFLAGS` in SCons build
 - SCons shortcut to create code archive
 - Test of the elastic-plastic contact solver
 - Paraview data dumper (`.pvd` files)
 - Compression for UVW dumper
 - `__contains__` and `__iter__` Python bindings of model
 - Warning message of possible overflow in Kelvin
 
 ### Changed
 
 - Simplified `tamaas_info.cpp`, particularly the diff part
 - Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are
   re-used as much as possible with different data with the same shape. This is
   in view of future MPI developments
 - Redirecting I/O streams in solve functions so they can be used from Python
   (e.g. in Jupyter notebooks)
 - Calling `initialize()` and `finalize()` is no longer necessary
 
 ### Fixed
 
 - Convergence issue with non-linear solvers
 - Memory error in volume potentials
 
 
 ## v2.1.1 -- 2020-04-22
 
 ### Added
 
 - SCons shortcut to run tests
 
 ### Fixed
 
 - Correct `RPATH` for shared libraries
 - Issues with SCons commands introduced in v2.1.0
 - Tests with Python 2.7
 
 
 ## v2.1.0 -- 2020-04-17
 
 ### Added
 
 - SCons shortcuts to build/install Tamaas and its components
 - Selection of integration method for Kelvin operator
 - Compilation option to remove the legacy part of Tamaas
 - NetCDF dumper
 
 ### Fixed
 
 - Link bug with clang
 - NaNs in Kato saturated solver
 
 ## v2.0.0 -- 2019-11-11
 
 First public release. Contains relatively mature elastic-plastic contact code.
diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index bbac664c..03f94165 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,159 +1,159 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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 "kelvin.hh"
 #include "logger.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* Constructors */
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 Kelvin<type, derivative>::Kelvin(Model* model) : VolumePotential<type>(model) {
   setIntegrationMethod(integration_method::linear, 1e-12);
 }
 
 template <model_type type, UInt derivative>
 void Kelvin<type, derivative>::setIntegrationMethod(integration_method method,
                                                     Real cutoff) {
   this->method = method;
   this->cutoff = cutoff;
 
   Logger logger;
 
   if (this->method == integration_method::linear) {
     logger.get(LogLevel::debug)
         << TAMAAS_MSG("Setting linear integration method");
     this->initialize(dtrait::template source_components<type>,
                      dtrait::template out_components<type>,
                      this->model->getDiscretization()[0]);
   }
 
   else {
     logger.get(LogLevel::debug) << TAMAAS_MSG(
         "Setting cutoff integration method (cutoff ", this->cutoff, ')');
     this->initialize(dtrait::template source_components<type>,
                      dtrait::template out_components<type>, 1);
   }
 
   auto max_q = Loop::reduce<operation::max>(
       [] CUDA_LAMBDA(VectorProxy<const Real, trait::boundary_dimension> qv) {
         return qv.l2norm();
       },
       range<VectorProxy<const Real, trait::boundary_dimension>>(
           this->wavevectors));
 
   if (this->method == integration_method::linear and
       not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
     logger.get(LogLevel::warning)
         << "Probable overflow of integral computation (consider "
            "changing integration method to integration_method::cutoff or "
            "compiling with real_type='long double')\n";
 }
 
 /* -------------------------------------------------------------------------- */
 /* Operator implementation */
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 void Kelvin<type, derivative>::applyIf(GridBase<Real>& source,
                                        GridBase<Real>& out,
                                        filter_t pred) const {
   KelvinInfluence kelvin(this->model->getShearModulus(),
                          this->model->getPoissonRatio());
 
   parent::transformSource(source, pred);
 
   // Reset buffer values
   for (auto&& layer : this->out_buffer)
     layer = 0;
 
   if (method == integration_method::linear) {
     linearIntegral(out, kelvin);
   } else {
     cutoffIntegral(out, kelvin);
   }
 }
 
 template <model_type type, UInt derivative>
 void Kelvin<type, derivative>::linearIntegral(GridBase<Real>& out,
                                               KelvinInfluence& kelvin) const {
   detail::KelvinHelper<type, KelvinInfluence> helper;
 
   helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
                        this->model->getSystemSize().front(), kelvin);
   helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
   helper.makeFundamentalGreatAgain(this->out_buffer);
 
   parent::transformOutput(
       [](auto&& out_buffer, auto layer) ->
       typename parent::BufferType& { return out_buffer[layer]; },
       out);
 }
 
 template <model_type type, UInt derivative>
 void Kelvin<type, derivative>::cutoffIntegral(GridBase<Real>& out,
                                               KelvinInfluence& kelvin) const {
   detail::KelvinHelper<type, KelvinInfluence> helper;
 
   auto func = [&](auto&& out_buffer,
                   auto layer) -> typename parent::BufferType& {
     auto&& out = out_buffer.front();
     out = 0;
     helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors,
                          this->model->getSystemSize().front(), cutoff, kelvin);
     helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
     helper.makeFundamentalGreatAgain(out);
     return out;
   };
   parent::transformOutput(func, out);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 std::pair<UInt, UInt> Kelvin<type, derivative>::matvecShape() const {
   const auto N = this->model->getDisplacement().getGlobalNbPoints();
   return std::make_pair(dtrait::template out_components<type> * N,
                         dtrait::template source_components<type> * N);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 GridBase<Real> Kelvin<type, derivative>::matvec(GridBase<Real>& X) const {
   constexpr auto dim = trait::dimension;
   const auto shape = this->model->getDiscretization();
   Grid<Real, dim> x{shape, dtrait::template source_components<type>, X.view()};
   Grid<Real, dim> y{shape, dtrait::template out_components<type>};
   this->apply(x, y);
-  return std::move(y);
+  return y;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation */
 /* -------------------------------------------------------------------------- */
 template class Kelvin<model_type::volume_2d, 0>;
 template class Kelvin<model_type::volume_2d, 1>;
 template class Kelvin<model_type::volume_2d, 2>;
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 02678659..24f68c9e 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,309 +1,309 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2023 Lucas Frérot
  *
  *  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_hermitian.hh"
 #include "grid_view.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "model_type.hh"
 #include "static_types.hh"
 
 #include "westergaard.hh"
 
 #include <functional>
 #include <numeric>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <model_type mtype, IntegralOperator::kind otype>
 Westergaard<mtype, otype>::Westergaard(Model* model)
     : IntegralOperator(model), buffer(), engine(FFTEngine::makeEngine()) {
   // 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 = allocateGrid<mtype, true, Real, GridHermitian>(
       boundary_hermitian_sizes, nb_components * nb_components);
   (*this)["influence"] = influence;
   initInfluence();
 }
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <UInt m, UInt j>
 struct boundary_fft_helper {
   template <typename Buffer, typename Out>
   static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
     auto view = make_view(out, 0);
     e.backward(view, buffer);
   }
 };
 
 template <UInt m>
 struct boundary_fft_helper<m, m> {
   template <typename Buffer, typename Out>
   static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
     e.backward(out, buffer);
   }
 };
 }  // 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 {
   auto& i = dynamic_cast<Grid<Real, bdim>&>(in);
   auto& full_out = dynamic_cast<Grid<Real, dim>&>(out);
 
   engine->forward(i, buffer);
 
   /// Applying influence
   func(buffer, *influence);
 
   /// Backward fourier transform on boundary only
   detail::boundary_fft_helper<bdim, dim>::backwardTransform(*engine, buffer,
                                                             full_out);
 
 } catch (const std::bad_cast& e) {
   throw model_type_error{
       TAMAAS_MSG("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 = FFTEngine::template computeFrequencies<Real, bdim, 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));
 
   if (mpi::rank() == 0) {
     // Set fundamental mode to zero
     MatrixProxy<Complex, comp, comp> mat((*influence)(0));
     mat = 0;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 /// \cond DO_NOT_DOCUMENT
 #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() {}
 /// \endcond
 
 /* ------------------------------------------------------------------------ */
 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 <model_type type, IntegralOperator::kind otype>
 std::pair<UInt, UInt> Westergaard<type, otype>::matvecShape() const {
   auto n = model->getBoundaryDiscretization();
   auto N = std::accumulate(n.begin(), n.end(), comp, std::multiplies<>());
   return std::make_pair(N, N);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, IntegralOperator::kind otype>
 GridBase<Real> Westergaard<type, otype>::matvec(GridBase<Real>& X) const {
   auto n = model->getBoundaryDiscretization();
   // Creating correct input
   Grid<Real, bdim> x{n, comp, X.view()};
 
   // Creating correct output
   if (bdim != dim)
     n.insert(n.begin(), 1);
 
   Grid<Real, dim> y(n, comp);
   apply(x, y);
-  return std::move(y);
+  return y;
 }
 
 /* -------------------------------------------------------------------------- */
 /* 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
 /* -------------------------------------------------------------------------- */