diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 431a4e5..a8e2b51 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,254 +1,252 @@
 /**
  * @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 "fft_plan_manager.hh"
 #include "grid_view.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "static_types.hh"
 
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type mtype, IntegralOperator::kind otype>
 Westergaard<mtype, otype>::Westergaard(Model* model) : IntegralOperator(model) {
   // Copy sizes
   std::array<UInt, trait::boundary_dimension> sizes;
   auto bdisc = model->getBoundaryDiscretization();
 
   // Ignoring first dimension if model is volumetric
   std::copy(bdisc.begin(), bdisc.end(), sizes.begin());
 
   auto boundary_hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           sizes);
 
   constexpr UInt nb_components = trait::components;
 
   buffer.setNbComponents(nb_components);
   buffer.resize(boundary_hermitian_sizes);
   influence.setNbComponents(nb_components * nb_components);
   influence.resize(boundary_hermitian_sizes);
   initInfluence();
 }
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 template <UInt m, UInt j>
 struct boundary_fft_helper {
   template <typename Buffer, typename Out>
   static void backwardTransform(Buffer&& buffer, Out&& out) {
     auto view = make_view(out, 0);
     FFTPlanManager::get().createPlan(view, buffer).backwardTransform();
   }
 };
 
 template <UInt m>
 struct boundary_fft_helper<m, m> {
   template <typename Buffer, typename Out>
   static void backwardTransform(Buffer&& buffer, Out&& out) {
     FFTPlanManager::get().createPlan(out, buffer).backwardTransform();
   }
 };
 }
 
 template <model_type mtype, IntegralOperator::kind otype>
 template <typename Functor>
 void Westergaard<mtype, otype>::fourierApply(Functor func, GridBase<Real>& in,
                                              GridBase<Real>& out) const try {
   Grid<Real, bdim>& i = dynamic_cast<decltype(i)>(in);
   Grid<Real, dim>& full_out = dynamic_cast<decltype(full_out)>(out);
 
   FFTPlanManager::get().createPlan(i, buffer).forwardTransform();
 
   /// Applying influence
   func(buffer, influence);
 
   /// Backward fourier transform on boundary only
   detail::boundary_fft_helper<bdim, dim>::backwardTransform(buffer, full_out);
 
 } catch (const std::bad_cast& e) {
   TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type mtype, IntegralOperator::kind otype>
 template <typename Functor>
 void Westergaard<mtype, otype>::initFromFunctor(Functor func) {
   // Compute wavevectors for influence
   auto wavevectors = FFTransform<Real, bdim>::template computeFrequencies<true>(
       influence.sizes());
   // Get boundary physical size
   auto system_size = this->model->getBoundarySystemSize();
   VectorProxy<const Real, bdim> domain(system_size[0]);
-  for (auto x : domain)
-    std::cout << x << std::endl;
 
   // Normalize wavevectors
   wavevectors *= 2 * M_PI;
   wavevectors /= domain;
 
   // Apply functor
   Loop::stridedLoop(func, wavevectors, influence);
 
   // Set fundamental mode to zero
   MatrixProxy<Complex, trait::components, trait::components> mat(influence(0));
   mat = 0;
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define NEUMANN_BASIC(type)                                                    \
   template <>                                                                  \
   void Westergaard<type, IntegralOperator::neumann>::initInfluence() {         \
     auto E_star = model->getHertzModulus();                                    \
     auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> && q,            \
                                       Complex & k) {                           \
       k = 2. / (E_star * q.l2norm());                                          \
     };                                                                         \
     initFromFunctor(basic);                                                    \
   }
 #define DIRICHLET_BASIC(type)                                                  \
   template <>                                                                  \
   void Westergaard<type, IntegralOperator::dirichlet>::initInfluence() {       \
     auto E_star = model->getHertzModulus();                                    \
     auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> && q,            \
                                       Complex & k) {                           \
       k = E_star * q.l2norm() / 2;                                             \
     };                                                                         \
     initFromFunctor(basic);                                                    \
   }
 
 NEUMANN_BASIC(model_type::basic_1d);
 NEUMANN_BASIC(model_type::basic_2d);
 DIRICHLET_BASIC(model_type::basic_1d);
 DIRICHLET_BASIC(model_type::basic_2d);
 
 #undef NEUMANN_BASIC
 #undef DIRICHLET_BASIC
 /* -------------------------------------------------------------------------- */
 
 template <>
 void Westergaard<model_type::surface_2d,
                  IntegralOperator::neumann>::initInfluence() {
   auto E = model->getYoungModulus();
   auto nu = model->getPoissonRatio();
   const Complex I(0, 1);
   auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> && q,
                                         MatrixProxy<Complex, comp, comp> && F) {
     Real q_norm = q.l2norm();
     q /= q_norm;
 
     F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0));
     F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1));
     F(2, 2) = 2 * (1 - nu * nu);
 
     F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu);
 
     F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu);
     F(1, 2) = I * q(1) * (1 + nu) * (1 - 2 * nu);
 
     F(2, 0) = -F(0, 2);
     F(2, 1) = -F(1, 2);
 
     F *= 1. / (E * q_norm);
   };
 
   initFromFunctor(surface);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type mtype, IntegralOperator::kind otype>
 void Westergaard<mtype, otype>::initInfluence() {}
 
 /* ------------------------------------------------------------------------ */
 template <model_type mtype, IntegralOperator::kind otype>
 void Westergaard<mtype, otype>::apply(GridBase<Real>& input,
                                       GridBase<Real>& output) const {
   auto apply = [](decltype(buffer)& buffer,
                   const decltype(influence)& influence) {
     Loop::stridedLoop([] CUDA_LAMBDA(VectorProxy<Complex, comp> && i,
                                      MatrixProxy<const Complex, comp, comp> &&
                                          F) { i = F * i; },
                       buffer, influence);
   };
   fourierApply(apply, input, output);
 }
 
 /* -------------------------------------------------------------------------- */
 template <model_type type, IntegralOperator::kind otype>
 Real Westergaard<type, otype>::getOperatorNorm() {
   constexpr UInt comp = model_type_traits<type>::components;
   Real _norm = 0.0;
 
   _norm = Loop::stridedReduce<operation::plus>(
     [] CUDA_LAMBDA(MatrixProxy<Complex, comp, comp>&& m) {
       Real n = thrust::norm(m.l2squared());
       return std::isnan(n) ? 0 : n;
     },
     influence);
 
   auto size = model->getSystemSize();
   auto disc = model->getDiscretization();
   auto dim = model_type_traits<type>::dimension;
 
   /// TODO: why?
   switch (dim) {
   case 1:
     _norm /= (size[0] / disc[0]) * (size[0] / disc[0]);
     break;
   default:
     for (UInt i = 0; i < dim; i++) {
       _norm /= size[i] / disc[i];
     }
   }
 
   return std::sqrt(_norm);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation */
 /* -------------------------------------------------------------------------- */
 template class Westergaard<model_type::basic_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::basic_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::basic_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::basic_2d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::surface_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::surface_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::surface_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::surface_2d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::volume_1d, IntegralOperator::neumann>;
 template class Westergaard<model_type::volume_2d, IntegralOperator::neumann>;
 template class Westergaard<model_type::volume_1d, IntegralOperator::dirichlet>;
 template class Westergaard<model_type::volume_2d, IntegralOperator::dirichlet>;
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */