diff --git a/src/core/grid.hh b/src/core/grid.hh
index 75db563..333fd01 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,230 +1,233 @@
 /*
  *  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
   template <typename RandomAccessIterator>
   Grid(RandomAccessIterator begin, RandomAccessIterator end,
        UInt nb_components);
+  template <typename RandomAccessIterator>
+  Grid(RandomAccessIterator begin, RandomAccessIterator end,
+       UInt nb_components, 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);
   /// 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);
   /// Move constructor (transfers data ownership)
   Grid(Grid&& o) noexcept;
   /// 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) {
     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) {
     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);
   }
 
   // Wrap memory
   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_tmpl.hh b/src/core/grid_tmpl.hh
index 12293b7..69aee35 100644
--- a/src/core/grid_tmpl.hh
+++ b/src/core/grid_tmpl.hh
@@ -1,153 +1,171 @@
 /*
  *  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)
+    : 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());
+}
+
+/* -------------------------------------------------------------------------- */
+
 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/model/kelvin.cpp b/src/model/kelvin.cpp
index fdf2bad..a286c6d 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,161 +1,161 @@
 /*
  *  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;
-  Grid<Real, dim> x;
-  Grid<Real, dim> y{this->model->getDiscretization(),
-                    dtrait::template out_components<type>};
-  X.setNbComponents(dtrait::template source_components<type>);
-  x.wrap(X, this->model->getDiscretization());
+  const auto shape = this->model->getDiscretization();
+  Grid<Real, dim> x{shape.begin(), shape.end(),
+                    dtrait::template source_components<type>,
+                    X.getInternalData()},
+      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