diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index 956649a..3e6e80f 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,114 +1,87 @@
 /**
  * @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"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* Constructors */
 /* -------------------------------------------------------------------------- */
 template <>
 Kelvin<model_type::volume_2d, 0>::Kelvin(Model* model)
     : VolumePotential(model) {
   initialize(trait::components, trait::components);
 }
 
 template <>
 Kelvin<model_type::volume_2d, 1>::Kelvin(Model* model)
     : VolumePotential(model) {
   initialize(trait::components * trait::components, trait::components);
 }
 
 template <>
 Kelvin<model_type::volume_2d, 2>::Kelvin(Model* model)
     : VolumePotential(model) {
   initialize(trait::components * trait::components,
              trait::components * trait::components);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Operator implementation */
 /* -------------------------------------------------------------------------- */
 template <model_type type, UInt derivative>
 void Kelvin<type, derivative>::applyIf(GridBase<Real>& source,
                                        GridBase<Real>& out,
                                        filter_t pred) const {
   Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
   influence::Kelvin<trait::dimension, derivative> kelvin(mu, nu);
   detail::KelvinHelper<type, decltype(kelvin)> helper;
 
   auto apply = [&](auto&& source_buffer, auto&& out_buffer) {
     // Setting initial layer to zero
     out_buffer.front() = 0;
 
     helper.applyIntegral(source_buffer, out_buffer, this->wavevectors,
                          this->model->getSystemSize().front(), kelvin);
     helper.applyFreeTerm(source_buffer, out_buffer, kelvin);
     helper.makeFundamentalGreatAgain(source_buffer, out_buffer, kelvin);
   };
 
   this->fourierApplyIf(apply, source, out, pred);
 }
 
-/* -------------------------------------------------------------------------- */
-template <model_type type, UInt derivative>
-void Kelvin<type, derivative>::initialize(UInt source_components,
-                                          UInt out_components) {
-  auto hermitian_sizes =
-      GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
-          this->model->getBoundaryDiscretization());
-
-  // Initializing buffers
-  this->source_buffer.resize(this->model->getDiscretization()[0]);
-  this->disp_buffer.resize(this->model->getDiscretization()[0]);
-
-  // Resizing source buffer
-  std::for_each(this->source_buffer.begin(), this->source_buffer.end(),
-                [&](auto& buffer) {
-                  buffer.setNbComponents(source_components);
-                  buffer.resize(hermitian_sizes);
-                });
-
-  // Resizing output buffer
-  std::for_each(this->disp_buffer.begin(), this->disp_buffer.end(),
-                [&](auto& buffer) {
-                  buffer.setNbComponents(out_components);
-                  buffer.resize(hermitian_sizes);
-                });
-}
-
 /* -------------------------------------------------------------------------- */
 /* 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>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
diff --git a/src/model/kelvin.hh b/src/model/kelvin.hh
index c578423..3da3cb1 100644
--- a/src/model/kelvin.hh
+++ b/src/model/kelvin.hh
@@ -1,71 +1,67 @@
 /**
  * @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 KELVIN_HH
 #define KELVIN_HH
 /* -------------------------------------------------------------------------- */
 #include "grid_hermitian.hh"
 #include "model_type.hh"
 #include "volume_potential.hh"
 #include "influence.hh"
 #include "kelvin_helper.hh"
 #include "integration/accumulator.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 /**
  * @brief Kelvin tensor
  */
 template <model_type type, UInt derivative>
 class Kelvin : public VolumePotential<type> {
   static_assert(type == model_type::volume_1d || type == model_type::volume_2d,
                 "Only volume types are supported");
   using trait = model_type_traits<type>;
   using parent = VolumePotential<type>;
 
 protected:
   using KelvinInfluence = influence::Kelvin<trait::dimension, derivative>;
   using ktrait = detail::KelvinTrait<KelvinInfluence>;
   using Source = typename ktrait::source_t;
   using Out = typename ktrait::out_t;
   using filter_t = typename parent::filter_t;
 
 public:
   /// Constructor
   Kelvin(Model* model);
 
-protected:
-  void initialize(UInt source_components, UInt out_components);
-
-public:
   /// Apply the Kelvin-tensor_order operator
   void applyIf(GridBase<Real>& source, GridBase<Real>& out,
                filter_t pred) const override;
 };
 
 __END_TAMAAS__
 
 #endif  // KELVIN_HH
diff --git a/src/model/kelvin_helper.hh b/src/model/kelvin_helper.hh
index 397af9d..724e1a4 100644
--- a/src/model/kelvin_helper.hh
+++ b/src/model/kelvin_helper.hh
@@ -1,227 +1,227 @@
 /**
  * @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 KELVIN_HELPER_HH
 #define KELVIN_HELPER_HH
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "influence.hh"
 #include "integration/accumulator.hh"
 #include "model.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 namespace detail {
 
 /* -------------------------------------------------------------------------- */
 /// Trait for kelvin local types
 template <typename T>
 struct KelvinTrait;
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 0>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 1>> {
   using source_t = MatrixProxy<Complex, dim, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Kelvin<dim, 2>> {
   using source_t = MatrixProxy<Complex, dim, dim>;
   using out_t = MatrixProxy<Complex, dim, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Boussinesq<dim, 0>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = VectorProxy<Complex, dim>;
 };
 
 template <UInt dim>
 struct KelvinTrait<influence::Boussinesq<dim, 1>> {
   using source_t = VectorProxy<Complex, dim>;
   using out_t = MatrixProxy<Complex, dim, dim>;
 };
 /* -------------------------------------------------------------------------- */
 
 /// Helper to apply integral representation on output
 template <model_type type, typename kelvin_t,
           typename = typename KelvinTrait<kelvin_t>::source_t>
 struct KelvinHelper {
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   static constexpr UInt bdim = trait::boundary_dimension;
 
   using BufferType = GridHermitian<Real, bdim>;
   using source_t = typename KelvinTrait<kelvin_t>::source_t;
   using out_t = typename KelvinTrait<kelvin_t>::out_t;
 
   /// Apply the regular part of Kelvin to source and sum into output
   void applyIntegral(std::vector<BufferType>& source,
                      std::vector<BufferType>& out,
                      const Grid<Real, bdim>& wavevectors, Real domain_size,
                      const kelvin_t& kelvin) {
     accumulator.makeUniformMesh(source.size(), domain_size);
 
     auto waverange = range<VectorProxy<const Real, bdim>>(wavevectors);
 
     for (auto [l, xl_, acc_g0, acc_g1] :
          accumulator.forward(source, wavevectors)) {
       Loop::loop(
           [xl = xl_, &kelvin](auto qv, auto u, auto acc_g0, auto acc_g1) {
             Real q = qv.l2norm();
             auto tmp = kelvin.template applyU1<false, true>(qv, acc_g0);
             tmp *= -q * xl;
             tmp += kelvin.template applyU0<false, true>(qv, acc_g0);
             tmp += kelvin.template applyU1<false, true>(qv, acc_g1);
             tmp *= std::exp(-q * xl);
             u += tmp;
           },
           waverange, range<out_t>(out[l]), range<source_t>(acc_g0),
           range<source_t>(acc_g1));
     }
 
     for (auto [l, xl_, acc_g0, acc_g1] :
          accumulator.backward(source, wavevectors)) {
       Loop::loop(
           [xl = xl_, &kelvin](auto qv, auto u, auto acc_g0, auto acc_g1) {
             Real q = qv.l2norm();
             auto tmp = kelvin.template applyU1<true, true>(qv, acc_g0);
             tmp *= -q * xl;
             tmp += kelvin.template applyU0<true, true>(qv, acc_g0);
             tmp += kelvin.template applyU1<true, true>(qv, acc_g1);
             tmp *= std::exp(q * xl);
             u += tmp;
           },
           waverange, range<out_t>(out[l]), range<source_t>(acc_g0),
           range<source_t>(acc_g1));
     }
   }
 
   /// Apply free term of distribution derivative
   void applyFreeTerm(std::vector<BufferType>& /*source*/,
                      std::vector<BufferType>& /*out*/,
                      const kelvin_t& /*kelvin*/) {}
 
   /// Making the output free of communist NaN
   /// TODO change interface
   void makeFundamentalGreatAgain(std::vector<BufferType>& /*source*/,
                                  std::vector<BufferType>& /*out*/,
                                  const kelvin_t& /*kelvin*/) {}
 
 protected:
   Accumulator<type, source_t> accumulator;
 };
 
 /// Applying free term for double gradient of Kelvin
 template <>
 inline void
 KelvinHelper<model_type::volume_2d, influence::Kelvin<3, 2>>::applyFreeTerm(
     std::vector<BufferType>& source, std::vector<BufferType>& out,
     const influence::Kelvin<3, 2>& kelvin) {
   for (UInt l : Loop::range(source.size())) {
     Loop::loop(
         [&kelvin](auto u, auto w) { u += kelvin.applyDiscontinuityTerm(w); },
         range<out_t>(out[l]), range<source_t>(source[l]));
   }
 }
 
 /// Setting fundamental mode to zero: we've got no other choice here
 template <>
 inline void KelvinHelper<model_type::volume_2d, influence::Kelvin<3, 0>>::
     makeFundamentalGreatAgain(std::vector<BufferType>& /*source*/,
                               std::vector<BufferType>& out,
                               const influence::Kelvin<3, 0>& /*kelvin*/) {
   for (auto&& layer : out) {
     out_t u(&layer(0));
     u = 0;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
-class SurfaceTractionHelper {
+struct SurfaceTractionHelper {
   using trait = model_type_traits<type>;
   static constexpr UInt dim = trait::dimension;
   static constexpr UInt bdim = trait::boundary_dimension;
 
   using BufferType = GridHermitian<Real, bdim>;
   using kelvin_t = influence::Kelvin<3, 2>;
   using source_t = typename KelvinTrait<kelvin_t>::source_t;
   using out_t = typename KelvinTrait<kelvin_t>::out_t;
 
   /// Compute surface tractions due to eigenstress distribution
   template <bool apply_q_power>
   void computeSurfaceTractions(std::vector<BufferType>& source,
-                               GridHermitian<Real, bdim>& tractions,
+                               BufferType& tractions,
                                const Grid<Real, bdim>& wavevectors,
                                Real domain_size, const kelvin_t& kelvin,
                                const influence::ElasticHelper<dim>& el) {
     accumulator.makeUniformMesh(source.size(), domain_size);
 
     // only doing backward loop because forward loop gives 0 for l = 0
     for (auto [l, xl_, acc_g0, acc_g1] :
          accumulator.backward(source, wavevectors)) {
       // Just skip all layers and let the accumulator accumulate
       if (l != 0)
         continue;
 
       Loop::loop(
           [xl = xl_, &kelvin, &el](auto qv, auto t, auto acc_g0, auto acc_g1) {
             // computing strains
             auto strains =
                 kelvin.template applyU0<true, apply_q_power>(qv, acc_g0);
             strains += kelvin.template applyU1<true, apply_q_power>(qv, acc_g1);
 
             // computing traction vector
             Vector<Real, dim> normal{{{0, 0, -1}}};
             t = el(strains) * normal;
           },
           range<VectorProxy<const Real, bdim>>(wavevectors),
           range<VectorProxy<Real, dim>>(tractions), range<source_t>(acc_g0),
           range<source_t>(
               acc_g1));  // no need for headless ranges: should not be singular
     }
   }
 
 protected:
   Accumulator<type, source_t> accumulator;
 };
 
 }  // namespace detail
 
 __END_TAMAAS__
 
 #endif  // KELVIN_HELPER_HH
diff --git a/src/model/volume_potential.cpp b/src/model/volume_potential.cpp
index 30a7cdd..d289e2d 100644
--- a/src/model/volume_potential.cpp
+++ b/src/model/volume_potential.cpp
@@ -1,68 +1,96 @@
 /**
  * @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 "volume_potential.hh"
 #include "model.hh"
 #include <algorithm>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 template <model_type type>
-VolumePotential<type>::VolumePotential(Model * model): IntegralOperator(model) {
-    // Copy horizontal sizes
-    std::array<UInt, trait::boundary_dimension> sizes;
-    std::copy(model->getDiscretization().begin() + 1,
-              model->getDiscretization().end(), sizes.begin());
+VolumePotential<type>::VolumePotential(Model* model) : IntegralOperator(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);
+  auto hermitian_sizes =
+      GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
+          sizes);
 
-    wavevectors = FFTransform<Real, trait::boundary_dimension>::
-        template computeFrequencies<true>(hermitian_sizes);
+  wavevectors =
+      FFTransform<Real, trait::boundary_dimension>::template computeFrequencies<
+          true>(hermitian_sizes);
 
-    // Normalize wavevectors
-    auto system_size = model->getBoundarySystemSize();
-    VectorProxy<const Real, trait::boundary_dimension> boundary_domain{
-        system_size[0]};
-    wavevectors *= 2 * M_PI;
-    wavevectors /= boundary_domain;
-    wavevectors *= -1.;  // < this is important for the convolution computation
+  // Normalize wavevectors
+  auto system_size = model->getBoundarySystemSize();
+  VectorProxy<const Real, trait::boundary_dimension> boundary_domain{
+      system_size[0]};
+  wavevectors *= 2 * M_PI;
+  wavevectors /= boundary_domain;
+  wavevectors *= -1.;  // < this is important for the convolution computation
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type>
 model_type VolumePotential<type>::getType() const {
   return this->model->getType();
 }
 
+/* -------------------------------------------------------------------------- */
+template <model_type type>
+void VolumePotential<type>::initialize(UInt source_components,
+                                       UInt out_components) {
+  auto hermitian_sizes =
+      GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
+          this->model->getBoundaryDiscretization());
+
+  // Initializing buffers
+  this->source_buffer.resize(this->model->getDiscretization()[0]);
+  this->disp_buffer.resize(this->model->getDiscretization()[0]);
+
+  // Resizing source buffer
+  std::for_each(this->source_buffer.begin(), this->source_buffer.end(),
+                [&](auto& buffer) {
+                  buffer.setNbComponents(source_components);
+                  buffer.resize(hermitian_sizes);
+                });
+
+  // Resizing output buffer
+  std::for_each(this->disp_buffer.begin(), this->disp_buffer.end(),
+                [&](auto& buffer) {
+                  buffer.setNbComponents(out_components);
+                  buffer.resize(hermitian_sizes);
+                });
+}
+
 /* -------------------------------------------------------------------------- */
 /* Template instanciation */
 /* -------------------------------------------------------------------------- */
 template class VolumePotential<model_type::volume_1d>;
 template class VolumePotential<model_type::volume_2d>;
 
 __END_TAMAAS__
diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh
index d54e6f9..f383d56 100644
--- a/src/model/volume_potential.hh
+++ b/src/model/volume_potential.hh
@@ -1,130 +1,134 @@
 /**
  * @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 VOLUME_POTENTIAL_HH
 #define VOLUME_POTENTIAL_HH
 /* -------------------------------------------------------------------------- */
 #include "fft_plan_manager.hh"
 #include "grid_hermitian.hh"
 #include "grid_view.hh"
 #include "integral_operator.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 #include <functional>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 /**
  * @brief Volume potential operator class. Applies the operators for computation
  * of displacements and strains due to residual/eigen strains
  */
 template <model_type type>
 class VolumePotential : public IntegralOperator {
   using trait = model_type_traits<type>;
+
 protected:
   using filter_t = std::function<bool(UInt)>;
 
 public:
   VolumePotential(Model* model);
 
   /// Update from model (does nothing)
   void updateFromModel() override {}
 
   /// Kind
   IntegralOperator::kind getKind() const override {
     return IntegralOperator::neumann;
   }
 
   /// Type
   model_type getType() const override;
 
   /// Applying to all source
   void apply(GridBase<Real>& input, GridBase<Real>& output) const override {
     applyIf(input, output, [](UInt) { return true; });
   }
 
 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 {
     fourierApplyIf(func, in, out, [](UInt) { return true; });
   }
 
   /// Function to handle layer-by-layer Fourier treatment of operators
   template <typename Func>
   void fourierApplyIf(Func func, GridBase<Real>& in, GridBase<Real>& out,
                       filter_t pred) const;
 
+  /// Initialize fourier buffers
+  void initialize(UInt source_components, UInt out_components);
+
 protected:
   Grid<Real, trait::boundary_dimension> wavevectors;
   using BufferType = GridHermitian<Real, trait::boundary_dimension>;
   mutable std::vector<BufferType> source_buffer;
   mutable std::vector<BufferType> disp_buffer;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation */
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 template <typename Func>
 void VolumePotential<type>::fourierApplyIf(Func func, GridBase<Real>& in,
                                            GridBase<Real>& out,
                                            filter_t pred) 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");
   // ^^^^ not applicable for Boussinesq
 
   // Transforming source
   for (UInt layer : Loop::range(i.sizes().front())) {
     if (not pred(layer))
       continue;
 
     auto in_layer = make_view(i, layer);
     FFTPlanManager::get()
         .createPlan(in_layer, source_buffer[layer])
         .forwardTransform();
   }
 
   // Applying operator
   func(source_buffer, disp_buffer);
 
   // Transforming output
   for (UInt layer : Loop::range(o.sizes().front())) {
     auto out_layer = make_view(o, layer);
     FFTPlanManager::get()
         .createPlan(out_layer, disp_buffer[layer])
         .backwardTransform();
   }
 }
 
 __END_TAMAAS__
 
 #endif  // VOLUME_POTENTIAL_HH