diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index ad709f8..7604742 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,238 +1,238 @@
 /**
  * @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 "westergaard.hh"
 #include "fft_plan_manager.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Westergaard<type>::Westergaard(Model* model) : BEEngine(model) {
   // Copy sizes
   std::array<UInt, trait::dimension> sizes;
   std::copy(model->getDiscretization().begin(),
             model->getDiscretization().end(), sizes.begin());
   auto hermitian_sizes =
       GridHermitian<Real, trait::dimension>::hermitianDimensions(sizes);
 
   // Copy boundary sizes
   auto sizes_it = sizes.begin();
   std::array<UInt, trait::boundary_dimension> boundary_hermitian_sizes;
   // Ignoring first dimension if model is volumetric
   if (trait::dimension > trait::boundary_dimension)
     ++sizes_it;
   std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin());
   boundary_hermitian_sizes =
       GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
           boundary_hermitian_sizes);
 
   constexpr UInt nb_components = trait::components;
 
   buffer.setNbComponents(nb_components);
   buffer.resize(boundary_hermitian_sizes);
   influence.setNbComponents(nb_components * nb_components);
   influence.resize(hermitian_sizes);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::initInfluence() {
   TAMAAS_EXCEPTION("Not implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define INFLUENCE_BASIC_MACRO                                                  \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E_star = model->getHertzModulus();                                    \
     VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
     Loop::stridedLoop(influence::Basic<bdim>(E_star, domain), wavevectors,     \
                       influence);                                              \
     influence(0) = 0;                                                          \
   } while (0)
 
 template <>
 void Westergaard<model_type::basic_2d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 template <>
 void Westergaard<model_type::basic_1d>::initInfluence() {
   INFLUENCE_BASIC_MACRO;
 }
 
 #undef INFLUENCE_BASIC_MACRO
 
 #define INFLUENCE_SURFACE_MACRO                                                \
   do {                                                                         \
     constexpr UInt bdim = trait::boundary_dimension;                           \
     auto wavevectors =                                                         \
         FFTransform<Real, bdim>::computeFrequencies<true>(influence.sizes());  \
     Real E = model->getYoungModulus();                                         \
     Real nu = model->getPoissonRatio();                                        \
     VectorProxy<const Real, bdim> domain(model->getSystemSize()[0]);          \
     Loop::stridedLoop(influence::Surface<bdim>(E, nu, domain), wavevectors,    \
                       influence);                                              \
   } while (0)
 
 template <>
 void Westergaard<model_type::surface_1d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 template <>
 void Westergaard<model_type::surface_2d>::initInfluence() {
   INFLUENCE_SURFACE_MACRO;
 }
 
 #undef INFLUENCE_SURFACE_MACRO
 
 /* ------------------------------------------------------------------------ */
 
 #define NEUMANN_BASIC(type)                                                    \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer *= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 #define DIRICHLET_BASIC(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveDirichlet(GridBase<Real>& in,                   \
                                          GridBase<Real>& out) const {          \
     auto apply = [](decltype(buffer)& buffer,                                  \
                     const decltype(influence)& influence) {                    \
       buffer /= influence;                                                     \
       buffer(0) = 0;                                                           \
     };                                                                         \
                                                                                \
     fourierApply(apply, in, out);                                              \
   }
 
 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
 
 /* -------------------------------------------------------------------------- */
 
 namespace detail {
 /// Class needed for cuda
 template <UInt d>
 struct SurfaceApply {
   void operator()(GridHermitian<Real, d>& buffer_,
                   const GridHermitian<Real, d>& influence) {
     Loop::stridedLoop(
         [] CUDA_LAMBDA(VectorProxy<Complex, d + 1> && buf,
                        MatrixProxy<const Complex, d + 1, d + 1> && inf) {
           Vector<Complex, d + 1> copy(buf);
           buf.mul(inf, copy);
         },
         buffer_, influence);
 
     VectorProxy<Complex, d + 1> buff_0(buffer_(0));
     buff_0 = 0;
   }
 };
 }  // namespace detail
 
 #define NEUMANN_SURFACE(type)                                                  \
   template <>                                                                  \
   void Westergaard<type>::solveNeumann(GridBase<Real>& in,                     \
                                        GridBase<Real>& out) const {            \
     fourierApply(detail::SurfaceApply<trait::boundary_dimension>(), in, out);  \
   }
 
 /// \TODO Find a way to make it work with cuda
 NEUMANN_SURFACE(model_type::surface_1d);
 NEUMANN_SURFACE(model_type::surface_2d);
 
 #undef NEUMANN_SURFACE
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveNeumann(GridBase<Real>& /*neumann*/,
                                      GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Westergaard<type>::solveDirichlet(GridBase<Real>& /*neumann*/,
                                        GridBase<Real>& /*dirichlet*/) const {
   TAMAAS_EXCEPTION("Not yet implemented");
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Real Westergaard<type>::getNeumannNorm() {
   constexpr UInt comp = model_type_traits<type>::components;
   Real _norm = 0.0;
 
   _norm = Loop::stridedReduce<op::plus>(
-    [] CUDA_LAMBDA(StaticVector<Complex, comp * comp>&& m) {
+    [] CUDA_LAMBDA(MatrixProxy<Complex, comp, comp>&& m) {
       Real n = thrust::norm(m.l2squared());
       return std::isnan(n) ? 0 : n;
     },
     influence);
 
   return std::sqrt(_norm);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class Westergaard<model_type::basic_1d>;
 template class Westergaard<model_type::basic_2d>;
 template class Westergaard<model_type::surface_1d>;
 template class Westergaard<model_type::surface_2d>;
 template class Westergaard<model_type::volume_1d>;
 template class Westergaard<model_type::volume_2d>;
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp
index d61ea66..7ec0b3d 100644
--- a/src/solvers/kato.cpp
+++ b/src/solvers/kato.cpp
@@ -1,366 +1,366 @@
 /**
  * @file
  *
  * @author Son Pham-Ba <son.phamba@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2018 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 "kato.hh"
 #include "elastic_functional.hh"
 #include "loop.hh"
 #include "fft_plan_manager.hh"
 #include <iomanip>
 #include <iterator>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Kato::Kato(Model& model, const GridBase<Real>& surface, Real tolerance,
            Real mu)
     : ContactSolver(model, surface, tolerance),
       engine(model.getBEEngine()), mu(mu) {
   if (model.getType() != model_type::surface_1d &&
       model.getType() != model_type::surface_2d) {
     TAMAAS_EXCEPTION(
         "Model type is not compatible with Kato solver");
   }
 
   gap = this->_gap.get(); // locally allocated
   pressure = &model.getTraction();
   comp = pressure->getNbComponents();
   N = pressure->getNbPoints();
 
   if (model.getType() == model_type::surface_1d) {
     initSurfaceWithComponents<model_type::surface_1d>();
   } else {
     initSurfaceWithComponents<model_type::surface_2d>();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Kato::solve(GridBase<Real>& p0) {
   if (p0.getNbPoints() != comp) {
     TAMAAS_EXCEPTION(
       "Target mean pressure does not have the right number of components");
   }
 
   Real cost = 0;
 
   switch (model.getType()) {
   case model_type::surface_1d:
     cost = solveTmpl<model_type::surface_1d>(p0);
     break;
   case model_type::surface_2d:
     cost = solveTmpl<model_type::surface_2d>(p0);
     break;
   default:
     break;
   }
 
   return cost;
 }
 
 template <model_type type>
 Real Kato::solveTmpl(GridBase<Real>& p0) {
   Real cost = 0;
   UInt n = 0;
 
   // Printing column headers
   std::cout << std::setw(5) << "Iter"
             << " " << std::setw(15) << "Cost_f"
             << " " << std::setw(15) << "Error" << '\n'
             << std::fixed;
 
   pressure->uniformSetComponents(p0);
 
   // setMaxIterations(10);
   // setDumpFrequency(1);
 
   do {
     engine.solveNeumann(*pressure, *gap);
     *gap -= surfaceComp;
     *pressure -= *gap;
 
     enforcePressureConstraints<type>(p0);
     cost = computeCost();
     printState(n, cost, cost);
   } while (cost > this->tolerance && n++ < this->max_iterations);
 
   Real beta = computeBeta<type>();
   computeFinalGap<type>(beta);
   model.getDisplacement() = *gap;
 
   return cost;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Kato::solveRelaxed(GridBase<Real>& g0) {
   if (g0.getNbPoints() != comp) {
     TAMAAS_EXCEPTION(
       "Target mean gap does not have the right number of components");
   }
 
   Real cost = 0;
 
   switch (model.getType()) {
   case model_type::surface_1d:
     cost = solveRelaxedTmpl<model_type::surface_1d>(g0);
     break;
   case model_type::surface_2d:
     cost = solveRelaxedTmpl<model_type::surface_2d>(g0);
     break;
   default:
     break;
   }
 
   return cost;
 }
 
 template <model_type type>
 Real Kato::solveRelaxedTmpl(GridBase<Real>& g0) {
   constexpr UInt comp = model_type_traits<type>::components;
   Real cost = 0;
   UInt n = 0;
 
   // Printing column headers
   std::cout << std::setw(5) << "Iter"
             << " " << std::setw(15) << "Cost_f"
             << " " << std::setw(15) << "Error" << '\n'
             << std::fixed;
 
   *pressure = 0;
 
   do {
     engine.solveNeumann(*pressure, *gap);
     addUniform<comp>(*gap, g0);
     *gap -= surfaceComp;
     *pressure -= *gap;
 
     enforcePressureCoulomb<type>();
     cost = computeCost();
     printState(n, cost, cost);
   } while (cost > this->tolerance && n++ < this->max_iterations);
 
   *gap += surfaceComp;
   model.getDisplacement() = *gap;
 
   return cost;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::initSurfaceWithComponents() {
   constexpr UInt comp = model_type_traits<type>::components;
   surfaceComp.setNbComponents(comp);
   surfaceComp.resize({N});
   surfaceComp = 0;
 
   Loop::stridedLoop(
-    [] CUDA_LAMBDA(Real& s, StaticVector<Real, comp>&& sc) {
+    [] CUDA_LAMBDA(Real& s, VectorProxy<Real, comp>&& sc) {
       sc(comp - 1) = s;
     },
     surface, surfaceComp);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Projects $\vec{p}$ on $\mathcal{C}$ and $\mathcal{D}$.
  */
 template <model_type type>
 void Kato::enforcePressureConstraints(GridBase<Real>& p0) {
   UInt n = 50;
 
   for (UInt i = 0; i < n; i++) {
     enforcePressureMean<type>(p0);
     enforcePressureCoulomb<type>();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::enforcePressureMean(GridBase<Real>& p0) {
   constexpr UInt comp = model_type_traits<type>::components;
   Grid<Real, 1> corr({1}, comp);
 
   computeMean<comp>(*pressure, corr);
   corr -= p0;
 
-  StaticVector<Real, comp> corr_(corr(0));
+  VectorProxy<Real, comp> corr_(corr(0));
   Loop::stridedLoop(
-    [corr_] CUDA_LAMBDA(StaticVector<Real, comp>&& p) {
+    [corr_] CUDA_LAMBDA(VectorProxy<Real, comp>&& p) {
       p -= corr_;
     },
     *pressure);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::enforcePressureCoulomb() {
   constexpr UInt comp = model_type_traits<type>::components;
 
   Loop::stridedLoop(
-    [this] CUDA_LAMBDA(StaticVector<Real, comp>&& p) {
-      StaticVector<Real, comp - 1> p_T(p(0));
+    [this] CUDA_LAMBDA(VectorProxy<Real, comp>&& p) {
+      VectorProxy<Real, comp - 1> p_T(p(0));
       Real p_N = p(comp - 1);
       Real p_T_sqrd= p_T.l2squared();
 
       bool cond1 = (p_N >= 0 && p_T_sqrd <= mu * mu * p_N * p_N);
       bool cond2 = (p_N <= 0 && p_T_sqrd <= p_N * p_N / mu / mu);
 
       if (cond2) {
         p_T = 0;
         p(comp - 1) = 0;
       } else if (!cond1) {
         Real p_T_norm = std::sqrt(p_T_sqrd);
         Real k = (p_N + mu * p_T_norm) / (1 + mu * mu);
         p_T *= k * mu / p_T_norm;
         p(comp - 1) = k;
       }
     },
     *pressure);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Compute mean of the field taking each component separately.
  */
 template <UInt comp>
 void Kato::computeMean(GridBase<Real>& field, GridBase<Real>& mean) {
+  auto _mean = Loop::stridedReduce<op::plus>(
+    [] CUDA_LAMBDA(VectorProxy<Real, comp>&& f) -> Vector<Real, comp> {
+      return f;
+    },
+    field);
+
   for (UInt i = 0; i < comp; i++) {
-    mean(i) = Loop::stridedReduce<op::plus>(
-      [i] CUDA_LAMBDA(StaticVector<Real, comp>&& f) {
-        return f(i);
-      },
-      field);
+    mean(i) = _mean(i) / N;
   }
-
-  mean /= N;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt comp>
 void Kato::addUniform(GridBase<Real>& field, GridBase<Real>& vec) {
-  for (UInt i = 0; i < comp; i++) {
-    Loop::stridedLoop(
-      [&] CUDA_LAMBDA(StaticVector<Real, comp>&& f) {
+  Loop::stridedLoop(
+    [&] CUDA_LAMBDA(VectorProxy<Real, comp>&& f) {
+      for (UInt i = 0; i < comp; i++) {
         f(i) += vec(i);
-      },
-      field);
-  }
+      }
+    },
+    field);
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Kato::computeCost() {
   UInt N = pressure->getNbPoints();
   Real beta = 0;
   Grid<Real, 1> lambda({N}, 1);
   Grid<Real, 1> eta({N}, 1);
   Grid<Real, 1> p_N({N}, 1);
   Grid<Real, 1> p_C({N}, 1);
 
   switch (model.getType()) {
   case model_type::surface_1d:
     beta = computeBeta<model_type::surface_1d>();
     computeValuesForCost<model_type::surface_1d>(beta, lambda, eta, p_N, p_C);
     break;
   case model_type::surface_2d:
     beta = computeBeta<model_type::surface_2d>();
     computeValuesForCost<model_type::surface_2d>(beta, lambda, eta, p_N, p_C);
     break;
   default:
     break;
   }
 
   return p_N.dot(lambda) + p_C.dot(eta);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 Real Kato::computeBeta() {
   constexpr UInt comp = model_type_traits<type>::components;
 
   return Loop::stridedReduce<op::max>(
-    [this] CUDA_LAMBDA(StaticVector<Real, comp>&& g) {
-      StaticVector<Real, comp - 1> g_T(g(0));
+    [this] CUDA_LAMBDA(VectorProxy<Real, comp>&& g) {
+      VectorProxy<Real, comp - 1> g_T(g(0));
       Real g_N = g(comp - 1);
       Real g_T_norm = g_T.l2norm();
       return mu * g_T_norm - g_N;
     },
     *gap);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::computeValuesForCost(Real beta, GridBase<Real>& lambda, GridBase<Real>& eta,
                          GridBase<Real>& p_N, GridBase<Real>& p_C) {
   constexpr UInt comp = model_type_traits<type>::components;
 
   Loop::stridedLoop(
-    [this, beta] CUDA_LAMBDA(StaticVector<Real, comp>&& p,
-                             StaticVector<Real, comp>&& g,
+    [this, beta] CUDA_LAMBDA(VectorProxy<Real, comp>&& p,
+                             VectorProxy<Real, comp>&& g,
                              Real& lambda_,
                              Real& eta_,
                              Real& p_N_,
                              Real& p_C_) {
-      StaticVector<Real, comp - 1> g_T(g(0));
+      VectorProxy<Real, comp - 1> g_T(g(0));
       Real g_N = g(comp - 1);
       Real g_T_norm = g_T.l2norm();
       lambda_ =  g_N - mu * g_T_norm + beta;
       eta_ = g_T_norm;
 
-      StaticVector<Real, comp - 1> p_T(p(0));
+      VectorProxy<Real, comp - 1> p_T(p(0));
       Real p_N = p(comp - 1);
       Real p_T_norm = p_T.l2norm();
       p_N_ = p(comp - 1);
       p_C_ = mu * p_N - p_T_norm;
     },
     *pressure, *gap, lambda, eta, p_N, p_C);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::computeFinalGap(Real beta) {
   constexpr UInt comp = model_type_traits<type>::components;
 
   Loop::stridedLoop(
-    [beta] CUDA_LAMBDA(StaticVector<Real, comp>&& g,
-                       StaticVector<Real, comp>&& h) {
+    [beta] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
+                       VectorProxy<Real, comp>&& h) {
       g(comp - 1) += h(comp - 1) + beta;
     },
     *gap, surfaceComp);
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */