diff --git a/python/tamaas/nonlinear_solvers/__init__.py b/python/tamaas/nonlinear_solvers/__init__.py
index b10ba88..d30b39a 100644
--- a/python/tamaas/nonlinear_solvers/__init__.py
+++ b/python/tamaas/nonlinear_solvers/__init__.py
@@ -1,167 +1,171 @@
 # -*- mode:python; coding: utf-8 -*-
 #
 # Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 # Copyright (©) 2020-2022 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/>.
 
 """Nonlinear solvers for plasticity problems.
 
 Solvers in this module use :py:mod:`scipy.optimize` to solve the implicit
 non-linear equation for plastic deformations with fixed contact pressures.
 """
 from functools import wraps
 
 from scipy.optimize import newton_krylov, root
 from scipy.optimize.nonlin import NoConvergence
 
 from .. import EPSolver, Logger, LogLevel, mpi
 from .._tamaas import _tolerance_manager
 from .._tamaas import _DFSANESolver as DFSANECXXSolver
 
 
 __all__ = ['NLNoConvergence',
            'DFSANESolver',
            'DFSANECXXSolver',
            'NewtonKrylovSolver',
            'ToleranceManager']
 
 
 class NLNoConvergence(Exception):
     """Convergence not reached exception."""
 
 
 class ScipySolver(EPSolver):
     """Base class for solvers wrapping SciPy routines."""
 
     def __init__(self, residual, model=None, callback=None):
         """Construct nonlinear solver with residual.
 
         :param residual: plasticity residual object
         :param model: Deprecated
         :param callback: passed on to the Scipy solver
         """
         super(ScipySolver, self).__init__(residual)
 
         if mpi.size() > 1:
             raise RuntimeError("Scipy solvers cannot be used with MPI; "
                                "DFSANECXXSolver can be used instead")
 
         self.callback = callback
         self._x = self.getStrainIncrement()
         self._residual = self.getResidual()
         self.options = {'ftol': 0, 'fatol': 1e-9}
 
     def solve(self):
         """Solve R(δε) = 0 using a Scipy function."""
         # For initial guess, compute the strain due to boundary tractions
         # self._residual.computeResidual(self._x)
         # self._x[...] = self._residual.getVector()
         EPSolver.beforeSolve(self)
 
         # Scipy root callback
         def compute_residual(vec):
             self._residual.computeResidual(vec)
             return self._residual.getVector().copy()
 
         # Solve
+        Logger().get(LogLevel.debug) << \
+            "Entering non-linear solve\n"
         self._x[...] = self._scipy_solve(compute_residual)
+        Logger().get(LogLevel.debug) << \
+            "Non-linear solve returned"
 
         # Computing displacements
         self._residual.computeResidualDisplacement(self._x)
 
     def reset(self):
         """Set solution vector to zero."""
         self._x[...] = 0
 
     def _scipy_solve(self, compute_residual):
         """Actually call the scipy solver.
 
         :param compute_residual: function returning residual for given variable
         """
         raise NotImplementedError()
 
 
 class NewtonKrylovSolver(ScipySolver):
     """Solve using a finite-difference Newton-Krylov method."""
 
     def _scipy_solve(self, compute_residual):
         try:
             return newton_krylov(compute_residual, self._x,
                                  f_tol=self.tolerance,
                                  verbose=True, callback=self.callback)
         except NoConvergence:
             raise NLNoConvergence("Newton-Krylov did not converge")
 
 
 class DFSANESolver(ScipySolver):
     """Solve using a spectral residual jacobianless method.
 
     See :doi:`10.1090/S0025-5718-06-01840-0` for details on method and
     the relevant Scipy `documentation
     <https://scipy.github.io/devdocs/reference/optimize.root-dfsane.html>`_ for
     details on parameters.
 
     """
 
     def _scipy_solve(self, compute_residual):
         solution = root(compute_residual,
                         self._x,
                         method='df-sane',
                         options={'ftol': 0, 'fatol': self.tolerance},
                         callback=self.callback)
         Logger().get(LogLevel.info) << \
             "DF-SANE/Scipy: {} ({} iterations, {})".format(
                 solution.message,
                 solution.nit,
                 self.tolerance)
 
         if not solution.success:
             raise NLNoConvergence("DF-SANE/Scipy did not converge")
         return solution.x.copy()
 
 
 def ToleranceManager(start, end, rate):
     """Decorate solver to manage tolerance of non-linear solve step."""
 
     def actual_decorator(cls):
         orig_init = cls.__init__
         orig_solve = cls.solve
         orig_update_state = cls.updateState
 
         @wraps(cls.__init__)
         def __init__(obj, *args, **kwargs):
             orig_init(obj, *args, **kwargs)
             obj.setToleranceManager(_tolerance_manager(start, end, rate))
 
         @wraps(cls.solve)
         def new_solve(obj, *args, **kwargs):
             ftol = obj.tolerance
             ftol *= rate
 
             obj.tolerance = max(ftol, end)
             return orig_solve(obj, *args, **kwargs)
 
         @wraps(cls.updateState)
         def updateState(obj, *args, **kwargs):
             obj.tolerance = start
             return orig_update_state(obj, *args, **kwargs)
 
         cls.__init__ = __init__
         # cls.solve = new_solve
         # cls.updateState = updateState
         return cls
 
     return actual_decorator
diff --git a/src/core/array.hh b/src/core/array.hh
index 5bc3124..c145d01 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,168 +1,179 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef ARRAY_HH
 #define ARRAY_HH
 /* -------------------------------------------------------------------------- */
 #include "logger.hh"
+#include "span.hh"
 #include "tamaas.hh"
 #include <memory>
 #include <thrust/copy.h>
 #include <thrust/fill.h>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Generic storage class with wrapping capacities
 template <typename T>
 struct Array final {
   /// Default
   Array() = default;
+
   /// Empty array of given size
-  Array(UInt size) : Array() { this->resize(size); }
+  Array(UInt size) : Array() { resize(size); }
+
   /// Copy constructor (deep)
   Array(const Array& v) : Array() {
-    this->resize(v._size);
-    thrust::copy_n(v._data, v._size, this->_data);
+    resize(v.size());
+    thrust::copy(v.view_.begin(), v.view_.end(), view_.begin());
   }
+
   /// 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);
+    view_ = std::exchange(v.view_, span<T>{});
+    reserved_ = std::exchange(v.reserved_, 0);
+    wrapped_ = std::exchange(v.wrapped_, false);
   }
+
   /// Wrap array on data
-  Array(T* data, UInt size) : Array() { this->wrap(data, size); }
+  Array(T* data, UInt size) noexcept : Array() { wrap(data, size); }
+
+  /// Wrap on span
+  Array(span<T> view) noexcept : Array() { wrap(view); }
+
   /// Destructor
   ~Array() {
-    if (not wrapped)
-      _alloc.deallocate(_data, _reserved);
+    if (not wrapped_)
+      alloc_.deallocate(view_);
   }
 
   /// 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);
+    wrapped_ = false;
+    resize(v.size());
+    thrust::copy(v.view_.begin(), v.view_.end(), view_.begin());
     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);
-    }
+    if (this == &v)
+      return *this;
+
+    if (not wrapped_)
+      alloc_.deallocate(view_);
+    view_ = std::exchange(v.view_, span<T>{});
+    reserved_ = std::exchange(v.reserved_, 0);
+    wrapped_ = std::exchange(v.wrapped_, false);
     return *this;
   }
 
-  void wrap(const Array& other) { this->wrap(other._data, other._size); }
+  /// Wrap on view
+  Array& operator=(span<T> v) noexcept { wrap(v); }
 
-  /// Wrap a memory pointer
-  void wrap(T* data, UInt size) {
-    _data = data;
-    _size = size;
-    wrapped = true;
-    _reserved = 0;
+  /// Wrap array
+  void wrap(const Array& other) noexcept { wrap(other.view_); }
+
+  /// Wrap view
+  void wrap(span<T> view) noexcept {
+    view_ = view;
+    wrapped_ = true;
+    reserved_ = 0;
   }
 
+  /// Wrap a memory pointer
+  void wrap(T* data, UInt size) noexcept { wrap(span<T>{data, size}); }
+
   /// Data pointer access (const)
-  const T* data() const { return _data; }
+  const T* data() const { return view_.data(); }
+
   /// Data pointer access (non-const)
-  T* data() { return _data; }
+  T* data() { return view_.data(); }
 
   /// Resize array
-  void resize(UInt size, const T& value = T()) {
-    if (wrapped)
+  void resize(UInt new_size, const T& value = T()) {
+    if (wrapped_)
       TAMAAS_EXCEPTION("cannot resize wrapped array");
 
     // Erase array
-    if (size == 0) {
-      _alloc.deallocate(_data, _reserved);
-      _data = nullptr;
-      _size = 0;
-      _reserved = 0;
+    if (new_size == 0) {
+      alloc_.deallocate(view_);
+      view_ = span<T>{};
+      reserved_ = 0;
       return;
     }
 
     // Do nothing
-    if (size == this->_size)
+    if (new_size == size())
       return;
 
     // Allocate new data
-    _alloc.deallocate(_data, _reserved);
+    alloc_.deallocate(view_);
 
-    _data = _alloc.allocate(size);
-    _size = size;
-    _reserved = size;
+    view_ = alloc_.allocate(new_size);
+    reserved_ = new_size;
 
-    if (not wrapped)
-      thrust::fill(_data, _data + size, value);
+    if (not wrapped_)
+      thrust::fill(view_.begin(), view_.end(), value);
   }
 
   /// Reserve storage space
   void reserve(UInt size) {
-    if (_reserved >= 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;
+    auto new_view = alloc_.allocate(size);
+
+    if (new_view.data() != view_.data()) {
+      thrust::copy(view_.begin(), view_.end(), new_view.begin());
+
+      alloc_.deallocate(view_);
+      view_ = {new_view.data(), view_.size()};
+      reserved_ = view_.size();
+    } else {
+      reserved_ = size;
     }
   }
 
   /// Access operator
-  inline T& operator[](UInt i) {
-    TAMAAS_ASSERT(i < _size, "memory overflow");
-    return _data[i];
-  }
+  inline T& operator[](UInt i) { return view_[i]; }
 
   /// Access operator (const)
-  inline const T& operator[](UInt i) const {
-    TAMAAS_ASSERT(i < _size, "memory overflow");
-    return _data[i];
-  }
+  inline const T& operator[](UInt i) const { return view_[i]; }
 
   /// Get size of array
-  inline UInt size() const { return _size; }
+  inline UInt size() const { return view_.size(); }
+
+  span<T> view() const { return view_; }
 
 private:
-  T* _data = nullptr;
-  UInt _size = 0;
-  std::size_t _reserved = 0;
-  bool wrapped = false;
-  Allocator<T> _alloc;
+  span<T> view_;
+  typename span<T>::size_type reserved_ = 0;
+  bool wrapped_ = false;
+  Allocator<T> alloc_;
 };
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif /* ARRAY_HH */
diff --git a/src/core/cuda/unified_allocator.hh b/src/core/cuda/unified_allocator.hh
index 7630e9d..97f4983 100644
--- a/src/core/cuda/unified_allocator.hh
+++ b/src/core/cuda/unified_allocator.hh
@@ -1,51 +1,52 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef UNIFIED_ALLOCATOR_HH
 #define UNIFIED_ALLOCATOR_HH
 /* -------------------------------------------------------------------------- */
+#include "span.hh"
 #include <cuda_runtime_api.h>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Class allocating [unified
 /// memory](http://docs.nvidia.com/cuda/cuda-runtime-api/group__CUDART__MEMORY.html#group__CUDART__MEMORY_1gd228014f19cc0975ebe3e0dd2af6dd1b)
 template <typename T>
 struct UnifiedAllocator {
   /// Allocate memory
-  static T* allocate(std::size_t n) {
+  static span<T> allocate(typename span<T>::size_type n) noexcept {
     T* p = nullptr;
     cudaMallocManaged(&p, n * sizeof(T));
     return p;
   }
 
   /// Free memory
   static void deallocate(T* p, __attribute__((unused)) std::size_t n) {
     cudaFree(p);
   }
 };
 
 }  // namespace tamaas
 
 #endif  // UNIFIED_ALLOCATOR_HH
diff --git a/src/core/fftw/fftw_allocator.hh b/src/core/fftw/fftw_allocator.hh
index f0eab15..5d6f793 100644
--- a/src/core/fftw/fftw_allocator.hh
+++ b/src/core/fftw/fftw_allocator.hh
@@ -1,52 +1,52 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FFTW_ALLOCATOR_HH
 #define FFTW_ALLOCATOR_HH
 /* -------------------------------------------------------------------------- */
+#include "span.hh"
 #include <fftw3.h>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Class allocating
 /// [SIMD](http://www.fftw.org/fftw3_doc/SIMD-alignment-and-fftw_005fmalloc.html#SIMD-alignment-and-fftw_005fmalloc)
 /// aligned memory
 template <typename T>
 struct FFTWAllocator {
   /// Allocate memory
-  static T* allocate(std::size_t n) noexcept {
+  static span<T> allocate(typename span<T>::size_type n) noexcept {
     T* p = nullptr;
-    p = (T*)fftw_malloc(sizeof(T) * n);
-    return p;
+    p = static_cast<T*>(fftw_malloc(sizeof(T) * n));
+    span<T> view{p, n};
+    return view;
   }
 
   /// Free memory
-  static void deallocate(T* p, __attribute__((unused)) std::size_t n) noexcept {
-    fftw_free(p);
-  }
+  static void deallocate(span<T> view) noexcept { fftw_free(view.data()); }
 };
 
 }  // namespace tamaas
 
 #endif  // FFTW_ALLOCATOR_HH
diff --git a/src/core/grid.cpp b/src/core/grid.cpp
index 55f8204..8086e65 100644
--- a/src/core/grid.cpp
+++ b/src/core/grid.cpp
@@ -1,122 +1,100 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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.hh"
 #include "tamaas.hh"
 #include <algorithm>
 #include <complex>
 #include <cstring>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid() : GridBase<T>() {
   this->n.fill(0);
   this->strides.fill(1);
   this->nb_components = 1;
 }
 
 /* -------------------------------------------------------------------------- */
 
-template <typename T, UInt dim>
-Grid<T, dim>::Grid(const std::array<UInt, dim>& n, UInt nb_components,
-                   value_type* data)
-    : GridBase<T>() {
-  this->n = n;
-  this->nb_components = nb_components;
-  UInt size = this->computeSize();
-  this->data.wrap(data, size);
-  this->computeStrides();
-}
-
-template <typename T, UInt dim>
-Grid<T, dim>::Grid(const Grid<T, dim>& o)
-    : GridBase<T>(o), n(o.n), strides(o.strides) {}
-
-template <typename T, UInt dim>
-Grid<T, dim>::Grid(Grid<T, dim>&& o) noexcept
-    : GridBase<T>(std::move(o)), n(std::move(o.n)),
-      strides(std::move(o.strides)) {}
-
-/* -------------------------------------------------------------------------- */
-
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const std::array<UInt, dim>& n) {
   this->resize(n.begin(), n.end());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const std::vector<UInt>& n) {
   TAMAAS_ASSERT(n.size() == dim, "Shape vector not matching grid dimensions");
   this->resize(n.begin(), n.end());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(std::initializer_list<UInt> n) {
   TAMAAS_ASSERT(n.size() == dim,
                 "Shape initializer list not matching grid dimensions");
   this->resize(std::begin(n), std::end(n));
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::computeStrides() {
   std::copy(n.begin() + 1, n.end(), strides.begin());
   strides[dim] = 1;
   strides[dim - 1] = this->nb_components;
   std::partial_sum(strides.rbegin(), strides.rend(), strides.rbegin(),
                    std::multiplies<UInt>());
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::printself(std::ostream& str) const {
   str << "Grid(" << dim << ", " << this->nb_components << ") {";
   for (auto& val : *this) {
     str << val << ", ";
   }
   str << "\b\b}";
 }
 
 /* -------------------------------------------------------------------------- */
 
 /// Class instanciation
 #define GRID_INSTANCIATE_TYPE(type)                                            \
   template class Grid<type, 1>;                                                \
   template class Grid<type, 2>;                                                \
   template class Grid<type, 3>
 
 GRID_INSTANCIATE_TYPE(Real);
 GRID_INSTANCIATE_TYPE(UInt);
 GRID_INSTANCIATE_TYPE(Complex);
 GRID_INSTANCIATE_TYPE(Int);
 GRID_INSTANCIATE_TYPE(bool);
 #undef GRID_INSTANCIATE_TYPE
 
 }  // namespace tamaas
diff --git a/src/core/grid.hh b/src/core/grid.hh
index 333fd01..4229715 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,233 +1,234 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef GRID_HH
 #define GRID_HH
 /* -------------------------------------------------------------------------- */
 #include "grid_base.hh"
 #include "tamaas.hh"
 #include <array>
 #include <numeric>
 #include <utility>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @brief Multi-dimensional & multi-component array class
  *
  * This class is a container for multi-component data stored on a multi-
  * dimensional grid.
  *
  * The access function is the parenthesis operator. For a grid of dimension d,
  * the operator takes d+1 arguments: the first d arguments are the position on
  * the grid and the last one is the component of the stored data.
  *
  * It is also possible to use the access operator with only one argument, it is
  * then considering the grid as a flat array, accessing the given cell of the
  * array.
  */
 template <typename T, UInt dim>
 class Grid : public GridBase<T> {
   template <typename... D>
   using is_valid_index = fold_trait<std::is_integral, D...>;
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Types */
   /* ------------------------------------------------------------------------ */
   using value_type = T;
   using reference = value_type&;
   static constexpr UInt dimension = dim;
 
   /* ------------------------------------------------------------------------ */
   /* Constructors */
   /* ------------------------------------------------------------------------ */
 public:
   /// Constructor by default (empty array)
   Grid();
-  /// Constructor from iterators
+  /// Constructor with shape from iterators
   template <typename RandomAccessIterator>
   Grid(RandomAccessIterator begin, RandomAccessIterator end,
        UInt nb_components);
+  /// Construct with shape from iterators on data view
   template <typename RandomAccessIterator>
-  Grid(RandomAccessIterator begin, RandomAccessIterator end,
-       UInt nb_components, value_type* data);
+  Grid(RandomAccessIterator begin, RandomAccessIterator end, UInt nb_components,
+       span<value_type> data);
   /// Constructor
   Grid(const std::array<UInt, dim>& n, UInt nb_components)
       : Grid(std::begin(n), std::end(n), nb_components) {}
   /// Wrap on given data
-  Grid(const std::array<UInt, dim>& n, UInt nb_components, value_type* data);
+  Grid(const std::array<UInt, dim>& n, UInt nb_components,
+       span<value_type> data)
+      : Grid(std::begin(n), std::end(n), nb_components, data) {}
   /// Constructor with vectors
   Grid(const std::vector<UInt>& n, UInt nb_components)
       : Grid(std::begin(n), std::end(n), nb_components) {}
   /// Constructor with initializer list
   Grid(const std::initializer_list<UInt>& n, UInt nb_components)
       : Grid(std::begin(n), std::end(n), nb_components) {}
   /// Copy constructor
-  Grid(const Grid& o);
+  Grid(const Grid& o) : GridBase<T>(), n(o.n), strides(o.strides) {}
   /// Move constructor (transfers data ownership)
-  Grid(Grid&& o) noexcept;
+  Grid(Grid&& o) noexcept
+      : GridBase<T>(std::move(o)), n(std::move(o.n)),
+        strides(std::move(o.strides)) {}
   /// Destructor
   ~Grid() override = default;
 
 private:
   /// Init from standard container
   template <typename Container>
   void init(const Container& n, UInt nb_components);
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Common operations */
   /* ------------------------------------------------------------------------ */
   /// Resize array
   void resize(const std::array<UInt, dim>& n);
   /// Resize array (from std::vector)
   void resize(const std::vector<UInt>& n);
   /// Resize array (from initializer list)
   void resize(std::initializer_list<UInt> n);
   /// Resize array (from iterators)
   template <typename ForwardIt>
   void resize(ForwardIt begin, ForwardIt end);
   /// Compute size
   inline UInt computeSize() const;
   /// Get grid dimension
   inline UInt getDimension() const override { return dim; }
   /// Compute strides
   void computeStrides();
   /// Print
   virtual void printself(std::ostream& str) const;
   /// Get sizes
   const std::array<UInt, dim>& sizes() const { return n; }
   /// Get strides
   const std::array<UInt, dim + 1>& getStrides() const { return this->strides; }
 
   /* ------------------------------------------------------------------------ */
   /* Access operators (these cannot be moved outside the class) */
   /* ------------------------------------------------------------------------ */
   /// Variadic access operator (non-const)
   template <typename... T1>
   inline std::enable_if_t<is_valid_index<T1...>::value, T&>
   operator()(T1... args) {
     /// Checking access integrity
     constexpr UInt nargs = sizeof...(T1);
     static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
                   "number of arguments in operator() does not match dimension");
 
     constexpr UInt start = (nargs == 1) ? dim : 0;
     UInt offset = unpackOffset(0, start, args...);
     return this->data[offset];
   }
 
   /// Variadic access operator
   template <typename... T1>
   inline std::enable_if_t<is_valid_index<T1...>::value, const T&>
   operator()(T1... args) const {
     /// Checking access integrity
     constexpr UInt nargs = sizeof...(T1);
     static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
                   "number of arguments in operator() does not match dimension");
 
     constexpr UInt start = (nargs == 1) ? dim : 0;
     UInt offset = unpackOffset(0, start, args...);
     return this->data[offset];
   }
 
   /// Tuple index access operator
   template <std::size_t tdim>
   inline T& operator()(std::array<UInt, tdim> tuple);
   template <std::size_t tdim>
   inline const T& operator()(std::array<UInt, tdim> tuple) const;
 
 private:
   /// Unpacking the arguments of variadic ()
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index,
                            T1... rest) const;
   /// End case for recursion
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index) const;
   /// Computing offset for a tuple index
   template <std::size_t tdim>
   inline UInt computeOffset(std::array<UInt, tdim> tuple) const;
 
   /* ------------------------------------------------------------------------ */
   /* Move/Copy operators */
   /* ------------------------------------------------------------------------ */
 public:
   using GridBase<T>::operator=;
   // = operator
   Grid& operator=(const Grid& other);
   // = operator (move)
   Grid& operator=(Grid&& other) noexcept;
 
   // Copy data from another grid
   template <typename T1>
   void copy(const Grid<T1, dim>& other);
   // Move data from another grid
   template <typename T1>
   void move(Grid<T1, dim>&& other) noexcept;
 
   template <typename Container>
-  void wrap(GridBase<T> & other, Container&& n) {
+  void wrap(GridBase<T>& other, Container&& n) {
     GridBase<T>::wrap(other);
     std::copy(n.begin(), n.end(), this->n.begin());
     this->computeStrides();
   }
 
   template <typename Container>
-  void wrap(const GridBase<T> & other, Container&& n) {
+  void wrap(const GridBase<T>& other, Container&& n) {
     GridBase<T>::wrap(other);
     std::copy(n.begin(), n.end(), this->n.begin());
     this->computeStrides();
   }
 
   // Wrap memory (non-const)
-  void wrap(Grid& other) {
-    wrap(other, other.n);
-  }
+  void wrap(Grid& other) { wrap(other, other.n); }
 
   // Wrap memory
-  void wrap(const Grid& other) {
-    wrap(other, other.n);
-  }
+  void wrap(const Grid& other) { wrap(other, other.n); }
 
   /* ------------------------------------------------------------------------ */
   /* Member variables */
   /* ------------------------------------------------------------------------ */
 protected:
   std::array<UInt, dim> n;            ///< shape of grid: size per dimension
   std::array<UInt, dim + 1> strides;  ///< strides for access
 };
 
 }  // namespace tamaas
 
 /* -------------------------------------------------------------------------- */
 /* Inline/template function definitions                                       */
 /* -------------------------------------------------------------------------- */
 
 #include "grid_tmpl.hh"
 
 /* -------------------------------------------------------------------------- */
 
 #endif  // GRID_HH
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 5a7e2ac..fc96e29 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,342 +1,344 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #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 */
   /* ------------------------------------------------------------------------ */
   using iterator = iterator_::iterator<T>;
   using const_iterator = iterator_::iterator<const T>;
   using value_type = T;
   using reference = value_type&;
 
   /// 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.resize(size); }
   /// 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);
   }
 
+  decltype(auto) view() const { return data.view(); }
+
   /* ------------------------------------------------------------------------ */
   /* 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 */
   /* ------------------------------------------------------------------------ */
 
   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/grid_tmpl.hh b/src/core/grid_tmpl.hh
index 69aee35..499a661 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,171 +1,176 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef GRID_TMPL_HH
 #define GRID_TMPL_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 
 namespace tamaas {
 
 template <typename T, UInt dim>
 template <typename RandomAccessIterator>
 Grid<T, dim>::Grid(RandomAccessIterator begin, RandomAccessIterator end,
                    UInt nb_components)
     : GridBase<T>() {
   if (end - begin != dim) {
     TAMAAS_EXCEPTION("Provided sizes (" << end - begin
                                         << ") for grid do not match dimension ("
                                         << dim << ")");
   }
 
   this->nb_components = nb_components;
   this->resize(begin, end);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename RandomAccessIterator>
 Grid<T, dim>::Grid(RandomAccessIterator begin, RandomAccessIterator end,
-                   UInt nb_components, value_type* data)
+                   UInt nb_components, span<value_type> data)
     : GridBase<T>() {
   if (end - begin != dim) {
     TAMAAS_EXCEPTION("Provided sizes (" << end - begin
                                         << ") for grid do not match dimension ("
                                         << dim << ")");
   }
 
   std::copy(begin, end, this->n.begin());
   this->nb_components = nb_components;
-  this->data.wrap(data, this->computeSize());
+
+  if (data.size() != computeSize())
+    TAMAAS_EXCEPTION("incompatible wrap span size");
+
+  this->data.wrap(data);
+  this->computeStrides();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeSize() const {
   UInt size = 1;
   for (UInt i = 0; i < dim; i++)
     size *= n[i];
   size *= this->nb_components;
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename ForwardIt>
 void Grid<T, dim>::resize(ForwardIt begin, ForwardIt end) {
   std::copy(begin, end, this->n.begin());
   UInt size = this->computeSize();
   GridBase<T>::resize(size);
   this->computeStrides();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos, UInt index,
                                        T1... rest) const {
   offset += index * strides[index_pos];
   return unpackOffset(offset, index_pos + 1, rest...);  // tail-rec bb
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index_pos,
                                        UInt index) const {
   return offset + index * strides[index_pos];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <std::size_t tdim>
 inline UInt Grid<T, dim>::computeOffset(std::array<UInt, tdim> tuple) const {
   static_assert(tdim == dim or tdim == dim + 1, "Tuple dimension is invalid");
   return std::inner_product(tuple.begin(), tuple.end(), strides.begin(), 0);
 }
 
 template <typename T, UInt dim>
 template <std::size_t tdim>
 inline T& Grid<T, dim>::operator()(std::array<UInt, tdim> tuple) {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 template <std::size_t tdim>
 inline const T& Grid<T, dim>::operator()(std::array<UInt, tdim> tuple) const {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>& Grid<T, dim>::operator=(const Grid<T, dim>& other) {
   this->copy(other);
   return *this;
 }
 
 template <typename T, UInt dim>
 Grid<T, dim>& Grid<T, dim>::operator=(Grid<T, dim>&& other) noexcept {
   this->move(std::move(other));
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::copy(const Grid<T1, dim>& other) {
   GridBase<T>::copy(other);
   this->n = other.n;
   this->strides = other.strides;
 }
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::move(Grid<T1, dim>&& other) noexcept {
   GridBase<T>::move(std::move(other));
   this->n = std::move(other.n);
   this->strides = std::move(other.strides);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Stream operator                                                            */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline std::ostream& operator<<(std::ostream& stream,
                                 const Grid<T, dim>& _this) {
   _this.printself(stream);
   return stream;
 }
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif  // GRID_TMPL_HH
diff --git a/src/core/cuda/unified_allocator.hh b/src/core/span.hh
similarity index 50%
copy from src/core/cuda/unified_allocator.hh
copy to src/core/span.hh
index 7630e9d..a4622de 100644
--- a/src/core/cuda/unified_allocator.hh
+++ b/src/core/span.hh
@@ -1,51 +1,71 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
-#ifndef UNIFIED_ALLOCATOR_HH
-#define UNIFIED_ALLOCATOR_HH
+#ifndef SPAN_HH
+#define SPAN_HH
 /* -------------------------------------------------------------------------- */
-#include <cuda_runtime_api.h>
-#include <memory>
+#include "tamaas.hh"
+#include <cstddef>
+#include <iterator>
+#include <type_traits>
 /* -------------------------------------------------------------------------- */
-
 namespace tamaas {
 
-/// Class allocating [unified
-/// memory](http://docs.nvidia.com/cuda/cuda-runtime-api/group__CUDART__MEMORY.html#group__CUDART__MEMORY_1gd228014f19cc0975ebe3e0dd2af6dd1b)
 template <typename T>
-struct UnifiedAllocator {
-  /// Allocate memory
-  static T* allocate(std::size_t n) {
-    T* p = nullptr;
-    cudaMallocManaged(&p, n * sizeof(T));
-    return p;
+struct span {
+  using element_type = T;
+  using value_type = std::remove_cv_t<T>;
+  using size_type = std::size_t;
+  using difference_type = std::ptrdiff_t;
+  using pointer = T*;
+  using const_pointer = const T*;
+  using reference = T&;
+  using const_reference = const T&;
+  using iterator = T*;
+  using reverse_iterator = std::reverse_iterator<iterator>;
+
+  constexpr size_type size() const noexcept { return size_; }
+  constexpr pointer data() const noexcept { return data_; }
+
+  constexpr iterator begin() const noexcept { return data_; }
+  constexpr iterator end() const noexcept { return data_ + size_; }
+  constexpr iterator begin() noexcept { return data_; }
+  constexpr iterator end() noexcept { return data_ + size_; }
+
+  reference operator[](size_type idx) {
+    TAMAAS_ASSERT((idx < size_), "index out of span range");
+    return data_[idx];
   }
 
-  /// Free memory
-  static void deallocate(T* p, __attribute__((unused)) std::size_t n) {
-    cudaFree(p);
+  const_reference operator[](size_type idx) const {
+    TAMAAS_ASSERT((idx < size_), "index out of span range");
+    return data_[idx];
   }
+
+  pointer data_ = nullptr;
+  size_type size_ = 0;
 };
 
 }  // namespace tamaas
 
-#endif  // UNIFIED_ALLOCATOR_HH
+/* -------------------------------------------------------------------------- */
+#endif
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index 1f15e61..836c0bf 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,199 +1,200 @@
 /**
  *  @mainpage Tamaas - A high-performance periodic contact library
  *
  *  @section Introduction
  *  Tamaas is a spectral-integral-equation based contact library. It is made
  *  with love to be fast and friendly!
  *
  *  @author Lucas Frérot <lfrerot1@jhu.edu>
  *  @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *  @author Valentine Rey <valentine.rey@univ-nantes.fr>
  *  @author Son Pham-Ba <son.phamba@epfl.ch>
  *  @author Jean-François Molinari <jean-francois.molinari@epfl.ch>
  *
  *  @section License
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef TAMAAS_HH
 #define TAMAAS_HH
 /* -------------------------------------------------------------------------- */
 
 #ifndef TAMAAS_USE_CUDA
 #define TAMAAS_USE_FFTW
 #endif
 
 // Values for fftw backends
 #define TAMAAS_FFTW_BACKEND_OMP 2
 #define TAMAAS_FFTW_BACKEND_THREADS 2
 #define TAMAAS_FFTW_BACKEND_NONE 3
 // Values for thrust backends
 #define TAMAAS_LOOP_BACKEND_OMP 1
 #define TAMAAS_LOOP_BACKEND_TBB 2
 #define TAMAAS_LOOP_BACKEND_CPP 3
 #define TAMAAS_LOOP_BACKEND_CUDA 4
 // Default loop backend is OpenMP
 #ifndef TAMAAS_LOOP_BACKEND
 #define TAMAAS_LOOP_BACKEND TAMAAS_LOOP_BACKEND_OMP
 #endif
 // Default FFTW backend is none
 #ifndef TAMAAS_FFTW_BACKEND
 #define TAMAAS_FFTW_BACKEND TAMAAS_FFTW_BACKEND_NONE
 #endif
 // If the thrust device hasn't been set, set OpenMP
 #ifndef THRUST_DEVICE_SYSTEM
 #define THRUST_DEVICE_SYSTEM THRUST_DEVICE_SYSTEM_OMP
 #endif
+
+/// Convenience macros
+#define TAMAAS_DEBUG_MSG(mesg)                                                 \
+  __FILE__ << ':' << __LINE__ << ": " << mesg << '\n'
+#define TAMAAS_EXCEPTION(mesg)                                                 \
+  {                                                                            \
+    std::stringstream sstr;                                                    \
+    sstr << TAMAAS_DEBUG_MSG("FATAL: " << mesg);                               \
+    throw ::tamaas::Exception(sstr.str());                                     \
+  }
+
+#define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
+
+#if defined(TAMAAS_DEBUG)
+#define TAMAAS_ASSERT(cond, reason)                                            \
+  do {                                                                         \
+    if (not(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; }
+
 /* -------------------------------------------------------------------------- */
 // Standard includes
 #include <exception>
 #include <iostream>
 #include <memory>
 #include <string>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 // Special thrust includes
 #include <thrust/complex.h>
 #include <thrust/random.h>
 
 #ifdef TAMAAS_USE_CUDA
 #include "cuda/unified_allocator.hh"
 #else
 #include "fftw/fftw_allocator.hh"
 #endif
 
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /// Cuda specific definitions
 #define CUDA_LAMBDA __device__ __host__
 
 #ifdef TAMAAS_USE_CUDA
 template <typename T>
 using Allocator = UnifiedAllocator<T>;
 #else
 template <typename T>
 using Allocator = FFTWAllocator<T>;
 #endif
 
 /// Common types definitions
 // If type macros have not been set, put default values
 #ifndef TAMAAS_REAL_TYPE
 #define TAMAAS_REAL_TYPE double
 #endif
 #ifndef TAMAAS_INT_TYPE
 #define TAMAAS_INT_TYPE int
 #endif
 
 using Real = TAMAAS_REAL_TYPE;           ///< default floating point type
 using Int = TAMAAS_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
 
 /// 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();
 /* -------------------------------------------------------------------------- */
 
 /// CUDA-compatible exchange function
 template <class U, class V = U>
 __device__ __host__ U exchange(U& obj, V&& new_value) {
   U old_value = std::move(obj);
   obj = std::forward<V>(new_value);
   return old_value;
 }
 
 /// Generic exception class
 class Exception : public std::exception {
 public:
   /// Constructor
   Exception(std::string mesg) : msg(std::move(mesg)) {}
 
   const char* what() const noexcept override { return msg.c_str(); }
 
   ~Exception() override = default;
 
 private:
   std::string msg;  ///< message of exception
 };
+
 /* -------------------------------------------------------------------------- */
 /// Enumeration of reduction operations
 enum class operation { plus, times, min, max };
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 
-/// Convenience macros
-#define TAMAAS_DEBUG_MSG(mesg)                                                 \
-  __FILE__ << ':' << __LINE__ << ": " << mesg << '\n'
-#define TAMAAS_EXCEPTION(mesg)                                                 \
-  {                                                                            \
-    std::stringstream sstr;                                                    \
-    sstr << TAMAAS_DEBUG_MSG("FATAL: " << mesg);                               \
-    throw ::tamaas::Exception(sstr.str());                                     \
-  }
-
-#define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
-
-#if defined(TAMAAS_DEBUG)
-#define TAMAAS_ASSERT(cond, reason)                                            \
-  do {                                                                         \
-    if (not(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 9395ec2..9d8d568 100644
--- a/src/model/elasto_plastic/residual.cpp
+++ b/src/model/elasto_plastic/residual.cpp
@@ -1,160 +1,160 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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 "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}) {
     *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());
+                      strain_increment.view());
   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/kelvin.cpp b/src/model/kelvin.cpp
index a286c6d..2e8c2e7 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,161 +1,160 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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_DEBUG_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_DEBUG_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)["displacement"].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 {
   const auto dim = trait::dimension;
   const auto shape = this->model->getDiscretization();
   Grid<Real, dim> x{shape.begin(), shape.end(),
-                    dtrait::template source_components<type>,
-                    X.getInternalData()},
+                    dtrait::template source_components<type>, X.view()},
       y{shape, dtrait::template out_components<type>};
   this->apply(x, 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/tests/test_fftw.cpp b/tests/test_fftw.cpp
index 3316933..cbbd340 100644
--- a/tests/test_fftw.cpp
+++ b/tests/test_fftw.cpp
@@ -1,237 +1,240 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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 "fftw/fftw_engine.hh"
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "grid_view.hh"
 #include "test.hh"
 
 using namespace tamaas;
 using fft = fftw::helper<Real>;
 
 /* -------------------------------------------------------------------------- */
 
+namespace fftw_impl {
 template <typename T>
 struct span {
   T* ptr;
   std::size_t size;
 
   ~span() { fftw::free(ptr); }
   const T* begin() const { return ptr; }
   const T* end() const { return ptr + size; }
   T* begin() { return ptr; }
   T* end() { return ptr + size; }
 
   operator T*() { return ptr; }
 };
+}  // namespace fftw_impl
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTEngine, FFT1D) {
   constexpr UInt size = 1000;
   FFTWEngine engine;
 
-  span<Real> data{fft::alloc_real(size), size};
-  span<fft::complex> solution{fft::alloc_complex(size / 2 + 1), size / 2 + 1};
+  fftw::span<Real> data{fft::alloc_real(size), size};
+  fftw::span<fft::complex> solution{fft::alloc_complex(size / 2 + 1),
+                                    size / 2 + 1};
 
   fftw::plan<Real> solution_plan{
       fftw::plan_1d_forward(size, data, solution, engine.flags())};
 
   std::iota(data.begin(), data.end(), 0);
   fftw::execute(solution_plan);
 
   Grid<Real, 1> grid({size}, 1);
   GridHermitian<Real, 1> result({size / 2 + 1}, 1);
 
   std::iota(grid.begin(), grid.end(), 0);
   engine.forward(grid, result);
 
   ASSERT_TRUE(compare(result, solution, AreComplexEqual()))
       << "1D FFTW transform failed";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFT2D) {
   constexpr UInt size = 100;
   constexpr UInt rsize = size * size;
   constexpr UInt csize = size * (size / 2 + 1);
   FFTWEngine engine;
 
-  span<Real> data{fft::alloc_real(rsize), rsize};
-  span<fft::complex> solution{fft::alloc_complex(csize), csize};
+  fftw::span<Real> data{fft::alloc_real(rsize), rsize};
+  fftw::span<fft::complex> solution{fft::alloc_complex(csize), csize};
 
   fftw::plan<Real> solution_plan{
       fftw::plan_2d_forward(size, size, data, solution, engine.flags())};
 
   std::iota(data.begin(), data.end(), 0);
   fftw::execute(solution_plan);
 
   Grid<Real, 2> grid({size, size}, 1);
   GridHermitian<Real, 2> result({size, size / 2 + 1}, 1);
 
   std::iota(grid.begin(), grid.end(), 0);
   engine.forward(grid, result);
 
   ASSERT_TRUE(compare(result, solution, AreComplexEqual()))
       << "2D FFTW transform failed";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFT2DBackwards) {
   const std::ptrdiff_t N0 = 20, N1 = 20;
 
   Grid<Real, 2> real({N0, N1}, 1);
   GridHermitian<Real, 2> spectral({N0, N1 / 2 + 1}, 1);
 
   real = 1.;
 
   FFTWEngine engine;
 
   engine.forward(real, spectral);
   real = 0;
   engine.backward(real, spectral);
 
   Grid<Real, 2> reference({N0, N1}, 1);
 
   reference = 1.;
 
   ASSERT_TRUE(compare(real, reference, AreFloatEqual()));
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFT1D2Comp) {
   constexpr UInt size = 20;
 
   /// 1D single component FFT should be working here
   Grid<Real, 1> grid({size}, 2), data({size}, 1);
   std::iota(grid.begin(), grid.end(), 0);
   std::iota(data.begin(), data.end(), 0);
   GridHermitian<Real, 1> result({size / 2 + 1}, 2), solution({size / 2 + 1}, 1);
 
   FFTWEngine engine;
 
   engine.forward(grid, result);
 
   std::iota(data.begin(), data.end(), 0);
   data *= 2;
   engine.forward(data, solution);
 
   const Real tol = 200 * std::numeric_limits<Real>::epsilon();
 
   ASSERT_TRUE(
       compare(make_component_view(result, 0), solution, AreComplexEqual{tol}))
       << "1D FFTW transform with 2 components failed on 1st component";
 
   data += 1;
   engine.forward(data, solution);
 
   ASSERT_TRUE(
       compare(make_component_view(result, 1), solution, AreComplexEqual{tol}))
       << "1D FFTW transform with 2 components failed on 2nd component";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFT2D3Comp) {
   constexpr UInt size = 20;
 
   /// 2D single component FFT should be working here
   Grid<Real, 2> grid({size, size}, 3), data({size, size}, 1);
   std::iota(grid.begin(), grid.end(), 0);
   std::iota(data.begin(), data.end(), 0);
   data *= 3;
 
   GridHermitian<Real, 2> result({size, size / 2 + 1}, 3),
       solution({size, size / 2 + 1}, 1);
 
   FFTWEngine engine;
 
   engine.forward(grid, result);
 
   constexpr Real tol = 5000 * std::numeric_limits<Real>::epsilon();
 
   for (UInt i = 0; i < 3; ++i) {
     engine.forward(data, solution);
     ASSERT_TRUE(
         compare(make_component_view(result, i), solution, AreComplexEqual{tol}))
         << "2D FFTW transform with 3 components failed on " << i
         << "th component";
     data += 1;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFT2DViewTransform) {
   constexpr UInt size = 20;
 
   Grid<Real, 2> data({size, size}, 1);
   GridHermitian<Real, 2> solution({size, size / 2 + 1}, 1);
   std::iota(std::begin(data), std::end(data), 0);
 
   FFTWEngine engine;
   engine.forward(data, solution);
 
   Grid<Real, 2> grid({size, size}, 3);
   auto view = make_component_view(grid, 1);
   std::iota(view.begin(), view.end(), 0);
 
   GridHermitian<Real, 2> result({size, size / 2 + 1}, 1);
   engine.forward(view, result);
 
   constexpr Real tol = 5000 * std::numeric_limits<Real>::epsilon();
 
   ASSERT_TRUE(compare(result, solution, AreComplexEqual{tol}))
       << "Fourier transform on component view fail";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFTI1D2Comp) {
   constexpr UInt size = 20;
 
   Grid<Real, 1> grid({size}, 2);
   std::iota(grid.begin(), grid.end(), 0);
   GridHermitian<Real, 1> grid_hermitian({size / 2 + 1}, 2);
   Grid<Real, 1> result({size}, 2);
 
   FFTWEngine engine;
   engine.forward(grid, grid_hermitian);
   engine.backward(result, grid_hermitian);
 
   ASSERT_TRUE(compare(grid, result, AreFloatEqual()))
       << "1D FFTI transform with 2 components failed";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTWEngine, FFTI2D3Comp) {
   constexpr UInt size = 20;
 
   Grid<Real, 2> grid({size, size}, 3);
   std::iota(grid.begin(), grid.end(), 0);
   GridHermitian<Real, 2> grid_hermitian({size, size / 2 + 1}, 3);
   Grid<Real, 2> result({size, size}, 3);
 
   FFTWEngine engine;
   engine.forward(grid, grid_hermitian);
   engine.backward(result, grid_hermitian);
 
   ASSERT_TRUE(compare(grid, result, AreFloatEqual()))
       << "2D FFTI transform with 3 components failed";
 }
diff --git a/tests/test_grid.cpp b/tests/test_grid.cpp
index d1a4c12..1e9f0a1 100644
--- a/tests/test_grid.cpp
+++ b/tests/test_grid.cpp
@@ -1,287 +1,305 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 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.hh"
 #include "grid_view.hh"
 #include "model_type.hh"
 #include "test.hh"
 
 using namespace tamaas;
 
 /* -------------------------------------------------------------------------- */
 
 // Testing 1D creation
 TEST(TestGridCreation, Creation1d) {
   Grid<Real, 1> grid{{5}, 2};
   std::array<UInt, 1> correct_size{{5}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 2> correct_strides{{2, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 2D creation
 TEST(TestGridCreation, Creation2d) {
   Grid<Real, 2> grid({5, 2}, 3);
   std::array<UInt, 2> correct_size{{5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 3> correct_strides{{6, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
 // Testing 3D creation
 TEST(TestGridCreation, Creation3d) {
   Grid<Real, 3> grid({8, 5, 2}, 3);
   std::array<UInt, 3> correct_size{{8, 5, 2}};
   ASSERT_TRUE(grid.sizes() == correct_size) << "Wrong sizes";
 
   std::array<UInt, 4> correct_strides{{30, 6, 3, 1}};
   ASSERT_TRUE(grid.getStrides() == correct_strides) << "Wrong strides";
 }
 
+/* -------------------------------------------------------------------------- */
+TEST(TestWrap, Creation) {
+  std::array<Real, 10> buffer;
+  std::array<UInt, 2> shape{5, 2};
+  span<double> view{buffer.data(), 10};
+
+  // Fill values
+  std::iota(std::begin(buffer), std::end(buffer), 0);
+
+  // Wrap grid on view
+  Grid<Real, 2> grid(shape.begin(), shape.end(), 1, view);
+
+  grid(1, 1) = 1;
+
+  ASSERT_TRUE(compare(view, buffer));
+  ASSERT_TRUE(compare(grid, buffer));
+}
+
 /* -------------------------------------------------------------------------- */
 
 // Testing iterators in STL function iota
 TEST(TestGridIterators, Iota) {
   constexpr UInt N = 20;
   Grid<UInt, 1> grid({N}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   const UInt* p = grid.getInternalData();
   for (UInt i = 0; i < N; ++i)
     ASSERT_TRUE(p[i] == i) << "Iota fill failed";
 }
 
 // Testing filling grid with OpenMP loop
 TEST(TestGridIterators, OpenMPLoops) {
   Grid<UInt, 1> grid({20}, 1);
 
 #ifndef TAMAAS_USE_CUDA
 #pragma omp parallel for
 #endif
   for (auto it = grid.begin(); it < grid.end(); ++it) {
     UInt i = it - grid.begin();
     *it = i;
   }
 
   std::vector<UInt> correct(20);
   std::iota(correct.begin(), correct.end(), 0);
 
   ASSERT_TRUE(compare(grid, correct)) << "Fill using OpenMP loop failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing access operator
 TEST(TestGridAccess, 3D) {
   Grid<UInt, 3> grid({1, 3, 5}, 3);
   std::iota(grid.begin(), grid.end(), 0);
 
   for (UInt i = 0; i < grid.sizes()[0]; ++i)
     for (UInt j = 0; j < grid.sizes()[1]; ++j)
       for (UInt k = 0; k < grid.sizes()[2]; ++k)
         for (UInt m = 0; m < grid.getNbComponents(); ++m)
           ASSERT_EQ(grid(i, j, k, m),
                     m + 3 * k + (3 * 5) * j + (3 * 5 * 3) * i);
 
   Grid<UInt, 2> grid1({5, 3}, 1);
   grid1(2, 1, 0) = 1;
   ASSERT_EQ(grid1(2 * 3 + 1), 1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 // Testing scalar simple loop-based operators
 TEST(TestGridOperators, LoopOperators) {
   Grid<UInt, 1> grid({20}, 1);
   grid = 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 1;
   })) << "Assignement operator failed";
   grid += 1;
   EXPECT_TRUE(std::all_of(grid.begin(), grid.end(), [](UInt x) {
     return x == 2;
   })) << "Plus-equal operator failed";
 }
 
 // Testing loop-based functions with reductions
 TEST(TestGridOperators, Stats) {
   constexpr UInt n = 20;
   const UInt N = mpi::allreduce<operation::plus>(n);
 
   Grid<Real, 1> grid({n}, 1);
   std::iota(grid.begin(), grid.end(), 0);
 
   std::vector<Real> v(n);
   std::iota(v.begin(), v.end(), 0);
 
   const auto b = v.cbegin(), e = v.cend();
   Real min = mpi::allreduce<operation::min>(*std::min_element(b, e));
   Real max = mpi::allreduce<operation::max>(*std::max_element(b, e));
   Real sum = mpi::allreduce<operation::plus>(std::accumulate(b, e, 0.));
   Real mea = sum / N;
 
   std::transform(b, e, v.begin(),
                  [mea](Real v) { return (v - mea) * (v - mea); });
   Real acc = mpi::allreduce<operation::plus>(std::accumulate(b, e, 0.));
   Real var = acc / (N - 1);
 
   EXPECT_DOUBLE_EQ(grid.min(), min) << "Minimum function failed";
   EXPECT_DOUBLE_EQ(grid.max(), max) << "Maximum function failed";
   EXPECT_DOUBLE_EQ(grid.sum(), sum) << "Sum function failed";
   EXPECT_DOUBLE_EQ(grid.mean(), mea) << "Mean function failed";
   EXPECT_DOUBLE_EQ(grid.var(), var) << "Var function failed";
 }
 
 /* -------------------------------------------------------------------------- */
 
 TEST(TestGridView, View2D) {
   Grid<UInt, 3> grid3d({10, 10, 10}, 3);
   std::iota(grid3d.begin(), grid3d.end(), 0);
 
   for (UInt i : Loop::range(10)) {
     auto grid_view = make_view(grid3d, i);
     std::vector<UInt> reference(300);
     std::iota(reference.begin(), reference.end(), 300 * i);
     EXPECT_TRUE(
         std::equal(grid_view.begin(), grid_view.end(), reference.begin()))
         << "View on slice fail";
   }
 }
 
 TEST(TestGridView, ViewOnComponent) {
   Grid<UInt, 3> grid3d({10, 9, 8}, 3), solution({10, 9, 8}, 3);
   grid3d = solution = 0;
 
   auto view = make_component_view(grid3d, 2u);
   std::iota(view.begin(), view.end(), 0);
 
   for (UInt i = 0; i < solution.getNbPoints(); ++i)
     solution(3 * i + 2) = i;
 
   EXPECT_TRUE(std::equal(grid3d.begin(), grid3d.end(), solution.begin()))
       << "View on components fail";
 
   view(9, 8, 7) = 1;
   ASSERT_EQ(grid3d(9, 8, 7, 2), 1);
 }
 
 TEST(TestGridView, ViewCopy) {
   Grid<UInt, 3> grid3d({10, 10, 10}, 2);
   Grid<UInt, 1> copy, solution({10}, 1);
   auto view_on_slice = make_view(grid3d, 3, 4);
   Grid<UInt, 1>* view_on_slice_ptr = &view_on_slice;
   auto view = make_component_view(*view_on_slice_ptr, 1u);
   std::iota(view.begin(), view.end(), 0);
   std::iota(solution.begin(), solution.end(), 0);
   copy = view;
   EXPECT_TRUE(std::equal(copy.begin(), copy.end(), solution.begin()))
       << "Copy of view fail";
 }
 
 TEST(TestGridView, HigherOrderView) {
   Grid<UInt, 2> grid({10, 10}, 2);
   std::iota(grid.begin(), grid.end(), 0);
   GridView<Grid, UInt, 2, 3> view(grid, {}, -1);
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), view.begin()))
       << "Higher order view fail";
 }
 
 /* -------------------------------------------------------------------------- */
 
 struct BroadcastSet123 {
   CUDA_LAMBDA inline void operator()(VectorProxy<UInt, 3> v) {
     v(0) = 1;
     v(1) = 2;
     v(2) = 3;
   }
 };
 
 TEST(TestGridOperators, BroadcastOperators) {
   Grid<UInt, 2> grid({10, 10}, 3), solution({10, 10}, 3);
   auto set = BroadcastSet123();
 
   // Filling up solution
   Loop::loop(set, range<VectorProxy<UInt, 3>>(solution));
 
   Vector<UInt, 3> v;
   v(0) = 1;
   v(1) = 2;
   v(2) = 3;
 
   grid = v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast assignement failure";
 
   v(1) = 1;
   v(2) = 1;
 
   solution += 1;
   grid += v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast += failure";
 
   solution -= 1;
   grid -= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast -= failure";
 
   v(1) = 2;
   v(2) = 3;
 
   solution *= solution;
   grid *= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast *= failure";
 
   Loop::loop(set, range<VectorProxy<UInt, 3>>(solution));
   grid /= v;
 
   EXPECT_TRUE(std::equal(grid.begin(), grid.end(), solution.begin()))
       << "Broadcast /= failure";
 }
 
 TEST(TestAllocation, HardType) {
   std::array<UInt, 2> dims{4, 5};
   auto grid = allocateGrid<model_type::basic_2d, true, Real>(dims, 1);
   static_assert(
       std::is_same<decltype(grid), std::unique_ptr<Grid<Real, 2>>>::value,
       "allocation type incorrect");
   EXPECT_TRUE(std::equal(dims.begin(), dims.end(), grid->sizes().begin()));
   EXPECT_EQ(grid->getNbComponents(), 1);
 }
 
 TEST(TestAllocation, SoftType) {
   std::array<UInt, 2> dims{4, 5};
   auto grid = allocateGrid<true, Real>(model_type::basic_2d, dims, 1);
   static_assert(
       std::is_same<decltype(grid), std::unique_ptr<GridBase<Real>>>::value,
       "allocation type incorrect");
   EXPECT_EQ(grid->dataSize(), 20);
   EXPECT_EQ(grid->getNbComponents(), 1);
 }