diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index 63124a4..1d19532 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,134 +1,131 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "kelvin.hh"
 #include "elasto_plastic_model.hh"
 #include "influence.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* Volume 2D implementation */
 /* -------------------------------------------------------------------------- */
 template <>
 void Kelvin<model_type::volume_2d>::applyVolumeForcePotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   Real nu = model->getPoissonRatio(), mu = model->getShearModulus();
   VectorProxy<const Real, trait::dimension> domain(model->getSystemSize()[0]);
   influence::Kelvin<trait::dimension, 0> kelvin(mu, nu, domain);
 
-  auto apply = [this, &kelvin](decltype(source_buffers)& source_buffers,
-			       decltype(disp_buffers)& disp_buffers) {
+  auto apply = [this, &kelvin](UInt i, decltype(source_buffers)& source_buffers,
+                               decltype(disp_buffer)& displacement) {
     constexpr UInt dim = trait::dimension;
     const Real L = this->model->getSystemSize().front();
     const UInt N = this->model->getDiscretization().front();
     const Real dl = L / N;
-    // Do the matrix vector product
-    for (UInt i = 0; i < N; ++i) {
-      const Real xi = i * dl;
-      // Compute displacements u_i
-      auto& displacement = disp_buffers[i];
-      displacement = 0;
-      for (UInt j = 0; j < N; ++j) {
-	const Real yj = j * dl;
-	auto& source = source_buffers[j];
+    const Real xi = i * dl;
+
+    // Compute displacements u_i
+    displacement = 0;
+
+    for (UInt j : Loop::range(N)) {
+      const Real yj = j * dl;
+      auto& source = source_buffers[j];
 
 #define VOLUME_FORCE_POTENTIAL(yj_xi)                                          \
   Loop::stridedLoop(                                                           \
       [&kelvin, xi, yj, dl](VectorProxy<Complex, dim>&& u,                     \
-			    VectorProxy<Complex, dim>&& f,                     \
-			    VectorProxy<const Real, dim - 1>&& q) {            \
-	constexpr UInt dim = trait::dimension;                                 \
-	Matrix<Complex, dim, dim> F;                                           \
-	influence::KelvinIntegrator::integrate<dim, 0, yj_xi>(kelvin, yj, xi,  \
-							      dl, q, F);       \
-	u += F * f;                                                            \
+                            VectorProxy<Complex, dim>&& f,                     \
+                            VectorProxy<const Real, dim - 1>&& q) {            \
+        constexpr UInt dim = trait::dimension;                                 \
+        Matrix<Complex, dim, dim> F;                                           \
+        influence::KelvinIntegrator::integrate<dim, 0, yj_xi>(kelvin, yj, xi,  \
+                                                              dl, q, F);       \
+        u += F * f;                                                            \
       },                                                                       \
       displacement, source, this->wavevectors)
 
-	if (j > i) {
-	  VOLUME_FORCE_POTENTIAL(1);
-	} else if (j == i) {
-	  VOLUME_FORCE_POTENTIAL(0);
-	} else {
-	  VOLUME_FORCE_POTENTIAL(-1);
-	}
-
-#undef VOLUME_FORCE_POTENTIAL
+      if (j > i) {
+        VOLUME_FORCE_POTENTIAL(1);
+      } else if (j == i) {
+        VOLUME_FORCE_POTENTIAL(0);
+      } else {
+        VOLUME_FORCE_POTENTIAL(-1);
       }
-
-      // Setting fundamental frequency to zero
-      VectorProxy<Complex, dim> u_fundamental(displacement(0));
-      u_fundamental = 0;
+#undef VOLUME_FORCE_POTENTIAL
     }
+
+    // Setting fundamental frequency to zero
+    VectorProxy<Complex, dim> u_fundamental(displacement(0));
+    u_fundamental = 0;
   };
 
   this->fourierApply(apply, source, out);
 }
 
 template <>
 void Kelvin<model_type::volume_2d>::applyVolumeStressPotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   //
 }
 
 template <>
 void Kelvin<model_type::volume_2d>::applyHypersingularStressPotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   //
 }
 
 /* -------------------------------------------------------------------------- */
 /* Volume 1D implementation */
 /* -------------------------------------------------------------------------- */
 
 template <>
 void Kelvin<model_type::volume_1d>::applyVolumeForcePotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   //
 }
 
 template <>
 void Kelvin<model_type::volume_1d>::applyVolumeStressPotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   //
 }
 
 template <>
 void Kelvin<model_type::volume_1d>::applyHypersingularStressPotential(
     GridBase<Real>& source, GridBase<Real>& out) const {
   //
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Kelvin<model_type::volume_1d>;
 template class Kelvin<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
diff --git a/src/model/ve_engine.cpp b/src/model/ve_engine.cpp
index ad5bbde..3ea1497 100644
--- a/src/model/ve_engine.cpp
+++ b/src/model/ve_engine.cpp
@@ -1,73 +1,75 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "ve_engine.hh"
 #include "model.hh"
 #include <algorithm>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 template <model_type type>
 VEEngineTmpl<type>::VEEngineTmpl(Model * model): VEEngine(model) {
     // Copy horizontal sizes
     std::array<UInt, trait::boundary_dimension> sizes;
     std::copy(model->getDiscretization().begin() + 1,
               model->getDiscretization().end(), sizes.begin());
 
     auto hermitian_sizes =
         GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
             sizes);
 
     wavevectors = FFTransform<Real, trait::boundary_dimension>::
         template computeFrequencies<true>(hermitian_sizes);
 
     // Normalize wavevectors
     VectorProxy<const Real, trait::boundary_dimension> boundary_domain{
         model->getSystemSize()[1]};
     wavevectors *= 2 * M_PI;
     wavevectors /= boundary_domain;
 
     /// Initializing buffers
-    auto initialize = [this, &hermitian_sizes](std::vector<BufferType> & buffers) {
+    auto initialize = [this,
+                       &hermitian_sizes](std::vector<BufferType>& buffers) {
       buffers.resize(this->model->getDiscretization()[0]);
       std::for_each(buffers.begin(), buffers.end(),
                     [&hermitian_sizes](BufferType& buffer) {
-		      buffer.setNbComponents(trait::components);
+                      buffer.setNbComponents(trait::components);
                       buffer.resize(hermitian_sizes);
                     });
     };
 
     initialize(source_buffers);
-    initialize(disp_buffers);
+    disp_buffer.setNbComponents(trait::components);
+    disp_buffer.resize(hermitian_sizes);
 }
 
 /* ------------------------------------------------------------------------ */
 /* Template instanciation */
 /* ------------------------------------------------------------------------ */
 template class VEEngineTmpl<model_type::volume_1d>;
 template class VEEngineTmpl<model_type::volume_2d>;
 
 __END_TAMAAS__
diff --git a/src/model/ve_engine.hh b/src/model/ve_engine.hh
index 7c8e9f3..3a37347 100644
--- a/src/model/ve_engine.hh
+++ b/src/model/ve_engine.hh
@@ -1,125 +1,124 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2017 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __VE_ENGINE_HH__
 #define __VE_ENGINE_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid_base.hh"
 #include "grid_view.hh"
 #include "grid_hermitian.hh"
 #include "model_type.hh"
 #include "fft_plan_manager.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 /// Forward declaration
 class Model;
 
 /**
  * @brief Volume equation engine class. Applies the operators for computation
  * of displacements and strains due to residual/eigen strains
  */
 class VEEngine {
 public:
   /// Constructor
   VEEngine(Model* model) : model(model) {}
   /// Destructor
   virtual ~VEEngine() = default;
 
   // Methods to override
 public:
   /// Compute displacement field in volume to to a volume force distribution
   virtual void applyVolumeForcePotential(GridBase<Real>& source,
 					 GridBase<Real>& out) const = 0;
   /// Compute displacement field in volume due to an initial stress distribution
   virtual void applyVolumeStressPotential(GridBase<Real>& source,
                                           GridBase<Real>& out) const = 0;
   /// Compute strain field in volume due to an initial stress distribution
   virtual void applyHypersingularStressPotential(GridBase<Real>& source,
                                                  GridBase<Real>& out) const = 0;
 
 protected:
   Model* model;
 };
 
 /**
  * @brief Templated class to manage application of functor on fields in Fourier
  * space.
  */
 template <model_type type>
 class VEEngineTmpl : public VEEngine {
   using trait = model_type_traits<type>;
 public:
   VEEngineTmpl(Model * model);
 
 protected:
   /// Function to handle layer-by-layer Fourier treatment of operators
   template <typename Func>
   void fourierApply(Func func, GridBase<Real>& in, GridBase<Real>& out) const;
 
 protected:
   Grid<Real, trait::boundary_dimension> wavevectors;
   using BufferType = GridHermitian<Real, trait::boundary_dimension>;
   mutable std::vector<BufferType> source_buffers;
-  mutable std::vector<BufferType> disp_buffers;
+  mutable BufferType disp_buffer;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation */
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 template <typename Func>
 void VEEngineTmpl<type>::fourierApply(Func func, GridBase<Real>& in,
                                       GridBase<Real>& out) const {
   constexpr UInt dim = trait::dimension;
   Grid<Real, dim>& i = dynamic_cast<decltype(i)>(in);
   Grid<Real, dim>& o = dynamic_cast<decltype(o)>(out);
 
   TAMAAS_ASSERT(i.sizes().front() == o.sizes().front(),
                 "Number of layers does not match");
 
   for (UInt layer : Loop::range(i.sizes().front())) {
     auto in_layer = make_view(i, layer);
     FFTPlanManager::get()
         .createPlan(in_layer, source_buffers[layer])
         .forwardTransform();
   }
 
-  func(source_buffers, disp_buffers);
-
   for (UInt layer : Loop::range(i.sizes().front())) {
+    func(layer, source_buffers, disp_buffer);
     auto out_layer = make_view(o, layer);
     FFTPlanManager::get()
-        .createPlan(out_layer, disp_buffers[layer])
+        .createPlan(out_layer, disp_buffer)
         .backwardTransform();
   }
 }
 
 __END_TAMAAS__
 
 #endif  // __VE_ENGINE_HH__